The Geomechanical Interaction of Multiple Hydraulic Fractures in Horizontal Wells

The technology of multiple hydraulic fracture stimulation in horizontal wells has transformed the business of oil and gas exploitation from extremely tight, unconventional hydrocarbon bearing rock formations. The fracture stimulation process typically involves placing multiple fractures stage-by-stage along the horizontal well using diverse well completion technologies. The effective design of such massive fracture stimulation requires an understanding of how multiple hydraulic fractures would grow and interact with each other in heterogeneous formations. This is especially challenging as the interaction of these fractures are subject to the dynamic process of subsurface geomechanical stress changes induced by the fracture treatment itself.

This paper consists of two parts. Firstly, an idealised analytical model is used to highlight some key features of multiple hydraulic fractures interaction, and to provide a quantification of ‘stress shadow’. Secondly, a new non-planar three dimensional (3D) hydraulic fracturing numerical model is used to provide an insight into the growth of multiple fractures under the influence of subsurface geomechanical stress shadows. Attention is given to studying the height growth of multiple fractures.

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

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 developed 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 non-planar 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

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 in-situ 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-hand side of Figure 1.

Figure 1.

“Compact” and “Diffuse” fracture growth patterns

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 formation. 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}{12L_n}$$E1

where L_n is spacing between the fractures in the “n” zone, and average aperture of fractures in the zone “n” is:

$$w_n=a\frac{p_n}{E}{R}_n$$E2

Here, p_n is net pressure, E is Young’s modulus, R_n 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:

$$m=\frac{w_n\rho }{L_n}$$E3

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

$$\frac{\partial m}{\partial t}+div\rho \overrightarrow{v}=0$$E4

and Darcy’s law:

$$\overrightarrow{v_n}=-\frac{w_n^3}{12\mu L_n}\cdot grad(p_n)$$E5

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 grad(\frac{w_n}{R_n})\}=0$$E6

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 σ_x 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 termination of some fractures.

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

Figure 2.

Interaction between two fractures

Figure 3.

Fracture spacing versus stress contrast

3. Non-planar 3D model of hydraulic fractures

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 power-law fluid commonly characterized by parameters K and n:

$$\tau =K{\dot{\gamma }}^n$$E7

Where τ is the fluid shear stress and $\dot{\gamma }$ 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 [28]. The fluid loss to the formation is described by a Carter type leak-off model [25].

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 the viscosity is low).

Although the proppant transport is calculated at each time step based on a volume conservation 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 employed time step is generally small.

3.1. Non-planar fracture propagation

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 stress state of the leading element at the fracture tip. These assumptions exclude the application for certain geological formations. In particular, geological stress regimes as a result of reverse faulting may present difficulties in modelling as the hydraulic fracture would most likely grow in a horizontal plane. However, the assumptions made are sufficient and adequate for typical deep unconventional resource formations where the vertical stress is generally the maximum principal stress, or at least greater than the minimum horizontal principal stress – resulting in vertical fracture growth.

Fractures advance when the maximum tensile stress ahead of the crack tip exceeds the intrinsic tensile strength of the rock. In linear elastic fracture mechanics, the stress at the crack tip is singular, therefore the stress intensity factor, which correlates to the strength of such a singular stress field, is generally used to determine fracture tip propagation [29]. In practice, hydraulic fracturing processes take place under relatively large compressive in-situ earth stress that is normally several orders of magnitude higher than the rock strength. Therefore, the process is dominated by the balance of the compressive stress and fluid pressure at the fracture tip region. A group of ‘virtual’ elements is placed along the crack front and at each time step a check is made on the stress status of these elements. The virtual elements are allowed to be active as part of the new fracture surface if the potential fluid pressure can overcome the compressive stress plus the strength of the rock.

Importantly, the interaction of multiple fractures may result in significant shearing displacements along the fracture surface, causing the fracture growth to follow a curved path, and the computational technique leads these virtual elements to be oriented according to the stress field so that the local minimum principal stress is perpendicular to the new fracture surface as the fracture front advances.

An outline of the mathematical description of the fracture stiffness matrix, the displacement discontinuity method is given in [19], and will not be repeated here. The numerical approach is robust and efficient but less accurate than the variational boundary integral method which employs second order elements [30, 31]. For the current application, computational efficiency is deemed to outweigh the importance of a moderate improvement in accuracy. [19] also describes the coupling of derived non-linear fracture growth to flow equation and mass conservation.

3.2. Wellbore hydraulic model

