Numerical Simulation of Hydraulic Fracture Propagation on Multilayered Formation Using Limited Entry Fracturing Technique

: Hydraulic fracturing is one of the most effective stimulation methods for unconsolidated sandstone reservoirs. However, the design of hydraulic fracturing must take into account the mechanical and stress properties of different geological formations between layers. In this paper, a three-dimensional coupled fluid-solid model using the finite element method is developed to investigate multiple vertical fractures at different depths along a vertical wellbore under different geological and ge-omechanical conditions. The finite element model does not require further refinement of any new cracks, requiring much smaller degrees of freedom and higher computational efficiency. In addition, new elements were used to account for local pressure drop due to perforation entry friction along the vertical wellbore. Numerical simulation results indicate that hydraulic fracture connections are observed from adjacent layers. Furthermore, the low stress contrast and high Young’s modulus between the layers increases the likelihood of multiple fracture connections. Higher fluid leakage rates increase the likelihood of fracture branching, but decrease the area of fracture coverage near the wellbore. Increasing ﬂuid viscosity is eﬀective in improving the area of fracture coverage near the wellbore. These ﬁndings are useful for the design of hydraulic fracturing in multi -layered formations in un-consolidated sandstone formations.


Introduction
Multilayered formations are found rich in both offshore and onshore reservoirs.Multilayer reservoirs have been thought to benefit much from hydraulic fracturing stimulation.However, the development of multilayer reservoirs remains a major challenge in the oil and gas industry due to the presence of multiple thin vertical pay zones rather than successive oil and gas accumulation.In order to maintain uniform hydraulic fracture propagation in multilayered formations, separate hydraulic fracturing treatments can be performed in each layer with the help of packers, but the high cost and long operation time limit individual hydraulic fracturing treatments for each layer.Therefore, a limited fracturing technique has been proposed by adjusting the drilling friction.Compared to individual fracturing, limited-entry fracturing provides an opportunity for more uniform hydraulic fracture propagation.However, the design of fracturing with limited entry in multilayered formations must take into account the characteristics of the multilayered formations and the in situ stress contrast.Figure 1 shows the stress distribution in multilayered formations in the depth direction.Therefore, a better understanding of the fracturing mechanism of hydraulic fracturing with limited entry in multilayered formations is needed [1,2].
Citation: Liu, H.; Ji, W.; Huang, Y.; Zhang, W.; Yang, J.; Xu, J.; Mei, M.  To better understand the mechanisms of hydraulic fracturing in multilayered formations, several authors once studied this problem using numerical and experimental methods [3][4][5].Athavale and Miskimins, 2008 performed triaxial fracturing experiments on laminated cement blocks and found that interlayer rock characteristics and stress differences lead to complex fracture propagation [6].Goyal et al., 2017 presented a case study of limited-entry technology in multistage fracturing in tight gas volcanic reservoirs in western India.They found that more 10 stages are needed to effectively cover the available net pay [7].Gao et al., 2017 built a 3D finite element model to simulate the propagation of multiple hydraulic fractures in a wellbore to a tight reservoir.In addition, they analyzed and compared the effects of different fracturing techniques and formation properties on fracture morphology [8].Li et al. 2019 used the numerical simulator RFPA-Petrol (Rock-Failure-Process-Analysis on Petroleum-problems) linked hydraulic-mechanical-damage (HMD) model to investigate fracture paths.The results indicate that under given geological conditions, hydraulic fracturing cannot cross barriers between adjacent layers.Previous studies have found that the presence of layers makes the fracture geometry more complex and thus affects the fracturing effect [9].However, many studies fail to take into account the effects of limited press-in techniques.

