Investigation on the Propagation Mechanisms of a Hydraulic Fracture in Glutenite Reservoirs Using DEM

: The geometry heterogeneity induced by embedded gravel can cause severe stress heterogeneity and strength heterogeneity in glutenite reservoirs, and subsequently affect the initiation and propagation of hydraulic fractures. Since the discrete element method (DEM) can accurately describe the inter-particle interactions, the macromechanical behavior of glutenite specimen can be preciously represented by DEM. Therefore, the initiation and propagation mechanisms of hydraulic fractures were investigated using a coupling seepage-DEM approach, the terminal fracture morphologies of hydraulic fractures were determined, and the effects of stress differences, permeability, and gravel strength were studied. The results show that the initiation and propagation of hydraulic fractures are signiﬁcantly inﬂuenced by embedded gravel. In addition, the stress heterogeneity and strength heterogeneity induced by the gravel embedded near the wellbore increase local initiation points, causing a complicated fracture network nearby. Moreover, due to the effect of local stress heterogeneity, gravel strength, and energy concentration near the fracture tip, four interactions of attraction, deﬂection, penetration, and termination between propagating fractures and encountering gravel were observed.


Introduction
In recent years, the utilization of oil and gas in the unconventional reservoirs [1,2], including gas shale [3,4], oil shale, and glutenite has attracted increased attention. Hydraulic fracturing is essential to increase production, which can form fractures with high conductivity, thus enhancing the ability of oil or gas to flow into the wells and increasing their productivity for cost-effective mining. Therefore, understanding the initiation mechanism and final patterns of fractures is of great significance for hydraulic fracturing design and treatment. However, due to severe geometric, stress, and strength heterogeneities, the initiation mechanism and propagation morphologies of hydraulic fractures are complicated [5].
Hydraulic fracturing involves complex seepage-stress-damage coupling processes. The geometric, stress, and strength heterogeneities inside the heterogeneous reservoirs seriously affect the fracture initiation and propagation of hydraulic fractures, and ultimately lead to complex fracture morphologies. Microseism, acoustic emission (AE), CT scanning, and nuclear magnetic resonance [6][7][8] have been widely used in the in situ monitoring and laboratory experiments to explore initiation and propagation laws of fractures. However, those techniques cannot explain the nucleation and evolution mechanism of fractures from the micro scale. Alternatively, numerical simulation technology can accurately reveal the micro-initiation and macro-propagation behaviors of hydraulic fractures [9]. In addition, the finite element method (FEM), boundary element method (BEM), and extended finite element method (XFEM) are broadly used. As the finite element cannot be re-meshed during simulation, only planar fractures are predicted [10]. Although BEM and XFEM can simulate the complex fracture morphologies, DEM cannot capture multiple simultaneous initiation points and the interaction between fractures in the propagation process [11].
In recent years, DEM has been broadly applied in the numerical simulation of hydraulic fracturing due to its unique advantages in describing the micromechanical response. Compared with continuous methods, DEM can accurately reveal the mechanics of microscopic particle by tracking and calculating particle interactions, in order that the nucleation, initiation, and propagation of hydraulic fractures can be truly reappeared. Meanwhile, DEM can also describe the complex fracture morphologies and accurately capture the multi-point initiation and propagation [12][13][14].
Glutenites are generally distributed in tight oil and gas reservoirs, and since glutenite reservoirs are usually deeply buried, their permeability and porosity are extremely low. For glutenite, the embedded gravel increases the complexity of initiation and propagation of hydraulic fractures in glutenite. True tri-axial experiments show that although the principal stress state dominates the propagation of hydraulic fractures in glutenite, gravel plays a key role in this process. When fractures encounter gravel, there are four interaction patterns: Attraction, deflection, penetration, and termination [15][16][17]. The geometric heterogeneity induced by gravel greatly adds the difficulty of numerical simulation of hydraulic fracturing in glutenite. The FSD model [5,18,19] is mainly used to characterize the influence of gravel on the initiation and propagation of hydraulic fractures. Although the true tri-axial experiment and numerical simulation explore the laws of fracture propagation in glutenite from different angles, it is still difficult to accurately characterize the rock mechanics behavior of glutenite and the seepage-stress-damage coupling problem involved in hydraulic fracturing.
Since DEM can directly represent the microstructure and mechanical characteristics of rock particles, the application of DEM to simulate the fracture initiation and propagation mechanism in glutenite is very adaptive and robust [20,21]. Nevertheless, few studies on DEM have been reported. Based on DEM, this paper explored the rock mechanics response from the micro scale, simulated the seepage process with the network flow model, and coupled the seepage with the discrete element process to investigate the microscopic mechanism of fracture initiation and propagation in glutenite. In addition, the effects of stress difference, permeability, gravel strength, injection rate, and well diameter were studied in this work.

