# Numerical Investigation of Hydraulic Fracture Propagation Morphology in the Conglomerate Reservoir

Authors: Zehao Xu; Xiangjun Liu; Lixi Liang
Journal: Geofluids
 (2022)

---

## Abstract

The random distribution of gravels makes the conglomerate reservoir highly heterogeneous. A stress concentration occurs at the gravel-matrix interfaces owing to the embedded gravel and affects the local mechanical response significantly, making it difficult to control and predict hydraulic fracture (HF) propagation. The mechanism of HF propagation in conglomerate reservoirs remains unclear; thus, it is difficult to effectively design and treat hydraulic fracturing. Based on the global pore-pressure cohesive zone element (GPPCZ) model method, a two-dimensional (2D) fracture propagation model with flow-stress-damage (FSD) coupling was established to investigate HF nucleation, propagation, and coalescence in conglomerate reservoirs. This model was experimentally verified, and fractal theory was introduced to quantify the complexity of fracture morphology. The microscale interactions of the gravel, matrix, and interface have been taken into consideration during simulating HF propagation accurately in macroscale. The influence of the mechanical properties of gravel, matrix, matrix-gravel interface, and reservoir stress distribution state, on HF morphology (HF length, stimulated reservoir square, and HF complexity morphology), was investigated. Finally, the main factors affecting fracture propagation were analyzed. It was revealed that the difference between the mechanical properties of the gravel and the matrix in the conglomerate rock will affect the geometry of HF to varying degrees. The local behavior of fracture propagation is obviously dominated by the elastic modulus, tensile strength, and the strength for the matrix-gravel interface. However, the propagation of HF at the whole scale is mainly dominated by the horizontal stress state, including the minimum horizontal stress and horizontal stress difference. In addition, the difference in horizontal stress significantly affects the fracturing patterns (deflection, bifurcation, and penetration) when HF encounters gravel. In this study, a simulation method of HF propagation in conglomerate reservoirs is introduced, and the results provide theoretical support for the prediction of HF propagation morphology and plan design of hydraulic fracturing in conglomerate reservoirs.


---

## Body

## 1. Introduction

With the further exploitation of oil and gas reservoirs, the exploration and development of conglomerate reservoirs have received increasing attention [1]. Conglomerate reservoirs are generally deeply buried, resulting in poor physical properties of conglomerate and low porosity and permeability [1–3]. As the core analysis of conglomerate reservoirs show, the pore space in rock is intergranular, and the radius of the pore throat is far smaller than that of the conventional reservoir [4–7]. With all these reservoir characteristics, a complex fracture network is necessary to form to improve economic productivity using hydraulic fracturing in the development of conglomerate reservoirs [8]. The propagation trajectory of HF during the hydraulic fracturing process is mainly controlled by the stress state [9]. However, because of the presence of randomly embedded gravel in conglomerate reservoirs, the gravel may cause stress interference, deflect the propagation of the original fracture, and form more complex fractures [10]. Compared with conventional reservoirs, it is more difficult to predict the fracture morphology of fracturing in conglomerate reservoirs, and the effective implementation of fracturing design is full of challenges [11]. As the fracturing results of stimulation wells hinge upon the geometry of fractures, it is important to reveal the formation mechanism of HFs in conglomerate reservoirs.Hydraulic fracturing is a complex process involving seepage-stress-damage coupling, characterized by the initiation, expansion, penetration, and connection of new fractures. For sandstone reservoirs, HF generally propagates along the direction of maximum in situ stress and presents a double-wing shape [12]. However, for conglomerate reservoirs, high heterogeneity leads to an uneven strength and stress distribution, resulting in uncertainty in the formation and propagation of cracks [13]. Some techniques have been used in field monitoring and experiments to investigate the propagation law and terminal morphology of HF, such as X-ray CT scans, acoustic emission (AE), and digital image correlation (DIC) [14–21]. However, HF initiation and expansion cannot be clarified in detail owing to the limitation of accuracy. Ma et al. [22] conducted triaxial fracturing experiments and found that for isotropic horizontal stress, the interface properties have a great effect on the growth path of HF, and four interaction modes between fracture and gravel were observed: termination, penetration, deflection, and attraction. Large-scale true triaxial hydraulic fracturing equipment was used to study the propagation of HF under different conditions and combined with the fracturing treatment curves to quantitatively analyze the impact of gravel on fracture propagation [23]. Liu et al. [24] applied high-resolution microimaging technology in a hydraulic fracturing true triaxial physics experiment, obtaining a visual description of glutenite fractures. Intuitive data can be obtained through physical experiments. However, owing to the uncertainty in the content, size, properties, and distribution characteristics of gravel in the reservoir, the results of the physics experiment are generally highly discrete. The propagation law of fractures encountering randomly distributed gravel is still undefined, due to the limitation of the number of experimental specimens. In addition, fracture morphology under different conditions has not been quantitatively analyzed.In the past several years, based on basic physical and mechanical parameters, numerical simulation has become one of the most effective methods for revealing the initiation and propagation mechanisms of HF under complex conditions. The following methods are conducted for numerical simulation of complex HF propagation: (1) finite element method (FEM): in the process of simulating HF propagation, the FEM method cannot remesh the grid; thus, it can only predict planar fractures [25]; (2) conventional finite element method (CFEM): the accuracy of CFEM is limited by the element size in the model [15, 16, 18, 26, 27]; (3) extended finite element method (XFEM): the XFEM is used to simulate HF propagation with the characteristics of high precision and easy operation. However, when there were abundant of heterogeneous media (natural fractures, weak surface), the calculation may not converge. In addition, the simultaneous propagation of fracture bifurcation and multiple fractures cannot be simulated with this method, and the interaction between the HFs cannot be described [28]. The heterogeneity of the rock has the potential to cause problems in solving as well. (4) The boundary element method (BEM): even though the BEM has advantages in dealing with the statistical propagation of HFs, it has limitations in handling the load and flux generated in the fracture plane [29]. For more complex media conditions, such as increasing the number of layers for heterogeneous media, the solution process of this method will become complicated, which severely limits the scale of the model. (5) Nonelement method (NEM): the method is proposed to simulate HF propagation, which can better deal with HF statistical propagation. However, the theoretical basis and computational efficiency of this method are low and need to be further improved [30]. (6) Analytical model method: the judgment of the HF propagation direction is based on analytical models and equations, and low simulation accuracy makes it difficult to meet the needs [31]. (7) The cohesive zone element method (CZM): this method is implemented to model HF propagation effectively, which does not have the above limitations. The cohesive element was first introduced by Barenblatt [32, 33] and used in brittle materials to simulate crack propagation. Vyacheslav [34] presented the CZM for fracturing soft rock, which led to a more accurate fitting of the pressure log. T.K. Guo et al. [5] and J. Guo et al. [4] used CZM to study the influence of geomechanical parameters and natural fractures on HF in layered oil reservoirs. Considering the respective advantages of XFEM and CZM, Wang [35] combined the two methods to study the fracture reorientation caused by the difference in perforation and analyzed the interference and merging of multiple fractures in the multistage fracturing process. Manchanda et al. [36] proposed 3D-CZM to study the influence of heterogeneity of rock characteristics on multiple HF propagation. Based on CZM, C. Guo et al. [15], T. Guo et al. [16], and J. Guo et al. [18] simulated the interaction of HF with a single natural fracture, and Arash et al. [37] extended their models and introduced natural fracture networks into the simulation. Li et al. [38, 39] presented a new pore pressure cohesive zone element (PPCZ) for modeling the propagation of hydraulically induced fracture. Wang [40] developed a new model using PPCZ for poroelastic and porous plastic formations to simulate the propagation of HF. Baykin and Golovin [41] proposed a fully coupled poroelastic cohesive model to study the influence of penetration on HF propagation. The PPCZ method has been widely used in hydraulic fracturing simulations, but little attention has been paid to the HF propagation mechanism in conglomerate reservoirs. The random distribution of gravels and the heterogeneity of conglomerates significantly affect the formation and propagation of fractures, leading to a complex fracture morphology [42, 43]. The influence of rock mechanical properties, gravel content, size, roundness, reservoir in situ stress state, and hydraulic fracturing treatment parameters on HF propagation had been studied [15, 16, 18, 38, 39, 44, 45]. However, most existing models do not account for the influence of rock heterogeneity on the HF pattern. The weak matrix-gravel interface or the superposition of multilevel gravel are not considered. Such assumptions are inconsistent with the actual conglomerate reservoir characteristics.In this paper, the PPCZ elements are embedded between any two adjacent quadrilateral finite elements by an ABAQUS plugin written in the Fortran programming language, to obtain a global pore pressure cohesive zone element (GPPCZ). Through the GPPCZ method, the finite element method and the discrete element method could be coupled to simulate the initiation, propagation, and intersection of fluid-driven fractures in rock. The solid element is used to represent the continuous properties of the material, and the PPCZ element is used to represent the discontinuous properties of the material. The interfaces between multiple solid elements are connected by PPCZ elements, which can represent weak surfaces in the rock. The fracture of the solid elements is achieved by the fracture of the PPCZ element, providing a potential path for the failure.In this study, in order to reflect the strong heterogeneity of conglomerate, the conglomerate rock is regarded as a three-phase composite material composed of gravel, matrix, and gravel-matrix interface on a microscopic scale. A 2D FSD fully coupled finite element model is established based on GPPCZ, to simulate HF propagation in conglomerate reservoirs under multifield coupling. Compared with the published physical test results and numerical simulation results, the HF propagation morphology is analyzed at the microscale and macroscale to verify the effectiveness of the method. Fractal theory is introduced to quantitatively describe the geometry of fractures in this paper. Finally, taking the conglomerate reservoir in the Mahu Sag, Junggar Basin, as an example, considering the structural characteristics and local material properties of the conglomerate, the microscale interactions between particles-gravel and interparticles are studied, to predict fracture morphology in macroscale accurately. A series of numerical simulations under multifield coupling conditions is carried out. The influence of the gravel, matrix, and interface mechanics properties; horizontal stress difference; and minimal horizontal stress on HF propagation is studied, aiming at capturing the profound understanding of HF propagation morphology prediction in conglomerate reservoirs.


