Numerical Study on Hydraulic Fracture Propagation in a Layered Continental Shale Reservoir

: The distribution of beddings varies greatly in shale reservoirs. The inﬂuence of beddings on hydraulic fracture propagation has often been studied using simpliﬁed geological models, i


Introduction
Multistage hydraulic fracturing technology plays an important role in shale gas production.Engineering practice finds that geological factors, such as bedding distributions, interface properties and in situ stresses, and engineering factors, such as fluid injection rate and fluid viscosity, all affect fracture complexity [1][2][3][4][5].Field hydraulic fracturing tests reveal that it is difficult to obtain expected fracture heights through increasing fluid injection rates due to the blockage of beddings [6,7].The distribution of beddings varies greatly in the vertical direction, leading to the heterogeneity and mechanical anisotropy of shale reservoirs [8][9][10][11][12].Fracture growth in layered shale reservoirs has the following characteristics [13][14][15]: (1) they are affected by beddings and fracture height growth is restricted; (2) fracture geometry is complex, resulting in uncertainty in multistage hydraulic fracturing design.
Fracture propagation in layered reservoirs has been often studied by hydraulic fracturing experiments.Warpinski et al. [16] investigated the effect of stress variations on fracture growth in layered rock samples.Olson et al. [17], Liu et al. [18] and Huang et al. [19] investigated fracture height growth in layered rocks based on experimental analysis.Ranjith et al. [20] found that a low foam injection flow rate could improve fracture complexity.Triaxial fracturing tests provide an effective way to analyze the effect of geological and engineering parameters on fracture penetration [21][22][23][24][25]. Zhang et al. [26] and Zhao et al. [15] analyzed the role of beddings in facture growth, and concluded that bedding interfaces Energies 2022, 15, 8840 2 of 18 were the main cause of fracture network patterns.The interaction mechanisms between hydraulic fractures and bedding interfaces are well understood through laboratory experiments.However, limited by sample size, knowledge of the actual propagation process on a reservoir scale remains unclear.
Numerical simulations provide an insight into fracture growth in layered reservoirs [27].Commonly used methods for hydraulic fracturing include the finite element method (FEM), displacement discontinuity method (DDM), discrete element method (DEM) and extended finite element method (XFEM).For example, Fu et al. [28] investigated the influence of interlayer permeability and stress difference on fracture penetration using the FEM.Gu et al. [29] built a three-dimensional model, considering the influence of bedding plane interfacial slip and filtration loss based on the DDM; however, the interaction between bedding planes and hydraulic fractures was solved by a simplified analytical method.Zhang et al. [30] established a quasi-3D displacement discontinuity model to study fracture height growth in layered rocks.Tang et al. [31][32][33] and Xie et al. [34] used the DDM to study the limiting effect of horizontal weak interfaces on fracture height.Dahi-Taleghani and Olson [35], and Klimenko and Dahi-Taleghani [36] realized the simulation of asymmetric fracture wings using the XFEM.Settgast et al. [37] presented a fully coupled finite element/finite volume method to realize a filed-scale reservoir simulation, and found that under low vertical stress difference, T-shaped or horizontal fractures were easily formed, and the opening of bedding interfaces inhibited the growth of fracture height.Tan et al. [38] modeled the turning, twisting and expanding behavior of hydraulic fractures in a layered shale formation using the XFEM-based CZM method.The discrete element method (DEM) offers an effective tool to model arbitrary propagation of hydraulic fractures, and is more suitable for analyzing discontinuous media, such as fractured reservoirs [39,40].Currently, the DEM has been successfully applied to study the role of preexisting discontinuities on fracture extension.For example, Zhou et al. [41] and Zhang et al. [42] employed a smooth joint model to study interactions between hydraulic and natural fractures; Yao et al. [43] investigated the role of natural fracture characteristics in hydraulic fracture propagation using the DEM; Li et al. [44] used a grain-based DEM to model hydraulic fracture propagation.Although current progress in hydraulic fracturing modeling has been achieved, knowledge of the actual propagation process remains unclear, since the geological models used in the simulations cannot reflect the geological structure of shale reservoirs [3].Fracture behavior cannot be predicted or controlled unless an accurate geological model is built.
The Ordos Basin is the second largest sedimentary basin in China, where multiscale beddings are pervasive in the continental shale formations [45,46].At present, field fracturing tests in this area show that injected fluid is prone to enter bedding interfaces during hydraulic fracturing, resulting in unsatisfactory production and low economic benefits [6].Although the combination of high-viscosity fluids which is used to generate hydraulic fractures and low-viscosity slick water which is used to disrupt beddings has been used, the reconstruction effect of different blocks is quite different [6,47].To effectively guide the hydraulic-fracturing design, it is urgent to determine fracture growth in layered shale rock, and clarify the influencing factors.In this research, an outcrop-based bedding model of the continental shale formation in the Ordos Basin is built.Reservoir-scale simulations using the discrete element method are employed to investigate the influence of geological factors and engineering factors on fracture geometry.

Geological Background
The Ordos Basin is a multi-cycle sedimentary basin that developed on the North China Platform.The structural framework of the basin was formed in the Mesozoic period.The basin is a gentle monoclinic structure, low in the west and high in the east, and the east-to-west dip is smaller than 1 • [48].It is composed of six secondary structural units [14], including the Yimeng uplift belt, Yishan slope, Tianhuan depression, western edge thrusting belt, Jinxi flexural fold belt and Weibei uplift belt (Figure 1a).[14], including the Yimeng uplift belt, Yishan slope, Tianhuan depression, western edge thrusting belt, Jinxi flexural fold belt and Weibei uplift belt (Figure 1a).In the late Triassic period, as the material supply in the basin was sufficient, a series of terrigenous clastic deposits were formed.The Chang 7 member of the Triassic Yanchang Formation was formed in the heyday of the lake basin expansion, which is divided into three sub-members from top to bottom, namely Chang 7-1, Chang 7-2 and Chang 7-3.The Chang 7-2 shale, with a thickness varying from 10 m to 40 m, is the main source rock, which is deemed to be the key target formation for gas exploration.As a sedimentary rock, the structure of the Chang 7-2 shale is thin-layered or laminar, resulting in strong heterogeneity in its mechanical properties.The porosity of shale is approximately 10%, and the permeability is less than 1 × 10 −3 μm 2 .Filed investigation was carried out in the south margin of the Ordos Basin, as shown in Figure 1b.

Geological Model of the Shale Gas Reservoir
Field investigation, borehole imaging and logging were conducted to obtain the bedding data in the study area (Figure 2a,b).Figure 2c shows that the beddings have an average spacing of 20.9 cm.As shown in Figure 2d, the distribution of the spacing obeys a negative exponential distribution.Bedding planes are horizontally distributed in the shale formations in the study area.In the model, bedding planes are horizontal planes and the persistence of the bedding planes is assumed to be larger than the model size.A threedimensional (3D) bedding model of the Chang 7-2 shale is built using the Monte Carlo method based on the statistical parameters of bedding spacing, including the mean of 20.9 cm and standard deviation of 25.2 cm (Figure 3).The average spacing of the simulated beddings is 20.7 cm.The relative difference between the simulated and actual beddings is only 0.96%, indicating that the bedding model is reliable.In the late Triassic period, as the material supply in the basin was sufficient, a series of terrigenous clastic deposits were formed.The Chang 7 member of the Triassic Yanchang Formation was formed in the heyday of the lake basin expansion, which is divided into three sub-members from top to bottom, namely Chang 7-1, Chang 7-2 and Chang 7-3.The Chang 7-2 shale, with a thickness varying from 10 m to 40 m, is the main source rock, which is deemed to be the key target formation for gas exploration.As a sedimentary rock, the structure of the Chang 7-2 shale is thin-layered or laminar, resulting in strong heterogeneity in its mechanical properties.The porosity of shale is approximately 10%, and the permeability is less than 1 × 10 −3 µm 2 .Filed investigation was carried out in the south margin of the Ordos Basin, as shown in Figure 1b.

Geological Model of the Shale Gas Reservoir
Field investigation, borehole imaging and logging were conducted to obtain the bedding data in the study area (Figure 2a,b).Figure 2c shows that the beddings have an average spacing of 20.9 cm.As shown in Figure 2d, the distribution of the spacing obeys a negative exponential distribution.Bedding planes are horizontally distributed in the shale formations in the study area.In the model, bedding planes are horizontal planes and the persistence of the bedding planes is assumed to be larger than the model size.A three-dimensional (3D) bedding model of the Chang 7-2 shale is built using the Monte Carlo method based on the statistical parameters of bedding spacing, including the mean of 20.9 cm and standard deviation of 25.2 cm (Figure 3).The average spacing of the simulated beddings is 20.7 cm.The relative difference between the simulated and actual beddings is only 0.96%, indicating that the bedding model is reliable.