Parallel Bond Model
The DEM simulates rock by bonding microparticles and characterizes the macromechanical behavior of the rock based on the interaction of microparticles. In this paper, a parallel bond model [22,23] is used to generate a numerical specimen of glutenite. After bonding, the particles can withstand normal force, shear force, and bending moment, which are calculated according to the following formula: where F n is the normal force between particles (in N); F n 0 is the normal force at last time step (in N); k n is the normal stiffness (in MPa/m); A is the bonding area between adjacent particles (in m 2 ); ∆δ n is the relative normal-displacement increment at the contact of adjacent particles (in m); F s is the shear force between the particles (in N); F s 0 is the shear force at last time step (in N); k s is the shear stiffness (in MPa/m); ∆δ s is the relative shear-displacement increment at the contact of adjacent particles (in m); M t is the inter-particle bending moment (in N · m); M t 0 is the bending moment at last time step (in N · m); J is the polar moment of inertia at the contact of adjacent particles (in m 4 ); and ∆θ t is the relative twist-rotation increment ( • ).

Failure Criterions
DEM can accurately represent the mechanical interaction between microparticles. According to the bonding relationship between the particles, it can be judged whether cracks are generated. In other words, when the tensile stress between the particles exceeds the tensile strength limit, the normal bonding between the two particles is broken, thereby forming tensile micro-crack; when the shear stress between the particles exceeds the shear strength limit, the two particles relatively slip, and shear micro-crack is formed. Therefore, the initiation criterion can be expressed by the following formula [21,24], and the microcrack whose direction is perpendicular to the particle contact direction and length is equal to the average of the adjacent particle sizes, is located at the failure of the bonding: where σ t is the tensile stress (in MPa); pb_ten is the tensile strength limit (in MPa); τ is the shear stress (in MPa); σ n is the normal compressive stress (in MPa); ϕ is the internal friction angle ( • ); and pb_coh is the cohesive (in MPa).

Implementation of Flow-DEM Coupling
Seepage involves a complex flow-stress coupling process. This paper simulates the seepage of fracturing fluid in porous media based on the modified network flow model [24] and couples it with the stress module characterized by discrete elements. Adjacent contact particles enclose a domain ( Figure 1, green lines), and neighboring domains are connected by channels whose positions are at the contact of adjacent particles and directions are perpendicular to the contact direction between the particles, as shown in Figure 1. Due to the pressure difference between adjacent domains, the fluid flows from one domain to another, causing its pressure rises. Accordingly, the particles start to move under the influence of the increased force, leading to the changes in local stress, which in turn changes the volume of the domain and its internal fluid pressure, and further affects the fluid flow. The fluid flow through the flow tube can be calculated according to the Hagen-Poiseuille formula [25]: where Q is the volume flow through the channel (in m 3 /s); w is the channel aperture (in m); µ is the fluid viscosity (in Pa · s); ∆p is the pressure difference between two adjacent domains (in Pa); and L is the channel length (in m). The amount of fluid pressure change in the domain can be calculated according to the following formula [26]: where ∆p is the pressure change caused by fluid exchange between the domains (in Pa); K is the fluid bulk modulus (in Pa); V d is the domain volume (in m 3 ); ∆t is the time step (in s); and ∆V d is the volume change of the domain (in m 3 ). The amount of fluid pressure change in the domain can be calculated according to the following formula [26]: where Δp is the pressure change caused by fluid exchange between the domains (in Pa); K is the fluid bulk modulus (in Pa); Vd is the domain volume (in m 3 ); Δt is the time step (in s); and ΔVd is the volume change of the domain (in m 3 ). The aperture of the channel is closely related to the bond between particles. In essence, the seepage-stress coupling causes the particles to move, which in turn changes the aperture of the channel. When the bond is broken, the channel aperture at the contact is equal to the actual distance between the particles. Meanwhile, when the particles are compressed, the channel aperture can be calculated according to the residual aperture (the width at which the particles only contact), which can be expressed as follows [24]: where w is the aperture of the channel under the compression of the particles (in m); w0 is the residual aperture of the channel (in m); F0 is the contact force of the particle when the aperture of the channel is equal to half of the residual aperture (in N); and F is the normal contact force of the particle (in N). Figure 1b shows the numerical model and boundary conditions of the glutenite. Considering the calculation cost, a laboratory-scale numerical model (20 × 20 cm) was built, which was filled with 10,394 particles, and adjacent particles were bonded to form rock  The aperture of the channel is closely related to the bond between particles. In essence, the seepage-stress coupling causes the particles to move, which in turn changes the aperture of the channel. When the bond is broken, the channel aperture at the contact is equal to the actual distance between the particles. Meanwhile, when the particles are compressed, the channel aperture can be calculated according to the residual aperture (the width at which the particles only contact), which can be expressed as follows [24]:

Simulation Scheme
where w is the aperture of the channel under the compression of the particles (in m); w 0 is the residual aperture of the channel (in m); F 0 is the contact force of the particle when the aperture of the channel is equal to half of the residual aperture (in N); and F is the normal contact force of the particle (in N). Figure 1b shows the numerical model and boundary conditions of the glutenite. Considering the calculation cost, a laboratory-scale numerical model (20 × 20 cm) was built, which was filled with 10,394 particles, and adjacent particles were bonded to form rock through the parallel bond model. The particle radii were 0.6 mm ≤ R ≤ 1.5 mm, according to the Rosin-Rammler distribution [27]. The gravel was created by the cutting method. In this method, polygons with different sizes, shapes, and positions were generated randomly and the particles surrounded by the polygons constitute gravel which was endowed with high mechanical parameters. A wellbore with a radius of 3 mm was established in the center of the model, and the boundary of constant injection rate was set to the wellbore. According to the servo principle, the external boundary motion of the glutenite was controlled to add confining pressure to the numerical model, and the horizontal principal stresses were set to σ 1 and σ 2 , respectively.

Simulation Scheme
In this work, the microscopic parameters used in the DEM simulation refer to those calibrated by Ju et al. [21], as shown in Table 1. The confining stresses are applied on the boundaries of the model under a servo-mechanism using four movable walls surrounding the rock specimen. The primary far-field horizontal geo-stress in x-direction is defined as σ 1 , while σ 2 refers to the secondary far-field horizontal geo-stress in y-direction. In addition, the validation of DEM approach on the simulation of fracture propagation in glutenite reservoirs was carried out in our previous work [26], and a good agreement was observed. Moreover, as a typical DEM solver, particle flow code (PFC2D) was applied to implement the simulation, while the seepage module was self-defined in the code. The geo-stress state is generally considered to control the direction of hydraulic fracture propagation. In accordance with the linear elastic theory, the hydraulic fracture mainly extends along the direction of the maximum principal stress, thus forming a bi-wing fracture [5]. However, for glutenite reservoirs, the presence of gravel leads to severe geometric and stress heterogeneities, affecting the initiation and propagation of fractures. To clarify this effect, the hydraulic fracturing processes in glutenite under four kinds of geo-stress differences were simulated. The applied geo-stress states are shown in Table 2. In addition, four permeabilities of 0.0005 × 10 −3 0.005 × 10 −3 , 0.05 × 10 −3 , and 0.1 × 10 −3 µm 2 were applied to evaluate its influence on the initiation and propagation of hydraulic fractures. The matrix-gravel strength difference would affect the interaction between hydraulic fractures and encountering gravel, and subsequently change the terminal fracture morphologies, thus two gravel strength ratios of 1.6 and 1.8 were used. Moreover, the influence of injection rate and well diameter was investigated, and two injection rates of 2 × 10 −4 m 3 /s and 4 × 10 −4 m 3 /s were used, while the initiation and propagation of hydraulic fractures under two well diameters (2 and 4 mm) were evaluated.

