**Abstract**

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.

### **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

© 2013 Sesetty and Ghassemi; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

**2. Model development**

these model components are briefly described below.

**3. Displacement discontinuity method**

s

s

displacement discontinuity of each fracture element.

*and Ann ij*

crack surfaces is given by *uy*(*x*)=*Dy*(*<sup>x</sup>*

special element used herein is given in [12].

*Ass ij* , *Asn ij* , *Ans ij*

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

Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

http://dx.doi.org/10.5772/56309

681

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

(for 1 to N)

are the influence coefficients, representing the stresses due to constant

where a is half length of the crack tip element, Dy is

(1)

stresses on each fracture element is found by using superposition [11]:

*N N ij j ij j i s ss s sn n j j N N ij j ij j i n ns s nn n j j*

*AD AD*

*a* ) 1/2

*AD AD i*

shear and normal DD elements. The above system of linear equations can be solved for

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

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

=+ =

1 1

= =

å å

= +

1 1

= =

å å

In order to determine the optimum spacing and optimum staging between fractures produc‐ tion 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 simula‐ tions 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 impor‐ tance 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**

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,

In order to determine the optimum spacing and optimum staging between fractures produc‐ tion 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

[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 simula‐ tions 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 impor‐

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

which reduced the possibility of having improper fracture propagation.

fracturing process.

680 Effective and Sustainable Hydraulic Fracturing

tance of spacing optimization.

proppant [1].

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.

### **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 *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]:

$$\begin{aligned} \sigma\_s^i = \sum\_{j=1}^N A\_{ss}^{ij} \stackrel{j}{D}\_s + \sum\_{j=1}^N A\_{sn}^{ij} \stackrel{j}{D}\_n\\ \sigma\_n^i = \sum\_{j=1}^N A\_{ns}^{ij} \stackrel{j}{D}\_s + \sum\_{j=1}^N A\_{nn}^{ij} \stackrel{j}{D}\_n \quad \text{(for } i = 1 \text{ to N)} \end{aligned} \tag{1}$$

*Ass ij* , *Asn ij* , *Ans ij and Ann ij* are the influence coefficients, representing the stresses due to constant 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 *uy*(*x*)=*Dy*(*<sup>x</sup> a* ) 1/2 where a is half length of the crack tip element, Dy 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].

#### **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 non-linearly, here it is assumed to deform linear (linear model in Figure 1) with the stress for simplicity.

Where N is the total number of elements and M is the number of normal elements. The *Kn*, *Ks* are shear and normal stiffness's of a joint element and *Xn*, *Xs* 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 simu‐

Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

http://dx.doi.org/10.5772/56309

683

The fracture tips are allowed to propagate when mode-I stress intensity factor KI is equal to fracture toughness KIC according to LEFM [13]. KI, KII are calculated using displace‐ ment discontinuity obtained at the center of the crack tip element [12]. Continued injec‐ tion 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 KI=KIC since more than one fracture is growing at a given time. In this case zero net pressure boundary condi‐ tion is used at the fracture tip [15, 16]. The fluid is injected until KI=KIC 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

2

*II*

*if K*

*II*

(4)

*if K*

0 0

4

*I I II*

<sup>ï</sup> - + ç ÷ ç ÷ ç ÷ <sup>=</sup> <sup>í</sup> è ø <sup>¹</sup> <sup>ï</sup>

The fluid flow inside the fracture assumed to be laminar and is modeled using flow through

*K K <sup>K</sup>*

<sup>ì</sup> <sup>=</sup> <sup>ï</sup> æ ö <sup>ï</sup> æ ö ç ÷ <sup>ï</sup>

sgn( ) 8 (, ) 2arctan <sup>0</sup>

ïî è ø

*I II II II*

*K K K K*

lations shown in this paper.

**5. Fracture propagation**

to compute the angle at which the crack will grow.

ï ï

parallel plates equation often called as cubic law.

q

**6. Fluid flow**

**Figure 1.** Goodman Joint model and a linear joint model. In Goodman model the closure reaches an asymptotic value at high values of normal stress

Given the far field stresses (*σij* )0 *<sup>∞</sup>*and stresses acting on the joint element, the total joint defor‐ mation 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].

