The Responses of Injection Pressure and Fracture Width during Height Extension in Sand-Mud Interbed Reservoirs in the Dongsheng Gas Field

The Dongsheng gas field is characterized by low porosity and low permeability. Its principal producing reservoir is the H-1 zone, belonging to the Lower Shihezi Formation. Sand is the main lithology in the H-1 zone, while mud interlayers are also well developed in a vertical direction. As a result, the reservoir is a sand-mud interbed, which brings difficulty to fracture height extension. In order to understand the process of fracture height growth in a sand-mud interbed reservoir and obtain the responses of injection pressure and fracture width during a hydraulic fracturing, a hydromechanical-coupled model is established. Mud interlayers are fully considered and a cohesive zone model is adopted to deal with fracture propagation. Numerical results show that the fracture extends quickly to the sand-mud interface after initiation and breaks through rather than propagating along the interfaces. Pressure and width both increase continuously when fracture propagates in the mud interlayer. High-viscosity and high-injection rates are helpful for the fracture to break through the mud interlayer, especially at an early period. When the mud interlayers are asymmetric, pressure and width fluctuate several times once fracture propagates inside and breaks through the mud interlayer. Perforations close to the thinner mud interlayer can increase the fracture width and reduce fracturing risks.


Introduction
The Dongsheng gas field belongs to the Erdos basin, which locates in North China. The principal producing reservoir of the Dongsheng gas field is Lower Shihezi Formation, especially the H-1 zone. Without too many complex structures, the formation is flat overall and has a good continuity horizontally. The thickness of the H-1 zone ranges from 10 m to 50 m and there are multiple gas formations in the vertical direction.
The main sedimentary system of the H-1 zone includes the alluvial fan and braided stream. As a result, the lithology is very complicated, including pebbly sandstone, coarse sandstone, medium sandstone, fine sandstone, siltstone, and dark mud. Particularly, the coarse sandstone, medium sandstone, and fine sandstone are the main reservoirs.
Due to the influence of river channel migration and the variation of the flood scale, barely permeable and impermeable interlayers are well developed vertically in the H-1 zone [1]. The most common interlayers in the H-1 zone are mud and property interlayers [2], as shown in Figure 1. The lithologies of mud interlayers have the characteristic of high clay content, such as pure mudstone, silty mudstone, and argillaceous siltstone, while the lithologies of the property interlayer are primarily siltstone or argillaceous fine siltstone. Compared with the mud interlayer, the property interlayer has higher porosity and permeability but they are still below the lower limit of the effective reservoirs. In this paper, only mud interlayers are considered to investigate its influence on fracture height growth. The thickness of the mud interlayer varies widely from several centimeters to 25 m. The thick mud interlayers can be treated as caps which can prevent fracture propagation.
The lithologies of gas reservoirs are also complex, including coarse sandstone, medium sandstone, and fine sandstone. The porosity ranges from 5.0% to 16.9% and the permeability is mainly in the range from 0.3 mD to 0.9 mD [4,5]. Both porosity and permeability are very low, and hydraulic fracturing is the most effective way to develop the Dongsheng gas field. Therefore, it is especially important to understand fracture height growth and specify the responses of injection pressure and fracture width during fracture propagation. In this paper, a hydromechanical coupled numerical model about fracture height extension is established based on the reservoir characteristics in the Dongsheng gas field taking into account of mud interlayers and stress difference. The numerical results illustrate the process of fracture height growth and are instructive to the design of a hydraulic fracturing.

Model
Classic linear elastic fracture mechanics (LEFM) is a widely used method to evaluate fracture propagation. When the stress intensity factor (SIF) at fracture tip is larger than fracture toughness which is an inherent property of the solid material, fracture initiates and propagates [6,7]. At the fracture tip, there is a stress singularity which is described by SIF in LEFM. But the SIF is usually difficult to obtain and this stress singularity bothers numerical simulations at times.
Instead of using the classic LEFM to deal with an elastic crack tip region where stress singularity exists, the cohesive zone model [8][9][10][11] assumes the existence of a fracture process zone (where the rock has yielded or experienced microcracking) in front of the material crack tip, which is governed by a traction-separation law. The stress singularity at the crack tip is avoided in cohesive zone models through this constitutive law. In this way, the cohesive zone model provides an alternative approach to explicitly simulate the fracture process near the tip and is often applied in modeling hydraulic fracturing.