Numerical Simulation of Hydraulic
Limited entry technique has proved to be a cost effective method of increasing net pay coverage and EUR per well with minimum number of frac stages.While using the Limitedentry method, the diameter or quantity of perforations is reduced to enhance the fluid friction through the perforations.Thus, the difference in fracture propagation resistance can be balanced.However, the design of Limited-entry technique is based on the try and error method, and the success rate of Limited-entry technique cannot be guaranteed.In recent years, many researches focus on the multiple fracture propagation in the horizontal well.Prada et al., 2012 presented the frac-pack method for the multilayered formation with low stress contrast.In addition, they suggested a perforation interval to maximize reservoir exposure without compromising the frac pack quality [10].Yang, 2012 applied the pseudo three dimensional models to enhance the well production for multilayered formation.In addition, she also applied the unified fracture design for the fracture geometry optimization [11].Li et al., 2017 simulated multiple fracture propagation along the horizontal well, they demonstrated that a sufficiently large perforation pressure drop can counteract the stress shadowing among the multiple fractures and can help to equalize the length of multiple fractures [12].However, the fracture mechanisms under the limited-entry technique for vertical well in multilayered formation is relatively rare, although limited entry fracturing is a great stimulation technique for vertical wells with a manageable degree of cost and risk [13][14][15].It is worth noting that the fracture toughness parameter is an important parameter for hydraulic fracturing design.Li et al., 2023 developed a planar three-dimensional model to study the influence of rock fracture toughness heterogeneity on hydraulic fracture propagation [16].Zvyagin et al., 2021 proposed a numerical method for solving spatial problems in fracture mechanics, which performed well in both quantitative and qualitative display [17].Pestov et al., 2023 conducted a comparative study on cracks in one matrix and two matrices, and defined the influence coefficient as the research object [18].Shamina et al., 2024 proposed a study on a single polygonal crack under a constant load perpendicular to the crack plane.The boundary element method was used to numerically construct the solution [19].As for the multiple fracture propagation vertically, it's necessary to consider the combining effects of in-situ stress, mechanical and perforation friction difference.For the three-dimensional model of hydraulic fractures in heterogeneous fractured reservoirs, some scholars have started to study it.Li et al., 2020 andKiselev et al., 2021, etc. obtained numerical asymptotic solutions for the initial condition problem (early solution) of plane strain hydraulic fractures with fluid hysteresis.And complete the planar three-dimensional simulation of fluid flow [20,21].
The numerical simulation of hydraulic fracturing is hot topic in oil and gas industry.There are now three different types of numerical simulation methods: discrete element method (DEM), boundary element method (BEM), and finite element method (FEM).Because only fracture faces are discretized, BEM offers the benefit of an efficient engineering community, albeit its use is typically limited to linearly elastic materials.DEM is more costly to compute and takes longer to calibrate material constants than FEM and BEM, but it can handle intricate fracture patterns and provide precise physics about the interactions between rock particles and fluids during hydraulic fracturing.Comparing to traditional hydraulic fracturing methods, CZM simulations benefit from relatively short simulation time and the ability to capture the mechanical behavior of rock.In addition, stable fracture growth in which both fracture energy and maximum displacement for de-bonding are constant.CZM provides useful force-displacement relationships for fracture mechanics studies when the shear and tensile stress of the interface zones exceed the strength of rock.Several authors once investigated the hydraulic fracturing using CZM.Studies on the behavior of interactions between multiple vertical fractures at varying depths along vertical wells are common; however, it is uncommon to find studies on the 3D hydraulic fracture propagation for multilayered unconsolidated sandstone reservoirs.
The motivation of this paper is to present a 3D model to simulate the crack propagation problem in hydraulic fracturing with degradable diverting materials via the CZMbased finite element method (FEM).Spring elements are adopted to bridge off the opened HFs, thus better representing net pressure increases in sealed fractures.Cohesive elements can be automatically generated anywhere in the model using Python scripting.This adaptive model is useful for simulating different potential crack propagation paths in TPDF.In this paper, a 3D multilayer hydraulic fracturing model using finite element method is built to investigate the multiple vertical fractures at different depths along the vertical wells considering different formation mechanical and stress properties.In addition, a new element to account for the local pressure drop due to perforation entry friction, which can consider the fluid distribution using limit entry technique.The numerical simulations are beneficial for hydraulic fracturing design in multilayered formation.

Mathematical Methods
The governing equations for the multiple fracture propagation along the multiple perforations includes the fluid flow, rock deformation and pressure drop equations.

Fluid Flow and Rock Deformation
The fluid continuity equation can be described as: where J denotes the rate of volume change of porous media, dimensionless; ρw is the fluid density, kg/m 3 ; nw is the pore ratio, dimensionless; vw is the fluid percolation velocity, m/s; and X is the space vector, m.
The rock deformation equation can be characterized by Biot's law of effective stress.In the current setup, the virtual work principle is defined as follows [22]: where ′ and ′ are effective stress and virtual rate of deformation, respectively.I is the unit matrix, and t and f are surface traction per unit area and body force per unit volume, respectively.The nodal variables in this equation are displacements, which are discretized using a Lagrangian formulation.

