**2. Simulation of fluid flow and heat transfer**

Three distinct approaches exist in the literature to simulate the fluid flow in naturally fractured reservoirs namely: single continuum, dual continuum and discrete fracture approach. In single continuum, the fractured medium is represented by an equivalent homogeneous system using a specific permeability tensor. In dual continuum approach the whole domain is divided into two interacting domains: fractures and matrix where by matrix (represented by sugar cubes) provides the storage and fractures (having regular pattern) the permeability. In discrete fracture approach, fractures are explicitly discretized in the domain. These approaches are briefly discussed below followed by the proposed methodology which is used in this study.

#### **2.1. Hybrid of single continuum and discrete fracture**

Different approaches have been used in the literature to incorporate the fractures into the flow modeling. Each of these techniques has its own drawbacks and benefits. In this study a hybrid methodology combining the single continuum and discrete fracture networks model is used to increase accuracy and efficiency of the fluid flow simulation. In the proposed methodology a threshold value is defined for the fracture length. Fractures which are smaller than the threshold value are used to generate the grid based permeability tensor using boundary element technique. Fluid flow simulation is carried out by using the single continuum approach in the nominated blocks. Fractures which are equal to and longer than the threshold value are explicitly discretized in the domain using appropriate elements and the fluid flow is modeled using the discrete fracture approach. Such an approach provides a more accurate and realistic framework to consider the effect of long fractures on the fluid flow in fractured medium.

#### **2.2. Domain discretization using the hybrid methodology**

In this study the medium and long fractures (*l ≥* 50m) are discretized using triangular elements and the contribution of flow by fractures (*l* < 50m) are taken into account by calculating permeability tensor for each discretized element. A schematic representation of the domain discretization for a fractured reservoir is shown in Fig 3 (a) and (b).

Permeability tensor for each block is expressed as:

where *σ<sup>y</sup>* '

576 Effective and Sustainable Hydraulic Fracturing

reservoir scale.

medium.

**2. Simulation of fluid flow and heat transfer**

**2.1. Hybrid of single continuum and discrete fracture**

**2.2. Domain discretization using the hybrid methodology**

discretization for a fractured reservoir is shown in Fig 3 (a) and (b).

is the effective normal stress exerted on the fracture surfaces, *k* is the spring constant,

*p* is the pore pressure, *Δ* is the characteristic height of the fracture, *δ<sup>y</sup>* is the displacement of the fracture surfaces, *E* is the Modulus of elasticity, *τn* is the shear stress exerted on the fracture surfaces, *τ*0 is the threshold stress requires to start the shear displacement of the fracture surfaces and u is the displacement. *As* mentioned above, the aperture distribution along the fracture surface is calculated based on an analytical methodology in which fracture geometry, stress distribution and fluid pressure inside the fracture are needed to be known as a priori. For this purpose a thermo-poro-elastic model is developed to simulate the fluid flow in the

Three distinct approaches exist in the literature to simulate the fluid flow in naturally fractured reservoirs namely: single continuum, dual continuum and discrete fracture approach. In single continuum, the fractured medium is represented by an equivalent homogeneous system using a specific permeability tensor. In dual continuum approach the whole domain is divided into two interacting domains: fractures and matrix where by matrix (represented by sugar cubes) provides the storage and fractures (having regular pattern) the permeability. In discrete fracture approach, fractures are explicitly discretized in the domain. These approaches are briefly discussed below followed by the proposed methodology which is used in this study.

Different approaches have been used in the literature to incorporate the fractures into the flow modeling. Each of these techniques has its own drawbacks and benefits. In this study a hybrid methodology combining the single continuum and discrete fracture networks model is used to increase accuracy and efficiency of the fluid flow simulation. In the proposed methodology a threshold value is defined for the fracture length. Fractures which are smaller than the threshold value are used to generate the grid based permeability tensor using boundary element technique. Fluid flow simulation is carried out by using the single continuum approach in the nominated blocks. Fractures which are equal to and longer than the threshold value are explicitly discretized in the domain using appropriate elements and the fluid flow is modeled using the discrete fracture approach. Such an approach provides a more accurate and realistic framework to consider the effect of long fractures on the fluid flow in fractured

In this study the medium and long fractures (*l ≥* 50m) are discretized using triangular elements and the contribution of flow by fractures (*l* < 50m) are taken into account by calculating permeability tensor for each discretized element. A schematic representation of the domain

$$\mathbf{K} = \begin{bmatrix} k\_{xx} & k\_{xy} \\ k\_{yx} & k\_{yy} \end{bmatrix} \tag{13}$$

Permeability tensors are calculated by simulating fluid flow in individual fractures in each element. The concept of permeability tensor was first introduced by [18] by considering a set of parallel fractures in a Representative Elementary Volume (REV) with zero matrix permea‐ bility [18]. In another attempt [19] developed a methodology for calculation of permeability tensor for arbitrary oriented fractures using superposition technique [19].

**Figure 3.** Domain discretization by using the hybrid of the single continuum and discrete fracture approach. (a) frac‐ tures equal to and longer than 50 m are explicitly discretized in the reservoir domain by using the triangular elements. (b) after the discretization of the long fractures, the effect of short fractures (<50m) are taken into account by calcula‐ tion of the permeability tensor of the corresponding blocks which are cut by the fractures.

In this study the authors have considered interconnected fractures with fracture surface as infinite plate without roughness. In another approach [20] estimated permeability tensor by assuming fractures as a planar sink/source term [20]. Also [21] extended the approach and studied the effect of vertical fracture/ matrix permeability ratio on the permeability tensor. In a separate study, [22] used a numerical technique (BEM) to calculate the permeability tensor of the REV containing medium sized fractures considering fractures as a sink/source term [22]. Following this work [23] presented an analytical model to calculate the permeability tensor of the blocks containing infinite parallel fracture sets [23]. Also [24] improved the efficiency of their previous approach by considering the effect of short fractures using the analytical method proposed by [24]. In another approach [25] presented the first comprehensive methodology to calculate the permeability tensor for arbitrary oriented fractures in different length scales. In this study permeability tensor was determined by discretizing the solution domain into different subdomains depending on the length of the fractures using BEM [25]. Short fractures are considered as part of matrix porosity to improve the matrix permeability inhomogeneity. However, medium and long fractures are discretized explicitly in the domain and fluid flow is simulated using BEM. Then [26] extended [25] by increasing the efficiency of the BEM so that fluid flow in greater number of fractures can be simulated. The authors also presented for the first time effective permeability tensor calculation for the fractured REV by using the BEM. The effective permeability model was validated using laboratory derived data.

with the mass conservation equation to consider the effect of the fracture on the flow of fluid in the region close to the fractures. Size of this region depends on the size of the fracture. Also

A New Approach to Hydraulic Stimulation of Geothermal Reservoirs by Roughness Induced Fracture Opening

() 0 *<sup>f</sup> f ff <sup>p</sup> k Qq L L*

++=

fractures which are considered as part of the matrix porosity, the Laplace equation is solved

Where, pmi is the matrix pressure and pfi is the fracture pressure at the matrix/fracture interface and vmi is the normal fluid velocity at the ith fracture node along the fracture surface. Since the pressure on the matrix fracture interface is unknown, periodic boundary condition is applied

Fluid flow in long fractures (l>50m) is coupled with discretized element based permeability

Different numerical techniques have been used to model thermo-poro-elastic phenomena in fractured porous media. To have a detailed understanding of the complex geomechanical aspects of the fractured rocks and the induced perturbation, such as thermal drawdown caused by the cold injection fluid in geothermal reservoirs an appropriate numerical technique should be used which is capable of (a) adequately applying the boundary and initial conditions and (b) accurately representing the system geometry. In order to take the aforementioned issues

Weighted residual method and the Green's theorem are applied to discretize the mass, momentum and energy conservations equations [28]. As mentioned before, the finite element method is used in this study for the numerical simulation purpose. Therefore the state variables namely: displacement, pore pressure and temperature are defined using proper shape

tensor in poro-thermo-elastic environment by using local-thermal non-equilibrium.

is the fluid pressure inside the fractures and p is the pore fluid pressure. For the short

*V KP f ff* =- Ñ (14)

http://dx.doi.org/10.5772/56447

579

*mi fi p p* = (16)

*mi fi v v* = (17)

¶ ¶ (15)

fluid flow simulation in the matrix (region 1) is described as follow:

¶ ¶

where pf

using the following boundary conditions:

in an iterative scheme to calculate the pressure values.

**2.4. Reservoir scale fluid flow simulation**

into account, FEM is used in the current study.

functions as:

#### **2.3. REV discretization for permeability tensor calculation**

To calculate the effective permeability tensor, the fractured REV is divided into three distinct regions: matrix (region 1), fracture (region 2) and region around the fractures (region 3) as shown in Fig. 4.

**Figure 4.** Domain discretization based on different fracture lengths

Flow inside the fractures (region 2) is modeled using the cubic law. With the assumption of smooth fracture surfaces, cubic law can accurately simulate the flow inside the fractures [19, 27]. In matrix regions close to the fractures (region 3), the Darcy equation Eq. (14) is coupled with the mass conservation equation to consider the effect of the fracture on the flow of fluid in the region close to the fractures. Size of this region depends on the size of the fracture. Also fluid flow simulation in the matrix (region 1) is described as follow:

$$V\_f = -K\_f \nabla P\_f \tag{14}$$

$$\frac{\partial}{\partial \mathcal{L}} (\mathbf{k}\_f \frac{\partial p}{\partial \mathcal{L}}) + \mathcal{Q}\_f + q\_{ff} = 0 \tag{15}$$

where pf is the fluid pressure inside the fractures and p is the pore fluid pressure. For the short fractures which are considered as part of the matrix porosity, the Laplace equation is solved using the following boundary conditions:

$$p\_{mi} = p\_{ft} \tag{16}$$

$$
\boldsymbol{\upsilon}\_{mi} = \boldsymbol{\upsilon}\_{ft} \tag{17}
$$

Where, pmi is the matrix pressure and pfi is the fracture pressure at the matrix/fracture interface and vmi is the normal fluid velocity at the ith fracture node along the fracture surface. Since the pressure on the matrix fracture interface is unknown, periodic boundary condition is applied in an iterative scheme to calculate the pressure values.

#### **2.4. Reservoir scale fluid flow simulation**

proposed by [24]. In another approach [25] presented the first comprehensive methodology to calculate the permeability tensor for arbitrary oriented fractures in different length scales. In this study permeability tensor was determined by discretizing the solution domain into different subdomains depending on the length of the fractures using BEM [25]. Short fractures are considered as part of matrix porosity to improve the matrix permeability inhomogeneity. However, medium and long fractures are discretized explicitly in the domain and fluid flow is simulated using BEM. Then [26] extended [25] by increasing the efficiency of the BEM so that fluid flow in greater number of fractures can be simulated. The authors also presented for the first time effective permeability tensor calculation for the fractured REV by using the BEM.

To calculate the effective permeability tensor, the fractured REV is divided into three distinct regions: matrix (region 1), fracture (region 2) and region around the fractures (region 3) as

Flow inside the fractures (region 2) is modeled using the cubic law. With the assumption of smooth fracture surfaces, cubic law can accurately simulate the flow inside the fractures [19, 27]. In matrix regions close to the fractures (region 3), the Darcy equation Eq. (14) is coupled

The effective permeability model was validated using laboratory derived data.

**2.3. REV discretization for permeability tensor calculation**

**Figure 4.** Domain discretization based on different fracture lengths

shown in Fig. 4.

578 Effective and Sustainable Hydraulic Fracturing

Fluid flow in long fractures (l>50m) is coupled with discretized element based permeability tensor in poro-thermo-elastic environment by using local-thermal non-equilibrium.

Different numerical techniques have been used to model thermo-poro-elastic phenomena in fractured porous media. To have a detailed understanding of the complex geomechanical aspects of the fractured rocks and the induced perturbation, such as thermal drawdown caused by the cold injection fluid in geothermal reservoirs an appropriate numerical technique should be used which is capable of (a) adequately applying the boundary and initial conditions and (b) accurately representing the system geometry. In order to take the aforementioned issues into account, FEM is used in the current study.

Weighted residual method and the Green's theorem are applied to discretize the mass, momentum and energy conservations equations [28]. As mentioned before, the finite element method is used in this study for the numerical simulation purpose. Therefore the state variables namely: displacement, pore pressure and temperature are defined using proper shape functions as:

$$
\mu = \mathcal{N}\_u \overline{u} \tag{18}
$$

**4. Results and discussions**

**Mean**

**Fracture set**

**Distributio n Law**

discrete fracture map is shown in Table 1.

**Half-Width**

The proposed methodology is used to generate the discrete fracture map of the Soultz geothermal reservoir at the depth of 3650 m. the statistical parameters used to generate the

A New Approach to Hydraulic Stimulation of Geothermal Reservoirs by Roughness Induced Fracture Opening

**Half-Width**

F1 Normal 2 16 normal 70 7 NW 1.3E-7 187 6E-6 F2 Normal 162 19 normal 70 7 NE 3E-9 150 6E-6 F3 Normal 42 6 normal 74 3 NW 1.76E-8 95 4E-6 F4 Normal 129 6 normal 68 3 SW 3.3E-8 112 2E-6 F5 Uniform 0 180 normal 70 9 - 1E-8 100 5E-7

The discrete fracture map, the corresponding mesh generated for the reservoir domain and the permeability tensors for each triangular element (a sample region which is cut by a fracture

Also the reservoir properties used for the stimulation purpose are shown in Table 2. The reservoir is pressurized by injecting fluid through the injection well (GPK2). The pressurization was carried out over a period of 52 weeks. During the pressurization, the change in fracture width for each individual natural fracture and the resulting permeability tensor were calcu‐ lated. Following stimulation of the reservoir, a flow test was carried out over a period of 14 years. During the flow test, changes in fracture apertures due to thermo-poro-elastic stresses and the consequent changes in permeability were determined. Also estimated were the thermal drawdown, produced fluid temperature and production rate of the Soultz EGS.

Results of shear dilation are presented as average percentage increase in fracture aperture (see Fig. 6). From Fig. 6, it can be seen that there exists three distinct aperture histories: 0-40 weeks, 40-50 weeks and 50 weeks and above. Until about 40 weeks, a slow but linear increase in occurrence of dilation events due to induced fluid pressure of 51.7 MPa (bottom hole) and reaches a value of about 18% (average increase in aperture). Following this time, the rate of occurrence of dilation events increases sharply until about 50 weeks, thus reaching 60% increase in average fracture aperture. After which, no significant dilation events can be observed (a plateau of events is reached). When compared with previous study [29], in which shear dilation events are estimated based on a semi-empirical model (Willis-Richards et al,

**Dip Directi on**

**Fracture No.**

**Radius (m)**

http://dx.doi.org/10.5772/56447

**Transmissivit y (m2\s)**

581

**Azimuth Dip**

**Mean**

**Distribu tion Law**

**Table 1.** Statistical data used for the discrete fracture network generation. After [32]

of length<50m) are shown in Fig. 5 (a), (b) and (c) respectively.

$$p = N\_p \overline{p} \tag{19}$$

$$T = N\_T \overline{T} \tag{20}$$

Where N is the corresponding shape function and *u*¯, *p*¯ and *t* ¯ are the nodal values of the corresponding state variable. By applying the Galerkin's method and replacing the weighting functions by the corresponding variables' shape functions, the discretized form of the conser‐ vation equations can be written as follow [29, 30]:

$$(K + \frac{G}{3})\nabla(\nabla \cdot \mathbf{u}) + G \cdot \nabla^2 \mathbf{u} - \alpha \nabla p - \gamma\_1 \nabla T\_m = 0 \tag{21}$$

$$a(\nabla \cdot \mathbf{u}) + \beta \not{p} - \frac{k}{\mu} \nabla^2 p - \gamma\_2 \not{T} = 0 \tag{22}$$

$$\stackrel{\bullet}{T} + v(\nabla T) - c^T \nabla^2 T = 0\tag{23}$$

where, K is the bulk modulus of elasticity, G is the shear modulus, *γ*1 and *γ*<sup>2</sup> are the thermal expansion coefficient of the fluid and solid respectively; k is the permeability Tm is the matrix temperature, T is the fluid temperature and *μ* is the fluid viscosity.

#### **3. Fracture network generation**

Simulation of naturally fractured reservoirs offers significant challenges due to the lack of a methodology that can utilize field data. To date several methods have been proposed in the literature to characterize naturally fractured reservoirs. In this study a hybrid tectonostochastic simulation is proposed to characterize a naturally fractured reservoir [31]. A finite element based model is used to simulate the tectonic event of folding and unfolding of a geological structure. A nested neuro-stochastic technique is used to develop the interrelationship between different sources of data (seismic attributes, borehole images, core description, well logs etc.) and at the same time the sequential Gaussian approach is utilised to analyze field data along with fracture probability data. This approach has the ability to overcome commonly experienced discontinuity of the data in both horizontal and vertical directions.