Effect of Geo-Stress Difference
According to the circumferential stress theory [21], when the circumferential stress on the well wall is greater than the bond strength limit, the bond between the particles breaks, thereby forming micro-cracks, namely, crack initiation. However, due to the local stress heterogeneity induced by the geometrical heterogeneity near the wellbore, it is easy to form multiple initiation points at the well wall, and the number of initiation points increases as the principal stress difference decreases. Figure 2 shows the fracture propagation patterns under different principal stress differences. The results show that crack initiation near the wellbore wall is controlled by the stress heterogeneity and the principal stress state. The gravel around the wellbore leads to local stress concentration, which enhances the stress heterogeneity, thus affecting the initiation and propagation of fractures, and it is dominant in influencing the crack initiation. It is clear from Figure 2 that near the wellbore wall, the fracture mainly extends to the concentrated area of gravel, but its macro propagation direction is still controlled by the principal stress state. Near the wellbore wall, the injected fluid pressure and fracture spread synchronously. The energy of fracture initiation and propagation mainly comes from the fluid pressure [20], which is sufficient to form the multi-branch cracks and even the broken zones, as shown in Figure 2c,d. As the injection pressure increases, the crack gradually extends forward. At a distance from the wellbore region, the fluid pressure spread lags far behind the fracture propagation, thus the local stress concentration controls the initiation and propagation of fractures, allowing for the fractures to mainly extend along the direction of the maximum principal stress.  Figure 3 shows the initiation and penetration times of fractures under different g stress differences. It is clear from the figure that as the stress difference increases, the i tiation and penetration times are greatly reduced. When σ1 = 3 MPa and σ2 = 5 MPa, gravel-induced stress concentration is mainly distributed in the direction of minimu principal stress, while the stress concentration caused by the principal stress is mai along the direction of the maximum principal stress. Therefore, the injected fluid simu When the propagating fracture encounters gravel, a complex interaction between the two occurs, thus affecting the propagation direction and final morphology of the fracture. Initially, there is a serious stress concentration at the gravel boundary. When the fracture approaches the gravel, micro-cracks appear at the gravel boundary, thereby attracting the main fracture to the gravel boundary, namely, the attraction. Subsequently, the propagating fracture interacts with the gravel. If the gravel strength is low, the fracture can continue to penetrate through the gravel; if the gravel strength is high, the fracture will extend along the gravel boundary, resulting in deflection; if the energy at the tip of the fracture is insufficient, the fracture is trapped by the gravel and the termination appears, as shown in Figure 2b. In addition, the interaction between the fracture and the gravel is affected by the approximate relation of the fracture. When the fracture deviates from the gravel, it tends to be attracted by the gravel. Nevertheless, when the fracture hits the gravel on-head during propagation, the gravel can prevent the fracture from propagating further. Moreover, when the energy of the fracture tip is high, the fracture can deflect from or penetrate the gravel to continue extending. Figure 3 shows the initiation and penetration times of fractures under different geostress differences. It is clear from the figure that as the stress difference increases, the initiation and penetration times are greatly reduced. When σ 1 = 3 MPa and σ 2 = 5 MPa, the gravel-induced stress concentration is mainly distributed in the direction of minimum principal stress, while the stress concentration caused by the principal stress is mainly along the direction of the maximum principal stress. Therefore, the injected fluid simultaneously percolates in multiple directions, which causes the injection pressure in the wellbore to increase at a slower speed, as shown in Figure 4, thus the fracture takes a long time to accumulate energy to initiate and propagate. As the principal stress difference increases, the range of the fluid percolation in the wellbore decreases, thereby accelerating the rate of pressure increase. It is clear from Figure 4 that as the stress difference increases, the wellbore injection pressure reaches the peak earlier, and the peak pressure is higher. During the fluid injection pressurization process, due to the frequent interaction of fractures and gravel, the injection pressure fluctuates and decreases rapidly after reaching the peak pressure. However, since the fracture propagates along the path of least resistance, even if the highstrength gravel blocks the fracture propagation, the hydraulic fracture can always find the minimum resistance path and continue propagating in the matrix, thus the initiation pressure is less affected by the stress difference.

Effect of Permeability
The permeability can directly reflect the fluid percolation and pressure Figure 5 shows the fracture propagation patterns under different permea the permeability is small, it is difficult for the injection fluid to seep into Therefore, the fluid pressure in the wellbore rises rapidly, and the spread of sure lags far behind, causing the surrounding stress to rapidly increas strength limit. Accordingly, the micro-cracks appear in a short time, as show Meanwhile, since the high energy accumulated near the wellbore is hard short time, the stresses of multiple points near the wellbore reach the stren chronously, resulting in multiple initiation points and even broken zones Figure 5a. Then, when the fracture encounters the gravel, it is more likely to the matrix since the energy spread is blocked, thus the deflection mainly permeability increases, more fluid can percolate into the reservoir, the pres tion rate of the fluid increases, and its effect for crack initiation starts to resulting in a single initiation point near the wellbore. The greater the per more the fracture is inclined to extend in the direction of the maximum pr forming a more regular fracture, namely, straighter bi-wing fracture, as sh 5d. The increase in permeability enhances the fluid percolation energy an rate of stress spread, thus the energy can be transmitted forward evenly. A hydraulic fracture has sufficient energy to penetrate the gravel to continue thus the penetration is increased. Additionally, high permeability reduces th centration near the wellbore, thereby increasing the initiation and penetra shown in Figure 6.

Effect of Permeability
The permeability can directly reflect the fluid percolation and pressure increase rates. Figure 5 shows the fracture propagation patterns under different permeabilities. When the permeability is small, it is difficult for the injection fluid to seep into the reservoir. Therefore, the fluid pressure in the wellbore rises rapidly, and the spread of the fluid pressure lags far behind, causing the surrounding stress to rapidly increase to the rock strength limit. Accordingly, the micro-cracks appear in a short time, as shown in Figure 6. Meanwhile, since the high energy accumulated near the wellbore is hard to release in a short time, the stresses of multiple points near the wellbore reach the strength limit synchronously, resulting in multiple initiation points and even broken zones, as shown in Figure 5a. Then, when the fracture encounters the gravel, it is more likely to extend inside the matrix since the energy spread is blocked, thus the deflection mainly occurs. As the permeability increases, more fluid can percolate into the reservoir, the pressure propagation rate of the fluid increases, and its effect for crack initiation starts to be dominant, resulting in a single initiation point near the wellbore. The greater the permeability, the more the fracture is inclined to extend in the direction of the maximum principal stress, forming a more regular fracture, namely, straighter bi-wing fracture, as shown in Figure 5d. The increase in permeability enhances the fluid percolation energy and reduces the rate of stress spread, thus the energy can be transmitted forward evenly. As a result, the hydraulic fracture has sufficient energy to penetrate the gravel to continue propagating, thus the penetration is increased. Additionally, high permeability reduces the energy concentration near the wellbore, thereby increasing the initiation and penetration times, as shown in Figure 6.

Effect of Gravel Strength
Gravel strength mainly affects the interaction between the fractures and the gravel, which affects the initiation and propagation of fractures. Figure 7a,b shows the propagation patterns when the gravel strength is 1.6 and 1.8 times the initial strength, respectively. The results show that as the gravel strength increases, the penetration decreases. When the fracture encounters the gravel, although the gravel strength is significantly larger than the matrix, the hydraulic fracture can always find the least resistance channel in the matrix to extend, thus it mainly appears in the form of deflection. When the ratio of gravel strength to initial gravel strength is less than 1.6, the hydraulic fracture can easily penetrate the gravel, and the influence of gravel on the fracture propagation direction is not clear; when the strength ratio is higher (≥1.8), the deflection occurs distinctly. After encountering the gravel, the fracture deflects multiple times along the boundary of the gravel, and at the same time forms a series of branch cracks, thus forming a complex fracture network. The greater the gravel strength, the more clear the fracture deflects, and the more complicated the final distribution morphology.

Effect of Gravel Strength
Gravel strength mainly affects the interaction between the fractures and the gravel, which affects the initiation and propagation of fractures. Figure 7a,b shows the propagation patterns when the gravel strength is 1.6 and 1.8 times the initial strength, respectively. The results show that as the gravel strength increases, the penetration decreases. When the fracture encounters the gravel, although the gravel strength is significantly larger than the matrix, the hydraulic fracture can always find the least resistance channel in the matrix to extend, thus it mainly appears in the form of deflection. When the ratio of gravel strength to initial gravel strength is less than 1.6, the hydraulic fracture can easily penetrate the gravel, and the influence of gravel on the fracture propagation direction is not clear; when the strength ratio is higher (≥1.8), the deflection occurs distinctly. After encountering the gravel, the fracture deflects multiple times along the boundary of the gravel, and at the same time forms a series of branch cracks, thus forming a complex fracture network. The greater the gravel strength, the more clear the fracture deflects, and the more complicated the final distribution morphology.