Numerical Method 4.1. Governing Equation
In the numerical method used in this study, rock mass is represented by an assembly of intact rock blocks separated by a fracture network.Crack propagation is achieved by contact breakage [50].The contact constitutive behavior is described by [51,52].
where σ n and τ s are the normal and shear stresses, respectively; u n and u s are the normal and shear displacements, respectively; ∆u s is the displacement increment; k n and k s are the normal and shear stiffnesses, respectively; ϕ and c are the friction angle and cohesive strength, respectively; and L is the length of the contact surface.
The rock blocks are assumed to be impermeable, and fluid infiltrates the rock only through fractures.The flow rate q f along a fracture is calculated by the following cubic law [43,53,54]: where e is the aperture, ∆p is the pressure differential, and l and w are the length and width of a crack segment, respectively.e is calculated by [50].
where e 0 is the initial domain aperture, and ∆u n is the domain normal displacement increment.The domain pressure p is calculated as shown in [52].
where p 0 is the initial domain pressure, Q is the sum of flow rates into the domain, K w is the fluid bulk modulus, ∆t is the time step, and V and V 0 are the new and old domain volumes, respectively.

Validation
We choose a homogeneous and elastic model with a Young's modulus of 10.3 GPa and a Poisson's ratio of 0.22 (Figure 4a) to verify the numerical method through a PKN solution [55].A fluid with a 0.001 Pa•s viscosity is injected through the injection point at a rate of 0.002 m 3 /s.The fracture aperture distribution calculated by the simulation matched well with that obtained by the PKN solution (Figure 4b,c).The maximum fracture apertures obtained by the numerical and PKN solutions are 0.98 mm and 0.94 mm, respectively.

Numerical Model and Parameter Settings
Numerical simulations of fracture propagation in the layered shale reservoir are conducted based on the 3D bedding model.To improve computational efficiency, a cubic model with a length of 10 m (Figure 5a) is extracted from the bedding model and used as the computational model in this work.In the simulations, to avoid the boundary effect, a larger model region with the dimensions of 20 m × 20 m × 20 m that is centered at (10, 10, 10) is built, while the inner box (i.e., the area of interest in the bedding model) with the dimensions of 10 m × 10 m × 10 m, which is also centered at (10, 10, 10), is built.A vertical profile is chosen to monitor the fracture growth process, as shown in Figure 5b.The isotropic elastic model is used for the rock matrix, and the Coulomb slip model is used for the beddings.The parameters used in the simulations are obtained from laboratory experiments on shale rock carried out by previous studies [11,56], as listed in Table 1.The injection well with a radius of 0.1 m is placed in the model center, and the injection time is 0.7 s.The model boundaries are fixed by setting the velocity to zero.In the simulation, the rock matrix is deemed to be impermeable, and fluid only flows along the fractures.In addition, one must note that fluid leak-off is not considered and the computation is timeconsuming.5a) is extracted from the bedding model and used as the computational model in this work.In the simulations, to avoid the boundary effect, a larger model region with the dimensions of 20 m × 20 m × 20 m that is centered at (10, 10, 10) is built, while the inner box (i.e., the area of interest in the bedding model) with the dimensions of 10 m × 10 m × 10 m, which is also centered at (10, 10, 10), is built.A vertical profile is chosen to monitor the fracture growth process, as shown in Figure 5b.The isotropic elastic model is used for the rock matrix, and the Coulomb slip model is used for the beddings.The parameters used in the simulations are obtained from laboratory experiments on shale rock carried out by previous studies [11,56], as listed in Table 1.The injection well with a radius of 0.1 m is placed in the model center, and the injection time is 0.7 s.The model boundaries are fixed by setting the velocity to zero.In the simulation, the rock matrix is deemed to be impermeable, and fluid only flows along the fractures.In addition, one must note that fluid leak-off is not considered and the computation is time-consuming.The vertical stress difference strongly influences hydraulic fracture geometry, shown in Figures 6 and 7. When Δσ = 0, the vertical growth of hydraulic fractures is r stricted by bedding interfaces (Figure 6a).The opening of the bedding interfaces strong inhibits the fracture height growth.Hydraulic fractures preferentially expand horizo tally along the beddings when Δσ = 0, resulting in an ellipse-shaped fracture geometry, shown in Figure 6a.Hydraulic fractures are more likely to be captured by bedding inte faces when Δσ is low.When Δσ increases from 0 MPa to 10 MPa, hydraulic fractures co tinue to expand upward, but fail to break through the densely distributed beddings in th lower part of the model (Figure 6b).When Δσ = 20 MPa and Δσ = 30 MPa, hydraulic fra tures propagate through the beddings in both upward and downward directions drive by the vertical stress; the fracture height increases rapidly, and the fracture shape chang from an ellipse to a circle (Figure 6c,d), indicating that high vertical stress difference favorable for fracture penetration through bedding interfaces.As shown in Figure 8, hig vertical stress difference can improve fracture height and fracture area.The vertical stress difference strongly influences hydraulic fracture geometry, as shown in Figures 6 and 7. When ∆σ = 0, the vertical growth of hydraulic fractures is restricted by bedding interfaces (Figure 6a).The opening of the bedding interfaces strongly inhibits the fracture height growth.Hydraulic fractures preferentially expand horizontally along the beddings when ∆σ = 0, resulting in an ellipse-shaped fracture geometry, as shown in Figure 6a.Hydraulic fractures are more likely to be captured by bedding interfaces when ∆σ is low.When ∆σ increases from 0 MPa to 10 MPa, hydraulic fractures continue to expand upward, but fail to break through the densely distributed beddings in the lower part of the model (Figure 6b).When ∆σ = 20 MPa and ∆σ = 30 MPa, hydraulic fractures propagate through the beddings in both upward and downward directions driven by the vertical stress; the fracture height increases rapidly, and the fracture shape changes from an ellipse to a circle (Figure 6c,d), indicating that high vertical stress difference is favorable for fracture penetration through bedding interfaces.As shown in Figure 8, high vertical stress difference can improve fracture height and fracture area.

Effect of Fluid Viscosity on Fracture Complexity
Four kinds of fluid viscosity, i.e., μ = 0.001 Pa•s, 0.01 Pa•s, 0.1 Pa•s and 1 Pa•s, are applied in the simulations.The vertical/horizontal stress is 30/20 MPa, and the injection rate is 0.05 m 3 /s, and the bedding aperture is 0.1 mm.The results suggest that fracture height and fracture area both decrease with the increase in fluid viscosity (Figures 9 and  10).Fracture propagation is significantly affected by beddings for the case of μ= 0.001 Pa•s, which is characterized by horizontal expansion along the beddings.The hydraulic fracture height is controlled by the densely distributed beddings, and the fracture shape is an ellipse (Figure 9a).With the increase in fluid viscosity, fracture propagation in the vertical and horizontal directions is blocked (Figure 9b-d); at this time, fluid viscosity controls the propagation pattern of hydraulic fractures.According to Figure 11, it seems that increasing fluid viscosity can restrict fracture growth.In conclusion, high fluid viscosity increases the resistance of hydraulic fractures to propagate, thereby reducing the height and length of hydraulic fractures.

Effect of Fluid Viscosity on Fracture Complexity
Four kinds of fluid viscosity, i.e., μ = 0.001 Pa•s, 0.01 Pa•s, 0.1 Pa•s and 1 Pa•s, are applied in the simulations.The vertical/horizontal stress is 30/20 MPa, and the injection rate is 0.05 m 3 /s, and the bedding aperture is 0.1 mm.The results suggest that fracture height and fracture area both decrease with the increase in fluid viscosity (Figures 9 and  10).Fracture propagation is significantly affected by beddings for the case of μ= 0.001 Pa•s, which is characterized by horizontal expansion along the beddings.The hydraulic fracture height is controlled by the densely distributed beddings, and the fracture shape is an ellipse (Figure 9a).With the increase in fluid viscosity, fracture propagation in the vertical and horizontal directions is blocked (Figure 9b-d); at this time, fluid viscosity controls the propagation pattern of hydraulic fractures.According to Figure 11, it seems that increasing fluid viscosity can restrict fracture growth.In conclusion, high fluid viscosity increases the resistance of hydraulic fractures to propagate, thereby reducing the height and length of hydraulic fractures.

Effect of Fluid Viscosity on Fracture Complexity
Four kinds of fluid viscosity, i.e., µ = 0.001 Pa•s, 0.01 Pa•s, 0.1 Pa•s and 1 Pa•s, are applied in the simulations.The vertical/horizontal stress is 30/20 MPa, and the injection rate is 0.05 m 3 /s, and the bedding aperture is 0.1 mm.The results suggest that fracture height and fracture area both decrease with the increase in fluid viscosity (Figures 9 and 10).Fracture propagation is significantly affected by beddings for the case of µ= 0.001 Pa•s, which is characterized by horizontal expansion along the beddings.The hydraulic fracture height is controlled by the densely distributed beddings, and the fracture shape is an ellipse (Figure 9a).With the increase in fluid viscosity, fracture propagation in the vertical and horizontal directions is blocked (Figure 9b-d); at this time, fluid viscosity controls the propagation pattern of hydraulic fractures.According to Figure 11, it seems that increasing fluid viscosity can restrict fracture growth.In conclusion, high fluid viscosity increases the resistance of hydraulic fractures to propagate, thereby reducing the height and length of hydraulic fractures.

Effect of Injection Rate on Fracture Complexity
In the simulation, the injection rate q is set as 0.02 m 3 /s, 0.04 m 3 /s, 0.06 m 3 /s and 0.08 m 3 /s, respectively, to investigate its influence on fracture complexity.The bedding aperture is 0.1 mm and the vertical/horizontal stress is 30/20 MPa.The fluid viscosity is set as 0.001 Pa•s.
When q = 0.02 m 3 /s, fracture height growth is strongly restricted by the beddings, and the fracture shape is elliptical (Figure 12a).When q increases from 0.04 to 0.08 m 3 /s, the fracture height increases (Figure 12b-d

Effect of Injection Rate on Fracture Complexity
In the simulation, the injection rate q is set as 0.02 m 3 /s, 0.04 m 3 /s, 0.06 m 3 /s and 0.08 m 3 /s, respectively, to investigate its influence on fracture complexity.The bedding aperture is 0.1 mm and the vertical/horizontal stress is 30/20 MPa.The fluid viscosity is set as 0.001 Pa•s.
When q = 0.02 m 3 /s, fracture height growth is strongly restricted by the beddings, and the fracture shape is elliptical (Figure 12a).When q increases from 0.04 to 0.08 m 3 /s, the fracture height increases (Figure 12b-d

Effect of Injection Rate on Fracture Complexity
In the simulation, the injection rate q is set as 0.02 m 3 /s, 0.04 m 3 /s, 0.06 m 3 /s and 0.08 m 3 /s, respectively, to investigate its influence on fracture complexity.The bedding aperture is 0.1 mm and the vertical/horizontal stress is 30/20 MPa.The fluid viscosity is set as 0.001 Pa•s.
When q = 0.02 m 3 /s, fracture height growth is strongly restricted by the beddings, and the fracture shape is elliptical (Figure 12a).When q increases from 0.04 to 0.08 m 3 /s, the fracture height increases (Figure 12b-d) because the hydraulic fractures break through the beddings in the lower part of the model and expand downward.However, hydraulic fractures cannot penetrate the densely distributed beddings in the upper part of the model.Fluid pressure distributions in the model with different injection rates are shown in Figure 13.q = 0.02 m 3 /s q = 0.04 m 3 /s q = 0.06 m 3 /s q = 0.08 m 3 /s q = 0.02 m 3 /s q = 0.04 m 3 /s q = 0.02 m 3 /s q = 0.04 m 3 /s q = 0.06 m 3 /s q = 0.08 m 3 /s q = 0.02 m 3 /s q = 0.04 m 3 /s  The injection rate can improve the fracture height and the fracture area (Figure 14).Usually, fluid pressure has a positive relationship with the fluid injection rate.Fluid pressure can generate induced tensile stress at the tip of a hydraulic fracture, which can cause the failure of rock matrices and shearing of bedding interfaces.Therefore, hydraulic fractures can easily penetrate through bedding interfaces under the condition of a high injection rate.In summary, the injection rate significantly affects fracture complexity, and a higher injection rate can yield a complex hydraulic fracture network.The simulation results are shown in Figures 15 and 16.Affected by the beddings, the vertical expansion of hydraulic fractures is constrained, and they mainly propagate along the bedding planes in the horizontal direction, resulting in an ellipse-shaped fracture geometry (Figure 15a).As the bedding aperture increases, the fracture height decreases but the fracture width increases (Figure 15b-d), resulting in a slight increase in the fracture area (Figure 17).The injection rate can improve the fracture height and the fracture area (Figure 14).Usually, fluid pressure has a positive relationship with the fluid injection rate.Fluid pressure can generate induced tensile stress at the tip of a hydraulic fracture, which can cause the failure of rock matrices and shearing of bedding interfaces.Therefore, hydraulic fractures can easily penetrate through bedding interfaces under the condition of a high injection rate.In summary, the injection rate significantly affects fracture complexity, and a higher injection rate can yield a complex hydraulic fracture network.The injection rate can improve the fracture height and the fracture area (Figure 14).Usually, fluid pressure has a positive relationship with the fluid injection rate.Fluid pressure can generate induced tensile stress at the tip of a hydraulic fracture, which can cause the failure of rock matrices and shearing of bedding interfaces.Therefore, hydraulic fractures can easily penetrate through bedding interfaces under the condition of a high injection rate.In summary, the injection rate significantly affects fracture complexity, and a higher injection rate can yield a complex hydraulic fracture network.The simulation results are shown in Figures 15 and 16.Affected by the beddings, the vertical expansion of hydraulic fractures is constrained, and they mainly propagate along the bedding planes in the horizontal direction, resulting in an ellipse-shaped fracture geometry (Figure 15a).As the bedding aperture increases, the fracture height decreases but the fracture width increases (Figure 15b-d), resulting in a slight increase in the fracture area (Figure 17).The simulation results are shown in Figures 15 and 16.Affected by the beddings, the vertical expansion of hydraulic fractures is constrained, and they mainly propagate along the bedding planes in the horizontal direction, resulting in an ellipse-shaped fracture geometry (Figure 15a).As the bedding aperture increases, the fracture height decreases but the fracture width increases (Figure 15b-d), resulting in a slight increase in the fracture area (Figure 17).

Conclusions
In this work, an outcrop-data-based bedding model of a continental shale formation in the Ordos Basin, China, is built.The discrete element method is used to carry out a series of numerical simulations of hydraulic fracture propagation to investigate the influence of geological factors and engineering factors on hydraulic fracture complexity.The main conclusions are as follows: (1) Beddings have significant inhibitory effects on fracture height growth.It is difficult for hydraulic fractures to penetrate zones with more densely distributed beddings.When a hydraulic fracture encounters a bedding with a larger aperture, it is more likely to be captured and expand along the weak interface.
(2) High vertical stress difference is favorable for the penetration of hydraulic fractures through bedding interfaces.When the vertical stress difference is Δσ =0 MPa and Δσ =10 MPa, the vertical growth of hydraulic fractures is strongly restricted by bedding interfaces.When Δσ = 20 MPa and Δσ = 30 MPa, hydraulic fractures propagate through the beddings in both upward and downward directions driven by the vertical stress.Thus, increasing vertical stress difference can improve fracture complexity.
(3) The properties of the injected fluid significantly affect the fracture propagation process.When the injection rate increases from 0.02 m 3 /s to 0.08 m 3 /s, the fracture height increases, suggesting that a higher injection rate can produce a more complex hydraulic fracture geometry.The fracture height and fracture area both decrease when the fluid

Conclusions
In this work, an outcrop-data-based bedding model of a continental shale formation in the Ordos Basin, China, is built.The discrete element method is used to carry out a series of numerical simulations of hydraulic fracture propagation to investigate the influence of geological factors and engineering factors on hydraulic fracture complexity.The main conclusions are as follows: (1) Beddings have significant inhibitory effects on fracture height growth.It is difficult for hydraulic fractures to penetrate zones with more densely distributed beddings.When a hydraulic fracture encounters a bedding with a larger aperture, it is more likely to be captured and expand along the weak interface.
(2) High vertical stress difference is favorable for the penetration of hydraulic fractures through bedding interfaces.When the vertical stress difference is Δσ =0 MPa and Δσ =10 MPa, the vertical growth of hydraulic fractures is strongly restricted by bedding interfaces.When Δσ = 20 MPa and Δσ = 30 MPa, hydraulic fractures propagate through the beddings in both upward and downward directions driven by the vertical stress.Thus, increasing vertical stress difference can improve fracture complexity.
(3) The properties of the injected fluid significantly affect the fracture propagation process.When the injection rate increases from 0.02 m 3 /s to 0.08 m 3 /s, the fracture height increases, suggesting that a higher injection rate can produce a more complex hydraulic fracture geometry.The fracture height and fracture area both decrease when the fluid

Conclusions
In this work, an outcrop-data-based bedding model of a continental shale formation in the Ordos Basin, China, is built.The discrete element method is used to carry out a series of numerical simulations of hydraulic fracture propagation to investigate the influence of geological factors and engineering factors on hydraulic fracture complexity.The main conclusions are as follows: (1) Beddings have significant inhibitory effects on fracture height growth.It is difficult for hydraulic fractures to penetrate zones with more densely distributed beddings.When a hydraulic fracture encounters a bedding with a larger aperture, it is more likely to be captured and expand along the weak interface.
(2) High vertical stress difference is favorable for the penetration of hydraulic fractures through bedding interfaces.When the vertical stress difference is ∆σ = 0 MPa and ∆σ = 10 MPa, the vertical growth of hydraulic fractures is strongly restricted by bedding interfaces.When ∆σ = 20 MPa and ∆σ = 30 MPa, hydraulic fractures propagate through the beddings in both upward and downward directions driven by the vertical stress.Thus, increasing vertical stress difference can improve fracture complexity.
(3) The properties of the injected fluid significantly affect the fracture propagation process.When the injection rate increases from 0.02 m 3 /s to 0.08 m 3 /s, the fracture height increases, suggesting that a higher injection rate can produce a more complex hydraulic fracture geometry.The fracture height and fracture area both decrease when the fluid

Figure 2 .
Figure 2. Characteristics of beddings in the Chang 7-2 shale: (a) distribution of beddings at an outcrop; (b) downhole camera image from a well [45]; (c) spacing of beddings based on a total of 90 mapped layers; (d) distribution of spacing of beddings.

Figure 2 .Figure 2 .
Figure 2. Characteristics of beddings in the Chang 7-2 shale: (a) distribution of beddings at an outcrop; (b) downhole camera image from a well [45]; (c) spacing of beddings based on a total of 90 mapped layers; (d) distribution of spacing of beddings.

Figure 4 .
Figure 4. (a) The validation model; (b) fracture aperture distributions calculated by the numerical and PKN solutions; (c) fracture comparison between the fracture apertures calculated by the numerical and PKN solutions.

Figure 4 .
Figure 4. (a) The validation model; (b) fracture aperture distributions calculated by the numerical and PKN solutions; (c) fracture comparison between the fracture apertures calculated by the numerical and PKN solutions.

5 .
Hydraulic Fracture Propagation in the Layered Shale Reservoir 5.1.Numerical Model and Parameter Settings Numerical simulations of fracture propagation in the layered shale reservoir are conducted based on the 3D bedding model.To improve computational efficiency, a cubic model with a length of 10 m (Figure

Figure 5 .
Figure 5. Numerical model used for hydraulic fracturing simulations: (a) 3D bedding model wi an edge length of 10 m; (b) x-z plane with y = 10 m cut from the bedding model.

Figure 5 .
Figure 5. Numerical model used for hydraulic fracturing simulations: (a) 3D bedding model with an edge length of 10 m; (b) x-z plane with y = 10 m cut from the bedding model.

Figure 8 .
Figure 8. Variation in fracture height and fracture area versus vertical stress.

Figure 7 . 18 Figure 7 .
Figure 7. Fluid pressure distributions in different cases of stress configurations.

Figure 8 .
Figure 8. Variation in fracture height and fracture area versus vertical stress.

Figure 8 .
Figure 8. Variation in fracture height and fracture area versus vertical stress.

Figure 11 .
Figure 11.Variation in fracture height and fracture area versus fluid viscosity.

Figure 11 .
Figure 11.Variation in fracture height and fracture area versus fluid viscosity.

Figure 11 .
Figure 11.Variation in fracture height and fracture area versus fluid viscosity.

Figure 13 .
Figure 13.Fluid pressure distributions in different cases, with the injection rate varying from q = 0.02 m 3 /s to 0.04 m 3 /s, 0.06 m 3 /s and 0.08 m 3 /s.

Figure 14 .
Figure 14.Variation in fracture height and fracture area versus injection rate.5.2.4.Effect of Bedding Aperture on Fracture Complexity Fluid flow in shale reservoirs is closely related to bedding aperture.Four kinds of bedding aperture, including w = 0.1 mm, 0.2 mm, 0.3 mm and 0.4 mm, are applied to the model.The vertical/horizontal stress is 30/20 MPa, the fluid viscosity is 0.001 Pa•s, and the injection rate is 0.05 m 3 /s.The simulation results are shown in Figures15 and 16.Affected by the beddings, the vertical expansion of hydraulic fractures is constrained, and they mainly propagate along the bedding planes in the horizontal direction, resulting in an ellipse-shaped fracture geometry (Figure15a).As the bedding aperture increases, the fracture height decreases but the fracture width increases (Figure15b-d), resulting in a slight increase in the fracture area (Figure17).

Figure 13 .
Figure 13.Fluid pressure distributions in different cases, with the injection rate varying from q = 0.02 m 3 /s to 0.04 m 3 /s, 0.06 m 3 /s and 0.08 m 3 /s.

Energies 2022 , 18 Figure 13 .
Figure 13.Fluid pressure distributions in different cases, with the injection rate varying from q = 0.02 m 3 /s to 0.04 m 3 /s, 0.06 m 3 /s and 0.08 m 3 /s.

Figure 14 .
Figure 14.Variation in fracture height and fracture area versus injection rate.5.2.4.Effect of Bedding Aperture on Fracture Complexity Fluid flow in shale reservoirs is closely related to bedding aperture.Four kinds of bedding aperture, including w = 0.1 mm, 0.2 mm, 0.3 mm and 0.4 mm, are applied to the model.The vertical/horizontal stress is 30/20 MPa, the fluid viscosity is 0.001 Pa•s, and the injection rate is 0.05 m 3 /s.The simulation results are shown in Figures15 and 16.Affected by the beddings, the vertical expansion of hydraulic fractures is constrained, and they mainly propagate along the bedding planes in the horizontal direction, resulting in an ellipse-shaped fracture geometry (Figure15a).As the bedding aperture increases, the fracture height decreases but the fracture width increases (Figure15b-d), resulting in a slight increase in the fracture area (Figure17).

Figure 14 .
Figure 14.Variation in fracture height and fracture area versus injection rate.5.2.4.Effect of Bedding Aperture on Fracture Complexity Fluid flow in shale reservoirs is closely related to bedding aperture.Four kinds of bedding aperture, including w = 0.1 mm, 0.2 mm, 0.3 mm and 0.4 mm, are applied to the model.The vertical/horizontal stress is 30/20 MPa, the fluid viscosity is 0.001 Pa•s, and the injection rate is 0.05 m 3 /s.The simulation results are shown in Figures15 and 16.Affected by the beddings, the vertical expansion of hydraulic fractures is constrained, and they mainly propagate along the bedding planes in the horizontal direction, resulting in an ellipse-shaped fracture geometry (Figure15a).As the bedding aperture increases, the fracture height decreases but the fracture width increases (Figure15b-d), resulting in a slight increase in the fracture area (Figure17).

Figure 16 .
Figure 16.Fluid pressure distributions in different cases, with the bedding aperture varying from w = 0.1 mm to 0.2 mm, 0.3 mm and 0.4 mm.

Figure 17 .
Figure 17.Variation in fracture height and fracture area versus bedding aperture.

Figure 16 . 18 Figure 16 .
Figure 16.Fluid pressure distributions in different cases, with the bedding aperture varying from w = 0.1 mm to 0.2 mm, 0.3 mm and 0.4 mm.

Figure 17 .
Figure 17.Variation in fracture height and fracture area versus bedding aperture.

Figure 17 .
Figure 17.Variation in fracture height and fracture area versus bedding aperture.