InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Energy Engineering » "Effective and Sustainable Hydraulic Fracturing", book edited by Andrew P. Bunger, John McLennan and Rob Jeffrey, ISBN 978-953-51-1137-5, Published: May 17, 2013 under CC BY 3.0 license. © The Author(s).

# An Efficient and Accurate Approach for Studying the Heat Extraction from Multiple Recharge and Discharge Wells

By Bisheng Wu, Xi Zhang, Andrew Bunger and Rob Jeffrey
DOI: 10.5772/56373

Article top

## Overview

Figure 1. (a) Multi-well geothermal reservoir model, (b) Geometry for Habanero project with 5 injection wells and 4 production wells (red denotes production well, blue denotes injection wells).

Figure 2. Case studies with different injection wells and extraction wells. Blue circle denotes injection well and red circle denotes production well.

Figure 3. Comparison between the present approach and the analytical solutions: (a) Temperature profile along the radial direction (physical time t=115 days, 3.17 years and 31.7 years) for the case with only a single injection well at the origin. Parameters are listed in the paper by Ghassemi et al. [6]; (b) Temperature distribtion along the line connecting a recharge well (left) and a discharge well (right). Blue circle denotes injection well and red circle denotes extraction well.

Figure 4. Normalized change in production temperature from cases (b) to (e) under the same total injection and production rates Q=0.06m3/s. The average value of the temperatures at all the points on the small circle around the production well is chosen as the production temperature.

Figure 5. Normalized temperature change for the production wells while keeping the controlling parameter set to χ=4.76e-3.

Figure 6. Isothermal lines of the normalized temperature change at (a) time=5 years (τ=5.77e-4) and (b) time=21 years (τ=2.42e-3) for case (d). Blue circle denotes injection well and red circle denotes production well.

Figure 7. (a) Streamlines profiles and (b) normalized production temperature change in case (6) when Q in =Q/4=0.015 m3/s and Q out =Q/5=0.012 m3/s. Plot (a) shows a quarter of the whole streamline profiles which are symmetric across the plot axes.

# An Efficient and Accurate Approach for Studying the Heat Extraction from Multiple Recharge and Discharge Wells

Bisheng Wu1, Xi Zhang1, Andrew Bunger1, 2 and Rob Jeffrey1

## 1. Introduction

Heat contained in the upper 10 km of the Earth’s crust represents a large, accessible, low-emission energy source that can substitute for other energy sources that produce significantly more greenhouse gas. Geothermal energy has become one of the most promising energy alternatives in the future. For example, Australia has a large volume of identified high heat producing granites within 3 to 5km of the surface. In some places at 5kms the temperature is more than 2500 C. One cubic kilometer of hot granite at 2500 C has the stored energy equivalent of 40 million barrels of oil (Geodynamics website). Therefore, developing methods to capture this resource of clean energy will help realize the potential of commercial EGS. In addition to the physical experiments and field exploitation, theoretical studies such as mathematical modeling are important in developing ways to maximize heat extraction from geothermal reservoirs treated by hydraulic fracturing.

Mathematical modeling of thermal problems associated with oil recovery or geothermal heat extraction from reservoirs containing hydraulic or natural fractures has been studied extensively. The models can be classified into three types based on the dimension of the heat transfer problem. The first type is based on one-dimensional (1-D) heat diffusion in the fracture and it mainly contains two cases. The first one considers 1D fluid flow in a 2D fracture. In this case, the fluid velocity is uniform and no singular point in pressure exists at the injection or production well [1-3]. The second one considers fluid flow that is radial and axisymmetric, for example the fluid is injected into a single well and produced from a ring of production wells that are all at the same radial distance from the injection well [4]. Due to the symmetry, the problem can be treated as one-dimensional based on the radial distance. The fluid velocity can be obtained either by considering the effect of the wellbore size or neglecting it. When the injection or pumping at the well is regarded as a point source or sink, the velocity will have a singularity of 1/r theoretically at the injection and production points. Fortunately, this singularity is overcome by the geometrical symmetry [5, 6] for the case with a single well.

