Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

Hydraulic fracturing of horizontal well hydraulic fracturing technology can help develop unconventional geothermal and petroleum resources. Today, industry uses simultaneous and sequential (also known as zipper) fracturing in horizontal petroleum well stimulations. To achieve successful and desired stimulated rock volumes and fracture networks, one must understand the effect of various rock and fluid properties on stimulation to minimize the risk of unwanted fracture geometries. This paper describes the development of a 2D coupled displacement discontinuity numerical model for simulating fracture propagation in simulta‐ neous and sequential hydraulic fracture operations. The sequential fracturing model consid‐ ers different boundary conditions for the previously created fractures (constant pressure restricting the flow back between stages and proppant-filled). A series of examples are pre‐ sented to study the effect of fracture spacing to show the importance of spacing optimiza‐ tion. The results show the fracture path is not only affected by fracture spacing but also by the boundary conditions on the previously created fractures.


Introduction
Increased interest in exploration and production of low permeability reservoirs presents new challenges in design and evaluation of hydraulic stimulation treatment. Hydraulic fracturing of unconventional petroleum resources (oil and gas shales) relies on multiple transverse hydraulic fracturing of horizontal wells. Each treatment stage in a well is designed to generate a stimulated volume defined as the rock volume contacted by treatment fluid with a desired enhancement to permeability. The collective network of stimulations should affect the maximum volume with minimal overlap of adjacent treatment stages. Usually, HF treatment of horizontal wells is carried out using two schemes namely, Simul frac and Sequel Frac. In simultaneous fracturing multiple fractures are created and propagated at same time whereas in sequential fracturing, fractures are created one after another usually by keeping the previously created fracture propped [1] or pressurized with fluid [2]. In both cases, the perforation clusters should be placed appropriately to reduce stress-shadow effects. By reducing the number of clusters per stage, they were able to minimized stress interference, which reduced the possibility of having improper fracture propagation.
In order to determine the optimum spacing and optimum staging between fractures production forecasting analysis is used by assuming simple straight lined fractures, but in reality fractures may propagate in complex manner when they are closely spaced or are near preexisting fractures as they will interfere to repel or attract each other [3]. In simultaneous fracturing closely spaced fracture interferences causes some of the fractures to stop in between or some may not even initiate due to the stress shadow between them [4]. The design of efficient systems can benefit from hydraulic fracture simulations that couple fluid flow to fracture deformation and fracture mechanics principles. Since the fracturing itself is too complicated these problems are difficult to analyze using laboratory experiments. Numerical method that can accurately model 2D or 3D fracture propagation can help to understand and improve the fracturing process.
[5] solved the growth of multiple simultaneous fractures assuming no fluid-flow inside the fractures; [6] simulated the sequential fracturing with no explicit fluid flow and assuming the previously created fracture dimensions are constant. In [3] previously created fracture in sequential fracturing is assumed to be propped and having a shape elliptical fracture similar to the fracture geometry formed from uniform pressure distributed fracture. The fracture curving is attributed to opening and sliding of previously created fracture. [7] Reproduced to similar results by considering the previously created fracture is uniformly pressurized instead of propped while stating the reason behind the fracture curving is unclear. In this paper a fully coupled DD-based sequential fracturing model is presented which considers previously created fracture as uniformly pressurized and also propped. A linear joint model is used to model the propped fracture. This allows the propped fracture to open/close and shear as the next fracture propagates. This paper also includes the simulation of simultaneous propagation of multiple fractures spaced at different distances. The fracture curving observed in simulations are explained using the stress distribution plots around the fractures. The model can be used to study the effect of parameters such as differential stress, Young's modulus, Poisson's ratio, viscosity of the fluid on fracture propagation. The model calculates the flow rate and pressure within each fracture as they propagate with injection onto the wellbore. Currently, we use the model to analyze propagation of multiple hydraulic fractures to show the importance of spacing optimization.
In our simulations, we consider two different scenarios for sequential fracturing, one scenario is where the previously created fractures remain pressurized by restricting the flow back between stages [8] and the other is where the previously created fractures are filled with proppant [1].

Model development
The model developed for the research is based on 2D plane strain and uses the displacement discontinuity method (DDM) to calculate fracture deformation and propagation. The fluid flow inside the fracture network is governed by Lubrication equation [9]. The hydraulic fracture model fluid flow and fracture deformation through an iterative scheme between fracture aperture along the fracture length and fluid pressure. This is a non-linear problem that is solved using Newton-Raphson method. The fracture propagation scheme for sequential hydraulic fractures propagation employs an iterative scheme to find the pressure at fracture tip required to meet the propagation criterion. Joint deformation model is used to simulate the propped fractures by specifying the proppant properties in terms of stiffness. Finally, fracture propagation path is determined using the maximum tensile-stress criterion of [10]. Each of these model components are briefly described below.