## 2. Modeling Approach

### 2.1. Formulations and Numerical Implementation

As shown in Figure1(a), the hydraulic fracturing model in the formation is under the action of the fracturing fluid. The y-axis is a schematic diagram of a vertical wellbore, the left side is a schematic diagram of the reservoir fractures, and the right side is a fracturing fracture simulation element. The fracture element is filled with pressurized fluid (blue). Under the action of fluid pressure, different parts of the fracture element have different responses: the blue elements are the areas that have been ruptured, and the fracture tip A on the right is the rupture process area that is about to rupture. In the GPPCZ method, the main theory used for calculation includes two parts: fluid flow equations in pores and fractures and traction separation criterion.Figure 1
Fluid flow in pores and fractures during hydraulic fracturing.The fractured element is filled with pressurized fluid, which has a normal flow perpendicular to the top and bottom planes of the element, and a tangential flow parallel to the plane of the element. The normal flow indicates that the fracturing fluid has leaked into the formation, and the tangential flow forces the fracture to propagate. According to Newtonian rheological theory, the fluid in the fracture element is considered incompressible. The displacement at any time isq, and the tangential flow at the fracture plane satisfies the Newtonian fluid pressure conduction equation [46]:
(1)q=−kt∇pf,where q is the tangential flow, which is equal to the average velocity of the tangential flow multiplied by the crack width; kt is the flow coefficient; and ∇pf is the fluid pressure. According to the Reynolds equation, kt can be expressed as
(2)kt=w312μ,where w is the crack width and μ is the fracturing fluid viscosity.For tangential flow, a power-law model is used to characterize the fluid shear stress [47], and the constitutive relationship is
(3)τ=μγ̇n′,where τ is the fluid shear stress, γ̇ is the tangential strain rate, and n′ is the power-law coefficient.The normal fluid loss in the fracture is expressed as the normal flow of the fluid in the element [48], that is, the fluid loss along the top and bottom planes of the cohesive element shown in Figure 1(b). The calculation equation is
(4)qt=ctPf−Pt,qb=cbPf−Pb,where qt and qb are the volume flow rates of fracturing fluid on the top and bottom planes of the element, respectively; ct and cb are the leak-off coefficients of the top and bottom planes of the element, respectively; Pf is the pressure of the fluid in the element; and Pt and Pb are the pore pressures in the bulk element immediately adjacent to the cohesive element.Currently, in the study of HF propagation using the displacement discontinuity method (DDM) and the XFEM, to strengthen the calculation convergence, most of the reservoirs are considered as elastic entities regardless of porosity or permeability [49]. Such assumptions are inconsistent with the actual reservoir characteristics. In the numerical model used in this paper, the conglomerate rock is regarded as a porous continuous medium. The fracturing fluid exchanges with the fluid in the pores along the solid phase framework during the fracturing process, which is characterized by the fluid-solid coupling equation. Based on the Biot pore elastic theory, the equilibrium equation of the rock skeleton stress is expressed in the form of a virtual work principle [50]:
(5)∫Vσ¯+PpI:εδdV=∫StvδdS+∫VfvδdV,where σ¯ is the effective stress matrix; I is the second-order element tensor; Pp is the pore pressure; εδ is the virtual strain matrix; and t, f, and vδ are the plane force, volume force, and virtual velocity vectors, respectively.The seepage process of fluid in rock pores obeys Darcy’s law, and the differential form of the fluid continuity equation is [50]
(6)ddt∫VρwnwdV+∫Sρwnwn⋅−1nwgρwk⋅∇Pw−ρwgdS=0,where ρw is the density of the fluid in the pores, nw is the porosity, Pw is the fluid pressure, k is the permeability matrix, g is the acceleration due to gravity, and n is the element vector perpendicular to plane S.


### 2.2. Traction Separation Criterion

