Introduction

The Hydraulic Fracturing (HF) is a process through which a fluid is injected into a well under high pressure and, as a result, long cracks are created in the reservoir rock at the injection site. These cracks create paths in the rock with high permeability that enhances the amount of oil production from the well (Fig. 1).

Fig. 1
figure 1

Hydraulic fracturing (HF) process (Veatch Jr. and Moschovidis 1986)

Full size image

Although the main application of this operation is in the oil and gas industry, it has various applications in other areas of engineering. In civil engineering, this method is used to measure in-situ stresses of rocks in dam and tunnel construction projects. In environmental engineering, HF is used to increase soil permeability and enhance the effectiveness of decontamination methods for polluted soils. Geothermal energy extraction also benefits from HF in hot and dry rock masses. From 1960 to the end of the 2012s, approximately one million hydraulic fracturing operations in the United States and approximately 2.5 million Hydraulic Fracturing worldwide have been recorded. Since 1993, nearly 40% of oil wells and 70% of gas wells in the United States have been supplemented by Hydraulic Fracking (Veatch and Moschovidis 1986; King 2012; Pak 1997).

The first simplified analytical models for analyzing and designing HF initiated in 1950s. Among the earliest published papers in this area were Zheltov (1955), and Perkins and Kern (1961). Later, Geertsma and De Klerk (1969) and Nordgren (1972) continued these early works. These researchers assumed that the induced cracks propagate in one direction only (plane strain conditions) with no fluid flow perpendicular to the crack wall. They obtained two different analytical solutions for maximum crack opening (w), crack length (L) and injection pressure (p) which are called KGD and PKN models. These models will be briefly introduced later in this paper. The economics of Hydraulic Fracturing is entirely dependent on accurate prediction of crack parameters. For example, under-prediction of the crack length renders using higher pressure and higher flow rate for injection. On the other hand, predicting longer length in the design, causes low effectiveness of real hydraulic fracture. If the pressure considered in the HF design is below the real required value, the hydraulic fracturing operation cannot induce cracks with the desired dimensions. On the other hand, if the design pressure is greater than the real required pressure, the crack dimensions will be more than needed (Ershaghi and Milligan 2004). Therefore, accurate estimation of HF characteristics can determine the real performance of the well after the operation (Fjaer et al. 2008).

There exists a large body of the literature about HF since its invention. Some of the recent works are reviewed here. Fisher et al. (2004) studied the effect of fracture length on reservoir productivity. Lolon et al. (2009) studied the effect of placement and fracture spacing in HF process on the productivity of the reservoir. Chen et al. (2009) simulated HF process using FEM in ABAQUS and numerically studied the effect of cohesive zone length, tensile strength of the reservoir rock, and viscosity of the injected fluid. Zhang et al. (2010) studied the effect of in-situ stresses, rock elastic modulus, rock tensile strength, and the viscosity of the injected fluid on the HF process, by three-dimensional modeling in the ABAQUS software. Sadrnezhad (1999) studied the influence of initial effective principal stresses on HF pressures in a typical porous medium. They found that the ratio of initial effective stresses has a considerable impact on the hydraulic fracturing pressure. Nasehi and Mortazavi (2013) studied the effect of in-situ stress ratios on HF using numerical modelling by employing two-dimensional UDEC4.0 software. Chiti and Pak (2017) investigated the effect of fracture fluid viscosity and compressibility, initial pore pressure, hydraulic conductivity, and the tensile strength of the rock on HF assuming no leak-off using XFEM method in ABAQUS. Han et al. (2018) investigated the effect of heterogeneity on the size of the constituent grains and the in-situ stress ratio on the HF process using PFC2D software. Salighedoost and Pak (2019) investigated the effect of the initial crack length, fracture energy, in-situ stress and the tensile strength of the reservoir rock, leak-off effect and fluid pressure inside the crack on the HF process using ABAQUS. Mohammadnejad and Khoei (2013) studied the effect of permeability of the reservoir rock, flow rate, and fluid viscosity on HF using numerical modelling. Haddad and Sepehrnoori (2014) studied the effect of the injection flow rate, fluid viscosity, fluid leak-off coefficient, elasticity modulus and Poisson’s ratio on HF using the FEM and CZM technique. Later Haddad and Sepehrnoori (2015) studied the effect of the maximum horizontal stress and stress shadow using the XFEM, cohesive zone and phantom node technique. Wang (2015) studied the non-planar HF in the permeable media using XFEM and considering Mohr-Columb plastic model along with the cohesive zone. Wang et al. (2016) used XFEM and CZM to study the effect of the existing cracks on the creation of new cracks (stress shadow effect). Vahab et al. (2018) studied the propagation of HF in a naturally layered domain using XFEM. Saberhosseini et al (2017) studied the effect of the fracture toughness, rock tensile strength, modulus of elasticity, fluid viscosity, Poisson's ratio, and leak-off coefficient on the characteristics of HF using XFEM. Mehrgini et al. (2017) studied the effect of the elasticity modulus, Poisson's ratio, cohesion, friction angle and tensile strength of the HF using FEM. Golovin and Baykin (2018) presented a new numerical model for simulation of propagation of planar cracks and investigated the effect of the pore pressure on HF process.

