A Numerical Study of Hydraulic Fracture Propagation Geometry in a Layered Shale Reservoir

An unconventional shale reservoir commonly develops multiple layers and shows strong anisotropy in mechanical properties, which has a great effect on hydraulic fracture propagation geometry. The most common mechanical properties are elastic modulus, Poisson’s ratio, tensile strength, and formation permeability. Therefore, the extended finite element method(XFEM) based cohesive zone model (CZM) is applied to analyze the effect of these mechanical properties on both the fracture propagation geometry and breakdown pressure in a layered shale reservoir after verifying the present numerical method by published analytical results. The parametric analysis indicates that the stiff or soft outer layers limit the fracture propagation width, while promoting the fracture propagation length. Higher Poisson’s ratio and formation permeability in outer layers narrow the fracture propagation width in middle layer. Poisson’s ratio contrast between different layers almost has no significant effect on the fracture propagation length and breakdown pressure. The hydraulic fracture propagation geometry presents the trend of “shorter and wider” with the increase of the tensile strength in outer layer. For asymmetric specimen with different mechanical parameters in each layer, hydraulic fracture shows an asymmetric propagation behavior, and hydraulic fracture preferentially propagates to the layer with higher elastic modulus, Poisson’s ratio, formation permeability, and low resistance.


Introduction
Hydraulic fracturing is the major technique used in the oil industry to simulate a low-porosity and low-permeability reservoir to enhance oil recovery. High-pressure fracturing fluid is pumped into boreholes to break down oil and gas reservoirs during hydraulic fracturing. The fracture propagates under an injected fluid pressure, forming a path for oil and gas to flow into the well from reservoirs. To acquire the highest possible postfracturing productivity, hydraulic fracturing is generally used to form long and penetrating fractures and increase the stimulated reservoir volume to establish complex fracture networks. Hydraulic fracture propagation geometry (length, width, and height) is an important factor in evaluating fracture complexity. As shown in Figure 1 [1,2], an unconventional shale reservoir, as a layered sedimentary rock, commonly develops multiple layers and shows strong mechanical heterogeneity. The mechanical heterogeneity of shale reservoir rock makes the stress and displacement distribution quite different from a homogeneous reservoir, thus affecting the hydraulic fracture propagation characteristics. Hence, determining the influencing factors of mechanical heterogeneity, such as elastic modulus, Poisson's ratio, tensile strength, and formation permeability on the hydraulic fracture propagation geometry, is critical to understanding the fracture complexity in shale reservoirs.
Considerable efforts have been devoted to establishing various analytical models to investigate the fracture propagation geometry during hydraulic fracturing. For example, Khristianovic et al. [3,4] have proposed and developed the KGD hydraulic fracture model, which is often used to describe the hydraulic fracture geometry initiated from a wellbore. Nordgren [5] developed the PKN hydraulic fracture model proposed by Perkins and Kern [6], which included fluid loss and fracture volume change in the continuity equation. In addition, for multilayered formations with arbitrary in situ horizontal stress distribution, a general semianalytical solution to obtain the vertical fracture geometry was proposed by Fung et al. [7]. Those analytical models provide important help for investigating the hydraulic fracture propagation geometry. However, these analytical models are only limited to studying some simple factors. Investigating the influence of mechanical heterogeneity on the hydraulic fracture propagation geometry is challenging.
Some experiments were also conducted to study the effect of mechanical heterogeneity on hydraulic fracture propagation behaviors. The laboratory experiments and analytical and numerical studies were carried out by Llanos et al. [8] to investigate how hydraulic fracture propagates through an orthogonal discontinuity. Teufel et al. [9] carried out laboratory experiments in which three-layer rock specimens were hydraulically fractured under a uniaxial compression stress state and demonstrated that two distinct geologic conditions could inhibit or contain hydraulic fracture growth in layered rock. Altammar et al. [10] performed an experimental study in layered specimens that directly observed fracture growth near the material interface. The results provide fundamental insights on fracture propagation behavior in heterogeneous porous media. Fisher and Warpinski [11] presented the data from microseismic and tiltmeter monitoring technologies that the height-growth mechanisms would likely be dominated by the layering and heterogeneities of both properties and stress. These experimental results provide an important understanding of some basic hydraulic fracture propagation behaviors. However, experimental investigations of the fundamental mechanisms of hydraulic fracture behavior are complicated and expensive, and many tests are required to obtain useful results. In addition, it is generally difficult to repeat the experiments [12][13][14]. These limitations strongly restrict the usefulness of results obtained from experimental studies.
Numerical methods have been predominant in modeling crack propagation because of their high generality and flexibility in modeling structures with complex geometries and various boundary and loading conditions [15][16][17]. Therefore, numerical simulation techniques were applied to study hydraulic fracturing. For example, Liu et al. [18] investigated the effects of elastic modulus, Poisson's ratio, in situ stress difference, injection rate, and perforation parameters on the fracture propagation length and width using the extended finite element method (XFEM). However, their model is a homogeneous numerical model rather than a heterogeneous model. For a layered unconventional gas reservoir, many researchers have also studied hydraulic fracture propagation behaviors using the XFEM [19][20][21][22]. In addition, Wu et al. [23] studied the influence of heterogeneity on hydraulic fracturing using the combined finite-discrete method. They mainly studied the effect of modulus contrast on hydraulic fracture propagation behavior in layered reservoirs without considering other factors, such as tensile strength, Poisson's ratio, and formation permeability. Behnia et al. [24] presented a boundary element method based on the higher order displacement discontinuity formulation to solve the general problem of hydraulic fracture propagation in layered formations. Guo et al. [25] studied the hydraulic fracture propagation path in heterogeneous formation under anisotropic in situ stress using the phase field method, in which the influence of reservoir elastic modulus was mainly discussed. Aimene et al. [26] proposed the material point method (MPM) geomechanical tool that combines an anisotropic damage mechanics model and interface model to describe the hydraulic fracture propagation height in layered shale and tight gas rock. Most studies focus on the effect of mechanical heterogeneity on hydraulic fracture propagation path, while few studies focus on fracture propagation geometry.
The mechanical heterogeneity of reservoir rock affects the hydraulic fracture propagation path and the fracture propagation geometry. This study established a fully fluidsolid coupling model of layered shale using the XFEMbased cohesive zone model (CZM) to simulate the hydraulic fracture propagation behaviors in a layered shale reservoir. Some typical examples were used to illustrate the rationality of the presented methods, and the obtained results were compared with the asymptotic analytical KGD solution. The effects of elastic modulus, Poisson's ratio, tensile strength, and formation permeability on the fracture propagation geometry during hydraulic fracturing were analyzed. The study is organized as follows. In Section 2, the numerical results of the homogeneous reservoir are compared with the KGD solution to verify the reliability of the model. The fracture propagation geometry characteristic under different influencing factors is discussed in detail in Section 3. Finally, the conclusions of this study are presented in Section 4.

