A Fully Coupled Discontinuous Deformation Analysis Model for Simulating Hydromechanical Processes in Fractured Porous Media

Abstract

Numerical simulations play a key role in the optimization of fracturing operation designs for unconventional reservoirs. Because of the presence of numerous natural discontinuities and pores, the rock masses of reservoirs can be regarded as fractured porous media. In this paper, a fully coupled discontinuous deformation analysis model is newly developed to simulate the hydromechanical processes in fractured and porous media. The coupling of fracture seepage, pore seepage, and fracture network propagation is realized under the framework of DDA. The developed model is verified with several examples. Then, the developed DDA model is applied to simulate the hydraulic fracturing processes in fractured porous rock masses, and the effects of rock mass permeability on fracturing are investigated. Our findings suggest that high rock permeability may inhibit the stimulation of fracture networks, while increasing the viscosity of fracturing fluids can enhance the fracturing efficiency. This study provides a valuable numerical tool for simulating hydromechanical processes in fractured and porous media and can be used to analyze various geo-mechanical problems related to fluid interactions.

1. Introduction

Unconventional reservoirs including shale oil and shale gas reservoirs require special recovery operations to improve the relatively low permeability of the rock masses. The trapped hydrocarbons in the reservoirs can be exploited through reservoir fracturing techniques for the generation of fracture networks with high permeability. While water-based hydraulic fracturing is the most commonly used technique for stimulating these reservoirs, it necessitates significant water consumption and may lead to the swelling of clay minerals. Consequently, waterless fracturing techniques have been proposed to improve the fracturing efficiency and reduce environmental concerns, including foam fracturing, CO₂ fracturing, explosive fracturing, electric fracturing, etc. [1]. Recently, a novel method of cryogenic fracturing using liquid nitrogen has garnered increasing attention from researchers due to its advantages in creating complex fractures, preventing fracturing closure, and offering environmental benefits [2,3,4,5]. Although these fracturing techniques differ in the choice of fracturing fluids and energy sources, they all involve similar physical processes, including fracture propagation, fluid (water, gas, oil, etc.) flow, and rock deformation, which represent a typical hydromechanical coupling problem.

Because of the presence of numerous natural discontinuities in rock masses of shale reservoirs, complex hydraulic fracture networks instead of bi-wing fractures will be stimulated by fracturing operations [6]. The classical fracturing models (KGD, PKN, etc.) [7,8,9,10] are only applicable for modeling fractures with simple geometry. Fracturing experiments have been widely used for investigating the evolution mechanism of hydraulic fractures [11,12], but the evolution processes of fracture networks are still difficult to capture because of the limitations of rock sample sizes. Thus, investigating complex hydraulic fracture networks is challenging, and numerical methods have been regarded as essential tools. Commonly used numerical methods include the finite element method (FEM) [13], extended finite element method (XFEM) [14], displacement discontinuity method (DDM) [15], discrete element method (DEM) [16], etc. Some emerging numerical methods have drawn increasing attention, such as the finite–discrete element method (FDEM) [17], phase-field method [18], and peridynamic method [19]. In most existing fracturing models, the shale reservoir is typically treated as continuous and impermeable media. Under such an assumption, the effects of pore seepage on fracturing are neglected and the modeling efficiency could be improved. However, multi-scale natural fractures have been observed in shales [20], and it has been demonstrated that these fractures significantly enhance the permeability of rock masses [21]. Field applications have shown that the flowback of injected fluids is quite challenging, with a significant portion of the fluids remaining trapped within the rock formation [22]. This suggests the presence of highly permeable flow channels and fluid storage spaces within the reservoirs. Thus, although the intrinsic permeability of shale is relatively low, the effects of pore seepage on fracturing still cannot be neglected, or the stimulated volume may be overestimated.

While shale reservoirs should be regarded as fractured porous media rather than continuous and impermeable media, developing effective numerical models for simulating fracturing in such media remains a challenge. One numerical approach with the potential for simulating the propagation of complex fracture networks in fractured porous media is the discontinuous deformation analysis (DDA) method [23]. DDA is known for its theoretical robustness and has gained considerable attention in the field of discontinuous computing since its inception [24]. DDA employs an energy-based method to calculate the motion of independent blocks. Based on DDA, the discontinuous displacement of fracture surfaces can be easily captured, and no mathematical problems exist when solving the propagation and interaction of multiple fractures. Previous studies have demonstrated that DDA is highly effective for simulating hydraulic fracturing, particularly for the simulation of the formation of dense fracture networks [25,26,27]: Morgan and Aral [28,29] have attempted to simulate hydraulic fracture networks using DDA; Hu et al. [30,31] proposed an improved DDA-based fracturing model and investigated the formation mechanism of dense fracture networks in shale reservoirs. However, previous DDA-based models have neglected the effects of reservoir permeability on fracture propagation and fluid flow.