Although many investigators have studied HF and effects of different factors on its characteristics, none of them has focused on the well-known classical methods such as KGD and PKN.

In this study the effect of the in-situ stress field in the reservoir on the fracture width (w), length (L), and fluid pressure at the crack mouth (P) have been studied using XFEM. The outcome of this research is a number of correction (modification) factors that can be used for improvement of the results of KGD and PKN models for a more accurate design of Hydraulic Fracture process in oil reservoirs.

Formulation

Governing equations in hydraulic fracture modeling

In isothermal hydraulic fracture modeling, the following equations are considered (Zielonka 2014):

  1. 1.

    Equilibrium equation in porous medium (conservation of linear momentum)

  2. 2.

    Porous medium behavioral equation (elasticity in porous medium using Biot theory)

  3. 3.

    Fluid continuity equation in porous medium

  4. 4.

    Fluid continuity equation inside the crack

  5. 5.

    Momentum equation for the pore fluid (Darcy’s Law)

  6. 6.

    Momentum equation for the fracturing (Lubrication Equation)

These equations will be briefly discussed in the following sections.

Porous medium deformation

The deformation of the porous medium of the reservoir can be modeled as the quasi-static deformation of isotropic poroelastic materials. So, the equilibrium equation without considering the body forces is as follows:

$$\sigma_{ij,j} = 0$$
(1)

Considering small strains, the elastic behavioral relationship in a porous medium is as follows:

$$\sigma_{ij} - \sigma_{ij}^{0} = 2G\varepsilon_{ij} + \left( {K - \frac{2}{3}G} \right)\varepsilon_{kk} - \alpha \left( {p - p_{0} } \right)\delta_{ij}$$
(2)

In this equation:

$$\begin{aligned} G & = \frac{E}{{2\left( {1 + \upsilon } \right)}} \\ K & = \frac{E}{{3\left( {1 - 2\upsilon } \right)}} \\ \end{aligned}$$
(3)

in which \(\sigma_{ij}\) and \(\sigma_{ij}^{0}\) are current total stress and initial (in-situ) total stress tensors, respectively. α is Biot’s coefficient, G and K are the elastic shear and bulk moduli, E is Young’s modulus, and \(\upsilon\) is Poisson’s ratio. Abaqus is formulated in terms of Terzaghi effective stresses σ′, defined for fully saturated media as (Zielonka 2014)

$${\sigma }_{ij}^{^{\prime}}={\sigma }_{ij}+p{\delta }_{ij}$$
(4)

In terms of the latter, the constitutive relation takes the form

$$\sigma_{ij}^{^{\prime}} - \sigma_{ij}^{^{\prime}0} = 2G\varepsilon_{ij} + \left( {K - \frac{2}{3}G} \right)\varepsilon_{kk} \delta_{ij} - \left( {\alpha - 1} \right)\left( {p - p_{0} } \right)\delta_{ij}$$
(5)

Defining effective strains as:

$${\varepsilon }_{ij}^{^{\prime}}={\varepsilon }_{ij}-\frac{\alpha -1}{3K}\left(p-{p}_{0}\right){\delta }_{ij}$$
(6)

the constitutive relation simplifies to (Zielonka 2014)

$$\sigma_{ij}^{^{\prime}} - \sigma_{ij}^{^{\prime}0} = 2G\varepsilon_{ij}^{^{\prime}} + \left( {K - \frac{2}{3}G} \right)\varepsilon_{kk}^{^{\prime}} \delta_{ij}$$
(7)

Pore fluid flow

The continuity equation for pore fluid, assuming small volumetric strains, is given by

$$\frac{1}{M}\dot{p}+\alpha {\dot{\varepsilon }}_{kk}+{v}_{k,k}=0$$
(8)

