Simulation of hydraulic fracture initiation and propagation for glutenite formations based on discrete element method

The existence of gravels in the glutenite formations leads to the complex geometries of hydraulic fracturing propagation and difficult construction in fracturing engineering. To study the hydraulic fracturing propagation law of glutenite formations, this paper establishes a fracture propagation model for the heterogeneous glutenite formations based on discrete element method, and analyzes the effects of gravel content, particle size, distribution, horizontal stress difference, fracturing fluid viscosity and flow rate on hydraulic fracturing propagation behavior. Results show that the complex geometries of hydraulic fractures in glutenite formations can lead to the generation of branched fractures and fracture bifurcation. Small-sized gravels have little effect on the fracture propagation shape which leads to a single main fracture with a flat fracture surface, on the contrary, large-sized gravels may induce hydraulic fractures to deflect along the gravel interface and form branched fractures with distorted fracture surfaces. Hydraulic fractures can propagate around gravels under the condition of high stress difference, high viscosity and medium flow rate. Gravels can prevent the propagation of hydraulic fractures under low stress difference, low viscosity and small flow rate. Hydraulic fracture bifurcation can occur when encountering gravels under high stress difference and large displacement. Properly increasing the high viscosity of fracturing fluids can effectively promote the main hydraulic fracture propagation and reduce the fracture tortuosity, thereby avoiding sand up.


Introduction
Due to lithology change, low porosity, poor permeability and strong heterogeneity, fracturing design and parameters optimization of glutenite formations face great challenge (Jia et al., 2017;Wang et al., 2017). For the gravels in glutenite formation, fracture distortion and multi-fracture propagation are prone to occur in the hydraulic fracturing process, resulting in sand up and affecting hydraulic fracturing effects (Meng et al., 2010;Li et al., 2013;Feng et al., 2016;Ma et al., 2017). Therefore, studies of the fracture propagation law is of great significance for the stimulation and reconstruction in glutenite formation.

OPEN ACCESS EDITED BY
A large number of numerical simulations (Feng et al., 2017;Zhang et al., 2017;Zhang et al., 2018;Lu et al., 2020;Wang et al., 2021;Wei et al., 2021) and experimental investigations (Fan et al., 2014;Ma et al., 2017;Moghadasi et al., 2019;Tan et al., 2019;Xie et al., 2022) have been carried out for the hydraulic fractures propagation of different lithologies. Tan et al. (2017), Tan et al. (2020) summarized fracture morphology types in mid-depth and deep shale formation. For middepth shales, fracture geometries were characterized by simple fracture, fishbone-like fracture, fishbone-like fracture with fissure opening, and multilateral fishbone-like fracture network (2017). For deep shales, fracture geometries were characterized by horizontal hydraulic fracture, transverse fracture, step-shaped hydraulic fracture with fissure opening, and multilateral step-shaped hydraulic fracture network (2020). Wan et al. (2019); Tan et al. (2021) came up with the concept of "lithology transition zone" for coal measure strata and studied the fracture propagation and penetration behavior; Zhao et al. (2007) compared hydraulic fractures propagation law in different lithologies such as basalt, conglomerate, marl. For glutenite formations; Luo et al. (2013) considered that particle size, content of gravel and the properties of the interface are the internal reason for irregular hydraulic fractures and pressure fluctuations based on the numerical simulation results; Meng et al. (2010) analyzed the influence of gravel particle size and horizontal stress difference on fracture propagation morphology and pressure curve using artificial samples; Li et al. (2013) analyzed the effects of the in situ stress and gravels on fracture propagation with RFPA-Flow, and summed up four typical fracture propagation modes at the gravel, including fracture arrest, deflection, penetration and adsorption. Based on the three-dimensional reconstructed heterogeneous glutenite model; Ju et al. (2016) studied the initiation and propagation of hydraulic fractures under different horizontal stress differences by using the discrete element method, and compared with the physical simulation results. Although this research reveals the fracture propagation law of glutenite formation to some extent, there are still many issues worthy of further discussion. On the one hand, the glutenite formation is different from the sandstone reservoir which has fast phase transformation, large changes in lithology and permeability, and poor sandstone connectivity and strong heterogeneity. The brittle fracture morphology of the glutenite formation is obviously different from that of the normal sandstone reservoir especially when glutenite formation contains high gravel content and large gravel particles, the fracture surface is extremely irregular and prone to multiple fractures (Ma et al., 2017). Also, multiple fractures will increase the filtration of the fracturing fluid, thus affecting the fracture and sand carrying capacity, resulting in difficult large-scale fracturing treatment. At present, many studies has been carried out on the physical characteristics of glutenite formations, and these studies focus mainly on reservoir characteristics and flow problems in the development process Wang et al., 2020), However, the research on fracture propagation law in glutenite formation is relatively less. On the other hand, it is impossible to obtain directly theoretical solution of the fracture propagation law under the restriction of the factors such as the great lithology changes and the strong heterogeneity for glutenite formations. The laboratory experiment is one of the methods that can be tried, and some intuitive observations can be obtained through the laboratory experiments (Ma et al., 2017). However, due to the uncertainty of the content, size, properties and distribution characteristics of the gravels in the glutenite formations, the laboratory results are generally very discrete, and it is difficult to conclude a repeatable conclusion by a limited number of tests.
Based on the physical and mechanical parameters, the discrete element method is also an effective way to study the hydraulic fracturing propagation law under complex conditions Tan et al., 2018;Huang et al., 2019;Huang et al., 2023). Therefore, this paper uses 2D Particle Flow Code (PFC2D, 2006) to analyze the hydraulic fracturing propagation law in glutenite formations, in order to guide the fracturing design and construction for glutenite formations.