Effect of Injection Rate
The effect of injection rate on the initiation and propagation is similar to the effect of permeability. When the injection rate is high, the pressure in the wellbore increases rapidly due to the fact that the pressure spread lags behind, which causes energy concentration around the wellbore wall and multiple simultaneous initiation points, as shown in Figure 8; when injected at ultra-high speed, the wellbore wall even collapses and a broken zone is formed. At this time, the initiation and propagation are mainly affected by three factors: (1) Principal stress state; (2) random distribution of gravel and stress concentration caused by gravel irregularity; (3) energy concentration near the wellbore caused by highspeed injection, whose effect is mainly in the vicinity of the wellbore. In the early stage of hydraulic fracturing, the initiation and propagation of fractures are mainly affected by the latter two factors. As the fractures extend forward, the fluid pressure is difficult to spread to the far wellbore area in time, thus the initiation and propagation are mainly influenced by the principal stress. The fracture mainly extends along the direction of the maximum principal stress, thereby forming a regular bi-wing fracture. It is clear from Figure 8 that the local energy is difficult to release in time due to the high injection speed, thus the hydraulic fracture rapidly extends along the minimum resistance channel in the matrix.

Effect of Injection Rate
The effect of injection rate on the initiation and propagation is similar to the effect of permeability. When the injection rate is high, the pressure in the wellbore increases rapidly due to the fact that the pressure spread lags behind, which causes energy concentration around the wellbore wall and multiple simultaneous initiation points, as shown in Figure 8; when injected at ultra-high speed, the wellbore wall even collapses and a broken zone is formed. At this time, the initiation and propagation are mainly affected by three factors: (1) Principal stress state; (2) random distribution of gravel and stress concentration caused by gravel irregularity; (3) energy concentration near the wellbore caused by high-speed injection, whose effect is mainly in the vicinity of the wellbore. In the early stage of hydraulic fracturing, the initiation and propagation of fractures are mainly affected by the latter two factors. As the fractures extend forward, the fluid pressure is difficult to spread to the far wellbore area in time, thus the initiation and propagation are mainly influenced by the principal stress. The fracture mainly extends along the direction of the maximum principal stress, thereby forming a regular bi-wing fracture. It is clear from Figure 8 that the local energy is difficult to release in time due to the high injection speed, thus the hydraulic fracture rapidly extends along the minimum resistance channel in the matrix. As a result, with the increase in the injection speed, the relationship between the fracture and gravel is traversed from penetration to deflection.
As a result, with the increase in the injection speed, the relationship between the fracture and gravel is traversed from penetration to deflection.

Effect of Well Diameter
The size of the wellbore seriously affects the surrounding geometric heterogeneity, causing stress heterogeneity, which in turn affects the initiation and propagation of fractures. Figure 9 shows the fracture propagation morphologies at different wellbore radii. The results show that when the diameter of the wellbore is small, the internal energy of the wellbore is sufficient, and it is mainly concentrated in the direction of the maximum principal stress. As a result, the fracture starts to initiate on the left and right sides of the wellbore and gradually extends outward, eventually forming a straight bi-wing shape. When the diameter of the wellbore is large, the heterogeneity around the wellbore is enhanced, and the energy distribution around the wellbore is dispersed, in order that multiple points reach the strength limit synchronously, leading to multi-point initiation. Moreover, the fracture slightly deflects from the horizontal direction when extending forward, which finally turns into a bi-wing shape with branched cracks. Compared with the large-size wellbore, the energy concentration is increased when the wellbore is small. At a distance from the wellbore there is still sufficient energy at the fracture tip to penetrate the gravel, while for the larger wellbore, the energy of the fracture tip is insufficient, it mainly extends forward in the form of deflection.

