Input parameters used in example1
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 simultaneous and sequential hydraulic fracture operations. The sequential fracturing model considers different boundary conditions for the previously created fractures (constant pressure restricting the flow back between stages and proppantfilled). A series of examples are presented to study the effect of fracture spacing to show the importance of spacing optimization. The results show the fracture path is not only affected by fracture spacing but also by the boundary conditions on the previously created fractures.
1. 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 stressshadow 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 DDbased 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].
2. 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 nonlinear problem that is solved using NewtonRaphson 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 tensilestress criterion of [10]. Each of these model components are briefly described below.
3. 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
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
4. 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 nonlinearly, here it is assumed to deform linear (linear model in Figure 1) with the stress for simplicity.
Given the far field stresses
Where N is the total number of elements and M is the number of normal elements. The
5. Fracture propagation
The fracture tips are allowed to propagate when modeI 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.
6. 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.
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 crosssectional area of the fracture and q_{L} is fluid leakoff volume rate per unit length of the fracture. Because of the ultralow permeability nature of shale reservoir matrix, leafoff 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.
At the well, the injection rate is specified
At the fracture tip finite pressure P_{tip} is specified.
7. Simulation examples
In this example we consider sequential fracturing of a horizontal well. Fracture that is generated first (i.e., Stage1 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 example1.



Young’s modulus  27  GPa 
Poisson’s ratio  0.25  
σ_{H/} σ_{h} (Max/Min Insitu horizontal stress)  5/4  MPa 
Injection rate (stage1)/(stage2)  20/40  bpm 
Viscosity  1  cP 
Fracture height  30  ft 
Fracture toughness  2  MPa.m^{1/2} 
Since zero fluid leakoff is assumed into the formation, the entire fluid inside the stage1 is assumed to redistribute 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 Stage2 transverse fracture is created at a distance of 3m from the stage1 fracture. Figure 2 show the stage2 fracture turns towards the stage1 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 Stage2 fracture appears to be oriented in the direction of maximum principal compressive stresses. The opening of stage2 fracture causes a reduction in width along the center of Stage1 fracture (Figure 3 left side.). The compressive stress due to opening of Stage1 fracture has not been reduced by the Stage 2 crack near the tips of the former, causing the Stage2 fracture to curve towards the tips of Stage1 fracture for a few steps as it eventually follows the maximum insitu compression direction. The widths of stage2 fracture (Figure 3 right side) increased as it grows in length.
This example considers sequential fracturing keeping the previously created fracture propped. The fluid and proppant properties for this simulation are given in Table 2. The insitu stress and rock properties are kept same as in example1. Keeping the spacing between the two fractures 3 m (9.84 ft.) same as in example1, the simulation results in this example shows the Stage2 fracture curve away from Stage1 fracture (Figure 4). The maximum principal compressive stress diagram shows (Figure 4 right side) a huge compression zone near the center of both fractures and towards the right of Stage2 fracture which caused it to curve away from Stage1 fracture. The creation of Stage2 fracture caused the tips of Stage1 fracture to close (Figure 5) due to the stress shadow between them. The center part of Stage1 fracture remained open and contributed to large compressive stresses around the center of fractures. Since the tips of Stage1 fracture are closed in this case there is no attraction towards tips (higher tensile zones) in this case.



Injection rate (stage1)/(stage2)  20/20  bpm 
Proppant normal stiffness/shear stiffness  15/15  GPa/m 
This example compares sequential fracturing while keeping previously created fracture pressurized and propped. The rock and proppant properties are given in 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 Stage2 fracture. The results from Figure 6 shows Stage2 fracture curves away from Stage1 fracture in both scenarios (i.e. Stage1 fracture pressurized/propped). The curving of Stage2 fracture near propped Stage1 fracture is stronger than near pressurized Stage1 fracture. This phenomenon is expected as from Figure 7 we can see the tips of pressurized Stage1 fracture remained open after initiation of Stage2 fracture. The opening of Stage1 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 Stage2 fracture. This phenomenon also leads to straight Stage2 fracture when the repulsion between fractures and attracion between the tips balances out. On the other hand Stage2 fracture curve away from propped Stage1 fracture as there is no attraction between the tips (Satge1 tips closed) from Figure 7.



Injection rate (stage1)/(stage2)  20/20  bpm 
Viscosity  1  cP 



Young’s modulus  30  GPa 
Poisson’s ratio  0.22  
σ_{H/} σ_{h} (Max/Min Insitu horizontal stress)  5/4  MPa 
Injection rate  80  bpm 
Viscosity  7  cP 
Fracture height  100  ft 
Fracture toughness  2  MPa.m^{1/2} 
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., Frac1 & Frac5) 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., Frac3) is more than the fractures Frac2 & Frac4. 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 Frac2 & Frac4) which will inhibit the growth of these fractures more than the center fracture (i.e.Frac3). As the outer fractures grows large enough the stress shadow over the remaining fractures increase and completely suppress their growth after sometime.
8. 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.
References
 1.