In general, the hydraulic stimulation of shale reservoirs represents a typical hydromechanical coupling problem, where the reservoir can be treated as fractured porous media. Numerical simulations serve as effective tools for investigating hydromechanical behaviors and optimizing fracturing operations. However, previous hydromechanical models often neglect the effects of pore seepage, which can significantly influence the simulation results. DDA shows great potential for simulating the fracturing processes in fractured porous media. In this work, a fully coupled hydromechanical model is developed based on DDA, which couples the fracture seepage, real pore seepage, and fracture network propagation. Several examples are presented to verify the model. In addition, comparative simulations are conducted to investigate the impact of rock mass permeability on hydraulic fracturing. Furthermore, critical discussions and future applications are also presented.

2. Brief Introduction to DDA Method

The discontinuous deformation analysis (DDA) method, which was originally proposed by Shi [23], is a fully implicit numerical method for simulating the discontinuous displacement of block systems. In DDA, there are no constraints on the shapes of blocks and every block has six independent variables:

$$\left[D_i\right]={\left(u_0{v}_0{r}_0{\epsilon }_x{\epsilon }_y{\epsilon }_{xy}\right)}^T,$$

(1)

where $u_0$ and $v_0$ represent the translation of the centroid $\left(x_0,y_0\right)$ along the x-axis and y-axis, respectively, while $r_0$ indicates the rotation of the rigid body around $\left(x_0,y_0\right)$. The variables ${\epsilon }_x$, ${\epsilon }_y,$ and ${\epsilon }_{xy}$ represent the three strain components of block i.

According to these six variables, the displacement $\left(u,v\right)$ of any point $\left(x,y\right)$ within block i can be obtained using a function of first-order complete approximation:

$$\left(\begin{array}{c}u\\ v\end{array}\right)=\left[T_i\right]\left[D_i\right],$$

(2)

$$\left[T_i\right]=\left(\begin{array}{cc}\begin{array}{ccc}1& 0& -(y-y_0)\\ 0& 1& (x-x_0)\end{array}& \begin{array}{ccc}(x-x_0)& 0& (y-y_0)/2\\ 0& (y-y_0)& (x-x_0)/2\end{array}\end{array}\right),$$

(3)

where $\left[T_i\right]$ represents the displacement transformation matrix. To establish the equilibrium equations, DDA uses an energy-based principle to build a set of symmetric equations, which have the following form:

$$\left[K\right]\left[D\right]=\left[F\right],$$

(4)

where vector $\left[D\right]$ (6n × 1) denotes the total displacement variables of block systems, vector $\left[F\right]$ (6n × 1) represents the global force vector (6n × 1), matrix $\left[K\right]$ (6n × 6n) stands for the global stiffness matrix, and the variable n represents the number of blocks. $\left[K\right]$ is a sparse symmetric matrix, which can be calculated according to the material properties of blocks and contact relationships between blocks. According to Equation (4), the displacement variables of each block can be solved and obtained.

Based on the above equations, it can be inferred that DDA shares significant similarities with the traditional FEM. Consequently, DDA is equally rigorous and accurate in solving deformations and displacements, while DDA can also simulate discontinuous deformations.

3. Hydromechanical DDA Model

3.1. Basic Assumption

Hydraulic fracturing in fractured porous rock masses involves several physical processes, including rock deformation, fracture propagation, fluid seepage in fractures, fluid seepage in the rock matrix, and fluid–solid interaction. To simulate these processes based on DDA, it is assumed that (1) the matrix of rock masses is homogeneous, isotropic, and linear elastic; (2) fractures can only propagate along the block interfaces; and (3) fluid can flow in fractures and the rock matrix, exhibiting laminar flow behavior.

3.2. Fracture Initiation

The fracture initiation is simulated under the framework of DDA. Efficient open–close iteration of block contacts is one of the most important characteristics of DDA, which enables DDA to simultaneously simulate both the continuous deformation and discontinuous displacement (such as separation, sliding, and rotation) of block systems. Assuming that fractures can only propagate along the boundaries of blocks, this model simulates the fracture initiation and propagation based on the block contact relationships.

