Shale reservoirs are typically characterized by low porosity and ultra-low permeability. The key to effectively developing shale oil and gas is to form a large-scale complex fracture network by multi-stage and multi-cluster fracturing in horizontal wells. Different from conventional reservoirs, shale, as a type of layered sedimentary rock, has strong heterogeneity. Meanwhile, shale reservoirs have a large number of local discontinuities, such as faults, bedding planes and natural fractures, which result in significant structural anisotropy that exerts a certain impact on the hydraulic fracture propagation.

The research on the expansion mechanism of complex fracture networks in shale reservoirs has become an important topic in the field of unconventional oil and gas development. In recent years, a lot of theoretical research on the propagation of hydraulic fractures in shale reservoirs has been carried out, and a series of fracturing models have been established that can simulate propagation of multiple fractures, including the line network model^[1], discrete fracture network model^[2,3], unconventional fracture model^[4], and other models based on the finite element method^[5,6], boundary element method^[7,8,9,10], extended finite element method^[11,12], discrete element method^[13-16] and phase-field method^[17,18]. The above models simplify the rock as an isotropic medium and primarily focus on the fracture distribution under the influence of high-angle natural fractures. However, these models generally fail to consider the influences of typical characteristics of shale such as bedding and mechanical anisotropy on the morphology of fracture propagation. As demonstrated in the laboratory fracturing simulation and field fracture monitoring, bedding fractures have an important impact on the vertical extension and propagation morphology of hydraulic fractures^{[19,20,21,22,23]}. When encountering bedding surfaces, hydraulic fractures may result in penetration, diversion, termination, or stepped extension, which cause uncertainty of the shape of the final fracture network. As a result, most numerical models can not reflect the essential difference between hydraulic fracture propagation in laminated shale reservoirs and that in naturally fractured reservoirs.

Field coring data shows that shale reservoirs have highly dense bedding fractures (several to hundreds per meter) commonly. The strong anisotropy brings great challenges to the simulation of fracture propagation in shale reservoirs. It is accurate to analyze the influence of the weak surface of the natural fracture on fracture propagation using the cohesive element of the finite element^[24]. Nevertheless, when the natural fracture density is set too high, the calculation stability and convergence are significantly poor. For the shale with abundant weak bedding planes, the rock mass medium is more like a discrete body of discontinuous media. As a numerical simulation for discontinuous medium, the discrete element method has great advantages in dealing with large deformation problems of discontinuous medium such as rock and soil fractures (like highly complex fracture network)^[25]. Zhao et al.^[13] used the two-dimensional particle discrete element method to simulate the behavior of hydraulic and natural fractures, in which they simulated natural fractures by weakening the bond strength between particles. Zangeneh et al.^[14] used the two-dimensional discrete element method to simulate the hydraulic fracture network. In their model, the stratum was divided by multiple sets of joints, bounded by deformable rock blocks. Nagel et al.^[15] studied the types and influencing factors of hydraulic fractures in naturally fractured reservoirs using the 3D discrete element method. However, they didnot consider the influence of bedding fracture. Zou et al.^[26,27] established a three-dimensional discrete element fracture network model to explore the influence of bedding fractures on the propagation of hydraulic fractures in shale reservoirs at engineering scale, but they did not take into account the influence of the longitudinal reservoir stress heterogeneity on the hydraulic fracture propagation pattern.

Taking the Longmaxi shale reservoir in the main area of the Jiaoshiba anticline in the Fuling shale gas field as the research object, we carried out a series of experimental tests to investigate the mechanical anisotropy of the shale reservoir in this study. Then, a three-dimensional complex fracture propagation model was built and solved based on the discrete element method. Moreover, the propagation model was combined with the indoor experimental test to investigate the law of hydraulic fracture propagation in shale reservoirs under the influences of longitudinal stress difference and bedding.