Guimares de Carvalho. Horizontal Well Completion and Stimulation Techniques, SPE 108075, presented at the Latin American & Caribbean Petroleum Engineering Conference, Buenos Aires,Rodrigues V. F FNeumann L Torres D. S 15 18 April,2007  2.
Geomechanics aspects of multiple fracturing of horizontal and vertical wells. SPE Drilling& Completion,Soliman M. Y Loyd E David A 2008 217 228  3.
Parameters affecting the interaction among closely spaced hydraulic fractures, SPE 140426, SPE Hydraulic Fracturing Technology Conference and Exhibition, Woodlands, Texas, USA,Bunger A. P andX Zhang R. G Jeffery 24 26 January,2011  4.
Experimental study of hydraulic fracture geometry initiated from horizontal wells. SPE Annual Technical Conference and Exhibition. San Antonio, TX, USA,El Rabba W 1989  5.
Modeling Simultaneous Growth of Multiple Hydraulic Fractures and Their Interaction With Natural Fractures. Presented at SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, USA, andOlson J. E D. A Taleghani 19 21 January,2009  6.
Strategies to Minimize Frac Spacing and Simulate Natural Fractures in Horizontal Completions. SPE 146104, presented at SPE Annual Technical Conference and Exhibition held in Denver, Colorado, USA, 30 October2 November,Roussel N. P Sharma M. M 2011  7.
Hyunil Jo Optimizing fracture spacing to induce complex fractures in a hydraulically fractured horizontal wellbore.2012 SPE 154930.  8.
Waters, Barry K.D, Robert C.D, Kenneth J.K, Lance A., Simultaneous Hydraulic Fracturing of Adjacent Horizontal Wells in the Woodford Shale, SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas,George A 19 21 January2009  9.
An Introduction to Fluid Dynamics. Cambridge, UK: Cambridge University Press;Batchelor G. K 1967  10.
A numerical method with a posteriori error estimation for determining the path taken by a propagating crack, Computer Methods in Applied Mechanics and Engineering,Stone T. J Babuska I 1998 1998 160 245 271  11.
Boundary Element Methods in Solid Mechanics. London: George Allen & Unwin; andCrouch S. L Starfield A. M 1983  12.
A Special Crack Tip Displacement Discontinuity Element, Mechanics Research Communications,Yan X 2004  13.
Mathematical analysis in the mechanics of fracture. In Fracture, An Advanced Treatise, Leibowitz H (ed.). Academic Press: New York,Rice J. R 1968 172 308  14.
Analysis of a deformable fracture in permeable material. International Journal for Numerical and Analytical Methods in Geomechanics,Murdoch L. C L. N Germanovich 2006  15.
Propagation of a Vertical Hydraulic Fracture. Society of Petroleum Engineers Journal,Nordgren R. P 1972 1972 12 306 314  16.
Schlumberger. Modeling of Hydraulic Fracture Network Propagation in a Naturally Fractured Formation. SPE 140253, SPE Hydraulic Fracturing Technology Conference and Exhibition, Woodlands, Texas, USA,Weng X O Kresse C Cohen andR Wu H Gu 24 26 January,2011  17.
Propagation of a system of cracks under thermal stress. Presented at the 45^{th} US rock mechanics Symposium, San Fransisco, andTarasvo S A Ghassemi 2010