Hydraulic fracturing is the process of transmitting pressure by liquid or gas to create cracks or to open the existing cracks in hydrocarbon bearing rocks underground. The purpose of hydraulic fracturing in an oil or gas reservoir is to enable the oil or gas to flow more easily from the formation to the wellbore. Mathematical modeling of the hydraulic fracturing process is usually performed to predict different aspects and phenomena through the process. For this reason, many models have been developed over the past several decades. All of these solutions are approximated as they require assumptions about either the fracture opening or the pressure field. Such assumptions are necessary because of the difficulty in treatment of the complex fracture geometry growing under different stress and well conditions. This paper surveys a review of the most commonly used and accurate existing hydraulic fracturing prediction.Keywords: Hydraulic fracturing; modeling; two-dimensional; three-dimensional.
Hydraulic fracturing is a widely used technique for the exploitation of oil, gas, and geothermal resources [1, 2]. The technique of hydraulic fracturing was introduced to the petroleum industry during the 1930s when Dow Chemical Company discovered that by applying a large enough down-hole fluid pressure, it was possible to deform and fracture the rock formation to have a more effective acid stimulation, and is now a standard operating procedure . By 1981, more than 800,000 hydro fracturing treatments had been performed and recorded. Nowadays, about 40% of all currently drilled wells are hydraulically fractured . In some cases, hydraulic fracturing is used to improve the connection between a well and a productive reservoir . In other applications, hydraulic fracturing is used to engineer the reservoir, creating or stimulating fractures in a low permeability matrix. Two well-known examples of the latter case are enhanced geothermal systems (EGS)  and gas shale stimulation .
Hydraulic fracturing occurs when fractures initiate and propagate as a result of pressure applied by a fluid inside the fractures. Fractures in the earth's crust are desired for a variety of reasons, including enhanced oil and gas recovery, re-injection of drilling or other environmentally sensitive wastes, measurement of in situ stresses, geothermal energy recovery, and enhanced well water production. The size of fractures can vary from a few meters to hundreds of meters . A typical hydraulic fracturing process in the petroleum industry is shown in Figure 1. The process is to work as the following steps :
It is generally assumed that the induced fracture has two wings, which extend in opposite directions from the well and is oriented, more or less, in a vertical plane. Other fracture configurations, such as horizontal fractures, are also reported to occur, but they constitute a relatively low percentage of situations documented [4, 8].
A successful hydraulic fracturing may increase the production up to 10s of times, making the technique economically attractive . Experience indicates that at a depth of below 600 meters, fractures are usually oriented vertically. At shallow depths, horizontal fractures have been reported . The fracture pattern, however, may not be the same for different types of soil and rock.
The geometry of the induced fracture is significantly affected by the rock’s mechanical properties, in-situ stresses, the rheological properties of the fracturing fluid and local heterogeneities such as natural fractures and weak bedding planes. For performing hydraulic fracturing in an isotropic and homogeneous medium, the in-situ stress state is the controlling factor on fracture development . Generally, because of reservoirs do not satisfy these ideal conditions, it should understand the effect of other factors on hydraulic fractures .
Modeling of fracture propagation by hydraulic induced fractures is of great interest in order to define the required amount of fluid, injection pressure, and proppant volume and to predict the effectiveness and feasibility of the treatment. Hydraulic fracturing is intrinsically a three dimensional nonlinear coupled problem, where fluid flow and diffusion into rock formation, fracture propagation, and inelastic rock deformation are mechanisms to be described by the model . A hydraulic fracture will grow in the direction normal to the smallest of the three principal stresses as it tends to open in the direction of least resistance. For most reservoir depths of interest in the petroleum industry, the smallest principal stress is in the horizontal plane, which restricts fractures to the vertical plane . Therefore, one may map fractures on the horizontal plane. However, an ideal model should include the three-dimensional aspects of the problem; the available techniques in the literature for three-dimensional analysis of fractures are computationally so expensive that current computers can only handle one or two fractures. Therefore to save the computational efforts and avoid further complexities, the investigation was limited to the two-dimensional analysis. However, two dimensional models will provide a framework for further developments to three dimensional analyses. Despite the development of sophisticated 3D fracture modeling tools, hydraulic fracture behavior is difficult to predict with any degree of confidence.
This article presents a review of the most commonly used and accurate existing hydraulic fracturing prediction models for reservoir rocks with emphasis on their qualitative and quantitative predictive capability, their scientific basis, and their limitations.
An easily static worked procedure for calculating the fluid displacement that accompanies hydraulic fracturing was introduced by Piggott and Elsworth . The solution applied to the particular case of a Perkins and Kern, Nordgren (PKN) [13, 14] hydraulic fracture subject to fracturing fluid loss to the formation at a rate that is equal to the rate of fluid injection. This is referred to as a static solution, since it assumes that the fluid particle is static with respect to the calculation of velocity .
A schematic illustration of a PKN hydraulic fracture is shown in Figure 2. As it can be seen from this figure, two symmetrical fracture segments propagate away from the wellbore with a constant height, H, to a maximum length of L=Lp at the end of fluid injection time of t= tp.
Other parameters that characterize fracture extension and fluid displacement are the diffusivity of the formation (D), porosity of the formation (n), the fracturing fluid leak-off coefficient (CL), and the rate of fracturing fluid injection (Q). The extension of a PKN hydraulic fracture subject to fluid loss at the rate of fluid injection is given by [14, 15]:
where the local rate of fluid loss to the formation through the two opposing fracture walls is defined by:
in this equation, τ is defined as time at which the formation is first exposed to fracturing fluid. Fluid flow through the formation occurs in the plane of x-y and fluid displacement is indexed by the displacement components in the directions of coordinate axes as:
where x0 and y0 and xf and yf are the initial and final positions of a reference fluid particle, respectively. ∆xj is the increment in fracture length. Then, the total displacement of the particle is given by:
Here, T denotes the transmissivity of the formation. In this solution it is assumed that the applied velocity to the particle is unchanged with respect to the motion of the particle due to sufficiently small displacement of the reference fluid particles. This approximation would results in the following equations:
In this equation, geometry is expressed in dimensionless form relative to the length of the fracture at the end of fluid injection xd =x/Lp and yd =y/Lp. In Eq. (5), integration implies the superposition of displacement increments due to fluid loss along the length of the fracture with
The major limitation of the static model is the assumption that the velocity of the fluid particle is everywhere equal to the velocity calculated at the initial location of the particle. From the basic theory of the two-dimensional static model, it is expected that the accuracy of this model will degrade as the displacement of the particle increases, and that a dynamic solution that explicitly represents the motion of the particle is required for the case of large displacement magnitudes that is considered extensively in section 2.2.
2.2. Dynamic Model
In this section, a dynamic solution of a PKN hydraulic fracture subject to high fluid loss is presented . The dynamic model is derived through the superposition of velocity increments that correspond to the discretization of fluid loss to the formation with respect to time and position along the length of the fracture. To start developing of the hydraulic fracturing dynamic model, it should begin with the hydraulic head induced by continuous, point fluid injection as stated by the following equation (with using the same nomenclature to the section of 2.1) :
In this equation, Q denotes the rate of fluid injection, T the transmissivity of the formation, D the diffusivity of the formation, t time relative to the onset of fluid injection and u is defined as:
Where r is the distance between the fluid source and the observation location. From Eqs. (7) and (8), it can be concluded that the advective velocity induced by a finite injection event can be expressed as:
Here, Qi,j, ti, and ∆t, denote the rate, onset, and duration of the Injection event; H and n are the thickness and porosity of the formation; and rj is the distance between the injection and observation locations (subscripts i and j indicate the timing and position of the injection event, respectively). Schematic domain of the dynamic model is illustrated in Figure 3. Partitioning the radially directed velocity increment given by Eq. (9) into components in the directions of x- and y-axes yields:
From Eqs. (1) and (2), it may determine the variation of fluid injection with time and position during the period of fracture extension as:
where time is expressed in dimensionless form relative to the duration of fluid injection, td= t/tp, and ∆xj, is an increment of fracture length as shown in Figure 3. Substitution of Eq. (12) into Eq. (11) yields:
These equations describe the velocity contributed by fracturing fluid loss at a particular time and position along the length of the fracture. The velocity resulting from fluid loss over the duration of injection and length of the fracture is obtained by superposition of these velocity increments. In integral form, this can be written as follows:
At the end of fluid injection, fracturing fluid loss to the formation ceases since the high rate of fluid loss precludes the retention of the fluid within the fracture. Thus, the limits of integration with respect to time are ti,d = 0 and t*i,d = 0, where the t*i,d is the lesser of the current time and the time at the end of fluid injection. The limits of integration with respect to position are xj,d = - x*j,d and xj,d = x*j,d where x*j,d is the lesser of the current length of the fracture and the length of the fracture at the end of fluid injection.
Eq. (14) is an integral equivalent of a system of nonlinear ordinary differential equations and is subject to the initial conditions x(t) = x0 and y(t) = y0, at t = 0. The desired output of the dynamic solution are the values x(t) = xf and y(t) = yf, corresponding to t→∞.
Two-dimensional modeling of fluid displacement through a hydraulic fracturing process using the dynamic model is considerably more complicated and computationally time consuming than the analogous exercise using the static model discussed in the section of 2.1. This model is further complicated by the singular behavior of the integrands in Eq. (14) and by the lack of a definite upper bound for the solution of the system of differential equations. The singular behavior of the integrands may be managed through the application of Gauss-Chebyshev quadrature . An approximate upper bound for the solution of the system of differential equations may be obtained by monitoring the motion of the fluid particle in terms of the characteristic time as:
In this equation, tc is the time to peak velocity assuming instantaneous injection at the wellbore. In this approach, the motion of the particle is determined over increments of the characteristic time, and iteration is stoped when the displacement of the particle over an increment is less than a nominal fraction of the calculated, total displacement of the particle.
The procedure of obtaining the solution of the dynamic model would results in many orders of magnitude more than the static model. This supremacy of static solution indicates its value as an alternative to the dynamic solution. Indeed, the static calculation of fluid displacement for complex fracture and formation geometries is sufficiently demanding computationally that the additional demands imposed by the dynamic solution would be likely to preclude routine analysis.
2.3. Finite Element Cohesive Zone Model
One of other well-known and most applicable 2D models of hydraulic fracturing that takes into account the inelastic deformations that arise in the fracture process zone, is the cohesive zone model (CZM), proposed by Barrenblatt . With this method the deformations at the crack tip prior to fracture propagation are accounted for and energy dissipation occurs in a finite region. Shorter and wider fractures are thus computed with the CZM. The size of this region is a function of the rock material properties (Young’s Modulus and fracture toughness), the fluid’s viscosity and the in-situ stresses . As proposed by Chen et al. , the cohesive zone model not only gives more realistic results but also simplifies the numeric procedure for solving the model (finite element modeling), once pre-determination of the crack tip is not required. Crack initiation and propagation are natural outcomes of the solution of the CZM.
A finite element model based on the cohesive zone model for the analysis of vertical well hydraulic fracturing is presented by Carvalho et al. . In this method, the rock softening behavior is included via a scalar damage constitutive model in terms of a traction-separation law. This law defines the propagation criterion. Due to lack of available data a linear softening material is assumed.
This traction-separation law is presented in Figure 4. The area under the stress-displacement curve equals the strain energy release rate (GC) when the size of the cohesive zone is small compared to the crack length. For elastic material behavior, the fracture toughness, KIC is related to GC through the following relation:
In this equation, E denotes the Young’s modulus and V is the Poisson ratio. For given values of the fracture toughness, KIC and rock tensile strength, Tmax the stress-displacement law is uniquely defined. The opening displacement, σf for which the fracture traction falls to zero is calculated from the following equation:
Here, the fracture energy release rate (GC) is obtained from Eq. (16). The critical opening displacement to activate material damage (σ0) is related to σf through a so called separation coefficient of α:
To complete the model description, the initial cohesive stiffness, K, is given by:
It should be mentioned that both the critical separation (σ0) and the damage stiffness change according to the chosen value of α, GC kept constant. Small values of α imply that a large fraction of the cohesive energy is consumed in the damage phase, in other words, there is large fracture process zone and the fluid front penetrates significantly into it. On the other hand, for large values of α, most of the energy is stored in the cohesive element as elastic energy. In this case the fracture process zone is small and, correspondingly, so is the penetration distance of the fluid front inside de fracture process zone. A thorough discussion on the effect of the separation coefficient on fracture propagation is presented by Chen et al. .
For a cohesive frictional dilatant material, the inelastic rock material behavior follows the Mohr-Coulomb flow theory of plasticity. Moreover, while model development associative behavior with constant dilatation angle is considered. These assumptions are justified by the presence of high confining stresses prior to crack propagation and to a decrease in the initial in-situ mean pressure near the crack tip during propagation .
The model includes the effect of the viscous fluid flow, considering an incompressible, uniform, Newtonian fracturing fluid. The continuity equation and momentum balance will give the lubrication equation relating the pressure gradient to the fracture width. Constant flow rate at the wellbore and zero fluid pressure at the fluid front are the boundary conditions to be enforced in the model. One disadvantage of this model is that in this model, the effect of leak-off from the fracture surface into the rock formation is not considered and it may results in some errors. An application of the finite element cohesive zone model to a specific case is provided by Carvalho et al. . From their results, the model well predicted the investigated system behaviors.
3. Three-dimensional (3D) Hydraulic Fracturing Models
3.1. General 3D Hydraulic Fracturing Model
A simple and general model by considering the effects of fluid flow inside the fracture, fluid leak-off through fracture walls and elastic response of the surrounding porous media was developed and presented by Devloo et al. .
A fracture is a region in the tridimensional space varying with time (t). Based on the model presented by Devloo et al. , consider a planar configurations where the variations of fracture longitude is in the x direction, the fracture height is in the vertical axis z, and fracture opening is in the y direction, as indicated in Figure 5. Typically, hydraulically induced fractures are large but quite narrow. Moreover, assume that the fluid flow in the fracture follows the pattern of flow between parallel plates and that the pressure is uniform along the y axis. Thus, the main variables inside the fracture are pressure p and opening displacement w, which only depend on (x, z, t). Here, a fracture occupies a region defined by points (x, y, z) such that:
By assuming normal traction and using single and double layer elastic potentials, an equation can be derived for the pressure p(x, z, t) acting on the fracture surface in function of the surrounding stress distribution σ(x, y) and the fracture opening w(x, z, t). It is given by the singular integral, considered in the Cauchy principal value sense,
Where , G the shear modulus and ν is the Poisson’s ratio of the rock formation. The problem of a plane fracture of arbitrary shape in an infinite elastic medium has been considered in .
In this model, a bi-dimensional fluid flux q = q(x, z, t) was considered inside the fracture. Typically, the flux depends on the pressure gradient and the fracture opening. For Newtonian fluids, with constant viscosity of μ the relation can be written as:
On the other hand, there is also interest in non-Newtonian fluids satisfying a rheological power law behavior. In this cases, an equation of flux can be rewritten as:
here k is the consistency index and α is the power law exponent [19, 21]. In an ideal case in which we have a hydraulic fracturing process with 100% of efficiency, all injected fluid would contribute to the fracture opening. However, for fractures in a porous media, part of the injected fluid is filtered through the fracture walls. A usual model for this filtering effect is Carter’s model :
This equation express the leak-off rate (Q) per unit fracture area at a current time t, where τ = τ(x, z) denotes the instant at which the fracture opens at position (x, z). The parameter vsp is called spurt-loss since it represents the amount of fluid instantaneously filtered to the porous media, and C is a filtration parameter.
The conservation law is applied for modeling of fluid flow inside the fracture through the general hydraulic fracturing model:
This equation is obtained by expressing the variation of fluid content in a control cell as the result of the balance of fluid flow in and out the cell and the fluid lost to the porous media during a period of time. For all test functions v, solution of Eq. (25) in its variational form is supposed to verify:
Here ∂Ω denotes the boundary of the computational domain. A finite element approach was applied in order to evaluate the above mentioned model by Devloo et al. . For this reason, a rectangular computational region on the plane x × z was considered as Ω, where the opening or closing of nodes determines the propagation of the fracture, and results were obtained. From the results, the presented model provides very well estimation of hydraulic fracturing and it is nominated as a three dimensional powerful model. 3.2. Improved Flow-Stress-Damage (FSD) Model
A three-dimensional model (RFPA3D-Parallel or Improved Flow-Stress-Damage Model) that considers the coupled effects of seepage, damage, and the stress field is introduced by Li et al. . RFPA3D-Parallel is an extension of two-dimensional Rock Failure Process Analysis (RFPA2D) [23, 24]. RFPA2D is a 2D finite element code that can simulate the fracture and failure process of quasi-brittle materials such as rock. The constitutive law of this model considers strength and stiffness degradation, stress-dependent permeability for the pre-peak stage, and deformation dependent permeability for the post-peak stage. Three coupled processes of (1) mechanical deformation of the solid medium induced by the fluid pressure acting on the fracture surfaces and the rock skeleton, (2) fluid flow within the fracture, and (3) propagation of the fracture are considered through this model. Different assumptions were considered in developing of the FSD model including: