Hydraulic fracturing technology entails injecting fracturing fluid at high pressure into targeted formations through a wellbore, with the aim of generating highly conductive hydraulic fractures (Zhang et al. 2021). Over the past few decades, hydraulic fracturing has emerged as a highly effective stimulation technique for the development of unconventional petroleum reservoirs characterized by low permeability and porosity (Zheng et al. 2022). Predicting fracture propagation trajectories and fracture network morphology is crucial for optimizing the design of hydraulic fracturing operation schemes. Accurate predictions enable effective planning and control of the hydraulic fracturing process, leading to improved reservoir productivity and recovery. As the utilization of hydraulic fracturing technology in tight sandstone and shale gas reservoirs has become more widespread, there has been a growing focus on understanding the mechanisms of fracture propagation and the formation of fracture networks (Guo et al. 2015; Suboyin et al. 2020).

During hydraulic fracturing, the fracturing fluid seeps from the fracture surface into the rock pores, exerting the action of seepage force on the rock skeleton (Rozhko 2007). Indeed, in geotechnical mechanics, the seepage force is defined as the drag force exerted by flowing fluid on a unit volume of rock or soil (Taylor 1948). To maintain force balance, the seepage force acting on the rock element needs to be counterbalanced by a corresponding change in effective stress. This adjustment is necessary to ensure equilibrium within the system (Wang 2000). In the past decades, many numerical and experimental researches have been conducted to analyze the effect mechanism of seepage force on the rock or soil failure (AlKhafaji et al. 2020; Bruno and Nakagawa 1991; Cobbold and Rodrigues 2007; Mourgues and Cobbold 2003; Pan et al. 2017; Rozhko 2007; Saada et al. 2012; Zareifard and Fahimifar 2015; Zou et al. 2016; Zeng et al. 2021). The research findings have shown that the seepage force has a significant effect on the deformation characteristics and failure modes of the rock or soil. In a recent study, Zhou et al. (2020) developed an analytical solution to calculate the stress field induced by seepage force around a borehole during steady-state flow conditions. By superimposing the derived stress field, it was determined that the seepage force plays a crucial role in influencing the effective stress field around the wellbore and the initiation of hydraulic fractures. When fracturing fluid enters into the reservoir pores, the effects of seepage forces should not be neglected in the analysis of fracture initiation and propagation mechanisms (Rahimi-Aghdam et al. 2019; Wang and Zhou 2024; Zhou et al. 2021).

Numerous studies have been conducted to investigate the influence mechanism of fracturing fluid viscosity and injection rate on fracture propagation behaviors. Chen et al. (2015) investigated the effects of three different viscous fracturing fluids on anisotropic granite fracture propagation through true triaxial laboratory experiments. The results show that lower viscous fracturing fluids are beneficial to form complex fracture networks in the granite specimens, and supercritical CO₂ induces the most fracture branches from the primary fractures. Tan et al. (2017) investigated the fracture initiation and vertical propagation behavior and found that the product of fluid viscosity and injection rate has an optimal value favorable to the formation of complex fracture networks. Through fracturing experiments, Wanniarachchi et al. (2018) found that foam fracturing tends to induce complex twisted fractures with larger surface areas than water-based fracturing. In a study conducted by Hou et al. (2019), the behavior of HF propagation was analyzed in various lithologic rocks using larger rock specimens. The experimental findings revealed that less viscous fracturing fluid and a smaller injection rate were conducive to the formation of complex fractures. On the other hand, a higher injection rate was found to be favorable for fracture development. Liu et al. (2021) conducted a study to examine the impact of fracturing fluid viscosity on directional hydraulic fracturing. Their findings indicated that both fracture initiation pressure and fracture length exhibit an increase with an increase in fracturing fluid viscosity.

The studies mentioned above have emphasized two crucial factors: fracturing fluid viscosity and injection rate, which are closely linked to the seepage force and have a significant impact on HF propagation. However, it is noteworthy that most studies have overlooked the analysis of their effects on the seepage force and the role of this force in HF propagation.

This paper commences by introducing the mechanism of seepage force and subsequently presents a discrete element fluid–solid coupling algorithm that incorporates its effects. Subsequently, we construct a dynamic 2D discrete-element model using the particle flow code (PFC2D) to simulate fracture propagation under the influence of seepage force. A comprehensive series of simulations is conducted to investigate the impact of seepage force in hydraulic fracturing scenarios within both homogeneous and heterogeneous reservoirs. Finally, we conduct an analysis to examine the influence of fracturing fluid viscosity and injection rate on hydraulic fracture propagation while considering the presence of seepage force.