Nodal Enrichment Functions.
To obtain more accurate stress and displacement fields at the crack tip, singular elements are used at the crack tip when the finite element method is used to study the crack problem. The stress concentration at the crack tip needs to be processed by the element grid encryption, so preprocessing of the finite element method is relatively complicated [15]. In order to overcome the disadvantages of the finite element method, the XFEM was first introduced by Belytschko and Black [27] to remove the mesh refinement and allow discontinuities to cross elements. XFEM is an extension of the conventional finite element method according to the concept of partition of unity by Melenk and Babuska [28], which allows local enrichment functions to be easily incorporated into a finite element approximation. In commercial software ABAQUS [29], the hydraulic fracture propagation characteristics can be investigated by the cohesive element method, while the cohesive element method can only make the hydraulic fractures propagate along the cell boundary. The hydraulic fracture propagation path is complicated in the complex geo-stress field and heterogeneous reservoir. Therefore, combining the finite element and cohesive element with the virtual node method can simulate the arbitrary propagation path of hydraulic fracture [19,30]. Furthermore, virtual nodes representing the degree of freedom of pore water pressure are introduced at the edge of each enrichment element to simulate the fluid flow on the surface of the element. The virtual nodes are superimposed with the original real nodes to represent the displacement discontinuity and fluid pressure of the fracture element [20][21][22]. For fracture analysis, the special enriched function typically consists of three parts: a conventional shape function, a near-tip asymptotic function, and a discontinuous function. The conventional shape functions are used to determine the continuous displacement field. The near-tip asymptotic functions are used to capture the singularity around the crack tip, while the discontinuous functions represent the jump in displacement across the crack surfaces. Hence, the approximation for a displacement vector function u with the partition of unity enrichment can be given by Equation (1) [31].
where N IðxÞ is the usual nodal shape functions, u I is the usual nodal displacement vector associated with the continuous part of the finite element solution, a I is the nodal enriched degree of freedom vector, HðxÞ is the discontinuous jump function across the crack surfaces, b α I is the nodal enriched degree of freedom vector, and F α ðxÞ is the elastic asymptotic crack-tip function.
On the right-hand side of the Equation (1), the first term u I is applicable to all the nodes in the model, the second term HðxÞa I is valid for nodes whose shape function support is cut by the crack interior, and the third term ∑ 4 α=1 F α ðxÞb α I is used only for nodes whose shape function support is cut by the crack tip. The calculation of HðxÞ and F α ðxÞ could be given in the previous researches [29].

Cohesive Zone Method
2.2.1. Damage Initiation Criterion. Damage initiation refers to the onset of degradation of the cohesive response at an enriched element. The degradation begins when the stresses or the strains satisfy specified crack initiation criteria. The maximum principal stress criterion is used in this study as the criteria for the initial damage on fracture, and damage is assumed to initiate when the maximum principal stress ratio reaches 1.0. Thus, this criterion can be represented as where f is the maximum principal stress ratio, σ o max represents the maximum allowable principal stress, σ max represents the stress applied on the fracture surface, and the symbol < > represents the Macaulay bracket with the usual interpretation (i.e., hσ max i = 0 if σ max < 0 and hσ max i = σ max if σ max ≥ 0) and is used to signify that a purely compressive deformation or stress state does not initiate damage. In addition, the crack always propagates along the direction perpendicular to the maximum principal stress.