where \({v}_{k}\) is the pore fluid seepage velocity and \(M\) and \(\alpha\) are Biot’s modulus and Biot’s coefficient. \(\dot{p}\) and \(\dot{\varepsilon }_{kk}\) are rate of changes of the pressure and volumetric strain, respectively. These two poroelastic constants are defined by the following identities:

$$\begin{aligned} \frac{1}{M} & = \frac{{\phi_{0} }}{{K_{{\text{f}}} }} + \frac{{\alpha - \phi_{0} }}{{K_{{\text{s}}} }} \\ \frac{1}{{K_{{\text{s}}} }} & = \frac{1 - \alpha }{K} \\ \end{aligned}$$
(9)

where \({K}_{\mathrm{f}}\) is the pore fluid bulk modulus, \({K}_{\mathrm{s}}\) is the porous medium solid grain modulus, and \({\phi }_{0}\) is the initial porosity. The pore fluid is assumed to flow through an interconnected network of pores according to Darcy’s law:

$${v}_{i}=-\frac{ k }{\mu }{p}_{,i}= -\frac{ \overline{k} }{\gamma }{p }_{,i}$$
(10)

in which \(k\) is the permeability, \(\mu\) is the pore fluid viscosity, \(\overline{k }\) is the hydraulic conductivity and \(\gamma\) is the pore fluid specific weight. Combining (10) with the continuity equation, the pore fluid diffusion equation is obtained (Zielonka 2014)

$$\frac{1}{M}\dot{p}+\alpha \dot{{\varepsilon }_{kk}}=\frac{ \overline{k} }{\gamma }{p }_{,kk}$$
(11)

Fracture fluid flow

The longitudinal fluid flow within the fracture is governed by Reynold’s lubrication theory defined by the continuity equation

$$\dot{g}+\frac{\partial {q}_{f}}{\partial s}+{v}_{\mathrm{T}}+{v}_{\mathrm{B}}=0$$
(12)

and the momentum equation for incompressible flow of Newtonian fluids through narrow parallel plates:

$${q}_{f}=-\frac{{ g}^{3} }{ 12{\mu }_{\mathrm{f}} }\frac{\partial {p}_{\mathrm{f}}}{\partial s}$$
(13)

where \(g\) is the fracture gap (Fig. 2), \({q}_{\mathrm{f}}={v}_{\mathrm{f}}\) is the fracturing fluid flux (per unit width) across the fracture, \({v}_{\mathrm{T}}\) and \({v}_{\mathrm{B}}\) are the normal flow velocities of fracturing fluid leaking through the top (\({v}_{\mathrm{T}})\) and bottom (\({v}_{\mathrm{B}})\) faces of the fracture into the porous medium, \({\mu }_{\mathrm{f}}\) is the fracturing fluid viscosity, and \({p}_{\mathrm{f}}\) is the fracturing fluid pressure along the fracture surface parameterized with the curvilinear coordinate, \(S\) (Zielonka 2014).

Fig. 2
figure 2

Fracture aperture, width and fracturing fluid flow (Zielonka 2014)

Full size image

Abaqus computes the normal fracturing fluid velocities as:

$$\begin{aligned} v_{{\text{T}}} & = c_{{\text{T}}} \left( {p_{{\text{f}}} - p_{{\text{T}}} } \right) \\ v_{{\text{B}}} & = c_{{\text{B}}} \left( {p_{{\text{f}}} - p_{{\text{B}}} } \right) \\ \end{aligned}$$
(14)

where \(c_{{\text{T}}}\) and \(c_{{\text{B}}}\) are the normal fluid leak-off coefficients in the top and bottom faces of the crack.

KGD and PKN analytical models

KGD and PKN are the most widely used analytical models to estimate characteristics of the hydraulic fracture. In Table 1 the proposed formulas by KGD and PKN models for fracture width (w), fracture length (L), and fluid pressure at the crack mouth (p) are displayed.

Table 1 analytical formulas of KGD and PKN models
Full size table

As can be seen, there are some geomechanical and hydromechanical factors that have not been taken into account in the simplified solutions of PKN and KGD that render these closed form solutions to be approximate.

Numerical simulation of KGD and PKN models

Before evaluating the effect of reservoir in-situ stress on PKN and KGD formulae, numerical simulations have been carried out to validate the results of the numerical method employed for this research.

Rock specifications