Seepage force is a fundament concept in rock and soil mechanics, and it is defined as the drag force of the flowing fluid acting on the particles in the unit volume of rock or soil under the action of the head difference (Taylor 1948; Das 2021). Figure 1 shows the soil seepage failure test under a steady-state flow condition.

Schematic diagram of the effect of the seepage force

During the fluid flow, the fluid pressure will drop due to the drag force on the soil particles, resulting in a head loss. The equilibrium condition of the force of the water body can be described as:

where h₁ is the pressure head at the inlet, h₂ is the pressure head at the outlet, and L is the length of the water body. γ_w is the unit weight of the fluid γ_w = ρg, where ρ is the density of the fluid and g is the acceleration of gravity. j represents the seepage force, which is equal to the drag force applied to the fluid by the soil particles in the unit volume of soil.

From Eq. 1, the seepage force is defined as:

where i represents the hydraulic gradient, P denotes the pore pressure. The seepage force is aligned with the direction of fluid flow.

Based on the aforementioned analysis, it is evident that the seepage force directly impacts the soil or rock skeleton and necessitates an adjustment in effective stress to maintain stress equilibrium. The effective stress controls the deformation and failure behavior of the rock and soil, so the seepage force will have a significant impact on the HF propagation process.

A fully dynamic 2D discrete-element model with particle flow code (PFC2D) is developed for the simulations of HF propagation considering the seepage force. The particle flow code includes the discretization of rock material into many circular particles that are connected through different types of contact bonds. Based on the traditional fluid–solid coupling algorithm in bonded granular matrix proposed by Cundall (2000), the fully coupled hydro-mechanical simulation can be conducted in the particle flow model. In this study, we have enhanced the traditional algorithm to account for the influence of seepage force on rock deformation during fluid seepage through pore throats. Additionally, we have developed a dedicated discrete element method (DEM) program to investigate the process of fracture propagation.

3.1 Improved fluid–solid coupling algorithm

The fluid–solid coupling algorithm generates fluid domain networks by searching all pores of the bonded assembly of particles. The fluid domain networks act as flow channels and reservoirs for fluid storage in porous media during fluid injection. As depicted in Fig. 2a, the polygons connected by black lines represent the fluid domains. The blue line represents the fluid flow channel between adjacent fluid domains, while the blue square signifies the center of the fluid domain. The fluid domain is responsible for storing fluid and carrying pore pressure. In the fluid–solid coupling process, the fluid domain of the wellbore is first selected as the fluid injection point. Then, when the fracturing fluid is injected, the fluid will flow into the fluid domain around the injection point along the flow channel under the action of pressure difference and change the pore pressure of the fluid domain. The flowing fluid will exert water pressure and friction on the rock particles, causing the movement of the particles and the failure of the contact bond. The change in particle position and the failure of contact bonds will in turn affect the aperture of the flow channel and the fluid flow rate. Through alternating iterations, the coupling between fluid and solid particles will eventually be achieved. The governing equations for fluid–solid coupling are as follows: Eqs. (3) to (7) represent fluid governing equations, Eqs. (8) to (12) correspond to mechanical governing equations, and Eq. (13) outlines the contact bond fracture criterion.

Fluid domain networks modelled in PFC2D

As depicted in Fig. 2c, the fluid flow channel between two particles is considered to be a pair of parallel plates with an aperture w. The flow rate Q of the fluid within the flow channel follows the Poiseuille equation.

where w represents the aperture of the flow channel, L denotes the length of the flow channel, μ represents the viscosity of the fluid, and Δp signifies the pressure difference between the fluid domains.

Considering the radius difference between the two particles, L is obtained from the harmonic average of the radius of two particles (see Fig. 2c).

The aperture (w) of the flow channel serves as an indicator of the model's permeability and plays a crucial role in governing the process of fluid–solid coupling. When the rock particles are subjected to compressive stress, the pore throats are compacted and the fluid flow is restricted. While under the effect of tensile stress, the pore throats become larger and the fluid flow is enhanced. If the contact bond between particles breaks to form a fracture, the fluid flow capacity in the fractured channel should be much greater than that of the rock matrix. Based on this, the aperture w of the flow channel in three different cases is set according to the following criteria:

where w_o represents the initial aperture of the flow channel, λ is the dimensionless safety factor, and d is the distance between two particles when contact is not broken under tensile stress. w_gap signifies the distance between particles after contact bond failure, and w_HFo represents the initial aperture of the fracture flow channel. F denotes the current contact normal force, while F₀ represents the contact normal force when the aperture decreases to half of its original value.

The initial aperture of the flow channel is generally set by an empirical formula, which makes it difficult to calibrate the rock permeability. This study defined the initial aperture w_o by the maximum and minimum particle radii of the model.

where α is a dimensionless coefficient through which the reservoir permeability is controlled, r_max is the maximum particle radius, and r_min is the minimum particle radius.

The fluid volume stored in the fluid domain changes during fluid flow, which alters the pore pressure of the fluid domain. The pressure change can be obtained from the following continuity equation:

where K_f is the fluid bulk modulus; V_d is the apparent volume of the fluid domain;

ΔV_d is the volume change of the fluid domain; Δt is flow time step. ∑ is the summation of the volumes of all fluids flowing into or out of the fluid domain during Δt time step.

As shown in Fig. 2a, fluid will exert force on rock particles when flowing in the fluid domain network. The traditional fluid–solid coupling algorithm only considers the water pressure on the particle surface. As shown in the local enlarged view of the fluid domain network (Fig. 2b), fluid flowing into the fluid domain induces changes in pore pressure and imparts normal water pressure onto the surfaces of rock particles. By integrating the water pressure on the particle surface, the resultant force F_P acts on the particle center in the form of body force.

where P is the pore pressure of the fluid domain, r is the particle radius, n_i is the normal vector of the line connecting two adjacent contact points of the particle, and s is the length of the line.

Not only the water pressure but also a frictional force acts on the particles during the fluid flow (Shimizu et al. 2011; Liu et al. 2019). However, the traditional algorithm ignores the friction force on the particle surfaces, resulting in flawed solutions in simulating the seepage force. In this study, we incorporate the frictional effect of fluid flow into our algorithm to accurately simulate the effect of seepage force.

As shown in Fig. 2c, to calculate the friction force applied to the two parallel plates, the fluid velocity is first obtained by the Poiseuille formula (Munson et al. 2013):

Then take the derivative of velocity u with respect to y:

The shear stress exerted by the fluid is further obtained:

Set y = 0 or y = w, integrate the shear stress along the x-direction to get the friction force F_S applied by the fluid:

F_S is the frictional force exerted on the particle surface during the fluid flow. When the resultant force F_P of the water pressure and the friction force F_S acts together on the particle center in the form of the body force, the seepage force exerted by the fluid can be accurately simulated.

The particle flow method updates the contact forces between particles based on the Force–Displacement Law of the contact model, and then updates the positions of the particles using Newton's second law of motion. In this paper, the parallel bond model is employed to bond the circular particles, as it has proven successful in numerical simulations of rock deformation and failure (Potyondy and Cundall 2004; Shimizu et al. 2011; Wang et al. 2022; Qi et al. 2023). The mechanical criteria of the parallel bond model can be found in the PFC2D user manual. For the parallel bond model, the fracture tip extension criterion follows the pattern below:

where τ_max and σ_max are the current maximum shear stress and maximum tensile stress specified for the parallel bond, τ_c and σ_c are the shear strength and tensile strength assigned to this parallel bond. According to Eq. 13, bond failures and damages occur when the maximum tensile stress surpasses the tensile strength or when the maximum shear stress exceeds the shear strength.

3.2 Verification of Seepage force simulation

3.2.1 Analytical verification

To verify the accuracy of the seepage force simulation results, we developed a DEM program based on an improved fluid–solid coupling algorithm to simulate the Darcy steady-state seepage experiment. As depicted in Fig. 3a, a core model with length L and width W is established, and a fluid domain network is generated using the improved fluid–solid coupling algorithm. The fluid domain network with length ΔL at the left and right sides of the core was selected as the inlet and outlet in the process of fluid displacement. The pressure at the outlet was set to 0 MPa, and the constant pressure fluid displacement was carried out from the inlet. During fluid displacement, the flow rates at the inlet and outlet are monitored, and the permeability is calculated using Darcy's law. As shown in Fig. 3c, when the inlet and outlet rates are stable and close, the flow reaches a steady state. Figure 3a, b illustrate the corresponding pore pressure distribution and contact force chain under the steady state. The pore pressure evenly distributes along the horizontal direction, and the contact force chain gradually changes from compression to tension under the effect of seepage force.