Displacement discontinuity method
In this model the displacement discontinuity boundary element method is used to find fracture deformation. In implementing this method, a fracture is divided into n small elements and by specifying the normal and shear stress acting on each element, the resultant normal and shear stresses on each fracture element is found by using superposition [11]: shear and normal DD elements. The above system of linear equations can be solved for displacement discontinuity of each fracture element.
Using constant displacement discontinuity elements at the crack tips lead to inaccurate value of stress intensity factor [12]. In fracture mechanics it is very important to have an accurate value to stress intensity factor, as it decides the condition for propagation and crack paths. In order to calculate accurate displacement discontinuities at crack tips, this model incorporates a crack tip element [11] in which the relative normal displacement discontinuity between the crack surfaces is given by u y (x) = D y ( x a ) 1/2 where a is half length of the crack tip element, D y is the displacement discontinuity at the center of special element and x is the distance measured along the element from the tip of the crack. The influence coefficients and formulation for the special element used herein is given in [12].

Joint model
A joint model is useful to simulate propped fractures where one can model the width reduction of propped fractures (proppant closure) due to the creation of new fractures. In this model we assumed a propped fracture behaves like a joint (natural fracture). The proppant pack inside the fracture is assumed to be a compressible element and its displacements can be calculated using the DD method. The joint elements have normal and shear stiffness that represents the filling material characteristics. Though the joint filling material usually deforms non-linearly, here it is assumed to deform linear (linear model in Figure 1) with the stress for simplicity. Given the far field stresses (σ ij ) 0 ∞ and stresses acting on the joint element, the total joint deformation will be the sum of initial displacements (due to initial stresses on the joint) and induced displacements (due to induced stresses caused by the fracturing in the formation) can be calculated from the following set of equations [11].
Where N is the total number of elements and M is the number of normal elements. The K n , K s are shear and normal stiffness's of a joint element and X n , X s are the total joint shear and normal deformations respectively. The maximum deformation of a joint element is limited by its closure value (See Figure 1) which is the hydraulic conductivity (0.1 mm for all the simulations shown in this paper.

Fracture propagation
The fracture tips are allowed to propagate when mode-I stress intensity factor K I is equal to fracture toughness K IC according to LEFM [13]. K I , K II are calculated using displacement discontinuity obtained at the center of the crack tip element [12]. Continued injection of fluid into the fracture will increase the stress intensity at fracture tip and eventually cause it to propagate. In sequential fracturing this can be achieved by changing the pressure boundary condition at fracture tip iteratively till the propagation condition is satisfied. We recognize that fully fluid filled crack has a singular pressure at the fracture tip and requires and asymptotic analysis. For computational purposes we choose to have a finite pressure boundary condition at the last grid block of finite difference scheme for fluid flow inside the fracture which is fracture tip [14]. Whereas for simultaneous fracturing, it is not feasible to iteratively find the pressure distribution in the fracture to satisfy K I =K IC since more than one fracture is growing at a given time. In this case zero net pressure boundary condition is used at the fracture tip [15,16]. The fluid is injected until K I =K IC to satisfy the propagation condition. The crack propagation path is calculated using the method of [10] as implemented in [17] in which the crack propagation direction relies on the maximum principal tensile stress criterion so that one can use the ratio of the stress intensity factors to compute the angle at which the crack will grow.  II   I  I  II  II  II  I  II  II if K

Fluid flow
The fluid flow inside the fracture assumed to be laminar and is modeled using flow through parallel plates equation often called as cubic law. 3 12 Where q is the volumetric flow rate, µ is dynamic viscosity of the fluid and w is fracture width.
The fluid is assumed to be Newtonian and incompressible. The continuity equation (eq 6) along with cubic law governs the fluid flow inside the fracture.
Where A is the cross-sectional area of the fracture and q L is fluid leak-off volume rate per unit length of the fracture. Because of the ultra-low permeability nature of shale reservoir matrix, leaf-off is assumed to be zero in these calculations [7]. After every time step in the simulation global mass balance is satisfied. The partial differential equation 7 obtained from substituting eq 5 (cubic law) in to eq 6 (continuity) is solved using the finite difference approximation with the following boundary conditions. 3 12 At the well, the injection rate is specified At the fracture tip finite pressure P tip is specified.

