Morphology of multiple constant height hydraulic fractures versus propagation regime

Multistage hydraulic fracturing is often used to enhance recovery of hydrocarbons from low permeability formations [1]. To optimize the process, multiple hydraulic fractures are initiated simultaneously from a horizontal wellbore. In addition, this is usually supplemented by a limited entry design, which promotes uniform flow distribution between the fractures. It was shown, however, that even in the simplest case of no geological layers, such a problem can lead to the development of instabilities that can cause complex fracture shapes [2]. In addition, similar behavior was observed for field scale simulations of hydraulic fracturing [3]. This behavior is important from both fundamental and practical points of view. First, it allows to better understand processes happening underground, and second, it permits to control this behavior using operational procedures.

As was shown in [2], hydraulic fracture propagation regime is the key factor that determines morphology of simultaneously propagating hydraulic fractures. For simple fracture geometries, such as plane strain, radial, and constant height cracks, there are four limiting cases or propagation regimes [4]. These are determined by the competition between viscous dissipation and fracture toughness, as well as the fluid storage inside the fracture or the surrounding formation. It is important to note that these mathematical conclusions are based on the model that features linear elastic rock deformation, Newtonian fracturing fluid, linear elastic fracture mechanics, Carter’s leak-off, and uniform injection rate. Other important physical processes, such as proppant transport or the effect of hydrostatic pressure are ignored. Nevertheless, it was concluded in [2] that when viscosity dominates, then all the fractures are nearly identical. At the same time, when toughness dominates, then the fractures strongly compete with each other and develop petal-like shapes. This was the case for “radial” geometry, or when there were no geological layers. Qualitatively similar observations were made for field scale hydraulic fracture modeling in [3]. Namely, the increase of fracture toughness leads to more fracture asymmetry, while the increase of viscosity (and/or decrease of toughness) leads to more regular and symmetric fractures. One of the challenges, however, is to quantitatively determine the parameter that governs such a transition in fracture morphology. As a result, this study aims to extend the previous analysis [2] for the case of constant height hydraulic fractures. While such a geometric assumption is still far from a general field scale simulation, it is arguably a better approximation relative to the radial case, especially when strong barriers limit fracture propagation in the vertical direction. Finally, and most importantly, it is possible to perform a relatively simple analysis for the constant height hydraulic fractures, which is an essential component that is needed to quantify the fracture morphology.

Constant height hydraulic fractures are often referred to as Perkins–Kern–Nordgren (PKN) [5, 6] fractures. Many researchers have investigated propagation of PKN fractures from various perspectives. A few papers that are relevant to this work are [7, 8, 9]. The classical PKN model does not account for the effect of fracture toughness and therefore is unsuitable for analyzing the problem under consideration since the transition from toughness domination to viscosity domination presents the most interest. The first attempt to include the effect of toughness was done in [7], where the propagation condition in terms of pressure was taken from the uniformly pressurized radial fracture with the diameter equal to the reservoir height. This methodology was then revised in [8], in which the authors proposed a more accurate procedure that utilizes energy considerations. Finally, it was shown in [9] that the second energetic approach is indeed more accurate (albeit quite marginally) by comparing the results to the solution with non-local elasticity as well as to the simulations performed with a fully planar model.

One prerequisite for conducting the study on simultaneous propagation of multiple hydraulic fractures is the knowledge of the parametric space for the problem. For instance, for the case of radial geometry [2], the prerequisite paper is [10], which builds on top of many previous publications and outlines boundaries in the parametric space that quantify dominance of various physical processes. In the context of constant height or PKN fractures, such a parametric space was recently constructed in [11]. In particular, the model that accounts for toughness based on the results in [8] is used due to its simplicity and still admissible accuracy. Section 2 summarizes mathematical formulation for the model and describes the aforementioned parametric space. Equipped with this parametric space, several cases are constructed for multiple fractures that lie within different regions of the parametric space. Numerical simulations are performed using commercial code ResFrac¹¹1www.resfrac.com [12] and results of simulations for simultaneous propagation of multiple hydraulic fractures are presented in section 3 for different regimes. Finally, conclusions are summarized in section 4.

This section outlines the mathematical model for a single PKN fracture and also presents the parametric space for the problem. With the reference to Fig. 1, let the fracture be vertical and occupy the $(x,y)$ plane, where $y$ is the vertical coordinate and $x$ is the horizontal coordinate. The fracture height is denoted by $H$, while the half-length is $l(t)$. The model assumes that $H\ll l$, which ensures that the flow is predominantly horizontal and also confirms validity of the local elasticity assumption.