In the following numerical simulations, the reservoir rocks are assumed to be homogeneous and isotropic with elastic behavior during hydraulic loading until cracking occurs. The reservoir is considered to be completely saturated (Sr = 1). The effects of possible temperature change during the HF process have been ignored. The rock reservoir specifications for KGD and PKN models are given in Table 2.

Table 2 The rock specifications for KGD and PKN model
Full size table

The parameters in Table 2 are: \(E\) (modulus of elasticity), \(\upupsilon\) (Poisson’s ratio), \(e\) (porosity ratio), \(k\) (permeability coefficient of the reservoir rock), \({G}_{c}\) (reservoir rock fracture energy), \({\sigma }_{\mathrm{t}}\) (tensile strength of the reservoir rock), \(q\) (injection flow rate), \(\mu\) (injection fluid viscosity coefficient), \({c}_{\mathrm{T}}, {c}_{\mathrm{B}}\) (leak-off coefficient into the rock reservoir).

Injection fluid specifications

Fluid flow inside the crack can be divided into two components: (1) Tangential Flow (2) Normal Flow. In this study, only the tangential flow inside the crack is analyzed. The normal flow inside the crack is defined by introducing leak-off coefficients (\({c}_{\mathrm{T}}, {c}_{\mathrm{B}}\)).

The leak-off coefficient can be defined independently with no relation to the permeability coefficient of the reservoir rock. In this study, the amount of leak-off coefficient for the upper and lower crack surfaces is equal to 5 × 10−30 m3/N s. By introducing the pore fluid bulk modulus and the solid grain bulk modulus to the software, the software uses the Biot’s relationship to obtain the effective stress. In this study, the amount of the bulk modulus of the reservoir rock and pore fluid are 40 MPa and 2.15 Mpa, respectively.

Verification for KGD model

In this section a typical problem is analyzed using KGD model, and the numerical results are compared with the closed form solutions. After obtaining good results, the effects of the in-situ stresses on the outcome will be examined. The KGD model consists of a half-circle of 160 m in diameter, a rectangle of 10 m width and 160 m length (Fig. 3). Since the assumption of plane strain in the horizontal direction is applicable in KGD, the model is 2D, and its thickness is unity. Due to symmetry, only half of the reservoir is modeled. In the middle area, a finer mesh is used because the crack propagates in this area (zone 1) (Fig. 3). The model geometry is divided into five zones whose properties are given in Table 3. The type of elements used in this model, is CPE4P (4-node bilinear displacement and pore pressure) which is suitable for simulating Hydraulic Fracture. Size of the domain is considered such that the changes in the pore pressure at the outer boundary (boundary 2) is less than 0.5%.

Fig. 3
figure 3

KGD numerical model

Full size image
Table 3 Geometrical Properties of the KGD model
Full size table

Boundary conditions: In this model, there are two main boundaries, 1 and 2. Due to axial symmetry along the X-axis, boundary 1 is fixed for displacement along the X-axis and fixed for rotation around Y axis. Boundary 2 is far from the crack area, so the mechanical boundary conditions (displacement in the Y and X directions) are fixed and the hydraulic conditions (i.e., pore water pressure) are considered zero.

Initial conditions: In the KGD model, the initial conditions are: (1) The medium is completely saturated (Sr = 1), (2) The initial pore pressure is considered zero. Along boundary 2, after injection (due to the distance of this boundary from the injection zone), the pore pressure will remain zero, the injection rate is 0.1 L/s/m and the injection duration is 10 s (the injection zone is shown in Fig. 3).

In this section, the results of the KGD numerical model, which is in its original form (without in-situ stresses), are compared with the results of the closed form relations. The results of the KGD model are illustrated in Table 4 and Fig. 4.

Table 4 Comparison between the numerical and analytical results of KGD model
Full size table
Fig. 4
figure 4

Diagrams of the width (a), length (b) and mouth pressure (c) of the crack for the original KGD model

Full size image

As shown in Table 4 and Fig. 4, The numerical results of crack width, length, and mouth pressure are very close to the analytical results. The maximum error is related to fluid pressure at the crack mouth which is less than 15%.

Verification for PKN model

The PKN model consists of a half-circle of 1600 m in diameter and a rectangle of 300 m width and 1600 m length. Due to the basic configuration of PKN model which is 3D, the thickness of the mesh is considered 25 m (Fig. 5).

Fig. 5
figure 5

PKN numerical model

Full size image