Fluid Flow within Hydraulic Fracture
There are two parts for the fluid flow in the hydraulic fracture.One is the longitudinal flow of fluid along the fracture another is the normal flow of fluid from the fracture surface to the surrounding porous medium.The tangential flow rate q within the fracture can be written as: where qf is the flow rate along the fracture length, µ is the viscosity of the fluid in the fracture, and w is the fracture opening, and pf is the pressure in the fracture.The continuity equation of mass conservation can be expressed as: where q refers to the fluid flux of the tangential flow for each element, p represents fluid pressure inside the fracture, w is the hydraulic fracture width, µ is the viscosity of used hydraulic fracturing-fluid.
To explain the typical flow of hydraulic fracturing, the Carter leak-off model is also presented.The fully coupled poro elastic analysis is used in the proposed model to determine the leak-off of the fracturing fluid through the fracture faces into the surrounding medium.The regular flow in the fracture is written as follows: where qt and qb are the typical flow flux cross into the cohesive element's top and bottom surfaces, respectively.Pt and Pb are the pore pressures on the top and bottom surfaces, respectively, while Pi is the fluid pressure within the cohesive element gap.Ct and Cb are fluid leak-off parameters.

Fluid Flow within Multiple Perforations
It is possible to treat the fluid flow in the wellbore and through the perforations as a single-phase, one-dimensional pipe flow.
To build the integrated wellbore perforation flow model, denoted by Equation (6), Kirchhoff's law is introduced (1).n = , 1, , where QT is the total flow rate in the wellbore, m 3 /s; Qi is the flowrate of the ith fracture, m 3 /s; Pout,i is the pore pressure of the outflow node in the ith fracture, Pa; Pwellbore is the pore pressure in the wellbore, Pa; Ppf,i is the perforation friction of the ith fracture, Pa.
The limited entry approach is a fracturing method that promotes the formation of perforation friction pressure in multilayered reservoir stimulation treatment by limiting the number or size of perforation holes in a completion interval.This makes it possible to direct the fracturing fluid through the selected perforations into the formation of interest.The diameter and quantity of perforations affect the flow resistance.The perforation holes serve as chokes connecting the wellbore and the formation.Using pipe fluid elements (FPE) with two pore pressure nodes, fluid transport between multi-cluster fractures and perforation friction were both modelled.In this investigation, pressure node will receive the fluid as it enters, and the other will receive it as it exits.The link between perforation friction and flowrate during fluid flow through the perforation cluster is characterized as follows [23,24]: where ith is the number of fracture clusters, dimensionless; Δpifric is the perforation friction through the ith fracture, Pa; Qi is the fluid flow rate of the ith fracture, m 3 /s; ρ is the fluid density, kg/m 3 ; Np,i is the perforation number of the ith fracture before diversion, dimensionless; Dp,i is the perforation diameter of the ith fracture, m; C is a dimensionless discharge coefficient ranging between 0.60 (before erosion) and 0.90 (after erosion).

Cohesive Element Method (CZM)
An effective approach to the consideration of hydraulic fracture propagation is the cohesive zone model: the stress specificity at the fracture tip present in the (Linear elastic fracture mechanics) LEFM is avoided by assuming the existence of a fracture process zone at the fracture tip characterized by a traction separation constitutive law.The cohesive zone approach benefits from the simplicity with which it may be integrated into a conventional finite element framework.Figure 2 shows the schematic of hydraulic fracture propagation with the CZM approach.The cohesive element's damage initiation and evolution are described, respectively, by the quadratic nominal stress law and linear damage evolution law based on effective displacement.When a nominal stress ratio-based quadratic interaction function approaches the value of unity, damage is assumed to have started.The standard might be stated as where 0 s t and 0 n t stand for the normal and shear strengths of an undamaged cohesive element, respectively; the symbol stands for the Macaulay bracket; and tn and ts stand for the normal and shear stress components, respectively.
The normal and shear stress components can be defined as [25,26]: ( ) The elastic separation-traction behavior predicts Tn and Ts as the normal and shear stress components, respectively, for current displacement without damage.The complete fracturing procedure is described by the scalar damage variable D. Damage evolution, based on mix-mode energy theory, establishes the proportions of normal and shear deformation.The area under the traction separation curve is equivalent to the fracture energy.
D varies with the fracturing operation time and increases linearly from 0 to 1 after the damage begins.In this model, the overall damage variable D is determined by the traction separation law and calculated as follows: Among then, δ is a constant related to material proper- ties that can be measured experimentally.m δ is the effective displacement, defined as: Among then, n δ , s δ and t δ represent the normal, first, and second shear displace- ment components of the cohesive elements, respectively.
f m δ is defined as: ( ) G and t G respectively represent the work done by the traction force and its conjugate relative displacement in the normal, first shear direction, and second shear direction.These parameters can be obtained using traction separation curves.

