Numerical modelling and diagnostic techniques of hydraulic fractures based on their inlet behaviour

Modelamiento numérico y técnicas de diagnóstico de fracturas hidráulicas basados en su comportamiento en el punto de inyección

 

Cristóbal Valderrama and Esteban Sáez

Departamento de Ingeniería Estructural y Geotécnica, P. Universidad Católica de Chile, Avda. Vicuña Mackenna 4860, Macul, Santiago, Chile, cevalder@uc.cl, esaez@ing.puc.cl

---

In the current conditions of mines exploited by caving methods, pre-conditioning by hydraulic fracturing has proven to have positive impacts such as the decrease in the block’s size related to the primary fragmentation. An essential part of the pre-conditioning design is the estimation of the hydraulic fractures length. On the other hand, the energy dissipated by the viscous flow imide the fracture modifies the propagation characteristics of these fractures, making difficult the use of standard methods of fracture mechanics in numerical modelling. For this reason, for plain-strain and axisymmetric cases we propose a numerical resolution strategy, which can be used for any set of hydraulic fracturing parameters. These criteria are based on the pressure and opening values at the inlet, and additionally on their length (or radius). Furthermore, this last characteristic, allow us to modify the propagation criteria in order to generate diagnostic tools for the estimation of fractures dimensions in the field.

Keywords: hydraulic fracturing, numerical modelling, diagnostic technique, mine pre-conditioning

---

En las actuales condiciones de las minas explotadas por métodos de hundimiento, el pre-acondicionamiento con fracturamiento hidráulico ha probado tener impactos positivos como la disminución del tamaño de bloques asociados a la fragmentación primaria. Parte fundamental del diseño del preacondicionamiento es la estimación de la longitud de las fracturas hidráulicas. Por otro lado, la energía disipada por el flujo dentro de la fractura modifica las características de propagación de estas fracturas, haciendo difícil el uso de los métodos estándar de la mecánica de fracturas en la modelación numérica. Por esta razón, para los casos de deformaciones planas y de simetría de revolución se presentan estrategias de resolución numérica, que pueden ser utilizadas para cualquier conjunto de parámetros del fracturamiento hidráulico. Estos criterios utilizan los valores de la presión o de la apertura de la fractura en el punto de inyección, además de su longitud (o radio). Además, esta última característica permite modificar los criterios de propagación generando herramientas de diagnóstico para la estimación de las dimensiones de las fracturas hidráulicas en terreno.

Palabras clave: fracturamiento hidráulico, modelamiento numérico, técnica de diagnóstico, pre-acondicionamiento de minas

---

 

Introduction

Hydraulic fracturing is the process of fracture inducement and propagation, which exist in natural conditions, as intrusive dykes generated by magma driven flows (Spence and Turcotte, 1985); and artificially developed for industrial applications such as: increase of the permeability in oil and gas reservoirs (Valkó and Economides, 1995), stimulation of geothermal reservoirs, nuclear waste disposal (Abou-Sayed, 1994), in-situ stress measurements (Haimson and Cornet, 2003) and rock mass pre-conditioning in underground mines exploited by non-supported methods such as caving and longwall methods (Jeffrey et al., 2001; van As et al., 2004).

On the other hand, the current trend in mining is an increase in the importance of underground mines, because reserves are deeper and with lower grades (Moss, 2011), where a large number of caving mines are located in hard rock masses, with a high-stress environment and low density of discontinuities. These conditions involve an increment of seismicity, stagnation of caving and increase in the blocks size. To overcome these problems, pre-conditioning by hydraulic fracturing is being used to generate a controlled alteration of the rock mass previous to mining, which has induced a decrease in the magnitude of seismicity, an increase of the connection time to the crater and a reduction in the size of rock fragments (Araneda and Sougarret, 2007).

A proper estimation of fracture’s length evolution is an important part of the design of hydraulic fracturing, which is generally done by means of numerical modelling; hence, it is crucial to develop mathematical models able to predict the fracture growth (length, opening and internal fluid pressure variation). When the rock mass is impermeable, the problem involves three coupled phenomena: the mechanical deformation of fracture’s surfaces, the fluid flow inside the fracture and fracture’s propagation. From the point of view of the problem’s mathematical formulation, coupling of the aforementioned phenomena leads to a complex structure solution which is dominated by near tip processes (Adachi et al., 2007; Detournay, 2004). For the case of hydraulic fractures that propagates in a medium without fracture resistance, it has been shown that a stress singularity proportional to r^-1/3 exists (Desroches et al., 1994), where r is the coordinate ahead of the tip, which is different to the r^-1/2 classical singularity of Linear Elastic Fracture Mechanics LEFM. Furthermore, fracture toughness K_Ic is not necessarily a leading term in the solution. These facts complicate (and, in some cases prevent) the use of typical numerical methods to control fracture propagation; however, the problem has been solved imposing the propagation condition K_I = K_Ic (Akulich and Zvyagin, 2008; Desroches and Thiercelin, 1993), where the stress intensity factor K_I is calculated numerically using the weighted function method. Also, the imposition of proper asymptotic behaviour of tip openings has been used to solve some specific cases (Garagash and Detournay, 2007; Garagash, 2006). Nevertheless, these approaches are based on the complex tip behaviour and explicit energy concepts have not been used.

