A Study on Fracture Propagation of Hydraulic Fracturing in Oil Shale Reservoir Under the Synergistic Effect of Bedding Weak Plane–Discrete Fracture

Abstract

Hydraulic fracturing is a critical process in the development of oil shale reservoirs. The presence of widespread bedding planes and natural fractures significantly influences the propagation of hydraulic fractures. Additionally, the injection point density plays a crucial role in the effectiveness of reservoir reconstruction. The Global Embedded Cohesive Zone Method (FEM-CZM) was employed to model the initiation and propagation of fractures from perforation holes, considering the combined effects of bedding planes and natural fractures. The results indicate the following: (1) the initiation and propagation of fractures from perforation holes lead to the co-propagation of two to four asymmetric main fractures, alongside open bedding planes and natural fractures; (2) larger bedding plane thickness and smaller bedding plane spacing promote hydraulic fractures’ tendency to propagate along the bedding planes, resulting in longer fracture lengths and predominance of tensile failure; and (3) a higher in situ stress difference facilitates the fracture’s penetration of the bedding plane, leading to an initial increase and subsequent decrease in fracture length. Tensile failure remains dominant, while the proportion of shear failure increases. Based on these findings, it is recommended to select fracturing sites with thicker bedding planes, larger bedding plane spacing, and a smaller vertical in situ stress field. Additionally, a perforation scheme with six injection points should be adopted to enhance the formation of high-efficiency seepage and heat transfer channels between hydraulic fractures, bedding planes, and natural fractures.

1. Introduction

Oil shale is regarded as a significant source of alternative oil and is widely distributed across the globe [1,2]. China’s oil shale resources are vast, with considerable potential for exploitation. The development of oil shale is primarily divided into two approaches: surface retorting and in situ thermal conversion [3]. Surface retorting of oil shale, however, results in environmental pollution, which conflicts with China’s carbon peak and carbon neutrality goals. As a result, in situ underground conversion technology, characterized by lower pollution and promising development outcomes, has garnered considerable attention. Nevertheless, the low permeability, porosity, and thermal conductivity of oil shale present significant challenges for in situ conversion mining. Hydraulic fracturing serves as a core technology for the efficient development of oil shale. The goal is to establish an “efficient heating pyrolysis fracture network” that facilitates heat transfer within oil shale reservoirs and provides pathways for pyrolysis products to migrate to production wells [4]. Currently, Jilin University in China has implemented in situ hydraulic fracturing and heat injection mining in the Fuyu and Nong’an pilot demonstration projects, yielding positive experimental results. Taking the Ordos Basin oil shale as an example, the oil shale displays well-developed bedding planes and an irregular in situ stress field. Through the study of outcrops, core samples, and image log data, some researchers have identified that weak planes (such as natural fractures and bedding planes) are prevalent in shale. However, few studies have examined the potential significance of these weak planes [5]. Therefore, investigating the influence of bedding planes and natural fractures on fracture propagation during hydraulic fracturing in oil shale reservoirs is of great importance.

Perforation design plays a crucial role in the success of hydraulic fracturing. Proper perforation design can significantly enhance reservoir stimulation, while poor perforation design can lead to hydraulic fracturing failure and hinder the migration of proppant. The perforation factor is an important parameter that influences fracture initiation near the wellbore and affects the subsequent propagation of hydraulic fractures. Currently, most studies focus on the influence of perforation parameters on fracture initiation and propagation through large-scale laboratory experiments [6,7,8]. Other studies investigate hydraulic fracture propagation using acoustic emission technology [9], introducing preset fractures into cement during fracturing experiments [10,11,12], or generating random fractures through thermal cycling [13]. However, there is limited research on perforation and fracture initiation and propagation using numerical simulation methods. Ispas et al. [14] employed numerical simulations to examine the effects of wellbore offset and azimuth on fracture propagation in sandstone formations, as well as the transition of fractures from a single planar fracture to a non-planar transverse fracture in deep-water oil fields in the Gulf of Mexico. Zhang et al. [15] simulated the initiation and propagation of multiple hydraulic fractures in the near-wellbore region, highlighting the critical role of the coupling effect between fluid flow and rock deformation in fracture initiation, extension, and relocation. The study also observed competitive propagation between fractures when multiple fractures were propagating simultaneously. Zhao et al. [16] investigated the influence of perforation friction and perforation cluster spacing on multi-fracture propagation and proposed an optimization scheme for uniform fracture propagation. While numerous experimental and numerical studies have explored the internal mechanisms of perforation and fracture initiation and propagation, there remains a lack of detailed descriptions regarding the dynamic process of perforation density and its effect on fracture initiation and propagation at individual perforation holes.