Discretization of Coupling Equations
The weak form of the fluid seepage stress equation is: Among then, Among then, 0 w ρ represents the fluid density in the reference configuration.Use Newton's iterative method to solve Equations (18) and (19).

Numerical Model Setup
A 2-D CZM model combined with user-defined perforation elements was developed to investigate the complex fracture mechanisms of hydraulic fracturing in multi-layered formations.The same simulation steps were applied for each case considered.In general, there are four necessary steps.The geostatic equilibrium between the initial loads (stresses and pore pressure) and the boundary conditions should be set, which is called Geostatic step.In addition, the initial loads can be applied programming subroutines, which are match with the actual in-situ stress environment.In the second step, defined as consolidation analysis, a non-Newtonian fluid is injected at a given flow rate to investigate the influence of fracture fluid viscosity on hydraulic fracture growth.In the third steps, a temporary plugging agent mixed with the proppant is injected.At this stage, due to the filter cake, fluid pressure is not allowed to be transmitted to the tip of the crack, so it reduces fluid pressure transmission within each crack and increases pressure in the wellbore.In the fourth step, inject fracturing fluid again and observe crack propagation.No confining stresses and reservoir pressure are applied to the formation, and the outer normal displacement and pore pressure are fixed to zero during the entire process.Figure 3 shows the schematic of hydraulic fracture propagation for multilayer zones at the unconsolidated formation.The size of the model is 55 m × 100 m.The target formation is composed of five zones, and two pay zones are set, which are upper formation and lower formation.This study used a nonlinear FEM equation solver ABAQUS to solve the fracture initiation and propagation equations.The C3D8P element, which has three displacement nodes and eight seepage nodes, incorporates the equations for matrix seepage and rock deformation.The cohesive zone element COH3D4P, which has a pore pressure degree of freedom, is coupled to the fracture initiation and propagation equations to simulate the fracture propagation.The three-dimensional fluid pipe elements with two pore-pressure nodes FPC3D2 have the fluid distribution equations between many fractures coded into them.It should be noted that C3D8P: A three-dimensional porous elastic element used for the analysis of permeability and structural interactions.The FPC3D2 unit embeds a main unit with both translational and pore pressure degrees of freedom, and only the pore pressure degrees of freedom shared by the embedded unit nodes will be constrained.COH3D4P is a cohesive unit.The finite element framework may perform automatic computations by combining all these components into a single integrated model.In fact, this simulation is based on the treatment of ultra static water pressure, and the pore pressure boundary is defined as 0. The initial stress and pore stress of the second step are automatically balanced by the model after the geostress is balanced, without any additional settings.The input information for the simulation of hydraulic fracturing in multilayer zones for the unconsolidated formation is provided in Table 1.The input data for the numerical model are all from on-site actual relevant data and indoor experimental data [27].As shown in Figure 4, the distribution diagram of the cohesive force unit nodes in the model.The model size is 120 × 80 m.

The Influence of Young's Modulus on the Fracture Geometry
In order to quantify the effect of Young's modulus for the pay zone in different depths of limited entry fracturing, we compared the geometry in multilayered fracturing in Figure 5.We study the effect of Young's modulus on crack propagation by changing the magnitude of the Young's modulus in the production layer in Figure 6.The color contours represent the fracture width for each figure.It shows that the increase of Young's modulus can increase the fracture length for each pay zone.The fracture length in both two pay zones increase to far extent.That's because the increase of the Young's modulus increases the formation brittleness.The increase in rock brittleness index increases the fracture propagation capability.It can be seen that the fracture length increases more rapidly with the Young's modulus of the rock before 14 GPa, and the fracture length of the upper reservoir increases nearly 15 m between 12 GPa and 14 GPa.In the contrary, from the simulation results of the fracture width under different conditions, it can be obtained that the fracture width decreases more before 14 GPa.For example, in the upper reservoir, the fracture width decreases from about 5.54 mm to 4.53 mm.The reason for this is that the increase in elastic modulus makes the reservoir denser and more brittle, leading to the rapid expansion of the main fracturing fracture once it forms, resulting in a decrease in the expansion in the width direction.