$$\begin{aligned} -\left(\stackrel{i}{\sigma}\_{s}\right)\_{0}^{\circ} &= \sum\_{j=1}^{N} \{A\_{ss}\stackrel{ij}{X}\_{s} + A\_{sn}^{\circ}\stackrel{ij}{X}\_{n}\} \\ -\left(\stackrel{i}{\sigma}\_{n}\right)\_{0}^{\circ} &= \sum\_{j=1}^{N} \{A\_{ns}\stackrel{ij}{X}\_{s} + A\_{nn}^{\circ}\stackrel{ij}{X}\_{n}\} \text{ (1 \le i \le M\text{)} }\end{aligned} \tag{2}$$

$$\begin{aligned} -\left(\stackrel{i}{\sigma}\_{s}\right)\_{0}^{\alpha} &= K\_{s}X\_{s} + \sum\_{j=1}^{N} \left(A\_{ss}X\_{s} + A\_{sn}^{ij}X\_{n}\right) \\ -\left(\stackrel{i}{\sigma}\_{n}\right)\_{0}^{\alpha} &= K\_{n}X\_{n} + \sum\_{j=1}^{N} \left(A\_{ns}X\_{s} + A\_{nn}X\_{n}\right) \text{ ( $M+1 \le i \le N$ )} \end{aligned} \tag{3}$$

Where N is the total number of elements and M is the number of normal elements. The *Kn*, *Ks* are shear and normal stiffness's of a joint element and *Xn*, *Xs* 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 simu‐ lations shown in this paper.

#### **5. Fracture propagation**

**4. Joint model**

682 Effective and Sustainable Hydraulic Fracturing

at high values of normal stress

Given the far field stresses (*σij*

)0

1 0

¥

¥

1 0


1 0

*j*


=

å

1 0

*j*

=

*N ij j ij j i ii s s s ss s sn n j N ij j ij j i ii n n n ns s nn n j*

=

å

=

å

å

calculated from the following set of equations [11].

s

æ ö

ç ÷ è ø æ ö

s

s

æ ö

ç ÷ è ø æ ö

¥

¥

s

è ø

ç ÷ è ø

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.

**Figure 1.** Goodman Joint model and a linear joint model. In Goodman model the closure reaches an asymptotic value

mation 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