In this study, the main area of the Jiaoshiba anticline is taken as the research object. The organic-rich, gas-bearing shale intervals of the Upper Ordovician Wufeng Formation and the Lower Silurian Longmaxi Formation are mainly concentrated in the bottom of Wufeng Formation-Longmaxi Formation, below which is a set of light-gray nodular limestone of the Upper Ordovician Jiancaogou Formation. According to the reservoir physical properties, shale layers in the Longmaxi and Wufeng Formations are divided into 9 sublayers from bottom to top. Gray-black carbonaceous and siliceous shale take dominance in Wufeng Formation-Longmaxi Formation. The main gas-bearing layers have rich beddings, and the numbers of beddings decrease from the bottom to the top^[28]. Among them, the sublayer ① is developed strongly in bedding fractures, with hundreds of bedding fractures per meter. In contrast, the sublayers ⑧ and ⑨ at the top have fewer bedding fractures. The core observation results are shown in Table 1. It can be seen that the high-angle natural fractures in the reservoir are generally small in scale and not connected with each other.

The vertical stress profile of Well Jiaoye A was interpreted based on the corrected logging curve obtained during the Kaiser in-situ stress test (Fig. 1). The minimum horizontal principal stress varies in the longitudinal direction. The in-situ stresses of sublayers ①-④ fluctuate within a small range (only 1-2 MPa), with an average of 49 MPa; the in-situ stress of sublayer ⑤ is significantly higher (up to 53 MPa), and there is a stress difference of 3-4 MPa between sublayer ⑤ and sublayers ①-④; and the in-situ stress of sublayer ⑥ increases further and reaching 56 MPa in maximum. In addition, the Jiancaogou Formation that underlies the target layer is a good stress-shielding layer, resulting in a stress difference of more than 12 MPa.

Taking the Fuling shale gas field as the research object, triaxial rock mechanics tests of shale samples taken from different sublayers and different coring directions of the same well were carried out with a dynamic rock triaxial test system under high temperature and high pressure conditions, and the results are shown in Fig. 2. It can be seen that elastic modulus E_h measured in the direction parallel to the bedding (0°) is generally higher than E_v measured in the direction perpendicular to the bedding (90°). The E_h/E_v value decreases from 1.26 in sublayer ① with many bedding fractures to 1.06 in sublayer ⑧ with few bedding fractures (Fig. 2a). When the loading stress acts vertically on the bedding surface, the compaction effect of the external force will lead to the closure of the bedding microfractures. A larger axial strain generally corresponds to smaller elastic and deformation moduli. The more developed the bedding fractures,the more obvious the phenomenon will be. As the confining pressure increases, and the compaction effect of the interlayer fractures enhances, the elastic moduli measured in different directions tend to be more similar (Fig. 2b). In addition, due to the influence of bedding fractures, the fracturing patterns of core samples from different bedding directions are different from each other. Generally, longitudinal tensile fractures occur more frequently in core samples with parallel bedding, while conjugate shear failures with tensile and shear fractures are predominant in core samples with vertical bedding.

To study the mechanical properties of bedding and their influences on the shear strength anisotropy of shale, the shear strengths of standard samples parallel to and perpendicular to the bedding direction were tested using a direct shear testing machine. The internal friction angle and cohesion of the samples were approximately obtained by linear regression (Fig. 3). As shown in the Mohr-Coulomb strength envelope curves obtained by the direct shear test, the shear stress required for shear failure of the rock mass in the case that the direction of the shear stress is parallel to that of the bedding direction is lower than that in the case that the shear stress direction is perpendicular to the bedding direction, under the same normal stress. At the same time, more developed bedding corresponds to lower shear strength in the direction parallel to the bedding. When sheared in the direction parallel to the bedding, the cohesion of sublayer ① with extremely developed bedding is 5.24 MPa, while that of sublayer ⑧ with few bedding is 15.18 MPa. Besides, affected by the degree of bedding development, the fracture surface is flat and smooth when the sample is sheared in the direction parallel to the bedding. In contrast, when the shearing direction is perpendicular to the bedding direction, composite failure modes such as matrix tensile splitting and bedding plane shear slip occur, the fracture surface is unsmooth and some fragments are generated by local splitting and peeling.