Model calibration and schemes 2.1 Model calibration
In order to illustrate the reliability of PFC2D, simulation of uniaxial compression is carried out, and the simulation results are compared with the laboratory results. The size of the model is 25 mm × 50 mm. The top and bottom surfaces are set as rigid walls with constant displacement. The left and right sides are set to the constant pressure wall controlled by servo to simulate the loading of confining pressure. The confining pressure is set to 0.1 MPa. Since the pressure is very small compared with the fracture pressure of the samples, thus it can be ignored and can be approximated to uniaxial loading. The micro parameters in the model are shown in Table 1. The samples in the laboratory were collected from sandstone cores of a block in Changqing Oilfield. The rock mechanics parameters obtained by numerical simulation are close to those obtained by laboratory experiment, as shown in Table 2. The hydraulic fracture geometries in the laboratory and simulation are similar which are shear fractures, as shown in Figure 1. Therefore, it proves the reliability of the PFC 2D method to reconstruct the rock structure and simulate the hydraulic fracturing process.

Hydraulic coupling model
In this paper, PFC2D software was utilized to establish the fluid solid coupling model of hydraulic fracturing. In the model, each

Micromechanical parameters
Value Elastic modulus of the particle element, E mod (GPa)

1.2
Elastic modulus of the bonding element, P b _E mod (GPa)

17
Ratio of normal stiffness to tangential stiffness for particle element, k n /k 1.0 Ratio of normal stiffness to tangential stiffness for bonding element, k n /k s

1.2
Tensile strength for the bonding element, σ t (MPa) 11 Shear strength for the bonding element, τ (MPa) 21 Frontiers in Earth Science frontiersin.org particle contact is a flow channel to simulate the fluid flow, and that these channels contact up small "domains" that store some fluid pressure, as shown in Figure 2. The particles position can be changed by the change of contact force, and finally the change of channel pores is realized. The flow rate of fluid exchange can be expressed by Hagen Poiseuille equation (Hiroyuki et al., 2011): Where, ΔP is the pressure difference between the two pore basins, μ is the fluid viscosity, w is the aperture and L is the length of the fluid channel.
Each channel has an aperture associated with it. The aperture is calculated by:

FIGURE 2
The model of domain and flow channel in PFC2D (Tan et al., 2018).

FIGURE 3
The geometric model.

Frontiers in Earth Science frontiersin.org 03
Where, w 0 is the residual aperture for particles that are just touching, F is the normal stress at the contact and F 0 is the contact stress at which the aperture decreases to half of its initial aperture.
For each time increment, the change of fluid pressure ΔP is (Hiroyuki et al., 2011): Where, ΣQ is the total flow rate of fluids, Δt is the duration in one time step, K f is the compressive modulus of the fluid, V r is the "domain", that is the pore volume and ΔV r is the pore volume change.

Simulation scheme
A two-dimensional plane model was established by using PFC2D, as shown in Figure 3. The size of the model is 50 cm × 50 cm, and the diameter of the composed spherical particles is set to 1 mm, with a total of 6,547 in particles. The diameter of the middle circular hole is 3 cm, which is utilized to simulate the injection hole. In Figure 3, the maximum horizontal in situ stress is applied along the X-axis direction, and the minimum horizontal in situ stress is applied along the Y-direction. The simulation parameters are shown in Table 3.

Results and discussion
3.1 Effects of geological factors 1) In-situ stress Figure 4 shows the hydraulic fracture propagation patterns under different conditions of in situ stress differences for the two cases of whether there is gravel or not. The simulation results show that the presence of gravels increases the complexity of the hydraulic fracture propagation. When the horizontal stress difference is small, the fracture propagation is less restricted by the in situ stress, and radial fractures are easily formed along the wellbore direction. However, during the hydraulic fracture propagation, the branched fractures propagating in the direction of non-maximum in situ stress are prone to be arrested after encountering gravels, and the hydraulic fractures extending in the direction of maximum in situ stress become domination. When the horizontal stress difference is large, the hydraulic fractures tend to extend in the direction of the maximum horizontal stress, and there are few branched fractures. In this condition, the main hydraulic fracture will detour after encountering gravels during the propagation process.

