Open access peer-reviewed conference paper

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

Written By

Bisheng Wu, Xi Zhang, Andrew Bunger and Rob Jeffrey

Submitted: July 26th, 2012 Reviewed: March 13th, 2013 Published: May 17th, 2013

DOI: 10.5772/56373

Chapter metrics overview

2,553 Chapter Downloads

View Full Metrics


In order to understand the thermal recovery Behavior of an engineered geothermal system (EGS), this paper develops a model in which fluid circulates in a single, planar hydraulic fracture with a constant hydraulic aperture via multiple recharging and discharging wells. The coupled equations for heat convection in the fracture plane and heat transfer into the rock are provided for steady and irrotational fluid flow conditions. By using velocity potentials and streamline functions, the temperature along a streamline is found to be only a function of the potential. By utilizing the Laplace transformation, the analytical solutions in the Laplace space for the temperature field are found, which are numerically inverted for time-domain results. Several examples with different arrangements of injection and production wells are investigated and the comparison with other published results is provided. The semi-analytical results demonstrate that the proposed model provides an efficient and accurate approach for predicting the temperatures of a multi-well reservoir system.

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 2500C. One cubic kilometer of hot granite at 2500C 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 T0. When time t>0, a cold fluid with a constant injection rate Qi(i=1, 2...M) and constant temperature Tiniis injected from the injection points (xi, yi) (i=1,..M) and the heated fluid is pumped out at a constant production rate Qi(i=M+1, M+2...N), from the discharge wells (xi, yi) (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 xand y. Therfore, the governing equations for the fluid flow obeys the potential flow theory


where is the gradient operator, qfdenotes the discharge vector and ϕfis the velocity potential for the fluid.

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


where ·is the divergence operator, δdenotes the Dirac delta function, Qjis flow rate into or out of the jth well (xj, yj) (Qj>0 for injection and Qj<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


where Tfand Trdenote the fluid temperature and rock temperature, respectively, ρwand cware the mass density and specific heat, respectively, of the fluid, and λris the thermal conductivity of the rock.

For the surrounding rock, the heat conduction equation is


where κr= λr/ρrcr, ρrand crare 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,   E5

and the initial velocity potential in the fracture is

ϕf=ϕ0   on z=0.E6

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

Tr=Tf  on z=0.E7

The temperature at each injection well is fixed

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

where zdenotes 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.


4. Dimensionless formulation

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

Θrτ=2ΘrZ2, ·Qf=j=1NΩjδ(XXj,YYj),Qf=Φf,    Qf·Θf+χΘrz=0, on Z=0,Θf=Θinj,   on (X,Y)=(Xj,Yj), Θf=Θr   on Z=0,Θr=Θf=0,   on τ=0,E10

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 Lis 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


and the flow stream function ψ is expressed as


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


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

v2ΘfΦ+χΘrZ=0,  Z=0,E14



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 Θ^fis obtained as


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


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


We note that in equation (16), Θ incan 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


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.

Number of injection wellsMRock ther. conductivity λr(W/(m∙K))2.2
Number of total wellsNRock mass density ρr(Kg/m3)2700
Total injection rate Qtotal(m3/s)0.06Rock specific heat cr(J/(kg∙K))790
Total Extraction rate Qtotal(m3/s)0.06Fluid mass density ρw (Kg/m3)900
Injection temperature Tin90Fluid specific heat cw(J/(kg∙K))4200
Initial temperature T0 (º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 timet=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 ratesQ=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 cw Q) is an important parameter, where Lis the well separation and Qis 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 Qin=Q/4=0.015 m3/s and Qout=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 cw 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.


  1. 1. LauwerierH. AThe transport of heat in an oil layer caused by the injection of hot fluid, Applied Scientific Research,19555145150
  2. 2. BödvarssonG. SOn the Temperature of Water Flowing through Fractures. Journal of Geophysical Research,1969,74(8),1987-1992.
  3. 3. GringartenA. CTheory of heat extraction from fractured hot dry rock, Journal of Geophysical Research,1975
  4. 4. ZhangXJeffreyR. GWuBTwo basic problems for hot dry rock reservoir stimulation and production. Australian Geothermal Conference2009Brisbane, QLD, Australia,1113Nov., 2009.
  5. 5. MossopASeismicity, subsidence and strain at the Geysers geothermal field. PhD. Dissertation, Stanford University,2001
  6. 6. GhassemiATarasovsSandChengA. H. DAn integral equation solution for three-dimensional heat extraction from planar fracture in hot dry rock. International Journal for Numerical and Analytical Methods in Geomechanics,2003279891004
  7. 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. 8. GringartenA. CRameyH. JThe use of source and Green’s functions in solving unsteady-flow problems in reservoirs. Society of Petroleum Engineers Journal,1973285296
  9. 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,196052524536
  10. 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. 11. MuskatMThe flow of homogeneous fluids through porous media. MicGraw-Hill Book Company, Inc.,1937
  12. 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. 13. RodemannHAnalytical model calculations on heat exchange in a fracture. Urach Geothermal Project, Haenel R. (ed.), Stuttgart,1982
  14. 14. SchulzRAnalytical model calculations for heat exchange in a confined aquifer. J. Geophys. Int.,1997611220
  15. 15. HeuerNKopperTandWindelbergDMathematical model of a Hot Dry Rock system. Geophys. J. Int.,1991105659664
  16. 16. OginoFYamamuraMFukudaTHeat transfer from hot dry rock to water flowing through a circular fracture. Geothermics,1999282144
  17. 17. Harlow, F. H, & Pracht, W. Theoretical study of geothermal energy extraction. Journal of Geophysical Research, (1972).,77(35), 7038-7048.
  18. 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. 19. GhassemiATarasovsSandChengA. H. DIntegral equation solution of heat extraction-induced thermal stress in enhanced geothermal reservoirs, International Journal for numerical and Analytical Methods in Geomechanics,200529829844

Written By

Bisheng Wu, Xi Zhang, Andrew Bunger and Rob Jeffrey

Submitted: July 26th, 2012 Reviewed: March 13th, 2013 Published: May 17th, 2013