The second type is so-called 21/2 dimensional model [7] which results from 2-D heat transfer within fractures and 1-D heat transfer within the adjacent rocks, it becomes more difficult to find analytical solutions for these cases because of the coupling of the steady fluid flow and heat transfer. There are mainly two ways to obtain the fluid velocity analytically for the cases with single, dipole or multiple wells; one method calculates the pressure by using a Green’s function [8], and the other method calculates the velocity potential and stream functions by using the source solution for the 2-D Laplace equation [9, 10]. From the perspective of modeling the long lifetime of a geothermal reservoir, the above two approaches are functionally the same. Although the first method can obtain accurate results for the fluid pressure and thus the fluid velocity, there still exists two coordinates in the governing equation and the difficulty of solving this equation is not reduced. In addition, the singularity issues make it invalid to use the boundary condition (i.e. injection temperature) when solving the equations. For this case, solutions still require the use of numerical methods.

In order to overcome the above issues including the computational difficulty and singularity at the source or sink points, Muskat [11] proposed a method in which an orthogonal set of curvilinear coordinates is introduced that correspond to the permanent set of equipotential and streamline surfaces. By using this transformation, the terms involving cross products related to the fluid diffusion and heat transfer are greatly simplified. Gringarten and Sauty [12] were the first to apply this concept to solve the dipole well problem of fluid flow and heat diffusion in an infinite fracture by using velocity potentials and stream functions based on the results of Dacosta and Bennett [9]. In their model, each channel leaving a particular injection well is treated separately and the two-dimensional heat diffusion problem in the fracture plane is simplified to be one-dimensional. Later Rodemann [13], Schulz [14], Heuer et al. [15] and Ogino et al. [16] extended Gringarten and Sauty‘s [12] approach to the problems related to heat extraction from a finite fracture or multiple fractures with a recharge-discharge well pair and obtained good analytical solutions.

The third type is the complete 3D problem for the heat transfer in the fracture and the rock. In these more complicated cases, numerical schemes, such as finite element [7], Marker-and-Cell method [17] and finite difference approaches [18, 19], have to be used to obtain the thermal behavior of the fluid and reservoir. Normally, finding the solution is computationally demanding.

Besides applying the method to a particular set of boundary conditions for EGS reservoirs, the advance of the present model is in that the semi-analytical solutions for the rock and fluid temperatures are obtained, thus significantly reducing computational cost. It also provides an efficient and accurate way to further study the effect of some factors, such as the number of wells, the well spacing, the injection rates and production rates, to maximize the heat extraction from the EGS reservoirs.

## 2. Problem description

The geometry considered by the present model is shown in Fig. 1. There exist M (M≥1) vertical recharge well and N-M (N≥2) vertical discharge wells, which intercepts the fracture containing the injected fluid. The whole system (liquid and rock formation) is initially in an equilibrium state with the uniform temperature T 0. When time t>0, a cold fluid with a constant injection rate Q i (i=1, 2...M) and constant temperature Tini is injected from the injection points (x i, y i ) (i=1,..M) and the heated fluid is pumped out at a constant production rate Q i (i=M+1, M+2...N), from the discharge wells (x i, y i ) (i=M+1, 2...N).

It is assumed that after a hydraulic fracturing treatment, there is a connected flat fracture plane intercepted by all wells. The geometry of the fracture is defined by its radius In particular, the fracture radius is infinite in this paper.

Some assumptions are made for the present model:

1. The fluid is incompressible with Newtonian rheology and the rock is impermeable, homogenous and isotropic;

2. The stress-induced change of the fracture aperture is ignored;

3. The material properties (density, specific heat capacity and thermal conductivity) of the rock and the fluid are constant and independent of the temperature;

4. The heat condution in the fluid is neglected when the injection rate is sufficiently great.

#### Figure 1.

(a) Multi-well geothermal reservoir model, (b) Geometry for Habanero project with 5 injection wells and 4 production wells (red denotes production well, blue denotes injection wells).

## 3. Governing equations

### 3.1. Pressure diffusion in planar fractures