The linear elastic traction-separation criterion is used to judge the initial damage and final failure of HFs [51]. As shown in Figure 2, quantitative characterization is achieved through the degree of damage in the cohesive element, and the destruction process of the element consists of three parts: initial stage of loading: the normal displacement of the top and bottom planes of the cohesive element is smaller than the initial damage. The normal stress increases linearly with the increase of the displacement and reaches to the tensile strength. When the normal displacement reaches to the limit value of the element, the cohesive element is in the damage phase. The normal stress that the cohesive element can withstand decreases with an increase in displacement. When the normal displacement reaches to the limit failure displacement of the cohesive element, the element is destroyed, and artificial cracks begin to emerge.Figure 2
Bilinear traction-separation law.In complex fracture propagation, there are multiple extension modes after the HF intersects with the potential propagation paths (cohesive elements). The second nominal stress criterion is used to judge whether microcracks appeared, and the influence of the normal stress and shear stress of the cohesive element is analyzed. Assuming that the sum of the squares of the stress in the three directions of the element and its critical stress ratio reaches 1, the element starts to be damaged and artificial cracks begin to occur. The criterion can be expressed as [50]
(7)tntn02+tsts02+tttt02=1,where tn is the normal stress components; ts and tt are the first and second shear stress components, respectively; tn0, ts0, and tt0 are the normal, first, and second shear strengths of the complete cohesive element, respectively; and the symbol “” is a Macaulay bracket, which means that the cohesive element will not be damaged under pure extrusion deformation or stress state.To characterize the degree of damage of the cohesive element, a damage factorD is introduced into the model. The values of D are 0 and 1, as shown in Figure 2, indicating that the material is undamaged and completely damaged, respectively. According to the linear displacement expansion criterion, D is expressed as [52, 53]:
(8)D=δmfδmmax−δm0δmmaxδmf−δm0,(9)δm=δn2+δs2+δt2,where δmmax is the limit displacement value of the element, δmf is the displacement value when the element starts to crack, δm0 is the displacement value when the element starts to be damaged, δm represents the effective displacement, δn is the normal displacement component, and δs and δt are the first and second shear displacement components, respectively.δmf can be defined as
(10)δmf=2GcTmax,where Tmax is the tensile strength or shear strength of the material and Gc is the mixed-mode fracture energy.According to the Benzeggagh-Kenane fracture criterion, the damage degree of the cohesive element is evaluated, and the damage evolution of the element during the fracture propagation process is determined [15, 16, 18]:
(11)Gnc+Gsc−GncGsGTη=Gc,where Gn, Gs, and Gt are the work done by the traction force in the normal, first, and second shear directions, respectively.
(12)Gsh=Gsc+Gtc,(13)GT=Gsh+Gnc,where Gnc, Gsc, and Gtc are the fracture energy in the normal, first, and second shear directions, respectively.Using the damage factor, the cohesive element damage evolution model based on the traction-separation criterion is expressed as [15, 16, 18]
(14)Ts=1−DTs¯,(15)Tt=1−DTt¯,(16)Tt=1−DTn¯Tn¯,Tn¯≥0,,Tn¯<0,where T¯n, T¯s, and T¯t are the stresses calculated by the three-direction cohesive element during the elastic deformation.


### 2.3. Global Pore-Pressure Cohesive Zone Element Model

In this study, conglomerate rock is regarded as a three-phase composite material on a microscopic scale, which is composed of gravel, matrix, and the gravel-matrix interface. According to the Fuller grading curve and the Walraven plane transformation equation, a 2D gravel random distribution is realized by programming. The Fuller-Walraven plane transformation equation is expressed as follows [54]:
(17)A=D0Dmax,(18)PD<D0=Pk1.065A0.5−0.053A4−0.012A6−0.0045A8−0.0025A10,where P is the percentage of the gravel passing through the sieve with a diameter of D; D, D0, and Dmax are the mesh diameter and the minimum and maximum gravel particle sizes, respectively; and Pk represents the percentage of gravel volume to total volume.Programmatically insert PPCZ elements into any two adjacent quadrilateral finite elements to obtain the GPPCZ model, as shown in Figure3.Figure 3
The GPPCZ model.


## 2.1. Formulations and Numerical Implementation

As shown in Figure1(a), the hydraulic fracturing model in the formation is under the action of the fracturing fluid. The y-axis is a schematic diagram of a vertical wellbore, the left side is a schematic diagram of the reservoir fractures, and the right side is a fracturing fracture simulation element. The fracture element is filled with pressurized fluid (blue). Under the action of fluid pressure, different parts of the fracture element have different responses: the blue elements are the areas that have been ruptured, and the fracture tip A on the right is the rupture process area that is about to rupture. In the GPPCZ method, the main theory used for calculation includes two parts: fluid flow equations in pores and fractures and traction separation criterion.Figure 1
Fluid flow in pores and fractures during hydraulic fracturing.The fractured element is filled with pressurized fluid, which has a normal flow perpendicular to the top and bottom planes of the element, and a tangential flow parallel to the plane of the element. The normal flow indicates that the fracturing fluid has leaked into the formation, and the tangential flow forces the fracture to propagate. According to Newtonian rheological theory, the fluid in the fracture element is considered incompressible. The displacement at any time isq, and the tangential flow at the fracture plane satisfies the Newtonian fluid pressure conduction equation [46]:
(1)q=−kt∇pf,where q is the tangential flow, which is equal to the average velocity of the tangential flow multiplied by the crack width; kt is the flow coefficient; and ∇pf is the fluid pressure. According to the Reynolds equation, kt can be expressed as
(2)kt=w312μ,where w is the crack width and μ is the fracturing fluid viscosity.For tangential flow, a power-law model is used to characterize the fluid shear stress [47], and the constitutive relationship is
(3)τ=μγ̇n′,where τ is the fluid shear stress, γ̇ is the tangential strain rate, and n′ is the power-law coefficient.The normal fluid loss in the fracture is expressed as the normal flow of the fluid in the element [48], that is, the fluid loss along the top and bottom planes of the cohesive element shown in Figure 1(b). The calculation equation is
(4)qt=ctPf−Pt,qb=cbPf−Pb,where qt and qb are the volume flow rates of fracturing fluid on the top and bottom planes of the element, respectively; ct and cb are the leak-off coefficients of the top and bottom planes of the element, respectively; Pf is the pressure of the fluid in the element; and Pt and Pb are the pore pressures in the bulk element immediately adjacent to the cohesive element.Currently, in the study of HF propagation using the displacement discontinuity method (DDM) and the XFEM, to strengthen the calculation convergence, most of the reservoirs are considered as elastic entities regardless of porosity or permeability [49]. Such assumptions are inconsistent with the actual reservoir characteristics. In the numerical model used in this paper, the conglomerate rock is regarded as a porous continuous medium. The fracturing fluid exchanges with the fluid in the pores along the solid phase framework during the fracturing process, which is characterized by the fluid-solid coupling equation. Based on the Biot pore elastic theory, the equilibrium equation of the rock skeleton stress is expressed in the form of a virtual work principle [50]:
(5)∫Vσ¯+PpI:εδdV=∫StvδdS+∫VfvδdV,where σ¯ is the effective stress matrix; I is the second-order element tensor; Pp is the pore pressure; εδ is the virtual strain matrix; and t, f, and vδ are the plane force, volume force, and virtual velocity vectors, respectively.The seepage process of fluid in rock pores obeys Darcy’s law, and the differential form of the fluid continuity equation is [50]
(6)ddt∫VρwnwdV+∫Sρwnwn⋅−1nwgρwk⋅∇Pw−ρwgdS=0,where ρw is the density of the fluid in the pores, nw is the porosity, Pw is the fluid pressure, k is the permeability matrix, g is the acceleration due to gravity, and n is the element vector perpendicular to plane S.


## 2.2. Traction Separation Criterion