Darcy steady state seepage experiment simulation

In the steady state, according to Eqs. 1–2, the total seepage force in the stable seepage section is:

where p_inlet is the pressure of the inlet, p_outlet is the pressure of the outlet, and V_Total is the total volume of the stable seepage section.

The sum of the forces on all particles in the stable seepage section can be calculated, and the calculation results are compared with the traditional algorithm and the analytical solution in Fig. 4. Figure 4 shows that with the increase of the pressure difference between the inlet and the outlet, the total seepage force in the stable seepage section increases linearly. The seepage force in the steady state depends on the pressure difference. In comparison to the traditional algorithm, the improved algorithm can more accurately simulate the seepage force exerted on the rock particles during fluid seepage.

Validation of seepage force simulation results. (L = 5 cm; W = 2.5 cm; ΔL = 1 cm; μ = 1 mPa·s)

3.2.2 Experimental verification

In a previous study by Bruno and Nakagawa (1991), the propagation path of hydraulic fractures under non-uniform pore pressure in rock formations was investigated. Water injection holes were pre-set in the rock sample to simulate non-uniform pore pressure fields on both sides of the main fracture during hydraulic fracture propagation. The experimental results shows that the hydraulic fracture propagates in a straight line in the region where the pressure gradient induced by the injection hole has not spread, while in the range where the pressure gradient field created by the injection hole has spread, the hydraulic fracture deviates significantly towards the side with a high pore pressure gradient, exhibiting a curved pattern.

To compare our numerical simulation results with the experimental findings, we used an improved algorithm to establish a non-uniform permeability model that simulates hydraulic fracture propagation, as demonstrated in Fig. 5. The left side of the blue baseline represents the permeable reservoir with the fluid domain channel open, while the right side represents the impermeable reservoir with the fluid domain channel closed. The color scale caliper on the left represents the magnitude of the seepage force vector, while the color scale caliper on the right signifies dimensionless pore pressure. Our numerical simulation results demonstrate that during the injection process, a significant pore pressure gradient field forms on the permeable reservoir side. The pore pressure field generated by the injection point and the corresponding pore pressure gradient extend farther from the injection point, resulting in a higher seepage force field. As one moves away from the injection point, both the range of the pore pressure field and the pore pressure gradient gradually diminish, leading to a gradual attenuation of the seepage force field. Once the hydraulic fracture initiates from the injection point, it deflects toward the side characterized by a higher pore pressure gradient, displaying a curved pattern akin to experimental observations. As the pore pressure gradient decreases and the seepage force weakens with increasing distance from the injection point, the hydraulic fracture gradually transitions to a more linear expansion pattern.

Comparison of numerical simulation result and experimental result of hydraulic fracture propagation under non-uniform pore pressure field

4.1 Rock specimen model setup

A square-shaped rock specimen model is established to simulate the HF propagation in a homogeneous reservoir as shown in Fig. 6. The model dimensions are set at 0.25 m width and height, with a wellbore radius of 5 mm. The parallel bond model is used to assemble particles as it has been successfully applied in the simulations (Potyondy and Cundall 2004; Yoon 2007; Zhou et al. 2016). A uniform particle size distribution is adopted in this rock model, where the ratio of maximum to minimum particle radius is 1.66. The servo mechanism is used to apply confining pressure to the model including the maximum in-situ horizontal principal stress SH in the x-direction and the minimum in-situ horizontal principal stress Sh in the y-direction. The maximum horizontal principal stress SH is set to twice the minimum horizontal principal stress Sh to prevent the formation of branch fractures and the redirection of HFs. To ensure a smooth wellbore surface and avoid unnecessary stress concentration, particles of the same radius are arranged in a circular to form the wellbore for fracturing fluid injection. The microscopic input parameters adopted in this model are calibrated by numerical simulations of uniaxial compression, uniaxial tension, and Darcy seepage tests. In this paper, macroscopic mechanical properties of a tight sandstone reservoir are used to calibrate the microscopic input parameters of the model (Zhang et al. 2020). The microscopic input parameters of the model and the calibrated macroscopic mechanical properties are presented in Table 1.