( , (1 )

*AX AX i M*

( ) , ( 1 )

( )

( )


*KX A X A X M i N*

*AX AX*


*N ij j ij j i s ss s sn n*

*N ij j ij j i n ns s nn n*

*KX A X A X*

*<sup>∞</sup>*and stresses acting on the joint element, the total joint defor‐

(2)

(3)

The fracture tips are allowed to propagate when mode-I stress intensity factor KI is equal to fracture toughness KIC according to LEFM [13]. KI, KII are calculated using displace‐ ment discontinuity obtained at the center of the crack tip element [12]. Continued injec‐ tion 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 KI=KIC since more than one fracture is growing at a given time. In this case zero net pressure boundary condi‐ tion is used at the fracture tip [15, 16]. The fluid is injected until KI=KIC 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.

$$\theta(K\_I, K\_{II}) = \begin{cases} 0 & \text{if } K\_{II} = 0\\ 2 \arctan\left(\frac{\frac{K\_I}{K\_{II}} - \text{sgn}(K\_{II})\sqrt{\left(\frac{K\_I}{K\_{II}}\right)^2 + 8}}{4}\right) & \text{if } K\_{II} \neq 0 \\ \end{cases} \tag{4}$$

#### **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.

$$q = -\frac{w^3 \hbar}{12\,\mu} \frac{\partial p}{\partial x} \tag{5}$$

**Parameter Value Units** Young's modulus 27 GPa

Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

http://dx.doi.org/10.5772/56309

685

σH/ σh (Max/Min In-situ horizontal stress) 5/4 MPa Injection rate (stage-1)/(stage-2) 20/40 bpm Viscosity 1 cP Fracture height 30 ft

Fracture toughness 2 MPa.m1/2

**Figure 2.** Geometry of the Stage-2 fracture near pressurized Stage-1 fracture after they reached their target lengths

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 perpen‐ dicular 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

Poisson's ratio 0.25

and maximum principal stress distribution around them at an instant

of stage-2 fracture (Figure 3 right side) increased as it grows in length.

**Table 1.** Input parameters used in example-1

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.

$$\frac{\partial q}{\partial \mathbf{x}} + \frac{\partial A}{\partial t} + q\_L = 0 \tag{6}$$

Where A is the cross-sectional area of the fracture and qL 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.

$$\frac{\partial w}{\partial t} = \frac{\partial}{\partial x} \left( \frac{w^3}{12\,\mu} \frac{\partial p}{\partial x} \right) \tag{7}$$

At the well, the injection rate is specified

*q*(0, *t*)=*Q*<sup>0</sup>

At the fracture tip finite pressure Ptip is specified.

*P*(*L* , *t*)=*Ptip*

#### **7. 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.


3 12 *w h <sup>p</sup> <sup>q</sup>* m*x*

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

*q A <sup>q</sup> x t* ¶ ¶

+ +=

0 *<sup>L</sup>*

Where A is the cross-sectional area of the fracture and qL 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

> 3 12 *w w p tx x* m

è ø

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

¶ ¶ æ ö ¶ <sup>=</sup> ç ÷ ¶¶ ¶ ç ÷

with cubic law governs the fluid flow inside the fracture.

the following boundary conditions.

684 Effective and Sustainable Hydraulic Fracturing

At the well, the injection rate is specified

At the fracture tip finite pressure Ptip is specified.

**Example-1:** Sequential fracturing with fracture spacing 3 m (9.84 ft.)

*q*(0, *t*)=*Q*<sup>0</sup>

*P*(*L* , *t*)=*Ptip*

example-1.

**7. Simulation examples**

¶ = - ¶ (5)

¶ ¶ (6)

(7)

**Figure 2.** Geometry of the Stage-2 fracture near pressurized Stage-1 fracture after they reached their target lengths and maximum principal stress distribution around them at an instant

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 perpen‐ dicular 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.

**Parameter Value Units** Injection rate (stage-1)/(stage-2) 20/20 bpm Proppant normal stiffness/shear stiffness 15/15 GPa/m

**Example-3:** Sequential fracturing with fracture spacing 7 m (20 ft)

attraction between the tips (Satge-1 tips closed) from Figure 7.

**Parameter Value Units** Injection rate (stage-1)/(stage-2) 20/20 bpm Viscosity 1 cP

**Figure 5.** Widths of Stage-1 fracture at various instants

**Table 3.** Input parameters used in example-3

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 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

Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

http://dx.doi.org/10.5772/56309

687

**Table 2.** Input parameters used in example-2

**Figure 3.** Widths of Stage-1 and Stage-2 fractures at various instants

#### **Example-2:** Sequential fracturing while keeping previously created fracture propped

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 in-situ stress and rock properties are kept same as in example-1. Keeping the spacing between the two fractures 3 m (9.84 ft.) same as in example-1, the simulation results in this example shows the Stage-2 fracture curve away from Stage-1 fracture (Figure 4). The maximum principal com‐ pressive stress diagram shows (Figure 4 right side) a huge compression zone near the center of both fractures and towards the right of Stage-2 fracture which caused it to curve away from Stage-1 fracture. The creation of Stage-2 fracture caused the tips of Stage-1 fracture to close (Figure 5) due to the stress shadow between them. The center part of Stage-1 fracture remained open and contributed to large compressive stresses around the center of fractures. Since the tips of Stage-1 fracture are closed in this case there is no attraction towards tips (higher tensile zones) in this case.

**Figure 4.** Geometry of the Stage-2 fracture near propped Stage-1 fracture after they reached their target lengths and maximum principal stress distribution around them at an instant


**Table 2.** Input parameters used in example-2

**Figure 3.** Widths of Stage-1 and Stage-2 fractures at various instants

686 Effective and Sustainable Hydraulic Fracturing

maximum principal stress distribution around them at an instant

zones) in this case.

