Research on the Competition Mechanism of Fractures in Multi-Cluster Fracturing of Horizontal Wells: Dynamic Response and Influence of Engineering Parameters

1  Introduction

The development of unconventional oil and gas resources is of strategic significance to global energy security. Horizontal well multi-cluster fracturing is currently the key technology for exploiting low-permeability reservoirs such as shale oil and tight gas. At present, sequential fracturing is predominantly employed to generate multiple fractures, thereby maximizing the stimulated reservoir volume. However, downhole monitoring indicates that fractures generally exhibit non-uniform propagation, resulting in ineffective clusters in certain sections and significantly compromising stimulation outcomes. Unraveling the competitive mechanisms of fracture growth and achieving uniform fracture propagation represent the core challenges in optimizing current fracturing techniques.

In recent years, numerical simulation methods have played a crucial role in studying fracturing mechanisms. For instance, Yin et al. 2025 [1] developed a fluid-solid coupling model incorporating real gas characteristics to analyze induced stress evolution during hydraulic fracturing and production, highlighting the complex spatiotemporal changes in the stress field. Wang et al. 2025 [2] utilized the lattice method to investigate multilayer combined fracturing in coal-bearing strata, specifically analyzing the influence of perforation design on fracture propagation, demonstrating the method’s capability in handling complex initiation scenarios. These studies provide valuable insights into stress evolution and near-wellbore fracture complexity, which are relevant to the cluster interaction mechanisms explored in this work.

Currently, significant theoretical and modeling advancements have been achieved across various areas, including improved cohesive-zone models, coupled DDM-FEM algorithms, novel CFD-DEM models, and fracturing simulations under mining-induced stress fields. These advancements have enhanced the accuracy of fracture-propagation predictions, computational efficiency, proppant-transport modeling, and simulations in rotational stress fields [3,4]. Damjanac and Cundall 2016 [5,6], proposed a fully coupled hydro-mechanical particle-based model, validating its effectiveness through analytical solutions and field cases. Huang et al. 2024 [7], analyzed the competitive propagation of multi-cluster fractures using the finite element–discrete element method. Using this model, a stochastic assignment procedure was developed to establish a coupled distribution of mechanical parameters between matrix elements and discrete elements, characterizing reservoir heterogeneity. Huang et al. 2020 [8], introduced three-dimensional mesh modeling for hydraulic fracture initiation and near-wellbore propagation, considering different perforation models including helical perforation, directional perforation, and Tristim perforation. Guo et al. 2024 [9], introduced a modeling workflow based on the finite element method and displacement discontinuity method for simulating dynamic porous media flow, geomechanics, and hydraulic fracturing. Zhao et al. 2020 [10], revealed through fluid-solid coupling simulations, the combined effects of multiple factors on fracture interactions. Cong et al. 2022 [11], investigated fracture height growth in layered formations, emphasizing the controlling roles of in-situ stress and lithology. Suo et al. 2026 [12], integrated theoretical analysis with finite element–discrete element method (FDEM) simulations to investigate the mixed-mode mechanisms, fracture-tip plastic zone evolution, and anisotropic mechanical responses in shale formations.

Existing studies have deepened the understanding of the influencing factors, yet most have focused on analyzing static parameters such as final fracture geometry or breakdown pressure. There remains a relative lack of systematic research on the dynamic evolution of pressure behavior within each cluster during fracturing. Key questions—such as how cluster pressures sequentially respond to fracture initiation, how they differentiate under stress shadowing, and which critical parameters control their dynamic trajectories—reflect the essence of competitive fracture propagation but remain poorly elucidated. This knowledge gap limits the precision of optimization strategies. To address this, this study employs the discrete lattice method to simulate simultaneous multi-cluster fracturing in a horizontal well with three clusters. The research focus shifts from static outcomes to dynamic processes, enabling real-time tracking of inlet fluid pressure, flow rate, and fluid intake volume per perforation cluster. By specifically examining the effects of cluster spacing, injection rate, and horizontal stress difference, this work aims to decipher the dynamic pressure evolution and uncover the initiation and development mechanisms of fracture competition at a fundamental level. The findings are expected to provide a theoretical basis for establishing optimization methods grounded in dynamic pressure characteristics [13–17], offering practical guidance for improving the uniformity of reservoir stimulation.

2  Discrete Lattice Method Numerical Model

2.1 Mechanical Model