Rock specimen model for the hydraulic fracturing simulation

4.2 Effect of fracturing fluid viscosity

The effect of seepage force on rock skeleton mainly depends on rock permeability, fluid viscosity, and injection rate. Based on the calibrated permeability in the Darcy flow simulation experiment, we studied the fracture propagation behavior with three different formation permeabilities, 2 mD, 20 mD, and 200 mD, under different fluid viscosities and injection rates. The wellbore was considered as an injection point at constant injection rate. An impermeable rock model is generated by closing the flow channels for comparison study.

As shown in Fig. 7, the HFs propagate along the direction parallel to the maximum horizontal principal stress with different fracturing fluid viscosities, and the fracture geometries are slightly different. When the fracturing fluid infiltrates into the rock pores, the presence of the seepage force leads to a reduction in the fracture length compared to the case without permeability. The lower the fracturing fluid viscosity is, the shorter the fracture length is. Figure 8 shows the number of fractures formed by contact bond failure during the hydraulic fracturing. Since there are no natural cracks and faults in the rock model, there are no fractures formed by contact bond shear failure. Under the action of water pressure, all fractures shown in Fig. 8 are formed due to contact bond tensile failure. The fracture length increases with the number of contact bond failures. According to Fig. 8, in reservoirs with higher permeability values of 20 mD and 200 mD, if the viscosity of the fracturing fluid is excessively low, the fluid rapidly permeates from the wellbore into the reservoir. Consequently, HFs cannot form at a low injection rate. The number of fractures in three different permeability reservoirs increases nonlinearly with the increasing fracturing fluid viscosity. High-permeability reservoirs require more viscous fracturing fluid to mitigate the role of seepage force and form fractures with the same length as low-permeability reservoirs, while low-permeability reservoirs can generate relatively longer fractures with lower fracturing fluid viscosity.

Fracture propagation at the same time step using different fluid viscosities. (Reservoir permeability:2 mD; injection rate: 2.5 × 10⁻³ m³/s)

The number of contact bond failure in the rock model under the same time step and different fluid viscosity. (Injection rate: 2.5 × 10⁻³ m³/s)

To analyze the effect of fracturing fluid viscosity on the seepage force, the distribution field of fluid domain pressure, ball appliedforce, and induced stress field around HF are studied. Figure 9 illustrates the fluid domain pressure distribution around the HF. The value of the color scale caliper represents the dimensionless fluid pressure P_d, P_d = 0.2 × (fluid domain pressure/water injection point pressure). Figure 10 presents the magnitude of the force vector exerted by the fluid on the rock particles during fracture propagation. Meanwhile, Fig. 11 displays the induced stress field around HF, with positive values representing tensile stress.

Fluid domain pressure distribution field under different fluid viscosities. (Reservoir permeability: 20 mD; Injection rate: 2.5 × 10⁻³ m³/s)

Ball appliedforce distribution field under different fluid viscosities. (Reservoir permeability: 20 mD; Injection rate: 2.5 × 10⁻³ m³/s)

Induced stress field around HF under different fluid viscosities. (Reservoir permeability: 20 mD; Injection rate: 2.5 × 10⁻³ m³/s)

Figure 9 reveals that the fluid pressure within the HF exhibits a gradient distribution, with pressure increasing closer to the injection point of the wellbore. In contrast to impermeable reservoirs, the seepage of fluid into rock pores significantly alters the pore pressure distribution surrounding HFs. A reduction in the viscosity of the fracturing fluid leads to a broader impact on the pore pressure field surrounding the HF. In impermeable reservoirs, where fluids cannot penetrate rock pores, there is a higher fluid pressure within the HFs.

The magnitude of the force exerted by the fluid on the rock particles determines whether the contact bond can undergo tensile failure and whether HFs can initiate and propagate. Seepage force is a body force exerted on the rock skeleton when the fracturing fluid flows into the reservoir pores along the fracture surface. However, in impermeable reservoirs, where the fluid cannot penetrate the rock pores, the water pressure within the HF merely exerts a surface force on the fracture wall. As illustrated in Fig. 10, the surface force exerted by the fluid in impermeable reservoirs on the fracture wall surpasses the seepage force acting on the rock skeleton by the fluid in permeable reservoirs. As the viscosity of the fracturing fluid decreases, the primary effect zone of the seepage force around the HF expands. However, when the viscosity of the fracturing fluid is low, despite the extensive range of the seepage force, its actual magnitude remains limited. This limitation results in insufficient normal force required for the tensile fracture of the contact bond, consequently leading to shorter HF lengths.

