**1. Introduction**

Hydraulic fracturing or fracture stimulation establishes conductive fractures hydraulically from a horizontal well in the tight formation/reservoir. They provide large surface area contact with the formation and thus facilitate the production of oil and gas, as evident from the experience in North Americas [1]. Multiple fractures are now placed in sequential stages in

horizontal wells. Typically, four or more fractures are pumped simultaneously into a single frac stage, and it is not uncommon to place 20 to 40 frac stages in a single horizontal well. Attempts have also been made to simultaneously fracture two adjacent horizontal wells to generate complex fracture networks thought to be beneficial for production [2, 3].

fractures is captured meticulously by means of boundary integral formulation with dislocation segments solution techniques. In the model, fracture growth and fluid flow equations are solved in a coupled manner via a proprietary, robust and efficient algorithm where mass conservation is strictly observed. This new non-planar 3D model was first described in [19]. It substantially enhances the capability of a plane strain 2D model presented in [7]. The present paper complements the published results [7, 19] and extends to considering fluid viscosity and

The Geomechanical Interaction of Multiple Hydraulic Fractures in Horizontal Wells

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

663

The complete theory of fracture interaction for an arbitrary number of fractures is quite cumbersome. Therefore, a 'simplified' analytical model has been developed to describe a set of fractures growing in a viscous mass-transfer dominated regime. In this model, two types of fracture growth patterns are defined. First is the 'compact' fracture growth where interaction and interference between fractures are minimum. The second is 'diffuse' fracture growth where some fractures will either terminate or collapse into each other because of minimum insitu stress rotation (Figure 1). In diffuse growth, we distinguish different zones. The first zone is the closest to the well. The next zone, farther from the well, contain less fractures because some of them have either terminated or merged into other fractures, as shown at the right-

The transition from one growth pattern to another is delineated by a 'critical' fracture spacing, influenced by the original principal stress difference in the un-fractured forma‐ tion. In the case of diffuse growth, different zones are depicted by the dashed lines on the right-hand picture of Figure 1. The Average permeability of the fracture set in the zone

3

*n*

where *Ln* is spacing between the fractures in the "n" zone, and average aperture of fractures

Here, *pn* is net pressure, *E* is Young's modulus, *Rn* is a radius of a penny-shape fracture, and *a* is a numerical coefficient (of approximately 1) and depends on Poisson's ratio. For simplicity, it is assumed that there is no leak-off of fluid into the extremely tight formation. Therefore, the

*<sup>L</sup>* (1)

*<sup>E</sup>* (2)

=

*<sup>w</sup> <sup>k</sup>*

*n*

12 *n*

= *<sup>n</sup> n n <sup>p</sup> waR*

height growth of multiple fractures.

**2. Analytical model**

hand side of Figure 1.

in the zone "n" is:

"n" is given by (see e.g. [20]):

average fluid density, *m*, in the fractures becomes:

It is important to consider the changes to subsurface in-situ stresses induced by fracture stimulation. This fracturing induced formation stress perturbation is also called 'stress shadowing' and needs to be understood and quantified for optimization and to avoid problems such as fracture 'screen-out'. Added complexity is caused by the fact that the fractures are created in an extremely tight formation - where the fracture pumping operations are often conducted within a time span that is in the same order of magnitude needed for frac fluid to leak-off and the resultant stress perturbation to completely dissipate. Moreover, the stress perturbation caused by inclusion of millions of pounds of proppant in the formation does not dissipate.

Observations from substantial amounts of field data show that generally 25% or more of the fractures placed are ineffective. In fact, recent fracture stimulation surveillance with advance technologies such as microseismic monitoring and fibre optics temperature and strain sensing [4-6] have shown that multiple fractures do not grow and develop in the same way. Some of the fractures are pre-maturely 'terminated' during the treatment. Nonetheless, to increase production, there is a tendency to use longer horizontal wells, place more fractures per stage, longer fractures and more stages of fractures. Consequently, more fractures with closer spacing and larger volumes of frac fluid are employed. This translates to higher pumping rates, higher treating pressure and surface horse power, and consumes more materials (frac fluid, chemicals and proppant). The optimal design of hydraulic fractures clearly necessitates an engineering optimization, considering production, treatment cost and the feasibility of placing these closely spaced conductive fractures in the subsurface. The latter requires knowledge of key parameters such as in-situ stress, rock stiffness and strength, frac fluid rheology and leak-off behaviour, and importantly, the impact of stress shadows on multiple fracture growth in heterogeneous formations [7-13]. Some experts have aimed to provide general guidelines for the design of optimal fracture placement based on simplified analytical models whilst others have devel‐ oped numerical models of various degree of sophistication to account for multiple fracture interaction and design [14-18].