Fracture initiation is described by the transformation of contact relationships between blocks: In a closed contact scenario, no fractures are present, and blocks undergo continuous deformation. When the contact shifts from closed to open, Type I tensile fractures emerge between blocks, allowing for surface separation and sliding. Similarly, when the contact shifts from closed to sliding, Type II shear fractures emerge between blocks, allowing for surfaces sliding. Additionally, this model assumes that once fractures occur, the corresponding tensile strength and cohesive strength of block interfaces become zero, maintaining a minimum aperture of $a_0$.

The assessment of tensile failure of fractures is based on the maximum tensile strength criterion, while the assessment of shear failure of fractures is based on the Mohr–Coulomb criterion:

$${\sigma }_n=T_0$$

(5)

$${\sigma }_{\tau }=c+{\sigma }_n\tan \phi$$

(6)

where ${\sigma }_n$ denotes the normal stress exerted on the fracture surfaces, ${\sigma }_{\tau }$ denotes shear stress acting on fracture surfaces, $T_0$ represents the tensile strength, $c$ denotes the cohesive strength, and $\phi$ denotes the friction angle.

3.3. Fluid Flow in Fractured Porous Media

Fluid flow in fractured and porous rock involves seepage through fractures and the rock matrix, as well as fluid transfer between the rock matrix and fractures. Figure 1 illustrates the fluid flow network model used to simulate these complex seepage processes. The network model consists of fluid cavities and connected flow channels, including block cavities (1, 2, 3, 4) and fracture cavities (i, j, k). The flow channels can be classified into fracture seepage channels, pore seepage channels, and pore–fracture seepage channels based on the properties of the two connected cavities.

In this model, the volumetric flow rate $q$ within hydraulic fractures (ij, jk) obeys the cubic law:

$$q=-{\displaystyle \frac{{\omega }^3}{12\mu }}{\displaystyle \frac{\partial p_f}{\partial l}},$$

(7)

In this equation, the symbol $\mu$ represents the viscosity of the fracturing fluid, $p_f$ represents the fluid pressure within fractures, and $l$ represents the coordinates along the length of the fracture. $\omega$ represents the fracture aperture, which is calculated based on the relative normal displacements between blocks:

$${\omega }_{ij}=e_m{\left[{\displaystyle \frac{16r^2}{{\left(1+r\right)}^4}}\right]}^{{\displaystyle \frac{1}{3}}},$$

(8)

where ${\omega }_{ij}$ denotes the fracture aperture of channel ij. $e_m$ and $r$ are calculated by

$$e_m={\displaystyle \frac{e_i+e_j}{2}},$$

(9)

$$r={\displaystyle \frac{e_i}{e_j}}$$

(10)

where $e_i$ is the relative normal displacement at point i, and $e_j$ is the relative normal displacement at point j.

Given the assumption of incompressible fluid, the following equation can be used to express the conservation of fluid mass within the fractures:

$${\displaystyle \frac{\partial \omega }{\partial t}}+{\displaystyle \frac{\partial q}{\partial l}}+c_0=0,$$

(11)

where $t$ represents the time, and $c_0$ denotes the fluid injection rate.

Along the flow channel between blocks (1→2, 3→4), the fluid flow rate obeys Darcy’s law:

$$v_w=-{\displaystyle \frac{k}{\mu }}{\displaystyle \frac{\partial p_p}{\partial l}},$$

(12)

where $k$ denotes the rock mass permeability, $p_p$ is the pore pressure within rock masses. Similar to fluid flow in fracture seepage channels, fluid flow in pore seepage channels also follows the principle of mass conservation:

$${\displaystyle \frac{\partial \varnothing }{\partial t}}+{\displaystyle \frac{\partial \left(\varnothing v_w\right)}{\partial l}}=0,$$

(13)

where $\varnothing$ indicates the porosity of rock masses.

In the developed model, as shown in Figure 1, when the rock masses are initially continuous, there exists a pore seepage channel between every two blocks. Once the interface between blocks breaks, a fracture seepage channel is added along the interface, and a virtual node is added at the block interface, forming two pore–fracture seepage channels (o1, o3) from the fracture surfaces to the adjacent blocks. The pressure continuity across the fracture determines the fluid flow from the fracture to the surrounding porous rock matrix, which can be written as

$$p_{+}=p_f=p_{-},$$

(14)

where $p_f$ denotes the fluid pressure at the midpoint of the fracture channel, and $p_{+}$ and $p_{-}$ represent the pressures in the upper (+) and lower (−) surrounding medium, respectively.

According to the above governing equations, an iteration-based finite volume scheme [32] is employed to solve the fluid pressure distribution in porous and discontinuous rock masses.