Figure 11 reveals that the tensile stress concentration area resides around the crack tip, while the compressive stress concentration area surrounds the crack wall. In comparison to the impermeable reservoir in the control group, the seepage force exerted by the fracturing fluid infiltrating into the reservoir pores profoundly alters the induced stress field around the fracture. As the viscosity of the fracturing fluid decreases, the tensile stress at the fracture tip weakens. Upon comparing Figs. 9, 10, and 11, it becomes apparent that the influence range of the seepage force around the HF on the induced stress field greatly exceeds the seepage range of the fracturing fluid.

4.3 Effect of injection rate

Figure 12 represents the simulation results of fracture propagation in low-permeability reservoirs (2mD) and high-permeability reservoirs (200mD) under different injection rates.

Simulation results of fracture propagation under the same time step and different injection rates

The fracture length with different permeabilities increases with the increasing injection rate. To form fractures of the same length, the injection rate required for low-permeability reservoirs is significantly smaller than that of high-permeability reservoirs. Even if the reservoir is homogeneous, higher injection rates may produce small amounts of secondary fractures around HFs. Short branch fractures may form from the wellbore when high-permeability reservoirs utilize large injection rates to form long fractures, as shown in Fig. 12b.

As shown in Fig. 13, the pore pressure distribution field around the HF does not change significantly with the increase of injection rate. Compared with the effect of fracturing fluid viscosity, injection rate exerts a minimal influence on the seepage effect of fluid in the reservoir. Figure 14 demonstrates that with the increase of injection rate, the effect range of seepage force does not change significantly, and the main areas affected by seepage force around HFs are the same. However, as the injection rate increases, the value of the seepage force on the particles increases, and the possibility of contact bond tensile failure increases, so the HF length is longer at higher injection rates.

Pore pressure distribution around HFs under different injection rates. (Reservoir permeability: 20 mD; Fluid viscosity: 10 mPa·s)

Seepage force distribution around HFs under different injection rates. (Reservoir permeability: 20 mD; Fluid viscosity: 10 mPa·s)

During hydraulic fracturing, a large injection rate may cause the wellbore pressure to exceed the safety threshold, resulting in casing deformation. Therefore, during the hydraulic fracturing process, it is recommended to utilize the viscosity of the fracturing fluid to control the impact of seepage force and regulate the length of the HF. This approach has the potential to effectively mitigate construction risks associated with hydraulic fracturing operations.

It is well known that reservoir heterogeneity has a significant influence on the HF propagation behaviors (Ouchi et al. 2017; Zou et al. 2020). In the section analyzing the impact of seepage force on fracture propagation, we assessed the influence of rock cementation strength and permeability heterogeneity on fracture propagation using the rock sample model established in Sect. 4.1.

Figure 15 depicts the fracture propagation behaviors within a heterogeneously cemented rock. To ensure uniformity in modeling, all rock particles are bonded using the parallel bond model to avoid the influence of different bonding models on fracture propagation. The heterogeneity in rock cementation is manifested through variations in the tensile strength of parallel bonds across different regions. Furthermore, interfaces with different strengths are also bonded using the parallel bond model. The upper green is strong cementation, the lower red is weak cementation, and the thin blue ring in the middle is the compaction zone near the borehole. The tensile strength of the parallel bond was used to control the strength of rock cementation. In cases where the difference between the two horizontal principal stresses is significant (e.g., 10 MPa), the heterogeneity of reservoir cementation strength exhibits minimal influence on the trajectory of fracture propagation. The HFs propagate in the direction of the maximum horizontal principal stress consistent with that of the homogeneous reservoir, and the fracture propagation trajectory is controlled by the stress rather than the seepage force. As the difference between the two horizontal principal stresses decreases to 5 MPa, it is observed that one side of the HF within the impermeable reservoir starts to deflect towards the weakly cemented side. Considering the fluid seepage into the reservoir, the effect of seepage force significantly enhances the tendency of fracture deflection. When the difference of the two horizontal principal stresses is close (2 MPa, 0 MPa), branch fractures are formed from the wellbore. In comparison to impermeable reservoirs, when the viscosity of the fracturing fluid is lower, the branch fractures are more prone to deflect towards the weakly cemented side due to the stronger influence of the seepage force. The stronger the seepage force is, the easier for the HFs on the weakly cemented side to form branch fractures at the tip is. The effect of seepage force is conducive to the formation of complex fracture networks in heterogeneously cemented reservoirs under a low horizontal in-situ stress difference.