After hydraulic fracture is generated, an evaluation of the performance is necessary, which can be accomplished through indirect methods, such as well testing and net-pressure analysis, and by means of direct methods, such as tiltmeters and micro-seismic mapping (Moss and Maron, 1987). Normally, direct methods give better results than indirect ones, but they require specific equipment, and additionally its results may be affected by ambient noise, such as wind, mines induced seismicity and subsidence. Indirect methods are able to give a rough description of the geometry, which might be enough for design purposes; therefore, the focus of this work is on the development of a tool for fracture’s dimensions estimation using field injection-pressure records (net-pressure analysis). To achieve this objective, we re-examine the problem of hydraulic fracture propagation driven by a Newtonian fluid, propagating in an impermeable, linear-elastic medium for a plain-strain and axisymmetric geometries. For both cases, we develop mode I numerical resolution strategy defined in terms of inlet values, such as pressure and opening at the injection point, avoiding tip behaviour. Based on these propagation conditions, we propose a net-pressure analysis, which considers fracture resistance, viscous dissipation, local mass conservation, and solid-fluid coupling.

This paper is organized as follows. First, the mathematical formulation is described. Then, new parameters based on inlet values and fracture length (radius), are introduced and propagation conditions with these parameters are generated. Finally, results of these criteria in numerical modelling and the development of net-pressure analysis founded on the previously generated propagation criteria are shown.

Mathematical formulation

A symmetric fracture under plain-strain assumption and a radially symmetric fracture are considered, both have half-length (or radius) l(t), and propagate through an impermeable, linear-elastic solid, driven by the injection of a incompressible Newtonian fluid in the center of fractures. The medium is characterized by a plain-strain Young’s modulus E’, fracture toughness K_Ic, and a constant confining stress σ₀, while the injection is characterized by the fluid viscosity p and a constant injection rate Q₀ (see Figure 1).

 

Figure 1: Schematic representation of the problem

 

The problem consists in determining the evolution of the fracture’s half-length l(t), opening w(x,t) and the net-pressure p(x,t) = p_f(x,t) - σ₀(x), where p_f is the total fluid pressure, t is the injection time, and x is the coordinate from the inlet to the tip. The model assumptions are: the fracture is completely filled with fluid; therefore the solution does not depend on σ₀; the flow is unidirectional and laminar; the viscous shear stresses in fracture’s walls are negligible; and the fracture propagates in mobile equilibrium. To impose the quasi-static propagation (mobile equilibrium), a new strategy is introduced.

The opening of the fracture walls w(x,t) produced by the net-pressure p(x,t) can be obtained by means of integral equations, which shows the non-local character of the problem: 

where m = 1 and m = 2, for plane-strain and radially symmetric geometries, respectively. The elastic kernels G₁ and G₂ are given in Sneddon (1951) and Sneddon and Lowengrub (1969).

By combining the local fluid mass balance and Poiseuille law we obtain Reynolds’s equation, which is used to model the fluid flow inside the fracture:

where v_tip(t) is the propagation velocity and the right hand term comes from the assumption that fluid reaches the tip of the fracture, i.e. q(l^-,t) = v_tip(t)w(l^-,t), superscript ^-indicate left limit. Fluid injection is modeled by a source at the center of the fracture.

Also, the global mass conservation is considered:

The uniqueness of the solution is achieved applying a quasi-static propagation condition such as K_I = K_Ic; however, in this paper this condition is enforced in a different manner, which is presented and discussed in the next section.

Proposed method

The complex near-tip behaviour of the hydraulic fractures makes it difficult to impose the classical quasi-static criterion K_I = K_Ic. For this reason, we translate the evaluation of the propagation condition to the injection point, avoiding calculations in the tip. Furthermore, hydraulic fracture’s evolution depends on their energy dissipation processes, so, it seems reasonable to use an energy parameter in numerical modelling. Therefore, we developed new parameters, W_m, which depends on inlet values, such as pressure and opening at the injection point, and with the same units that the energy release rate (Joules/ m²):