3.4. Hydromechanical Coupling

During hydraulic fracturing, fluid pressure acting on fractures and the rock matrix will cause stress changes in rock masses, thus resulting in fracture propagation. It is essential to consider the influence of fluid pressure in the hydromechanical model. Here, we assume that within each fracture ij, the fluid pressure is distributed along a linear gradient. The potential energy induced by the fluid pressure within fractures is minimized, and a submatrix (6 × 1) can be obtained as follows:

$$f_r={\int }_0^1{\left[T_i\right]}^T\left(\begin{array}{c}{p}_x\left(\xi \right)\\ p_y\left(\xi \right)\end{array}\right)md\xi ⟶\left[F_i\right] r=1,\cdots ,6,$$

(15)

where $\xi$ represents the local coordinates of fracture ij, $0\le \xi \le 1$, $m$ denotes the fracture length. This submatrix is incorporated into the force vector of block i. Based on the above processes, the fluid pressure within fractures is implemented into the DDA calculations. Then, the consideration of pore pressure within rock masses will be introduced.

Based on the effective stress concept, the effective stress of block i could be mathematically represented as

$${\sigma }_{ij}^{\prime }={\sigma }_{ij}-\alpha P_p{\delta }_{ij}$$

(16)

where ${\sigma }_{ij}$ is the stress component of block i, $\alpha$ represents the Biot coefficient, $p_p$ is the pore pressure within block i, ${\delta }_{ij}$ denotes the Kronecker symbol, when $i=j$, ${\delta }_{ij}=1$; when $i\ne j$, ${\delta }_{ij}=0$. Hence, the stress increment induced by the pore pressure is articulated as follows:

$${∆\sigma }_{ij}=-\alpha P_p{\delta }_{ij}$$

(17)

Here, we assume that the pore pressure within each block is constant. The potential energy induced by the pore pressure is minimized, and a submatrix (6 × 1) can be obtained as follows:

$$f_r=-S\left[{∆\sigma }_{ij}\right]⟶\left[F_i\right] r=1,\cdots ,6,$$

(18)

where $S$ represents the area of block i. This submatrix is then added to the force vector of block i.

Based on the above methodologies, the influence of fluid pressure on rock mass deformation and fracture propagation is implemented in the model. Furthermore, the volume of fluid cavities and the fracture permeabilities are dynamically updated according to the rock mass deformation, thereby influencing the fluid pressure distribution.

An iteration scheme is developed to achieve the hydromechanical coupling, which is shown in Figure 2. Firstly, the fluid pressure distribution within the rock matrix and fractures is calculated. Then, the fluid pressures within fractures and the matrix are incorporated into the global force vector $\left[F\right]$, and DDA calculations are performed to obtain the total deformation variables $\left[D\right]$ of the block systems. Then, the contact relationships between blocks are checked to determine whether new fractures are generated. At the same time, the fracture aperture and the block volume are calculated to update the fluid cavity volume and fracture permeability. In addition, the fluid flow network including fracture seepage channels and pore seepage channels is also updated according to the fracture network. This iterative procedure continues until reaching the final time step.

4. Model Verification

In this part, several examples will be presented to verify the proposed hydromechanical DDA model. Since our previous research has already demonstrated the accuracy of our earlier model in addressing fracture deformation, fracture seepage, and fracture propagation [30,31], we will not include similar verifications in this paper. Instead, two new examples are provided here to illustrate the model’s capability of simulating fluid flow within fractured porous media.

4.1. Unsteady Pore Seepage

To assess the accuracy of the hydromechanical DDA model in modeling unsteady pore seepage, a verification example by Yan et al. [33] is employed here. The model is shown in Figure 3a. The length L is 10 m and the width is 1 m. The rock matrix is considered an isotropic porous medium. Constant pressure boundary conditions are applied to the model, with fluid pressure set to $p_1=10 \mathrm{M}\mathrm{P}\mathrm{a}$ on the left side and $p_2=0 \mathrm{M}\mathrm{P}\mathrm{a}$ on the right side. The analytical solution to this problem is mathematically expressed as follows:

$$p\left(x,t\right)=p_1+{\displaystyle \frac{x}{L}}\left(p_2-p_1\right)+{\displaystyle \frac{2}{\pi }}{\sum }_{n=1}^{\infty }{e}^{-\kappa n^2{\pi }^{2}t/L^2}\left({\displaystyle \frac{p_2\cos \left(n\pi \right)-p_1}{n}}\right)\sin {\displaystyle \frac{n\pi x}{L}},$$