The Discrete Lattice Method is an explicit numerical method capable of effectively simulating highly nonlinear processes, such as sliding, opening/closing, and fracturing, between discrete blocks. Its characteristic is to represent an object as an assemblage of blocks connected by nodes, with springs linking the nodes. Each node possesses six degrees of freedom, including three translational degrees governed by momentum conservation equations and three rotational degrees determined by Newton’s second law. The motion of translational degrees of freedom follows a central-difference formulation applied to each node [18,19]:

{vi(t+Δt2)=vi(t−Δt2)+∑Fi(t)Δtmui(t+Δt)=ui(t)+vi(t+Δt2)Δt(1)

where vit is the velocity component of node i (i = 1, 2, 3) at time t, m/s; uit is the position component of node i (i = 1, 2, 3) at time t, m; ∑Fi(t) is the sum of all force components acting on the node with mass m, N; Δt is the time step, s; m is the nodal mass, kg.

The central difference formula for the angular velocity in the rotational degrees of freedom of node i (1, 2, 3) at time t is as follows:

ωi(t+Δt2)=ωi(t−Δt2)+ΣMi(t)IΔt(2)

where ΣMi(t) is the sum of the moments of all i—components at node at time t, N/m; I is the moment of inertia, kg·m2; ωi is the angular velocity, rad/s.

For an unbroken spring, the relative velocity between the connected nodes A and B is:

vir=viA−viB(3)

The normal and tangential velocity components are:

{viN=virniviS=vir−viNni(4)

where N represents the normal direction; S represents the tangential direction; ni is the unit normal vector. The Einstein summation convention applies to repeated indices. The normal and tangential stresses of the spring are updated based on the relative displacement of the nodes:

{FiN←FiN+viNKNΔtFiS←FiS+viSKSΔt(5)

In the formula, KN and KS are, respectively, the normal stiffness and tangential stiffness of the spring, N/m; Fi(t)N and Fi(t)S represent the normal force and tangential force components (i = 1, 2, 3) of the node at time t, respectively, N; vi(t)N and vi(t)S represent the normal displacement and tangential displacement components (i = 1, 2, 3) of the node at time t, respectively, m.

The relationship between the spring stiffness, tensile strength and the bulk modulus, tensile strength of the rock mass is given by [20]:

{FHmax=αtσbR2FSmax=λFHmax+αsσcR2(6)

where FHmax and FSmax represent the tensile strength and shear strength of the spring, respectively, N; αt and αs are the tensile strength correction factor and shear strength correction factor, respectively; σb is the macroscopic rock mass tensile strength, Pa; σc is the macroscopic rock mass shear strength, Pa; R is the half-length of the fracture, m; λ is the friction coefficient.

2.2 Flow Field Model

This study utilizes the Discrete Lattice Method (DLM) implemented in the Xsite software platform for numerical simulation. DLM exhibits distinct advantages in modeling discontinuous medium deformation and the evolution of complex fracture networks. By discretizing the rock mass into a system of spring-connected mass points, the method inherently captures non-planar fracture initiation, multi-directional branching, and dynamic reorientation behaviors without requiring predefined fracture propagation paths. This capability makes DLM particularly well-suited for characterizing competitive inter-cluster propagation phenomena driven by stress shadow effects in multi-cluster fracturing. In contrast to continuum-based numerical approaches, DLM provides enhanced physical rationality in simulating fracture initiation and discontinuous deformation processes. The adoption of this method ensures that the numerical model accurately represents fracture interaction mechanisms, thereby establishing a reliable foundation for parameter sensitivity analysis.

When the force on a spring exceeds its predefined strength, the spring fractures and generates a micro-crack. Adjacent microcracks form flow channels, and the fluid units at the centers of the springs are interconnected through these channels. Newly formed micro-cracks are automatically integrated into the channel network, enabling real-time updates to the network [19]. A schematic diagram of the lattice network and fluid channels is shown in Fig. 1.

Figure 1: Schematic diagram of the lattice network and fluid channels.

The flow of fluid within the channels is described by the lubrication equation, assuming the channel length and width are equal. The volumetric flow rate of fluid moving from node A to node B through a channel is given by the expression:

q=βkrw312μ[pA−pB+ρwg(zA−zB)](7)

where q is fluid flow rate, m3/s; β is dimensionless coefficient; w is fracture width, m; μ is fluid viscosity, Pa·s; pA and pB represent the fluid pressures at nodes A and B, respectively, Pa; zA and zB represent the hydraulic heads at nodes A and B, respectively, m; ρw is fluid density, kg/m3; g is Gravitational acceleration, m/s2; kr is the relative permeability, which is a function of saturation s:

kr=s2(3−2s)(8)

The flow evolution model employs an explicit numerical solution method, where the fluid pressure increment per unit flow time step is given by [20]:

Δp=∑qiVKFΔtf(9)

In the formula, Δp is fluid pressure increment, Pa; qi is flow rate in the pipeline connected to node i, m3/s; V is node volume, m3; KF is apparent fluid bulk modulus, Pa; Δtf is unit fluid time step, s.

2.3 Hydraulic Fracture Propagation Criterion

The propagation of hydraulic fractures in this model is determined based on the J-integral criterion. Under quasi-static conditions, assuming the presence of body forces and tractions on the fracture surface, the J-integral formulation for calculating the local energy release rate along the unit vector at the fracture front is expressed as:

J=1R∫V[(σijuij−Wδki)q,ir^nk−Fiui,jq,ir^nj]dV−1R∫STjuj,kqr^nkdS(10)

{qr^=1−r^,|r^|<1qr^=0,     |r^|≥1(11)

{r^=x^12+x^22+x^32x^i=xi−xicR(12)

where W is strain energy density, J/m3; V is region containing the fracture front, m3; σ is cauchy stress tensor, Pa; u is displacement vector, m; S is Fracture surface, m2; qr^ is smooth function in V; F and T represent the body force and traction force, respectively, N; xi and xic represent the coordinates of an arbitrary point and the domain center, respectively.

The stress intensity factor based on the J-integral is expressed as:

KI=JE(13)

In the formula, KI is the stress intensity factor, Pa·m0.5; E is the Young’s modulus, Pa.

2.4 Numerical Model Validation

The model dimensions are 40 m (length) × 20 m (width) × 30 m (height). A horizontal wellbore is placed along the model’s longitudinal axis and includes three perforation clusters. The initial in-situ stress field consists of vertical stress, maximum horizontal principal stress, and minimum horizontal principal stress, with specific values listed in Table 1. Each perforation cluster inlet is defined as a constant-rate injection boundary to simulate fracturing fluid injection, and the initial pore pressure is assigned as the reservoir pore pressure [21].

To validate the accuracy of the discrete lattice method for numerical simulation, a comparison was conducted between the geometric dimensions of a simulated 3D penny-shaped fracture and the theoretical solution under identical geological conditions and engineering parameters [22–24]. Fig. 2a illustrates the simulated morphology of the penny-shaped fracture, while Fig. 2b compares the numerically simulated fracture geometry with the theoretical results, showing generally good agreement [25,26]. The maximum fracture width and length obtained from the numerical simulation are 0.00112 and 24.36 m, respectively, compared to the theoretical values of 0.00124 and 23.82 m. The relative errors for fracture width and length are 9.7% and 2.2%, respectively. At the fracture initiation point, the theoretical result is larger than the simulated value, whereas at the fracture tip, the theoretical result is smaller. This observation is consistent with findings from previous studies. In conclusion, the numerical model is suitable for subsequent research and evaluation [27].

Figure 2: (a) Simulated fracture geometry, (b) comparison of dimensions between numerical and theoretical results.

The model validation methodology employed in this study demonstrates both academic rigor and engineering applicability. Although the ultimate engineering application focuses on the competitive propagation of multi-cluster fractures, the use of a single-fracture model for validation represents a standard methodological approach and a necessary foundation in computational fracture mechanics [28]. The theoretical basis for this lies in the understanding that the complex behavior of multi-cluster fracture competition fundamentally arises from the coupling of two essential physical processes: the “single-fracture propagation mechanism” and the “stress interference effects between clusters”. By systematically comparing the numerical simulation results of a three-dimensional penny-shaped fracture with classical theoretical solutions, this study first verifies the accuracy of the model in characterizing the core physical process of single-fracture initiation and propagation within a homogeneous stress field [29,30]. The reliability of stress shadow effects, as an inherent coupling mechanism of the model, is based on the model’s ability to represent fracture deformation behavior and the resulting induced stress field. The single-fracture validation provides indirect yet robust mechanical evidence for this capability. Therefore, the validation at the single-fracture level sufficiently meets the requirements for establishing the fundamental reliability of the model in this study [31,32].

3  Simulation Results and Analysis

3.1 Effect of Cluster Spacing

By comparing the three-dimensional fracture geometries at cluster spacings of 4, 3, and 2 m (Fig. 3), the significant influence of cluster spacing on fracture dimensions, propagation paths, and uniformity is clearly evident. The underlying mechanism primarily stems from variations in the intensity of the “stress shadow effect.” During simultaneous multi-cluster fracturing, the initially initiated or dominantly propagating outer-cluster fractures generate a high-stress zone around them. The middle cluster, situated within the overlapping stress shadow regions of the left and right outer clusters, experiences a substantial increase in stress during its propagation [24].

Figure 3: Effect of cluster spacing on fracture geometry.

From the apparent characteristics of the cloud charts, the fracture geometry undergoes systematic changes as the cluster spacing decreases: the maximum fracture width is reduced by 3.92 mm, and the volume of the middle-cluster fracture decreases multiplicatively. At a 4 m spacing, all three clusters achieve practical propagation. Although there are specific differences in fracture length and width, the overall geometries remain relatively independent of each other. When the spacing is reduced to 3 m, the interaction between fractures becomes evident, with the propagation of the middle-cluster fracture significantly suppressed, manifesting as reduced length and width compared to the outer clusters. At the minimal spacing of 2 m, the competition is most intense: the middle-cluster fracture is severely inhibited, exhibiting minimal propagation, while the two outer-cluster fractures receive more fluid and develop more fully, resulting in a final fracture geometry that shows pronounced non-uniformity. Based on the above discussion, by influencing the intensity of the stress shadow effect, cluster spacing becomes a key factor controlling fracture competition and propagation uniformity.

Figs. 4–6 clearly reveal the systematic influence of cluster spacing on fracture propagation. As the cluster spacing decreases from 4 to 2 m, competition among fractures intensifies, as reflected in dynamic changes in three parameters: fluid pressure, flow rate, and fluid intake volume. The fluid pressure curves indicate that the initiation pressure of the middle cluster is significantly delayed, and its steady-state pressure remains consistently higher than that of the outer clusters. The maximum inter-cluster fluid pressure difference reaches 7.32 MPa, directly reflecting the enhancement of the stress shadow effect. As the spacing decreases, not only does the peak pressure difference increase, but its duration also extends noticeably, indicating that the middle cluster faces greater propagation resistance throughout the entire treatment. The variation in flow rate curves is even more pronounced, showing severe uneven distribution. The initial peak flow rate of the middle cluster is only 60% of that of the outer clusters and decays more rapidly. At a 2 m spacing, the instantaneous flow rate of the middle cluster drops to less than 30% of that of the outer clusters at its lowest point. The cumulative fluid intake volume curves intuitively reflect the outcome of this competition. The stimulated volume of the middle cluster is significantly reduced, with a final difference in fluid intake volume reaching 2.6 times the baseline. Comprehensive analysis indicates that reduced spacing, by enhancing the stress shadow effect, adversely affects all three dimensions—fluid pressure, flow rate, and fluid intake volume—making it a key parameter controlling the intensity of fracture competition.

Figure 4: Effect of cluster spacing on fluid pressure.

Figure 5: Effect of cluster spacing on flow rate.

Figure 6: Effect of cluster spacing on fluid intake volume.

3.2 Effect of Fracturing Fluid Injection Rate

A comparison of the three-dimensional fracture geometries at injection rates of 0.02, 0.04, and 0.06 m3/s (Fig. 7) reveals that increasing the injection rate directly enhances the overall scale and complexity of the fracture network. The maximum fracture widths at the three rates are 7.35 mm, 1.01 cm, and 15.1 mm, respectively, showing a significant increasing trend.

Figure 7: Effect of injection rate on fracture geometry.

Under low injection rate conditions, fracture propagation is limited, and the geometry tends to be relatively straightforward. As the injection rate increases, not only does the length of the main fracture significantly increase, but the activation of secondary fractures is also markedly enhanced, demonstrating improved interconnectivity within the fracture network. This phenomenon stems from the direct control of the injection rate over the fluid pressure field within the fractures and the intensity of stress interference. According to the principle of mass conservation, an increase in injection rate implies a larger volume of fluid injected into the fractures per unit time, thereby increasing the net pressure within the fractures. The resulting higher internal pressure provides greater driving force for the fractures to overcome in-situ stresses and propagate further, while simultaneously increasing the likelihood of fluid infiltrating natural weak planes (such as natural fractures and bedding planes), thereby promoting the formation of a complex fracture network.

Furthermore, an increased injection rate increases the fluid’s kinetic energy, leading to a more balanced distribution of fluid across multiple clusters. This counteracts inter-cluster competition caused by the stress shadow effect and effectively mitigates the negative impacts associated with close cluster spacing. while a higher injection rate promotes fracture uniformity and complexity, it is crucial to consider its potential side effects. The significant increase in net pressure within the fractures, while providing the driving force for propagation, may approach or exceed the pressure limits of surface equipment and wellbore integrity. Additionally, the elevated pressure and enhanced fluid energy might facilitate the activation of natural fractures or bedding planes, potentially promoting the formation of a more complex micro-fracture network. This requires a balanced design considering equipment capabilities and the desired fracture complexity. In summary, the injection rate, by governing intrafracture pressure and fluid distribution efficiency, dictates the macroscopic scale and overall stimulated volume of the fracture system [20], making it a core engineering parameter for controlling fracture extension dimensions and complexity. When optimizing treatment design, ensuring a sufficiently high injection rate within equipment capabilities should be prioritized to provide the fundamental energy required to form an ideal fracture network.

Increasing the injection rate effectively improves the uniformity of fracture propagation, as clearly demonstrated by all three monitoring parameters. The fluid pressure curves (Fig. 8) show that as the injection rate increases from 0.02 to 0.06 m3/s, the pressure in all clusters rises synchronously, and the inter-cluster pressure difference narrows significantly. This indicates that the additional energy provided by a higher injection rate can partially suppress the stress shadow effect, enabling more balanced propagation conditions across clusters. The flow rate curves (Fig. 9) also exhibit notable changes, with all clusters converging toward a similar trend. The initial flow rate of the middle cluster increases significantly, and the flow rates among clusters tend to achieve dynamic equilibrium. At an injection rate of 0.06 m3/s, the instantaneous flow-rate difference between clusters is controlled to within 20%, demonstrating improved uniformity. The fluid intake volume curves (Fig. 10) further confirm this effect. As the injection rate increases, the rate of change in fluid intake volume becomes more consistent across clusters, leading to a significant improvement in the final uniformity of fluid distribution. These data confirm that a higher injection rate, by enhancing fluid kinetic energy and promoting balanced fluid intake among clusters, is an effective means to achieve uniform reservoir stimulation.

Figure 8: Effect of injection rate on fluid pressure.

Figure 9: Effect of injection rate on flow rate.

Figure 10: Effect of injection rate on fluid intake volume.

3.3 Effect of Horizontal Stress Difference

Analysis of the three-dimensional fracture geometries under horizontal stress differences of 5, 3, and 1 Mpa (Fig. 11) indicates that a reduction in the horizontal stress difference significantly alters the fracture propagation trajectory and extension efficiency [15], with the fracture width also progressively increasing from 7.35 to 11.3 mm.

Figure 11: Effect of horizontal stress difference on fracture geometry.

Under a high horizontal stress difference of 5 MPa, fractures exhibit relatively low propagation resistance, facilitating the development of longer, more freely turning geometries. In contrast, when the stress difference decreases to 1 MPa, fracture propagation becomes significantly constrained, resulting in more complex and shorter geometries with extension directions strictly controlled by the stress field orientation. The underlying mechanism for this phenomenon lies in the stress field’s control over the stress intensity factor at the fracture tip. According to linear elastic fracture mechanics, fracture initiation and propagation depend on the degree of stress concentration at the tip. A reduced horizontal stress difference elevates the stress threshold that fractures must overcome for propagation, thereby diminishing the effective driving force at the fracture tip. Under such conditions, fractures tend to propagate along the direction of the maximum principal stress, following the path of least energy consumption, which suppresses the generation of branch fractures. Concurrently, a lower stress difference increases rock mechanical stiffness, making the stress shadow effect more pronounced and further restricting the balanced development of multiple fractures.

In summary, the horizontal stress difference establishes the energy threshold and directional constraints for fracture propagation, forming the mechanical foundation of fracture geometry. This mechanism fundamentally differs from the engineering parameters discussed earlier. While parameters such as cluster spacing and injection rate are adjusted within a given in-situ stress background, the in-situ stress field itself determines the effectiveness of these optimization measures [23]. Therefore, accurate characterization of the in-situ stress field is a prerequisite for achieving optimized fracturing design.

Variations in horizontal stress difference primarily affect the baseline pressure and propagation difficulty, exerting a relatively indirect influence on fracture competition. The fluid pressure curves (Fig. 12) indicate that as the stress difference decreases from 5 to 1 MPa, the initiation pressure of each cluster increases linearly by approximately 4 MPa, with a synchronous rise in steady-state pressure. This confirms the fundamental control of the stress field on propagation resistance. The flow rate curves (Fig. 13) show a synchronous decrease in flow rate across all clusters, while the distribution ratio remains stable. This suggests that although changes in stress difference increase the difficulty of propagation, they do not alter the fundamental competitive dynamics. The fluid intake volume curves (Fig. 14) reveal that, under low horizontal stress-difference conditions, the overall stimulated volume of each cluster shows slight variation in scale. At the same time, the uniformity index remains essentially unchanged. This phenomenon indicates that the stress difference, as a background field parameter, primarily controls the absolute scale of stimulation. In contrast, the relative competitive relationship between fractures is more significantly influenced by other engineering parameters. In practical engineering applications, adjustment strategies for other parameters must be tailored to the prevailing stress-difference conditions.

Figure 12: Effect of horizontal stress difference on fluid pressure.

Figure 13: Effect of horizontal stress difference on flow rate.

Figure 14: Effect of horizontal stress difference on fluid intake volume.

4  Conclusion

This study employs a discrete lattice method-based numerical simulation to systematically investigate the influence of cluster spacing, injection rate, and horizontal stress difference on fracture propagation by monitoring the dynamic responses of pressure, flow rate, and fluid intake volume in each cluster during multi-cluster fracturing. The main conclusions are as follows:

(1)   The stress shadow effect, induced by the propagation of leading fractures, significantly alters the distribution of the in-situ stress field around the wellbore. This is the fundamental cause of non-uniform fracture propagation and is directly reflected in the dynamic evolution of inter-cluster pressure differences.

(2)   Cluster spacing primarily governs the intensity of the stress shadow effect through a geometric mechanism. When the spacing is reduced from 4 to 2 m, Inter-fracture competition is most intense, the maximum inter-cluster pressure difference increases by 7.32 MPa, and the fluid intake volume decreases substantially. This leads to strong suppression of the middle-cluster fracture and results in severely non-uniform propagation.

(3)   Appropriately increasing the injection rate can significantly suppress the inter-cluster pressure difference, while also increasing the fracture width. The additional fluid kinetic energy provided by a higher injection rate helps overcome the rock stress barriers created by the stress shadow effect, making it the most direct and effective controllable means for optimizing uniform fracture propagation in field operations.

(4)   The horizontal stress difference linearly determines the overall pressure level required for fracture initiation and propagation, serving as the foundational background field for fracturing design. Its magnitude modulates the absolute intensity of other engineering parameters.

(5)   The present study primarily investigates the influence of individual parameter variations on fracture morphology in hydraulic fracturing simulations. To gain a more comprehensive understanding of fracture propagation mechanisms and enhance the generalizability of the findings, future research could incorporate dimensionless parameter analysis. This approach would allow for a more fundamental interpretation of the underlying physical processes and broaden the applicability of the results across different geological and engineering conditions.

Acknowledgement: The authors gratefully acknowledge the financial support provided by the Open Fund of the Hubei Key Laboratory of Oil and Gas Drilling and Production Engineering (Yangtze University), YQZC202406. Sincere gratitude is also extended to the Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, for granting the software license used in this research.

Funding Statement: The authors gratefully acknowledge the financial support provided by the Open Fund of the Hubei Key Laboratory of Oil and Gas Drilling and Production Engineering (Yangtze University), YQZC202406.

Author Contributions: The authors confirm their contributions to the paper as follows: conceptualization: Pujin Wang, Wenwei Zhao; methodology: Pujin Wang, Wenwei Zhao, Liangping Yi; software: Pujin Wang, Guofa Ji, Wenwei Zhao, Liangping Yi; validation: Pujin Wang, Guofa Ji, Wenwei Zhao, Liangping Yi; resources: Wenwei Zhao; data curation: Pujin Wang, Wenwei Zhao; writing—original draft preparation: Pujin Wang; writing—review and editing: Pujin Wang, Guofa Ji, Wenwei Zhao, Liangping Yi; supervision: Guofa Ji, Liangping Yi; funding acquisition: Guofa Ji. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data can be provided upon request from the authors. The data supporting the findings of this study are available from the corresponding author, Guofa Ji, upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Yin Z, Zhang Y, Zhang F, Wang X, Zhang L. Dynamic evolution mechanisms of induced stresses in hydraulically fractured wells: incorporating real gas characteristics. Int J Rock Mech Min Sci. 2025;195(1):106298. doi:10.1016/j.ijrmms.2025.106298. [Google Scholar] [CrossRef]

2. Wang X, Zhang F, Zhu H, Hou B. Numerical simulation on fracture propagation of multilayer combined fracturing in coal-bearing strata: analyzing perforation design effects. Eng Fract Mech. 2026;332:111807. doi:10.1016/j.engfracmech.2025.111807. [Google Scholar] [CrossRef]

3. Li L, Zeng J, Chen W, Guo T, Luo Z, Zhang Y. Numerical simulation of multi-cluster fracture competitive propagation and proppant placement. J Shenzhen Univ Sci Eng. 2025;42(2):183–91. doi:10.3724/sp.j.1249.2025.02183. [Google Scholar] [CrossRef]

4. Yang C, Yang Z, Wang H, Jiang L, Yi L, Cheng Y, et al. Mechanisms of non-uniform propagation of hydraulic fractures: a comprehensive numerical investigation. Eng Anal Bound Elem. 2025;178(5):106307. doi:10.1016/J.ENGANABOUND.2025.106307. [Google Scholar] [CrossRef]

5. Damjanac B, Cundall P. Application of distinct element methods to simulation of hydraulic fracturing in naturally fractured reservoirs. Comput Geotech. 2016;71(1):283–94. doi:10.1016/j.compgeo.2015.06.007. [Google Scholar] [CrossRef]

6. Damjanac B, Detournay C, Cundall PA. Application of particle and lattice codes to simulation of hydraulic fracturing. Comput Part Mech. 2016;3(2):249–61. doi:10.1007/s40571-015-0085-0. [Google Scholar] [CrossRef]

7. Huang T, Zhu H, Liu Q, Dai X, Tang X, Yi X, et al. Study on multi-cluster fracture interlaced competition propagation model of hydraulic fracturing in heterogeneous reservoir. Geoenergy Sci Eng. 2025;244(4):213415. doi:10.1016/J.GEOEN.2024.213415. [Google Scholar] [CrossRef]

8. Huang L, Liu J, Zhang F, Fu H, Zhu H, Damjanac B. 3D lattice modeling of hydraulic fracture initiation and near-wellbore propagation for different perforation models. J Pet Sci Eng. 2020;191:107169. doi:10.1016/j.petrol.2020.107169. [Google Scholar] [CrossRef]

9. Guo X, Aibaibu A, Wu Y, Chen B, Zhou H, Zhu B, et al. Numerical modeling of hydraulic fracturing interference in multi-layer shale oil wells. Processes. 2024;12(11):2370. doi:10.3390/PR12112370. [Google Scholar] [CrossRef]

10. Zhao K, Stead D, Kang H, Damjanac B, Donati D, Gao F. Investigating the interaction of hydraulic fracture with pre-existing joints based on lattice spring modeling. Comput Geotech. 2020;122:103534. doi:10.1016/j.compgeo.2020.103534. [Google Scholar] [CrossRef]

11. Cong Z, Li Y, Tang J, Martyushev DA, Hubuqin, Yang F. Numerical simulation of hydraulic fracture height layer-through propagation based on three-dimensional lattice method. Eng Fract Mech. 2022;264:108331. doi:10.1016/j.engfracmech.2022.108331. [Google Scholar] [CrossRef]

12. Suo Y, Wei HX, Zhao JY, Wei YJ, Miao S, Zhao Y, et al. Simulation analysis of hydraulic fracture initiation and propagation mechanisms under mixed-mode and reservoir fracturing. Pet Sci. 2026;23(1):220–35. doi:10.1016/J.PETSCI.2025.09.025. [Google Scholar] [CrossRef]

13. Li J, Zhu Z, Hua W, Yang Y, Dong S, Guo T. Numerical investigation of fracture competitive propagation mechanisms for temporary plugging staged fracturing (TPSF) in shale gas reservoirs. Comput Geotech. 2025;180(5):107109. doi:10.1016/J.COMPGEO.2025.107109. [Google Scholar] [CrossRef]

14. Paullo Munoz LF, Mejia C, Rueda J, Roehl D. Iterative coupling of DDM and FEM for planar hydraulic fracturing under non-uniform confining stresses. Comput Geotech. 2025;186(16):107432. doi:10.1016/j.compgeo.2025.107432. [Google Scholar] [CrossRef]

15. Dontsov EV. An approximate solution for a penny-shaped hydraulic fracture that accounts for fracture toughness, fluid viscosity and leak-off. R Soc Open Sci. 2016;3(12):160737. doi:10.1098/rsos.160737. [Google Scholar] [PubMed] [CrossRef]

16. Wan X, Rasouli V, Damjanac B, Pu H. Lattice simulation of hydraulic fracture containment in the North Perth Basin. J Pet Sci Eng. 2020;188(5):106904. doi:10.1016/j.petrol.2020.106904. [Google Scholar] [CrossRef]

17. Liu YH, Zhang JT, Bai J, Zhang FS, Tang JZ. Numerical study of hydraulic fracturing in the sectorial well-factory considering well interference and stress shadowing. Pet Sci. 2023;20(6):3567–81. doi:10.1016/j.petsci.2023.05.020. [Google Scholar] [CrossRef]

18. Liu X, Qu Z, Guo T, Sun Y, Wang Z, Bakhshi E. Numerical simulation of non-planar fracture propagation in multi-cluster fracturing with natural fractures based on Lattice methods. Eng Fract Mech. 2019;220(3):106625. doi:10.1016/j.engfracmech.2019.106625. [Google Scholar] [CrossRef]

19. Chen B, Li S, Tang D. Numerical simulation study on hydraulic fracture propagation of multi-cluster fracturing of horizontal well in deep fractured coal seams. Eng Fract Mech. 2025;318:110983. doi:10.1016/j.engfracmech.2025.110983. [Google Scholar] [CrossRef]

20. Lu Y, Shelepov VV, Yang Z, Liu J, Hana J, Li X, et al. Research on initial pressure of fracture in the technique of staged fracturing of a horizontal well. J Appl Mech Tech Phys. 2018;59(6):1111–7. doi:10.1134/S0021894418060172. [Google Scholar] [CrossRef]

21. Zhao Y, Wang L, Xu W, Wang Q, Qiu S, Xiong F. Study on the initiation and propagation law of hydraulic fractures in tight reservoirs with cemented natural fractures for clean natural gas production. Int J Hydrogen Energy. 2026;220(5):153912. doi:10.1016/J.IJHYDENE.2026.153912. [Google Scholar] [CrossRef]

22. Zhang Y, Guo T, Chen M, Qu Z, Cao J, Yang X, et al. Numerical simulation study of fracture height growth considering the influence of bedding planes. Comput Geotech. 2024;168(8):106144. doi:10.1016/j.compgeo.2024.106144. [Google Scholar] [CrossRef]

23. Zhang G, Guo T, Chen M, Qu Z. Investigation of hydraulic fracture initiation and propagation behavior in glutenite reservoir using an improved global cohesive zone model. Comput Geotech. 2025;186(4):107421. doi:10.1016/j.compgeo.2025.107421. [Google Scholar] [CrossRef]

24. Zhang B, Gamage RP, Kong L, Zhang C. Cross-scale proppants transport in the shale hydraulic fracture network: a hybrid CFD-DEM investigation. Eng Geol. 2025;354:108160. doi:10.1016/j.enggeo.2025.108160. [Google Scholar] [CrossRef]

25. Zheng Y, Wei H, Zhang A, Niu Q. Exploring the influence of fluid viscosity on the hydraulic fracture propagation in composite coal seams. Eng Fract Mech. 2025;323(3):111203. doi:10.1016/j.engfracmech.2025.111203. [Google Scholar] [CrossRef]

26. Dontsov EV, Peirce AP. A multiscale Implicit Level Set Algorithm (ILSA) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-off asymptotics. Comput Meth Appl Mech Eng. 2017;313:53–84. doi:10.1016/j.cma.2016.09.017. [Google Scholar] [CrossRef]

27. Adachi JI, Detournay E, Peirce AP. Analysis of the classical pseudo-3D model for hydraulic fracture with equilibrium height growth across stress barriers. Int J Rock Mech Min Sci. 2010;47(4):625–39. doi:10.1016/j.ijrmms.2010.03.008. [Google Scholar] [CrossRef]

28. Hou J. Numerical simulation studies on multi-cluster fracture propagation in tight-oil horizontal wells. Pet Geol Recovery Effic. 2024;31(3):137–46. (In Chinese). doi:10.13673/j.pgre.202311006. [Google Scholar] [CrossRef]

29. Li JL, Zhong Y, Yan YH, Lu T, Yang XL, Li GF. A study on non-uniform propagation of multi-cluster hydraulic fractures in deep shale reservoirs based on a two-level flow distribution model. J Eng Geol. 2025;33(5):1697–711. (In Chinese). doi:10.13544/j.cnki.jeg.2025-0146. [Google Scholar] [CrossRef]

30. Yang CJ, Yang XD, Wang YB, Zhou SY, Song CP. Study on hydraulic fracture propagation in layered rock mass. Coal Geol China. 2025;37(8):35–40. (In Chinese). [Google Scholar]

31. Bao X, Xiao Y, Ding L, Tang Y, Ai R. Numerical simulation of fracture propagation in Sand conglomerate reservoir under 3D lattice element method. Results Eng. 2025;28(7):107365. doi:10.1016/j.rineng.2025.107365. [Google Scholar] [CrossRef]

32. Li S, Li Y, Cheng W. Numerical simulation on the propagation behaviour of hydraulic fractures in sandstone-shale interbeds. Processes. 2025;13(10):3318. doi:10.3390/pr13103318. [Google Scholar] [CrossRef]