Hydro-Mechanical Coupled FDEM Modeling of Multi-Cluster Hydraulic Fracture Propagation and Interaction in Shale Reservoirs

1. Introduction

Hydraulic fracturing is a key stimulation process for improving fluid access to low-porosity and low-permeability unconventional reservoirs. Its effectiveness depends not only on creating individual fractures but also on controlling the geometry, connectivity, and interaction of multiple fractures initiated from closely spaced clusters. During treatment, fracture growth is governed by the coupled effects of in-situ stress, rock deformation, fluid pressure, leak-off, fracture interaction, and operational parameters. Consequently, understanding the mechanisms that control fracture deflection, branching, aperture, and stress-shadow interference is central to improving stimulated reservoir volume (SRV) and treatment efficiency [1,2,3,4,5,6,7,8].

Field observations and previous studies have shown that in-situ stress is a primary control on hydraulic-fracture orientation and extent [9,10,11]. Analytical models such as the PKN, KGD, and penny-shaped formulations provide important benchmarks for idealized fracture growth [12,13,14,15,16,17,18,19,20,21,22], but their simplifying assumptions limit their ability to represent branching, coalescence, contact, and multi-fracture competition. Conventional finite-element approaches can also require remeshing as cracks propagate, while enriched or discontinuum-based approaches are better suited to evolving fracture paths [23,24,25]. For shale stimulation, a numerical framework should therefore resolve both the continuous deformation of the rock matrix and the discontinuous evolution of fractures while maintaining fluid-solid coupling.

The finite-discrete element method (FDEM) provides such a framework by combining finite-element deformation with discrete fracture and contact mechanics. Zero-thickness cohesive interfaces can represent fracture initiation and progressive damage, whereas contact algorithms describe the mechanical response after interface failure. When coupled with matrix seepage and fracture flow, FDEM can explicitly represent the transition from intact porous rock to an evolving hydraulic-fracture network. This capability is particularly useful for investigating multi-cluster stimulation, where neighboring fractures interact through stress shadows and compete for fluid and propagation space.

Accordingly, this study develops a two-dimensional hydro-mechanical FDEM model for multi-cluster hydraulic fracturing in shale. The framework couples Biot poroelasticity, Darcy seepage, cubic-law fracture flow, Carter-type leak-off, cohesive-zone damage, and post-failure contact. The model is first verified against the classical KGD solution. It is then used to quantify the effects of four process variables—horizontal principal stress difference, cluster spacing, injection rate, and fracturing-fluid viscosity—on fracture morphology and propagation. Finally, the individual sensitivities are integrated into a combined stimulation strategy. The principal contribution is a unified process-level analysis linking controllable treatment parameters to stress-shadow interaction, fracture opening, and branching behavior within one verified FDEM framework.

The remainder of the paper is organized as follows. Section 2 presents the governing equations, constitutive assumptions, fracture-damage formulation, and hydro-mechanical coupling strategy. Section 3 describes model verification and the parametric simulations. Section 4 discusses the coupled mechanisms, engineering implications, and limitations of the model, and Section 5 summarizes the main conclusions.

2. Materials and Methods

2.1. Modeling Assumptions

A coupled hydro-mechanical FDEM framework is established to simulate hydraulic-fracture initiation and propagation in saturated shale. The following assumptions are introduced to retain the dominant physical processes while maintaining computational tractability:

(1) The intact rock matrix is simplified as fully saturated linear elastic porous media, and fluid migration within internal micropores complies with the classical Darcy seepage theory;

(2) Propagating hydraulic fractures are idealized as smooth parallel-plate flow channels, where low-Reynolds laminar transport of injected fracturing fluid follows the standard cubic flow law;

(3) Rock deformation generated during hydraulic fracturing is categorized as infinitesimal strain deformation, and the initiation and progressive expansion of fractures are governed by the irreversible damage evolution criterion of the cohesive zone model;