For the axisymmetric case we choose w(0,t) instead of p(0,t) because the last one is singular at the inlet. Based on existing semi-analytical solutions (Garagash, 2006; Garagash and Detournay, 2007; Savistski and Detournay, 2002), we described the behaviour of W_m for viscosity W_mV and toughness W_mT dominated regimes, i.e. when the energy is expended mainly in viscous flow or fracture creation, respectively. The transition between these two limiting cases is controlled by the dimensionless toughness T_m:

It is preferable to have an expression W_m in terms of dimensionless toughness T_m, because this parameter Controls the evolution of fracture. Hence, we define the dimensionless form of W_m as = W_m/W_mV:

The existing semi-analytical Solutions (Garagash, 2006; Garagash and Detournay, 2007; Savistski and Detournay, 2002) only apply for limited cases; then, we generated continuous functions for the entire range T_m = [0, ∞), which can be included as propagation condition for the hydraulic fracturing algorithm. The form of these functions is:

where W_^DmIN introduces the interaction and it is used to fit the function in the transition between viscosity and toughness regimes; W_^DmIN must have the following characteristics:

• Strictly positive in T_m = [0, ∞)

• lim_Tm→∞ W_mIN = 1 and lim_Tm→0 W_mIN = 0

In line with these requirements, we suggest a simple rational function of the form:

where for the plane-strain case a₁ = -0.4246, b₁ = 4.9317, and for the axisymmetric case a₂ = -0.0253, b₂ = 4.0940. These coefficients were obtained by imposing the continuity of the expression (11) for both cases. By combining the above equations, we obtain quasi-static criteria for the plane-strain case in expression (12), and for the axisymmetric case in expression (13):

Applications

Algorithm implemented

The solution strategies to impose the quasi-static propagation were implemented in an algorithm which solves the hydraulic fracture propagation’s problem in a coupled manner. For a given fracture half length (or radius) l(t), the algorithm solves iteratively the governing equations to achieve the time which meets quasi-static propagation according to the (12) and (13). The elasticity equation (1) and the fluid flow equations (2) and (3) are solved independently, and are iteratively coupled through a fixed-point method (Figure 2). A boundary element method is used to solve the elasticity equation (Gordeliy and Detournay, 2011; Crouch and Starfield, 1983); for the fluid flow a finite difference method is employed to calculate de pressure gradient ∂p/∂x(x,t), then, the net pressure p(x,t) is obtained by integrating ∂p/∂x(x,t) and by imposing the global mass conservation (3).

The coupling process starts with trial values for the opening w_k and net pressure p_k. A new pressure p_k+1/2 is obtained through the fluid flow equations, and afterwards the pressure is updated as p_k+1/2 = a p_k+1/2 + (1 - a) p_k (where 0 < a < 1/2). The updated pressure is used in the elasticity equation to obtain w_k+_1/2, and the updated opening is calculated as w_k+1 = aw_k+1/2 + (1 - a)w_k. The convergence of the iterative coupling is controlled by a tolerance on the norm of the pressure’s vector between successive steps. Once the convergence is reached, the quasi-static propagation is iteratively accomplished by a combined bisection-Müller method.

 

Figure 2: Algorithm implemented to solve the coupled problem of
hydraulic fracture propagation

 

Numerical modelling

In this section we show the results obtained with the proposed strategies and compare it against those obtained by existing semi-analytical Solutions and existing softwares. To compare the results, we use the dimensionless half-length in viscosity:

In our calculations, we used a uniform grid of boundary elements along the fracture, with the same pressure tolerance of 10^-3 for the coupling algorithm and for the propagation condition. For the plane-strain case, Figure 3 compares our results with those obtained numerically by Garagash (2006) and Garagash and Detournay (2007) and with the semi-analytical solutions of Adachi (2001). It is possible to see that our results agree with the existing solutions.

 

Figure 3: Comparison of computed plane-strain dimensionless half-length γ_1V by
the proposed implementation against semi-analytical solutions (Garagash,
2006; Garagash and Detournay, 2007) and a numerical solution (Adachi, 2001)

 

Figure 4 shows how the dimensionless radius γ_2V varies with the dimensionless toughness T₂ which is a description of the fracture’s evolution with the injection time. Our results were compared with those obtained by the Loramec software (Desroches and Thiercelin, 1993), and with the semi-analytical solution (Savistski and Detournay, 2002). The results agree with the Loramec results, however, for T₂ < 0.5 differences are observed. This discrepancy is due to the assumption of injection as a point source, while Loramec takes into account the borehole radius.

 