In-situ stress perturbation can potentially result in the growth of complex non-planar fractures. An analytical model is developed to highlight some of the salient features of multiple nonplanar hydraulic fractures interaction. The benefit of using an analytical model is that it provides immediate insights into the controlling parameters on fractures interaction, and guides further numerical analysis for stimulation optimization. Here, we present only the impact of fracture spacing and in-situ stress difference on fracture growth pattern. Other parameters that affects fracture growth pattern are evident from the model, but they will not be described here.

Emboldened by the results of our analytical model, a new numerical model based on rigorous mechanistic formulation has been devised to allow for more realistic solution of field problems. This is a non-planar 3D numerical model, where the interaction of multiple non-planar fractures is captured meticulously by means of boundary integral formulation with dislocation segments solution techniques. In the model, fracture growth and fluid flow equations are solved in a coupled manner via a proprietary, robust and efficient algorithm where mass conservation is strictly observed. This new non-planar 3D model was first described in [19]. It substantially enhances the capability of a plane strain 2D model presented in [7]. The present paper complements the published results [7, 19] and extends to considering fluid viscosity and height growth of multiple fractures.

#### **2. Analytical model**

horizontal wells. Typically, four or more fractures are pumped simultaneously into a single frac stage, and it is not uncommon to place 20 to 40 frac stages in a single horizontal well. Attempts have also been made to simultaneously fracture two adjacent horizontal wells to

It is important to consider the changes to subsurface in-situ stresses induced by fracture stimulation. This fracturing induced formation stress perturbation is also called 'stress shadowing' and needs to be understood and quantified for optimization and to avoid problems such as fracture 'screen-out'. Added complexity is caused by the fact that the fractures are created in an extremely tight formation - where the fracture pumping operations are often conducted within a time span that is in the same order of magnitude needed for frac fluid to leak-off and the resultant stress perturbation to completely dissipate. Moreover, the stress perturbation caused by inclusion of millions of pounds of proppant in the formation does not

Observations from substantial amounts of field data show that generally 25% or more of the fractures placed are ineffective. In fact, recent fracture stimulation surveillance with advance technologies such as microseismic monitoring and fibre optics temperature and strain sensing [4-6] have shown that multiple fractures do not grow and develop in the same way. Some of the fractures are pre-maturely 'terminated' during the treatment. Nonetheless, to increase production, there is a tendency to use longer horizontal wells, place more fractures per stage, longer fractures and more stages of fractures. Consequently, more fractures with closer spacing and larger volumes of frac fluid are employed. This translates to higher pumping rates, higher treating pressure and surface horse power, and consumes more materials (frac fluid, chemicals and proppant). The optimal design of hydraulic fractures clearly necessitates an engineering optimization, considering production, treatment cost and the feasibility of placing these closely spaced conductive fractures in the subsurface. The latter requires knowledge of key parameters such as in-situ stress, rock stiffness and strength, frac fluid rheology and leak-off behaviour, and importantly, the impact of stress shadows on multiple fracture growth in heterogeneous formations [7-13]. Some experts have aimed to provide general guidelines for the design of optimal fracture placement based on simplified analytical models whilst others have devel‐ oped numerical models of various degree of sophistication to account for multiple fracture

In-situ stress perturbation can potentially result in the growth of complex non-planar fractures. An analytical model is developed to highlight some of the salient features of multiple nonplanar hydraulic fractures interaction. The benefit of using an analytical model is that it provides immediate insights into the controlling parameters on fractures interaction, and guides further numerical analysis for stimulation optimization. Here, we present only the impact of fracture spacing and in-situ stress difference on fracture growth pattern. Other parameters that affects fracture growth pattern are evident from the model, but they will not

Emboldened by the results of our analytical model, a new numerical model based on rigorous mechanistic formulation has been devised to allow for more realistic solution of field problems. This is a non-planar 3D numerical model, where the interaction of multiple non-planar

generate complex fracture networks thought to be beneficial for production [2, 3].

dissipate.

662 Effective and Sustainable Hydraulic Fracturing

interaction and design [14-18].

be described here.

The complete theory of fracture interaction for an arbitrary number of fractures is quite cumbersome. Therefore, a 'simplified' analytical model has been developed to describe a set of fractures growing in a viscous mass-transfer dominated regime. In this model, two types of fracture growth patterns are defined. First is the 'compact' fracture growth where interaction and interference between fractures are minimum. The second is 'diffuse' fracture growth where some fractures will either terminate or collapse into each other because of minimum insitu stress rotation (Figure 1). In diffuse growth, we distinguish different zones. The first zone is the closest to the well. The next zone, farther from the well, contain less fractures because some of them have either terminated or merged into other fractures, as shown at the righthand side of Figure 1.

The transition from one growth pattern to another is delineated by a 'critical' fracture spacing, influenced by the original principal stress difference in the un-fractured forma‐ tion. In the case of diffuse growth, different zones are depicted by the dashed lines on the right-hand picture of Figure 1. The Average permeability of the fracture set in the zone "n" is given by (see e.g. [20]):

$$k\_n = \frac{w\_n^{-3}}{12 \, L\_n} \tag{1}$$

where *Ln* is spacing between the fractures in the "n" zone, and average aperture of fractures in the zone "n" is:

$$aw\_n = a \frac{p\_n}{E} R\_n \tag{2}$$

Here, *pn* is net pressure, *E* is Young's modulus, *Rn* is a radius of a penny-shape fracture, and *a* is a numerical coefficient (of approximately 1) and depends on Poisson's ratio. For simplicity, it is assumed that there is no leak-off of fluid into the extremely tight formation. Therefore, the average fluid density, *m*, in the fractures becomes:

**Figure 1.** "Compact" and "Diffuse" fracture growth patterns

$$m = \frac{w\_n \, \rho}{L\_n} \tag{3}$$

The above Eq. 6 is highly nonlinear and cannot be solved analytically, instead, it is solved by considering the integral relationships and using a modification of the so-called method of successive change of steady states [21]. The essential idea is to solve the steady-state version of Eq. 6 for a given size of the fracture (i.e. without the first term in Eq. 6). Subsequently, the time dependency of fracture radius growth is found from the integral form of mass balance equation, the differential form of which is given by Eqs. 4 and 6. The detail of this solution will

The Geomechanical Interaction of Multiple Hydraulic Fractures in Horizontal Wells

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

665

Interaction between fractures is governed by the induced stress due to the growing fractures. Solutions for various cases of interacting fractures have been published [22]. Since these solutions are cumbersome, several approximations have been proposed [23], which make the problem more tractable for our present purpose. Following these approximations we can obtain the stress state induced by two neighbouring growing fractures. The induced stress depends on the fracture aperture (or equivalently, net pressure) and a dimensionless ratio of

To appreciate the geomechanical interaction of fractures, see the left side of Figure 2, which depicts a top view of two simultaneously growing fractures. We choose to examine the incremental stress changes in the x and y directions at a mid point in between the fractures, as shown in the figure. The induced stress is obviously not constant in the space between the fractures. But for illustration purpose, we will simply consider the stress evolution at this chosen point to be representative of the stress state of the large area between the fractures. The analytical model on the right side of Figure 2 shows that σ<sup>x</sup> increases with fracture length. Although σy also increases, it does so at a much slower pace than σx. We can now define a deviatoric stress, Δσ = σx-σy, which exhibits a maximum value. It can be immediately observed that if the maximum value of Δσ is greater than the original in-situ stress difference (σH-σh), then rotation of minimum in-situ stress direction is feasible and the fracture(s) would grow in a non-planar fashion. For multiple fracture growth, this will result in the collapse or termina‐