Damage Evolution
Law. The damage evolution law describes the rate at which cohesive stiffness is degraded once the corresponding initiation criterion is stratified. To quantify damage evolution, a scalar damage variable D with an initial value of 0 is introduced, which represents the overall damage at the intersection between the crack surfaces and the edges of cracked elements. After the damage initiation criterion is reached, the value of D monotonically evolves from 0 to 1.0 as further loading. D = 0 means that the cohesive stiffness is not completely degraded, and no new cracks are generated. D = 1:0 means that the cohesive stiffness is completely degraded to 0, leading to the expansion of the crack or generation of an additional crack. The normal and shear stress components of the traction-separation model are affected by the damage according to where T n , T s , and T t are the normal and shear stress components predicted by the elastic traction-separation behavior for the current separations without damage. The BK law model is described by Benzeggagh and Kenane [32], which can combine energy release rates in 3 Geofluids Mode I, Mode II, and Mode III into a power-law relationship of a single scalar fracture criterion. This criterion is suitable for the critical fracture energy of similar rock material along the first and the second shear directions. The equivalent fracture-energy release rate can be expressed by where G c I and G c II denote the critical fracture energy in the normal and first shear directions, respectively; G I , G II , and G III denote the work done by the tractions and their conjugate relative displacements in the normal, first and second shear directions, respectively; η denotes the contribution of shear mode ratio to critical fracture energy; and the value of η is assumed to be 1.5 in this study.

Fluid
Flow inside the Fracture. The cohesive element method with phantom nodes can also be extended to model hydraulic fracture. An additional phantom node with pore pressure degrees of freedom is introduced on the edges of each enriched element to model the fluid flow within the crack element surfaces in conjunction with the phantom nodes that are superposed on the original real nodes to represent the discontinuities of displacement and fluid pressure in a cracked element. The flow patterns of the pore fluid in the cracked elements are shown in Figure 2(a) [29]. The fluid flow is continuous, indicating both tangential and normal flow within and across the cracked element surfaces and the rate of opening of the cracked element surfaces. The fluid pressure on the cracked element surfaces contributes to the traction-separation behavior of the cohesive segments in the enriched elements, which enables the modeling of hydraulic fracture. Many hydraulic fracturing fluids exhibit power-law rheological behavior and temperature-related properties. Slick water is often used as a fracturing fluid in most low permeability. In order to avoid additional complexity, the fluid is assumed to be incompressible with Newtonian rheology.
The tangential flow rate q inside the fracture can be determined by the pressure gradient to the fracture width for a Newtonian fluid of viscosity, which can be expressed as Equation (7). According to Equation (7), the pressure drop along the fracture can be determined by the local flow rate and local fracture width [33].
where d is the fracture opening width (as shown in Figure 2 (a)), μ is the injection fluid viscosity, and p f is the injection fluid pressure along the fracture surface. According to Reynold's lubrication theory, the conservation of fluid mass inside the fracture can be expressed as where q t and q b are the normal flow rate of injection fluid leaking through the top and bottom surface of the fracture into the porous medium, respectively. A fluid leak-off coefficient for the pore fluid material is defined to determine the normal flow. This coefficient defines a pressure-flow relationship between the phantom nodes located at the cracked element edges and cracked element surfaces. The fluid leak-off coefficients can be interpreted as the permeability of a finite material layer on the cracked element surfaces, as shown in Figure 2(b) [29].
The normal flow is defined as where c t and c b are the fluid leak-off coefficients for the top and bottom surface of the fracture, respectively, and p t and p b are the pore fluid pressures on the top and bottom surface of the fracture, respectively.

Fluid-Solid Coupling in a Porous
Medium. The rock matrix is conventionally considered a multiphase material, consisting of a solid skeleton and connected or unconnected pores saturated with formation fluid. Therefore, the total stresses include the effective stresses generated by solid skeleton and the pore pressure stored in the formation fluid.
Based on the pioneering work of Biot [34], in a porous medium filled with fluid, the relation of the total stresses 4 Geofluids and the effective stresses can be expressed by Equation (10), and the effective stresses govern the deformation and failure of the rock.
where σ ij is the total stress, σ ij ′ is the effective stress, α is Biot's coefficient which related to the rock property, and p m is the pore pressure.
According to the principle of virtual work, the stress equilibrium for the solid phase of rock matrix under its current configuration can be expressed by Equation (11) [29].
where σ ij ′ is the effective stress tensor, I is the unit matrix, V is the solid volume, t is the surface traction vector per unit area, f is the body force vector per unit volume, S is the surface area controlled by surface traction, δε is the virtual strain rate, and δv is the virtual velocity vector. The porous medium is modeled in ABAQUS by attaching the finite element mesh to the solid phase and then allowing liquid to flow through the mesh. A continuity equation is required for the fluid: the time rate of change of the total liquid mass in the control volume V is equal to the mass of liquid crossing the surface S per unit time. The liquid mass continuity equation is given [29].
where ρ f is the fluid density, n f is the porosity of the permeable formation, n is the outward normal to surface S, and q m is the fluid velocity vector in the porous medium. According to Darcy's Law, the fluid flow in a porous medium can be expressed by where g is the gravity acceleration vector, k is the matrix permeability tensor of the porous medium, and X is a spatial coordinate vector.