The linear elastic traction-separation criterion is used to judge the initial damage and final failure of HFs [51]. As shown in Figure 2, quantitative characterization is achieved through the degree of damage in the cohesive element, and the destruction process of the element consists of three parts: initial stage of loading: the normal displacement of the top and bottom planes of the cohesive element is smaller than the initial damage. The normal stress increases linearly with the increase of the displacement and reaches to the tensile strength. When the normal displacement reaches to the limit value of the element, the cohesive element is in the damage phase. The normal stress that the cohesive element can withstand decreases with an increase in displacement. When the normal displacement reaches to the limit failure displacement of the cohesive element, the element is destroyed, and artificial cracks begin to emerge.Figure 2
Bilinear traction-separation law.In complex fracture propagation, there are multiple extension modes after the HF intersects with the potential propagation paths (cohesive elements). The second nominal stress criterion is used to judge whether microcracks appeared, and the influence of the normal stress and shear stress of the cohesive element is analyzed. Assuming that the sum of the squares of the stress in the three directions of the element and its critical stress ratio reaches 1, the element starts to be damaged and artificial cracks begin to occur. The criterion can be expressed as [50]
(7)tntn02+tsts02+tttt02=1,where tn is the normal stress components; ts and tt are the first and second shear stress components, respectively; tn0, ts0, and tt0 are the normal, first, and second shear strengths of the complete cohesive element, respectively; and the symbol “” is a Macaulay bracket, which means that the cohesive element will not be damaged under pure extrusion deformation or stress state.To characterize the degree of damage of the cohesive element, a damage factorD is introduced into the model. The values of D are 0 and 1, as shown in Figure 2, indicating that the material is undamaged and completely damaged, respectively. According to the linear displacement expansion criterion, D is expressed as [52, 53]:
(8)D=δmfδmmax−δm0δmmaxδmf−δm0,(9)δm=δn2+δs2+δt2,where δmmax is the limit displacement value of the element, δmf is the displacement value when the element starts to crack, δm0 is the displacement value when the element starts to be damaged, δm represents the effective displacement, δn is the normal displacement component, and δs and δt are the first and second shear displacement components, respectively.δmf can be defined as
(10)δmf=2GcTmax,where Tmax is the tensile strength or shear strength of the material and Gc is the mixed-mode fracture energy.According to the Benzeggagh-Kenane fracture criterion, the damage degree of the cohesive element is evaluated, and the damage evolution of the element during the fracture propagation process is determined [15, 16, 18]:
(11)Gnc+Gsc−GncGsGTη=Gc,where Gn, Gs, and Gt are the work done by the traction force in the normal, first, and second shear directions, respectively.
(12)Gsh=Gsc+Gtc,(13)GT=Gsh+Gnc,where Gnc, Gsc, and Gtc are the fracture energy in the normal, first, and second shear directions, respectively.Using the damage factor, the cohesive element damage evolution model based on the traction-separation criterion is expressed as [15, 16, 18]
(14)Ts=1−DTs¯,(15)Tt=1−DTt¯,(16)Tt=1−DTn¯Tn¯,Tn¯≥0,,Tn¯<0,where T¯n, T¯s, and T¯t are the stresses calculated by the three-direction cohesive element during the elastic deformation.


## 2.3. Global Pore-Pressure Cohesive Zone Element Model

In this study, conglomerate rock is regarded as a three-phase composite material on a microscopic scale, which is composed of gravel, matrix, and the gravel-matrix interface. According to the Fuller grading curve and the Walraven plane transformation equation, a 2D gravel random distribution is realized by programming. The Fuller-Walraven plane transformation equation is expressed as follows [54]:
(17)A=D0Dmax,(18)PD<D0=Pk1.065A0.5−0.053A4−0.012A6−0.0045A8−0.0025A10,where P is the percentage of the gravel passing through the sieve with a diameter of D; D, D0, and Dmax are the mesh diameter and the minimum and maximum gravel particle sizes, respectively; and Pk represents the percentage of gravel volume to total volume.Programmatically insert PPCZ elements into any two adjacent quadrilateral finite elements to obtain the GPPCZ model, as shown in Figure3.Figure 3
The GPPCZ model.


## 3. Model Verification and Comparison

It is essential to test the mechanical and hydraulic behaviors of the GPPCZ method before its use. The reliability of this model had been verified by a number of researchers, and the simulation results based on the GPPCZ method were compared against the analytic solutions of the KGD problem [38, 39, 50, 55]. To further verify the reliability of the model, the simulation results were compared against the triaxial fracturing experiment results presented by Ma et al. [22]. A plane square model was established, with the same experimental parameters and similar gravel sizes and distribution morphologies. The accuracy of the model was validated by determining whether the HF propagation morphology and propagation law of the fracture tip encountering the gravel in the numerical simulation were the same as those in the physical simulation.The validated model with a size of0.3m×0.3m contains 48,269 quadrilateral finite elements and 98,188 PPCZ elements; the particle radius spread was 5mm<R<20mm, following a uniform distribution, and 223 gravels were randomly embedded in the specimen. The properties of the gravel, matrix, and matrix-gravel interface are provided in Table 1 and set to be the same as those presented by Ma et al. [22].Table 1
Experimental properties in the physical simulation and numerical validated model.
Model properties and unitValuePermeability of matrix (10−3 μm2)0.57Elastic modulus of matrix (GPa)36.7Poisson’s ratio of matrix0.248Tensile strength of matrix (MPa)8.1Shear strength of matrix (MPa)81Porosity of matrix0.1Critical tensile energy of matrix (J/m2)40Critical shear energy of matrix (J/m2)1800Tensile strength of interface (MPa)7.4Critical tensile energy of interface (J/m2)20Critical shear energy of interface (J/m2)900Maximum horizontal stress (MPa)15Minimum horizontal stress (MPa)5Permeability of gravel (10−3 μm2)0.12Elastic modulus of gravel (GPa)50.1Poisson’s ratio of gravel0.247Tensile strength of gravel (MPa)12.4Shear strength of gravel (MPa)124Porosity of gravel0.05Critical tensile energy of gravel (J/m2)100Critical shear energy of gravel (J/m2)4500Acceleration due to gravity (m/s2)9.8Fluid viscosity (mPa·s)120Leak-off coefficient (m3/s/Pa)1E-14Initial pore pressure (MPa)0Injection rate (ml/min)20The comparison indicates that the HF propagation from the numerical simulation results is consistent with that from the physical simulation, not only in terms of fracture global propagation morphology but also in concrete details. As per the results of the numerical simulation shown in Figure4, the global fracture geometry of this specimen is mainly governed by the horizontal stress state and the embedded gravel. The HF deflects locally after encountering the gravel, which makes the HF growth path tortuous and complex. However, with the influence of the horizontal stress difference Δσ=10MPa, the HF propagates along the direction perpendicular to the minimum horizontal stress, forming a relatively straight fracture on a macroscopic scale.Figure 4
Propagation morphology of HF in the specimen.These propagation morphology characteristics of the HF are consistent with those proposed by Ma (Figures5(a) and 5(b)) and consistent with the simulation results of Li et al. [38, 39], Rui et al. [44], Zhang et al. [45], and others. Previous studies have shown that the embedded gravel will have a substantial “shielding” effect on the HF propagation, resulting in slower energy release, which is a manifestation of the stress intensity factor and toughness locally [56]. When HF propagates to the vicinity of the gravel, the stress close by the matrix-gravel interface increases rapidly, causing failures at the interface; thus, HF is easily attracted to the interface (Figure 4), which is consistent with the results from a physical experiment, as shown in Figure 5(d). Interestingly, it can be noted from Figure 4 that if the HF initially propagates away from the gravel at a distance, the HF is usually attracted by the gravel owing to the interface failure formed earlier. However, when HF propagates and hits the gravel, it may deflect along the interface or damage and penetrate the gravel. Three types of HF gravel particle intersections (termination, penetration, and deflection) occur in this simulation, which were also observed in physical experiment results presented by Ma et al. [22]. In addition, in Figures 5(c) and 5(f), a process zone with multiple narrow fractures in the matrix is observed, and these fractures compete for propagation space. Figure 5(e) shows a large opening fracture in the gravel; dislodged particles inside the fracture can be observed. Figure 5(f) shows that the HF forms a process zone and is divided into multiple branches. These phenomena are reflected in the numerical simulation. Ultimately, the GPPCZ method can not only accurately simulate the global HF propagation in the conglomerate reservoir but also capture the microscopic interactions between propagating HF and gravel. It is foreseeable that the simulation results obtained in this work are highly acceptable.Figure 5
Physical experimental propagation morphology of HF presented by Ma et al.
(a)(b)(c)(d)(e)(f)


