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

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 2^{1}/_{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
*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:

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

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

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

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

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

where
**q**
_{
f
} denotes the discharge vector and

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

where
*δ* denotes the Dirac delta function, *Q*
_{
j
} is flow rate into or out of the *j*th 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

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

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

and the initial velocity potential in the fracture is

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

The temperature at each injection well is fixed

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.

## 4. Dimensionless formulation

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

(10) |

based on the following transformation

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

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

where

Appling Laplace transform to equation (14) and making use of the general Laplacian solution of equation (10) for

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

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

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

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.

### 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.

### 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 m^{3}/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 m^{3}/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 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 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.

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.06m^{3}/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:

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;

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