Currently, numerical simulations of hydraulic fracturing primarily involve four methods: the Displacement Discontinuity Method (DDM) [17,18], the Distinct Element Method (DEM) [19,20], the Extended Finite Element Method (XFEM) [21,22], and the Finite Element Method (FEM). Among these, the FEM has been widely used in hydraulic fracturing studies involving various types of rocks. The Cohesive Zone Method (CZM), which incorporates viscous, elastic, and plastic damage constitutive models, is increasingly applied in hydraulic fracturing simulations. For instance, the CZM has been successfully used to simulate the hydraulic fracturing process and the interaction between hydraulic fractures and natural fractures [23]. Additionally, based on the CZM, studies have explored the initiation and propagation of fractures under different geological and fracturing parameters [24], as well as the influence of geological factors and stress shadow on fracture propagation [25]. The CZM has also been employed to investigate the effects of viscosity parameters and fluid viscosity on fracture propagation [26]. In existing research, when studying the impact of bedding planes and natural fractures on hydraulic fracture propagation, a single cohesive element is often inserted into the reservoir rock to represent bedding planes or natural fractures in the simulation. However, this one-layer cohesive element fails to account for the thickness of the bedding plane and does not accurately reflect the expansion of hydraulic fractures along arbitrary paths.

Consequently, to address the challenges posed by the development of bedding planes and natural fractures in oil shale reservoirs, a stress-damage-flow model for hydraulic fracture propagation was established. Natural fractures were incorporated into the model using random functions to simulate the in situ conditions of the reservoir. Numerical simulations were performed based on the global FEM-CZM, enabling more accurate predictions of fracture behavior. This work analyzed the effects of bedding plane thickness, bedding plane spacing, and in situ stress field on the initiation and competitive propagation of hydraulic fractures under the combined influence of weak bedding planes and discrete natural fractures. Practical recommendations were provided for optimizing hydraulic fracturing in oil shale reservoirs, particularly in terms of injection point configuration. These insights offer valuable guidance for improving the hydraulic fracturing efficiency in unconventional reservoirs, such as oil shale.

2. Methods

This study uses Dassault Systèmes’ ABAQUS 2022a finite element software to establish the model. The core mathematical equations of the model include five key components: the linear elastic traction–separation criterion, the fracture initiation damage criterion, the fracture damage evolution criterion, the fluid flow model in hydraulic fractures, and the fluid flow–geomechanics coupling model [3,4,13,16].

2.1. Linear Elastic Traction–Separation Criterion

It is assumed that the stress and corresponding strain of the cohesive bonding element follows a linear elastic relationship at the onset of damage and prior to its evolution, as expressed in the following formula. The corresponding traction–separation curve is shown in Figure 1a.

$$\sigma =\left\{\begin{array}{c}{\sigma }_{\mathbf{n}}\\ {\sigma }_{\mathbf{s}}\\ {\sigma }_{\mathbf{t}}\end{array}\right\}=\mathit{K}·\epsilon =\left(\begin{array}{ccc}{k}_{\mathbf{nn}}& k_{\mathbf{ns}}& k_{\mathbf{nt}}\\ k_{\mathbf{ns}}& k_{\mathbf{ss}}& k_{\mathbf{st}}\\ k_{\mathbf{nt}}& k_{\mathbf{st}}& k_{\mathbf{tt}}\end{array}\right)·\left\{\begin{array}{c}{\epsilon }_{\mathbf{n}}\\ {\epsilon }_{\mathbf{s}}\\ {\epsilon }_{\mathbf{t}}\end{array}\right\}=\mathit{E}\epsilon$$

(1)

where $\sigma$ is the stress vector, MPa; ${\sigma }_{\mathrm{n}}$, ${\sigma }_{\mathrm{s}}$, and ${\sigma }_{\mathrm{t}}$ are normal stress, the first tangential stress, the second tangential stress, MPa; K is the elastic stiffness matrix, GPa; ε is the strain vector; ${\epsilon }_{\mathrm{n}}$, ${\epsilon }_{\mathrm{s}}$, and ${\epsilon }_{\mathrm{t}}$ are the normal strain, the all-directional strain, and the second tangential strain, respectively; E is the elastic modulus, GPa.

2.2. Initial Fracture Damage Criterion

Currently, four hydraulic fracturing damage criteria are commonly used in ABAQUS 2022a software: the maximum nominal stress damage criterion, the maximum nominal strain damage criterion, the secondary stress damage criterion, and the secondary strain damage criterion. This model adopts the maximum nominal stress damage criterion.

$${\left\{{\displaystyle \frac{〈{\sigma }_{\mathrm{n}}〉}{{\sigma }_{\mathrm{n}}^0}},{\displaystyle \frac{{\sigma }_{\mathrm{s}}}{{\sigma }_{\mathrm{s}}^0}},{\displaystyle \frac{{\sigma }_{\mathrm{t}}}{{\sigma }_{\mathrm{t}}^0}}\right\}}_{\max }=1$$

(2)

where ${\sigma }_{\mathrm{n}}^0$, ${\sigma }_{\mathrm{s}}^0$, and ${\sigma }_{\mathrm{t}}^0$ are the normal ultimate stress, the all-direction ultimate stress, and the second tangential ultimate stress, MPa; 〈〉 indicates that the cohesive element only bears tensile pressure.

2.3. Fracture Damage Evolution Criterion

When the cohesive bond element meets the initial damage criterion, the damage evolution criterion governs the stiffness degradation rate of the material. The damage variable D represents the overall damage of the material, accounting for the combined effects of all evolution modes. The initial damage value is zero, indicating that the material is undamaged. Once the damage evolution mode is defined, damage begins to develop. As the load continues to increase, the damage variable D gradually evolves towards 1, according to the specified damage mode. At this point, the material stiffness is fully degraded, and the damage is complete.