By following the classical PKN model [5, 6], fracture width in each vertical fracture cross-section is assumed to be elliptical, and the fracturing fluid pressure is determined based on the local plane strain elasticity assumption:

Here $w(x,y)$ is the fracture opening, $\bar{w}(x)$ is the effective or averaged width, $p(x)$ is the fluid pressure, and $E^{\prime}=E/(1\!-\!\nu^{2})$ is the plane strain Young’s modulus. The primary governing equation is the vertically-averaged volume balance, which reads:

In the above equation $\mu$ is the fluid viscosity, $\bar{q}_{x}$ denotes the horizontal flux, $C^{\prime}=2C_{l}$ is the scaled Carter’s leak-off coefficient, and $t_{0}(x)$ corresponds to time when the fracture front was located at the point $x$.

Propagation condition at the lateral fracture tip is based on the results from [8] and can be written in terms of the toughness dependent pressure condition as

where $K_{Ic}$ is fracture toughness. This propagation condition is supplemented by the zero flux at the tip, which translates into the following equation for fracture length:

By following the results in [11], one possibility to reduce the number of parameters in the problem is to introduce the following normalization:

Here the width, length, and time scales are given by

By adopting such a normalization, the governing equations (1)–(4) become

where the dimensionless parameter that quantifies leak-off is

Thus, the solution to the problem in terms of length and effective width depends on two dimensionless parameters, namely dimensionless time $\tau$ and dimensionless leak-off $\phi$. Conceptually, these parameters are similar to those for a radial fracture [10], but quantitatively they are very different.

Parametric space for this problem was constructed in [11] and is shown in Fig. 2. The blue, red, green, and magenta regions indicate zones of applicability of storage viscosity or $M$, storage toughness or $K$, leak-off viscosity or $\tilde{M}$, and leak-off toughness or $\tilde{K}$ solutions, respectively. These are the zones, in which the difference between the “true” solution and the limiting self-similar solution falls below a small threshold, see [11] for more details. Here the storage viscosity limit corresponds to dominance of viscosity over toughness and small leak-off, i.e. when most of the injected fluid stays inside the fracture. In contrast, storage toughness limit occurs when toughness dominates over the viscosity, while most of the injected fluid still stays inside the fracture. The leak-off regimes are defined in a similar fashion, but most of the injected fluid leaks into the formation in these limits. Equations that define boundaries of the transition regions are given by

In the context of the analysis of propagation of multiple hydraulic fractures, this parametric space allows to select problem parameters that correspond to dominance of either toughness, viscosity, or fall into the situation when both are important.

A series of numerical simulations for simultaneous propagation of multiple hydraulic fractures is performed in ResFrac, which is a fully coupled 3D hydraulic fracturing and reservoir simulator [12]. Three basic cases are considered: toughness domination (denoted by $K$), transition (denoted by $MK$), and viscosity (denoted by $M$) domination. The input parameters as well as the corresponding dimensionless parameters $\tau$ and $\phi$ are summarized in Tab. 1, see (5), (6), and (8) for the definitions of $\tau$ and $\phi$. Note that leak-off cases are not considered because previous study [2] showed that leak-off does not qualitatively affect the behavior. Uniform flux distribution is achieved by substantially increasing perforation friction, the effect of gravity is neglected, and the height constraint is introduced by having strong barriers above and below the primary layer of interest. All the injection points are located in the middle of the layer. Only 5 fractures are initiated, all of which are planar and perpendicular to the well in the base case. Also, the base simulations have 10 m uniform spacing between the fractures, while other values are used to investigate the effect of spacing, namely 20, 40, 80, and 160 m.

Fig. 3 shows results of numerical simulations for 10 m spacing. First of all, the parametric map (similar to the one shown in Fig. 2) is shown and three circular markers indicate relative location of the three cases under consideration. As can be seen, one point lies within the toughness domination zone, the second one is in between the toughness and viscosity domination zones, while the third point lies within the viscosity domination zone. Hydraulic fractures for all these three cases are shown in two views: original, in which physical proportions are preserved, and stretched, in which the result is stretched 10 times along the wellbore for better fracture visualization.