In the same way, the height growth of multiple fractures may also see similar interaction if Δσ is greater than the original in-situ stress difference between vertical stress and minimum horizontal stress (σV-σh). As fractures grow upward, assuming no stress barrier or contrast is encountered, both σV and σh would typically decrease in magnitude. Potentially, the stress shadowing effect is even more prominent for height growth of multiple fractures. It is interesting to note that field experience does indicate that it is challenging to grow in fracture height for multiple fracture stimulation. A plausible explanation is indeed the 'stress shadow'

Although the analytical solution is an approximation, it captures the foremost qualitative features of fracture interaction. One of the first results shows how the in-situ stress contrast influences the domains of fracture growth. Figure 3 shows the results based on the analytical model, where a distinct transition between the two regimes of fracture growth can be identi‐ fied. Note that in practice, typical fracture spacing is in the order of 70 to 100ft (20-30m). Therefore, according to this simplistic, homogeneous model, today's multiple hydraulic fracture stimulation practises can potentially experience 'strong interference', depending on

effect caused by geomechanical interaction of multiple fractures.

fracture length and the distance between fractures (i.e. the initial fracture spacing).

be given in a future paper.

tion of some fractures.

the stress contrast.

Where, *ρ* is fluid density. Applying the equation of fluid mass balance:

$$\frac{\partial m}{\partial t} + \operatorname{div} \rho \stackrel{\rightarrow}{\boldsymbol{v}} \boldsymbol{v} = \mathbf{0} \tag{4}$$

and Darcy's law:

$$\stackrel{\rightarrow}{w}\_{n} = -\frac{w\_{n}^{3}}{12\,\mu\text{L}\_{n}} \cdot \text{grad}(p\_{n})\tag{5}$$

We obtain an equation for fluid transport in the fractures:

$$\frac{\partial}{\partial t}(\frac{w\_n}{L\_n}) - \frac{E}{12\,\mu\,a}div \{\frac{w\_n^3}{L\_n} \cdot \text{grad}(\frac{w\_n}{R\_n})\} = 0\tag{6}$$

The above Eq. 6 is highly nonlinear and cannot be solved analytically, instead, it is solved by considering the integral relationships and using a modification of the so-called method of successive change of steady states [21]. The essential idea is to solve the steady-state version of Eq. 6 for a given size of the fracture (i.e. without the first term in Eq. 6). Subsequently, the time dependency of fracture radius growth is found from the integral form of mass balance equation, the differential form of which is given by Eqs. 4 and 6. The detail of this solution will be given in a future paper.

Interaction between fractures is governed by the induced stress due to the growing fractures. Solutions for various cases of interacting fractures have been published [22]. Since these solutions are cumbersome, several approximations have been proposed [23], which make the problem more tractable for our present purpose. Following these approximations we can obtain the stress state induced by two neighbouring growing fractures. The induced stress depends on the fracture aperture (or equivalently, net pressure) and a dimensionless ratio of fracture length and the distance between fractures (i.e. the initial fracture spacing).

To appreciate the geomechanical interaction of fractures, see the left side of Figure 2, which depicts a top view of two simultaneously growing fractures. We choose to examine the incremental stress changes in the x and y directions at a mid point in between the fractures, as shown in the figure. The induced stress is obviously not constant in the space between the fractures. But for illustration purpose, we will simply consider the stress evolution at this chosen point to be representative of the stress state of the large area between the fractures. The analytical model on the right side of Figure 2 shows that σ<sup>x</sup> increases with fracture length. Although σy also increases, it does so at a much slower pace than σx. We can now define a deviatoric stress, Δσ = σx-σy, which exhibits a maximum value. It can be immediately observed that if the maximum value of Δσ is greater than the original in-situ stress difference (σH-σh), then rotation of minimum in-situ stress direction is feasible and the fracture(s) would grow in a non-planar fashion. For multiple fracture growth, this will result in the collapse or termina‐ tion of some fractures.

r = *<sup>n</sup> n*

*<sup>L</sup>* (3)

*<sup>t</sup>* (4)

*<sup>v</sup> grad p <sup>L</sup>* (5)

*w ww <sup>E</sup> div grad tL a L R* (6)

*w m*

r

<sup>0</sup> *<sup>m</sup> div v*

® ¶ + =

m

¶ - ×= ¶ 3 ( ) { ( )} 0 <sup>12</sup> *n nn n nn*

( ) <sup>12</sup> *n n n n*

=- × 3

m

*w*

Where, *ρ* is fluid density. Applying the equation of fluid mass balance:

**Figure 1.** "Compact" and "Diffuse" fracture growth patterns

664 Effective and Sustainable Hydraulic Fracturing

¶

®

We obtain an equation for fluid transport in the fractures:

and Darcy's law:

In the same way, the height growth of multiple fractures may also see similar interaction if Δσ is greater than the original in-situ stress difference between vertical stress and minimum horizontal stress (σV-σh). As fractures grow upward, assuming no stress barrier or contrast is encountered, both σV and σh would typically decrease in magnitude. Potentially, the stress shadowing effect is even more prominent for height growth of multiple fractures. It is interesting to note that field experience does indicate that it is challenging to grow in fracture height for multiple fracture stimulation. A plausible explanation is indeed the 'stress shadow' effect caused by geomechanical interaction of multiple fractures.

Although the analytical solution is an approximation, it captures the foremost qualitative features of fracture interaction. One of the first results shows how the in-situ stress contrast influences the domains of fracture growth. Figure 3 shows the results based on the analytical model, where a distinct transition between the two regimes of fracture growth can be identi‐ fied. Note that in practice, typical fracture spacing is in the order of 70 to 100ft (20-30m). Therefore, according to this simplistic, homogeneous model, today's multiple hydraulic fracture stimulation practises can potentially experience 'strong interference', depending on the stress contrast.

**3. Non-planar 3D model of hydraulic fractures**

power-law fluid commonly characterized by parameters *K* and *n*:

the viscosity is low).

employed time step is generally small.

**3.1. Non-planar fracture propagation**

t

The numerical simulator essentially follows a well accepted numerical methodology [24, 25]. In this numerical framework, the rock mechanics of fracture growth takes place in a linear elastic medium and therefore can be modelled by boundary integral equations based on fundamental solutions of dislocation segments. The fracture growth is controlled by both fracture mechanics criteria and fluid volume balance at the fracture tip. The formulation adequately portrays both viscosity and toughness controlled fracture growth as well as the 'fluid lag' behaviour [25]. The proppant laden frac fluid is idealized as an incompressible

*n*

Where *τ* is the fluid shear stress and *γ*˙ the fluid shearing rate. *K* and *n* are termed the consistency index and power law index, respectively. They are naturally dependent on the proppant concentration, and this effect of proppant concentration on fluid rheology is accounted for in the simulator using Shah's model [27]. The proppant transport and settlement are computed via the transport equation and empirical equations of particles settling between parallel plates

The ensuing flow of proppant laden frac fluid within a fracture is captured according to the Reynolds lubrication theory, and the derived non-linear fracture growth and fluid flow equations are now solved in a coupled manner in which mass conservation (i.e., frac fluid and proppant) are strictly enforced. The equations are discretized on the same structured grid, and they form a system of moving boundary, transient coupled equations of fracture width and fluid pressure. A new adaptive time integration algorithm has been developed to provide numerically stable solutions for a wide range of power-law fluid properties (especially when

Although the proppant transport is calculated at each time step based on a volume conserva‐ tion equation, this computation is decoupled from the process of solving the fluid flow equation for the sake of computational efficiency. Specifically, the fluid pressure and fracture width are first obtained by 'freezing' the proppant concentration associated with each element. The contribution of proppant transport to the volume concentration is then updated based on the calculated fluid velocity field. This proves to be effective as well as efficient since the

In dealing with non-planar fracture growth, the heterogeneous formation is modelled as multiple isotropic parallel layers with heterogeneities characterized by variations in in-situ stress, elastic stiffness modulus, Poisson's ratio, fracture toughness, pore pressure, leak-off coefficients, and spurt-loss coefficients. All fractures are assumed to be vertical but they can turn in any horizontal direction. For each fracture, the growth direction is determined by the

= *K* & (7)

The Geomechanical Interaction of Multiple Hydraulic Fractures in Horizontal Wells

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

667

 g

[28]. The fluid loss to the formation is described by a Carter type leak-off model [25].

**Figure 2.** Interaction between two fractures

**Figure 3.** Fracture spacing versus stress contrast