(4) The fracturing fluid adopted in the simulation is defined as incompressible Newtonian fluid, and fluid percolation from fracture surfaces toward the surrounding rock matrix is incorporated into the fluid mass conservation equation;

(5) A bidirectional fluid-solid coupling algorithm is developed based on Biot poroelastic theory, and an explicit time-stepping integration scheme is utilized for iterative dynamic solving of the multi-field coupled system.

2.2. Solid Mechanical Governing Equations

The FDEM numerical technique discretizes fractured rock media into continuous finite element units to characterize matrix deformation responses. Zero-thickness cohesive interface units are embedded at the boundaries of adjacent rock blocks to capture the whole process of hydraulic fracture initiation and progressive expansion. Meanwhile, discrete element contact mechanics are coupled to quantify the discontinuous mechanical characteristics of fractured rock, including interface extrusion, relative sliding and post-failure block displacement.

2.2.1. Dynamic Equilibrium Equation of Solid Medium

In accordance with Newton’s second law of motion, a global dynamic equilibrium formula considering structural damping is constructed to describe the transient mechanical response of rock mass nodes under fracturing disturbance.

$$\begin{array}{c}\mathit{M}\ddot{\mathit{u}}+\mathit{C}\dot{\mathit{u}}={\mathit{F}}_{\mathrm{total}}\end{array}$$

(1)

where M denotes the nodal lumped mass matrix, C represents the global damping matrix of the rock structural system, and $\ddot{\mathit{u}}$, $\dot{\mathit{u}}$, u refer to nodal acceleration, velocity and displacement vectors respectively.

The total external load acting on rock blocks can be decomposed into five independent mechanical components:

$$\begin{array}{c}{\mathit{F}}_{\mathrm{total}}={\mathit{F}}_{\mathrm{el}}+{\mathit{F}}_{\mathrm{co}}+{\mathit{F}}_{\mathrm{ct}}+{\mathit{F}}_{\mathrm{fl}}+{\mathit{F}}_{\mathrm{g}}\end{array}$$

(2)

Specifically, F_el is the elastic internal force induced by finite element unit deformation; F_co corresponds to the bonding constraint force of damageable cohesive interfaces; F_ct is the frictional contact force between discrete rock blocks after fracture opening; F_fl stands for the equivalent load generated by pore and fracture fluid pressure; F_g is the body load induced by rock gravity acceleration.

2.2.2. Biot Poroelastic Constitutive Relation

To quantitatively characterize the fluid-solid interaction effect of saturated porous rock, the Biot poroelastic effective stress principle is adopted to establish the correlation between rock mechanical state and pore fluid pressure. The total stress tensor is superposed by the effective stress borne by rock skeleton and the additional stress induced by pore fluid pressure:

$$\begin{array}{c}\mathit{\sigma }={\mathit{\sigma }}^{\text{'}}+\alpha p_{\mathrm{m}}\mathit{I}\end{array}$$

(3)

where σ is the macroscopic total stress tensor of rock mass, σ^' is the effective stress carried by rock solid skeleton, α is the Biot effective stress coefficient, p_m represents matrix pore fluid pressure, and I is the second-order unit tensor.

Under infinitesimal deformation conditions, the rock skeleton exhibits linear elastic mechanical properties, and its constitutive relationship is expressed as:

$$\begin{array}{c}{\mathit{\sigma }}^{\text{'}}=\mathbb{D}:\mathit{\epsilon }\end{array}$$

(4)

where 𝔻 refers to the fourth-order elastic stiffness tensor of rock materials, and ε is the infinitesimal strain tensor calculated from displacement gradient:

$$\begin{array}{c}\mathit{\epsilon }={\displaystyle \frac{1}{2}}[\nabla \mathit{u}+(\nabla \mathit{u}{)}^{\mathrm{T}}]\end{array}$$

(5)

2.2.3. Cohesive Zone Damage and Fracture Evolution Model

Zero-thickness cohesive units are pre-laid along potential fracture propagation paths to reproduce the physical fracture process zone formed during hydraulic fracturing. A bilinear traction-separation constitutive model is employed to simulate the full-range mechanical evolution of rock interfaces, covering initial damage generation, progressive stiffness degradation and ultimate structural failure.

The mechanical response of cohesive interfaces is dominated by the quantitative relationship between interface traction and relative displacement:

$$\begin{array}{c}{\sigma }_{\mathrm{n}}=f\left({\delta }_{\mathrm{n}}\right),{\tau }_{\mathrm{s}}=f\left({\delta }_{\mathrm{s}}\right)\end{array}$$

(6)

where σ_n and τ_s denote normal tensile traction and tangential shear traction on cohesive interfaces, while δ_n andδ_s represent normal opening displacement and tangential sliding displacement of interface units respectively.

A dimensionless irreversible damage variable varying from 0 to 1 is defined to quantitatively assess the progressive damage degree of rock bonding interfaces. The dynamic evolution formula of damage is presented as follows:

$$\begin{array}{c}D={\displaystyle \frac{{\delta }_{\max }-{\delta }_0}{{\delta }_{\mathrm{f}}-{\delta }_0}}\end{array}$$

(7)

where δ₀ is the critical displacement for triggering initial interface damage, δ_f is the ultimate displacement causing complete interface fracture, and δ_max is the maximum interface separation displacement recorded in the loading process, which ensures the irreversibility of damage accumulation. The interface remains intact at D = 0 and completely loses mechanical bearing capacity at D = 1.

A mixed-mode fracture criterion based on critical energy release rate is adopted to evaluate the dynamic propagation status of hydraulic fractures:

$$\begin{array}{c}{\left({\displaystyle \frac{G_{\mathrm{I}}}{G_{\mathrm{Ic}}}}\right)}^m+{\left({\displaystyle \frac{G_{\mathrm{II}}}{G_{\mathrm{IIc}}}}\right)}^m\ge 1\end{array}$$

(8)

where G_I and G_II are transient energy release rates corresponding to mode-I tensile fracture and mode-II shear fracture, G_Ic and G_IIc are critical fracture toughness values for tensile and shear failure, and m is the calibration exponent for the mixed fracture criterion.

2.2.4. Contact Mechanical Equations for Discrete Rock Blocks

After cohesive units fail completely and hydraulic fractures open, the contact interaction algorithm between adjacent discrete rock blocks is activated. Inter-block contact force is divided into normal compression component and tangential friction component to accurately simulate the post-fracturing movement and mechanical response of rock blocks.

The normal contact interaction between rock blocks conforms to the linear stiffness mechanical model:

$$\begin{array}{c}{F}_{\mathrm{n}}=k_{\mathrm{n}}{\mathsf{\Delta }}_{\mathrm{n}}\end{array}$$

(9)

where k_n is the normal contact stiffness of block contact surfaces, and Δ_n represents the normal embedding depth between adjacent contact blocks.

The tangential frictional contact behavior obeys the classical Mohr-Coulomb shear strength criterion:

$$\begin{array}{c}{F}_{\mathrm{s}}=\mathit{min}(k_{\mathrm{s}}{\mathsf{\Delta }}_{\mathrm{s}},F_{\mathrm{n}}\tan \phi +c{A}_{\mathrm{c}})\end{array}$$

(10)

where k_s denotes tangential contact stiffness, Δ_s is the relative tangential sliding displacement of contact surfaces, φ is the internal friction angle of fractured rock mass, c is the residual cohesion of fracture surfaces, and A_c is the effective contact area between discrete rock blocks.

2.3. Multi-Field Seepage Governing Equations

2.3.1. Seepage Equation for the Rock Matrix

Combining Darcy seepage theory and Biot poroelastic mass conservation principle, a transient seepage governing formula for saturated rock matrix is established. This equation fully considers the coupling effects of rock volumetric deformation and fracture fluid percolation on dynamic variation of pore pressure:

$$\begin{array}{c}{\displaystyle \frac{1}{M_{\mathrm{b}}}}{\displaystyle \frac{\partial p_{\mathrm{m}}}{\partial t}}=\nabla \cdot \left({\displaystyle \frac{k_{\mathrm{m}}}{{\mu }_{\mathrm{l}}}}\nabla p_{\mathrm{m}}\right)-\alpha {\displaystyle \frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}}-q_{\mathrm{lk}}\end{array}$$

(11)

where M_b is the Biot modulus characterizing pore volume compressibility, k_m is the inherent permeability of rock matrix, μ_l is the dynamic viscosity of fracturing fluid, ε_v is the volumetric strain of rock skeleton, and q_lk is the unit-volume fluid percolation flux from hydraulic fractures to surrounding matrix. Volumetric strain is defined as the divergence of nodal displacement vector:

$$\begin{array}{c}{\epsilon }_{\mathrm{v}}=\nabla \cdot \mathit{u}\end{array}$$

(12)

The Biot modulus is determined by the comprehensive compression characteristics of rock solid skeleton and pore filling fluid:

$$\begin{array}{c}{M}_{\mathrm{b}}={\left({\displaystyle \frac{\varphi }{K_{\mathrm{l}}}}+{\displaystyle \frac{\alpha -\varphi }{K_{\mathrm{s}}}}\right)}^{-1}\end{array}$$

(13)

where ϕ represents matrix porosity, K_l and K_s are the bulk modulus of fracturing fluid and rock solid skeleton respectively.

2.3.2. Flow Governing Equation for Hydraulic Fractures

Fluid flow inside hydraulic fractures is simplified as parallel-plate laminar flow, and the flow conductivity of fractures is governed by the aperture-dependent cubic law. The transient mass conservation equation for internal fracture flow incorporates wellbore fluid injection and fracture wall fluid percolation effects:

$$\begin{array}{c}{\displaystyle \frac{\partial w}{\partial t}}=\nabla \cdot \left({\displaystyle \frac{w^3}{12{\mu }_{\mathrm{l}}}}\nabla p_{\mathrm{f}}\right)-q_{\mathrm{total}}+Q_{\mathrm{inj}}\end{array}$$

(14)

where w is the dynamic aperture of propagating hydraulic fractures, which is numerically equivalent to the normal opening displacement of cohesive interface units; p_f is the fluid pressure within fracture cavities; w³/12μ_l denotes the dynamic flow conductivity of hydraulic fractures; q_total is the cumulative fluid percolation flux on fracture surfaces; Q_inj represents the fluid injection source term at wellbore perforation positions.

The average transport velocity of fracturing fluid along fracture propagation direction is calculated via the cubic flow law:

$$\begin{array}{c}{\mathit{v}}_{\mathrm{f}}=-{\displaystyle \frac{w^2}{12{\mu }_{\mathrm{l}}}}\nabla p_{\mathrm{f}}\end{array}$$

(15)

The seepage velocity of fluid in rock matrix pores complies with the classical Darcy velocity formula:

$$\begin{array}{c}{\mathit{v}}_{\mathrm{m}}=-{\displaystyle \frac{k_{\mathrm{m}}}{{\mu }_{\mathrm{l}}}}\nabla p_{\mathrm{m}}\end{array}$$

(16)

The time-dependent fluid percolation behavior from fractures to matrix is quantified by the widely applied Carter percolation theoretical model:

$$\begin{array}{c}{q}_{\mathrm{total}}={\displaystyle \frac{2C_{\mathrm{l}}}{\sqrt{t-t_0}}}\end{array}$$

(17)

where C_l is the Carter percolation coefficient correlated with rock physical properties, t is the real-time simulation duration, and t₀ corresponds to the initiation moment of each independent hydraulic fracture.

The Biot effective stress coefficient is calculated based on the modulus ratio of dry rock matrix and pure rock skeleton:

$$\begin{array}{c}\alpha =1-{\displaystyle \frac{K_{\mathrm{dry}}}{K_{\mathrm{s}}}}\end{array}$$

(18)

where K_dry denotes the bulk modulus of dry rock matrix without pore fluid filling.

The critical normal opening displacement for tensile damage initiation of cohesive units is determined by rock tensile strength and interface elastic stiffness:

$$\begin{array}{c}{\delta }_0^{\mathrm{n}}={\displaystyle \frac{{\sigma }_{\mathrm{t}}}{k_{\mathrm{coh}}}}\end{array}$$

(19)

where σ_t is the uniaxial tensile strength of rock materials, and k_coh is the initial elastic stiffness of cohesive interface units.

The critical tangential sliding displacement triggering shear damage of cohesive interfaces is defined by the peak shear strength of rock internal bonding surfaces:

$$\begin{array}{c}{\delta }_0^{\mathrm{s}}={\displaystyle \frac{{\tau }_{\max }}{k_{\mathrm{coh}}}}\end{array}$$

(20)

where τ_max represents the peak shear strength of rock bonding interfaces.

In numerical simulation procedures, the macroscopic aperture of hydraulic fractures is quantitatively characterized by the normal opening displacement of embedded cohesive units:

$$\begin{array}{c}w={\delta }_{\mathrm{n}}\end{array}$$

(21)

3. Results

This section first evaluates the numerical accuracy of the hydro-mechanical FDEM framework using a classical analytical benchmark and then examines the effects of horizontal principal stress difference, cluster spacing, injection rate, and fracturing-fluid viscosity on multi-cluster fracture propagation.

3.1. Model Verification Against the KGD Solution

The coupled fluid-flow and solid-deformation formulation was verified using a simplified single-fracture problem corresponding to the classical two-dimensional KGD configuration. An incompressible Newtonian fluid was injected into the fracture using the parameters defined in the original model setup. Fracture half-length and maximum aperture at the fracture center were selected as verification metrics because they directly reflect the coupled pressure–deformation response.

Table 1 compares the temporal evolution predicted by the analytical KGD solution and the FDEM model. Across the five reported injection times, the relative errors in fracture half-length range from 1.17% to 2.73%, while the errors in maximum fracture aperture range from 1.26% to 2.68%. All deviations remain below 3%, supporting the ability of the numerical framework to reproduce the benchmark fracture-growth response before it is applied to the multi-cluster cases.

3.2. Parametric Analysis of Multi-Cluster Fracture Propagation

A two-dimensional domain measuring 400 m × 200 m with three injection points is used for the parametric simulations. Unless a parameter is varied in a sensitivity case, the baseline properties are those listed in Table 2. The analysis focuses on four variables that directly affect the fracturing process: horizontal principal stress difference, cluster spacing, injection rate, and fracturing-fluid viscosity.

3.2.1. Effect of Horizontal Principal Stress Difference

The horizontal principal stress difference was set to 2, 6, 10, and 14 MPa. At low stress differences, fracture paths exhibit stronger deflection and branching, producing Y-shaped and dendritic geometries and a larger stimulated fracture-network area. The reported total fracture lengths for the low-stress-difference cases reach 215.4 and 187.4 m. As the stress difference increases, the stress field increasingly constrains propagation to a preferred direction, and Mode-I opening becomes dominant. At 14 MPa, the fractures are comparatively straight and parallel with little lateral branching, and the total fracture-network length decreases to 142.1 m. These results indicate that increasing stress anisotropy reduces the opportunity for fracture-path reorientation and therefore lowers network complexity.

Figure 1. Fracture geometries under different horizontal principal stress differences. (a) 2 MPa; (b) 6 MPa; (c) 10 MPa; (d) 14 MPa.

Figure 1. Fracture geometries under different horizontal principal stress differences. (a) 2 MPa; (b) 6 MPa; (c) 10 MPa; (d) 14 MPa.

3.2.2. Effect of Cluster Spacing

For three simultaneously propagating fractures, cluster spacing was set to 5, 10, 15, and 20 m. Fluid pressure and fracture opening generate induced stresses around each fracture, producing a stress-shadow interaction that is strongest when adjacent clusters are close. As spacing increases, the compressive interference acting on the middle fracture weakens. The middle-fracture length consequently increases from 18.5 m at 5 m spacing to 39.8 m at 20 m spacing, while its share of the total three-fracture length rises from 11.78% to 24.43% (Table 3). At small spacing, the central fracture is strongly suppressed and may arrest, bifurcate, or redirect toward a lower-resistance path. Thus, cluster spacing controls not only individual fracture length but also the uniformity of propagation among simultaneously stimulated clusters.

Figure 2. Fracture geometries under different cluster spacings. (a) 5 m; (b) 10 m; (c) 15 m; (d) 20 m.

Figure 2. Fracture geometries under different cluster spacings. (a) 5 m; (b) 10 m; (c) 15 m; (d) 20 m.

3.2.3. Effect of Injection Rate

Injection rates of 3, 6, 9, and 12 m3/min were evaluated. According to fluid mass conservation and the cubic-flow relation, a higher injection rate increases the rate of pressure accumulation within the fractures and therefore raises the mechanical driving force for opening and propagation. The simulations show that increasing the injection rate from 6 to 12 m3/min increases the maximum dynamic fracture aperture from 1.4 to 3.2 mm. Higher-rate cases also develop more secondary branches near the primary-fracture tips, indicating that elevated net pressure activates a larger number of damageable interfaces and increases fracture-network complexity.

Figure 3. Fracture geometries under different injection rates. (a) 3 m3/min; (b) 6 m3/min; (c) 9 m3/min; (d) 12 m3/min.

Figure 3. Fracture geometries under different injection rates. (a) 3 m3/min; (b) 6 m3/min; (c) 9 m3/min; (d) 12 m3/min.

3.2.4. Effect of Fracturing-Fluid Viscosity

Fracturing-fluid viscosities of 1, 10, 50, and 100 mPa·s were investigated. Fluid viscosity alters both pressure transmission along the fracture and leak-off into the porous matrix. In the 1 mPa·s case, stronger fluid penetration into the matrix increases pore pressure around the fracture and reduces effective confinement, which is associated with more extensive shear-related branching and a fine, reticulated fracture pattern. At 100 mPa·s, leak-off is reduced and pressure is concentrated more strongly within the primary fracture, favoring Mode-I opening. The resulting fractures are wider but less branched; the reported average width reaches 2.6 mm, while the branching frequency is 72% lower than in the 1 mPa·s case. The results therefore demonstrate a trade-off between fracture width and geometric complexity as viscosity changes.

Figure 4. Fracture geometries under different fracturing-fluid viscosities. (a) 1 mPa·s; (b) 10 mPa·s; (c) 50 mPa·s; (d) 100 mPa·s.

Figure 4. Fracture geometries under different fracturing-fluid viscosities. (a) 1 mPa·s; (b) 10 mPa·s; (c) 50 mPa·s; (d) 100 mPa·s.

3.3. Integrated Process Optimization

The sensitivity results were combined to identify a process window that favors fracture-network development. A high injection rate increases fracture pressure and promotes opening and secondary damage, whereas excessively close cluster spacing intensifies stress-shadow suppression of the central fracture. Low-viscosity slickwater favors leak-off-assisted pressure diffusion and branching. On this basis, the simulation matrix in Table 4 indicates that Schedule 4—10 m cluster spacing, 12 m3/min injection rate, and 1 mPa·s fluid viscosity—produces the most developed branching pattern among the four tested schedules (Figure 5). This combination should be interpreted as the optimum within the investigated numerical cases rather than as a universal field optimum.

4. Discussion

4.1. Coupled Mechanisms Governing Multi-Cluster Fracture Growth

The simulations show that the four investigated variables act through two coupled mechanisms: the redistribution of stress around neighboring fractures and the redistribution of fluid pressure between fractures and the porous matrix. Horizontal stress anisotropy sets the background directional constraint on fracture growth. Cluster spacing then controls the magnitude of fracture-induced stress shadows superimposed on that background field. Injection rate governs the rate at which net pressure and fracture opening are generated, while viscosity influences how pressure is partitioned between along-fracture transport and leak-off into the matrix. The resulting fracture geometry is therefore not controlled by any single parameter; it emerges from the competition between directional stress confinement, fracture–fracture interaction, pressure buildup, and fluid penetration.

4.2. Implications for Hydraulic-Fracturing Design

From a process-design perspective, the results support coordinated rather than independent selection of cluster spacing, injection rate, and fluid rheology. Very small cluster spacing can increase the number of initiation sites but may reduce propagation uniformity because the central fracture experiences stronger stress-shadow suppression. Increasing injection rate can partly offset propagation resistance by raising net pressure, but the present simulations also show that fluid viscosity changes whether this additional energy produces wider primary fractures or a more highly branched network. The combined schedule analysis therefore favors high-rate injection with moderate spacing followed by low-viscosity slickwater when the design objective is fracture-network complexity and SRV. Field application should nevertheless account for pumping limits, proppant transport, formation heterogeneity, and containment, which are outside the present model scope.

4.3. Model Limitations and Future Work

Several limitations define the range of applicability of the present results. First, the model is two-dimensional and therefore does not represent vertical fracture growth, height containment, or three-dimensional stress interaction. Second, the rock matrix is treated as a saturated linear-elastic porous medium and the injected fluid as incompressible and Newtonian; shale heterogeneity, anisotropy, non-Newtonian rheology, and multiphase effects are not explicitly represented. Third, potential fracture paths are controlled by the embedded cohesive-interface discretization, and the current study does not report a systematic mesh-sensitivity or cohesive-parameter uncertainty analysis. Finally, the integrated optimization is based on the tested numerical schedules and is not calibrated to a specific field treatment. Future work should therefore include three-dimensional simulations, heterogeneous and anisotropic shale properties, explicit natural-fracture networks, mesh and time-step sensitivity analyses, and validation against laboratory or field fracture diagnostics.

5. Conclusions

A two-dimensional hydro-mechanical FDEM framework was developed and verified to investigate multi-cluster hydraulic-fracture propagation in shale. The model couples matrix deformation and seepage, fracture flow, cohesive damage, and post-failure contact. Comparison with the KGD benchmark yielded relative errors below 3% for fracture half-length and maximum aperture. The main findings are as follows:

(1) Horizontal principal stress difference strongly controls fracture-path complexity. Low stress anisotropy permits greater deflection and branching, whereas high stress anisotropy promotes straighter, more parallel Mode-I-dominated propagation; the reported total fracture-network length decreases to 142.1 m at a 14 MPa stress difference.

(2) Cluster spacing governs stress-shadow competition. Increasing spacing from 5 to 20 m increases the middle-fracture length from 18.5 to 39.8 m and raises its share of the total three-fracture length from 11.78% to 24.43%, demonstrating improved propagation uniformity as inter-fracture interference weakens.

(3) Injection rate controls pressure buildup and fracture opening. Increasing the rate from 6 to 12 m3/min increases the maximum dynamic aperture from 1.4 to 3.2 mm and promotes additional secondary branching.

(4) Fracturing-fluid viscosity produces a trade-off between branching complexity and primary-fracture width. The 1 mPa·s case favors leak-off and branching, whereas the 100 mPa·s case produces wider, more planar fractures and a 72% lower branching frequency.

(5) Within the investigated parameter space, the combined schedule using a high injection rate, moderate cluster spacing, and low-viscosity slickwater generates the most developed fracture network. This result provides a process-screening guideline, but field-scale optimization requires further calibration and three-dimensional validation.