Fracture propagation in heterogeneous cemented reservoirs under different horizontal in-situ stress differences and fluid viscosities. (SH: 20 MPa)

Figure 16 demonstrates the fracture propagation in heterogeneous permeability reservoirs. The blue line represents the baseline, and the initial aperture of the flow channel w_o was used to control the heterogeneity of the reservoir permeability. Above the baseline is the high-permeability fluid domain network generated by the improved fluid–solid algorithm, and below is the low-permeability fluid domain network, as observed in Figs. 16a, b. When the reservoir is impermeable or the fracturing fluid viscosity is high, the HFs mainly propagate along the direction parallel to the maximum horizontal principal stress, and the permeability heterogeneity has little effect on the fracture propagation trajectory. When the viscosity of the fracturing fluid reduces to 1 mPa·s, the HFs begin to deflect to the side of the high-permeability reservoir under the effect of seepage force. For ultra-low viscosity fracturing fluids (0.01 mPa·s), HFs begin to deflect significantly toward the high-permeability reservoir, and the smaller the horizontal in-situ stress difference, the greater the degree of HF deflection. Due to the strong seepage force, the ultra-low viscosity fracturing fluid forms many secondary fractures around the HFs of the reservoir. During the process of hydraulic crack propagation, a tensile stress concentration zone forms at the fracture tip, accompanied by a compressive stress concentration zone on both sides of the fracture wall. In regions characterized by concentrated compressive stress, the initiation and propagation of secondary fractures are constrained. Although the seepage force near the primary fracture wall is substantial, it is only after escaping the compressive stress concentration area that the potential for secondary fracture generation emerges. Consequently, as depicted in Fig. 16, the initiation of secondary fractures occurs at a certain distance from the primary hydraulic fracture.

Fracture propagation in heterogeneous permeability reservoirs under different horizontal in-situ stress differences and fluid viscosities. a High permeability fluid domain network, w_o = (r_max + r_min)*0.01; b Low permeability fluid domain network, w_o = (r_max + r_min)*0.001

The cementation strength of the rock is closely related to reservoir permeability. A rock with weaker cementation strength exhibits higher reservoir permeability (Kadhim et al. 2013; Lyu et al. 2019). Combining the simulation results in Figs. 14 and 15, it can be seen that the effect of seepage force is conducive to the formation of complex fracture networks in reservoirs with weak cementation and high permeability under low horizontal in-situ stress difference.

This paper developed a discrete element fluid–solid coupling algorithm capable of accurately simulating the seepage force exerted on the rock skeleton during the infiltration of fracturing fluid. The numerical modeling results exhibit satisfactory agreement with analytical solutions and previously published experimental findings, affirming the reliability of the model. Using this algorithm, the effect of seepage force on fracture propagation behavior in homogeneous and heterogeneous reservoirs was investigated. The simulation results show that:

1. (1)

   The viscosity of the fracturing fluid plays a significant role in determining the range and magnitude of the seepage force effect. When the viscosity of the fluid decreases, resulting in an increased seepage force, the length of fractures in homogeneous reservoirs tends to decrease.

2. (2)

   The fluid injection rate has a limited effect on the range of the seepage force effect, but does impact its magnitude. An increase in fluid injection rate results in a longer fracture length in homogeneous reservoirs.

3. (3)

   As the fracturing fluid seeps into the rock pores, the seepage force effect induces a phenomenon where the hydraulic fractures tend to deflect towards the side of the reservoir with weaker cementation strength and higher permeability.

4. (4)

   The influence of the seepage force promotes the development of intricate fracture networks within heterogeneous reservoirs characterized by varying cementation strength and permeability, particularly in scenarios where there are low horizontal in-situ stress differences.

Future work will focus on examining the interplay between naturally occurring fractures and hydraulic fractures, with a view to understanding the seepage force control mechanisms that govern the formation of complex fracture networks. The ultimate goal of this line of inquiry is to develop a comprehensive set of theoretical principles for controlling seepage force in hydraulic fracturing operations.