Figure 4: Comparison of the computed axisymmetric dimensionless radius γ_2V by
the proposed method against the Loramec results and the semi-analytical solution

 

The proposed strategies have the advantage that despite of the pressure singularities of the problem, they do not require a fine grid in the solution of the elasticity equation. In both cases only 10 boundary elements were used.

Diagnostic techniques

From the proposed criteria (12) and (13), it is possible to isolate the half-length l(t); then they could be used as an estimator of the fracture’s length. For the plane-strain case we can write (12) as:

where v_tip = dl/dt. Integration of (16) gives an expression for the half-length in terms of net-pressure:

In (17), p(0,t) is the inlet net-pressure obtained from the field recording, and Q₀ = Q_T/h_f, being Q_T (m³/s) the total injection rate, while h_f is the out-of-plane dimension of the 2D fracture. tI and tF are the times in which the propagation starts, and the time when injection ends, respectively. To obtain l(t) is necessary to set h_f and calcúlate the right-hand of (17).

The plane-strain geometry (KGD in the literature) is appropriate for fractures in which h_f > 2l, i.e. fractures in stratified rocks, such as sedimentary formations, or in treatments with a long pressurized interval (the zone where the fluid can enters to the rock mass).

In the axisymmetric case, it is possible to isolate l(t) directly from (15):

To the author’s knowledge, it is impossible to measure the inlet’s opening in the course of the injection; however, it is possible to identify generated fractures by performing borehole televiewer loggings before and after the hydraulic fracturing.

To obtain the final radius of the fracture by (18), it is necessary an aperture at the end of the injection w(0,tF); nevertheless, the after-fracturing measurement needs to be done once the fluid pressure has dissipated and fracture has closed (t = tR). Therefore, we are capable to measure the ‘residual opening’ w(0,tR) generated by the water remaining inside the fracture (into the voids between asperity-supported pinched regions).

To overcome this restriction, an elliptical shape (the shape produced by constant pressure) is assumed for both openings w(0,tF) and w(0,tR). Under this assumption, w(0,tF) and w(0,tR) can be written in terms of their respective inlet’s openings, and also the volume of fracture is lineal with respect to the inlet's openings.

On the other hand, in the field, it is possible to measure the volume of fluid that comes back to the pump system once the injection ends, V_rec. Then, by mass conservation we can obtain w(0,tF) as:

Finally, an estimation of the fracture’s radius is obtained replacing (19) in (18). Also, instead of using V_rec, it is possible to use the volume of fluid stored into the fracture, which can be measured through a tiltmeter analysis.

It is important to discuss the validity of the presented diagnostic techniques, mainly because a rock mass contains pre-existing discontinuities, material heterogeneities and variable stresses. For plane-strain and axisymmetric case, the solutions are valid for low variations in material properties and field stresses along the propagation path.

Since the solutions do not consider the interaction between pre-existing discontinuities and hydraulic fractures, they are suited for low-fractured rock masses, or rock masses in which their discontinuities have hard infill and high in-situ stresses. Also, attention must be paid to the spacing between hydraulic fractures; if it is small, interaction among fractures take a place and the solution is no longer valid (a measure of the "smallness" of the spacing is given in Bunger et al., 2012).

Conclusions

The hydraulic fracture’s propagation analysis depends strongly on the parameters, which define the dominant dissipation mechanism; even the singularity type changes depending on whether the fracture generation or the viscous flow is the controlling phenomena. For this reason, the conventional approach of linear elastic fracture mechanics may not be appropriate for this problem.

For the numerical modelling, we proposed solution strategies for the plane-strain and penny-shaped geometries. The strategies consist in the construction of continuous functions of inlet’s pressure and inlet’s opening (for the plane-strain and penny-shaped cases, respectively), which depend on the dimensionless toughness T_m. These functions can be used for any combination of hydraulic fracture’s parameters (0 ≤ T_m < ∞). Also, the developed method is efficient and with a coarse grid we can obtain satisfactory results, and it is not affected by inlet and tip singularities.

Further, we indicate the advantage of writing propagation conditions in terms of inlet parameters and fracture dimensions, because they can be converted in simple tools of fracture’s length assessment. These tools take into account; viscous dissipation, local fluid-mass conservation, fractures resistance and coupling effects.