Because the fluid flow along the narrow fracture is much faster than the reservoir heat temperature changes, it can be assumed to be a pseudo-steady state. This means that the fluid pressure or velocity potential is a function of only the in-plane coordinates x and y. Therfore, the governing equations for the fluid flow obeys the potential flow theory

 qf=−∇ϕf, (1)

where is the gradient operator, q f denotes the discharge vector and ϕf is the velocity potential for the fluid.

In the present model, there are a total of M recharge wells and N-M production wells to be considered. The continuity equation for the flow of an incompressible Newtonian fluid in the fracture plane is expressed as

 ∇·qf=∑j=1NQjδ(x−xj,y−yj), (2)

where · is the divergence operator, δ denotes the Dirac delta function, Q j is flow rate into or out of the jth well (x j , y j ) (Q j >0 for injection and Q j <0 for production).

### 3.2. Heat transport in closed fractures

According to Cheng et al. [18], the heat transport in the fractures can be described by the following equation when the effect of the heat storage and longitudinal dispersion can be ignored

 ρwcwqf·∇Tf+2λr∂Tr∂z=0, (3)

where T f and T r denote the fluid temperature and rock temperature, respectively, ρ w and c w are the mass density and specific heat, respectively, of the fluid, and λ r is the thermal conductivity of the rock.

For the surrounding rock, the heat conduction equation is

 κr∇2Tr=∂Tr∂t, (4)

where κ r= λ r/ρ r c r , ρ r and c r are the mass density and specific heat, respectively, of the rock formation.

### 3.3. Initial and boundary condition

The initial temperature distribution for the rock and fluid is assumed to be a constant

 Tf=Tr=T0, (5)

and the initial velocity potential in the fracture is

 ϕf=ϕ0   on z=0. (6)

It must be mentioned that for the temperature continuity along the fracture, we have

 Tr=Tf  on z=0. (7)

The temperature at each injection well is fixed

 Tf=Tinj,   at (xj,yj) (j=1,..M), (8)

where z denotes the vertical coordinate with z=0 being the fracture plane.

In order to avoid the existence of the volume associated with the infinite boundary of the fracture plane, the total injection rate is always equal to the total production rates, which is reasonable for long-term estimatation, i.e.

 ∑j=1NQj=0. (9)

## 4. Dimensionless formulation

After some manipulation, the above equations can be written in a simplified form

 ∂Θr∂τ=∂2Θr∂Z2, ∇·Qf=∑j=1NΩjδ(X−Xj,Y−Yj),Qf=−∇Φf,    Qf·∇Θf+χ∂Θr∂z=0, on Z=0,Θf=Θinj,   on (X,Y)=(Xj,Yj), Θf=Θr   on Z=0,Θr=Θf=0,   on τ=0, (10)

based on the following transformation

Θr=TrT0T0,  Θf=TfT0T0,  Θinj=TinjT0T0,  Φf=ϕfϕ0Q, Ωj=QjQ,Qf=LQqf, τ=tκrL2,  X=xL, Y=yL, Z=zL, χ=2λrLρwcwQ,

where L is the characteristic length to be chosen from the geometrical configuration of the wells (for example, the minimum distance between wells).

## 5. Velocity potential and stream function

According to the superposition principle, the flow velocity potential Ф for a sum of the injection and producing wells is

 Φ=−14π∑j=1NΩjln[(X−Xj)2+(Y−Yj)2], (11)

and the flow stream function ψ is expressed as

 Ψ=−12π∑j=1NΩjarctan(Y−Yj)(X−Xj). (12)

The following formulas are valid when the variables Ф and ψ denote the velocity potential function and stream function, respectively,

 VX=−∂Φ∂X=−∂Ψ∂Y, ∂Θf∂X=∂Θf∂Φ(−VX)+∂Θf∂ΨVY,VY=−∂Φ∂Y=∂Ψ∂X,    ∂Θf∂Y=∂Θf∂Φ(−VY)+∂Θf∂Ψ(−VX). (13)

Then equation (10) 3 after using the above identities, becomes

 −v2∂Θf∂Φ+χ∂Θr∂Z=0,  Z=0, (14)