(19)

where $t$ represents time and $x$ represents the coordinates. The parameter $\kappa$ is defined as $\kappa =\left(k/\mu \right)M,$ where $k$ is the rock permeability, $\mu$ is the dynamic viscosity, and $M$ represents the Biot modulus. $M$ is calculated as $M=K_w/n,$ where $K_w$ is the bulk modulus of the fluid and n is the rock porosity. Figure 3b illustrates the distribution of pore pressure ($t=0.1, 1, 3, 5 \mathrm{s}$), with $K_w=2.2 \mathrm{G}\mathrm{P}\mathrm{a}$, $\mu =1 \mathrm{m}\mathrm{P}\mathrm{a}·\mathrm{s},n=0.1$. The numerical and analytical results match well, demonstrating the accuracy of the proposed DDA model in modeling the unsteady pore seepage within the porous media.

4.2. Mixed Pore–Fracture Seepage

To examine the developed model in simulating mixed pore–fracture seepage, a standard example by Tatomir [34] is employed here. The model is shown in Figure 4a. The upper boundary is defined as a pressure boundary, and all other boundaries are impermeable. The dimensions of the model and pressure boundary settings, as well as the permeability value, are illustrated in Figure 4a. The model mesh is presented in Figure 4b.

The fluid pressure distribution obtained from the numerical simulations is depicted in Figure 5a. We compare the simulated hydraulic head $h$ ($h=p/\rho g$, where $\rho$ is water density, $\rho =1.0\times {10}^3 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$, and $g$ is gravity acceleration, $g=9.8 \mathrm{m}/{\mathrm{s}}^2$) at different elevations ($y=\mathrm{600,800,1000} \mathrm{m}$) with the results of Tatomir [34]. As depicted in Figure 5b–d, the numerical results align closely with those of the standard example. Thus, this example demonstrates that the proposed model can accurately simulate the fluid pressure distribution in fractured porous media.

5. Hydraulic Fracturing in Fractured Porous Rock Masses

This section employs the developed hydromechanical DDA model to explore the evolution processes of hydraulic fracture networks in fractured porous rock masses.

Figure 6 displays a typical model of porous and fractured rock masses, which was built based on the features of the Silurian Longmaxi shale formations in the Sichuan Basin, China [30]. The Longmaxi shale is usually regarded as a dense rock with relatively low permeability and small porosity. Its matrix permeability ranges from μD (${10}^{-18} {\mathrm{m}}^2$) to nD (${10}^{-21} {\mathrm{m}}^2$) based on permeability test results [21]. Previous studies on the fracturing of the Longmaxi shale reservoirs usually neglect the effect of pore seepage. However, numerous natural discontinuities such as beddings, high-angle fractures, and multi-scale laminations have been observed in the Longmaxi shale. The experimental results suggest that the presence of a single fracture can increase the macroscopic permeability of shale to ${10}^2$–${10}^3$μD, which is significantly higher than the permeability of the matrix itself [21]. Thus, it is necessary to consider the influence of these natural fractures on the permeability of rock masses. In deep shale reservoirs, the natural discontinuities close under in situ stresses, but highly permeable channels still exist, and the macroscopic permeability of the reservoir remains significantly larger than the matrix permeability. Previous simulation studies of hydraulic fracturing in the Longmaxi shale have not given sufficient attention to the influence of rock mass permeability. In addition, the inclusion of rock mass permeability increases the complexity of simulating fracturing processes. The hydromechanical DDA model developed in this paper can account for the effect of the permeability, making it suitable for investigating the impact of rock mass permeability on the Longmaxi shale reservoirs.

The model shown in Figure 6 represents a vertical section of the rock masses in the Longmaxi shale formation. The model considers both natural fractures and bedding planes, while the rock matrix is regarded as a porous and permeable medium. In this model, the natural fractures and bedding planes are reconstructed based on statistical measurement data obtained from an outcrop of the Longmaxi shale formation in the Sichuan Basin [30]. According to outcrop observations, the natural discontinuities in the Longmaxi shale reservoir mainly consist of one set of bedding planes and one set of high-angle fractures, which are approximately perpendicular to the bedding planes. In the reconstructed model illustrated in Figure 6, the bedding planes are uniformly distributed with a dip angle of 30°, and the average trace length of high-angle natural fractures is 15 m. To consider the influence of discontinuities on fracturing, natural fractures and bedding planes are initially treated as weak planes with very low strengths and seepage paths with relatively high permeabilities. Fracturing fluid can permeate into the surrounding rock masses, which may affect the fracturing processes and fracture configurations. The model was meshed using triangular elements with fine refinement in simulations. Fracturing fluid is injected into the center of the rock mass consistently.