**Example-2:** Sequential fracturing while keeping previously created fracture propped

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 in-situ stress and rock properties are kept same as in example-1. Keeping the spacing between the two fractures 3 m (9.84 ft.) same as in example-1, the simulation results in this example shows the Stage-2 fracture curve away from Stage-1 fracture (Figure 4). The maximum principal com‐ pressive stress diagram shows (Figure 4 right side) a huge compression zone near the center of both fractures and towards the right of Stage-2 fracture which caused it to curve away from Stage-1 fracture. The creation of Stage-2 fracture caused the tips of Stage-1 fracture to close (Figure 5) due to the stress shadow between them. The center part of Stage-1 fracture remained open and contributed to large compressive stresses around the center of fractures. Since the tips of Stage-1 fracture are closed in this case there is no attraction towards tips (higher tensile

**Figure 4.** Geometry of the Stage-2 fracture near propped Stage-1 fracture after they reached their target lengths and

#### **Example-3:** Sequential fracturing with fracture spacing 7 m (20 ft)

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 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.

**Figure 5.** Widths of Stage-1 fracture at various instants


**Table 3.** Input parameters used in example-3

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

Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

http://dx.doi.org/10.5772/56309

689

**Figure 8.** Fracture network obtained after outer fractures reached their target lengths and maximum principal stress

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

growth after sometime.

distribution around them

**8. Conclusions**

**Figure 6.** Comparison of Stage-2 fracture geometry after reaching its target length with different boundary condi‐ tions for Stage-1 fracture (i.e. pressurized/propped)

**Figure 7.** Comparison of pressurized/propped Stage-1 fracture widths after Stage-2 fracture reached target length

#### **Example-4:** Simultaneous propagation of hydraulic fractures.


**Table 4.** Input parameters used in example-4

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.

**Figure 8.** Fracture network obtained after outer fractures reached their target lengths and maximum principal stress distribution around them

### **8. Conclusions**

**Figure 6.** Comparison of Stage-2 fracture geometry after reaching its target length with different boundary condi‐

**Figure 7.** Comparison of pressurized/propped Stage-1 fracture widths after Stage-2 fracture reached target length

**Example-4:** Simultaneous propagation of hydraulic fractures.

Poisson's ratio 0.22

**Table 4.** Input parameters used in example-4

**Parameter Value Units** Young's modulus 30 GPa

σH/ σh (Max/Min In-situ horizontal stress) 5/4 MPa Injection rate 80 bpm Viscosity 7 cP Fracture height 100 ft

Fracture toughness 2 MPa.m1/2

tions for Stage-1 fracture (i.e. pressurized/propped)

688 Effective and Sustainable Hydraulic Fracturing

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.

draulic Fracturing Technology Conference, The Woodlands, Texas, January (2009). ,

Numerical Simulation of Sequential and Simultaneous Hydraulic Fracturing

http://dx.doi.org/10.5772/56309

691

[9] Batchelor, G. K. An Introduction to Fluid Dynamics. Cambridge, UK: Cambridge

[10] Stone, T. J, & Babuska, I. A numerical method with a posteriori error estimation for determining the path taken by a propagating crack, Computer Methods in Applied

[11] Crouch, S. L, & Starfield, A. M. Boundary Element Methods in Solid Mechanics. Lon‐

[12] Yan, X. A Special Crack Tip Displacement Discontinuity Element, Mechanics Re‐

[13] Rice, J. R. Mathematical analysis in the mechanics of fracture. In Fracture, An Ad‐ vanced Treatise, Leibowitz H (ed.). Academic Press: New York, (1968). , 172-308. [14] Murdoch, L. C, & Germanovich, L. N. Analysis of a deformable fracture in permea‐ ble material. International Journal for Numerical and Analytical Methods in Geome‐

[15] Nordgren, R. P. Propagation of a Vertical Hydraulic Fracture. Society of Petroleum

[16] Weng, X, Kresse, O, Cohen, C, Wu, R, & Gu, H. Schlumberger. Modeling of Hydraul‐ ic Fracture Network Propagation in a Naturally Fractured Formation. SPE 140253, SPE Hydraulic Fracturing Technology Conference and Exhibition, Woodlands,

[17] Tarasvo, S, & Ghassemi, A. Propagation of a system of cracks under thermal stress. Presented at the 45th US rock mechanics Symposium, San Fransisco, (2010).

Mechanics and Engineering, (1998). , 1998(160), 245-271.

19-21.

University Press; (1967).

don: George Allen & Unwin; (1983).

Engineers Journal, (1972). , 1972(12), 306-314.

Texas, USA, January, (2011). , 24-26.

search Communications, (2004).

chanics, (2006).