## 4. Hydraulic Fracturing Propagation Simulation in Conglomerate Reservoirs

In this study, the conglomerate reservoirs of well A in the Mahu Sag, Junggar Basin, NW China, were used as the modeling objects. Conglomerate reservoirs are developed prominently in Mahu Sag, and breakthroughs have been made in the exploration of conglomerate reservoirs in recent years [57, 58]. According to the data of well A, the average gravel volume content in the conglomerate reservoirs of this well is 63.5%. The diameter of gravel varies, and the maximum gravel diameter is approximately 10 cm, which is generally 0.5~2 cm. Three types of gravel sizes were designed in this specimen: φ=9~10cm (Figure 6(a)), φ=4~5cm (Figure 6(b)), and φ=1~2cm (Figure 6(c)). As shown in Figure 6(d), the specimen was of size 2m×2m, in which the diameter of the wellhole was 0.166 m, and gravels of various diameters were superimposed together. The specimen had a total of 49,985 quadrilateral finite elements and 99,548 PPCZ elements with a minimum mesh size of 0.01 m. The percentage of the total gravel area was 62%, and the percentage of the matrix area was 38%. The input parameters for the simulation are listed in Table 2. The HF propagation morphology after 1 s of injection and 50 times magnification is shown in Figure 6(d).Figure 6
HF propagation after 1 s of fracturing and 50 times magnification.Table 2
Input properties in the simulation.
Model properties and unitValuePermeability of matrix (10−3 μm2)10Elastic modulus of matrix (GPa)18Poisson’s ratio of matrix0.28Tensile strength of matrix (MPa)3Shear strength of matrix (MPa)30Porosity of matrix0.1Critical tensile energy of matrix (J/m2)100Critical shear energy of matrix (J/m2)1000Tensile strength of interface (MPa)3Critical tensile energy of interface (J/m2)80Critical shear energy of interface (J/m2)800Maximum horizontal stress (MPa)93Minimum horizontal stress (MPa)90Initial pore pressure (MPa)65Permeability of gravel (10−3 μm2)1Elastic modulus of gravel (GPa)22Poisson’s ratio of gravel0.26Tensile strength of gravel (MPa)6Shear strength of gravel (MPa)60Porosity of gravel0.05Critical tensile energy of gravel (J/m2)150Critical shear energy of gravel (J/m2)1500Acceleration due to gravity (m/s2)9.8Fluid viscosity (mPa·s)120Leak-off coefficient (m3/s/Pa)1E-14Percentage of the gravel62%Percentage of the matrix38%Injection rate (m3/s)0.001As the results show, the fractures on both sides of the wellbore developed significant fracture bifurcations. The HF propagation morphology is mainly governed by the in situ stress state on the whole, and the propagation trajectories are significantly affected by the embedded gravel, making the HF growth path highly tortuous and complex locally. Fractures are easier to branch from the matrix, which increases the complexity of the fracture network.To quantitatively describe the fracture geometry, the HF is extracted as shown in Figure7. The maximum length and width of the left fracture (Figure 7(a)), L1 and H1, are 0.750 m and 0.240 m, respectively; those of the right fracture, L2 and H2, are 0.650 m and 0.300 m, respectively. The total fracture length and the stimulated reservoir square, L and SRS, are 1.400 m and 0.375 m2, respectively.
(19)L=L1+L2,(20)SRS=L1×H1+L2×H2.Figure 7
Extraction of fracturing fractures.
(a)(b)The hydraulic fracturing results depend on the morphology of the HF. The more complex the fractures, the better the stimulation effect of hydraulic fracturing [59]. The fracture distribution after fracturing has a fractal structure with a statistical sense of self-similarity. The fractal dimension can be used for the quantitative evaluation of fracture network complexity [60]. In this study, the fractal dimension was determined based on the grid covering method of the box dimension method.Firstly, a fractal set is defined, and the equation is [60]
(21)Nλ=Aλ−D,where N is the number of grids, λ is the length of the grid, A is the ratio constant, and D is the fractal dimension of the fracture network.Take the logarithm of both sides of the equation:(22)lnNλ=lnA−Dlnλ.The method has three steps: (1) Cover the entire HF morphology with a grid with side lengthλ and count the number of square grids containing fractures, as shown in Figure 7(b). (2) Gradually change the side length of the square grid and count the corresponding number of grids. (3) Take the logarithm of the side length λ and the corresponding number of grids; afterward, use lnλ as the abscissa and lnNλ as the ordinate. Use the least square method to perform regression analysis on the statistical data. If the HF morphology has significant fractal characteristics, it will satisfy the linear relationship shown in Equation (22), and the slope will be the fractal dimension D.The fracture fractal statistics according to Equation (22) are plotted in Figure 8. The slope of the regression line is −0.8532 and R2 is 0.9966, indicating that it has good fractal characteristics. It can be obtained that the fracture fractal dimension of the fractured reservoir is 0.8532. The larger the fractal dimension, the more complex the fracture morphology.Figure 8
Calculation of fracture fractal dimension.


### 4.1. Effect of Elastic Modulus

Elastic modulus is expressed as the ability of a rock to resist deformation under stress. Rocks with a high elastic modulus have an ability to resist deformation, which are more rigid and brittle. Thus, there is a more possibility of brittle failure for these rocks. In contrast, rocks with a low elastic modulus have greater flexibility. Previous studies have shown that, compared with the elastic modulus of gravel or matrix, “the ratio of the elastic modulus of the matrix to gravel (Em/Eg)” has a larger influence on the propagation morphology of HF [44]. In this study, models with Em/Eg were set to 0.5, 0.7, 0.8, 1.0, 1.2, 1.5, and 2.0, and the gravel elastic modulus was 22 GPa, as shown in Figure 9.Figure 9
HF propagation morphology under different ratios of the elastic modulus of the matrix to gravel.The elastic modulus of gravel or matrix has a great influence on the propagation trajectory of HF in the conglomerate reservoir. As shown in Figure9, the fracture in the specimen with Em/Eg=0.5 extends for a short distance. When Em/Eg=0.5, i.e., the lower the elastic modulus of the matrix, the higher the probability of ductile failure. HF tends to turn, kink, and branch in the vicinity of gravel with lower elastic modulus, causing local bending of the main fracture, resulting in short fracture length. With increased elastic modulus of the matrix, the heterogeneity of the elastic modulus for the rock material is weakened, and the fracture morphology should become simple. However, a positive relationship between D and Em/Eg is shown in Figure 10(c), which means the increase of the matrix elastic modulus results in an increment of the complexity of the fracture morphology. The reasons behind this situation are that the increased matrix elastic modulus results in an increased matrix stiffness and, thus, higher possibility of brittle failure. The fracture morphology becomes straight and extends further (Figure 10(a)), which is due to the decrease of local bending of the fracture on a microscale. As the elastic modulus of the matrix increases, the SRS tends to increase (Figure 10(b)). In addition to the elastic modulus, the fracture propagation is also affected by other differences in strength (tensile strength, etc.), resulting in a macroscopically distinct branching of HF. As the elastic modulus increases, the influence on HF morphology is weakening, while the influence of other strength differences is strengthening. The rate of change of D gradually decreases with the increases of Em/Eg which illustrates this point.Figure 10
Geometric parameters of HF under different ratios of the elastic modulus of the matrix to gravel.
(a)(b)(c)


