Numerical Simulation of Fracture Propagation in Tight Formation Considering Natural Fractures Distributions

1  Introduction

China’s fundamental energy reserve pattern is characterized by abundant tight sandstone formations yet limited oil and gas reserves. Projections suggest that by 2050, tight sandstone is anticipated to account for no less than 50% of China’s primary energy consumption, underscoring its enduring dominance in the nation’s energy mix over the foreseeable future [1–3]. Currently, hydraulic fracturing technology is the most widely applied and effective technique for tight sandstone reservoir stimulation in China. Its primary objective is to enhance permeability by generating new fracture networks or connecting primary fractures in tight sandstones [4].

Comprehending the propagation characteristics of hydraulic fractures under varied conditions is paramount to the success of hydraulic fracturing operations. Scholars both domestically and internationally primarily investigate the propagation principles of these fractures under varied operational scenarios, leveraging physical experimentation, numerical simulations, and theoretical analyses. Fatah and Al-Yaseri [5] initially formulated the “criterion for tensile fracture induction via stress concentration along borehole walls,” postulating that hydraulic fractures initiated by hydraulic fracturing are oriented perpendicular to the direction of the minimum principal stress and propagate toward the vertical minimum principal stress axis. Uwakwe et al. [6] conducted a meticulous classification of the natural fracture network within tight sandstone formations, elucidating the morphological features of fractures observed during field fracturing operations. Dehghan [7] delved into the impact of intrinsic natural fractures on the propagation dynamics of hydraulic fractures in tight sandstones, unraveling the underlying mechanisms governing hydraulic fracture propagation in the presence of such fractures. Fan et al. and Li et al. [8,9] examined the role of cleats in initiating and directing the propagation of hydraulic fractures within tight sandstones. To mitigate the confounding effects of tight sandstone heterogeneity on hydraulic fracture propagation outcomes, researchers have resorted to employing homogeneous simulation materials or sophisticated numerical modeling techniques to conduct hydraulic fracturing experiments. Shrivastava et al. [10], leveraging numerical simulations, analyzed the impact of bedding plane thickness and spacing on hydraulic fracture length and connectivity. Tang and Yang [11] employed cement mortar-based tight sandstone samples to investigate the initiation and propagation characteristics of hydraulic fractures under varied true triaxial stress conditions. Edirisinghe and Perera [12], utilizing discrete element numerical simulation software, constructed a fluid-structure interaction (FSI) model to explore the effects of displacement, Poisson’s ratio, and stress on the evolution of hydraulic fractures during fracturing operations. Lastly, Li et al. [13] applied the displacement discontinuity method (DDM) to simulate hydraulic fracture propagation in layered and naturally fractured reservoirs, capturing the intricate interplay between hydraulic fractures, natural fractures, and bedding planes. Although numerous experiments and numerical simulations related to tight sandstone hydraulic fracturing have been conducted, relatively little research has focused on fracture propagation in tight sandstone with pre-existing fractures. It remains unclear how hydraulic fractures propagate, how they interact with pre-existing natural fractures (NFs), and how geological and engineering parameters influence fracture propagation. Among these robust fracture propagation models, the finite element method (FEM) implemented in ABAQUS does not require additional refinement for new fractures, demands fewer degrees of freedom, and exhibits higher computational efficiency. Additionally, the finite element method (FEM) facilitates the establishment of natural fracture models. By globally embedding zero-thickness elements within the finite element method (FEM), arbitrary fracture flow can be achieved, and the simulation results are more realistic.

To analyze hydraulic fracture propagation in tight sandstone with pre-existing fractures, numerical simulations were performed. Grasping the propagation characteristics of hydraulic fractures under diverse conditions is pivotal to the success of hydraulic fracturing endeavors. However, uncertainties persist regarding how hydraulic fractures propagate, their interplay with pre-existing natural fractures (NFs), and the impact of geological and engineering parameters on fracture propagation. Therefore, this study investigates hydraulic fracture propagation in tight sandstone formations with inherent fractures. Firstly, a fracture intersection model was established using the Cohesive module in ABAQUS, and validated via laboratory experiments. Afterwards, cohesive elements are globally embedded and integrated with Python to generate a conjugate joint distribution model, thereby establishing a dispersed mesh model. Then, based on the established model, a hydraulic fracturing simulation scheme is employed to conduct simulation research on fracture propagation. The findings can provide theoretical support for hydraulic fracture propagation in fractured tight sandstone reservoirs. Currently, hydraulic fracturing technology is the most prevalent and effective means for stimulating tight sandstone reservoirs in China. Mastering the propagation characteristics of hydraulic fractures under varied conditions is key to the success of hydraulic fracturing operations. However, uncertainty remains regarding how hydraulic fractures propagate, their interplay with pre-existing natural fractures (NFs), and the impact of geological and engineering parameters on fracture propagation. Therefore, four geological and engineering parameters were selected to explore their impact on the simulation of tight sandstone fracturing. This work provides theoretical support for simulating the fracturing of tight sandstone.

2  Mathematical Model

2.1 Seepage-Stress Coupling Model

Rock is a porous medium that includes both solid and fluid filled pore spaces. There is an instantaneous interaction between pore pressure and solid matrix strain. The linear pore elasticity equation is used to establish the stress-strain relationship as follows [14–16]:

σij=λεvδij+2Gεij−Cξδij(1)

p=Cεv−Mξ(2)

Among them: σij represents stress, MPa; λ and G are the Lame parameters of porous materials; C and M are the additional elastic moduli of the two-phase medium, MPa; εv is the volumetric strain; δij is the Kronecker function; ξ is the strain parameter of the fluid relative to the volume deformation of the solid; p is the pore pressure, MPa.

When the total stress of the rock is balanced, the component of stress εij on j should be 0, that is:

σij,j=0(3)

The rock pore fluid also satisfies the mass conservation equation:

∂ξ∂t+qi,j=0(4)

Among them, qi,j is the fluid flow rate, m3/s.

2.2 Cohesion Law

The Cohesive element model is composed of the viscous traction separation law that controls the process of fracture propagation and the fluid dynamics equations that describe the fluid flow in permeable rocks around fractures and hydraulic fractures [17,18]. As shown in Fig. 1, the Cohesive element model can be divided into cohesive zone and cohesive unbroken zone. In the viscous zone, a pre-defined element composed of viscous units is embedded in the model, and hydraulic fracturing grows along and within the pre-defined element. The viscous unbroken zone refers to the incomplete fracture zone where the internal surface tension is not zero during separation.

Figure 1: Schematic diagram of cohesive model

Fig. 2 shows the traction separation criterion for viscous units. Point A in the figure represents the initial damage location, which corresponds to the displacement δ0 and maximum traction force Tmax at the time of the initial damage, B represents the point of complete damage and failure, δf corresponds to the effective displacement at complete damage, and Gc represents the energy release rate caused by complete damage. As the fluid is injected, under the action of fluid pressure, the viscous traction force increases with the opening of the viscous surface. When the traction force reaches the cohesive strength Tmax, it reaches its cohesive strength δ0. When the separation strength exceeds δ0, irreversible linear damage occurs. When the separation strength reaches the critical value δf (at the fracture tip of the material), the viscous strength and cohesion disappear, and the fracture tip is completely fractured [19,20].

Figure 2: Damage diagram of cohesive model

The bilinear viscous traction separation rule was used in the model, as shown in Fig. 2. The bilinear viscous traction separation constitutive relationship can be derived from the following equation [21]:

T={k0δ0<δ≤δ0(1−D)kδδ0<δ<δf0δ≥δf(5)

Among than, k0 is the initial stiffness of the viscous interface, k=T0/δ0; δmax is the maximum separation value obtained during the load period; D is a damage variable that represents the overall degree of damage to the material.

D=δmf(δmmax−δm0)δmmax(δmf−δm0)(6)

Among then, δm0 represents the displacement of the element node at the beginning of the failure; δmf is D = 1 the displacement opening value of the unit node; δmmax is the maximum displacement at which the unit nodes open during the loading process.

The analysis in this article adopts the maximum nominal stress criterion:

f=max{|tn|tn0,tsts0,tttt0}(7)

Among then, Traction vector t={tn,ts,tt}, three components, tn is the traction force perpendicular to the possible fracture surface, and the other two are the tangential traction forces perpendicular to each other on the possible fracture surface; tn0, ts0, tt0 represents the peak value of nominal stress [22]. Once the cohesive elements (CE) fractures, damage can be evaluated based on fracture energy theory. The damage evolution of cohesive elements (CE) during fracture propagation is determined by the B-K fracture criterion, defined as:

Gnc+(Gsc−Gnc){GshGT}η=Gc(8)

Among then, Gn, Gs and Gt, respectively represent the work done by the traction force and its conjugate relative displacement in the normal, first shear direction, and second shear direction. These parameters can be obtained using traction separation curves. In order to more accurately characterize the material properties, the traction separation curve was obtained through experimental methods. There is a relationship between Gsh=Gsc+Gtc and GT=Gsh+Gnc. Gnc, Gsc and Gtc represent the fracture energy in the normal direction, first and second shear directions, respectively.

2.3 Flow Properties of Viscous Units

It is presupposed that the tangential fluid flowing parallel to the fracture surface behaves as an incompressible Newtonian fluid, adhering to the pressure transmission formula characteristic of such fluids. The tangential fluid dynamics within the fracture conform to both the principles of lubrication theory and continuity equations [23,24].

q=d312μ∇pi(9)

δdδt+∇⋅(−d312μ∇pi)+(qt+qb)=Qinj(10)

Among then, μ is the dynamic viscosity of fracturing water, cp; d is the width of the fracture opening, m; q is the tangential volumetric flow rate of fracturing water, m3/s; ∇pi is the fluid pressure gradient along the viscous region, Pa/m; qt and qb are the flow rate of fracturing water flowing into and out of the upper and lower surfaces of the viscous unit, m3/s; They respectively reflect the filtration of fluids from the surface of fractures to adjacent layers. For non-permeable rock layers, qt=qb=0.

In the viscous unit, the normal flow rate of fracturing water is [25,26]:

{qt=Ct(pi−pt)qb=Cb(pi−pb)(11)

Among then, pt and pb is the pore pressure in the adjacent porous elastic material at the top and bottom of the viscous unit surface, respectively, Pa; Ct and Cb is the corresponding fluid filtration coefficient, m3/(Pa·s); pi is the pressure on the middle surface of the unit, Pa.

3  Numerical Simulation Model

3.1 Numerical Model

This study is based on a two-dimensional plane model and establishes a 300 × 300 m fracturing area. The fracturing injection point is located for simulating the wellbore injection, and a 0.5 m long perforation is pre-set above and below the injection point [27]. Many natural fractured reservoirs have more than one set of natural fractures (i.e., a set of fractures with one dominant orientation), and reservoirs containing two sets of natural fractures with different orientations are common. Usually, the initial set of natural fractures is related to tectonic stress and develops parallel to the compression direction, while the second set of fractures is caused by viscoelastic effects or orthogonal or sub orthogonal loading. Therefore, we designed a conjugate natural fracture distribution model, as shown in Fig. 3. Two dimensional 2D models are usually simple and clear, with lower simulation costs, more suitable for crack propagation, and do not involve interlayer problems. The modeling of 3D models is very complex and suitable for studying crack propagation in problems such as layering.

Figure 3: Conceptual model diagram of hydraulic fracturing modeling

The fundamental segment of the model is discretized utilizing a free quadrilateral meshing technique, with a universal mesh size of 0.5 m. To mimic the irregular expansion of fractures, cohesive elements of zero thickness are seamlessly integrated into the global grid. Fig. 4 illustrates the nodal distribution at the intersections of these cohesive units, showcasing the shared nodes employed at the interfaces between viscous elements and tight sandstone rock elements, as well as among viscous elements themselves. Notably, at the injection point, both viscous elements and injection nodes are configured as shared entities. The model’s boundary conditions encompass zero displacement, a static water pressure boundary, and an initial geostress that mirrors the static water pressure, eliminating any geostress differential. The fracturing fluid is injected at a rate of 0.01 m3/s, and the overall simulation duration time last for 25 s. The presence of natural fractures increases the complexity of fractures after fracturing. The strength of its weak surface (natural fractures) is much lower than that of the reservoir rock. Two innovative numerical models were established with different parameters. By globally embedding zero thickness elements in the finite element method, arbitrary fracture flow can be achieved, and simulation results are more realistic.

Figure 4: Distribution diagram of Cohesive element unit nodes in the model (modified after [22])

The mechanical parameters of tight sandstone rock foundation are as follows: elastic modulus of 4.2 GPa, Poisson’s ratio of 0.26, permeability coefficient of 3.5 × 10−8 m/s, and porosity ratio of 0.1. Set two sets of viscous unit properties, representing the mechanical properties of ordinary tight sandstone and weak joint surfaces, respectively. The elastic modulus, Poisson’s ratio, etc. are the same, but the difference is the initial damage stress (i.e., strength). The initial damage stresses of ordinary rock cohesive units are 1.2, 0.8, and 0.8 MPa, respectively; The initial damage stresses of the cohesive unit corresponding to the joint face are 0.4, 0.25, and 0.25 MPa, respectively. The fracturing fluid is water, with a density of 9800 kg/m3 and a viscosity of 0.001 Pa·s. The detailed model parameter settings are shown in Table 1.

The interaction between hydraulic fractures (HFs) and natural fractures (NFs) represents a complex yet crucial research topic in fracturing technology. To demonstrate the complex fracture propagation, Table 2 shows the interaction between hydraulic fracture and natural fracture. This interplay significantly influences the formation of fracture networks and, consequently, the stimulation effectiveness in hydrocarbon reservoirs. As high-pressure fluids are injected to create hydraulic fractures, stress waves are generated within the surrounding rock, triggering the propagation of natural fractures. When encountering natural fractures, hydraulic fractures may penetrate, intercept, or deflect them. In unconventional hydrocarbon reservoirs, the presence of natural fractures substantially enhances the complexity of fracture networks, thereby boosting fracturing effectiveness. Numerical simulation methods are prevalent in analyzing the influence of various factors, including the approach angle between hydraulic and natural fractures, differential principal stress, and formation elastic properties, on their interactions. Experimental investigations have identified three primary modes of interaction between these fractures under different conditions: penetration, interception, and deflection. In scenarios of high horizontal differential principal stress, steep approach angles, and substantial interfacial friction, hydraulic fractures predominantly traverse natural fractures. Conversely, in environments characterized by low differential stress, shallow approach angles, and reduced friction, hydraulic fractures are more prone to being intercepted by natural fractures. Moreover, the angle and spacing of natural fractures significantly influence the propagation path and extent of hydraulic fractures. For example, when natural fractures are perpendicular to hydraulic fractures, the latter tend to traverse a majority of natural fractures, resulting in intricate fracture networks. Numerical simulations underscore the crucial role of fluid viscosity, injection duration, and fluid displacement rate in governing the propagation of hydraulic fractures. The interplay between hydraulic and natural fractures is vital in enhancing reservoir permeability and recovery efficiency. By meticulously crafting fracturing strategies, we can optimize the utilization of natural fracture networks, thereby facilitating efficient oil and gas extraction. In essence, the interaction between hydraulic and natural fractures represents a complex, multistage process. By integrating theoretical analysis, numerical simulations, and experimental validations, we can gain profound understanding of this process, furnishing scientific rationale and technical backing for practical engineering applications.

3.2 Numerical Simulation Scheme

The propagation of hydraulic fractures induced by hydraulic fracturing within tight sandstone reservoirs, particularly those with well-developed fracture networks, presents a highly intricate process. The joint fracture energy (JFE) inherent in these tight sandstone fractures plays a pivotal role in determining the extent and nature of fracture development. Furthermore, the morphology of these fractures is intricately intertwined with the elastic modulus (E) of the tight sandstones, the fracturing fluid injection rate (FFIR), and the viscosity of the fracturing fluid (FFV). Notably, the impact of Poisson’s ratio on fracture propagation, though present, is deemed relatively negligible and will not be the primary focus of this discussion. Regarding the influencing factors, the elastic modulus and the reservoir parameters themselves, such as the inherent properties of the tight sandstone, are fixed and not readily alterable. In contrast, the fracturing fluid injection rate and viscosity represent engineering parameters that offer a degree of flexibility and can be adjusted strategically based on specific reservoir conditions and operational objectives. By leveraging this understanding, engineers and geologists can meticulously tailor the fracturing process to optimize fracture growth, ensuring efficient utilization of the natural fracture network and maximizing the potential for oil and gas recovery. This approach necessitates a comprehensive assessment of the reservoir’s characteristics, coupled with the strategic adjustment of engineering parameters to achieve the desired fracture propagation patterns.

The specific simulation plan is shown in Table 3. It should be noted that when changing the FFIR, maintain the same amount of fracturing fluid injection, and the injection rates of fracturing fluid are 6.0, 8.0, 10.0, and 12.0 m3/min, respectively. The corresponding fracturing times are 25 s, 18.75 s, 15 s, and 12.5 s, respectively.

4  Simulation Results and Analysis

4.1 The Influence of Joint Fracture Energy (JFE) on Fracture Propagation

The fracture energy inherent within tight sandstone formations holds a profound influence on the expansion dynamics of hydraulic fractures within these materials. To delve deeper into this phenomenon, we will embark on a study that specifically examines the consequences of varying the magnitude of joint fracture energy on the propagation patterns of fractures, while maintaining a consistent set of key parameters. This includes fixing the fracturing simulation duration at 25 s, ensuring a horizontal stress difference (HSD) of zero, and adhering to the parameter values outlined in Table 3 for all other relevant factors. Then, we will increase the joint fracture energy JFE to 100, 150, 200, and 250, respectively. Fig. 5 illustrates the transformation in the geometric configuration of fracture propagation across varying levels of joint fracture energy (JFE). As the JFE escalates, a discernible trend emerges: the joint’s resistance to deformation intensifies, resulting in an augmentation of its strength. With all other conditions held steady, this enhanced strength leads to a reduction in both the overall length and width of the primary fracture. However, paradoxically, the increased joint strength prompts the fracturing fluid to disseminate more extensively along the preexisting fractures, thereby enhancing the intricacy of the fracture network. It is noteworthy that when JFE reaches a threshold of 200 J/m2, further increments do not significantly alter the qualitative patterns of fracture propagation. Additionally, it is imperative to consider that fractures situated remotely from the fracture tip may be susceptible to shear failure, a phenomenon that underscores the complexity and dynamic nature of fracture behavior in tight sandstone formations. The shear failure of these fractures can release acoustic energy into the entire reservoir and can be detected through microseismics during the production process [28]. Fig. 6 shows the Fracture width and length under different JFE conditions. It can see that the fracture width decrease with the increase of JFE. In addition, the fracture length also decrease with the increase of JFE. Overall, JFE has a significant influence on the hydraulic fractures propagation in tight sandstones, and the fracturing fracture paths in the near wellbore area are more complex than those in the far wellbore area. When JFE increased from 100 to 250 J/m2, the fracture width decreased by 50% from 0.42 to 0.21 mm, and the fracture length decreased by 20% from 135 to 108 m. JFE has a significant impact on fracture morphology.

Figure 5: Fracture propagation morphology under different JFE conditions, HSD = 5 MPa

Figure 6: Fracture width and length under different JFE conditions, HSD = 5 MPa

4.2 The Influence of Elastic Modulus (E) on Fracture Propagation

Previous mechanical tests have shown that elastic modulus is a key parameter in tight sandstone rock reservoirs. Establishing elastic model parameters for different tight sandstone rock reservoirs. The joint fracture energy is set to 100 J/m2, the fracturing simulation time is 25 s, the horizontal stress difference (HSD) is 0, and other parameters are shown in Table 3. The elastic modulus is set to 2, 3, 4, and 5 GPa, respectively. Fig. 7 shows the geometric shape of fracture propagation under different elastic model E conditions. The findings suggest that as the elastic modulus (E) increases, the primary fracture lengthens, while the number of smaller, bifurcated fractures diminishes. This aligns with experimental outcomes from relevant mechanical tests, wherein higher E values correlate with enhanced brittleness in the reservoir, facilitating the formation of longer fractures. Notably, at an E of 4 GPa, the fracture length achieves its maximum. However, as E further escalates, the number of branching fractures commences to augment. It is also worth mentioning that during fracturing operations, these slender, bifurcated fracture channels are prone to closure due to their narrow widths, thus contributing minimally to production. From a fracture morphology perspective, E exerts a pronounced effect on fracture propagation. Specifically, as E grows larger, hydraulic fractures tend to become longer and narrower, while the count of primary fractures remains constant, accompanied by an increase in the number of branch fractures. The elastic modulus (E) of tight sandstone and rock reservoirs significantly impacts the expansion of hydraulic fractures. Fig. 8 clearly depicts the inverse relationship between fracture width and E, as well as the negative correlation between fracture length and E. When E increased from 2 to 5 GPa, the fracture width decreased from 0.38 to 0.25 mm, a decrease of 34%; And the fracture length increased from 121 to 151 m, an increase of 25%. The simulation results are consistent with the results of indoor experiments. The results of rock mechanics fracture experiments show that the larger the elastic modulus, the greater the brittleness of the rock, and the easier it is to form long fractures [29].

Figure 7: Fracture propagation morphology under different E conditions, HSD = 5 MPa

Figure 8: Fracture width and length under different E conditions, HSD = 5 MPa

4.3 The Effect of Fracturing Fluid Injection Rate (FFIR) on Fracture Propagation

Maintain the same amount of fracturing fluid injection and change the fracturing fluid injection rate FFIR (corresponding to a change in fracturing time) to study the effect of FFIR on the expansion of tight sandstone hydraulic fracturing fractures. The fracturing simulation time is 25 s, the horizontal stress difference (HSD) is 0, and other parameters are shown in Table 3. The injection rates of fracturing fluid are 6.0, 8.0, 10.0, and 12.0 m3/min, respectively, with corresponding fracturing times of 25, 18.75, 15, and 12.5 s. Fig. 9 shows the geometric shape of fracture propagation under different FFIR conditions. The results show that, under the same amount of fracturing fluid injection, with the increase of fracturing fluid injection rate, the width and length of fractures both increase. Moreover, when the FFIR increases to 10.0 m3/min, the length of fractures is the longest and the number of branching fractures is the highest. The fracturing effect is best under this model. When the FFIR continues to increase, due to the occurrence of perforation, the fracturing fluid leaks out, resulting in a decrease in the complexity of the final fracture. Fig. 10 illustrates the fracture width and length under different FFIR conditions. Both the fracture width and fracture length can be seen for the increase of FFIR. Within the range set by the model, while maintaining the same amount of fracturing fluid injection, the complexity of fractures increases with the increase of FFIR, and the range of fracturing control increases. However, for fracturing construction, an increase in displacement requires higher requirements for wellhead equipment, which is also a problem that needs to be considered. When FFIR increased from 6 to 12 m3/min, the fracture width increased from 0.26 to 0.42 mm, an increase of 62%; The fracture length increased from 116 to 153 m, an increase of 32%. The increase in FFIR has a significant impact on the fracture morphology and is beneficial for fracturing operations. However, the increase in FFIR requires too much equipment and results in more leakage.

Figure 9: Fracture propagation morphology under different FFIR conditions, HSD = 5 MPa

Figure 10: Fracture width and length under different FFIR conditions, HSD = 5 MPa

4.4 The Effect of Fracturing Fluid Viscosity (FFV) on Fracture Propagation

The viscosity of the fracturing fluid serves as a pivotal parameter in hydraulic fracturing operations. To isolate its influence, we maintain all other variables constant and delve into how varying the fracturing fluid viscosity (FFV) impacts the expansion of hydraulic fractures in tight sandstone formations. Our simulation spans 25 s, assumes a horizontal stress difference (HSD) of zero, and adheres to the parameters outlined in Table 3. We set the FFV at 1, 10, 50, and 100 cp, respectively, to observe its effects. Fig. 11 visually represents the geometric configurations of fracture propagation under different FFV scenarios. The data reveals that as FFV escalates, the fracture length and complexity remain relatively stable, implying that for tight sandstones with pre-existing fractures, augmenting the viscosity of the fracturing fluid has a limited impact. Fig. 12 further illustrates the relationship between fracture width and length against varying FFV. Notably, while an increase in FFV leads to a narrowing of fracture length, it concurrently enhances fracture width. Conventionally, boosting fracturing fluid viscosity is considered advantageous for fracture propagation and extension. However, in practical scenarios featuring densely packed formation fractures, this effect is negligible within the FFV range of 1–100 cp. Overall, the influence of fracturing fluid viscosity is minimal. Therefore, during on-site construction, when using low viscosity fracturing fluid, the influence of fracturing fluid viscosity can be appropriately ignored. When FFV increased from 1 to 100 cp, the fracture width increased from 0.24 to 0.44 mm, an increase of 83%; The fracture length decreased from 131 to 122 m, a decrease of 7%. An increase in FFV is beneficial for fracture width, but the length will decrease. Research suggests that in general, increasing the viscosity of fracturing fluid is beneficial for providing crack width, but the length of crack propagation is significantly reduced. Low viscosity fracturing fluid has strong diffusion ability and is more likely to penetrate into microcracks and pores, making it easier to form complex fracture networks. Therefore, during on-site construction, when using low viscosity fracturing fluid, the influence of fracturing fluid viscosity can be appropriately ignored to achieve sand carrying effect.

Figure 11: Fracture propagation morphology under different FFV conditions, HSD = 5 MPa

Figure 12: Fracture width and length under different FFV conditions, HSD = 5 MPa

5  Field Applications

To valid the numerical simulations, a case study from the Shunbei oilfield is selected. The target formation of sandstone matrix is tight, with lithology mainly composed of sandstone—predominantly lithic sandstone and feldspathic lithic sandstone, and a small amount of quartz sandstone. The lithic content is high (32%–41%), mainly consisting of volcanic rock fragments (andesite and basalt fragments) and metamorphic rock fragments (gneiss and phyllite fragments); the feldspar content ranges from 17% to 21.8% (dominated by potassium feldspar), and the quartz content is 28%–47%, indicating an overall low compositional maturity. The reservoir is distributed in a “strip-like” pattern around strike-slip fault zones. The fault fracture zones and derived fractures generated by multi-stage fault activities serve as the core reservoir spaces of the tight sandstone reservoir. For stimulation, this well adopted the horizontal well staged fracturing technology. The pumping rate was 10 m3/min, and a total of 8 fracturing stages were completed along the horizontal wellbore. Each stage was perforated with 2–3 clusters, with a cluster spacing of 30–60 m. During the proppant injection process, parameters such as proppant particle size and proppant-liquid ratio were considered, and a pulsed proppant injection technique was employed to ensure effective fracture conductivity by the proppant. An overall volume stimulation mode with a proppant intensity of 2.5–3.5 m3/m was adopted, achieving favorable results. Numerical simulation results confirm that the generated fractures are highly consistent with those detected by wide-area electromagnetic monitoring. Fig. 13 shows numerical simulation for one stage of the horizontal well in the target formation.

Figure 13: Numerical simulation for one stage of the horizontal well in the target formation

The fractures on both sides extend more fully, with the left-side fracture reaching a length of 147 m and the right-side fracture 138 m. However, the middle fracture is affected by stress shadowing, resulting in limited overall extension, with a total length of approximately 125 m. The total error is about 12.11%. Further findings from the simulation reveal that natural fractures at local locations exhibit differences in activation. Specifically, natural fractures perpendicular to the fault strike open earlier than those parallel to the fault strike, and fracturing fluid leakage mainly occurs along the direction of the minimum principal stress. In terms of their respective roles, the wide-area electromagnetic method focuses more on the actual monitoring of the real conditions of the fractured network induced by fracturing, while numerical simulation emphasizes the analysis and prediction of the formation and propagation of the fractured network from a theoretical and modeling perspective. If the results of the two methods are compared, it may be found that the numerical simulation shows a certain degree of consistency with the actual monitoring results of the wide-area electromagnetic method in terms of predicting the propagation direction and morphology of fractures. These findings verify the reliability and accuracy of the model and the error mainly comes from the 2D (two dimensional) assumptions.

In real reservoirs, heterogeneities such as bedding planes, natural fractures, and stress anisotropy cause hydraulic fractures to frequently undergo diversion, branching, or curved propagation. However, 2D models force fractures to extend within a fixed plane, which can significantly underestimate or overestimate the actual extension range and morphology of fractures. Furthermore, 2D models are generally categorized into “plane strain models” (assuming a fixed fracture height, with only length and width calculated) and “plane stress models” (assuming a fixed fracture width). In actual fracturing operations, however, fracture height dynamically changes under the influence of factors like the stress difference between the upper and lower barrier layers, reservoir thickness, and fluid loss. Additionally, the “thickness” of fractures (i.e., fracture width) exhibits non-uniform distribution along the length and height in 3D space (e.g., wider fracture width near the wellbore and narrower width at the far end). In field applications, the assumption of fixed parameters in 2D models leads to deviations in the calculation of in-fracture fluid pressure distribution and net pressure. This, in turn, affects the prediction accuracy of fracture propagation speed and final scale. Particularly, 2D models struggle to handle well-fracture coupled flow, making it impossible to provide effective guidance for the optimization of complex fracturing schemes. Thus, the results of this model are only valid for cases where the fracture height is confined within the pay zone.

From the result analysis, the “sliding, termination, and crossing” phenomena between hydraulic fractures and natural fractures observed in the model are consistent with Dehgean’s (2020) experimental conclusions. The quantitative relationship between JFE and fracture complexity supplements the qualitative research by Fan et al. (2014), and the law that “increasing E makes fractures longer and narrower” is also consistent with the discrete element simulation results by Li et al. (2021). Additionally, the ABAQUS Cohesive element model used in this study does not require mesh refinement for newly generated fractures, balancing computational accuracy and efficiency. In terms of the priority of parameter adjustment, FFIR is the core adjustable parameter, directly determining the scale and complexity of the fracture network, but it is necessary to balance fracturing effect and wellhead equipment cost. FFV has a limited impact on fracture propagation, and low-viscosity fracturing fluid (1–50 cp) can meet the requirements. As fixed geological parameters, JFE and E can be compensated for their adverse effects on fractures by adjusting FFIR (e.g., FFIR can be increased to 10 m3/min in formations with high JFE).

In terms of engineering application, for formations with high JFE (e.g., deep tight sandstones with strong natural fractures), FFIR can be adjusted to 8–10 m3/min to ensure the complexity of the fracture network. For high-E brittle formations, the viscosity of the fracturing fluid needs to be controlled within 10–50 cp to avoid further narrowing of fractures. For formations with high natural fracture density (the focus of this study), low-viscosity fracturing fluid (1–50 cp) is used, which can not only reduce the cost of fluid preparation but also ensure the fracturing effect. Under real-world field circumstances, geological parameters are typically fixed, so adjustments must be made through engineering parameters, primarily the viscosity and displacement rate of fracturing fluid. A higher Fracturing Fluid Injection Rate (FFIR) exerts a notable influence on fracture geometry and is advantageous for fracturing operations. Nevertheless, increasing FFIR demands excessive equipment and leads to substantial leakage. The selection of this parameter should be based on actual on-site conditions. The viscosity of the fracturing fluid has a relatively minor impact on on-site operations: excessively high viscosity hinders fracture propagation, while overly low viscosity is unfavorable for fracture propping. These results are useful for hydraulic fracturing design in tight formations.

6  Conclusion

This paper investigates the problem of fracture propagation in tight sandstone rock with pre-existing fractures. Firstly, a fluid structure coupling model of fracture intersection was established using the Cohesive module of ABAQUS, and validated through indoor experiments. Afterwards, the cohesive unit is embedded globally and combined with Python to generate a conjugate joint distribution model to establish a dispersed mesh model. Then, based on the establishment of a model, a hydraulic fracturing simulation scheme is used to complete the simulation research on fracture propagation. In summary, the key findings of this study are as follows:

(1) Under normal stress conditions, the numerical simulation results of scenarios with different approach angles and stress differences are highly consistent with the experimental results of scholars such as Zhou and Blanton. In field fracturing operations, when induced fractures encounter natural fractures, phenomena such as sliding, termination, or crossing generally occur, which verifies the model’s ability to characterize actual fracture interaction behaviors.

(2) As joint fracture energy (JFE) increases, joint strength increases simultaneously. Under otherwise unchanged conditions, the length and width of the main fracture decrease significantly; however, the enhanced joint strength promotes the wider diffusion of fracturing fluid along pre-existing fractures, ultimately increasing the complexity of the fracture network. This provides a theoretical basis for fracture network regulation in high-JFE formations.

(3) As the elastic modulus (E) of tight sandstone reservoirs increases, hydraulic fractures exhibit a “long and narrow” characteristic—the number of main fractures remains unchanged while the number of branch fractures increases. Moreover, E has a significant impact on the propagation range of hydraulic fractures, clarifying the key to controlling fracture morphology in high-E brittle formations.

(4) Within the model’s parameter range, when the total fracturing fluid injection volume is constant, increasing the fracturing fluid injection rate (FFIR) can enhance fracture complexity and expand the fracturing control area. However, this requires supporting more robust wellhead equipment, providing a trade-off basis for on-site displacement optimization and equipment selection.

(5) When formation fractures are densely developed, the impact of increasing fracturing fluid viscosity (FFV) within the range of 1–100 cp on fracture propagation is negligible, and the overall effect of FFV is weak. This provides support for the application of low-viscosity fracturing fluid in tight fractured formations.

7  Challenges and Future Prospects

From the perspective of rock mechanics, the core logic behind the narrow width of hydraulic fracturing fractures under high Young’s modulus conditions can be explained through three dimensions: the elastic deformation characteristics of rocks, the matching relationship between the net pressure in the fracture and rock response, and the energy distribution law. Rocks with high Young’s modulus exhibit high rigidity, resulting in minimal elastic strain under the same stress. The “width” of a hydraulic fracture essentially refers to the “tensile elastic deformation” of the rock on the fracture walls under the action of fluid pressure inside the fracture: after fracturing fluid is injected into the fracture, the pressure within the fracture overcomes the in-situ stress of the reservoir to form “net pressure”. This net pressure exerts a “tensile force” on the upper/lower or left/right walls of the fracture, forcing the rock to open sideways, and the distance of this opening is the fracture width. The opening degree of fracture width is jointly controlled by the net pressure driving fracture propagation and the stiffness of rock resisting deformation, and its simplified correlation is:

w∝pnet⋅LE(12)

Among then, w is the average fracture width, m; E is the Young’s modulus of rock, GPa; pnet is the effective net pressure, MPa; L is the fracture length, mm; This proportional relationship indicates that the fracture width is inversely proportional to Young’s modulus—when Young’s modulus of rock increases, its ability to resist elastic deformation is enhanced, and the fracture is difficult to open sufficiently under the same net pressure.

Even after the fracture initiates, the “rigidity” of rocks with high Young’s modulus causes the “energy loss” of fluid pressure in the fracture to be more inclined towards “fracture extension” rather than “fracture opening”. A portion of the energy from the injected fracturing fluid is used to overcome the elastic resistance of the rock to open the fracture (converted into elastic potential energy of the rock), while the other portion is used to overcome the fracture toughness of the rock to expand the fracture (converted into surface energy of the new fracture surface). Rocks with high Young’s modulus have high elastic resistance, so more of the injected energy is consumed in “maintaining fracture extension” (e.g., overcoming the stress shadow effect and passing through porous media). This limits the increase in effective net pressure available for “fracture opening”, further restricting the fracture width. In hydraulic fracturing engineering, the core functions of fracturing fluid viscosity are to “maintain net pressure in the fracture, carry proppants, and reduce fluid loss”. Theoretically, it regulates the fracture propagation morphology by influencing the pressure distribution inside the fracture. However, in practical scenarios, when the dominant effects of inherent reservoir properties, construction parameters, or geological conditions far exceed the influence of viscosity, the impact of fracturing fluid viscosity on fracture propagation morphology decreases significantly. For the results of this study, the existence of natural fracture properties and stress characteristics has, to a certain extent, weakened the influence of fracturing fluid viscosity on the overall fracture morphology. In this case, the fracture morphology is dominated by “natural fracture activation”, and the influence of viscosity is diminished. This phenomenon is relatively common in unconventional reservoirs with well-developed natural fractures, such as shale and coal-measure reservoirs. The core goal of fracturing in such reservoirs is to “activate natural fractures to form complex fracture networks” rather than “rely on fracturing fluid viscosity to create new main fractures”. At this point, the “distribution density, orientation, and connectivity” of natural fractures become the core factors determining the fracture propagation morphology, and the role of viscosity is greatly reduced. In summary, when the effects of factors such as reservoir stress, natural fractures, construction displacement, and fluid loss far exceed the ability of viscosity to regulate net pressure, the influence of viscosity will be reduced.

Acknowledgement: The authors received funding from the National Natural Science Foundation of China. We gratefully acknowledge these contributions.

Funding Statement: This work was financially supported by the National Natural Science Foundation of China (Grant No. U1762212).

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Yujie Yan, Xiongwei Liu; data collection: Yanling Wang, Na An; model construction and calculation: Yujie Yan, Cheng Ji; analysis and interpretation of results: Yujie Yan, Xiongwei Liu; draft manuscript preparation: Yujie Yan, Yanling Wang; manuscript revision: Yujie Yan, Shu Jiang. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: The data supporting the conclusions of this study can be obtained from the corresponding authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

References

1. Zhao J, Wang F, Cai J. 3D tight sandstone digital rock reconstruction with deep learning. J Petrol Sci Eng. 2021;207:109020. doi:10.1016/j.petrol.2021.109020. [Google Scholar] [CrossRef]

2. Man K, Wang J, Su R, Zhou H. Theoretic research on seepage model of “water-rock-structural plane” under 3D stresses. In: Proceedings of the 2011 International Conference on Remote Sensing, Environment and Transportation Engineering; 2011 Jun 24–26; Nanjing, China. p. 1837–9. doi:10.1109/RSETE.2011.5964654. [Google Scholar] [CrossRef]

3. Iyare UC, Frash LP, KC B, Meng M, Li W, Madenova Y, et al. Experimental investigation of shear in granite fractures at Utah FORGE: implications for EGS reservoir stimulation. Geothermics. 2025;131:103344. doi:10.1016/j.geothermics.2025.103344. [Google Scholar] [CrossRef]

4. Ameli P, Elkhoury JE, Morris JP, Detwiler RL. Fracture permeability alteration due to chemical and mechanical processes: a coupled high-resolution model. Rock Mech Rock Eng. 2014;47(5):1563–73. doi:10.1007/s00603-014-0575-z. [Google Scholar] [CrossRef]

5. Fatah A, Al-Yaseri A. Geomechanical integrity and geochemical reactions of shale caprocks for hydrogen storage: a comprehensive review. Fuel. 2025;400:135728. doi:10.1016/j.fuel.2025.135728. [Google Scholar] [CrossRef]

6. Uwakwe OC, Riechelmann S, Mueller M, Reinsch T, Balcewicz M, Igbokwe OA, et al. Scaling in fractured geothermal carbonate reservoir rocks: an experimental approach. Geothermics. 2025;125:103199. doi:10.1016/j.geothermics.2024.103199. [Google Scholar] [CrossRef]

7. Dehghan AN. An experimental investigation into the influence of pre-existing natural fracture on the behavior and length of propagating hydraulic fracture. Eng Fract Mech. 2020;240:107330. doi:10.1016/j.engfracmech.2020.107330. [Google Scholar] [CrossRef]

8. Fan T, Zhang G, Cui J. The impact of cleats on hydraulic fracture initiation and propagation in coal seams. Pet Sci. 2014;11(4):532–9. doi:10.1007/s12182-014-0369-7. [Google Scholar] [CrossRef]

9. Li Y, Hubuqin, Wu J, Zhang J, Yang H, Zeng B, et al. Optimization method of oriented perforation parameters improving uneven fractures initiation for horizontal well fracturing. Fuel. 2023;349(5):128754. doi:10.1016/j.fuel.2023.128754. [Google Scholar] [CrossRef]

10. Shrivastava S, Banerjee A, Singh A, Singh MK, Pasupuleti S. Copula-based dependency modelling of hydraulic properties for non-linear filtration through porous media. Powder Technol. 2025;460:121069. doi:10.1016/j.powtec.2025.121069. [Google Scholar] [CrossRef]

11. Tang X, Yang H. Research on energy and instability criterion of sandstone hydraulic fracturing in the Three Gorges Reservoir area, Chongqing, China. Geofluids. 2023;2023:9985452. doi:10.1155/2023/9985452. [Google Scholar] [CrossRef]

12. Edirisinghe EAAV, Perera MSA. Review on the impact of fluid inertia effect on hydraulic fracturing and controlling factors in porous and fractured media. Acta Geotech. 2024;19(12):7923–65. doi:10.1007/s11440-024-02389-7. [Google Scholar] [CrossRef]

13. Li Y, Hu W, Zhang Z, Zhang Z, Shang Y, Han L, et al. Numerical simulation of hydraulic fracturing process in a naturally fractured reservoir based on a discrete fracture network model. J Struct Geol. 2021;147:104331. doi:10.1016/j.jsg.2021.104331. [Google Scholar] [CrossRef]

14. Zhu H, Deng J, Jin X, Hu L, Luo B. Hydraulic fracture initiation and propagation from wellbore with oriented perforation. Rock Mech Rock Eng. 2015;48(2):585–601. doi:10.1007/s00603-014-0608-7. [Google Scholar] [CrossRef]

15. Zhang J, Li X, Wang Q, Zheng J, Huai Q. Experimental study on permeability enhancement of tight sandstone by hydraulic fracturing under the true tri-axial system. Petrol Sci Technol. 2025;26:1–20. doi:10.1080/10916466.2025.2465848. [Google Scholar] [CrossRef]

16. Li J, Dong S, Hua W, Li X, Guo T. Numerical simulation of temporarily plugging staged fracturing (TPSF) based on cohesive zone method. Comput Geotech. 2020;121:103453. doi:10.1016/j.compgeo.2020.103453. [Google Scholar] [CrossRef]

17. Ma T, Jiang L, Shen W, Cao W, Guo C, Nick HM. Fully coupled hydro-mechanical modeling of two-phase flow in deformable fractured porous media with discontinuous and continuous Galerkin method. Comput Geotech. 2023;164:105823. doi:10.1016/j.compgeo.2023.105823. [Google Scholar] [CrossRef]

18. Zhang H, Chen J, Gong D, Liu H, Ouyang W. Effects of fracturing parameters on fracture network evolution during multicluster fracturing in a heterogeneous reservoir. Comput Geotech. 2023;159:105474. doi:10.1016/j.compgeo.2023.105474. [Google Scholar] [CrossRef]

19. Carrier B, Granet S. Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model. Eng Fract Mech. 2012;79:312–28. doi:10.1016/j.engfracmech.2011.11.012. [Google Scholar] [CrossRef]

20. Liao J, Zhang Z, Tang H, Yang J, Zhang X, Wang B. Simulation of fracturing and well pattern optimization of fractured tight sandstone reservoirs. Front Earth Sci. 2022;10:873617. doi:10.3389/feart.2022.873617. [Google Scholar] [CrossRef]

21. Gao Q, Han S, Cheng Y, Yan C, Sun Y, Han Z. Effects of non-uniform pore pressure field on hydraulic fracture propagation behaviors. Eng Fract Mech. 2019;221:106682. doi:10.1016/j.engfracmech.2019.106682. [Google Scholar] [CrossRef]

22. Wang H, Wang GL, Yue GF, Gan HN. Numerical simulation of granite hydraulic fracture propagation under the influence of natural fractures. Acta Geol Sin. 2020;94(7):2124–30. [Google Scholar]

23. Zimmerman RW, Bodvarsson GS. Hydraulic conductivity of rock fractures. Transp Porous Medium. 1996;23(1):1–30. doi:10.1007/BF00145263. [Google Scholar] [CrossRef]

24. Zhang P, Pu C, Shi X, Xu Z, Ye Z. The numerical simulation and characterization of complex fracture network propagation in multistage fracturing with fractal theory. Minerals. 2022;12(8):955. doi:10.3390/min12080955. [Google Scholar] [CrossRef]

25. Zhou J, Chen M, Jin Y, Zhang GQ. Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs. Int J Rock Mech Min Sci. 2008;45(7):1143–52. doi:10.1016/j.ijrmms.2008.01.001. [Google Scholar] [CrossRef]

26. Blanton TL. Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs. In: Proceedings of the SPE Unconventional Gas Technology Symposium; 1986 May 18–21; Louisville, Kentucky. p. 613–27. doi:10.2118/15261-ms. [Google Scholar] [CrossRef]

27. Távara L, Mantič V, Salvadori A, Gray LJ, París F. Cohesive-zone-model formulation and implementation using the symmetric Galerkin boundary element method for homogeneous solids. Comput Mech. 2013;51(4):535–51. doi:10.1007/s00466-012-0808-5. [Google Scholar] [CrossRef]

28. Wang H. Hydraulic fracture propagation in naturally fractured reservoirs: complex fracture or fracture networks. J Nat Gas Sci Eng. 2019;68:102911. doi:10.1016/j.jngse.2019.102911. [Google Scholar] [CrossRef]

29. Daneshy A. Three-dimensional analysis of interactions between hydraulic and natural fractures. In: Proceedings of the SPE Hydraulic Fracturing Technology Conference and Exhibition; 2019 Feb 5–7; The Woodlands, TX, USA. doi:10.2118/194335-ms. [Google Scholar] [CrossRef]