Due to symmetry, only half of the reservoir is modeled. In the middle area, a finer mesh is used because the crack propagates in this area (zone 1). The model geometry is divided into five zones (Fig. 5) and their properties are given in Table 5.

Table 5 Geometrical properties of PKN model
Full size table

The type of elements used for this model in Abaqus, is C3D8RP (8-node trilinear displacement and pore pressure, reduced integration) which has the efficiency for modeling HF. Mesh size are considered such that the changes in the pore pressure in the remote boundary (boundary 2) of the model are less than 0.5%.

Boundary conditions: In PKN model, the boundary conditions are applied to the surface nodes of the boundaries 1 and 2. In this model, the boundary conditions are: (1) Due to symmetry along the X axis for boundary 1, the displacements along the X-axis and the rotation around Y and Z axes are fixed, (2) Boundary 2 is far from the crack area so the mechanical boundary conditions (displacement in the X and Y directions) and the hydraulic conditions (pore water pressure) at this boundary are zero, (3) To prevent crack mouth nodes from moving in the X and Z directions, their phantom nodes must be fixed in both directions (Fig. 6a), (4) Given that in the PKN model the upper and lower edges of the crack have a width of zero and the crack geometry is elliptical in the vertical direction, in zone 1 at top and bottom edges of the crack, the phantom nodes are fixed in the direction of crack growth along the y-axis (Fig. 6b), and (5) All nodes in the model must be fixed in the Z direction.

Fig. 6
figure 6

a Prevent crack mouth nodes from moving in the X and Z directions, b upper and lower crack nodes have a width of zero, c fluid injection area in PKN model

Full size image

Initial conditions: In PKN model, the initial conditions are: (1) The medium is completely saturated (Sr = 1), (2) The initial pore pressure is considered zero, everywhere.

Along boundary 2, after injection, the pore pressure remains the same as the initial pore pressure. For PKN model the injection rate is 3 (L/s)/m and the injection time is 140 s. In this model,

for creating the same injection flow rate along the crack height, the fluid needs to be injected between the nodes in such a way that the injection flow rate at the side nodes is half of that in the internal nodes (Fig. 6c).

$$\begin{aligned} Q_{{\text{T}}} & = 3 \times 25 = 75\, \frac{{\text{L}}}{{\text{S}}} \\ Q_{1} & = \frac{75}{{\left( {2 \times 2.5} \right) + \left( {4 \times 5} \right)}} \times 1 = 3 \,\frac{{\text{L}}}{{\text{S}}} \\ Q_{2} & = \frac{75}{{\left( {2 \times 2.5} \right) + \left( {4 \times 5} \right)}} \times 2 = 6 \,\frac{{\text{L}}}{{\text{S}}} \\ \end{aligned}$$

The results of the PKN numerical model, which is in its original form (without in-situ stresses), are obtained and compared with the results of the closed form relations. The results of PKN model are illustrated in Table 6 and Fig. 7.

Table 6 Comparison of the numerical and analytical results of PKN model
Full size table
Fig. 7
figure 7

Diagrams of width (a), length (b) and mouth pressure (c) of the crack for original PKN model

Full size image

As shown in Table 6 and Fig. 7, The numerical results of crack width, crack length and mouth pressure are close to the analytical results.

Evaluation of the effects of in-situ stresses

In-situ stresses are one of the most effective parameters on HF operations, Valko and Economides (1996) state that it is very complex to deal with the exact effect of the stress state on HF. They have used the Principle of Superposition to solve this problem and considered the pressure inside the crack as the net pressure. Application of the superposition law requires pure elastic behavior. Despite existence of the plastic zone at the crack tip, due to the small processing area at the crack tip compared to the large size of the hydraulic cracks, the principle of superposition can be used. This simplification has caused the in-situ stress enters only into the equation of fluid pressure at the crack opening in the analytical relations of KGD and PKN models. In this study, the effect of this important parameter on Hydraulic Fracture characteristics is studied in detail. To obtain the in-situ stress, the specific weight of the rock of the layers above the crack level is considered to be 21 KN/m3. Assuming the absence of faults in the area, the vertical stress is calculated as follows:

$${\sigma }_{v}={\int }_{0}^{z}\gamma \left(z\right).dz$$
(15)

In which \(z\) is the depth of crack location.

Using the experimental results of Haimson and Fairhurst (1969), horizontal in-situ stresses can be estimated based on the following relations:

$$\begin{aligned} \sigma_{{{\text{h}},\max }} & = \frac{2}{3}\sigma_{{\text{v}}} \\ \sigma_{{{\text{h}},\min }} & = \frac{1}{3}\sigma_{{\text{v}}} \\ \end{aligned}$$
(16)

In ABAQUS software, to apply the in-situ stresses in the “Load” module, one should introduce “Predefined field” and "Stress” type. In “Edit Predefined field” window, the stress values should be entered as follows:

$$\begin{array}{*{20}l} { \sigma_{{\text{v}}} = \sigma_{Z} = \sigma_{33} = sigma\;33} \hfill & {\quad \sigma_{{{\text{h}},\max }} = \sigma_{X} = \sigma_{11} = sigma\; 11} \hfill \\ {\sigma_{{{\text{h}},\min }} = \sigma_{Y} = \sigma_{22} = sigma \;22} \hfill & {\quad \tau_{XY} = \sigma_{12} = sigma\; 12 = 0} \hfill \\ \end{array}$$

The overburden pressure at depths of 0.25, 0.5, 1, 1.5, 2 and 3 km are calculated and then applied to KGD and PKN numerical models to determine the effect of in-situ stresses on the results. The results obtained from the KGD numerical simulation and closed form solutions are shown in Table 7 and Fig. 8.

Table 7 Results of the numerical simulations considering the in-situ stresses effect in KGD model
Full size table
Fig. 8
figure 8

Diagrams of the crack width (a), crack length (b) and fluid pressure (c) at crack mouth at different depths (different in-situ stresses) in KGD model

Full size image

The results obtained from the PKN numerical simulation and closed form solutions are shown in Table 8 and Fig. 9.

Table 8 Results of the numerical simulations considering the in-situ stresses effect in PKN model
Full size table
Fig. 9
figure 9

Diagrams of the crack width (a), crack length (b) and fluid pressure (c) at crack mouth at different depths (different in-situ stresses) in PKN model

Full size image

Comparing Fig. 8 with 9, reveals the fact that simulation of PKN model is generally more difficult than simulation of KGD model and the number of discrepancies is higher. This is because of two reasons: First is that the assumptions regarding shape of fracture in PKN model is different from reality, especially at the first steps of fracture propagation. The second is that since the crack length to width (\(\frac{L}{W}\)) ratio for PKN is greater than that of KGD, the effect of fluid viscosity (and loss of energy) in PKN model is more than KGD.

In the first study step, in both simulations, it is observed that with increasing the reservoir depth (increasing the in-situ stresses), the width of the crack opening increases and it is almost greater than the value obtained from the closed form solutions. Since the leak-off coefficient is small, the volume of fluid injected into the crack is almost constant, and as the width of the crack increases, its length must decrease. The amount of crack length obtained from numerical results is generally less than the analytical values, as expected. Also, with increasing the reservoir depth (increasing in-situ stresses), to overcome the in-situ stresses to make the cracking possible, the pressure at the crack mouth must increase, which is also shown by the numerical results. The crack mouth pressure obtained from the numerical results is generally greater than the analytical values.

As can be seen, the effects of in-situ stress state on the crack characteristics are considerable. For example, if the depth of the reservoir increases from 1.0 to 3.0 km, the crack width increases up to 2.3 times, the crack length decreases up to 3 times and the crack mouth pressure increases up to 3.25 times.

Depending on the type of reservoir rock (different Poison’s ratios), the lateral pressure coefficients can get values other than \(\frac{2}{3}\) and\(\frac{1}{3}\). Therefore, the maximum and minimum horizontal in-situ stresses (\({\sigma }_{\mathrm{h},\mathrm{max}} , {\sigma }_{\mathrm{h},\mathrm{min}})\) can be different from what have already been assumed. In the following, the effect of the ratios of the maximum and minimum horizontal in-situ stresses on the hydraulic crack characteristics are studied. In the second study step, two different cases for the maximum and minimum horizontal in-situ stresses are considered:

  1. 1.

    The minimum horizontal in-situ stress is kept constant and the maximum horizontal in-situ stress increases with the ratio of 1.1, 1.2, 1.3, 1.5, 2 and 3.

  2. 2.

    The maximum horizontal in-situ stress is constant and the minimum horizontal in-situ stress decreases with the ratio of 1.1, 1.2, 1.3, 1.5 and 3. The constant depth for KGD simulation is 500 m and for PKN is 1000 m.

Case 1: the minimum horizontal in-situ stress is constant and the maximum horizontal in-situ stress increases: For this case, the results obtained from the KGD numerical simulations for 20 s injection are shown in Table 9 and Fig. 10.