The mechanical parameters of shale, such as elastic parameters and shear strength, are anisotropic due to the existence of weak planes or bedding. For the shale gas reservoir with developed bedding, accurately describing the development degree of bedding and the anisotropic characteristics under its influence are the key to reflecting the difference of hydraulic fracture propagation in shale reservoirs and naturally fractured reservoirs.

The governing equation of the model is mainly composed of the fracturing fluid flow equation, the rock mass deformation equation, and the fracture failure criterion. This governing equation is solved by combining the finite element and discrete element methods^[26,27]. According to the discrete element method^[29], the solution domain of the stratum model (Ω_f) is discretized into several triangular-prism block elements that are linked by virtual springs to transfer the interaction force. The break of spring represents the rock fracture. The joint elements between all the contact block elements constitute the connected fracture network for fracturing fluid flow. The distribution of the fluid pressure in the contact block is calculated using the finite element method. The pressure is taken as external load acting on the fracture surface (i.e. the contact surface of the blocks), and the deformation of the block and the stress state of the spring are calculated. The spring break (fracture propagation) is determined by the maximum tensile stress criterion and the Mohr-Coulomb criterion. At the same time, the constitutive equation of transversely isotropic linear elastic material is used to replace the isotropic constitutive equation to study the influence of the mechanical anisotropy of shale rock on the morphology of fracture propagation^[26,27].

The fluid flow in the fracture is regarded as the plane plate flow of incompressible Newtonian fluid, which satisfies the continuity equation without considering the effect of gravity^[28]:

Due to extremely low permeability of the shale matrix, the matrix block is regarded impermeable and the filtration of the fracturing fluid is ignored, i.e., q_l=0. Then, the global mass balance equation in the fracture network is as follows:

When N cracks initiate simultaneously, and the total fluid rate entering all the fractures is Q, the fluid flow into each fracture depends on the width and pressure of the fracture in the process of its extension. The model assumes the fractures are completely filled with fluid, and there is no flow at the end of the hydraulic fractures.

The linear elastic dynamic equilibrium equation is as follows^[30]:

The boundary of the model block is fixed, and u_i=0; in the condition that a contact force is applied on the contact surface of the block, the fluid pressure p_i ($p_{i}=\sigma_{ij}\times n_{i}$) is applied on the fracture wall when the hydraulic fracture occurs. Since the stress-strain relationship conforms to the linear elastic constitutive equation^[30], then:

As a type of laminated sedimentary rock, shale can be regarded as a type of transversely isotropic material, namely, its elastic characteristics are the same in the bedding planes, but different in the direction perpendicular to the bedding plane^[31]. Five elastic constants, E_h, E_v, υ_h, υ_v, and G_v, are used to characterize the linear elastic characteristics of transversely isotropic rock. If the shale layer is horizontal, its flexibility coefficient matrix is as follows:

where the shear modulus G_v is^[32]:

There are normal and tangential springs between adjacent blocks, and the springs can be broken by tensile or shear slip, i.e. tensile failure or shear failure. Whether a spring is broken is judged according to the maximum tensile stress criterion and the Mohr-Coulomb criterion. When -F_n⁽ⁿ⁾<F_n,max (the tensile stress is negative) or $\bracevert F_{s}^{(n)}\bracevert$<F_s,max , t.here is no failure between adjacent contact points. The normal stress and shear stress in step n are solved with equations (7) and (8).

When F_n⁽ⁿ⁾≥F_n,max=AT₀, tensile failure occurs, the spring breaks, adjacent blocks are separated (Δu_n>0), and the interaction force disappears, namely:

When ${{F}_{\text{s}}}^{\left( n \right)}\]≥\[{{F}_{\text{s,max}}}=A{{S}_{0}}+{{F}_{n}}^{\left( n \right)}\text{tan}\varphi$, shear failure occurs, accompanied by shear displacement between adjacent blocks (Δu_n≤0, Δu_s>0). At this moment, different from the tensile failure, the shear failure is associated with compressive stress and friction resistance between blocks:

Under the action of the fluid pressure in the fracture, after the matrix block deforms, the width of the fracture changes. Meanwhile, the fluid rate in the fracture also changes with the width, which affects the fluid pressure in the fracture. In other words, the fluid pressure in the fracture and the fracture width affect each other. The weak coupling method is used to realize the iterative process of the fracturing fluid flow and solid deformation in fracture, while the flow equation and rock mass deformation equation are discretized respectively, and solved sequentially and iteratively^[26,27]. Before the tensile or shear failure occurs in the matrix block, the permeability of the initial joint unit is equal to that of the matrix block. The permeability of the matrix block is equivalent to the initial width of the fracture, that is, w₀=(12K₀)^1/2, and the original fracture width is designed for the initial flow of the fluid^[33]. Then, in the current time step, it is necessary to select appropriate experimental solutions p_m and w_m to solve the pressure p_m+1 in the next time step, then w_m+1 is calculated with equation (11). Generally, 0<β≤0.5 is taken, and the time step is small enough to facilitate the iterative convergence of the equation^[34].

The propagation patterns of symmetrical vertical fracture and horizontal radial fracture were simulated by the model, and compared with the analytical solutions of the classical fracturing models. The main input parameters were E=35 GPa, ν=0.2, μ=5 MPa·s, and Q=5 m³/min. The results are shown in Fig. 4. In the classical PKN model and the radial model, the energy consumption during fracture propagation is assumed to be mainly used for fluid flow in the fracture, and the influences of mechanical fracture parameters, such as fracture toughness and tensile strength of the rock, are not considered. The tensile strength is set at zero in the numerical model. The results of the numerical model are basically consistent with those of the PKN and radial models, with a difference of 3.6% and 4.9% respectively, which verifies the reliability of the numerical model.

Based on the reservoir parameters of the Jiaoshiba anticline in the Fuling shale gas field, a layered formation model was established. Considering that the formation is of horizontal occurrence, the boundary of the triangular-prism grid element is along the preset trace line of the interlayer interface (weak bedding plane). In the early stage of simulation, it was found that the fractures were not affected by the upper sublayers ⑦-⑨ in the simulation of sublayers ①-⑨. Out of the consideration of computational efficiency of the model, the influences of sublayers ⑦-⑨ are not considered in the following modeling. The model was 200 m wide, 500 m long and 50 m thick in X, Y, and Z directions, respectively. According to the results of laboratory experiments and in-situ stress interpretation from well logging, the 9 sublayers were divided into 3 intervals with different stresses vertically (Table 2). Among them, the total thickness of sublayers ①-④ at the bottom was 28 m; the thickness of sublayer ⑤ in the middle part was 10 m, and the stress difference with sublayer ⑥ was 3 MPa; the thickness of sublayer ⑥ in the upper part was 12 m. The horizontal well passed through the middle of sublayers ①-④, and the grid model is shown in Fig. 5. To simulate 3-cluster fracturing in a single stage, it was assumed that only one hydraulic fracture was produced in each perforation cluster. The basic fracturing simulation parameters were 3 clusters at the spacing of 25 m in one stage, displacement of 14 m³/min, and fracturing fluid viscosity of 2.5 MPa·s.

In the study area, bedding fractures are developed in general, especially in the sublayers ①-④ at the bottom of the Wufeng-Longmaxi Formations, bedding fractures reach the density of over one hundred per meter. During the simulation of the fracture propagation, it is necessary to use an equivalent method to approximate the bedding fractures to "explicit" weak bedding planes with the interval of 2-10 m before fracturing. At the same time, in order to study the influence of natural fractures on the morphology of the fracture growth, discrete high angle natural fractures 20 m long and 10 m high were set up in random distribution in the model. Given the fact that high-angle natural fractures in the reservoir are of small scales and not connected with each other, their linear density was set at 0.05/m. In reference to results of the laboratory rock mechanics test, the mechanical parameters of the shale matrix, natural fracture and weak bedding plane such as tensile strengths, cohesion and internal friction angles were determined (Table 3).

To illustrate the influence of the weak bedding plane on the distribution of artificial fractures in the heterogeneous interval, the simulated morphology of fracture propagation in the homogeneous reservoir interval of a single stage of horizontal well (without considering bedding and natural fractures) is taken as a comparison (Fig. 6a). The simulation with the above fracturing parameters showed three single-wing main fractures 450-480 m long were formed after fracturing in the homogeneous interval. Under the influence of the stress shadow, the fractures in the outer clusters deflected to the outside; the fractures in the middle clusters extended along the fracture height direction more than the outer cluster fractures under the compression of the outer fractures. The fracture near the wellbore passed through the interface of sublayers ⑤ and ⑥, and was up to 50 m high. With the increase of the distance to the wellbore, the height of the artificial fracture decreased gradually. Restricted by the stress difference of the interlayer, it can be seen that the fracture height at the fracture tip was limited to the top of sublayer ④, and couldn't enter sublayer ⑤, with the height of 28 m. In the case considering the influence of natural fractures (without considering bedding), with the opening of natural fractures, the complexity of the hydraulic fractures increased, and the main fractures decreased in length to about 330 m. Due to the existence of local high-angle cross-layer natural fractures, the fracture height was very large at local positions (Fig. 6b). In general, the near-wellbore fractures passed through two stress-shielding layers in the simulation not considering the influence of bedding.

4.2.1. Influence of the density of weak bedding plane

According to the physical modeling of fracturing in laboratory, a "fence-shaped" complex fracture network interwoven by bedding and hydraulic fractures is generated in the shale after fracturing^[16], and the density of the weak bedding plane would significantly affect the morphology of the propagation of artificial fractures. Under the basic fracturing operation parameters (N_C=3, Q=14 m³/min, μ=2.5 mPa·s) and medium strength (T_BP=4 MPa, S_BP=15 MPa, and φ_BP=25°) and 2 m of spacing between weak bedding planes, the fracture propagation morphologies at different injection times are shown in Fig. 7. In the early stage of expansion (5 min into injection in the simulation), the hydraulic fracture extended in the sublayers ①-④, accompanied with bedding opening (Fig. 7a); after 15 min of injection, due to the influence of the interlayer stress difference, the fracture height propagation was obstructed and stabilized at 28 m. The hydraulic fracture extend along the direction of the fracture length and the bedding opened continuously (Fig. 7b); after injection for 45 min, the hydraulic fractures in the middle cluster broke through sublayer ④ and expanded to sublayer ⑤. Finally, they were cut off by bedding fractures, and the fracture height was stable at 38 m (Fig. 7c).

Fig. 8 shows the numerical simulation results of the artificial fracture propagation at different weak plane densities. As a large-area, continuous weak plane, beddings in shale can increase the fracture density during the hydraulic fracturing, and thus improve the adequacy and effect of the reservoir stimulation. The larger the bedding density, the higher the hydraulic fracture density in unit volume will be. However, the opening of too many beddings will severely restrict the expansion of hydraulic fractures in the direction of length and height, consequently, the stimulated reservoir volume will greatly reduce. Under different bedding densities, the relationships between the length and height of artificial fractures and the stimulation time are shown in Fig. 9. When d_BP was 8 m (weakly developed), the fracture length and height expanded rapidly, and the fracture near the well went through the interface of sublayers ⑤ and ⑥ at 15 min. Then, the overall artificial fracture height was stabilized at 40 m, and the final simulated fracture length was 400 m. When d_BP was 2 m, the fracture length was less than 300 m, the hydraulic fractures ended at the interface of sublayers ⑤ and ⑥, and the fracture near the wellbore reached a maximum height of 37 m. In general, the fractures near the well mostly ended at the interface of sublayers ⑤ and ⑥, with a fracture height of about 38 m, which is obviously smaller than that in the model without considering the bedding effect. In contrast, the artificial fractures far from the wellbore terminated at the interface of sublayers ④ and ⑤, with a fracture height of about 28 m.

4.2.2. Influence of the strength of the weak bedding plane

The tensile strength and shear strength (including cohesion and internal friction angle) of beddings are used to describe whether the bedding is easy to open or not. To study the influence of the bedding strength on the fracture propagation morphology, four groups of experiment schemes were set up, in which different values of T_BP, S_BP, and φ_BP were input to simulate weak bedding planes with different strengths. The strength grading is shown in Table 4.

The bedding spacing was 2 m. The dynamic parameters of the fracture propagation and the fracture morphology at the end of the simulation are shown in Figs. 10 and 11, respectively. It can be seen that the artificial fractures in the shale reservoir with low-strength weak bedding planes were cut off by the bedding fractures near the well and only propagated horizontally along the bedding, and were less than 5 m high (Fig. 11a). In contrast, when the mechanical strength of the bedding is moderate to high, the fractures height expanded smoothly. When the hydraulic fractures passed through sublayer ④, they stopped expanding in height for a period of time due to the stress difference between sublayer ④ and ⑤. At the later stage of fracturing simulation (40 min), the middle cluster fractures penetrated through the stress shielding layer under the condition of medium-high to high bedding strength (Fig. 10). As the bedding strength increases, it becomes more difficult to open the bedding, while the fracturing fluid is more efficient in creating main fractures, thereby the fractures increase in length (Fig. 11d). When the bedding strength is weak, the bedding is easy to be opened by hydraulic fracture, then the hydraulic fracture would divert along the bedding, ending with much smaller fracture height.

4.2.3. Influences of fracturing engineering parameters

To investigate the correlation between fracturing parameters and fracture propagation morphology, under constant geological parameters, medium bedding strength (T_BP=4 MPa, S_BP=15 MPa, and φ_BP=25°), and bedding spacing of 2 m, simulations at different cluster numbers, displacements and fracturing fluid viscosities were conducted, and the results are shown in Fig. 12.

The same injection parameters, the less the clusters in one stage, the better it is to increase the fluid pressure in the fracture, thereby promote the expansion of the fracture height. In the simulation under 2 clusters in one stage and the cluster spacing of 25 m, two transverse fractures 400 m long were formed; the near-well artificial fractures extended to the top of sublayer ⑥, with a height of 48 m, while those far from the wellbore stopped growth at the interface of sublayers ⑤ and ⑥, with a height of 38 m (Fig. 12a). As the cluster number increases, the fracture height decreases. In the simulation under 3 (Fig. 7c) and 5 clusters (Fig. 12b) per stage, the artificial fractures near wellbore stopped growth at the interface of sublayers ⑤ and ⑥ longitudinally, with fracture heights of 38-40 m, while those far from the wellbore terminated at the interface of sublayers ④ and ⑤, with a height of 28 m and length of about 300 m. Compared with the 3 clusters fracturing, the 5 clusters fracturing was weaker in maintaining the fracture height along the fracture length direction, the fractures far away from the wellbore dropped in height quickly to 28 m, and all fractures reduced in height.

Injection displacement and fracturing fluid viscosity are also important engineering parameters to be considered in shale reservoir volume fracturing. As the displacement was increased from 12 m³/min to 16 m³/min, the propagation trend of artificial fractures was obvious, and the stimulation scope of sublayer ⑤ increased (Fig. 12c, 12d). Increasing the viscosity of the fracturing fluid can significantly reduce the restriction of the bedding and longitudinal interlayer stress difference on the fracture height propagation. When the viscosity was 1.0 mPa·s, the hydraulic fracture fully opened the bedding near the injection point, and the artificial fracture vertically ended at the interface of sublayers ④ and ⑤, with a fracture height of 28 m; meanwhile, due to the significant diversion effect of the bedding fracture, the length of the main fracture was only 168 m, as a result, the fracture-controlled reserves of the single well will drop significantly, which is unfavorable for the long-term stable production (Fig. 12e). As the viscosity of the fracturing fluid increased to 2.5 mPa·s, the fracture height reached sublayer ⑤, and the beddings at different positions were opened to various degrees. As the viscosity of fracturing fluid was further increased to 25 mPa·s, the hydraulic fractures propagated along the direction of the maximum principal stress (vertical fracture), and the fractures broke through to sublayer ⑥ (Fig. 12f); however, increasing the viscosity of the fracturing fluid greatly reduced the opening degree of the bedding fractures and the complexity of the fractures.

The aim of the hydraulic fracturing in the shale reservoir is to improve the complexity of fractures and increase the volume of stimulated reservoir, so as to maximize the productivity. According to the simulation results, highly developed bedding fractures would increase the complexity of hydraulic fractures, but also restrict the extension of the fracture height and reduce the volume of stimulated reservoir. For the interval with developed bedding and easy to open, it is necessary to reduce the cluster number, increase the displacement, and increase the amount and proportion of high-viscosity fracturing fluid in the prepad fluid stage, to reduce the restriction of the bedding on fracture height expansion; for the interval with undeveloped bedding and difficult to open (upper gas layer with undeveloped bedding or deep shale gas reservoirs with large overlying stress), the target is to increase the fracture complexity and enlarge the gas drainage fracture area. In this context, it is necessary to increase the cluster number, increase the displacement, and reduce the viscosity of the fracturing fluid (using low-viscosity slick water).

Taking Well Jiaoye A in the Fuling shale gas field as an example, the horizontal section of the well mainly goes through sublayers ①-③. The well was fractured in a total of 28 sections (16 with 2 clusters and 12 with 3 clusters of perforations), at a cluster spacing of 17-25 m, average fluid consumption of single stage of 1965 m³ and proppant consumption of 50.7 m³. Micro seismic monitoring was carried out in the 28 sections during fracturing, and the number of micro seismic events in each stage was 5-135. The monitoring results show that the fractures are 229-483 m long, the fracture network is 34-204 m wide, and the fractures are 16-49 m high (52% of the fractures are 30-40 m high). The three-dimensional numerical model of fracture propagation was used to simulate each stage, and the results were compared with the results of micro-seismic monitoring (Fig. 13, the third section of monitoring is taken as an example for display). From the simulation, the fractures in the sections were 250-430 m long and 25-50 m high, and fracture network was 76-167 m wide, which are generally consistent with the results of micro-seismic monitoring.

Based on the discrete element method, a fracture propagation model of shale reservoir with developed bedding which considers the influence of the weak bedding plane and the longitudinal stress difference is established. The actual parameters of the horizontal shale gas well in the Longmaxi Formation of the main area of the Jiaoshiba anticline in the Fuling shale gas field were used for numerical simulation analysis. It is found that high-density beddings can increase the fracture complexity of the stimulated volume, but inhibit the fracture height and fracture length obviously. When the low-strength bedding is opened, artificial fractures change to expand horizontally, resulting in low bottom hole pressure and limited fracture height expansion. When the bedding effect is not considered, and there is a 7 MPa interlayer stress difference, some local hydraulic fractures in the middle cluster can pass through the stress shielding layer, reaching a height of close to 50 m. In contrast, when the bedding effect is considered, there is a 3-4 MPa interlayer stress difference, and the fractures are controlled in the stress shielding layer, and less than 38 m high. The influences of weak bedding plane and longitudinal stress difference should be considered when predicting fracture propagation in shale reservoirs, otherwise, the predicted result of fracture height would be quite different from the actual height. In fracturing of the interval with abundant bedding, the cluster number should be reduced, the displacement and the amount and proportion of high-viscosity fracturing fluid in the prepad fluid stage should be increased to reduce the restriction of bedding on the fracture height expansion, so as to improve the vertical extension of fractures and increase the volume of stimulated reservoirs.

a_i,t—acceleration of gravity, m/s²;

A—contact area, m²;

A—flexibility coefficient, Pa^-1;

b_i—force per unit volume, N/m³;

d_BP—weak plane spacing, m;

D_ijkl—elastic tensor, Pa;

E—elastic modulus in isotropic and homogeneous reservoir, Pa;

E_h, E_v—elasticity moduli parallel to and perpendicular to the bedding direction, Pa;

f—right vector, m/s;

F_n—normal force, N;

F_{n, max}—the normal force required for tension break of spring between block elements, Pa;

F_s—tangential force, N;

F_{s, max}—tangential force required for shear break of spring between block elements, Pa;

F_{s, re}—residual shear resistance after spring break, N;

G_v—shear modulus, Pa;

GR—gamma ray, API;

i, j, k, l—tensor indexes, dimensionless;

k_n and k_s—normal and tangential spring stiffness, N/m;

K₀—matrix permeability, m²;

m, m+1—current iteration and next iteration;

n-1, n—current time step and next time step;

n_j—unit vector normal to the fracture surface, dimensionless;

N—number of fractures;

N_c—number of clusters;

p—pressure, Pa;

p_i—fluid pressure exerted on the fracture wall, Pa;

p_m—experimental solution of pressure from the current iteration step, Pa;

p_m+1/2—experimental solution of pressure from the current iteration step to the next iteration step, Pa;

p₀—pore pressure, Pa;

q_l—filtration rate, m/s;

q_x, q_y—velocity in x and y directions, m²/s;

Q—total injection volume, m³/s;

R_LLD—deep lateral resistivity, Ω·m;

R_LLS—shallow lateral resistivity, Ω·m;

s—fracture area, m²;

S_BP—shear strength of the bedding fracture, MPa;

S₀—matrix shear strength, Pa;

t—time, s;

T_BP—tensile strength of the bedding fracture, MPa;

T₀—tensile strength of the matrix, Pa;

w—dynamic crack width, m;

w₀—initial fracture width, m;

w_m—experimental solution of fracture width of current iteration step, m;

x, y—two directions of the coordinate system, m;

u_i—displacement, m;

u_i,t—velocity, m/s;

α—damping per unit volume, (kg·s)/m³;

β—iterative correction coefficient;

γ—experimental solution coefficient, m/(Pa·s);

ε_kl—strain tensor, dimensionless;

μ—fluid viscosity, Pa·s;

υ—Poisson's ratio, dimensionless;

υ_h, υ_v—the Poisson's ratios in the directions of parallel to bedding and vertical to bedding, dimensionless;

ρ—rock density, kg/m³;

σ_h—minimum horizontal principal stress, Pa;

σ_H—maximum horizontal principal stress, Pa;

σ_ij—Cauchy tensor, Pa;

σ_{ij, j}—Cauchy Tensor derivative, N/m³;

σ_v—vertical stress, Pa;

φ—internal friction angle, (°);

φ_BP—internal friction angle of bedding fracture, (°);

Ω_f—solution domain;

Δt_m—size of an adaptive time step in Step m, s;

Δw_m—fracture width change of current iteration step, m;

Δu_n, Δu_s—normal and tangential relative displacement between adjacent nodes, m.