Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Sci Rep . 2026 Mar 4;16:11984. doi: 10.1038/s41598-026-43148-w Search in PMC Search in PubMed View in NLM Catalog Add to search Investigation of the propagation behavior of hydraulic fractures and its influencing mechanisms in fractured reservoirs based on a hydromechanical coupling numerical model Yuyang Liu Yuyang Liu 1 Research Institute of Petroleum Exploration and Development (RIPED), PetroChina, Beijing, 100083 China 2 National Shale Gas Research & Development (Experiment) Center, LangFang, 065007 China 3 Key Laboratory of Coal-rock Gas, CNPC, Langfang, 065007 China Find articles by Yuyang Liu 1, 2, 3 , Xun Gong Xun Gong 4 School of Space and Earth Science, Peking University, Beijing, 100871 China 5 Institute of Energy, Peking University, Beijing, 100871 China Find articles by Xun Gong 4, 5, ✉ , Xinhua Ma Xinhua Ma 1 Research Institute of Petroleum Exploration and Development (RIPED), PetroChina, Beijing, 100083 China 2 National Shale Gas Research & Development (Experiment) Center, LangFang, 065007 China 3 Key Laboratory of Coal-rock Gas, CNPC, Langfang, 065007 China Find articles by Xinhua Ma 1, 2, 3 Author information Article notes Copyright and License information 1 Research Institute of Petroleum Exploration and Development (RIPED), PetroChina, Beijing, 100083 China 2 National Shale Gas Research & Development (Experiment) Center, LangFang, 065007 China 3 Key Laboratory of Coal-rock Gas, CNPC, Langfang, 065007 China 4 School of Space and Earth Science, Peking University, Beijing, 100871 China 5 Institute of Energy, Peking University, Beijing, 100871 China ✉ Corresponding author. Received 2025 Nov 26; Accepted 2026 Mar 2; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13068976 PMID: 41781489 Abstract Numerous discontinuities, such as bedding planes and natural fractures, occur in reservoir rocks and significantly influence the propagation behavior of hydraulic fractures during reservoir stimulation. To elucidate the mechanisms underlying the influence of natural fracture networks in reservoir rocks on the propagation of hydraulic fractures, a finite element–discrete fracture model is employed to establish a fluid–solid coupled finite element–discrete fracture model. This model is used to investigate the propagation behavior of hydraulic fractures in fractured reservoirs and their underlying mechanisms. The results indicate that when distant from natural fracture networks, hydraulic fractures typically propagate along the direction of the maximum principal stress. Upon approaching natural fracture networks, the propagation path of hydraulic fractures is altered, leading to localized deflection. The mechanical properties of the rock matrix versus those of natural fracture networks, in situ stress, fracturing fluid viscosity, and injection rate significantly influence the propagation of hydraulic fractures in fractured reservoirs. Increased mechanical disparity between the rock matrix and natural fractures promotes deflection along natural fracture networks, resulting in the formation of complex fracture networks. However, increased in situ stress, fracturing fluid viscosity, and injection rate facilitate direct penetration of natural fractures by hydraulic fractures, yielding the formation of simple, long, straight primary fractures. Furthermore, the propagation distance of hydraulic fractures along the direction of the maximum principal stress is positively correlated with the in situ stress, fracturing fluid viscosity, and injection rate. The findings of this study provide theoretical guidance for optimizing fracturing design. Keywords: Natural fractures, Reservoir stimulation, Hydromechanical coupling, Hydraulic fracturing, Fracturing fluid Subject terms: Energy science and technology, Engineering, Solid Earth sciences Introduction Hydraulic fracturing technology, as a key method for enhancing oil and gas recovery, has become the cornerstone of unconventional reservoir stimulation in the exploitation of unconventional hydrocarbon resources 1 – 5 . This technique involves injecting high-pressure fluid into the formation to fracture the reservoir rock and activate preexisting natural fractures, thereby creating a complex fracture network that significantly enhances the reservoir permeability and connectivity 1 – 3 , 6 , 7 . In naturally fractured reservoirs, the propagation of hydraulic fractures is influenced by numerous factors, including the in situ stress, rock mechanical properties, temperature, and operational parameters 2 , 3 , 8 – 10 . Therefore, an in-depth investigation into the propagation mechanisms of hydraulic fractures in naturally fractured reservoirs and their influencing factors is important for optimizing fracturing designs and increasing the stimulated reservoir volume (SRV). The propagation behavior of hydraulic fractures in naturally fractured reservoirs has been extensively investigated through laboratory experiments and numerical simulations. As a fundamental research approach, laboratory-scale hydraulic fracturing experiments provide relatively reliable results and have become a critical method for validating numerical simulation outcomes. Blanton 11 , 12 and Warpinski et al. 13 investigated the interaction mechanisms between hydraulic fractures and natural fractures/discontinuities through hydraulic fracturing experiments. These authors identified stress anisotropy and the approach angle (the intersection angle between hydraulic and natural fractures) as the dominant factors governing fracture interactions, subsequently establishing interaction criteria for hydraulic–natural fracture systems. Zhou et al. 14 conducted true triaxial hydraulic fracturing experiments to analyze interaction mechanisms, and they revealed that hydraulic fractures approaching natural fractures may exhibit complex behaviors, including arrest, diversion, or direct crossing. These behaviors are primarily controlled by the strength of natural fractures, in situ stress conditions, and approach angle, findings that have been widely accepted in the research community 15 – 18 . Furthermore, Liu et al. 19 demonstrated that the distance between natural fractures and wellbores significantly influences the geometry of hydraulic fracture. While laboratory experiments have substantially advanced our understanding of hydraulic–natural fracture interactions, several limitations persist: (1) the scale effects inherent in experimental samples may reduce the reliability of the results, and (2) the constrained dimensions of laboratory setups cannot fully replicate field conditions. These inherent constraints consequently restrict the direct application of laboratory findings to field operations. In recent years, with the advancement of computational technologies, numerical simulation methods such as the finite element method (FEM), extended finite element method (XFEM), boundary element method (BEM), discrete element method (DEM), phase-field method, and peridynamics have been increasingly applied in the oil and gas industry and are emerging as critical tools for studying the propagation of hydraulic fractures 20 – 28 . Lei et al. 29 employed a fluid‒solid coupled FEM model to simulate the propagation of hydraulic fractures at the field scale in naturally fractured reservoirs. Their findings highlighted that well-connected fracture systems facilitate the formation of complex fracture networks. Li et al. 30 developed a numerical model based on the FEM to investigate the interaction between hydraulic and natural fractures and reported that higher natural fracture density and lower stress anisotropy levels are favorable for the generation of complex networks. Zhang et al. 31 applied an XFEM-based numerical model to simulate the interaction between hydraulic fractures and shale bedding planes and demonstrated that the bedding orientation and fracture toughness significantly influence fracture propagation. Weng et al. 9 employed an unconventional fracture model to study the propagation of multiple hydraulic fractures in the presence of natural fractures and reported that reduced stress anisotropy or interfacial friction promotes the formation of complex networks. Gong et al. 32 adopted a phase-field model to analyze the interaction between hydraulic fractures and lithological interfaces and identified the in situ stress and rock mechanical properties as key factors controlling fracture behavior at discontinuities. Yan et al. 33 successfully captured the hydraulic fracturing process by establishing a numerical model through the coupled application of FEM and DEM, utilizing novel algorithms for solution. Wu et al. 34 simulated the interaction between hydraulic fractures and shale bedding planes during fracturing, employing a numerical model developed using FEM and DEM to account for the extensive bedding structures prevalent in shales. The study indicates that the slip characteristics of bedding planes critically determine whether hydraulic fractures penetrate these planes. Long et al. 35 developed a thermo-hydro-mechanical coupled numerical model by integrating FEM and DEM to analyse the interaction between hydraulic fractures and natural fractures during hot dry rock fracturing. Their research highlights that the horizontal principal stress difference is the primary factor governing fracture propagation. Each numerical method exhibits distinct advantages and limitations, rendering the selection of an appropriate approach a critical consideration in hydraulic fracturing research. Among these methods, the FEM is the earliest and most widely adopted technique, with its extensive applications and abundant case studies serving as key factors influencing its selection in this research. This study addresses the practical challenges associated with the propagation of hydraulic fracture behaviors and its influencing mechanisms in fractured reservoirs, which restrict field fracturing design and optimization. To address these challenges, a finite element–discrete fracture model for hydraulic fracturing is established. In this model, the fundamental theories of elasticity, damage mechanics, seepage mechanics, and fracture mechanics are integrated with a discrete fracture model using the FEM. By linking the FEM with a discrete fracture model, the approach aims to simulate the interaction between hydraulic fractures and natural fractures in naturally fractured reservoir rocks. A numerical model is then employed to simulate the interaction between hydraulic fractures and natural fractures within fractured reservoirs, and the influencing factors are analyzed. This study reveals the mechanisms that govern this interaction process, thereby providing guidance for predicting hydraulic fracture paths and optimizing fracturing design and implementation. The remainder of this article is organized as follows: Sect. 2 presents the mathematical model formulation, including governing equations and constitutive relationships; Sect. 3 provides a description of the numerical solution methodology, detailing discretization and solver implementation processes; Sect. 4 outlines the model validation process, thereby comparing the numerical results with experimental and analytical benchmarks; Sect. 5 offers an analysis of the propagation behavior of hydraulic fractures in naturally fractured reservoirs and its key influencing factors; Sect. 6 provides a comprehensive examination of the findings within the context of existing theories and field applications; and Sect. 7 outlines the major conclusions and their implications for hydraulic fracturing design optimization. Mathematical model The proposed hydromechanical coupling model is established on the basis of the following key assumptions 24 , 29 , 36 – 41 : (1) The reservoir rock is represented as a dual-porosity/dual-permeability system comprising a rock matrix and a fracture network; (2) compared with the intact rock matrix, natural fractures are treated as weak interfaces with inferior mechanical properties (e.g., elastic modulus and strength). Both the matrix and the fracture network obey Darcy’s law for fluid filtration; no chemical reactions occur between the fracturing fluid and rock formation in the simulation. (3) The temperature influences the properties of the fracturing fluid and the thermodynamic characteristics of reservoir rocks, thereby controlling the effectiveness of in situ fracturing simulations. Given the complex nature of its impacts on the fracturing process, temperature effects are neglected, and a constant ambient temperature is assumed throughout the entire fracturing operation. (4) The rock material follows small deformation theory, and no phase transformations occur throughout the entire simulation process. Rock deformation The strain‒displacement relationship during rock loading can be expressed as follows: 1 where denotes the strain component, dimensionless; and and denote the displacement components, m. The mechanical equilibrium equation for rock during quasistatic loading is as follows: 2 where is the total stress tensor, and is the body force per unit volume of rock. Assuming that the reservoir rock behaves as a linear elastic porous medium under the influence of the pore fluid pressure, the effective stress expression can be derived as follows 37 – 41 : 3 where G is the shear modulus, Pa; p is the pore pressure, Pa; is the Biot coefficient ( and denote the proportions of the rock matrix and natural fractures, respectively), which changes with increasing rock damage degree; is the dimensionless Poisson’s ratio; and is the Kronecker number. The substitution of Eq. (2) into Eqs. (1) and (3) yields the following: 4 Fluid flow In the hydraulic fracturing process, the reservoir rock is treated as a porous medium, and the mass conservation equation governing fluid flow can be expressed as follows 37 – 41 : 5 where ρ is the fluid density; is the storage coefficient of the porous medium ( and denote the storage coefficients of the matrix and natural fractures, respectively); is the fluid velocity vector; is the sink term ( and denote the source–sink terms for the rock matrix and natural fractures, respectively); and is the volumetric strain. Ignoring the effects of gravity, Darcy’s velocity can be calculated as follows: 6 where k is the permeability of the porous medium ( and denote the permeabilities of the matrix and natural fractures, respectively), m 2 ; and is the fluid dynamic viscosity, Pa*s. The water storage coefficient of porous media can be expressed as follows: 7 where is the fluid compression coefficient, and is the bulk modulus of the porous medium ( and denote the bulk moduli of the rock matrix and natural fractures, respectively). The porosity of the porous medium can be calculated as follows: 8 where is the residual porosity of the porous medium ( and denote the residual porosities of the matrix and natural fractures, respectively); is the initial porosity of the porous medium ( and denote the initial porosities of the matrix and natural fractures, respectively); is the pore–stress correlation coefficient of the porous medium ( and denote the pore–stress correlation coefficients of the rock matrix and natural fractures, respectively) 29 ; and is the average effective stress, which can be calculated as follows: 9 where , and denote the stresses along the three principal directions, Pa. The permeability of porous media can be calculated as follows: 10 where is the initial permeability of the porous medium ( and denote the initial permeabilities of the rock matrix and fractures, respectively); D is the damage factor; and is the permeability‒stress correlation coefficient of the porous medium ( and denote the permeability‒stress correlation coefficients of the rock matrix and natural fractures, respectively), with different values for the matrix and natural fracture. Damage evolution equation In this study, the fractured reservoir is treated as a brittle rock formation. Accordingly, the maximum tensile stress criterion and the Mohr‒Coulomb criterion are employed as rock failure criteria. First, it is determined whether is greater than or equal to 0 and whether is greater than or equal to 0. As long as one criterion is satisfied, rock rupture occurs ( and can be obtained by Eq. (11)) 37 – 41 . 11 where and denote the threshold functions for tensile and shear damage, respectively, Pa; and denote the maximum and minimum principal stresses, respectively, Pa; and denote the tensile and compressive strengths of the rock, respectively, Pa; and denotes the angle of internal friction ( and denote the angles of internal friction of the rock matrix and natural fractures, respectively), °. Rock damage during hydraulic fracturing is irreversible. Consequently, the mechanical properties of rock progressively decrease with increasing damage accumulation. The damage state can be quantified using damage variable D, where D = 0 indicates intact rock (undamaged state) and D = 1 indicates complete rock failure (fully damaged state). The decrease in rock elastic modulus during damage progression can be expressed as follows: 12 where denotes the initial modulus of elasticity of the rock. Under uniaxial tensile stress conditions, rock damage evolution can be derived from the maximum tensile stress criterion as follows: 13 where is the residual strength of the reservoir rock, Pa; denotes the initial damage threshold strain (elastic limit under tension); is the ultimate tensile strain, with = ; denotes the ultimate strain coefficient, dimensionless; and is the residual strength coefficient, dimensionless. Both and are critical parameters for characterizing the constitutive behavior of reservoir rocks. Under uniaxial tensile stress conditions, the following relationships can be established: 14 In Eq. (13), parameter can be defined as follows 42 : 15 where , , and denote the three principal strains. Moreover, 〈x〉 can be expressed as follows: 16 Similarly, the damage constitutive relationship for the reservoir rock under shear loading can be derived as follows: 17 where denotes the compressive strain corresponding to the uniaxial compressive strength of the rock, which can be calculated as follows: 18 Rock heterogeneity The reservoir rock comprises minerals such as quartz, feldspar, carbonate, dolomite, and pyrite. Different minerals exhibit distinct chemical and mechanical properties, resulting in notable rock heterogeneity. Studies have confirmed that rock heterogeneity significantly influences the propagation of hydraulic fractures. Therefore, in this study, the Weibull distribution is employed to characterize the heterogeneity in the mechanical properties of the reservoir rock, including the elastic modulus and strength. The distribution function is as follows: 19 where is a mechanical parameter; is a scale parameter of the rock sample cell, which is related to the average shear rupture strength of all the cells; and m is the heterogeneity coefficient; the lower the value is, the more heterogeneous the rock 29 . Numerical solution The hydraulic fracturing simulation process involves highly nonlinear characteristics in both the temporal and spatial domains for fluid flow and rock deformation, thus significantly increasing the computational complexity. Commercial software COMSOL Multiphysics, which relies on the FEM, demonstrates superior capabilities in simulating coupled multiphysics fields, including temperature, fluid flow, rock deformation, and chemical reactions. Furthermore, this software provides a series of solvers and external interfaces for MATLAB, which enables users to customize programming and achieve effective integration between MATLAB and COMSOL Multiphysics. On the basis of this framework, to address the differences in mechanical parameters between the rock matrix and natural fractures in fractured reservoirs, MATLAB is first employed to generate heterogeneous rock mechanical parameters, thereby ensuring that the fracture strength is less than that of the rock matrix. A custom rock damage criterion is subsequently defined under the model component subnode, which serves as the damage evaluation criterion for the subsequent calculations. In the model, the governing equations for rock deformation and fluid flow are first initialized using a steady-state solver, followed by the determination of a simultaneous solution with a transient solver. The specific settings for the transient solver include temporal discretization on the basis of an implicit solution method, a relative tolerance of 0.001, a time step size of 1, and a full coupling algorithm selected for the solution. The detailed computational procedure is as follows (Fig. 1 ): Heterogeneous mechanical parameters such as Young’s modulus, compressive strength, and tensile strength are generated within the range of model scales using MATLAB. The solid mechanics and Darcy’s law module in COMSOL Multiphysics is selected to establish a geometric model, input model parameters, and define model boundary conditions and initial conditions. Moreover, tensile and compressive damage criteria for the rock are defined under the component definition subnode, which are employed as damage criteria in the subsequent calculation process. The steady-state solution of the model is first used. In one time step, the steady-state solver is chosen to calculate rock deformation and fluid flow to initialize the model. A fully coupled algorithm is adopted to solve the model within the framework of the transient solver. The tensile and shear damage criteria are employed to assess the nodal units. The tensile stress is evaluated first, while the shear stress is assessed if the damage conditions are not satisfied. Thereafter, the parameters are processed according to the results. The time step is increased sequentially, and the above steps are repeated until the maximum time step is reached. Fig. 1. Open in a new tab Flowchart of the model solution process. Model validation In this section, the reliability of the model is validated. Specifically, the model is validated on the basis of hydraulic fracturing experimental results and simulation results. Comparison with experimental results First, the established hydromechanical coupling model was validated against physical experimental data. Guo et al. 43 conducted large-scale true triaxial hydraulic fracturing experiments of sandstone outcrops to investigate the influences of stress conditions and operational parameters on the propagation of hydraulic fractures. In this study, a 500 × 500 mm two-dimensional model is established in accordance with the experimental setup of Guo et al. 43 , with boundary conditions and operational parameters configured to match their physical experiment. As shown in Fig. 2 , the simulation results indicate that the hydraulic fracture initiates from the wellbore and propagates along the direction of the maximum horizontal principal stress, which agrees well with the experimental observations of Guo et al. 43 . However, minor differences exist between the experimental and numerical simulation results because of slight curvatures in the fracture path. These differences can be attributed to rock heterogeneity and imperfect parameter matching in the simulation process. Fig. 2. Open in a new tab Results of the physical modeling experiments and numerical simulations. (a1–a4) Experimental results of Guo et al. 43 . (b) Numerical simulation results. (c) Pressure variation with time. Second, by analyzing the pump pressure curve in the fracturing process (Fig. 2 c), it can be observed that as the fracturing fluid is continuously injected, the fluid pressure in the wellbore and its rate of increase gradually increase until the fracture pressure is reached, at which point rock fracturing occurs. Owing to high brittleness of the rock, the fracturing process is concluded within a relatively short period. Overall, the numerical simulation results highly agree with the experimental results. The existing discrepancies can be attributed to rock heterogeneity and the incomplete matching of the simulation parameters with those of the physical model experimental rock samples. Interaction of hydraulic fractures with natural fractures In the reservoir stimulation process using hydraulic fracturing technology, the propagation behavior of hydraulic fractures becomes complex when they approach natural fractures, bedding planes, or other discontinuities within the rock. Existing research suggests that when a hydraulic fracture encounters a natural fracture, it may exhibit complex behaviors, such as diverting along the natural fracture, directly crossing the natural fracture, or both crossing and diverting along the natural fracture (Fig. 3 ) 14 , 17 . The specific behavior depends on the geological conditions and engineering factors. On this basis, in this study, the interactions between hydraulic and natural fractures and their influencing factors were analyzed, and the propagation behaviors of hydraulic fractures under different conditions were compared with that reported in previous research to validate the reliability of the established model. The geometric parameters and boundary conditions of the model are shown in Fig. 4 a: notably, the model is square with a side length of 300 mm, featuring a central borehole with a diameter of 25 mm. Two natural fractures are distributed at the top and bottom, each with a horizontal projection length of 100 mm. The lower and left boundaries are roller-supported, whereas the upper and right boundaries are subjected to the maximum and minimum principal stresses, respectively. During fracturing, the fracturing fluid is injected into the wellbore at a constant rate. During the simulation, the interaction process between hydraulic and natural fractures was investigated by varying the difference in principal stress (the difference between the maximum and minimum principal stresses) and the approach angle (the intersection angle between hydraulic and natural fractures). The simulation results (Fig. 3 ) revealed that when the difference in principal stress is 4 MPa and the approach angle is 45° or 60°, the hydraulic fracture is diverted along the natural fracture upon encountering it. However, when the approach angle increased to 90° while maintaining the same difference in principal stress, the hydraulic fracture directly crossed the natural fracture (Fig. 3 c). A comparison of the results of previous studies with the simulation results (Fig. 4 b) indicated that when the approach angle is low, hydraulic fractures typically propagate along natural fractures. With increasing approach angle, the behavior transitioned from diversion along the natural fracture to directly crossing it. Additionally, a larger difference in principal stress promoted the crossing of the natural fracture by the hydraulic fracture. The findings suitably agree, with minor differences attributable to rock heterogeneity and certain simulation parameters that require manual calibration. Fig. 3. Open in a new tab Numerical simulation results for a difference in principal stress of 4 MPa. ( a ) Approach angle of 45°, turning. ( b ) Approach angle of 60°, turning. ( c ) Approach angle of 90°, direct crossing. Fig. 4. Open in a new tab Geometric boundary and simulation results. ( a ) Geometric boundary of the model. ( b ) Experimental and simulation results 11 – 13 , 37 , 44 . Propagation behavior of hydraulic fractures in fractured reservoirs and its influencing factors In this section, the propagation behavior of hydraulic fractures in fractured reservoirs is investigated. On the basis of the results of indoor physical simulation experiments involving hydraulic fracturing, the geometric shape and boundary conditions of the numerical model were established. The model is square with a side length of 300 mm, and a borehole with a diameter of 25 mm is located at the center of the model. During fracturing, the fracturing fluid is injected primarily along the wellbore at a constant rate. The left and lower boundaries of the model are defined as roller-supported boundaries, and the upper boundary is subjected to a vertical stress of 18 MPa. Moreover, the right boundary is subjected to a minimum horizontal stress of 10 MPa (Fig. 5 a). Within the model, 100 line segments are assigned to simulate natural fractures. Among them, 50 natural fractures exhibit a dip angle of 45°, and the remaining 50 natural fractures exhibit a dip angle of 135°, rendering them orthogonal to each other. The minimum length of these natural fractures is 50 mm, and the maximum length is 200 mm, with their distribution following a power-law distribution. Fig. 5. Open in a new tab Model geometry boundary and mesh analysis results. ( a ) Model geometry boundary. ( b ) Mesh division. ( c ) Relationship between the grid cell size and mass. Next, the established model is subjected to mesh generation. In this study, the entire model domain is discretized using free triangular mesh elements (Fig. 5 ). The mesh cell sizes are 0.0008, 0.0009, 0.001, 0.002, and 0.003 m, yielding corresponding average mesh qualities of 0.9123, 0.9127, 0.9132, 0.8863, and 0.88, respectively. The closer the average mesh cell mass is to 1, the higher the mesh quality (Fig. 5 c). The analysis indicated that the highest mesh quality is achieved with a cell size of 0.001 m. Consequently, the model domain mesh cell size was set to 0.001 m. Specifically, the mesh element size for natural fractures and the borehole boundary was set to 0.001 m. In the remaining domain, the minimum mesh element size was set to 0.001 m, the maximum mesh element size was also set to 0.003 m, and the element growth rate was set to 1.05 (Fig. 5 b). This configuration ensures that the simulation results are minimally affected by the mesh, thereby increasing the model simulation accuracy Table 1 . Table 1. Simulation parameters[ 36 , 37 ]. Parameter Magnitude Parameter Magnitude Rock matrix Natural fracture Elastic modulus (E m ) 30 GPa Elastic modulus (E f ) 10 GPa Poisson’s ratio (ν m ) 0.2 Poisson’s ratio (ν f ) 0.23 Density (ρ m ) 2700 kg/m 3 Density (ρ b ) 2000 kg/m 3 Compressive strength (f Cm ) 120 MPa Compressive strength (f Cb ) 30 MPa Tensile strength (f Tm ) 8 MPa Tensile strength (f Tb ) 2 MPa Permeability (k m ) 1 × 10 − 18 m 2 Permeability (k f ) 1 × 10 − 15 m 2 Initial porosity (ε m0 ) 0.03 Initial porosity (ε f0 ) 0.06 Residual porosity (ε mr ) 0.005 Residual porosity (ε fr ) 0.003 Internal friction angle (θ m ) 30° Internal friction angle (θ f ) 15° Fluid properties Viscosity (µ) 0.05 Pa·s Density (ρ) 1000 kg/m 3 Heat capacity (C pf ) 4200 J/(kg*K) Fluid temperature (T f ) 293.15 K Other Vertical stress (σ V ) 18 MPa Horizontal stress (σ h ) 10 MPa Open in a new tab Propagation of hydraulic fractures in fractured reservoirs The propagation behavior of hydraulic fractures in fractured reservoirs was subsequently simulated. As shown in Fig. 6 , the fracturing fluid is injected along the wellbore and gradually extends along the direction of the maximum principal stress. With increasing number of time steps, the hydraulic fracture approaches natural fractures and is subsequently diverted, leading to complex propagation behavior. In other words, the natural fractures present within the reservoir rock locally alter the propagation path of the hydraulic fracture. Additionally, as time progresses, an increasing amount of fracturing fluid flows along the direction of the maximum principal stress, resulting in a gradual increase in the fracture length along this direction, thereby increasing the SRV. However, the presence of natural fracture networks causes changes in the propagation direction of the hydraulic fracture in certain regions. Discontinuous interfaces become connected, thereby significantly increasing the SRV. The simulation results also revealed that factors such as in situ stress conditions, rock mechanical properties, and operational parameters influence the process of hydraulic fracture propagation in fractured reservoirs. On the basis of these findings, the factors influencing hydraulic fracture behavior were analyzed in detail, and the underlying mechanisms were elucidated to provide guidance for fracture design and optimization. Fig. 6. Open in a new tab Propagation of hydraulic fractures in naturally fractured reservoirs. ( a ) and (d) Propagation at 20 s. ( b ) and (e) Propagation at 50 s. ( c ) and (f) Propagation at 80 s. Influence of the elastic modulus ratio (E r /E f ) of the rock matrix on natural fractures The mechanical properties of reservoir rocks are typically characterized by the elastic modulus, compressive strength, and tensile strength. Natural fractures, as weak planes within reservoir rocks, generally exhibit less favorable mechanical properties (e.g., lower elastic modulus and compressive strength) than the rock matrix does. On this basis, the influence of rock mechanical properties on the propagation behavior of hydraulic fractures in naturally fractured reservoirs was investigated by varying the elastic modulus ratio between the rock matrix and natural fractures (E r /E f ) while maintaining the other parameters constant. During the simulation, the Er/Ef ratio was sequentially set to 2, 3, 4, and 5 to examine the process of hydraulic fracture propagation under different combinations of rock matrix strength and natural fracture strength. The simulation results (Fig. 7 ) revealed that when the other parameters remained unchanged and the E r /E f ratio was set to 2, 3, 4, or 5, before the natural fractures were approached, the hydraulic fractures propagated along the direction of the maximum principal stress. Upon encountering natural fractures, the fracturing fluid was diverted, causing the hydraulic fractures to deflect along the orientation of the natural fractures. As the E r /E f ratio was increased from 2 to 5, the hydraulic fractures propagated more notably along natural fractures. This behavior could be attributed to the fact that while the strength of natural fractures remains constant, a higher rock matrix strength requires more energy to fracture. Since natural fractures function as weak planes, significantly less energy is needed to activate them than to fracture the rock matrix. Consequently, the hydraulic fractures preferentially propagate along natural fractures. However, the overall propagation direction remains governed by the maximum principal stress because of in situ stress effects. Thus, with increasing degree of mechanical contrast between the rock matrix and natural fractures, the hydraulic fractures increasingly follow the natural fractures while maintaining a dominant propagation direction aligned with the maximum principal stress. Additionally, this tendency for hydraulic fractures to propagate along natural fractures results in larger fracture lengths, thereby increasing the SRV. Fig. 7. Open in a new tab Propagation behavior of hydraulic fractures under different E r /E f ratios. ( a ) and (e) E r /E f =2. ( b ) and (f) E r /E f =3. ( c ) and (g) E r /E f =4. ( d ) and (h) E r /E f =5. Effect of the geostress conditions Rock reservoirs are generally located deep underground, and the in situ stress state of rock reservoirs significantly influences hydraulic fracture propagation. Currently, the difference in principal stress (the difference between the maximum and minimum principal stresses) is typically employed to characterize variations in the in situ stress state. By analyzing the propagation behavior of hydraulic fractures under different principal stresses, the mechanism underlying the influence of the stress state on hydraulic fracture propagation can be revealed. On this basis, the propagation behavior of hydraulic fractures in naturally fractured reservoirs under different stress states was investigated by varying the maximum and minimum principal stresses while maintaining the other parameters constant. The difference in principal stress was sequentially set to 2, 4, 6, and 8 MPa to explore the mechanisms underlying the influence of the stress state on hydraulic fracture propagation. The simulation results (Fig. 8 ) indicated that when the other parameters were held constant, as the difference in principal stress was increased from 2 to 8 MPa, the propagation distance of hydraulic fractures along the direction of the maximum principal stress also increased. Additionally, when the hydraulic fractures encountered natural fractures, their behavior transitioned from diverting along the natural fractures to directly crossing them (highlighted by red dashed circles in Fig. 8 ). Thus, with the other parameters unchanged, with increasing difference in principal stress, the behavior of hydraulic fractures in naturally fractured reservoirs transitions from propagating along natural fracture networks to directly crossing them, ultimately aligning with the direction of the maximum principal stress. Furthermore, a greater difference in principal stress results in larger hydraulic fracture propagation distances, fewer connected natural fractures, and a significantly reduced SRV. Therefore, the influence of the in situ stress must be comprehensively considered during fracturing design optimization. Fig. 8. Open in a new tab Expansion behavior of hydraulic fractures in bedding under varying differences in principal stress. ( a ) Principal stress difference of 2 MPa. ( b ) Principal stress difference of 4 MPa. ( c ) Principal stress difference of 6 MPa. ( d ) Principal stress difference of 8 MPa. Next, the fluid pressure variation over time at the red point marked in Fig. 8 was analyzed. The results are shown in Fig. 9 . Adopting the case of an in situ stress difference of 8 MPa as an example, the hydraulic fracturing process can be divided into three stages, i.e., AB, BC, and CD, on the basis of the injection pressure curve. Stage AB (fluid injection stage): At this stage, as the fracturing fluid is injected into the wellbore, the wellbore pressure continuously increases until point B is reached, at which point the rock is fractured. The fluid pressure at point B corresponds to the fracture pressure of the rock. Stage BC (fluid pressure decrease stage): At this stage, as the fracturing fluid flows through the rock, the accumulated energy in the wellbore rapidly decreases until point C is reached, where the energy is relatively low. Stage CD (stable fracture propagation stage): At this stage, with the continuous injection of the fracturing fluid, the hydraulic fracture propagates at a relatively stable rate. The fluctuations in the fluid pressure curve can be attributed to rock heterogeneity and the presence of natural fractures. Additionally, as the in situ stress difference is increased from 2 to 8 MPa, the fracture pressure of the rock (the fluid pressure at point B) gradually decreases, as does the corresponding fracture time. This finding indicates a negative correlation between the in situ stress difference and the fracture pressure of the rock. Fig. 9. Open in a new tab Variations in the fluid pressure with time. Effect of the fracturing fluid viscosity Extensive field operations and laboratory experiments have demonstrated that the viscosity of the injected fluids significantly influences the fracturing performance. Accordingly, in this study, the propagation of hydraulic fractures in naturally fractured reservoirs was investigated by varying the fracturing fluid viscosity (25, 50, 75, and 100 mPa·s) while keeping the other parameters unchanged. As shown in Fig. 10 , when the hydraulic fractures were not located near natural fractures, they propagated along the direction of the maximum principal stress. As the hydraulic fractures approached the natural fractures, the fracturing fluid began to flow along the natural fracture orientation, thereby reactivating preexisting fractures. Consequently, the hydraulic fractures were subsequently extended along the direction of natural fractures. As the viscosity was increased from 25 to 100 mPa·s, the fracture behavior transitioned from diversion along natural fracture networks to directly crossing them, ultimately aligning with the maximum principal stress direction. The analysis indicated that a higher fracturing fluid viscosity corresponds to higher energy, enabling more effective rock fracturing during propagation. This in turn resulted in the hydraulic fractures preferentially extending along the direction of the maximum principal stress, yielding longer and straighter primary fractures, which reduced the degree of natural fracture connectivity. Additionally, with increasing fracturing fluid viscosity, the propagation length of the hydraulic fractures at the same time step increased accordingly. It can be inferred that higher-viscosity fluids can carry more proppant and possess more energy, causing the hydraulic fractures to directly cross natural fractures rather than being diverted along them, resulting in the formation of simpler and straighter fracture geometries. Consequently, this limits the SRV. In other words, the fracturing fluid viscosity indirectly influences the propagation behavior of hydraulic fractures in naturally fractured reservoirs by altering its proppant-carrying capacity. Fig. 10. Open in a new tab Hydraulic fracture propagation behavior under different fracturing fluid viscosities. ( a ) and (e) 25 mPa·s. ( b ) and (f) 50 mPa·s. ( c ) and (g) 75 mPa·s. ( d ) and (h) 100 mPa·s. Effect of the fracturing fluid injection rate In addition to the fracturing fluid viscosity, the fracturing fluid injection rate significantly influences the reservoir stimulation effect. On this basis, the propagation behavior of hydraulic fractures in naturally fractured reservoirs was investigated under different injection rates (0.00009, 0.0001, 0.0002, and 0.0003 kg/s) while maintaining the other parameters unchanged, aiming to elucidate the mechanism underlying the influence of the injection rate on hydraulic fracture propagation. The simulation results (Fig. 11 ) demonstrated that as the injection rate was increased from 0.00009 to 0.0003 kg/s, the propagation length of hydraulic fractures gradually increased. Additionally, when the hydraulic fractures encountered natural fractures, their behavior transitioned from diversion along natural fractures to directly crossing them. This transition could be attributed to the fact that higher injection rates provide higher fluid energy, enabling the fractures to propagate across larger distances within the same time step. Moreover, the increased energy causes the hydraulic fractures to preferentially cross natural fractures rather than being diverted along them. It can be inferred that higher injection rates promote the formation of simpler, longer, and straighter fractures, which limits the activation of natural fracture networks and consequently restricts the SRV. Furthermore, higher injection rates yield larger fracture propagation distances. In other words, the fracturing fluid injection rate is positively correlated with both the near-wellbore fracture extension range and the length of straight fractures. Fig. 11. Open in a new tab Hydraulic fracture propagation behavior under different fracturing fluid injection rates. ( a ) and (e) 0.00009 kg/s. ( b ) and (f) 0.0001 kg/s. ( c ) and (g) 0.0002 kg/s. ( d ) and (h) 0.0003 kg/s. Discussion As described in Sect. 5, the propagation behavior of hydraulic fractures in naturally fractured reservoirs is influenced by factors such as the mechanical properties of the reservoir rock, in situ stress, fracturing fluid viscosity, and injection rate. Via the use of the control variable method, the process of hydraulic fracture propagation was investigated under single-variable conditions. The results are summarized in Table 2 . An analysis of the fracture propagation behavior under varying values of the elastic modulus ratio, difference in principal stress, fracturing fluid viscosity, and injection rate revealed that as the difference in principal stress was increased from 2 to 8 MPa, as the fracturing fluid viscosity was increased from 25 to 100 mPa·s, and as the injection rate was increased from 0.00009 kg/s to 0.003 kg/s, the propagation of hydraulic fractures along the direction of the maximum principal stress was promoted. This resulted in the formation of simple, long, straight primary fractures, with the length of the generated fractures gradually increasing. When the hydraulic fractures approached natural fracture networks, increasing the difference in principal stress, fracturing fluid viscosity, and injection rate caused the hydraulic fractures to traverse natural fractures directly, ultimately propagating along the direction of the maximum principal stress. As the elastic modulus ratio was increased from 2 to 5, the hydraulic fractures were encouraged to turn along natural fractures, thereby activating more natural fractures and increasing the complexity of the resulting fracture network. Combining existing research, analysis indicates that an increase in principal stress differential reduces rock fracture pressure, while higher fracturing fluid viscosity and injection rates confer greater kinetic energy to the fluid. Consequently, when hydraulic fractures propagate towards natural fractures, they are more readily able to penetrate these natural fractures. Conversely, an increase in rock matrix and elastic modulus ratio renders the rock matrix more resistant to fracturing. When hydraulic fractures approach natural fractures, they tend to propagate along these pre-existing fractures, thereby increasing the complexity of the resulting fracture network. Furthermore, the fracturing process within reservoir rocks involves complex multiphysics interactions including temperature, rock deformation, fluid flow, and chemical reactions. Consequently, studying this phenomenon necessitates a detailed consideration of the influence of each individual factor. Secondly, when optimizing fracturing designs, it is essential to analyze the stress state, rock mechanical properties, and the degree of natural fracture development of the target reservoir in detail, among other geological conditions. On the basis of the geological conditions, appropriate engineering parameters such as fracturing techniques, fracturing fluids, and proppants can subsequently be selected. This approach ensures that the hydraulic fractures generated during fracturing are connected as extensively as possible with natural fractures, bedding planes, and other discontinuities within the reservoir rock. Ultimately, this process significantly increases the SRV, thereby maximizing the effective fracture zone. Table 2. Parametric conditions during fracturing operations. Influencing factors Level 1 2 3 4 Elastic modulus ratio 2 3 4 5 Principal stress difference (MPa) 2 4 6 8 Fluid viscosity (mPa·s) 25 50 75 100 Fluid injection rate (kg/s) 0.00009 0.0001 0.0002 0.0003 Open in a new tab On the basis of the numerical simulation results, for naturally fractured reservoir rock matrices with elastic modulus ratios ranging from 2 to 5 and principal stress differences ranging from 2 to 8 MPa, the optimal range for maximizing the SRV can be achieved by selecting fracturing fluid viscosities between 25 and 50 mPa·s and injection rates between 0.00009 and 0.0001 kg/s. Under the target reservoir conditions, reducing the fracturing fluid viscosity and injection rate increases the likelihood that the fracturing fluid will reach more discontinuities within the reservoir. This approach facilitates the formation of a complex fracture network, thereby enhancing the reservoir stimulation effect. Conclusion and outlook In this paper, a fluid‒solid coupled hydraulic fracturing numerical model is established experimentally, incorporating the fundamental theories of rock mechanics, elasticity theory, seepage mechanics, and damage mechanics. In this model, the FEM is integrated with a discrete fracture model. Through the numerical simulation of fracturing processes in reservoir rocks with well-developed natural fractures, the conclusions are as follows: In this paper, the results of prior hydraulic fracturing experiments and the numerical simulation results are first compared with those of the established finite element–discrete fracture model, thereby demonstrating the reliability of the model. During fracturing operations, hydraulic fractures primarily propagate along the direction of the maximum principal stress. While the presence of natural fractures locally alters the fracture propagation direction, they do not affect the overall trend of propagation along the direction of the maximum principal stress. When the hydraulic fractures approach natural fractures, the fracturing fluid is diverted along the natural fractures, thereby activating them and increasing the complexity of the resulting fracture network. The behavior of hydraulic fractures in fractured reservoirs depends on the ratio of the elastic modulus of the rock matrix to that of natural fractures, difference in principal stress, fracturing fluid viscosity, and injection rate. The distance traversed by the hydraulic fractures along the direction of the maximum principal stress increases with these parameters. A larger difference in principal stress, a higher fluid viscosity, and a higher injection rate help hydraulic fractures cross natural fractures. Moreover, the fracture pressure decreases with increasing difference in principal stress. A higher elastic modulus ratio causes the hydraulic fractures to turn and follow natural fractures. The numerical model established in this study primarily accounts for rock deformation and damage and fracturing fluid flow in the fracturing process under field conditions. However, phenomena such as heat exchange or chemical reactions between the fracturing fluid and the reservoir are not incorporated. Hence, future work should aim to integrate these aspects (heat exchange and chemical reactions) into the model, thereby enhancing the approximation of real-world fracturing conditions by the simulation results. This study provides a theoretical reference for field fracturing design and optimization. Acknowledgements This work has received funding from the Oil & Gas Major Project of China (2025ZD1401403-04) and the PetroChina’s Fundamental Prospective Project (No. 2024DJ8705; No. 2023ZZ08). We extend our sincere gratitude for this support. Author contributions Conceptualization, Yuyang Liu and Xun Gong; Formal Analysis, Yuyang Liu and Xun Gong; Resources, Xinhua Ma; Data Curation, Yuyang Liu and Xun Gong; Writing-Original Draft Preparation, Yuyang Liu; Writing-Review and Editing, Yuyang Liu and Xun Gong; Visualization, Xinhua Ma; Supervision, Xinhua Ma; Project Administration, Xinhua Ma. Funding This work has received funding from the Oil & Gas Major of China of China (2025ZD1401403-04) and the PetroChina’s Fundamental Prospective Project (No. 2024DJ8705; No. 2023ZZ08). We extend our sincere gratitude for this support. Data availability Some or all data, models, or code generated or used during the study are available from the corresponding author by request. Declarations Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References 1. Hubbert, M. K. & Willis, D. G. Mechanics of hydraulic fracturing. Trans. AIME 210 , 153–168. 10.2118/686-G (1957). [ Google Scholar ] 2. Adachi, J. et al. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 44 (5), 739–757. 10.1016/j.ijrmms.2006.11.006 (2007). [ Google Scholar ] 3. Vengosh, A. et al. A critical review of the risks to water resources from unconventional shale gas development and hydraulic fracturing in the United States. Environ. Sci. Technol. 48 (15), 8334–8348. 10.1021/es405118y (2014). [ DOI ] [ PubMed ] [ Google Scholar ] 4. Gong, X., Ma, X. & Liu, Y. Analysis of geological factors affecting propagation behavior of fracture during hydraulic fracturing shale formation. Geomech. Geophys. Geo-Energy Geo-Resour. 10 , 102. 10.1007/s40948-024-00819-0 (2024). [ Google Scholar ] 5. Miao, H. et al. Hydrocarbon source and relationship between hydrocarbon charging process and reservoir tight period of the Denglouku Formation tight sandstone gas reservoirs in the Xujiaweizi Fault Depression, Songliao Basin. Nat. Resour. Res. 33 , 1657–1684. 10.1007/s11053-024-10359-9 (2024). [ Google Scholar ] 6. Detournay, E. Mechanics of hydraulic fractures. Annu. Rev. Fluid Mech. 48 , 311–339. 10.1146/annurev-fluid-010814-014736 (2016). [ Google Scholar ] 7. Miao, H. et al. The shale gas migration capacity of the Qiongzhusi Formation: Implications for its enrichment model. Phys. Fluids 37 (5), 056603. 10.1063/5.0267742 (2025). [ Google Scholar ] 8. Dahi-Taleghani, A. & Olson, J. E. Numerical modeling of multistranded-hydraulic-fracture propagation: Accounting for the interaction between induced and natural fractures. SPE J. 16 (3), 575–581. 10.2118/124884-PA (2011). [ Google Scholar ] 9. Weng, X. et al. Modeling of hydraulic-fracture-network propagation in a naturally fractured formation. SPE Prod. Oper. 26 (4), 368–380. 10.2118/140253-PA (2011). [ Google Scholar ] 10. Miao, H. et al. Hydrocarbon generation potential and organic matter accumulation patterns in organic-rich shale during the Mesoproterozoic oxygenation event: evidence from the Xiamaling Formation shale. Geomech. Geophys. Geo-Energy Geo-Resour . 9 , 134. 10.1007/s40948-023-00638-9 (2023). [ Google Scholar ] 11. Blanton, T. L. An experimental study of interaction between hydraulically induced and pre-existing fractures. SPE/DOE Unconventional Gas Recovery Symposium, SPE-10847-MS. Pittsburgh, PA, USA. (1982). 10.2118/10847-MS 12. Blanton, T. L. & Louisville Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs. SPE Unconventional Resources Conference/Gas Technology Symposium, SPE-15261-MS. KY, USA. (1986). 10.2118/15261-MS 13. Warpinski, N. R. & Teufel, L. W. Influence of geologic discontinuities on hydraulic fracture propagation. J. Pet. Technol. 39 (2), 209–220. 10.2118/13224-PA (1987). [ Google Scholar ] 14. Zhou, J. et al. Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs. Int. J. Rock. Mech. Min. Sci. 45 (7), 1143–1152. 10.1016/j.ijrmms.2007.11.011 (2008). [ Google Scholar ] 15. Cheng, W. et al. A criterion for identifying hydraulic fractures crossing natural fractures in 3D space. Pet. Explor. Dev. 41 (3), 371–376. 10.1016/S1876-3804(14)60043-6 (2014). [ Google Scholar ] 16. Zou, J. et al. Complex hydraulic-fracture-network propagation in a naturally fractured reservoir. Comput. Geotech. 135 , 104165. 10.1016/j.compgeo.2021.104165 (2021). [ Google Scholar ] 17. Xiong, D. & Ma, X. Influence of natural fractures on hydraulic fracture propagation behaviour. Eng. Fract. Mech. 276 , 108932. 10.1016/j.engfracmech.2022.108932 (2022). [ Google Scholar ] 18. Hu, Y. et al. Investigation of coupled hydro-mechanical modelling of hydraulic fracture propagation and interaction with natural fractures. Int. J. Rock Mech. Min. Sci. 169 , 105418. 10.1016/j.ijrmms.2023.105418 (2023). [ Google Scholar ] 19. Liu, Y. et al. Influence of natural fractures on propagation of hydraulic fractures in tight reservoirs during hydraulic fracturing. Mar. Pet. Geol. 138 , 105505. 10.1016/j.marpetgeo.2022.105505 (2022). [ Google Scholar ] 20. Cundall, P. A. & Strack, O. A discrete numerical model for granular assemblies. Géotechnique 30 (3), 331–336. 10.1680/geot.1980.30.3.331 (2008). [ Google Scholar ] 21. Ouchi, H. et al. A peridynamics model for the propagation of hydraulic fractures in naturally fractured reservoirs. Spe J. 22 (3), 1082–1102. 10.2118/182594-PA (2017). [ Google Scholar ] 22. Lecampion, B., Bunger, A. & Zhang, X. Numerical methods for hydraulic fracture propagation: A review of recent trends. J. Nat. Gas Sci. Eng. 49 , 66–83. 10.1016/j.jngse.2017.10.012 (2018). [ Google Scholar ] 23. Zhou, S., Zhuang, X. & Rabczuk, T. A phase-field modeling approach of fracture propagation in poroelastic media. Eng. Geol. 240 , 189–203. 10.1016/j.enggeo.2018.04.015 (2018). [ Google Scholar ] 24. Ju, Y. et al. Numerical analysis of the effects of bedded interfaces on hydraulic fracture propagation in tight multilayered reservoirs considering hydro-mechanical coupling. J. Pet. Sci. Eng. 178 , 356–375. 10.1016/j.petrol.2019.03.042 (2019). [ Google Scholar ] 25. Zhou, S. & Zhuang, X. Phase field modeling of hydraulic fracture propagation in transversely isotropic poroelastic media. Acta Geotech. 15 , 2599–2618. 10.1007/s11440-020-00960-6 (2020). [ Google Scholar ] 26. Chen, B. et al. A review of hydraulic fracturing simulation. Arch. Comput. Methods Eng. 28 , 1–58. 10.1007/s11831-020-09468-4 (2021). [ Google Scholar ] 27. Qin, M. & Yang, D. Numerical investigation of hydraulic fracture height growth in layered rock based on peridynamics. Theor. Appl. Fract. Mech. 125 , 103885. 10.1016/j.tafmec.2023.103885 (2023). [ Google Scholar ] 28. Shentu, J. et al. Investigation of hydraulic fracture propagation in conglomerate rock using discrete element method and explainable machine learning framework. Acta Geotech. 19 , 3837–3862. 10.1007/s11440-024-02290-3 (2024). [ Google Scholar ] 29. Lei, Q., Doonechaly, N. G. & Tsang, C. F. Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks. Int. J. Rock. Mech. Min. Sci. 138 , 104598. 10.1016/j.ijrmms.2021.104598 (2021). [ Google Scholar ] 30. Li, Z. et al. Numerical investigation on the propagation behavior of hydraulic fractures in shale reservoir based on the DIP technique. J. Pet. Sci. Eng. 154 , 302–314. 10.1016/j.petrol.2017.04.026 (2017). [ Google Scholar ] 31. Zhang, J. N. et al. Hydraulic fracture propagation at weak interfaces between contrasting layers in shale using XFEM with energy-based criterion. J. Nat. Gas Sci. Eng. 101 , 104502. 10.1016/j.jngse.2022.104502 (2022). [ Google Scholar ] 32. Gong, X. et al. Simulation investigation of hydraulic fracture propagation patterns at lithological interfaces based on the phase-field method. Rock. Mech. Rock. Eng. 58 , 2803–2828. 10.1007/s00603-024-04258-x (2025). [ Google Scholar ] 33. Yan, C. et al. Combined finite-discrete element method for simulation of hydraulic fracturing. Rock. Mech. Rock. Eng. 49 (4), 1389–1410 (2016). [ Google Scholar ] 34. Wu, S. et al. Influence of slip and permeability of bedding interface on hydraulic fracturing: A numerical study using combined finite-discrete element method. Comput. Geotech. 148 , 104801 (2022). [ Google Scholar ] 35. Long, T. et al. Hydraulic fracturing of naturally fractured hot dry rock based on a coupled thermo-hydro-mechanical model. Geoenergy Sci Eng 2025: 214163. (2025). 36. Gong, X. et al. Effect of bedding structures on hydraulic fracture propagation behavior investigated using a coupled thermo-hydraulic-mechanical numerical model based on the phase-field cohesive zone method. Comput. Geotech. 186 , 107427 (2025). [ Google Scholar ] 37. Zhang, B. et al. Numerical simulation of fracture propagation and production performance in a fractured geothermal reservoir using a 2D FEM-based THMD coupling model. Energy 273 , 127175. 10.1016/j.energy.2023.127175 (2023). [ Google Scholar ] 38. Wang, F. & Kobina, F. The influence of geological factors and transmission fluids on the exploitation of reservoir geothermal resources: Factor discussion and mechanism analysis. Res. Sci. 1 (1), 3–18 (2025). [ Google Scholar ] 39. Wu, J. & Ansari, U. From CO2 sequestration to hydrogen storage: Further utilization of depleted gas reservoirs. Re Sci. 1 (1), 19–35 (2025). [ Google Scholar ] 40. Li, M., Liu, J. & Xia, Y. Risk prediction of gas hydrate formation in the wellbore and subsea gathering system of deep-water turbidite reservoirs: Case analysis from the South China Sea. Res. Sci. 1 (1), 52–72 (2025). [ Google Scholar ] 41. Zheng, P. et al. Formation mechanisms of hydraulic fracture network based on fracture interaction. Energy 243 , 123057 (2022). [ Google Scholar ] 42. Mazars, J. & Pijaudier-Cabot, G. Continuum damage theory—Application to concrete. J. Eng. Mech. 115 (2), 345–365. 10.1061/(ASCE)0733-9399 (1989). [ Google Scholar ] 43. Guo, T. et al. Physical simulation of hydraulic fracturing of large-sized tight sandstone outcrops. Spe J. 26 (1), 372–393. 10.2118/205003-PA (2021). [ Google Scholar ] 44. Gu, H. et al. Hydraulic fracture crossing natural fracture at nonorthogonal angles: a criterion and its validation. Spe Prod. Oper. 27 (1), 20–26. 10.2118/139984-PA (2012). [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Data Availability Statement Some or all data, models, or code generated or used during the study are available from the corresponding author by request. Articles from Scientific Reports are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (3.3 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top