Cohesive Zone Model.
In this work, we use a hydromechanical model utilizing the finite element method with a special zero-thickness interface element based on the cohesive zone model which has shown its effectiveness in previous work [12][13][14][15].
Indeed, there are numerous kinds of T-S models [16,17] and the classic bilinear traction-separation criteria [18][19][20] are used in this paper, as shown in Figure 2. These criteria evaluate the stress ratios between a given stress value and the peak nominal stress value in each of the three directions and is based on a quadratic combination of all three ratios [21], as shown in equation (1) as follows:

Geofluids
where t n , t s1 , t s2 are nominal stress, first direction shear stress, and second direction shear stress. t 0 n , t 0 s1 , t 0 s2 represent the peak values of the nominal stress when the deformation is either purely normal to the interface or purely in the first or the second shear direction. hi is the Macaulay bracket with the usual interpretation ht n i = t n if t n ≥ 0 and ht n i = 0 if t n < 0.
The fracture initiation and propagation are natural outcomes of a cohesive zone model. After cohesive element debonds, fluid flows in the hydromechanical interface and transmit pressure on the cohesive element, as shown in Figure 3.

Fluid Flow.
We assume that the rock surrounding the hydraulic fractures is permeable. Leakoff and poroelastic effects are considered during fracture propagation. Therefore, the fluid flow can be divided into three processes: flow inside the fracture, leakage through fracture, and seepage in the matrix.
As to incompressible fluid, flow in the hydromechanical cohesive interface is derived from the fluid mass balance equation: where w is normal opening, which equals to the displacement jump in the normal direction of the interface, q is the longitudinal fluid flux, q L is fluid leakoff, QðtÞ is the source term, and δðx, yÞ is the Dirac delta function. The lubrication equation for fluid flow between parallel plates is used to simulate flow in the fracture. It is given as follows: where v is the average fluid velocity over the cross-section of the hydraulic fracture, μ is the viscosity of the Newtonian fluid, and p is fluid pressure inside the fracture. Combining equations (2) and (3), we obtain the equation for fluid flow in a hydraulic fracture: For the seepage in the matrix, according to the theory of poroelasticity, the continuity equation utilized in numerical formulation is as follows [22]: where k is matrix permeability coefficient, γ l is specific weight of liquid, p 0 is the pore pressure, b l is the body force  In addition, fluid leakage through the fracture wall can be expressed by the following: where C L is leakoff coefficient.

Rock Deformation.
The governing equations of rock deformation include the equilibrium equation, geometric equation, and constitutive equation [23,24]. The equilibrium equation is written as follows: where σ is stress tensor and b is body force vector of the rock. According to Boit's theory, effective stress can be expressed as follows: where σ′ is the effective stress tensor, α is the Biot coefficient, and I is the Kronecker delta.
The geometric equation is written as follows: The constitutive relationship can be expressed as where λ and G are the Lamé parameters of the porous material, while C and M are additional elastic moduli required to describe a two-phase medium. ε vol is the volumetric strain. ζ is a strain parameter describing the volumetric deformation of the fluid relative to that of the solid. δ is the Kronecker function.

Model Construction.
In this study, a sand-mud interbed model is established. In Figure 4 Geofluids interfaces in order to investigate fracture height growth behavior at sand-mud interfaces. Rock mechanical parameters adopted in numerical simulations are shown in Table 1, and they are all obtained from experimental measurements using downhole cores from the Dongsheng gas field. In order to investigate the influence of rock mechanics on fracture height growth, both fine sand and coarse sand are tested. Table 2 includes basic parameters except for rock mechanics. The reservoir parameters are all used according to the Dongsheng gas field so that the numerical models are consistent with the actual reservoir conditions as much as possible.

Results
In order to clearly illustrate about fracture height growth in sand, we first study the cases neglecting mud interlayers. Afterwards, we consider two symmetric mud interlayers whose thicknesses are both 3 m. In the end, two asymmetric mud interlayers as shown in Figure 4 are taken into account which is in accordance with the characteristics of the Dongsheng gas field. Figure 5 shows fracture height extension in sand. At the sand-mud interface, the fracture passes through directly and extends to the cap. This is because the reservoir shear strength is much larger than tensile strength and tensile failure dominates. Besides, when fluid wants to flow into the interface, it needs to overcome overburden stress or vertical stress, which is very difficult because the reservoir is deeper than 3000 m. But at the interface, it still can be observed that stress discontinuity exists. Figure 6 clearly shows the fracture width and pressure responses during propagation. When fracture propagates in sand, pressure drops rapidly after initiation and keeps being smooth, while the width increases slowly. When the fracture reaches the interface and extends in the cap, propagation slows down. As a result, pressure increases and width     Figure 10: The influence of minimum horizontal stress difference on injection pressure and fracture width. 6 Geofluids expands. In reality, when pressure is large enough, fracture will propagate horizontally in sand rather than keep the extension in the cap. The zigzags shown in Figure 6 are due to the mesh size used in the numerical model. With smaller meshes, the lines will be much more smooth but the calculation will be tremendous. Therefore, the element size along fracture is 0.25 m in this model and calculation accuracy is not affected.