### 4.2. Effect of Tensile Strength

During the hydraulic fracturing process, the HF generation is dominated by the tensile failure in the conglomerate rock. The tensile strength of the rock is one of the important factors that affect the propagation of HF. Eight cases with different ratios of matrix tensile strength to gravel (Tm/Tg=0.5,0.7,0.8,0.9,1.0,1.2,1.5,2.0 and Tg=6.0MPa) were studied.It can be seen from Figure11 that, with the difference of Tm/Tg, the HF morphology on both sides of the borehole is different. Figure 12 shows that as Tm/Tg increases, L first increases and then decreases (Figure 12(a)); SRS gradually decreases (Figure 12(b)). Taking Tm/Tg=0.5 as an example, the tensile strength of the matrix is small enough to form a weak surface. Compared with the energy of fracturing the matrix, the energy accumulated by the injection fluid is large enough, and the fractures propagate strongly in the matrix. At this time, affected by the random gravel, the HF propagation is blocked, passivated, built up, and detoured. The weak surface produces multiple initiation points, resulting in multibranched fractures; the HF propagation behavior becomes unpredictable. As the Tm/Tg increases, D decreases from 0.9329 to 0.8177 (Figure 12(c)), indicating that the complexity of the HF morphology tends to decrease. Similar to the law of elastic modulus, the rate of change of D gradually decreases as Tm/Tg increases indicating that the influence of tensile strength on HF morphology is weakening.Figure 11
HF propagation morphology under different ratios of the tensile strength of the matrix to gravel.Figure 12
Geometric parameters of HF under different ratios of the tensile strength of the matrix to gravel.
(a)(b)(c)


### 4.3. Effect of Shear Strength

A simulation of the HF propagation trajectory in the conglomerate reservoir with different ratios of matrix shear strength to gravel (Sm/Sg=0.5,0.7,0.8,0.9,1.0,1.2,1.5,2.0 and Sg=60.0MPa) was performed, as shown in Figure 13.Figure 13
HF propagation morphology under different ratios of the shear strength of the matrix to gravel.As shown in Figures14(a)–14(c), the values of geometric parameters of HF under different Sm/Sg are discrete. The correlations between Sm/Sg and geometric parameters are poor. The HF morphology has no obvious changing rule affected by the shear strength. This result is obtained because the HF propagation was dominated by tensile failure. The shear strength of the matrix has only a slight influence on the behavior of HF. This is consistent with the results presented by previous studies [61].Figure 14
Geometric parameters of HF under different ratios of the shear strength of the matrix to gravel.
(a)(b)(c)


### 4.4. Effect of the Strength for the Matrix-Gravel Interface

Due to the difference in composition, the mechanical properties of the matrix, gravel, and matrix-gravel interface are significantly different. In general, quartz-rich gravel particles demonstrate the highest strength, and the interfaces exhibit the lowest strength in most conglomerate rocks. Due to stress concentration, the HF initiated and propagated along the low-strength interface generally [62]. HF was also attracted by the gravel owing to the “shielding” effect and finally deflected into the interface [22]. The matrix-gravel interface strength has an impact on HF morphology in conglomerate reservoirs. A simulation of the HF propagation trajectory in the conglomerate reservoir was operated with different ratios of matrix-gravel interface tensile strength to matrix (Tint/Tm=0.4,0.5,0.6,0.7,0.8,1.0 and Tm=6.0MPa).The interface is the instability region. Affected by the fracturing fluid pressure, gravel embedded in the conglomerate causes strong stress concentrations at the matrix-gravel interface and significantly affects the local mechanical response. Under the impact of stress disturbance, a new microfracture is easy to form at the microscopic level. For the case ofTint/Tm=0.4 shown in Figure 15, the tensile strength of the matrix-gravel interface is low and HF mainly deflects and grows along the low-strength interface. The HF propagation path becomes very tortuous and complicated in a local area. However, due to the influence of the stress distribution state, the fractures expand along the direction of the maximum stress at the whole scale. When Tint/Tm increases from 0.4 to 1.0, it becomes difficult for HF to propagate farther, and L becomes smaller (Figure 16(a)), due to the tensile strength of the matrix-gravel interface which is higher. With the increase in the tensile strength of the matrix-gravel interface, a high-strength protective layer is formed on the gravel surface; HF deterministically chooses the path with the least resistance to propagate. When the fracture encounters the gravel, HF is prone to bifurcation in the matrix and develops multibranched fractures, resulting in an increase in SRS (Figure 16(b)). Owing to the dual factors of the decrease of L and the increase of SRS, the D of the fracture network decreased from 0.9101 to 0.8532 with a small change of 6.25%, as shown in Figure 16(c). In addition, it can be seen from the geometric parameters of HF under different Tint/Tm plotted in Figure 16, Tint/Tm has a good correlation with L and SRS. The R-square is 0.9304 and 0.9134, indicating that the propagation morphology of HF is affected by the strength of the matrix-gravel interface significantly.Figure 15
HF propagation morphology under different ratios of matrix-gravel interface tensile strength to matrix.Figure 16
Geometric parameters of HF under different ratios of matrix-gravel interface tensile strength to matrix.
(a)(b)(c)


### 4.5. Effect of Minimal Horizontal Stress

Governed by the horizontal stress state, the HF always deterministically chooses the path with the least resistance to propagate [9]. To study the influence of the minimal horizontal in situ stress (σh) on the HF propagation, the models with σh as 10, 15, 20, 25, 30, 35, and 40 MPa were established. To avoid the influence of the horizontal stress difference, keeping the horizontal stress difference at 3 MPa unchanged.Geometric parameters of HF have a good correlation withσh, as shown in Figures 17(a)–17(c). The R-square between L and σh is as high as 0.9356. This shows that the fracture propagation is greatly affected by the minimum horizontal stress state. The fracture initiation pressure can be expressed as [63]
(23)pb=3σh−σH−pp+kT,where pb is the fracture initiation pressure; σh,σH are the minimum and maximum horizontal stresses, respectively; T is the tensile strength of the rock; pp is the pore pressure; and k is the coefficient of tensile strength.Figure 17
Geometric parameters of HF under different minimal horizontal stresses.According to Equation (23), the fracture initiation pressure increases linearly by a factor of 3 compared to the minimum horizontal stress. Therefore, HF tends to initiate and propagate in formations with lower minimum horizontal principal stress and lower tensile strength. When the minimum stress is lower, it is easier to fracture on weak surfaces due to the lower fracture initiation pressure. The HF is greatly affected by the gravel, making the HF growth path highly tortuous and complex (Figure 18). For the case of σh=10MPa, the D of the fracture network reaches the maximum value, as shown in Figure 17(c). Fracture initiation requires higher pressure as the minimum horizontal stress increases from 10 MPa to 40 MPa. HF in the formation is difficult to initiate, resulting in lower L and SRS. The D of the fracture network decreases from 1.0388 to 0.7692; that is, the complexity of the HF morphology decreases as a result of the significant influence of minimal stress on fracture deflection.Figure 18
HF propagation morphology under different minimal horizontal stresses.
(a)(b)(c)