Simulation examples
Example-1: Sequential fracturing with fracture spacing 3 m (9.84 ft.) In this example we consider sequential fracturing of a horizontal well. Fracture that is generated first (i.e., Stage-1 fracture) is subjected to a constant pressure while injecting into subsequent fractures. Figure 2 shows the geometry of horizontal wellbore and transverse fractures from top view (the simulations in this paper considers a 2D plane strain similar to KGD model. Thus in its current form, the model cannot consider the stress shadow between the fractures due to their height). The properties used in simulation are given in Table 1 for example-1.  Since zero fluid leak-off is assumed into the formation, the entire fluid inside the stage-1 is assumed to re-distribute after injection is ceased into it and is set to a constant average pressure (assuming no gravity effect) corresponding to the target length. Then the Stage-2 transverse fracture is created at a distance of 3m from the stage-1 fracture. Figure 2 show the stage-2 fracture turns towards the stage-1 fracture. This is due to the altered stress distribution in the rock. From the fracture propagation criterion, the fracture propagates in the direction perpendicular to the maximum principal tensile stress. Figure 2 (right side) shows the distribution of maximum principal compressive stresses. The Stage-2 fracture appears to be oriented in the direction of maximum principal compressive stresses. The opening of stage-2 fracture causes a reduction in width along the center of Stage-1 fracture (Figure 3 left side.). The compressive stress due to opening of Stage-1 fracture has not been reduced by the Stage 2 crack near the tips of the former, causing the Stage-2 fracture to curve towards the tips of Stage-1 fracture for a few steps as it eventually follows the maximum in-situ compression direction. The widths of stage-2 fracture (Figure 3 right side) increased as it grows in length.   Tables 1&2 and fluid properties are given in Table 3. The spacing between the fractures is increased to 7 m (more than twice of previous examples) to observe the cuving behavior of Stage-2 fracture. The results from Figure 6 shows Stage-2 fracture curves away from Stage-1 fracture in both scenarios (i.e. Stage-1 fracture pressurized/propped). The curving of Stage-2 fracture near propped Stage-1 fracture is stronger than near pressurized Stage-1 fracture. This phenomenon is expected as from Figure 7 we can see the tips of pressurized Stage-1 fracture remained open after initiation of Stage-2 fracture. The opening of Stage-1 fracture near its tips created attraction towards its tips while its opening near its center created repulsion between the fractures. In this case the repulsion between the fractures slightly dominated the attraction between their tips, which lead to a slightly curve away Stage-2 fracture. This phenomenon also leads to straight Stage-2 fracture when the repulsion between fractures and attracion between the tips balances out. On the other hand Stage-2 fracture curve away from propped Stage-1 fracture as there is no attraction between the tips (Satge-1 tips closed) from Figure 7.   In this example simultaneous propagation of five hydraulic fractures with initial lengths of 5 m each and 10 m spacing between them is considered. The input parameters used in this simulation are shown in Table 4. The results from Figure 8 show that when fluid is injected into the wellbore at the injection point (Figure 8 left side) the outer fractures (i.e., Frac-1 & Frac-5) tend to grow more than the rest of the fractures. This behavior is expected because the stress shadow on the outer fractures will be less than the fractures inside. On the other hand the results show growth of center fracture (i.e., Frac-3) is more than the fractures Frac-2 & Frac-4. This effect is seen when the outer fractures start to propagate more rapidly than the remaining fractures in the fracture network. The larger outer fractures exert more stress shadow on the fracture near to them (in this case Frac-2 & Frac-4) which will inhibit the growth of these fractures more than the center fracture (i.e.Frac-3). As the outer fractures grows large enough the stress shadow over the remaining fractures increase and completely suppress their growth after sometime.

Conclusions
Several numerical examples are considered to show the behavior of closely spaced hydraulic fractures in sequential and simultaneous fracturing. Two models are presented in sequential fracturing; one considering the previously created fracture is kept pressurized with fluid and other considering the previously created fracture filled with proppant. Numerical results show the fracture geometries of later fractures are dependent on boundary conditions of previous fracture (i.e., pressurized fracture or propped fracture) and injection conditions. With variation in spacing and boundary condition on the previous fracture, the later fracture is observed to curve in/out from the previously created fractures. Later fractures curving from the previous fracture depends on the attraction forces between tips and repelling forces between the centers of the fractures. It was observed that fracture created near to propped fracture tend to curve away more as there is no attraction between the fractures' tips whereas fracture created near pressurized fracture tend to curve. In simultaneous propagation of hydraulic fractures, the outer fractures dominate the inner fracture in growth. The center fractures usually stop after they reached certain length due to the stress shadow between them. These simulations are useful for horizontal wellbore stimulation design and required spacing conditions to acquire the desired fracture lengths, proppant placement, and production rates.

Author details
Varahanaresh Sesetty and Ahmad Ghassemi Mewbourne School of Petroleum and Geological Engineering, The University of Oklahoma, USA