Reservoirs without Mud Interlayers.
3.1.1. The Effect of Sand Property. Figures 7 and 8 show that the fracture propagates in fine and coarse sands, respectively. With a lower tensile strength, the fracture propagates much easier in coarse sand. As a result, pressure is lower and width is smaller. However, when the fracture extends to the interface, pressure in coarse sand increases rapidly and surpasses that in fine sand. This is because the elastic modulus of coarse sand is lower. Under the same condition, fracture width is smaller and the friction along fracture is higher. If the fracture continuously propagates in the cap, it requires higher pressure to overcome friction.

The Effect of Minimum Horizontal Stress Difference between Cap and Sand.
Before the fracture reaches the interface, the responses of pressure and width are exactly the same. But after fracture extends into cap with higher stress, pressure increases higher and fracture opens wider, as shown in Figures 9 and 10.

The Effects of Fluid Viscosity and Injection Rate.
We consider four cases to investigate the influence of fluid viscosity and injection rate on height growth, and they are the low viscosity and low injection rate (LVLR), low viscosity and high injection rate (LVHR), high viscosity and low injection rate (HVLR), and high viscosity and high injection rate (HVHR). Compared with viscosity, it seems that the injection rate has the most influence on fracture width.   7 Geofluids Figure 11 clearly shows that with the high injection rate, fracture width is larger. At the same time, when the fracture extends to the interface, pressure increases rapidly. With higher viscosity, fluid friction along the fracture is higher, as shown in Figures 12(c) and 12(d).

Reservoirs within Symmetric Mud Interlayers.
In this section, two symmetric mud interlayers are inserted in sand. Figure 13(a) shows that fracture extends to the first interface which is corresponding to ① in Figure 14 when fracture pressure and width are about to increase rapidly. The pressure keeps increasing until ② when the fracture tip extends into the mud interlayer. While the fracture propagates in the mud interlayer, it seems that pressure no longer increases. Due to the high in situ stress of the mud interlayer, fracture propagation velocity slows down at the moment and friction along the fracture reduces. Therefore, pressure inside the fracture is much more uniform and pressure at the injection point tends to be similar to the pressure at the tip of fracture. As a result, fracture width keeps growing until ③, when the fracture breaks through the second interface and exits the mud interlayer. Afterwards, fracture extends to the interface of sand and cap. But it is a very short period from ③ to ④. This is because when the fracture propagates inside the mud, it requires a higher pressure. When the fracture leaves the mud and goes into the sand, the high   Figure 14: The responses of injection pressure and fracture width during height growth in the reservoir with symmetric mud interlayers. 8 Geofluids pressure will accelerate fracture growth. Therefore, the fracture reaches the sand and cap interface quickly as shown in Figure 13(c). In the end, the fracture extends inside the cap and the pressure increases integrally although the pressure at the injection point does not rise up obviously, which results in fracture width increasing continuously.

Effect of Fluid Viscosity and Injection
Rate. When the mud interlayers are imbedded in sand, the treatment parameters affect height growth differently. With a high injection rate, the fracture propagates to the first interface and breaks through the mud interlayer quickly. As shown in Figure 15, it seems that under the low fluid viscosity condition, the pressure required to break through mud interlayer is a little bit smaller. This may be resulted from the lower friction when fluid viscosity is low. Under the high viscosity and high injection rate conditions, pressure inside the fracture is the highest. Therefore, when we hope that the fracture breaks through mud interlayers in sand, high viscosity and high injection rate are the best choices especially at the early period.

Effect of Stress Difference between Sand and Mud
Interlayers. As shown in Figure 16, before the fracture reaches the first interface, the responses of pressure and width are exactly the same. But when the stress difference between sand and mud interlayers increases, the pressure required to break through the mud interlayer is higher and the fracture opens wider. Besides, it takes more time to break through the mud interlayers. As shown in Figure 17, when the stress difference increases from 2 MPa to 7 MPa, the net pressure required to break through the mud interlayer rises up linearly from 8.59 MPa to 12.42 MPa. The increment of net pressure,   Figure 17: The relationship between required net pressure and stress difference. 9 Geofluids 3.83 MPa, is smaller than the increment of stress difference, 5 MPa. This is because when the minimum horizontal stress is large, it becomes much more difficult for fracture to break through. Fracture propagation slows down and fluid friction along the fracture drops. As a result, pressure inside the fracture trends to be uniform and only a small increment of net pressure can keep the fracture propagating.