where

v2=VX2+VY2.

Appling Laplace transform to equation (14) and making use of the general Laplacian solution of equation (10) for Θ^r , the analytical solution for the Laplace tranform of the fluid temperature Θ^f is obtained as

 Θ^f=Θinse−χsg(Φ,Ψ), (15)

where the cap ^ denote the Laplace tranform and the function g(Ф,ψ) is defined as

g(Φ,Ψ)=Φ+dζv2(ζ,Ψ).

The solution for Θ f in the time domain is obtained by the inverse transformation of Eq. (15)

 Θf=Θinerfc[χg(Φ,Ψ)2τ]. (16)

We note that in equation (16), Θ in can be measured and χ is a constant. Therefore, the key to the solutions for the temperature is to calculate the value of the function g(Ф, ψ).

## 6. Methodology, verification and numerical results

Generally, when the number of the wells is larger than 2, the analytical solutions for g(Ф, ψ) are difficult and are not available in the literature. Then the following computational procedures are adopted. Start with the points which are very close to one of the injection wells and the velocity, velocity potential Ф0 and stream function ψ0 are calculated with the known coordinates. We then move along the streamline, and the velocity potential is increased by an increment ΔФ. In addition the coordinates of the next point on the same streamline can be evaluated via the following equations

 Ψ0=Ψ(X,Y),Φ0−ΔΦ=Φ(X,Y), (17)

which can be solved by using the Newton-Raphson method. In the same way, the information of the next point on the same streamline is obtained until the new point reaches a region which is regarded as being within the range of a production well or as infinity. There are no well defined criteria to define the region which is regarded as being within the range of a production well or as infinity. In the present calculations, we have compared cases using several small radial distances (such as 0.002m, 0.001m and 0.0001m) from the production wells and different far-field distances (such as 6km, 8km and 10km) from all the wells. When we find little difference between the results for these cases, the solutions are taken as being accurate. It should be noted that the position of the new point relative to the total wells should be compared with that of the starting point on the same streamline in each step. If the relative position changes, the stream function should be plus or minus PI which depends on the position change.

For verification and testing of the proposed model, numerical results for several well arrangements as shown in Fig. 2 are provided below. In particular, the required parameters are listed in Table 1.

 Parameters Value Parameters Value Number of injection wells M Rock ther. conductivity λ r (W/(m∙K)) 2.2 Number of total wells N Rock mass density ρr (Kg/m3) 2700 Total injection rate Qtotal (m3/s) 0.06 Rock specific heat cr (J/(kg∙K)) 790 Total Extraction rate Qtotal (m3/s) 0.06 Fluid mass density ρ w (Kg/m3) 900 Injection temperature T in 90 Fluid specific heat cw (J/(kg∙K)) 4200 Initial temperature T 0 (ºC) 270

#### Table 1.

Parameters for the present calculations

#### Figure 2.

Case studies with different injection wells and extraction wells. Blue circle denotes injection well and red circle denotes production well.

### 6.1. Cases with one single well and dipole wells

In Figure 3(a), the temperature evolution along the radial direction for the case with one single injection well is compared between the present approach and the analytical solution [5]. In Figure 3(b), the temperature profiles along the line connecting the injection well and the production well are also compared for the case with dipole wells[14] It is found that the results for the above two cases from the current method and analytical solutions are in excellent agreement.

### Figure 3.

Comparison between the present approach and the analytical solutions: (a) Temperature profile along the radial direction (physical time t=115 days, 3.17 years and 31.7 years) for the case with only a single injection well at the origin. Parameters are listed in the paper by Ghassemi et al. [6]; (b) Temperature distribtion along the line connecting a recharge well (left) and a discharge well (right). Blue circle denotes injection well and red circle denotes extraction well.

### 6.2. Cases with two injection wells or multiple wells