### 4.6. Effect of Horizontal Stress Difference

According to the rock failure theory, cracking starts at the point where the circumferential stress on the rock exceeds the tensile strength of the rock [45]. When the rock mechanical properties are constant, the HF propagation is mainly dominated by the maximum and minimum horizontal stress states. In this paper, six model cases with different in situ stress differences between the maximum and minimum horizontal stresses (Δσ=σH−σh) were established, where σh=25MPa is set as a constant and Δσ is set as 0, 3, 6, 8, 10, and 12 MPa.In the formation with low horizontal stress difference, the strength difference has a significant effect on the propagation trajectory of HF, resulting in more complex fracture morphology. As shown in Figure19, multiple branch fractures appear on the right side of the wellbore for the case of Δσ=0MPa. When a fracture encounters gravel, multiple fractures branch at the gravel interface, and complex fractures are more likely to occur in reservoirs with lower stress difference. As Δσ increases, the diversion and branching of HF are inhibited. HFs tend to propagate in the direction of the maximum stress, and straight biwing fractures are more likely to form. When Δσ=12MPa, the HF propagation is suppressed, which seriously affects the fracturing results. The values of geometric parameters of HF under different Δσ are shown in Figures 20(a)–20(c). With the increase of Δσ, SRS and D of the fractures decrease significantly. This is consistent with the conclusions obtained from the above analysis.Figure 19
Propagation morphology of HF under different horizontal stress differences.Figure 20
Geometric parameters of HF under different horizontal stress differences.
(a)(b)(c)


## 4.1. Effect of Elastic Modulus

Elastic modulus is expressed as the ability of a rock to resist deformation under stress. Rocks with a high elastic modulus have an ability to resist deformation, which are more rigid and brittle. Thus, there is a more possibility of brittle failure for these rocks. In contrast, rocks with a low elastic modulus have greater flexibility. Previous studies have shown that, compared with the elastic modulus of gravel or matrix, “the ratio of the elastic modulus of the matrix to gravel (Em/Eg)” has a larger influence on the propagation morphology of HF [44]. In this study, models with Em/Eg were set to 0.5, 0.7, 0.8, 1.0, 1.2, 1.5, and 2.0, and the gravel elastic modulus was 22 GPa, as shown in Figure 9.Figure 9
HF propagation morphology under different ratios of the elastic modulus of the matrix to gravel.The elastic modulus of gravel or matrix has a great influence on the propagation trajectory of HF in the conglomerate reservoir. As shown in Figure9, the fracture in the specimen with Em/Eg=0.5 extends for a short distance. When Em/Eg=0.5, i.e., the lower the elastic modulus of the matrix, the higher the probability of ductile failure. HF tends to turn, kink, and branch in the vicinity of gravel with lower elastic modulus, causing local bending of the main fracture, resulting in short fracture length. With increased elastic modulus of the matrix, the heterogeneity of the elastic modulus for the rock material is weakened, and the fracture morphology should become simple. However, a positive relationship between D and Em/Eg is shown in Figure 10(c), which means the increase of the matrix elastic modulus results in an increment of the complexity of the fracture morphology. The reasons behind this situation are that the increased matrix elastic modulus results in an increased matrix stiffness and, thus, higher possibility of brittle failure. The fracture morphology becomes straight and extends further (Figure 10(a)), which is due to the decrease of local bending of the fracture on a microscale. As the elastic modulus of the matrix increases, the SRS tends to increase (Figure 10(b)). In addition to the elastic modulus, the fracture propagation is also affected by other differences in strength (tensile strength, etc.), resulting in a macroscopically distinct branching of HF. As the elastic modulus increases, the influence on HF morphology is weakening, while the influence of other strength differences is strengthening. The rate of change of D gradually decreases with the increases of Em/Eg which illustrates this point.Figure 10
Geometric parameters of HF under different ratios of the elastic modulus of the matrix to gravel.
(a)(b)(c)


## 4.2. Effect of Tensile Strength

During the hydraulic fracturing process, the HF generation is dominated by the tensile failure in the conglomerate rock. The tensile strength of the rock is one of the important factors that affect the propagation of HF. Eight cases with different ratios of matrix tensile strength to gravel (Tm/Tg=0.5,0.7,0.8,0.9,1.0,1.2,1.5,2.0 and Tg=6.0MPa) were studied.It can be seen from Figure11 that, with the difference of Tm/Tg, the HF morphology on both sides of the borehole is different. Figure 12 shows that as Tm/Tg increases, L first increases and then decreases (Figure 12(a)); SRS gradually decreases (Figure 12(b)). Taking Tm/Tg=0.5 as an example, the tensile strength of the matrix is small enough to form a weak surface. Compared with the energy of fracturing the matrix, the energy accumulated by the injection fluid is large enough, and the fractures propagate strongly in the matrix. At this time, affected by the random gravel, the HF propagation is blocked, passivated, built up, and detoured. The weak surface produces multiple initiation points, resulting in multibranched fractures; the HF propagation behavior becomes unpredictable. As the Tm/Tg increases, D decreases from 0.9329 to 0.8177 (Figure 12(c)), indicating that the complexity of the HF morphology tends to decrease. Similar to the law of elastic modulus, the rate of change of D gradually decreases as Tm/Tg increases indicating that the influence of tensile strength on HF morphology is weakening.Figure 11
HF propagation morphology under different ratios of the tensile strength of the matrix to gravel.Figure 12
Geometric parameters of HF under different ratios of the tensile strength of the matrix to gravel.
(a)(b)(c)


## 4.3. Effect of Shear Strength

A simulation of the HF propagation trajectory in the conglomerate reservoir with different ratios of matrix shear strength to gravel (Sm/Sg=0.5,0.7,0.8,0.9,1.0,1.2,1.5,2.0 and Sg=60.0MPa) was performed, as shown in Figure 13.Figure 13
HF propagation morphology under different ratios of the shear strength of the matrix to gravel.As shown in Figures14(a)–14(c), the values of geometric parameters of HF under different Sm/Sg are discrete. The correlations between Sm/Sg and geometric parameters are poor. The HF morphology has no obvious changing rule affected by the shear strength. This result is obtained because the HF propagation was dominated by tensile failure. The shear strength of the matrix has only a slight influence on the behavior of HF. This is consistent with the results presented by previous studies [61].Figure 14
Geometric parameters of HF under different ratios of the shear strength of the matrix to gravel.
(a)(b)(c)


## 4.4. Effect of the Strength for the Matrix-Gravel Interface

Due to the difference in composition, the mechanical properties of the matrix, gravel, and matrix-gravel interface are significantly different. In general, quartz-rich gravel particles demonstrate the highest strength, and the interfaces exhibit the lowest strength in most conglomerate rocks. Due to stress concentration, the HF initiated and propagated along the low-strength interface generally [62]. HF was also attracted by the gravel owing to the “shielding” effect and finally deflected into the interface [22]. The matrix-gravel interface strength has an impact on HF morphology in conglomerate reservoirs. A simulation of the HF propagation trajectory in the conglomerate reservoir was operated with different ratios of matrix-gravel interface tensile strength to matrix (Tint/Tm=0.4,0.5,0.6,0.7,0.8,1.0 and Tm=6.0MPa).The interface is the instability region. Affected by the fracturing fluid pressure, gravel embedded in the conglomerate causes strong stress concentrations at the matrix-gravel interface and significantly affects the local mechanical response. Under the impact of stress disturbance, a new microfracture is easy to form at the microscopic level. For the case ofTint/Tm=0.4 shown in Figure 15, the tensile strength of the matrix-gravel interface is low and HF mainly deflects and grows along the low-strength interface. The HF propagation path becomes very tortuous and complicated in a local area. However, due to the influence of the stress distribution state, the fractures expand along the direction of the maximum stress at the whole scale. When Tint/Tm increases from 0.4 to 1.0, it becomes difficult for HF to propagate farther, and L becomes smaller (Figure 16(a)), due to the tensile strength of the matrix-gravel interface which is higher. With the increase in the tensile strength of the matrix-gravel interface, a high-strength protective layer is formed on the gravel surface; HF deterministically chooses the path with the least resistance to propagate. When the fracture encounters the gravel, HF is prone to bifurcation in the matrix and develops multibranched fractures, resulting in an increase in SRS (Figure 16(b)). Owing to the dual factors of the decrease of L and the increase of SRS, the D of the fracture network decreased from 0.9101 to 0.8532 with a small change of 6.25%, as shown in Figure 16(c). In addition, it can be seen from the geometric parameters of HF under different Tint/Tm plotted in Figure 16, Tint/Tm has a good correlation with L and SRS. The R-square is 0.9304 and 0.9134, indicating that the propagation morphology of HF is affected by the strength of the matrix-gravel interface significantly.Figure 15
HF propagation morphology under different ratios of matrix-gravel interface tensile strength to matrix.Figure 16
Geometric parameters of HF under different ratios of matrix-gravel interface tensile strength to matrix.
(a)(b)(c)


## 4.5. Effect of Minimal Horizontal Stress

Governed by the horizontal stress state, the HF always deterministically chooses the path with the least resistance to propagate [9]. To study the influence of the minimal horizontal in situ stress (σh) on the HF propagation, the models with σh as 10, 15, 20, 25, 30, 35, and 40 MPa were established. To avoid the influence of the horizontal stress difference, keeping the horizontal stress difference at 3 MPa unchanged.Geometric parameters of HF have a good correlation withσh, as shown in Figures 17(a)–17(c). The R-square between L and σh is as high as 0.9356. This shows that the fracture propagation is greatly affected by the minimum horizontal stress state. The fracture initiation pressure can be expressed as [63]
(23)pb=3σh−σH−pp+kT,where pb is the fracture initiation pressure; σh,σH are the minimum and maximum horizontal stresses, respectively; T is the tensile strength of the rock; pp is the pore pressure; and k is the coefficient of tensile strength.Figure 17
Geometric parameters of HF under different minimal horizontal stresses.According to Equation (23), the fracture initiation pressure increases linearly by a factor of 3 compared to the minimum horizontal stress. Therefore, HF tends to initiate and propagate in formations with lower minimum horizontal principal stress and lower tensile strength. When the minimum stress is lower, it is easier to fracture on weak surfaces due to the lower fracture initiation pressure. The HF is greatly affected by the gravel, making the HF growth path highly tortuous and complex (Figure 18). For the case of σh=10MPa, the D of the fracture network reaches the maximum value, as shown in Figure 17(c). Fracture initiation requires higher pressure as the minimum horizontal stress increases from 10 MPa to 40 MPa. HF in the formation is difficult to initiate, resulting in lower L and SRS. The D of the fracture network decreases from 1.0388 to 0.7692; that is, the complexity of the HF morphology decreases as a result of the significant influence of minimal stress on fracture deflection.Figure 18
HF propagation morphology under different minimal horizontal stresses.
(a)(b)(c)


## 4.6. Effect of Horizontal Stress Difference

According to the rock failure theory, cracking starts at the point where the circumferential stress on the rock exceeds the tensile strength of the rock [45]. When the rock mechanical properties are constant, the HF propagation is mainly dominated by the maximum and minimum horizontal stress states. In this paper, six model cases with different in situ stress differences between the maximum and minimum horizontal stresses (Δσ=σH−σh) were established, where σh=25MPa is set as a constant and Δσ is set as 0, 3, 6, 8, 10, and 12 MPa.In the formation with low horizontal stress difference, the strength difference has a significant effect on the propagation trajectory of HF, resulting in more complex fracture morphology. As shown in Figure19, multiple branch fractures appear on the right side of the wellbore for the case of Δσ=0MPa. When a fracture encounters gravel, multiple fractures branch at the gravel interface, and complex fractures are more likely to occur in reservoirs with lower stress difference. As Δσ increases, the diversion and branching of HF are inhibited. HFs tend to propagate in the direction of the maximum stress, and straight biwing fractures are more likely to form. When Δσ=12MPa, the HF propagation is suppressed, which seriously affects the fracturing results. The values of geometric parameters of HF under different Δσ are shown in Figures 20(a)–20(c). With the increase of Δσ, SRS and D of the fractures decrease significantly. This is consistent with the conclusions obtained from the above analysis.Figure 19
Propagation morphology of HF under different horizontal stress differences.Figure 20
Geometric parameters of HF under different horizontal stress differences.
(a)(b)(c)


## 5. Conclusion

(1)
Based on the GPPCZ method, a 2D FSD fully coupled finite element model is established to simulate HF propagation in conglomerate reservoirs under multifield coupling(2)
Mesostructural characteristics of conglomerate (volume ratio of components, number, gradation, and spatial distribution of gravel) and local material properties of each component (elastic modulus, Poisson’s ratio, permeability properties, strength properties, etc.) have been taken into consideration during simulations to predict fracture morphology accurately. Compared with the published physical test results and numerical simulation results, the propagation behavior of fractures encountering gravel can be captured at both the microscale and macroscale, which is effective for the simulation of hydraulic fractures in conglomerate reservoirs.L and SRS are used to describe the geometry of HF propagation effectively, and the fractal theory is introduced to successfully quantify the complexity of fracture morphology(3)
It is revealed that the difference between the mechanical properties of the gravel and the matrix in the conglomerate rock will affect the geometry of HF to varying degrees. The local behavior of fracture propagation is obviously dominated by the elastic modulus, tensile strength, and the strength for the matrix-gravel interface. The shear strength of the matrix has only a slight influence on the behavior of HF. The propagation of HF at the whole scale is mainly dominated by the horizontal stress distribution state, including the minimum principal stress, horizontal stress difference. In addition, the difference in horizontal stress significantly affects the fracturing patterns (deflection, bifurcation, and penetration) when it encounters gravel. The results provide theoretical support for the prediction of HF propagation morphology and plan design of hydraulic fracturing in conglomerate reservoirs



---

*Source: 6811300-2022-07-23.xml*