Table 1. Parameters used for simulations of fracturing in fractured porous rock masses.

Parameters	Value	Parameters	Value
Elastic modulus, $E$	$38 \mathrm{G}\mathrm{P}\mathrm{a}$	Fluid bulk modulus, $K_w$	$2.2 \mathrm{G}\mathrm{P}\mathrm{a}$
Poisson ratio, $v$	0.20	Fluid viscosity, $\mu$	$1 \mathrm{m}\mathrm{P}\mathrm{a}·\mathrm{s}$
Tensile strength, $T_0$	$3 \mathrm{M}\mathrm{P}\mathrm{a}$	Fluid density, ${\rho }_f$	$1000 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
Cohesive strength, $c$	$10 \mathrm{M}\mathrm{P}\mathrm{a}$	Injection rate, $Q_0$	$0.005 {\mathrm{m}}^3/\mathrm{m}/\mathrm{s}$
Friction angle, $\phi$	45°	Injection volume, $V$	$0.5 {\mathrm{m}}^3/\mathrm{m}$
Porosity, $n$	0.1	${\sigma }_y$	$27 \mathrm{M}\mathrm{P}\mathrm{a}$
Permeability, $k$	$0.5 \mathrm{m}\mathrm{D}$	${\sigma }_x$	$30 \mathrm{M}\mathrm{P}\mathrm{a}$

presents the mechanical parameters of the Longmaxi shale and the parameters used in the fracturing process simulations.

To explore the influence of rock permeability on fracturing processes, various values of rock permeability k ($0 \mathrm{m}\mathrm{D}$, $0.01 \mathrm{m}\mathrm{D}$, $0.1 \mathrm{m}\mathrm{D}$, $1 \mathrm{m}\mathrm{D}$) were employed in the simulations, with all other parameters held constant. The simulation results of the fracture network configurations are shown in Figure 7. It is shown that with the increase in rock permeability, the stimulated reservoir area is significantly reduced and the fracture network configurations become simpler. To obtain a more in-depth analysis of the impact of permeabilities, the total length and maximum aperture of stimulated fractures are statistically measured and compared. As illustrated in Figure 8, the length of the stimulated hydraulic fractures decreases roughly linearly with the increasing permeability, while the width of the fractures exhibits a non-linear decreasing trend. In other words, the increase in permeability is unfavorable for the formation of hydraulic fracture networks.

As previously discussed, higher rock permeability induces significant alterations in the hydraulic fracture network configurations. To further explore the extension processes of hydraulic fracture networks, an in-depth analysis is undertaken, focusing on the permeability condition of $k=0.50 \mathrm{m}\mathrm{D}$. Figure 9 illustrates the fluid pressure evolution processes of hydraulic fracture networks during fluid injection, while Figure 10 illustrates the fracture aperture evolution processes. In addition, Figure 11 delineates the relation between the fracture length and fluid injection volume under the condition of $k=0.50 \mathrm{m}\mathrm{D}$, alongside a comparative depiction of the fracture length evolution at $k=0.01 \mathrm{m}\mathrm{D}$. Generally, hydraulic fracture networks expand rapidly during the initial phase of fluid injection, then the expansion decelerates notably. Specifically, the main fracture initially propagates along the adjacent bedding planes towards the injection point, extending rapidly during the initial injection phase. At this stage, most of the injected fluids are contained within the main fracture, resulting in a significant increase in the fracture aperture due to high fluid pressure. Subsequently, surrounding natural fractures are activated, and the main fracture gradually reorients towards the horizontal direction, aligning with the orientation of the maximum principal stress. While the main fracture length increases, the aperture of the main fractures reaches its maximum value, and the aperture of the activated natural fractures remains relatively small. The overall morphology of the fracture network conforms to the orientation of bedding planes, diverging from the morphology observed at $k=0.01 \mathrm{m}\mathrm{D}$ (Figure 7b). According to the above simulation results, it can be deduced that because of the presence of relatively high permeability, fractures tend to propagate along preferential structural planes, thereby diminishing the influence of in situ stress.