As can be seen from the fracture geometries in Fig. 3, toughness domination (case 1) leads to unstable fracture growth, whereby the fractures tend to avoid each other and form irregular shapes. Some fractures propagate predominantly to the left from the well and some propagate to the right. Moreover, the “local” fracture height of each individual fracture is not equal to the reservoir height (or height of the primary layer). Once the viscosity is increased and toughness is reduced (case 2), the fractures become more symmetric and have nearly identical lengths, even though the “local” fracture height is still variable along each individual fracture. Finally, when viscosity dominates (case 3), all the fractures become nearly identical and the fracture height of each fracture is equal to the reservoir height in most regions of the fracture. The exception is the near tip region, whose size is on the order of fracture height $O(H)$, where the non-uniform behavior still persists. Thus, these results demonstrate that the behavior is qualitatively similar to that observed for “radial” geometry in [2], i.e. the case without the barriers restricting vertical growth. One very big difference, however, is the parameter that quantifies the transition from one limiting case to another. For the radial case, the transition from the viscosity to the toughness regime is determined by the parameter

where the lower bound corresponds to the onset of the viscosity regime, while the upper bound corresponds to the onset of the toughness regime. For constant height fractures, on the other hand, the result is

where the lower bound corresponds to the onset of the toughness regime, while the upper bound corresponds to the onset of the viscosity regime. Clearly, these equations (10) and (11) are very different and lead to very different results for these geometries.

Fig. 4 illustrates sensitivity of the results to fracture spacing. Numerical simulations are performed for the base simulations (see Tab. 1), but for various spacing ranging from 10 m to 160 m. Given that the reservoir height is 80 m, then the ratios between the spacing and height are from 1/8 to 2/1. The same number of fractures, i.e. 5, is used in all simulations. The vertical axis corresponds to the approximate value of $\log_{10}(\tau)$, but essentially shows the results for toughness (lowest $\tau$), transition (intermediate $\tau$), and viscosity (highest $\tau$) cases. As can be seen from the results, the normalized spacing significantly affects the fracture morphology, as expected. Once the spacing is twice as much as the fracture height, there is almost no interaction between the fractures, which is consistent with the fact that the stress shadow decays quickly at distances commensurate with fracture height (for constant height fractures). Once the fractures are brought closer together, stress interaction starts to distort the fracture shapes if toughness dominates or sufficiently large. Viscosity dominated fractures, on the other hand, remain very similar. These results illustrate how fracture height and spacing, as well as the fracture regime affects the resulting fracture morphology.

Finally, Fig. 5 presents sensitivity of the results to stress orientation. In particular, the change of the stress orientation angle allows to effectively adjust the angle between the well and the direction of fracture propagation. The zero angle corresponds to fractures propagating perpendicular to the well, while other cases correspond to the rotation (measured in degrees) relative to this baseline position. Note that the fracture spacing along the well is kept fixed, so that the physical spacing between the fractures, measured perpendicularly to the fracture surface, is effectively multiplied by the cosine of the rotation angle. The input parameters for the “base” simulation, i.e. the one corresponding to zero degree stress rotation, are given in Tab. 1. Results shown in Fig. 5 demonstrate that changing the stress orientation angle noticeably affects the results. First of all, once the symmetry is broken, the fractures that have an initial lateral shift tend to propagate predominantly in the direction of this shift. This applies to all cases, albeit to a different degree. The viscosity dominated cases still have significant fracture overlap, and this overlap drastically reduces towards lower values of $\tau$ or toughness domination. Another interesting observation is that the height of each individual fracture tends to be equal to the reservoir height for larger stress angles, even in the toughness domination regime. That is, the complexity of fracture morphology is reduced once the stress orientation angle increases. This is caused perhaps by the initial bias that fractures have, which significantly affects the stress interaction at early times of propagation. Note that such a behavior persists even if the fluid viscosity is further reduced, i.e. when the system is brought further inside the toughness domination zone.

In this work, hydraulic fracture morphology for the case of simultaneous propagation of several fractures is investigated in view of dominance of either fluid viscosity or toughness. In particular, the results are obtained for the case of height contained fractures. First, the dimensionless parameters that define the parametric space are outlined and a single parameter that determines the transition from toughness to viscosity domination is defined. Numerical results for 3 base cases (toughness domination, viscosity domination, and transition) are performed for the case of 5 uniformly spaced fractures. It is observed that toughness domination leads to strong stress interaction between the fractures, which leads to the creation of complex fracture geometries. At the same time, relatively large viscosity tends to regularize the fracture morphology and leads to the formation of nearly identical fractures. Numerical simulations with different fracture spacings indicate that the ratio between the fracture spacing to the reservoir height plays a crucial role, since the stress interaction decays quickly as the spacing to height ratio is increased. Also, the influence of stress orientation on the results is investigated, which showed that it also affects the fracture morphology, but to a smaller degree. Since the results are presented with the reference to the dimensionless parametric space, these findings can help to quickly identify the fracture regime and morphology for any input parameters and can be used as a guideline to reach a desired fracture morphology configuration.

The author would like to acknowledge Mark McClure for useful recommendations regarding this work.