Model Construction and Verification
3.1. Model Settings. According to the hydraulic fracturing of laboratory and field observations, the hydraulic fracture will propagate in the three-dimensional direction (length, width, and height) inside the reservoir under pump fracturing fluid pressure. However, the effect of different factors on the hydraulic fracture propagation length and width is studied in this study without considering the fracture height various. A two-dimensional stress-pore pressure fully coupled plain strain numerical model was established in this study, as shown in Figure 3.
The model size was assumed to be 90 m × 90 m, which was large enough to minimize boundary effects, and an original horizontal hydraulic fracture was predefined in the model center. The numerical model was divided into two outer layers (about 30 m) and one middle layer (about 30 m). The effect of mechanical properties on the hydraulic fracture propagation behaviors was explored by assigning different mechanical properties to different layers. The maximum horizontal stress σ H and minimum horizontal stress σ h were applied to the numerical model along the x and y directions, respectively. Displacement of all the outer boundaries was equal to zero along the direction that was perpendicular to its surface. Additionally, constant pore pressure was imposed on the outer boundaries, and incompressible Newtonian fluid was injected into the predefined hydraulic fracture at a constant rate.
All the basic input parameters for the simulation cases are summarized in Table 1. The reservoir rock was modeled with elements of two-dimensional CPE4P (4-node bilinear displacement and pore pressure) [29]. The numerical simulation of hydraulic fracturing can be divided into two-step. The first step is to obtain the equilibrium state by applying boundary displacement, pore pressure, and external load to the numerical model. The second step is to inject a constant rate of fracturing fluid into the predefined hydraulic fracture to propagate the hydraulic fracture.

Numerical Model Verification.
The capabilities of the proposed method for modeling hydraulic fractures were verified in this section. A numerical example was used to prove the correctness and accuracy of the proposed model, and then, the results were compared with the analytical results. The maximum principal stress and the vertical (y direction) displacement distributions around the hydraulic fracture in homogeneous reservoir are shown in Figure 4. In a homogeneous reservoir, the outer and middle layers have the same mechanical properties (as shown in Figure 3). It can be seen  Figure 4(a) that stress concentration occurs at crack tips under the joint effect of in situ stress, pore pressure, and fluid injection pressure. As the fracturing fluid pressure increases, the predefined hydraulic fracture was initiated from tips and then propagates symmetrically near the injection point. In addition, it can be observed from Figure 4(b) that with the increase of fracturing fluid pressure, the upper and lower hydraulic fracture surfaces are separated from each other in the y direction, and the displacement on both sides of the fracture is symmetrically distributed. The fracture opening is the largest at the injection point, and its value decreases with a long distance from the injection point. The above-mentioned analysis indicates that the simulation method proposed in this study is suitable to describe the stress and displacement distributions during hydraulic fracturing.
The fracture propagation width and injection pressure at the injection point with the injection time were recorded during the numerical simulations. The comparison of numerical result with the analytical solutions from the KGD model [3,4] is shown in Figure 5. It can be found that the numerical results of the fracture propagation width are agreed well with the analytical results, and the injection pressure decay after breakdown pressure is slightly larger than that of the analytical results. In addition, the injection pressure increases sharply in the early stage of hydraulic fracturing, and the fracture propagation width at the injection point increases gradually with the injection time. At this stage, the hydraulic fracture cannot initiate because the injection pressure does not reach the tensile strength of reservoir rock. The hydraulic fracture begins to initiate when its value reaches the tensile strength of reservoir rock. Furthermore, in addition to maintaining the continuous expansion of hydraulic fracture, the fracturing fluid injection pressure also needs to support the fracture width and overcome the viscous resistance of fracturing fluid, resulting in a gradual decrease in the injection pressure after the breakdown pres-sure. At the steady fracture propagation stage, the injection pressure remains unchanged with the injection time. During hydraulic fracturing, the fracture propagation width and length will maintain an increasing trend. When the injection rate of fracturing fluid is constant, the increase of fracture propagation width will result in decreased injection pressure. However, when the fracturing fluid is full of the fracture again, the injection pressure will gradually increase until it propagates hydraulic fracture. Therefore, the injection pressure and fracture propagation width show a "sawtooth" trend during the steady fracture propagation stage, as shown in Figure 5.