Effect of Well Diameter
The size of the wellbore seriously affects the surrounding geometric heterogeneity, causing stress heterogeneity, which in turn affects the initiation and propagation of fractures. Figure 9 shows the fracture propagation morphologies at different wellbore radii. The results show that when the diameter of the wellbore is small, the internal energy of the wellbore is sufficient, and it is mainly concentrated in the direction of the maximum principal stress. As a result, the fracture starts to initiate on the left and right sides of the wellbore and gradually extends outward, eventually forming a straight bi-wing shape. When the diameter of the wellbore is large, the heterogeneity around the wellbore is enhanced, and the energy distribution around the wellbore is dispersed, in order that multiple points reach the strength limit synchronously, leading to multi-point initiation. Moreover, the fracture slightly deflects from the horizontal direction when extending forward, which finally turns into a bi-wing shape with branched cracks. Compared with the large-size wellbore, the energy concentration is increased when the wellbore is small. At a distance from the wellbore there is still sufficient energy at the fracture tip to penetrate the gravel, while for the larger wellbore, the energy of the fracture tip is insufficient, it mainly extends forward in the form of deflection.

Discussion
For glutenite reservoirs, strong heterogeneity is induced by embedded gravel, and this causes the initiation and propagation of hydraulic fractures to be more complicated. It was reported that a substantial "shielding" effect on near fracture can be produced by

Discussion
For glutenite reservoirs, strong heterogeneity is induced by embedded gravel, and this causes the initiation and propagation of hydraulic fractures to be more complicated. It was reported that a substantial "shielding" effect on near fracture can be produced by embedded gravel, and this would result in a reduction in the energy release rate [28,29]. Therefore, the stress near the interfaces between gravel and matrix quickly increases as hydraulic fractures approach the gravel, and this change subsequently causes local failures, which gradually attract the propagating hydraulic fractures. It was found that a hydraulic fracture is anticipated to be attracted by the encountering gravel due to interface failures as it deviates from the centerline and propagates the gravel. However, when a fracture propagates along the centerline of the gravel, it may deflect or the gravel is penetrated depending on local stress distribution. Due to the limitation of true tri-axial experiments, a physical experiment on hydraulic fracturing in a glutenite specimen is rare, and the only one was presented by Ma et al. [15] to the best of our knowledge. Here, the authors found three interaction patterns between propagating fractures and encountering gravel, which were defined as attraction, deflection, and penetration [15]. Moreover, we observed these interaction patterns in this work, suggesting excellent robustness of the proposed approach on simulating this complicated initiation and propagation of hydraulic fractures. Furthermore, the interaction patterns are significantly influenced by permeability, gravel strength, and injection rate, ect., thus complicated morphologies of hydraulic fractures are usually anticipated in glutenite reservoirs.
Finally, geological stress controls the global propagation of hydraulic fracture, while local stress state significantly influences its initiation and propagation, and their combination gives rise to complicated fracture morphologies. The pore pressure is affected by the permeability and injection rate, which changes the local stress state, and is significantly influenced by gravel strength and well diameter, thus different fracture morphologies were observed under different parameters. In addition, compared with geological stress which is almost constant, local stress distribution can be changed by regulating the working parameters during the hydraulic fracturing treatment, thus these findings are helpful in increasing the stimulated reservoir volume in fields. However, since this work is only a laboratory scale simulation, the upscale of the simulation is required for better understanding of the mechanisms behind the initiation and propagation of hydraulic fracture in glutenite reservoirs.

Conclusions
(1) Based on the coupled seepage-DEM method, the mechanism of hydraulic fracture initiation and propagation in glutenite was explored, the interaction between fractures and gravel was revealed, and the effects of geo-stress difference, permeability, gravel strength, injection rate, and well diameter were analyzed.
(2) When hydraulic fractures encounter gravel, attraction, deflection, penetration, and termination are generally observed. The interaction mechanism between hydraulic fractures and gravel is affected by local stress heterogeneity, fracture tip energy, and gravel strength.
(3) The initiation of hydraulic fracture near the wellbore wall is controlled by local stress heterogeneity, while the existence of gravel increases the stress heterogeneity, which leads to simultaneous multi-point initiation. Moreover, local stress concentration dominates the propagation of hydraulic fractures, and when they enter into the far wellbore area, the influence of the geological stress state is further highlighted.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.