The presence of pores and micro-fractures may cause fracturing fluids to leak into rock masses, hindering the formation of hydraulic fracture networks. Figure 10 illustrates the variation in the fluid volume within fractures and the fluid leakoff volume. When $k=0.01 \mathrm{m}\mathrm{D}$, the fluid leakoff volume is minimal, and most of the injected fluid volume remains within the fractures, resulting in a higher efficiency of hydraulic fracturing. Significant fluid leakoff is observed when $k=0.50 \mathrm{m}\mathrm{D}$. The length of the stimulated fractures gradually increases with the increase in the fluid injection volume (Figure 11), resulting in an escalation in the fluid leakoff rate (Figure 12). As fluid injection continues, the extension rate of hydraulic fractures decreases. Therefore, fluid leakoff becomes the primary factor governing the expansion of stimulated fracture networks in the later stages of fracturing.

According to the above numerical results, it can be concluded that rock permeability plays a critical role in the fracturing processes. In reservoirs characterized by higher permeability, it is vital to reduce the fluid leakoff volume during the fracturing process to enhance the overall fracturing efficiency. Improving the viscosity of fracturing fluids is commonly used as a strategy to decrease the volume of fluid leakoff. Here, the effectiveness of augmenting the fluid viscosity in minimizing fluid leakoff is investigated through numerical simulations. Figure 13 shows the simulated fracture network configurations under different viscosities of fracturing fluids, while Figure 14 illustrates the variation in the stimulated fracture length and fluid leakoff volume. It can be observed that when the fluid viscosity is improved, the volume of fluid leakoff significantly decreases (Figure 14b), while the stimulated fracture length obviously increases (Figure 14a). Furthermore, as shown in Figure 13, the improvement in the fluid viscosity enhances the fluid net pressure inside the fractures, leading to a more complex fracture network morphology and a larger stimulated reservoir area. Hence, improving the fracturing fluid viscosity is useful in reducing the fluid leakoff volume and increasing the fracturing efficiency.

6. Discussions

In this work, a new fully coupled DDA model is developed for simulating the hydromechanical behavior of fractured porous media. As a discontinuum-based numerical method, DDA holds great potential for solving complex discontinuous geo-mechanical problems including hydraulic fracturing. Several DDA-based models have been previously developed to simulate the hydromechanical coupling processes in fractured media [25,26,27,28,29,30,31]. However, previous models mainly focus on the fluid seepage within fractures and their effects on rock deformation or fracture propagation; the fluid seepage within the rock matrix and fluid transfer between fractures and the matrix are always neglected. As fluid leakoff is a vital process for fracturing processes, some models used Carter’s model to consider the effects of leakoff:

$$u={\displaystyle \frac{C}{\sqrt{t}}},$$

(20)

where $u$ is the leakoff rate, $C$ is the leakoff coefficient, and $t^{\prime }$ is the time. Since Carter’s leakoff model is derived originally from Darcy’s law, $C$ can be written with

$$C={\left({\displaystyle \frac{k∆p\phi }{\mu }}\right)}^{0.5},$$

(21)

where $k$ is the rock permeability, $\phi$ is the rock porosity, $\mu$ is the viscosity of fluid, and $\mathsf{\Delta }p$ represents the fluid pressure differences between the fracture and the matrix. According to the above two equations, the fluid leakoff rate is not only related to the intrinsic rock parameters $\left(k,\phi \right)$ and fluid viscosity $\mu$, but also related to the fluid pressure distribution in porous rock masses. Thus, the classical Carter’s leakoff model with a constant coefficient C cannot reflect the effects of pressure evolution, which may cause errors in calculating the fluid flux and fluid pressure. In this work, the fluid transfer between fractures and the rock matrix is considered by simulating the real seepage processes instead of using Carter’s model. In addition, the fluid pressure evolution in the rock matrix can also be obtained by this model, and the effect of pore pressure on rock deformation is also incorporated into the coupled processes. Therefore, compared to previous DDA-based models, this work shows a significant improvement in the consideration of real pore seepage and coupling processes in fractured porous media.