Numerical Results and Analysis
This section uses several parameters (elastic modulus, Poisson's ratio, tensile strength, and formation permeability) to analyze the fracture propagation geometry. However, the difference in mechanical properties between different layers results in complex fracture propagation characteristics [35,36]. Daneshy [37] has indicated that fracture slides along the interface because of existing shear stresses (caused by material difference) when the fracture encounters a weak interface. In this study, to better understand the effect of mechanical properties on the fracture propagation geometry, it is assumed that the interface has good cohesive properties while ignoring fracture sliding along the interfaces. Hence, the hydraulic fracture will directly propagate through the interface without bifurcation or the "T" type failure model.

Effect of Elastic Modulus
Contrast. Many studies have investigated the effect of elastic modulus on hydraulic fracture behaviors by using numerical or experimental methods. However, these studies all focus on the hydraulic propagation path, and few studies have studied the influence of modulus contrast in different layers on the hydraulic fracture geometry [10,21,23,25,38]. To better understand the impact of multilayered systems on fracture propagation geometry, the impact of modulus contrasts on fracture propagation geometry has been investigated using the following specimens and strategies. Case 1 is a specimen with a soft middle layer bounded by two stiff outer layers, case 2 is a specimen with a stiff middle layer bounded by two soft outer layers, and case 3 is a specimen with a middle layer bounded by a stiffer outer layer on one side and a softer outer layer on the other side. To control the variables of the numerical model, other parameters except elastic modulus in this section are held constant. The fracturing fluid injection time is 250 sec.  Geofluids distribution has no obvious change during the hydraulic fracture propagation. Figures 6(b) and 6(c) show the evolution of stress when the ratio of elastic modulus of the middle layer to that of outer layers is 1 : 4. The results indicate that the maximum principal stress distribution is significantly affected by modulus contrast and thus affecting the hydraulic fracture propagation characteristics. Therefore, studying the effect of modulus contrast between different layers is necessary on the fracture propagation geometry.
The injection pressure and fracture propagation width at the injection point in a heterogeneous shale reservoir are shown in Figure 6(d). It can be found that the variation of the injection pressure and fracture propagation width is similar to that in a homogeneous reservoir when the hydraulic fracture propagates in layer 2 ( Figure 5). However, when the hydraulic fracture gradually propagates to outer layers, the injection pressure and fracture propagation width shows a decreasing trend with the injection time, and then, the injection pressure tends to be stable. The reason can be explained by the results shown in Figures 6(a) and 6(c). When the hydraulic fracture encounters an interface with different material properties (fracture propagation halflength is 15 m), the maximum principal stress in heterogeneous reservoir is higher than that in a homogeneous reservoir. Thus, a lower injection pressure is needed to facilitate fracture propagate. Furthermore, hydraulic fracture with lower injection pressure results in a narrower fracture propagation width.
The above analysis shows that the ratio of elastic modulus of the middle layer to that of outer layers significantly influences the stress distribution, injection pressure, and fracture propagation geometry during hydraulic fracturing. Therefore, the effects of modulus contrast on fracture propagation geometry under different conditions will be analyzed in detail in the following subsections.     Figure 7 shows the influence of modulus contrast on the fracture propagation geometry and injection pressure under a heterogeneous shale reservoir. As shown in Figure 7(a), the fracture propagation width gradually narrows as the elastic modulus of the outer layer increases when hydraulic fracture propagates in layer 2. The result shows that the outer layers with higher elastic modulus can inhibit the fracture propagation width of the middle layer. When hydraulic fracture propagates in outer layers, the greater ratio of elastic modulus of the middle layer to that of outer layers will cause a more obvious change of fracture propagation width at the interface. In addition, the fracture propagation width in the outer layers is narrower than that in the middle layer. According to the cubic law, the fracturing fluid pressure within the crack surface will increase as the fracture propagation width decreases under the constant injection rate, as shown in Figure 7(b). High fracturing fluid pressure will cause the continuous forward propagation of hydraulic fracture, as shown in Figure 7(a). Under the same injection time, the fracture propagation length increases gradually as the elastic modulus in outer layers increase.
The effect of modulus contrast on the breakdown pressure is presented in Figure 8. Modulus contrast is important for the stress and displacement distribution during hydraulic fracturing in a heterogeneous shale reservoir (as shown in Figures 6(b) and 6(c)). Thus, different modulus contrast has various breakdown pressure. As shown in Figure 8, the breakdown pressure increases as the ratio of elastic modulus of the middle layer to that of the outer layers increases.

Case 2: A Specimen with a Stiff Middle Layer Bounded by Two Soft Outer Layers.
In this case, a hydraulic fracture will propagate in a three-layer specimen from the stiff middle layer to the soft outer layer. The elastic modulus of the middle layer (layer 2) remains unchanged (E 2 = 50 GPa), and the elastic modulus of outer layers is changed to investigate the effect of modulus contrast on the fracture propagation geometry. The elastic modulus of outer layers (layers 1 and 3) gradually decreases (10 GPa and 20 GPa).

Geofluids
When the elastic modulus of the middle layer is higher than that of the outer layers, the influence of modulus contrast on the fracture propagation geometry is shown in Figure 9. The fracture propagation half-width curves no longer present an inverted U shape like a homogeneous shale reservoir. In a heterogeneous shale reservoir, the increase of elastic modulus in the middle layer will increase the antideformation ability in the middle layer, making the hydraulic propagation width in the middle layer narrower than that in a homogeneous shale reservoir. In addition, when the hydraulic fractures penetrate the middle layer and propagate in the outer layers, the fracture propagation half-width will first increase and then decreases gradually under the influence of modulus contrast. According to the analysis in case 1, the fracturing fluid pressure within the crack surface will increase as the fracture propagation width decreases, resulting in a longer hydraulic fracture propagation path. These results indicate that the outer layers with lower elastic modulus inhibit the fracture propagation width and promote the increase of fracture propagation length. In addition, the modulus contrast inhibits the development of fracture propagation width compared with the results of case 1 and case 2. In case 1, the effect of modulus contrast on fracture propagation width is more significant in outer layers (layers 1 and 3). However, in case 2, the effect of modulus contrast on Interface     9 Geofluids fracture propagation width is more significant in the middle layer (layer 2).

Case 3: A Specimen with a Middle Layer Bounded by a Stiffer Outer Layer on One Side and a Softer Outer Layer on the Other Side.
To investigate the hydraulic fracture propagation characteristic in the spatially asymmetric layered model, an asymmetric specimen with a middle layer bounded by a stiffer outer layer on one side and a softer outer layer on the other side is also considered. In this case, the elastic modulus of outer layers (layers 1 and 3) is assumed to be 10 GPa and 20 GPa, respectively. The elastic modulus of the middle layer (layer 2) is assumed to be 15 GPa. Figure 10(a) shows the fracture propagation geometry for the asymmetric specimen and the homogeneous shale reservoir specimen. The two cases exhibit noticeable differences fracture propagation length and width. Furthermore, the simulation results in two cases are compared to better understand the effect of modulus contrast on the fracture propagation geometry. Under the influence of modulus contrast, the stress distribution of asymmetric specimens will change significantly during hydraulic fracturing, resulting in the asymmetric-induced strain in reservoir rock. Compared with the results in homogeneous shale reservoir specimen, the fracture propagation length and width curves are asymmetrically distributed in relative to the model center. In addition, the softer and stiffer outer layers (layers 1 and 3) can inhibit the increase of fracture propagation width and promote the increase of fracture propagation length. Figure 10(b) shows the fracture propagation geometry under different injection time. In the early stage of hydraulic fracturing (t = 49:93 sec), hydraulic fracture propagates symmetrically to both sides near the model center. The stress intensity factor (SIF) increases to infinity as the hydraulic fracture approaches to the lower elastic modulus layers.
Hence, the hydraulic fracture will preferentially propagate to layer 1 with injection time. After the hydraulic fracture tip is located in layer 1, the SIF decreases and the propagation speed of hydraulic fracture gradually decreases. Then, the hydraulic fracture begins to propagate to layer 3 and gradually penetrates the layer interface between layer 2 and layer 3. According to the linear elastic fracture mechanics theory, the energy required for fracture propagation in layers with higher elastic modulus is lower than that in layers with lower elastic modulus. Hence, the fracture propagation speed in layer 3 is larger than that in layer 1 after the hydraulic fracture propagation in the outer layers. The numerical results have also been demonstrated by Hisanao Ouchi [39,40], indicating that the fracture in layered materials propagates preferentially towards a stiffer layer. Therefore, this analysis suggests that the modulus contrast between different layers is an important parameter for designing hydraulic fracturing.

Effect of Poisson's Ratio Contrast.
Unconventional reservoirs including coal seams and tight sandstones typically have low permeability and are located in heterogeneous geological environments. Hydraulic fracturing is an important method to enhance unconventional reservoirs production. However, hydraulic fracture is quite complicated due to the mechanical parameters (e.g., elastic modulus and Poisson's ratio), loading conditions, and drilling technology. Because Poisson's ratio of rocks spans a narrow range, particularly in a layered unconventional gas reservoir, its effect on the fracture propagation geometry has not received much attention. Poisson's ratio is assumed to be constant, which is unreasonable for heterogeneous shale reservoirs. Xie et al. [41] have reported that with the increase of angle between the loading direction and bedding direction, Poisson's ratio of Longmaxi laminated shale displays a "U-shaped" anisotropic model characterized by a first decrease and then  10 Geofluids increases. The effect of Poisson's ratio on the fracture propagation geometry has not been systematically studied. Hence, to better understand the impact of Poisson's ratio contrast on fracture propagation geometry, two potential cases of the fracture propagation geometry are investigated on trilayer specimens in this study. Case 1 is a specimen with a middle layer (lower Poisson's ratio) bounded by two higher Poisson's ratio outer layers, and case 2 is a specimen with a middle layer bounded by a lower Poisson's ratio outer layer on one side and a higher Poisson's ratio outer layer on the other side. Additionally, the other parameters except Poisson's ratio in this section are held constant.

Case 1: A Specimen with a Middle Layer (Lower Poisson's Ratio) Bounded by Two Higher Poisson's Ratio
Outer Layers. In this case, a hydraulic fracture will propagate from a middle layer with lower Poisson's ratio to an outer layer with higher Poisson's ratio. Gercek [42] has reported that Poisson's ratio for rocks varies considerably but it is never greater than 0.5. Poisson's ratio of the middle layer (layer 2) remains unchanged (υ 2 = 0:1), and Poisson's ratio of outer layers (layers 1 and 3) is gradually increased (0.1, 0.2, and 0.3). The fracturing fluid injection time is 250 sec. Figure 11(a) shows the influence of Poisson's ratio on fracture propagation geometry. Figure 11(a) shows that when the hydraulic fracture propagates in the middle layer, the fracture propagation half-width gradually decreases as Poisson's ratio in outer layers increases. Therefore, the outer layers with a higher Poisson's ratio can inhibit the fracture propagation width of the middle layer. After the hydraulic fracture propagation in outer layers, Poisson's ratio almost has no significant effect on fracture propagation geometry. Anisotropic Poisson's ratio can change the fracture propagation geometry. An accurate laboratory measurement of Poisson's ratio is applied to accurately predict the fracture propagation geometry, especially for hydraulic fracturing applications in a layered shale reservoir.
The effect of Poisson's ratio on the breakdown pressure is shown in Figure 11(b). As illustrated in Figure 11(a), the breakdown pressure increases as Poisson's ratio of outer layers increases, indicating that outer layers with a higher Poisson's ratio will need a greater fracturing fluid injection pressure to form a new fracture. Eaton [43] has also investigated these numerical results and indicated that the fracture breakdown pressure and Poisson's ratio increase with depth.

Case 2: A Specimen with a Middle Layer Bounded by Lower Poisson's Ratio Outer Layer on One Side and Higher
Poisson's Ratio Outer Layer on the Other Side. In this case, a hydraulic fracture will propagation in a three-layer specimen with different Poisson's ratios. Poisson's ratio of the outer layers (layers 1 and 3) is assumed to be 0.1 and 0.3, respectively. Furthermore, Poisson's ratio of the middle layer is assumed to be 0.2 to investigate the influence of Poisson's ratio on the fracture propagation geometry. The fracturing fluid injection is assumed to be 300 sec. Figure 12 illustrates the fracture propagation geometry for the asymmetric specimen (with different Poisson's ratios in each layer) and the homogeneous shale reservoir specimen. The fracture propagation geometry presents the same inverted U shape in these two cases. However, in the asymmetric specimen, the fracture propagation geometry curves are asymmetrically distributed in relative to the model center under the influence of Poisson's ratio contrast. In addition, it can be seen that Poisson's ratio contrast almost has no significant effect on the fracture propagation width, as illustrated in Figure 12(a). Figure 12(b) shows the fracture propagation geometry for asymmetric specimen under different injection times. When the injection time is smaller than 150.4 sec, the hydraulic fracture propagates symmetrically to both sides  11 Geofluids near the model center. When hydraulic fracture gradually penetrates layer 1, its propagation speed in layer 1 will decrease with the injection time. However, the hydraulic fracture will gradually propagate to layer 3 and finally penetrate the interface between layer 2 and layer 3. In addition, the hydraulic fracture propagation speed and length in layer 3 are greater than those in layer 1. This conclusion is also consistent with the numerical results obtained by the study of Li et al. [44], who have indicated that the fracture propagation length increases as Poisson's ratio increases.

Effect of Tensile Strength
Contrast. In this section, two different cases are selected to better understand the influence of reservoir rock tensile strength contrast on the fracture propagation geometry. Case 1 is a specimen with a middle layer (lower tensile strength) bounded by two outer layers with the higher tensile strength, and case 2 is a specimen with a middle layer bounded by an outer layer with the lower tensile strength on one side and a higher tensile strength outer layer on the other side. In addition, the other parameters except tensile strength in this section are held constant. The influence of tensile strength on the fracture propagation geometry in heterogeneous shale reservoirs is presented in Figure 13. It can be seen from Figure 13 that the hydraulic fracture propagates symmetrically to both sides near the model center during the total hydraulic fracturing. When hydraulic fracture propagates in layer 2, fracture half-width gradually increases as the tensile strength in outer layers increases. However, when hydraulic fracture propagates in outer layers, the fracture propagation width is narrower than under the homogeneous shale reservoir. In addition, the fracture propagation length gradually decreases slightly as the tensile strength in outer layers increases. According to the maximum principal stress criterion, outer layers with higher tensile strength will require greater energy for fracture initiation and propagation, resulting in a greater fracturing fluid injection pressure required to form a new  12 Geofluids fracture. Therefore, the outer layer with higher tensile strength will limit the increase of fracture propagation length. As a result, when hydraulic fracture propagates from the layer with lower tensile strength to the layer with higher tensile strength, the fracture propagation geometry curves present the feature of "shorter and wider" as the tensile strength in the outer layer increases.

Case 2: A Specimen with a Middle Layer Bounded by an
Outer Layer with the Lower Tensile Strength on One Side and a Higher Tensile Strength Outer Layer on the Other Side. In this case, the tensile strength of the middle layer (layer 2) is assumed to be 2.0 MPa, and the tensile strength of outer layers (layers 1 and 3) is assumed to be 1.0 MPa and 4.0 MPa, respectively. The fracturing fluid injection time is 300 sec. The fracture propagation geometry for the asymmetric specimen (with different tensile strengths in each layer) and the homogeneous shale reservoir specimen are presented in Figure 14. Figure 14(a) shows that the fracture propagation geometry presents inverted U shape in the two cases. However, due to the influence of tensile strength contrast between the middle and outer layers, the hydraulic fracture no longer extends symmetrically from the injection point to both sides in the asymmetric specimen. Additionally, the fracture propagation geometry in asymmetric specimens does not change significantly compared with that under a homogeneous shale reservoir. Figure 14(b) presents the fracture propagation geometry of asymmetric specimens under different injection times. When the injection time is smaller than 151.3 sec, the hydraulic fracture gradually propagates symmetrically to both sides near the model center in layer 2. When the hydraulic fracture gradually reaches the interface between the middle and outer layers, the hydraulic fracture will preferentially propagate to the layer of smaller resistance with the injection time. This result indicates that the fracture propagation length in layer 1 is greater than that in layer 3.

Effect of Formation Permeability
Contrast. Different from the conventional reservoirs, shale, as a type of layered sedimentary, is typically characterized by low porosity and extremely-low permeability. Due to the presence of bedding in a shale reservoir, the permeability of shale is strongly anisotropic between directions parallel and perpendicular to bedding. Furthermore, shale permeability anisotropy is also found between parallel to bedding directions because different horizontal directions may have experienced different tectonic movements [2]. This section investigates the influence of formation permeability contrasts on the fracture propagation geometry by using the following specimens and strategies. Case 1 is a specimen with a middle layer (lower formation permeability) bounded by two outer layers with higher formation permeability, and case 2 is a specimen with a middle layer bounded by an outer layer with lower formation permeability on one side and an outer layer with higher formation permeability on the other side. In addition, the other parameters except formation permeability in this section are held constant.

Case 1: A Specimen with a Middle Layer (Lower Formation Permeability) Bounded by Two Outer Layers
with Higher Formation Permeability. In this case, the hydraulic fracture will propagate from the middle layer with lower formation permeability to an outer layer with higher formation permeability. Therefore, the formation permeability of the middle layer (layer 2) remains unchanged (K 2 = 2:5 mD), and the formation permeability of outer layers (layers 1 and 3) is assumed to be 2.5 mD, 6.0 mD, and 10.0 mD, respectively. The fracturing fluid injection time is 300 sec. The influence of formation permeability on   13 Geofluids the fracture propagation geometry is shown in Figure 15. It can be found that when hydraulic fracture propagates in layer 2, the fracture propagation width decreases as the formation permeability in outer layers increases and is narrower than that in a homogeneous shale reservoir. When hydraulic fracture propagates in outer layers, the fracture propagation width is wider than in a homogeneous reservoir. The fracture propagation length increases slightly as the formation permeability in outer layers increases, indicating that the fracture propagation length is not sensitive to the formation permeability.

Case 2: A Specimen with a Middle Layer Bounded by a
Lower Formation Permeability Outer Layer on One Side and a Higher Formation Permeability Outer Layer on the Other Side. In this case, the formation permeability in the middle layer (layer 2) is assumed to be 2.5 mD and the formation permeability in outer layers (layers 1 and 3) is assumed to be 0.01 mD and 6.0 mD, respectively. The fracturing fluid injection time is 300 sec. The specimen of formation permeability variation can influence the fluid leak-off. Furthermore, the pore pressure will change in the vicinity of the fractures as the fluid leaks into the formation during hydraulic fracturing. The pore pressure variations will change the effective stresses and fracture propagation geometry. Figure 16 presents the fracture propagation geometry for the asymmetric specimen (with different formation permeabilities in each layer) and the homogeneous shale reservoir specimen, respectively. The fracture propagation geometry curves in these two specimens have the same shape (inverted U shape). In contrast, the hydraulic fracture in asymmetric specimen no longer presents symmetric distribution near the model center under the formation permeability contrast between different layers. Furthermore, it can be observed from Figure 16 that after the hydraulic fracture propagates to the interface between layer 2 and layer 1, it will stop at the interface with lower formation permeability. However, the hydraulic fracture will continue propagating towards layer 3 with the injection time. Higher formation permeability can induce greater fluid leak-off and increase the pressure around the fracture tip [45]. Higher pressure induces more tensile effective stress and hence more tensile effective traction at the fracture tip. According to the maximum principal stress damage initiation criterion, meeting the criteria of fracture propagation in a higher formation permeability layer is easy. Hence, hydraulic fracture preferentially propagates to the layer with higher formation permeability. These numerical results were verified by the study of Haimson [46], who has indicated that the increased pore pressure will generate additional stress and displacement in reservoir rock and decrease the critical injection pressure required for fracture initiation. They also explained why the hydraulic fracture propagates preferentially to the layer with higher formation permeability.

Conclusions
In this study, a 2D fully coupled hydraulic fracturing simulation model is developed using the XFEM-based CZM to investigate the hydraulic fracture propagation geometry in a layered heterogeneous shale reservoir. A parametric analysis of the elastic modulus, Poisson's ratio, tensile strength, and formation permeability is carried out. The results in this study provide a reference for completing the design and optimization of hydraulic fracture treatments. The following conclusions can be obtained: (1) In the specimen with a soft middle layer bounded by stiff outer layers or a stiff middle layer bounded by soft outer layers, the stiffer or softer outer layers  14 Geofluids decrease the fracture propagation width and increase the fracture propagation length. However, in an asymmetric specimen (with different elastic modulus in each layer), hydraulic fracture shows an asymmetric propagation behavior and the hydraulic fracture propagation rate in the stiffer outer layer are higher than those in the softer outer layer (2) The outer layers with higher Poisson's ratio decrease the fracture propagation width of the middle layer (lower Poisson's ratio). However, Poisson's ratio contrast does not significantly affect the fracture propagation length and breakdown pressure under the same injection time. For the asymmetric specimen with different Poisson's ratios in each layer, hydraulic fracture propagation length is longer in the layer with higher Poisson's ratio (3) When hydraulic fracture propagates from the middle layer with low tensile strength to the outer layer with high tensile strength, hydraulic fracture geometry presents the feature of "shorter and wider" as the tensile strength in the outer layer increases. The hydraulic fracture will preferentially propagate to the layer with lower tensile strength for an asymmetric specimen. In addition, hydraulic fractures tend to propagate preferentially to the layers with higher permeability

Data Availability
The data used to support the findings of this study are included within the article.