The Influence of Fluid Leak-Off Rate on the Fracture Geometry
In order to quantify the effect of flow leak-off rate in different depths of limited entry fracturing, we compared the geometry in multilayered fracturing in Figure 7.The color contours represent the fracture width for each figure.From the hydraulic fracturing simulation results, it can be seen that when the formation fluid leak-off rate increases, the fracture width and fracture length are undergoing a significant decrease.In the upper reservoir, when the fluid leak-off rate ranges from 10 −6 m/s 0.5 to 10 −4 m/s 0.5 , the width of the fracture decreases abruptly from about 5.55 mm to about 4 mm, and when the fluid leak-off rate is larger than 10 −4 m/s 0.5 , the fracture width does not decrease significantly.Fracture length in 10 −4 m/s 0.5 at the same rate of decline rate inflection point, but the fracture length will continue to decrease because of the increase in the fluid leak-off rate.From Figure 8, it can be seen that the increase of the fluid leak-off rate will cause a large amount of fracturing fluid to be lost into the formation, and the fracturing fluid will be present in a large amount in the near-well formation around the wellbore, and the energy of fracture initiation and expansion can not be maintained, which not only seriously affects the reservoir stimulation effect, but also possibly induce formation damage.The reason for this is that the greater the fluid leakage rate, the more severe the fracturing fluid loss, ultimately leading to a decrease in both crack length and width, which is not conducive to fracturing operations.

The Influence of Fracture Viscosity on the Fracture Geometry
The changes of fracture morphology under fracturing fluid viscosities of 10 mPa•s, 50 mPa•s, 100 mPa•s, and 150 mPa•s were compared in Figure 9, respectively.The color contours represent the fracture width for each figure.Numerical simulation results show that fracture length and width in unconsolidated sandstone reservoirs decrease and increase, respectively, with increasing fracture fluid viscosity.The increase in fracture fluid viscosity is mainly due to the increase in resistance to fluid flow in the formation, which significantly reduces the rate of fluid leakage around hydraulic fracturing.Moreover, the increase in the fracturing fluid viscosity can also increase the fracturing energy, thus more energy can be consumed for the fracture open.

The Influence of In-Situ Stress Difference between Layers on the Fracture Geometry
The direction of the current maximum horizontal stress affects the direction of hydraulic fracturing and influences fracture geometry.The color contours represent the fracture width, see Figure 11