Table 9 Results of the KGD model for case 1 (the minimum horizontal in-situ stress is constant and the maximum horizontal in-situ stress increases)
Full size table
Fig. 10
figure 10

Diagrams of variation of crack width (a),crack length (b) and fluid pressure (c) for KGD model in case 1 (the minimum horizontal in-situ stress is constant and the maximum horizontal in-situ stress increases)

Full size image

For case 1, the results obtained from the PKN numerical simulations for 140 s injection are shown in Table 10 and Fig. 11.

Table 10 Results of the PKN model for case 1 (the minimum horizontal in-situ stress is constant and the maximum horizontal in-situ stress increases)
Full size table
Fig. 11
figure 11

Diagrams of the variation of crack width (a),crack length (b) and fluid pressure (c) for PKN model in case (the minimum horizontal in-situ stress is constant and the maximum horizontal in-situ stress increases)

Full size image

Case 2: the maximum horizontal in-situ stress is constant and the minimum horizontal in-situ stress decreases:

In this case, the results obtained from the KGD numerical simulations for 20 s injection are shown in Table 11 and Fig. 12.

Table 11 Results of the KGD model for case 2 (the maximum horizontal in-situ stress is constant and the minimum horizontal in-situ stress decreases)
Full size table
Fig. 12
figure 12

Diagrams of the crack width (a), crack length (b) and fluid pressure (c) for KGD model in case 2 (the maximum horizontal in-situ stress is constant and the minimum horizontal in-situ stress decreases)

Full size image

For case 2, the results obtained from the PKN numerical simulations for 140 s injection are shown in Table 12 and Fig. 13.

Table 12 Results of the PKN model for case 2 (the maximum horizontal in-situ stress is constant and the minimum horizontal in-situ stress decreases)
Full size table
Fig. 13
figure 13

Diagrams of the crack width (a), crack length (b) and fluid pressure (c) for PKN model in case 2 (the maximum horizontal in-situ stress is constant and the minimum horizontal in-situ stress decreases)

Full size image

According to the results obtained in this section, for both models, when the minimum horizontal in-situ stress is kept constant and the maximum horizontal stress is changed, the values of the crack width, crack length, and mouth pressure are almost constant. But when the maximum horizontal in-situ stress is kept constant and the minimum horizontal in-situ stress changes, the values of the crack width, crack length, and mouth pressure change significantly. So, the minimum horizontal in-situ stress (\({\sigma }_{\mathrm{hmin}} )\) is the most effective parameter in studying of the in-situ stresses effect on HF. This was expected because the required fluid pressure for creating crack in rock is dependent to the in-situ-minimum principal stress. This study shows that with increasing the stress ratio (\(\frac{{\sigma }_{\mathrm{hmax}}}{{\sigma }_{\mathrm{hmin}}})\) from 1.0 to 3.0, the crack width decreases up to 1.5 times for KGD model and 2 times for PKN model, the crack length increases up to 1.4 times for KGD model and 3 times for PKN model, and the fluid pressure decreases up to 2.8 times for both models.

Modification of KGD and PKN analytical formulas

Considering that the classical analytical relationships in KGD and PKN models presented in Table 1, predict crack characteristics without considering the in-situ stresses, in this study, these relationships are modified based on the results obtained from numerical simulations, to enhance the efficiency of these formulas. These modifications are applied by introducing coefficients for analytical formulas. The coefficients \({C}_{1}\) to \({C}_{6}\) applied to the analytical formulas of the KGD and PKN models are given in Table 13.

Table 13 The constant coefficients applied to the analytical formulas of the KGD and PKN models
Full size table

In this study, 2 numerical simulations for initial validation of the numerical models, 13 numerical simulations to study of the effect of the in-situ stresses at different depths and 21 numerical simulations to study of the effect of the ratio of the in-situ stresses have been carried out. According to the results obtained from Sect. 5, it was shown that the maximum horizontal in-situ stress does not affect the values of the crack width, crack length and pressure of the crack mouth and only minimum horizontal in-situ stress is effective. Therefore, to modify the closed form solutions, normalized minimum horizontal in-situ stress (\(\frac{{\sigma }_{\mathrm{hmin}}}{{P}_{\mathrm{a}}} )\) is used, where \({P}_{\mathrm{a}}\) is the atmospheric pressure.

In this section, to obtain the modified coefficients of the crack width, crack length and crack mouth pressure, the following steps are used:

Step 1 The value of the crack parameter is obtained by numerical modeling for different \(\frac{{\sigma }_{\mathrm{hmin}}}{{P}_{\mathrm{a}}}\) (using section “Investigation on the effects of in-situ stresses”, numerical results).

Step 2 The value of the crack parameter is obtained by analytical formulas for different \(\frac{{\sigma }_{\mathrm{hmin}}}{{P}_{\mathrm{a}}}\). (using section “Investigation on the effects of in-situ stresses”, analytical results).

Step 3 For different \(\frac{{\sigma }_{\mathrm{hmin}}}{{P}_{\mathrm{a}}}\), the values obtained from step 1 are divided by the values obtained from step 2 and the ratio (\(C\)) is obtained.

Step 4 The data points of (\(\frac{{\sigma }_{\mathrm{hmin}}}{{P}_{\mathrm{a}}}\),\(C\)) are plotted as point diagrams and then a suitable curve is fitted to the data points using the least squares method. The equations of the fitted curves are the same as the modified coefficients of the analytical formulas.

The equations of the modified coefficients C1, C2, C3, C4, C5, and C6 are shown in Table 14 and Fig. 14.

Table 14 Modification coefficients for KGD and PKN analytical formulas
Full size table
Fig. 14
figure 14

Coefficients for considering the effect of initial in-situ stresses a \({C}_{1}\) b \({C}_{2}\) c \({C}_{3}\) d \({C}_{4}\) e \({C}_{5}\) f \({C}_{6}\)

Full size image

These coefficients which are functions of in-situ minimum horizontal stress, should be substituted for original constant values in KGD and PKN models. In Table 14, \({P}_{\mathrm{a}}\) is the atmospheric pressure (Pa = 100 kPa).

Summary and conclusions

One of the common methods for stimulating oil and gas wells and increasing their productivity is Hydraulic Fracturing. For design of Hydraulic Fracturing operations, a number of models have been developed. KGD and PKN are the mostly used models in this area. In both models the characteristics of the fracture namely: fracture length (\(L\)), fracture width (\(w\)), and fluid pressure at the crack opening (\(p\)) are determined based on closed form relations. Because of the simplifying assumptions that are used in KGD and PKN models, the effects of some geomechanical and hydraulic factors are ignored in these classical analytical formulae. In this paper, by conducting a series of numerical simulations using XFEM, the accuracy of the analytical formulas of KGD and PKN models were examined for cases in which the reservoir in-situ stresses are taken into account. The results demonstrate that incorporating the in-situ stresses in the analysis of HF can greatly change the characteristics of the induced cracks. In this paper based on the results of the numerical simulations, the effect of in-situ stresses on Hydraulic Fracture characteristics (width, length, and crack mouth fluid pressure) are taken into account using a series of correction factors that should be applied to the original KGD and PKN formulae to modify them and provide more accurate predictions. Based on the obtained results the following conclusions can be drawn:

  1. 1.

    The results of the numerical simulations showed that when the minimum horizontal in-situ stress is kept constant and the maximum horizontal stress is changed, the values of the crack width, crack length and mouth pressure are almost constant.

  2. 2.

    When the maximum horizontal in-situ stress is kept constant and the minimum horizontal in-situ stress changes, the values of the crack width, crack length and mouth pressure change significantly.

  3. 3.

    The minimum horizontal in-situ stress (\({\sigma }_{\mathrm{hmin}}\)) is the most effective parameter in studying of the in-situ stresses effect on HF because the required fluid pressure for creating crack in rock is dependent on the in-situ-minimum principal stress.

  4. 4.

    This study showed that with increasing the stress ratio (\(\frac{{\sigma }_{\mathrm{hmax}}}{{\sigma }_{\mathrm{hmin}}}\)) from 1.0 to 3.0, the crack width decreases up to 1.5 times for KGD model and 2 times for PKN model, the crack length increases up to 1.4 times for KGD model and 3 times for PKN model, and the fluid pressure decreases up to 2.8 times for both models.

  5. 5.

    The obtained results and the proposed correction factors that modify the original formulae of PKN and KGD can be used in the design of HF process for the cases where the in-situ stresses in the reservoir is of importance.

Although the advantage of using the proposed correction factors in the design of hydraulic fracture operations is acquiring a more realistic and economical design, application of the results of the current study is limited to the cases where the real reservoir conditions are similar to the reservoir conditions considered in this study.