The fractures are assumed to be initiated at the specified locations of a cluster of perforations. Near wellbore formation stress concentration is not considered. The downhole pressure is defined as the fluid pressure just upstream of all injection points. It is obvious that the distribution of the corresponding injection at each fracture follows wellbore fluid mechanics, depending on the injection area connected with the fracture and the fluid pressure within the fracture. An approximation is made by using Bernoulli equation to capture the fluid distribution. This wellbore hydraulic model is able to account directly for the empirically derived or calibrated frictional pressure drop at the perforation and along the wellbore. Therefore, it is possible to account for limited entry perforation designs commonly employed in multiple hydraulic fracture design.

4. Comparing simulation results with analytical solution and experiment

The complete validation of a hydraulic fracturing simulator is enormously challenging, because of the difficulty in obtaining the general analytical solutions or constructing realistic laboratory experiments. A carefully conducted laboratory hydraulic fracturing experiment has been published [32]. Figure 4 shows our simulation of the experiment. The comparison is excellent not only in the created fracture geometry, but also the pumping parameters.

Figure 4.

Comparing with experimental result of [32]

5. Frac fluid viscosity and fractures height growth

A single stage of four hydraulic fractures is considered. Parameters chosen broadly resemble a typical shale gas development field case. The horizontal wellbore is placed at the relative depth of 0m as shown in Figure 5. Fracture injection points are spaced at 25m (82ft) apart and each fracture has 12 perforations. Fluid is injected into the wellbore at constant injection rate of 0.2 m³/s (~75 bbl/min). Formation heterogeneities are limited to only in-situ stress variation. In this instance, the minimum horizontal stress gradient is set to 14kPa/m (0.62psi/ft) with no upper stress barrier to restrict height growth. The difference between minimum and maximum horizontal stress is chosen to be more than 20%, so that we do not cause the fracture direction to re-orient drastically. The key parameters are shown in Figure 5. The fluid leak-off factor is set to 1 with no spurt loss.

We consider two cases of distinctly different frac fluids, one with low viscosity (slick water) and the other a very high viscosity (gel frac). This corresponds to setting n = 0.8 and K = 0.05 for the low viscosity fluid case, and n = 0.8 and K = 0.5 for the high viscosity fluid in the power-law fluid model of Equation 7.

The injection parameters and the calculated total fracture volumes are shown in Figure 6. The simulated results are shown in Figure 7. As expected, the low viscosity fluid creates larger fractures with more fracture surface area at the expense of fracture aperture/width. Significantly, the low viscosity fluid creates a much larger fracture height. It is now interesting to see how the multiple fractures evolve in length and height as they compete in the presence of stress perturbation. Figure 8 shows fracture development at the end of the pumping of high viscosity fluid. The outer two fractures grew slightly outward and are longer than the inner fractures. This is due to stress shadowing. Interestingly, the inner two fractures developed more height, and this behaviour will become even more pronounced when we look at low viscosity fluid.

Figure 5.

Formation input parameters for the case studies

Figure 6.

Injection volume and calculated fracture volumes

Figure 7.

Calculated fracture parameters

Figure 8.

Fracture geometry with high viscosity frac fluid at 150 mins

Figure 9.

Fracture geometry of low viscosity frac fluid at 150 mins

Figure 9 shows the final fracture geometry of pumping low viscosity fluid (note that the scale in Figure 9 is different from that in Figure 8). Perhaps against intuition, the inner fractures grew appreciably in height. To understand this, we look at the evolution of fracture length and height growth over time as depicted in Figures 10 and 11. It becomes clear that as the inner fractures are constrained from growing in length because of competition with the outer fractures, they found relative freedom to grow upward instead.

Figure 10.

Fracture geometry of low viscosity frac fluid after 20 mins

Figure 11.

Fracture geometry of low viscosity frac fluid after 40 mins

6. Discussion

At present, multiple fracture stimulations in horizontal wells often assume a uniform growth of all fractures, with the total volume of frac fluid more or less evenly distributed amongst the fractures. We have pointed out some of the potential impacts of stress shadows, which can affect treatment cost and production. Crucially, it has consequences in the understanding of subsurface development. For instance, if, instead of stimulating four fractures uniformly, only two or three of the fractures grow disproportionally, then some of the fractures may be unknowingly over-stimulated, resulting in excessive length or height growth. Normally this does not pose any major issue unless there are near-by wells or geological faults.

The non-planar 3D simulator presented here is capable of capturing the influence of key parameters such as injection rate, fracture spacing, formation properties, etc. on fracture design. Attempts are being made to extend the present numerical framework to include the interaction of hydraulic fractures with natural fractures.

Acknowledgments

The authors recognize the valuable input from Alexei Savitski, Anastasia Dobroskok, Mauricio Farinas, Ernesto Fonseca, Kaiming Xia, and Zongyu Zhai. Permission by Shell International Exploration and Production Inc. to publish is gratefully acknowledged.