$${\sigma }_{\mathrm{n}}=\left\{\begin{array}{l}\left(1-D\right){\overline{\sigma }}_{\mathrm{n}},{\sigma }_{\mathrm{n}}\ge 0\\ {\overline{\sigma }}_{\mathrm{n}},{\sigma }_{\mathrm{n}}\prec 0\end{array}\right.$$

(3)

$$\begin{array}{l}{\sigma }_{\mathrm{s}}=\left(1-D\right){\overline{\sigma }}_{\mathrm{s}}\\ {\sigma }_{\mathrm{t}}=\left(1-D\right){\overline{\sigma }}_{\mathrm{t}}\end{array}$$

(4)

where ${\overline{\sigma }}_{\mathrm{n}}$, ${\overline{\sigma }}_{\mathrm{s}}$, and ${\overline{\sigma }}_{\mathrm{t}}$ are the pressure predicted in the normal direction according to the classification criterion of undamaged frontline elastic traction, the pressure predicted in the first direction according to the classification criterion of undamaged frontline elastic traction, and the pressure predicted in the second tangential direction according to the classification criterion of undamaged frontline elastic traction, respectively, MPa.

If the linear displacement damage evolution model is considered, the damage evolution in ABAQUS software can be simplified as follows:

$$D={\displaystyle \frac{d_{\mathrm{m}}^{\mathrm{f}}\left(d_{\mathrm{m}}^{\max }-d_{\mathrm{m}}^0\right)}{d_{\mathrm{m}}^{\max }\left(d_{\mathrm{m}}^{\mathrm{f}}-d_{\mathrm{m}}^0\right)}}$$

(5)

where $d_{\mathrm{m}}^{\max }$ is the maximum displacement of the element, m; $d_{\mathrm{m}}^{\mathrm{f}}$ is the displacement when the element is opened, m; and $d_{\mathrm{m}}^0$ is the displacement when the element begins to damage, m.

2.4. Fluid Flow Model in Hydraulic Fracture

The fluid flow within the cohesive unit fracture is primarily divided into radial flow along the fracture direction and normal flow, as shown in Figure 1b. The formula for calculating the radial fluid flow along the cohesive element is as follows:

$$q={\displaystyle \frac{d^3}{12\mu }}\Delta p$$

(6)

The formula for calculating the normal fluid flow along the cohesive element is as follows:

$$\begin{array}{l}{q}_{\mathrm{t}}=c_{\mathrm{t}}\left(p_{\mathrm{f}}-p_{\mathrm{t}}\right)\\ q_{\mathrm{b}}=c_{\mathrm{b}}\left(p_{\mathrm{f}}-p_{\mathrm{b}}\right)\end{array}$$

(7)

where q is the flow rate in the fracture, m³/s; q_t is the flow into the top surface, m³/s; q_b is the flow into the bottom surface, m³/s; d is the width of the cohesive element, m; μ is the viscosity of the fracturing fluid, Pa·s; Δp is the pressure gradient of the cohesive unit, MPa; c_t and c_b are the filtration coefficients of the upper and lower surfaces of the cohesive unit, m²/s; p_f is the fracturing fluid pressure in cohesive unit, MPa; p_t is the pore pressure on the upper surface of cohesive element, MPa; and p_b is the surface pore pressure under the cohesive element, MPa.

3. Establishment and Verification of Fracture Propagation Numerical Model

3.1. Numerical Modeling

To investigate the influence of multiple injection points on fracture initiation and propagation, this study establishes five hydraulic fracturing scenarios with single clusters and multiple injection points.

Table 1. Model parameters for the five research examples.

Instantiation	Perforation Hole Density/Number	Perforation Range/m	Perforation Hole Spacing/10 m	Perforation Depth/m	Model Scale (Length × Width/m)
T1	2	10	—	0.2	100 × 100
T2	4		5.0
T3	8		3.3
T4	12		2.0
T5	16		1.43

presents the main parameters for each example. Each scenario consists of a single cluster within a section, with a perforation range of 1 m. The hole length is 0.2 m, and the injection points are evenly and symmetrically distributed within the 1 m range. The number of perforations in the 5 examples is 2, 4, 8, 12, and 16, respectively, as illustrated in Figure 2a. To facilitate the analysis of perforation initiation and expansion, each perforation is assigned a number, as shown in Figure 2a. In the actual fracturing process, the flow rate of the fracturing fluid may vary across the perforation holes. However, it is assumed that the fluid pressure along the wellbore remains relatively constant, so the flow rate assigned to each perforation hole is considered equal. The flow distribution is represented by the following relations: q₁ + q₂ = q₃ + q₄ = …; q₁ + q₂ + q₃ + q₄ + … = Q, where q represents the flow into each perforation hole, and Q denotes the total flow in the cluster.

The cohesive model is primarily established using the FEM-CZM, with the global cohesive element forming the potential path for hydraulic fracture propagation. The global embedded cohesive element defines three distinct sections representing cohesive bedding planes, fractures, and the matrix. The model size is 100 m × 100 m, with a bedding plane dip angle of 0°. The bedding plane thickness is set to 0, 0.1, 0.2, and 0.3 m, and the bedding plane spacing is varied at 0.8, 1.0, 1.5, and 2.0 m. Natural fractures are randomly generated intermittent fractures that follow a Weibull distribution. The length of these natural fractures ranges from 0.5 to 2.5 m, and the initial fracture width is 0.001 m. As shown in Figure 3, the X-axis is aligned with the minimum horizontal principal stress, the Y-axis is aligned with the maximum horizontal principal stress, and the Z-axis is aligned with the vertical principal stress. The in situ stress field data are based on the in situ stress test results for continental Mesozoic shale gas reservoirs in the Ordos Basin by Wang et al. [27]. The magnitudes of the in situ stress field (${\sigma }_{\mathrm{h}}$, ${\sigma }_{\mathrm{H}}$, ${\sigma }_{\mathrm{V}}$) are (8, 9, 9 MPa), (8, 9, 10 MPa), (8, 9, 12 MPa), and (8, 9, 13 MPa), respectively. The model follows the principle of effective stress and super-hydrostatic pressure, with a corresponding pore pressure boundary set to 0. Fixed boundary conditions are applied in the X and Z directions, where the degrees of freedom are constrained to 0. In the XOZ plane, the degrees of freedom are unconstrained.

The mechanical parameters of the rock matrix, bedding plane, and natural fractures in the numerical model are based on the mechanical properties of oil shale in the Ordos Basin [28], as summarized in

Table 2. Related parameters of the hydraulic fracturing model.

Object	Parameter	Value
Oil shale matrix	Elastic modulus/GPa	6.66
	Poisson’s ratio	0.28
	Permeability/10−3 μm2	0.013
	Void ratio	0.02149
Oil shale matrix CZM unit	Elastic modulus/GPa	6.66
	$\mathrm{Normal} \mathrm{stress} {\sigma }_{\mathrm{n}}$/MPa	6
	$\mathrm{The} \mathrm{first} \tan \mathrm{gential} \mathrm{stress} {\sigma }_{\mathrm{s}}$/MPa	9
	$\mathrm{The} \sec \mathrm{ond} \tan \mathrm{gential} \mathrm{stress} {\sigma }_{\mathrm{t}}$/MPa	9
Oil shale bedding plane and natural fracture CZM unit	Elastic modulus/GPa	6.66
	$\mathrm{Normal} \mathrm{stress} {\sigma }_{\mathrm{n}}$/MPa	1
	$\mathrm{The} \mathrm{first} \tan \mathrm{gential} \mathrm{stress} {\sigma }_{\mathrm{s}}$/MPa	5
	$\mathrm{The} \sec \mathrm{ond} \tan \mathrm{gential} \mathrm{stress} {\sigma }_{\mathrm{t}}$/MPa	5
	Damage evolution mode displacement/mm	0.03
Fluid properties	Filtration coefficient/m2·s−1	$1 \times$ 10−14
	Fracturing fluid viscosity/Pa·s	0.001
	Fracturing fluid density/kg·m−3	9800
	Injection speed/m3·s−1	0.001
	Injection time/s	10

.

Figure 3. Comparison of fracture propagation paths derived from numerical simulations and experimental results: (a) numerical simulation results; (b) experimental result [29].

Figure 3. Comparison of fracture propagation paths derived from numerical simulations and experimental results: (a) numerical simulation results; (b) experimental result [29].

3.2. Numerical Model Validation

To verify the reliability of the numerical model, the fracture propagation path under the same conditions was simulated, based on the experiment conducted by Sun et al. [29] on the influence of shale bedding plane orientation and strength on hydraulic fracturing. The size of the numerical model matches that of the original physical simulation test, which is 400 mm × 400 mm. The maximum horizontal principal stress is applied in the X direction (4 MPa), the minimum horizontal principal stress is applied in the Y direction (2 MPa), and the bedding plane dip angle is 20°. The parameters used for model verification are the matrix tensile strength is 0.775 MPa, and the bedding plane tensile strength is 0.344 MPa.

Figure 3 compares the fracture propagation paths between the numerical simulation results and the experimental results. Figure 3a shows the numerical simulation results, while Figure 3b presents the experimental data from the literature. The fracture propagation paths in both the numerical simulation and experimental results are highly consistent, demonstrating similar propagation behaviors. This confirms the reliability of using the global FEM-CZM to simulate hydraulic fracture propagation paths.

4. Analysis of Modeling Results

4.1. Simulation Result Analysis

4.1.1. Effect of Multiple Injection Points on Fracture Initiation

Figure 4 illustrates the fracture opening distribution at different simulation times, magnified 100 times, for a scenario with 8 injection points. In this case, the bedding plane thickness is 0 m, the bedding plane spacing is 1.0 m, and the stress field (${\sigma }_{\mathrm{h}}$, ${\sigma }_{\mathrm{H}}$, and ${\sigma }_{\mathrm{V}}$) is (8, 9, 12 MPa). Based on the interaction between hydraulic fractures in the matrix, bedding planes, and natural fractures, the fracture propagation process can be divided into three stages.

Stage 1: From 0 to 0.344 s, the hydraulic fracture initiates and propagates within the matrix between two bedding planes. As fracturing fluid continues to be injected, the fracture expands vertically, with both its length and width increasing, as shown in Figure 4a.

Stage 2: From 0.344 to 0.853 s, the hydraulic fracture reaches the bedding plane, and the fracture penetrates the bedding plane, expanding within the matrix. During this stage, the length and width of the fracture increase more slowly, as seen in Figure 4b.

Stage 3: From 0.853 to 10.000 s, the hydraulic fracture progressively opens the bedding plane and propagates along it. As injection pressure increases, the fracture penetrates the bedding plane, continuing to propagate vertically, and eventually interacts with natural fractures. This results in a fracture propagation pattern where the hydraulic fracture directly penetrates and opens natural fractures.

Over time, more complex fracture networks emerge. Notably, at the start of the second stage, fractures at some injection points that were initially dormant begin to initiate and propagate due to the coupling effects of fluid injection and in situ stress. Consequently, some initial fractures are opened while others are suppressed. Ultimately, two main fractures penetrate the bedding plane and connect to natural fractures, with additional fractures propagating along the bedding plane.

Figure 5 illustrates the initiation time of each perforation at different injection points and whether the bedding plane is stimulated by hydraulic fracturing. In Figure 5, a value of “1” indicates that the perforation hole has initiated or the bedding plane has opened, while “0” indicates that it is either closed or has not initiated. From Figure 5, it is evident that there is a competitive interaction between perforation holes during the fracturing process. Due to stress interference and the balance between fluid and stress, the initiation of fractures at the bottom of the perforation holes occurs at different times. This behavior can be classified into the following three main competitive fracture initiation modes:

Mode 1: Perforation holes initiate early and continue to open and expand, as observed for No. 3 in T2, T3, and T4, and No. 5 in T5.

Mode 2: The perforation hole does not initiate at any point, such as No. 3, No. 7, and No. 8 in T5.

Mode 3: Fractures initiate initially and then close, as seen for No. 11 and No. 12 in T4.

The fracture initiation rate, defined as the ratio of initiated perforation holes to the total number of perforation holes, is examined at a fracturing time of 3 s. To ensure the accuracy and statistical significance of the data analysis, the numerical simulation results were examined for any anomalies or outliers [30,31]. The results are shown in Figure 6a. As depicted in Figure 6a, under the same total injection flow rate, the initiation rate of perforation holes decreases with an increasing number of injection points. This is attributed to a reduction in the flow rate per perforation hole and the increased stress interference among the perforation holes.

4.1.2. Effect of Multiple Injection Points on Fracture Propagation

Figure 7 compares the fracture opening distribution of T3, T4, and T5 at 10 s, under the conditions where the bedding plane thickness is 0 m, the bedding plane spacing is 1.0 m, and the stress field (${\sigma }_{\mathrm{h}}$, ${\sigma }_{\mathrm{H}}$, and ${\sigma }_{\mathrm{V}}$) is (8, 9, 12 MPa). The total length of the hydraulic fractures is shown in Figure 6b. By combining Figure 5 and Figure 7, it can be observed that, although complex fractures form in the early stages, these micro-fractures evolve into simpler main fractures in the middle and late stages due to the stress shadow effect. This results in a propagation mode characterized by the co-propagation of asymmetric main fractures, open bedding planes, and natural fractures.

Due to the limitation on the total amount of fracturing fluid and the competitive nature of fracture expansion, the number of fractures is highest in the early stages, gradually decreasing in the middle and late stages. Ultimately, a fracture pattern dominated by two to four main fractures forms. From Figure 5 and Figure 7, it is evident that as the number of injection points increases, the bedding plane becomes more easily opened, and the total length of the formed fractures increases. This indicates that increasing the number of injection points within a certain range can enhance both the total length and complexity of the fractures.

4.2. Analysis of Influencing Factors

To further investigate the influence of the synergistic effect between bedding weak planes and discrete fractures on the initiation and propagation of hydraulic fractures at multiple injection points, three experimental schemes were designed. These schemes aim to examine the effects of bedding plane thickness, bedding plane spacing, and in situ stress field on hydraulic fracture propagation.

4.2.1. The Influence of Bedding Plane Thickness

Based on the hydraulic fracturing model established in this study, the propagation behavior of hydraulic fractures was investigated under the conditions of eight injection points, an in situ stress field of (8, 9, 9 MPa), bedding plane spacing of 1 m, and varying bedding plane thicknesses of 0, 0.1, 0.2, and 0.3 m, while keeping other parameters constant. A bedding plane thickness of 0 m indicates the absence of a physical thickness, though the plane still has attributes like cementation strength. The simulation results are shown in Figure 8. The key findings are as follows:

Bedding plane thickness = 0 m: Hydraulic fractures intersect with and propagate along the bedding plane, reaching the maximum fracture width. Perforation holes No. 1, No. 4, and No. 6 close after connecting to the bedding plane, while No. 5 and No. 14 primarily propagate along the bedding plane with limited vertical extension.

Bedding plane thickness = 0.1 m: Hydraulic fractures continue to propagate along the bedding plane after intersection, significantly increasing fracture length. Perforation hole No. 4 closes after connecting to the bedding plane, while No. 12, No. 14, and No. 16 intersect and propagate along the bedding plane. Perforation holes No. 3 and No. 8 mainly propagate along the bedding plane with limited vertical extension.

Bedding plane thickness = 0.2 m: After intersecting the bedding plane, hydraulic fractures show slight vertical propagation. Only No. 3 and No. 4 propagate along the bedding plane, while No. 12, No. 14, and No. 16 intersect, connect to the bedding plane, and close.

Bedding plane thickness = 0.3 m: Hydraulic fractures propagate along the bedding plane after intersection, with limited vertical extension observed primarily in perforation No. 8. The fracture width significantly decreases. Perforation holes No. 12, No. 14, and No. 16 intersect, connect to the bedding plane, and then close.

Consequently, increasing bedding plane thickness results in a stronger tendency for hydraulic fractures to propagate along the bedding plane, while vertical propagation weakens. The stress shadow effect further encourages a dominant mode of simple, asymmetric main fractures propagating along the bedding plane.

As shown in Figure 9a, the fracture failure mode is primarily tensile failure. When the bedding plane thickness is 0 m, the proportion of tensile failure is highest. As the bedding plane thickness increases (from 0.1 m to 0.3 m), the proportion of tensile failure also increases, indicating that the in situ stress significantly influences the compaction effect of the bedding plane. However, this effect gradually diminishes as the bedding thickness increases. When the bedding plane is thin, the relatively weak cohesion between rock particles makes the bedding more prone to reaching its tensile strength, leading to tensile failure. Furthermore, the in situ stress field plays a crucial role in the propagation of fractures within the bedding plane. For thinner bedding planes, there is a higher stress concentration near the fracture tip, which may facilitate the initiation and propagation of fractures along the bedding plane. In contrast, for thicker bedding planes, the stress is more uniformly distributed, allowing fractures to propagate in more complex patterns. This explains why the proportion of tensile failure increases with bedding thickness: as the bedding plane becomes thicker, its stress distribution becomes less concentrated and its tensile resistance weakens, making tensile failure more likely to occur.

4.2.2. The Influence of Bedding Plane Spacing

The propagation behavior of hydraulic fractures under different bedding plane spacings (0.8, 1, 1.5, and 2 m) and an in situ stress field of (8, 9, and 9 MPa) was further studied, with the bedding plane thickness fixed at 0 m. The results are shown in Figure 10, as follows:

Bedding plane spacing = 0.8 m: The hydraulic fractures propagate entirely along the bedding plane with no vertical propagation. Perforation holes No. 12 and No. 14 initiate, form an intersection, and propagate along the bedding plane. Similarly, No. 5 and No. 13 initiate, connect to the bedding plane, and propagate along it, with No. 5 closing gradually after connection.

Bedding plane spacing = 1.0 m: Perforation holes No. 1, No. 4, No. 5, No. 6, and No. 14 initiate and propagate along the bedding plane after connecting to it. No. 1, No. 4, and No. 6 close after propagating along the bedding plane, while No. 5 and No. 14 mainly propagate along the bedding plane, with slight vertical propagation.

Bedding plane spacing = 1.5 m: Perforation holes No. 1 and No. 3 initiate, form an intersection, and propagate along the bedding plane. Holes No. 4, No. 15, and No. 16 initiate, connect to the bedding plane, and propagate along it, with slight vertical propagation in No. 4 and No. 15.

Bedding plane spacing = 2.0 m: Perforation holes No. 2, No. 11, and No. 16 close after connecting to the bedding plane. No. 11 and No. 16 exhibit significant vertical propagation and tend to open natural fractures upon intersection. No. 1 initiates, and shorter fractures form in the matrix at the bedding plane’s center, closing during the later stages.

These results show that fracture competition plays a critical role in propagation. When the bedding plane spacing is small, hydraulic fractures tend to propagate along the bedding plane, which restricts vertical propagation. However, as the bedding plane spacing increases, the fracture faces less resistance in reaching the bedding plane, allowing for more vertical propagation.

To further analyze the fracture failure patterns presented in Figure 9b, it can be observed that when the bedding spacing is 0.8 m, the proportion of tensile failure is relatively high. This is primarily because, at smaller bedding spacings, hydraulic fractures predominantly propagate along the bedding planes, leading to the occurrence of tensile failure. As the bedding spacing increases (from 1.0 m to 2.0 m), the proportion of tensile failure gradually decreases. This is because, over time, hydraulic fractures tend to traverse more of the rock matrix rather than propagate along the bedding planes. Furthermore, larger bedding spacings reduce the energy of the hydraulic fracturing fluid upon reaching the bedding planes, thereby weakening the driving force for fractures to extend along the bedding. As a result, the fracture propagation path becomes more complex, and the extension of hydraulic fractures is no longer confined to the bedding plane direction, leading to a reduction in tensile failure. With the increase in bedding spacing, the fracture length also increases. The increase in fracture length is influenced by the following two factors: first, the increase in bedding spacing itself; and second, as the bedding spacing increases, the difficulty of vertical propagation of hydraulic fractures decreases, making it easier for hydraulic fractures to connect with natural fractures. The connection with natural fractures facilitates smoother communication between the hydraulic fractures and bedding planes, which, in turn, aids in volume conversion and effective fluid permeation.

4.2.3. The Influence of In Situ Stress Field

To analyze the influence of the in situ stress field on hydraulic fracture propagation, dynamic simulations were performed under the conditions of eight injection points, 0 m bedding plane thickness, and 1 m bedding plane spacing. The in situ stress fields were set as (8, 9, 9 MPa), (8, 9, 10 MPa), (8, 9, 12 MPa), and (8, 9, 13 MPa), corresponding to vertical stress differences ($\Delta \sigma ={\sigma }_{\mathrm{V}}-{\sigma }_{\mathrm{h}}$) of 1.0, 2.0, 4.0, and 5.0 MPa. The results are shown in Figure 11. The following observations were made:

When the vertical in situ stress difference is 1.0 MPa, the hydraulic fracture primarily propagates along the bedding plane, with No. 5 and No. 14 fractures exhibiting slight vertical propagation.

When the vertical in situ stress difference is 2.0 MPa, fractures No. 5 and No. 14 propagate mainly along the bedding plane after connecting the bedding plane, but show significant vertical propagation. No. 1 and No. 10 propagate along the bedding plane after connection and gradually close.

When the vertical in situ stress difference is 4.0 MPa, fractures No. 5 and No. 14 primarily propagate along the bedding plane after the bedding plane connection, with significant vertical propagation. Additionally, No. 6 and No. 13 propagate along the bedding plane.

When the vertical in situ stress difference is 5.0 MPa, the hydraulic fracture directly penetrates the bedding plane and propagates in the matrix. Notably, fracture No. 10 propagates in the vertical direction, opening a natural fracture, and then gradually closing.

With the increase in vertical in situ stress difference, the hydraulic fracture exhibits greater opening in the vertical direction. This is because the in situ stress field influences the fracture propagation direction. A larger vertical in situ stress difference results in a compaction effect on the bedding plane, increasing the cohesion of the bedding plane. This increased cohesion raises the resistance to fracture propagation along the bedding plane, making it easier for fractures to penetrate the bedding plane.

Figure 9c shows that as the vertical in situ stress difference increases from 1.0 to 4.0 MPa, the proportion of tensile failure decreases. This is because, as the in situ stress difference increases, hydraulic fractures primarily propagate along the bedding planes. When the vertical in situ stress difference is large, the compaction effect on the bedding planes is enhanced, thereby increasing the resistance to fracture propagation along the bedding planes. As a result, although the fracturing fluid can penetrate the bedding plane, the higher resistance reduces the occurrence of tensile failure, promoting shear failure instead. This also forces the fractures to propagate more vertically, leading to an increase in fracture length. In cases with a significant vertical stress difference, the fracture propagation path becomes more complex, and the vertical extension of the fracture becomes more pronounced. When the vertical in situ stress difference reaches 5.0 MPa, the hydraulic fracture directly penetrates the bedding plane, resulting in tensile failure and the shortest fracture length. This suggests that a higher in situ stress difference overcomes the resistance of the bedding plane, causing the fracture to propagate rapidly in the vertical direction, leading to a reduction in fracture length. Once the hydraulic fracture intersects with natural fractures, it will open the natural fractures and continue propagating, causing the fracture mode to transition from tensile failure to shear failure.

4.2.4. Injection Point Density Optimization

As observed from Section 4.2.1, Section 4.2.2 and Section 4.2.3 hydraulic fractures tend to expand along the bedding plane with increasing bedding plane thickness. Conversely, smaller bedding plane spacing promotes the fractures to propagate toward the bedding plane, while larger bedding plane spacing reduces its obstruction. A smaller vertical in situ stress difference facilitates the opening of the bedding plane by hydraulic fractures. To enhance the communication between the bedding plane and natural fractures, improve reservoir stimulation, and establish efficient seepage and heat transfer channels, the injection point density was optimized under the conditions of bedding plane thickness of 0.3 m, bedding plane spacing of 2.0 m, and in situ stress fields values of (8, 9, and 10 MPa). The results are shown in Figure 12, where the number of primary fractures, fracture length, and fracture communication area (SRA) were used as indicators.

From Figure 12, it is evident that the fracturing ultimately resulted in an expansion mode dominated by two to four primary fractures. As the number of injection points increases, the number of main fractures initially increases, then decreases. The trends for fracture length and SRA are similar, both reaching their maximum values with six injection points. It is noteworthy that the fracture length and SRA at two injection points are both greater than those at four injection points. This suggests that increasing the number of injection points helps to increase the number of fractures. However, the increase in injection points does not significantly improve fracture propagation, and may even suppress effective fracture extension due to interference between injection points.

Analyzing the final number of main fractures, the number of main fractures is four when there are two and four injection points, three when there are six injection points, and two when there are eight injection points. When the number of main fractures is lower (six and eight injection points), fracture propagation seems to be more effective. This is because the main fractures formed after initiation of the fractures are concentrated and controlled in their dominant propagation direction. Although the number of fractures is smaller, the fracture length and SRA are larger.

From the analysis of fracture length and SRA, it is evident that fracture length shows a significant increasing trend as the number of injection points increases. Particularly, with six injection points, the fracture length reaches 50.68 m, the largest value among all injection point configurations. The variation in fracture area is consistent with the change in fracture length, with the maximum SRA (0.2016 m²) observed at six injection points. This indicates that an increase in injection point number promotes lateral fracture extension and the formation of effective seepage channels.

Considering both fracture length and SRA maximization, the number of injection points should be six, as this results in the highest fracture length and SRA, indicating the most effective fracture propagation. Additionally, although the number of main fractures is maximized (four) at two injection points, the fracture length and area are smaller, implying that fracture extension is not as effective as with six injection points. Therefore, three to four main fractures are an ideal choice, as they effectively avoid the competitive effects between fractures while expanding the extension area of individual fractures.

By comprehensively considering fracture length, SRA, and the number of main fractures, six injection points are the most optimal choice. This configuration not only maximizes fracture extension and forms more effective seepage channels, but also maintains an appropriate number of main fractures.

5. Discussion

5.1. Research Implications

In this study, a hydraulic fracturing model for oil shale was developed using the FEM-CZM, which simultaneously considers the synergistic effects of bedding planes, natural fractures, and in situ stress. This model’s accuracy was validated by comparing the results with experimental data from existing studies. According to the findings, the thickness of the bedding plane has a complex influence on the fracture propagation in multi-injection-point hydraulic fracturing. Specifically, as the bedding thickness increases, it becomes more difficult for the hydraulic fracture to penetrate the bedding plane, which is consistent with the results of Li et al. [4,32]. Therefore, fracturing could be more effectively conducted in reservoirs with larger bedding spacing and significant in situ stress differences. Injection point density, as one of the key parameters in field fracturing operations, plays a critical role. Selecting an appropriate injection point density can increase the number of main fractures and mitigate the impact of the stress shadow effect, thereby more efficiently enhancing oil shale reservoirs. Similar conclusions were also drawn by Ran et al. [33,34].

5.2. Model Limitation

Although the study reveals the impacts of bedding planes, natural fractures, multi-injection-point density, and in situ stress differences on the initiation and propagation of hydraulic fractures in oil shale reservoirs, it is important to note that the initiation and propagation of hydraulic fractures in oil shale are complex, multi-scale geo-mechanical problems. Several other critical factors, such as perforation orientation, injection volume, injection pressure, and the complexity of the bedding planes, also affect fracturing effectiveness. Therefore, there are still some unresolved issues in this study.

In the numerical model, to clearly demonstrate the initiation and propagation of fractures, multiple injection points were set only within a single cluster. However, when multi-cluster fracturing is used, hydraulic fracturing becomes more susceptible to stress shadow effects [35,36]. Thus, the findings and numerical model presented in this study lay the foundation for future research in this field, emphasizing the need for further investigation.

6. Conclusions

This study investigates the initiation and propagation of hydraulic fractures in oil shale reservoirs by thoroughly considering the synergistic effects of bedding planes, natural fractures, in situ stress differences, and injection point density. The following conclusions are drawn:

(1)

During the fracturing process, the initiation of fractures at the perforation holes follows three competitive modes: initial initiation and continued propagation; perforation hole opens but does not initiate; initial initiation, propagation, and subsequent closure.

Initially, the number of fractures formed is highest, with fracture numbers gradually decreasing in the middle and late stages. Ultimately, a propagation mode dominated by two to four asymmetric main fractures, along with open bedding planes and natural fractures, is established. As the number of injection points increases, the initiation rate of the perforation holes decreases due to several factors: reduced flow rate into the perforation holes, increased friction at the perforation points, and the influence of stress interference.

(2)

As the thickness of the bedding plane increases, the hydraulic fracture shows a stronger tendency to propagate along the bedding plane, with reduced vertical propagation. The fracture failure is predominantly tensile, and the stress shadow effect results in the formation of a simple main fracture along the bedding plane. Regarding fracture competition, when the bedding plane spacing is smaller, the hydraulic fracture is more likely to propagate along the bedding plane. Conversely, larger bedding plane spacing reduces the hindrance to vertical propagation, allowing the hydraulic fracture to connect with natural fractures, leading to longer fracture lengths. In these cases, fracture failure is mainly tensile. As the vertical in situ stress difference increases, fractures are more likely to penetrate the bedding plane. The fracture failure remains mainly tensile, but the proportion of shear fractures increases. The fracture length initially increases with the stress difference before decreasing.

(3)

Based on the above analysis, areas with large bedding plane thickness, wide bedding plane spacing, and a relatively small in situ stress field are more favorable for hydraulic fracture propagation. These conditions facilitate better communication between hydraulic fractures, bedding planes, and natural fractures. Additionally, selecting six injection points optimizes the reservoir transformation, enhancing the development of a complex fracture network and improving oil shale recovery.

(4)

This paper explores the competitive initiation and propagation behaviors of hydraulic fractures with multiple injection points under the combined influence of bedding planes and natural fractures. Future studies could focus on the effects of factors such as changes in bedding plane dip angle and fracturing fluid displacement on hydraulic fracture propagation. Additionally, evaluating the impact of complex fracture networks on in situ oil shale transformation would provide valuable insights for optimizing actual production processes.