Figure 4 displays the extraction temperature change at different times for the cases (b) to (e) (Figure 2). As the geometrical configurations of all cases are symmetrical, the output temperature for only one production well in each case is plotted. The total injection and production rates in all cases are kept the same, i.e. 0.06 m3/s, and for the purpose of simplicity, the individual injection rate for each case is assumed to be the same. For example, when there are two injection/production wells, the individual injection/production rate is 0.06/2=0.03 m3/s. These results show that case (b) has the least temperature decrease, the reason being that the well separation is greater than the other cases. However, this greater separation and the fact that only 2 wells are involved will undoubtedly lead to higher pressure differences in order to obtain the same total flow rate. Case (d), with 2 injectors and 2 producers, is the next best performing well from the standpoint of minimizing temperature falloff. Interestingly, case (c) and case (e) perform almost equally, but case (c) only requires 3 wells while case (e) requires 5.

### Figure 4.

Normalized change in production temperature from cases (b) to (e) under the same total injection and production rates Q=0.06m3/s. The average value of the temperatures at all the points on the small circle around the production well is chosen as the production temperature.

Figure 5 shows the output temperature change for different geometrical sizes and injection rates are plotted for case (c) (Figure 2). Χ=2λ r L/(ρ w c w Q) is an important parameter, where L is the well separation and Q is the total flow rate. When the value of χ is kept constant, the fluid temperature at the output will evolve equivalently.

### Figure 5.

Normalized temperature change for the production wells while keeping the controlling parameter set to χ=4.76e-3.

Figure 6 presents isothermal lines at 5 years and 21 years for case (d). We can see clearly the contours at the same temperature propagating away from the injection well with time. As there are associated with the singularity in flow velocity at the production well whereby the temperature there cannot be obtained analytically and is calculated approximately. In order to make this approximation, a small circle around the well is chosen and calculations are not carried out inside this region. The number of the streamlines used in the computation which flows into this well can be counted and the temperature averaged over all used streamlines can also be computed and chosen as the production well temperature. In Figure 6, because of the symmetry, each production well accepts flow from each injection well. However, which injection well the influx at a specific part of a well comes from depends on the angle at which it faces toward an injection well. Therefore, the temperature fields around the production wells are very complicated and require an averaging method.

### Figure 6.

Isothermal lines of the normalized temperature change at (a) time=5 years (τ=5.77e-4) and (b) time=21 years (τ=2.42e-3) for case (d). Blue circle denotes injection well and red circle denotes production well.

### Figure 7.

(a) Streamlines profiles and (b) normalized production temperature change in case (6) when Q in =Q/4=0.015 m3/s and Q out =Q/5=0.012 m3/s. Plot (a) shows a quarter of the whole streamline profiles which are symmetric across the plot axes.

The streamline profiles and production temperature in Figure 7 are for case (f) (Figure 2) where there are four injection wells and five production wells. The total flow rate is also 0.06m3/s. Due to the symmetry, the temperature distributions of injection wells 1, 2, 3 and 4 are similar, as are the temperature fields between wells 5 and 9 and between production wells 6 and 8. Figure 7(b) shows the trend in the average normalized temperature variations in time. On the other hand, we can find that the contribution to the normalized temperature change of a particular well may strongly depend on the well arrangement, as well as on the imposed outflow rate at the production well, based on the temperature changes in Figure 7(b). For instance, all the four injection wells contribute to outflux at well 7 (the middle one in case (f)), but its output temperature undergoes the biggest decrease as show Figure 7(b). This is attributed to a larger flow rate in this central region, defined by four injection wells. The heat exchange from the rock to the fluid is not sufficient to heat the fluid to the same degree as occurs at other production wells. In the end, this will affect the output temperature at well 7. It must be remembered that all production wells are assumed to have the same outflux. Therefore, some fluid heated in the central region is forced to move outward and eventually reaches the outer production wells. The direct result of the heated fluid movement is to increase the output temperature at the outer production wells.

## 7. Conclusions

By using the Laplace transformation, the analytical form of the solutions in the Laplace space and thus, by inversion, in the time domain are obtained. The approach provides an efficient and accurate way to calculate the rock and fluid temperatures. Through several preliminary case studies, some brief conclusions are obtained:

1. Even for a multi-well geothermal reservoir, we obtain analytical or semi-analytical solutions that provide an efficient and accurate way for predicting the rock and fluid temperature;

2. The flow rates and the relative locations of the wells determine the flow path of the fluid. The injection rate also determines the thermal behavoir of the fluid and the parammeter Χ=2λ r L/(ρ w c w Q) reflects the rate of heat extracted by the cold injection fluid.

3. The efficiency of heat extraction from the EGS reservoir studied here depends on the layout and spacing of the injection and production wells. In fact, our model showed some contrasting cases where fewer wells outperformed a greater number of wells. Ongoing work will be aimed at finding optimal well layouts and flow rates.

## References

1 - H. A Lauwerier, The transport of heat in an oil layer caused by the injection of hot fluid, Applied Scientific Research, 1955 5 145 150
2 - G. S Bödvarsson, On the Temperature of Water Flowing through Fractures. Journal of Geophysical Research, 1969,74(8),1987-1992.
3 - A. C Gringarten, Theory of heat extraction from fractured hot dry rock, Journal of Geophysical Research, 1975
4 - X Zhang, R. G Jeffrey, B Wu, Two basic problems for hot dry rock reservoir stimulation and production. Australian Geothermal Conference 2009Brisbane, QLD, Australia, 11 13Nov., 2009.
5 - A Mossop, Seismicity, subsidence and strain at the Geysers geothermal field. PhD. Dissertation, Stanford University, 2001
6 - A Ghassemi, S Tarasovs, and A. H. D Cheng, An integral equation solution for three-dimensional heat extraction from planar fracture in hot dry rock. International Journal for Numerical and Analytical Methods in Geomechanics, 2003 27 989 1004
7 - Kolditz, O. Modelling flow and heat transfer in fractured rocks: conceptual model of a 3-D deterministic fracture network. Geothermics, (1995). , 24(3),421-437.
8 - A. C Gringarten, H. J Ramey, The use of source and Green’s functions in solving unsteady-flow problems in reservoirs. Society of Petroleum Engineers Journal, 1973 285 296
9 - DaCosta, J.A, & Bennett R.R., The pattern of flow in the vicinity of a recharging and discharging pair of wells in an aquifer having areal parallel flow. International Association of Scientific Hydrology, 1960 52 524 536
10 - Grove, D. B, Beetem, W. A, & Sower, F. B. Fluid travel time between a recharging and discharging well pair in an aquifer having a uniform regional flow field, Water Resources Research,(1970)., 6(5), 1404-1410.
11 - M Muskat, The flow of homogeneous fluids through porous media. MicGraw-Hill Book Company, Inc., 1937
12 - Gringarten, A. C, & Sauty, J. P. A Theoretical Study of Heat Extraction from aquifers with uniform regional flow. J. Geophys. Res., (1975)., 80(35),4956-4962.
13 - H Rodemann, Analytical model calculations on heat exchange in a fracture. Urach Geothermal Project, Haenel R. (ed.), Stuttgart, 1982
14 - R Schulz, Analytical model calculations for heat exchange in a confined aquifer. J. Geophys. Int., 1997 61 12 20
15 - N Heuer, T Kopper, and D Windelberg, Mathematical model of a Hot Dry Rock system. Geophys. J. Int., 1991 105 659 664
16 - F Ogino, M Yamamura, T Fukuda, Heat transfer from hot dry rock to water flowing through a circular fracture. Geothermics, 1999 28 21 44
17 - Harlow, F. H, & Pracht, W. Theoretical study of geothermal energy extraction. Journal of Geophysical Research, (1972).,77(35), 7038-7048.
18 - Cheng, A. H. D, Ghassemi, A, & Detournay, E. Integral equation solution of heat extraction from a fracture in hot dry rock. International Journal for Numerical and Analytical Method in Geomechanics, (2001)., 25, 1327.
19 - A Ghassemi, S Tarasovs, and A. H. D Cheng, Integral equation solution of heat extraction-induced thermal stress in enhanced geothermal reservoirs, International Journal for numerical and Analytical Methods in Geomechanics, 2005 29 829 844