In addition to the simplicity of the developed tools, they could be useful when the typical tools for fracture’s length estimation, such as tiltmeters or micro-seismicity mapping, are exposed to noisy conditions such as mining-induced seismicity and subsidence, which can affect the results.

 

References

Abou-Sayed, A.S. (1994). Safe injections pressures for disposing of liquid wastes: a case study for deep well injection. Proceedings of the Second SPE/ISRM Rock Mechanics in Petroleum Engineering. Balkema, 769-776

Adachi, J. (2001). Fluid driven fracture in permeable rock. PhD thesis, University of Minnesota

Adachi, J., Siebrits, E., Peirce, A. and Desroches, J. (2007). Computer simulation of hydraulic fractures. International Journal of Rock Mechanics and Mining Sciences 44: 739-757

Akulich, A.V and Zvyagin, A.V. (2008). Numerical simulation of hydraulic fracture crack propagation. Moscow University Mechanics Bulletin 63(1), 6-12

Araneda, O. and Sougarret, A. (2007). Lessons learned in cave mining: El Teniente 1997-2007. Proceedings ofthe First International Symposium on Block and Sub-level Caving. Cape Town

Bunger, A., Zhang, X. and Jeffrey, R. (2012). Parameters affecting the interaction among closely spaced hydraulic fractures. SPE Journal 17(1), 292-306

Crouch, S.L. and Starfield, A.M. (1983). Boundary element methods in solid mechanics. London, George Allen & Unwin

Desroches, J., Detournay, E., Lenoach, B., Papanastasiou, P., Pearson, J.R.A., Thiercelin, M. and Cheng, A. (1994). The crack tip region in hydraulic fracturing. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 447(1929), 39-48

Desroches, J. and Thiercelin, M. (1993). Modelling the propagation and closure of micro-hydraulic fractures. International Journal of Rock Mechanics Mining Sciences and Geomechanics Abstracts 30(7), 1231-1234

Detournay, E. (2004). Propagation regimes of fluid-driven fractures in impermeable rocks. International Journal of Geomechanics 4(1), 35-45

Garagash, D. and Detournay, E. (2007). Plane-strain propagation of a fluid-driven fracture: small toughness solution. Journal of Applied Mechanics 72(6), 916-928

Garagash, D. (2006). Plane-strain propagation of a fluid-driven fracture during injection and shut-in: asymptotics of large toughness. Engineering Fracture Mechanics 73(4), 456-481

Gordeliy, E. and Detournay, E. (2011). Displacement discontinuity method for modeling axisymmetric cracks in an elastic half-space. International Journal of Solids and Structures 48(19): 2614-2629

Haimson, B.C. and Cornet, F.H. (2003). ISRM suggested methods for rock stress estimation - Part 3: hydraulic fracturing (HF) and/or hydraulic testing of pre-existing fractures (HTPF). International Journal of Rock Mechanics & Mining Sciences 40(7-8), 1011-1020

Jeffrey, R.G., Settari, A., Mills, K.W., Zhang, X. and Detournay, E. (2001). Hydraulic fracturing to induce caving: fracture model development and comparison to field data. Proceedings of the 38th US Rock Mechanics Symposium. Washington DC, USA, 251-259

Moss, A. (2011). Block caving. Presentation at the 20^th BMO Capital Markets Global Metals & Mining Conference. Miami, USA

Moss, A.V and Maron, I.A. (1987). Computational Mathematics: worked examples and problems with elements of theory. 4^th ed. Moscow: Mir Publishers

Savistski, A. and Detournay, E. (2002). Propagation of a penny-shaped fluid driven fracture in an impermeable rock: asymptotic solutions. International Journal of Solids and Structures 39(26), 6311-6337

Sneddon, I.N. (1951). Fourier Transforms. New York, McGraw-Hill

Sneddon, I.N. and Lowengrub, M. (1969). Crack problems in the classical theory of elasticity. New York, John Wiley & Sons

Spence, D.A. and Turcotte, D.L. (1985). Magma-driven propagation of cracks. Journal of Geophysical Research: Solid Earth 90(B1), 575-580

Valkó, P. and Economides, M.J. (1995). Hydraulic fracture mechanics. John Wiley & Sons, Chichester, UK

van As, A., Jeffrey, R.G., Chacón, E. and Barrera, V (2004). Preconditioning by hydraulic fracturing for block caving in a moderately stressed naturally fractured orebody. Proceedings of MassMin 2004. Santiago, Chile, 535-541

 

---

Fecha de entrega: 5 de enero 2015
Fecha de aceptación: 22 de septiembre 2015