Based on the newly developed DDA model, the authors investigated the fracturing processes in a fractured porous shale reservoir. The fracture network configurations under different reservoir permeabilities were obtained, and the fluid leakoff processes were analyzed. The comprehensive simulation results revealed several notable findings. Firstly, although the intrinsic permeability of shale is relatively low, the effects of pore seepage on the fracturing processes cannot be overlooked due to the presence of natural discontinuities within rock masses. As the reservoir permeability increases, both the stimulated fracture length and fracture aperture decrease significantly. Secondly, the results show that most fracturing fluids leak off into the rock masses, with only a small amount remaining in activated fractures, which explains the difficulties associated with the flowback of fracturing fluids observed in practical operations. Thirdly, the leakoff rate increases with the extension of the hydraulic fracture network until the fluid injection and leakoff rates reach a balance, at which point the fracture networks cease to propagate. In other words, reservoir permeability and fluid leakoff effects are the primary factors influencing fracture propagation in the later stages of fracturing. These findings may be helpful for the optimization of fracturing operations. It is essential to use rock mass permeability rather than intrinsic permeability in the design of fracturing operations; otherwise, the stimulated reservoir volume may be overestimated. To reduce the fluid leakoff volume and enhance the fracturing efficiency, increasing the viscosity of the fracturing fluid is beneficial. However, it is believed that high-viscosity fracturing fluids may hinder the activation of natural fractures due to increased flow resistance. Therefore, using variable viscosity fracturing fluids may provide an alternative solution. In the early stages of fracturing, relatively low viscosity fluids facilitate fluid flow and the formation of new fractures. Subsequently, the viscosity of the fracturing fluids should be gradually increased to slow the leakoff rate. This increase in flow resistance will enhance fracture apertures, allowing proppants to enter fractures more easily, while the injected low-viscosity fluids will penetrate the rock masses, promoting the activation of natural fractures.

In addition to considering the effects of fluid leakoff, another objective of integrating real pore seepage into the DDA model is to equivalently simulate the evolution of sub-scale fracture networks. Multi-scale natural fractures observed in shale will consequently activate similar multi-scale hydraulic fractures during fluid injection, which are quite essential for enhancing reservoir permeability. Previous studies [35] have shown that natural fractures at different scales exhibit similar geometric characteristics. Specifically, their occurrence, defined by the strike, dip angle, and dip direction, is fundamentally consistent, while fracture spacing and trace length display fractal characteristics. However, due to the vast number of multi-scale hydraulic fractures, it is impractical to consider all the fractures in a single model. Therefore, when modeling reservoir-scale main fractures, it is common practice to account for the evolution of sub-scale fractures using an equivalent approach. Our previous research [36] has numerically investigated the permeability of hydraulic fracture networks at the laboratory scale, yielding permeability values under varying fluid pressures and in situ stresses. The results indicated a strong anisotropy in permeability, with the permeability curves displaying similar shapes. This study primarily focuses on the laboratory-scale fracturing of the Silurian Longmaxi shale. From a macro perspective, the propagation of hydraulic fracture networks can be treated as the dynamic evolution of the anisotropic permeability of the rock masses. By incorporating dynamic anisotropic seepage into the reservoir-scale fracturing model, the impacts of multi-scale fracture networks can be assessed, allowing for the determination of fracture network conductivity distribution within the reservoir. Additionally, fitting equations for anisotropic permeability under different stresses and fluid pressures have been provided for convenient application in our prior research. Thus, building upon our previous findings regarding fracture network permeability and utilizing the newly developed porous DDA model from this work, the equivalent evolution processes of multi-scale fracture networks can be simulated without encountering numerical challenges associated with numerous fractures.

It should be noted that the newly developed hydromechanical DDA model in this work can also be used in various geo-mechanical problems, such as the rock slope induced by rainfall, dam foundation stability analysis, rock cavern stability analysis under hydrogeological conditions, geological CO₂ storage, etc. In addition, previous studies have demonstrated that thermal effects significantly impact rock mass deformation and fluid seepage in fractured porous media, thereby influencing fracture propagation processes [37]. In addition, chemical reaction processes are also important in some geo-mechanical problems. Thus, to meet the requirement of rising geo-mechanical problems, our future work also includes implementing heat transfer, heat conduction, and chemical reaction processes into this hydromechanical DDA model, which enables the development of a thermo-hydromechanical–chemical coupling DDA model.

7. Conclusions

This work developed a fully coupled DDA model for simulating the hydromechanical processes in fractured porous media. Compared to previous DDA-based models, this work shows a significant improvement in the consideration of real pore seepage and coupling processes in fractured porous media. Several verification examples were given to show the reliability of the developed model.

Comparative simulations were conducted to investigate the influence of rock permeabilities on the fracturing processes. The numerical results revealed that as the rock mass permeability increases, the fluid leakoff volume increases, resulting in a smaller length and aperture of hydraulic fractures. In highly permeable rock masses, hydraulic fractures exhibit a propensity to propagate along weak discontinuities, making it difficult to form new fractures and complex fracture networks. In addition, increasing the fracturing fluid viscosity is useful for reducing the fluid leakoff volume and improving the fracturing efficiency.

This study presented a viable numerical approach for examining hydromechanical processes in rock masses characterized by natural discontinuities and high permeabilities, which can be used in fracturing optimization designs and various geo-mechanical problems.