Reservoirs within Asymmetric Mud
Interlayers. Asymmetric mud interlayers are more common in reality and are considered in this section. The numerical model is the same as shown in Figure 4. The distances between the injection point to the upper and lower first interface are the same.
It can be observed that fracture extension becomes much more complex from that in Figure 18 and the process is divided into several stages, as shown in Figure 19. The upper and lower fracture branches reach ① and breaks through the first interface ② at the same time. But the stress distribution is no longer in symmetry from the reservoir top to bottom due to the influence of asymmetric mud interlayers. Afterwards, both upper and lower branches extend in the mud interlayer. But the lower branch quickly reaches the second   Figure 19: The responses of injection pressure and fracture width during height growth in the reservoir with asymmetric mud interlayers. 10 Geofluids interface and the fracture width keeps growing ③. Right after the lower branch breaks through the lower thin mud interlayer, it goes quickly to the interface between sand and lower cap ④. Figure 18(c) shows at the moment when the lower branch reaches the cap interface, the upper branch is still inside the upper mud interlayer and does not break through. As both branches extend in mud, the pressure increases continuously. Finally, the upper branch breaks through the second interface and the width decreases immediately once more ⑤. As the number of mud asymmetric interlayers increases, the fluctuation of pressure and width will be more frequent. Every time the fracture breaks through the mud interlayer, it requires a higher pressure. Therefore, in reality, it may be difficult for the fracture to run through all interlayers. But at the early period, fluid with high viscosity and high rate is helpful for fracturing all sand reservoirs. Once the height extends into the multiple sand layers, it will become much easier for the fracture to propagate inside the sand. In order to understand the influence of the perforation position on fracture height growth when the mud interlayers are asymmetric, two cases are studied. The distance between perforation to the center of the model in both Figures 20(a) and 20(b) is 1 m. In Figure 20(a), perforation is above the center of model. In Figure 20(b), perforation is below the center of the model. It can be observed that after both branches extend into the cap, these two cases have no difference from the perspective of stress filed. But during the propagation process, it can be observed in Figure 21 that when the perforation is below the center point, which means the position of perforation is close to the thinner mud interlayer, fracture width is much larger during fracture propagation. This is because the upper thick mud layer is difficult to break through and the fracture is asymmetric. That is to say, perforating close to the thinner mud layer can enlarge fracture width and make sand transport inside the fracture smoothly, which can reduce the risks of sand plugging.

Conclusions
This paper focuses on fracture height growth during hydraulic fracturing in the Dongsheng gas field. A hydromechanical sand-mud interbed model is established taking into full consideration reservoir characteristics. Besides, the parameters used in the simulations are from the experimental tests on downhole cores. Based on the numerical simulations, it can be concluded that fluid is difficult to flow inside the interfaces between sand and mud as it requires a much higher pressure to overcome overburden stress. As a result, the fracture extends vertically and breaks through mud interlayers without orientation. With a larger Young's modulus, fracture width in coarse sand is smaller than that in fine sand.
When there are mud interlayers within sand, the process of fracture height growth is complicated especially when the mud interlayers are asymmetric. Injection pressure and fracture width vibrate as fracture extends to and breaks through  11 Geofluids the mud interface. As the minimum horizontal stress difference between the mud interlayer and the sand layer increases, the net pressure required to break through the mud interlayer increases linearly. What is more, when the mud interlayers are asymmetric, the fracture branch near the thinner mud interlayer breaks through the interlayer and reaches cap interface first. Afterwards, the other side branch propagates to the cap interface as the pressure inside the fracture increases. What needs to be stressed is that what we studied in this paper is a 2D height growth problem. In a hydraulic fracturing, when fracture height extends to the mud interlayers, it may propagate continuously in sand rather than break through the mud interlayer. This is because the stress of mud is larger than that of sand under normal circumstance. Fracture propagates much easier horizontally. If we want the fracture height to break through multiple mud interlayers, a high net pressure inside the fracture is required. Therefore, according to the simulation results in this paper, high viscosity fluid and high injection rate are beneficial to build a high net pressure inside the fracture, especially at the early period. In addition, when the mud interlayers are asymmetric, perforation close to the thinner mud interlayer can enlarge fracture width and make sand transport inside the fracture smoothly, which can reduce the risk of sand plug.

Data Availability
The data presented in this paper are available upon request from the corresponding author.

Conflicts of Interest
The authors declare no conflict of interest.