Discussions
Comparing to single layer formation, it is well-known that the multilayered formations can exhibit an in-situ tectonic stress and mechanical properties distribution between pay layer and boundary layer.In order to study the effect of in-situ stress and mechanical properties on fracture propagation in multilayer formation, several examples are studied.Because the power of the perforation charge on separate layer, two hydraulic fractures can be successfully initiated.However, the detail observation of two fractures in different layers, there are still some differences in the fracture geometry.In general, the fracture stimulated zones in the upper zone and lower zone are different, and larger fracture stimulated zone can be observed in the upper zone.One important reason is the injection perforation interval in the upper zone can receive larger fluid injection volume than that in the lower zone.Compared to the fracturing simulation injected into the casing element, simulating the injection interface element again adds some initial damage, making the model easier to converge.Moreover, the flow frictional resistance can also reduce the fracturing energy from different perforations, thus reducing the fracture propagation length in the lower formation of multiple layers.
Although there have been many successful cases of in-situ fracturing treatment for multilayer reservoirs, the basic mechanism of fracture propagation is still not fully understood.The geometry of multiple fractures can be significantly affected by different fracture schemes (geomechanics and different fracture treatment parameters.The numerical results show that the Young's modulus, leak-off coefficient, stress contrast between layers and fracturing fluid viscosity have a negligible impact on the concurrent growth patterns of multiple hydraulic fractures.Due to limit-entry fracturing and interlayer stress shadows between fractures, fractures do not form obvious expansion modes from convergence to divergence.Thick barriers can also increase the fracture propagation lengths and stop fractures from communicating in the direction of height.The total fracture height in the stage dramatically drops when the layer's modulus lowers, whereas the total fracture length generally increases, according to quantitative study of the layered modulus effect.Additionally, it appears that the increase the layer stress difference is beneficial for fracture containment in the target layer.One important is the stress limit the fracture height increase.Considering the high permeability of unconsolidated sandstone, the influence of the fracturing fluid leakoff rate cannot be neglected.Branch fractures can be observed near the main fracture in the multilayered formation stimulation.The fracturing fluid leak-off volume increase can result in larger wellbore coverage, while the total fracture volume near wellbore increases.As a result, change of the fracturing fluid viscosity is a useful tool in the fracturing fluid leak-off control, the increase of the fracturing fluid viscosity not only improve the fracture coverage near wellbore, but also result in fracture extent area.Therefore, the selection the high efficient fracturing fluid types is necessary for the stimulation for the multilayer unconsolidated formation In practical engineering applications, the Young's modulus and interlayer geostress of reservoirs belong to geological conditions and cannot be changed.For the viscosity of fracturing fluid, in practical field applications, high viscosity fracturing fluid is usually used to control the loss of fracturing fluid, which is more important in high permeability loose sandstone and also reduces the fluid leakage rate.However, the viscosity should not be too high, otherwise it is not conducive to the expansion of fractures.For perforation, increasing the number of perforations helps to open more cracks, but it also increases fluid leakage, so perforation design needs to be optimized.

Conclusions
Understanding the hydraulic fracturing mechanism in multilayer formation using the limited entry fracturing technique need consider specific multilayered properties.In this paper, a 3D fluid-solid damage coupling finite element model was established to investigate the fracture configurations under different fracturing treatment parameters and formation properties in multilayered unconsolidated sandstone reservoirs.From the numerical simulation results, we can get the following conclusions.
1.In vertical wells, interference between neighboring vertical fracture tips caused by neighboring perforation intervals, increasing the potential of multiple fractures communication process.In particular, the hydraulic fracture in the upper formation has a larger reservoir stimulated area because of relatively low flow resistance.

Limit entry fracturing reduces the likelihood of communication between multiple
fractures from different layers.Because of this, activating multilayered reservoirs by distinct perforation design on each layer is preferred with respect to the fluid distribution capability.3. Neighboring pay zones and barriers with strong in situ stress contrast, high fracturing fluid viscosity, and low fluid leak-off rate and elastic modulus contrast can prevent hydraulic fractures from communicating along the height direction and length.Therefore, limit entry fracturing method is preferred in the above formations with large possibilities of fracture communication.
Propagation on Multilayered Formation Using Limited Entry Fracturing Technique.

δ
represent the displacement when the traction force reaches max T and the fracture is completely damaged; 0 m

)
Among then, max T represents the tensile or shear strength, and is a constant calculated from the cohesive element material parameters; c G represents the mixed mode frac- ture energy.Once the cohesive elements fractures, damage can be evaluated based on fracture energy theory.The damage evolution of cohesive elements during fracture propagation is determined by the Benzeggagh-Kenane (B-K) fracture criterion, defined as:

Figure 2 .
Figure 2. Schematic of hydraulic fracture propagation with the CZM approach.
in pore fluid; J represents the volume change rate of porous media; v represents the seepage velocity after the equation becomes weak.Through the divergence theorem, Equation (18) becomes:

Figure 3 .
Figure 3. Schematic of hydraulic fracture propagation for multilayer zones at the unconsolidated formation.

Figure 4 .
Figure 4. Distribution diagram of cohesive force unit nodes in the model.

Figure 9 .
Figure 9. Influence of different fracturing fluid viscosity on fracture morphology.

Figure 10 .
Figure 10.Relationship between Fracturing fluid viscosity and fracture length and width under separate layers.
. The magnitude of the in situ stress can be determined by leak-off test, and the orientation of the maximum horizontal stress can be determined by resistivity image logging to determine the distribution and change in the depth and lateral directions.The change of crack morphology with stress difference is simulated for interlayer stress differences of 1 MPa, 2 MPa, 3 MPa, and 4 MPa, and the change of hydraulic fractures with in situ horizontal stress difference is shown in Figure12.The numerical simulation results show that the fracture length increases linearly with increasing interlaminar stress.When the stress difference increases from 1 MPa to 4 MPa, the fracture length increases from 35.3 m to 48.2 m.The crack width does not change with respect to the in situ stress difference.Since the crack initiation zone tends to originate from the low stress location, the increase in stress difference limits crack propagation in the target zone.

(Figure 11 .
Figure 11.Influence of different interlayer in-situ stress on fracture morphology.

Figure 12 .
Figure 12.Relationship between interlayer in-situ stress and fracture length and width under separate layers.

Table 1 .
Input for the hydraulic fracturing simulation in multilayer zones for the unconsolidated formation.