2) Average gravel size
To study the influence of gravel particle size on hydraulic fracture propagation, five groups of calculation examples are designed. In each case, the principal stress difference is 2 MPa, and the gravel content is 30%. The simulation results are shown in Figure 5.
Simulation results show that, the distribution of gravels becomes more dispersed with the increase of the gravel size if the gravel content is same. The main fracture has enough spaces to propagate around the gravels, and the probability for encountering gravels decreases. However, under the control of the principal stresses, once the main fracture encounters the larger-sized gravel, more energy is required to bypass the gravel, and the hydraulic fracture will have a greater turning angle. When the particle size of the gravel is small, the number of the gravel increases, and the main fracture is more likely to encounter the gravels during propagation. However, the degree of stress concentration has been homogenized in the situation, and the hydraulic fracture can easily continue to propagate around the gravel. Hence, the inhibition and shielding effect of gravels on hydraulic fractures are relatively weak, and the fractures are relatively straight. Figure 6 shows the relationship between different gravel sizes and fracture initiation pressure. Results show that under the same gravel content, the fracture initiation pressure increases at initial stage and then decreases with the increase of gravel size. At a certain volume content, the size of the gravel has a probabilistic effect on the fracture propagation. When the gravel particle size increases, the sensitivity and uncertainty of the system stability increases, and the probability characteristics are obvious. The decreases in the gravel particle size contribute to the generation of micro-fractures clusters between the gravel, matrix and matrix, which reduces the rock strength and the fracture initiation pressure. (Khristianovic, 1955;Li et al., 2013).

3) Gravel content
In this section, to simulate the hydraulic fractures propagation under different gravel contents, six groups of calculation examples are designed with gravel contents of 0%, 10%, 20%, 30%, 40% and 50%. The average gravel size is 10 mm. The simulation results are shown in Figure 7.
Results show that the non-uniform distribution of gravels will lead to asymmetric hydraulic fracture propagation. With the increase of gravel content, micro-fractures at the end of the main fracture continue to initiate and propagate. The gravels distributed on the fracture propagation path will have an influence on the fracture propagation direction, which the gravels in other areas have little effect on the fracture propagation. With more gravel content on the propagation path, the fractures bypass and turn around the gravel to form S-shaped fractures.

1) Fracturing fluid viscosity
To study the influence of fracturing fluid viscosity on hydraulic fracture propagation, three groups of examples are designed in this section. The fracturing fluid viscosity is 10 mPa · s, 30 mPa · s and 100 mPa · s. The simulation results are shown in Figure 8.
Results show that the lower the viscosity of the fracturing fluid is, the easier the fracturing fluid is to loss along the fracture to the matrix, resulting in the generation of multi-branched fractures. Due to the joint distribution of hydraulic energy by the multi-branched hydraulic fractures, the length of the main hydraulic fracture becomes shorter. As the viscosity of fracturing fluid increases, the possibility of infiltration in all directions is reduced, as a result, the Frontiers in Earth Science frontiersin.org hydraulic fracture morphology is relatively straight and the length of the main fracture increases. Therefore, the length of the main hydraulic fracture can be increased by increasing the viscosity of fracturing fluids during constructions.

2) Flow rate of fracturing fluid
To study the effect of flow rate on the fracture propagation, two sets of calculation examples are designed in this section, and the Initiation pressures under different gravel sizes.

FIGURE 5
Fracture propagation geometry under different gravel sizes.

Frontiers in Earth Science
frontiersin.org 06 flow rates are 5 ml/ min and 50 ml/ min. Under each flow rates, two kinds of gravel contents and gravel sizes were set respectively to simulate the hydraulic fracture propagation process. The simulation results are shown in Figure 9.
The studies show that under small flow rate, the hydraulic fracture geometry is relatively simple, and the gravels can cause the arrest of the hydraulic fracture. Under large flow rate, after encountering gravels in the fracture propagation process, fluids in  Frontiers in Earth Science frontiersin.org 07 the fracture are likely to produce a "water hammer effect", which leads to the fracture bifurcation and multiple fractures. Contrastively, for homogeneous reservoirs, high flow rate also does not cause fracture bifurcation (Stanchits et al., 2014).

Conclusion
1) The gravels in the glutenite formation will have a significant impact on the initiation and propagation of hydraulic fractures. Compared with homogeneous sandstone formations, the glutenite formations own complex hydraulic fracture geometries, which could easily lead to the generation of branched fractures and bifurcation fractures, making nonplanar fracture type. 2) This study applied PFC to simulate the hydraulic fracturing propagation of glutenite formations, and there are three main behaviors for the hydraulic fracture as encountering gravels: bypass, stop and bifurcation. Under the high stress difference, high viscosity and moderate flow rate, hydraulic fractures can continue to propagate around the gravel. Under low in situ stress difference, lower viscosity and smaller flow rate, gravels can prevent the propagation of hydraulic fractures. Under the higher in situ stress difference and larger displacement, hydraulic fractures can bifurcate and continue to propagate.
3) The particle size and content of gravels have a certain influence on the fracture propagation. The larger gravel particle size makes it easier to cause the hydraulic fracture to stop or bypass at the gravel. Only the gravel distribution in the direction of hydraulic fracture propagation will cause the hydraulic fracture propagation changes, and the more gravels intensively distributes in the propagation direction, the more distorted the hydraulic fracture is, and which will result in a S-shaped fracture.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors. Frontiers in Earth Science frontiersin.org