# An Approximate Solution for Predicting the Heat Extraction and Preventing Heat Loss from a Closed-Loop Geothermal Reservoir

Authors: Bisheng Wu; Tianshou Ma; Guanhong Feng; Zuorong Chen; Xi Zhang
Journal: Geofluids
 (2017)

---

## Abstract

Approximate solutions are found for a mathematical model developed to predict the heat extraction from a closed-loop geothermal system which consists of two vertical wells (one for injection and the other for production) and one horizontal well which connects the two vertical wells. Based on the feature of slow heat conduction in rock formation, the fluid flow in the well is divided into three stages, that is, in the injection, horizontal, and production wells. The output temperature of each stage is regarded as the input of the next stage. The results from the present model are compared with those obtained from numerical simulator TOUGH2 and show first-order agreement with a temperature difference less than 4°C for the case where the fluid circulated for 2.74 years. In the end, a parametric study shows that (1) the injection rate plays dominant role in affecting the output performance, (2) higher injection temperature produces larger output temperature but decreases the total heat extracted given a specific time, (3) the output performance of geothermal reservoir is insensitive to fluid viscosity, and (4) there exists a critical point that indicates if the fluid releases heat into or absorbs heat from the surrounding formation.


---

## Body

## 1. Introduction

With the rapid increasing demand for clean and renewable energy, geothermal energy has become one of the most promising alternatives for energy supplies due to its many advantages. First, geothermal energy is inexhaustible and renewable due to radioactive decay from below the earth [1] and there is a wide range of geothermal resources in the world including in the USA, Mexico, Germany, Italy, Finland, Norway, Iceland, Sweden, New Zealand, Australia, Indonesia, and China [1, 2]. Second, geothermal energy is much cleaner and environmentally friendly compared to the conventional fossil fuels. It has been reported that the geothermal energy based plants produce 1/12 carbon dioxide of the coal based plant [3]; Third, up to date, as only small amount of geothermal resource has been utilized for heating or power generation compared to its large reserve, there is a large potential for extraction.Closed-loop geothermal system represents an important mean to extract heat from below the earth and study on maximizing its output performance is very important [4, 5]. The basic principle of a closed-loop system is simple: a cold fluid (CO2 or water) is injected into the well, absorbs heat from hotter surrounding rocks when flowing in a channel such as wellbore or fracture, and then is pumped out from the production well. In the present paper, the connecting channel is chosen to be a wellbore. In addition to extracting heat from the reservoir by fluid circulation, the closed-loop system can prevent contaminants in fluids from leaking into the reservoir. As the output performance, including temperature and thermal power, of a geothermal reservoir is affected by many factors such as number and location of wells and fractures, flow rates, and velocity and direction of areal flow [6, 7], an efficient mathematical model is required for studying the relationship between reservoir performance and these parameters depending on specific system configuration.Fluid flow and heat transfer in a wellbore/pipe have received a great deal of attention during the last several decades due to their importance in oil and gas industry. For example, the thermoelastic stress change is significant near the wellbore; it may change the fracture initiation and path during hydraulic fracturing [10–12]. The mathematical models can be mainly classified into two types, that is, analytical and numerical. Many studies belong to the first type. For example, Bullard [13] and Moss and White [14] used a line/source model to study the time for a wellbore to attain temperature equilibrium and the temperature evolution in a water-injection well, respectively. Recently, the infinite, cylindrical, and finite line source theories [15] have been widely used for geothermal simulation. For example, Eskilson [16], Zeng et al. [17], Sutton et al. [18], Bandos et al. [19], Michopoulos and Kyriakis [20], Molina-Giraldo et al. [21], Hecht-Méndez et al. [4], Rivera et al. [22], and Zhou et al. [23] used line source models to predict the temperature evolution in the ground source heat pumps. In addition, Ramey [24] presented approximate solutions for the wellbore/reservoir (W/R) system involved in injection of hot or cold fluid. Dowdle and Cobb [25] applied the Horner temperature plot method, which is similar to conventional pressure build-up, to predict the static wellbore temperature from well logs. Edwardson et al. [26] and Tragesser et al. [27] studied the wellbore temperature during mud circulation by using an exact method which is based on the solution of differential equation of heat conduction. Because the above models do not consider the heat exchange between the fluid inside and that outside the drill pipe, they cannot be used as a general tool for predicting the thermal behaviour in a wellbore/reservoir systems.As for fluid circulation in a wellbore/reservoir system, the work by Raymond has to be mentioned [28]. In addition to providing systematic derivation of the governing equations characterizing the transient fluid flow and heat transfer, Raymond [28] obtained the analytical solutions for the case with steady state heat transfer in the wellbore by using Laplace transformation and also numerically solved the problem with transient heat transfer in both wellbore and formation by using the finite difference method (FDM). Most of the models studied later on fluid circulation in a wellbore/reservoir system more or less borrow some concepts from Raymond’s work [28]. For example, Holmes and Swift [29] solved approximately similar governing equations under assumption of steady state heat transfer between the annular fluid and the formation. Keller et al. [30] considered the effects of heat generated by fluid friction, sliding casing strings, and other energy sources in the system. García et al. [31], Espinosa-Paredes et al. [32], Fomin et al. [33], and Izgec et al. [34] presented a fully transient FDM model by considering the transient heat behaviour during fluid circulation in a W/R system. Wu et al. [35] solved analytically the transient heat transfer both in the wellbore and formation during drilling with fluid circulation in the wellbore.In addition to the above-mentioned approaches, recently a large number of numerical open sources or commercial software are available for geothermal simulation. For example, HYDROTHERM by US Geological Survey [36] and TOUGH2 by Lawrence Berkeley National Laboratory [37] based on FDM, OpenGeoSys by Helmholtz Centre of Environmental Research, Germany [38], ROCKFLOW by Kolditz et al. [39], FEFLOW by Diersch et al. [40, 41], FEHM by Los Alamos National Laboratory [42] and COMSOL [43] based on finite element method (FEM), ECLIPSE by Schlumberger [44], and Fluent by ANSYS [45] based on finite volume method (FVM) have been applied successfully to model the heat extraction from geothermal reservoirs with complicated geometries. In addition, the boundary element method (BEM) was also used for modelling similar problems. For example, Ding [46], Ghassemi et al. [47], Kumar and Gutierrez [48], and McClure and Horne [49] used the BEM or displacement discontinuity element method (DDM) to simulate the heat and mass transfer in geothermal reservoir with single fracture or fracture networks.Although the above numerical methods can be used to solve complex problems, they sometimes are computationally costly. The analytical methods cannot be directly used to solve the present closed-loop geothermal model because of the complex model geometry. This is the main motivation for carrying out the present work.The objective of the present work is to provide an approximate approach to predict the heat extraction and thermal power from a closed-loop geothermal reservoir. This paper is organized as follows: in Section2, the problem description and formulation, including governing equations and boundary and initial conditions, are presented; Section 3 provides the dimensionless formulation, followed by solution method in Section 4. The model validation and numerical results are given in Section 5 and conclusions are presented in Section 6.


## 2. Problem Formulation

### 2.1. Problem Description and Assumptions

The model geometry studied is shown in Figure1. It contains one simplified well containing two vertical parts and one horizontal part. The initial temperature of the system is a function of depth; that is,T0 =A0z +B0, whereA0 is the geothermal gradient andB0 is the surface ground temperature. At time t>0, a fluid (water) with a constant temperature Tin is injected into the system at a constant volumetric rateQ1. The geometrical sizes of the wellbore shown in Figure 1 are defined as follows: the lengths of the injection, horizontal, and production well sections areL1,L2, andL3 (L3 =L1), respectively; the radius of the wellbore is rhi (i=1,2,3), possibly including all layers such as steel casing and cementation if they exist as shown in Figure 1. The origin of the general coordinate system (GCS), denoted by the coordinates (x,y,z), is located at the injection point of the vertical well and coincides with the origin of the local cylindrical coordinate system (LCCS), denoted by the coordinate (r1,z), for the injection well. The origin of the LCCS, denoted by the coordinate (r2,x), for the horizontal well is at the bottom of the injection well, that is, intersection A. In a similar way, the origin of the local cylindrical coordinate system (LCCS), denoted by the coordinate (r3,z1), for the production well is located at the bottom of the production well, that is, (L2, 0,L1) in the GCS.Figure 1
Model geometry. (a) Injection well system and (b) cross section of the well system.
(a)
(b)In order to make the problem tractable, some assumptions are made:(a)
The fluid is single phase, incompressible, and Newtonian, and the rock is impermeable.(b)
The material properties of the fluid and the rock are constants, independent of temperature.(c)
When there are multiple layers around the wellbore, the heat transfer through the layers is characterized by one overall heat transfer coefficient.(d)
The heat diffusion around the wellbore can be regarded to be axis-symmetric.Assumption (d) is the most important one based on which the present model is simplified. According to JiJi [50, page 242–246], for a semi-infinite domain with a homogenous initial temperature and one time-dependent surface heat flux imposed at the boundary at time t>0, the penetration depth of the thermal layer is in the order of 6κrt1/2, where κr=λr/(ρrcr) is the thermal diffusivity with λr, ρr, and cr being the thermal conductivity, mass density, and specific heat capacity, respectively, of the media. Take shale, for example, when λr=1.42 W/(m·K), ρr=2057 Kg/m3, and cr=2151 J/(Kg·K) (pp 106, [8]), the thermal diffusivity κr=3.209×10-7 m2/s. This means that, after 60 days and 30 years, the thermal layer is about 3.16 m and 42.68 m, respectively, illustrating the slow movement of the thermal front in rock. Because the reservoir size is large enough (say in the order of 1 km), the interaction of heat transfer in the horizontal and vertical directions around the bottom of the vertical wells can be neglected.Based on the above assumptions, the current model can be simplified greatly without causing large errors in predicting the temperature profiles.


### 2.2. Governing Equations

The governing equations can be written in the corresponding local cylindrical coordinate system (LCCS) for the vertical and horizontal parts of the well. Besides, temperature continuity is enforced at the intersection point at the heel of the horizontal well.


#### 2.2.1. Heat Exchange along the Wellbore

As the heat flow is uniform along the well, the transient energy equations for the fluid flow are written as follows, according to their LCCSs [28]:(1)ρfcfA1v1∂Tf1∂z+2πrt1U1Tf1-Tb1=-ρfcfA1∂Tf1∂t,(injection),ρlclA3v3∂Tf3∂z1+2πrt3U3Tb3-Tf3=ρfcfA3∂Tf3∂t,(production),for the vertical wells, and (2)ρfcfA2v2∂Tf2∂x+2πrt2U2Tf2-Tb2=-ρfcfA2∂Tf2∂t,for the horizontal well, where ρf and cf are the mass density and specific heat capacity, respectively, of the fluid and rtl, Al, Ul, vl, Tfl, and Tbl (l=1 for injection, l=2 for horizontal, and l=3 for production) denote the radius of fluid flow channel, the areas of the wellbore cross section, overall heat transfer coefficients (OHTCs), fluid velocity in the wellbore, fluid temperature, and the temperature at the wellbore wall, respectively. The value of the radius of fluid flow channel will be determined based on the condition as follows. If fluid flows in a pipe, it is equal to the pipe radius; if the fluid contacts directly with the rock formation, it is equal to the wellbore radius; that is, rtl=rhl.The OHTCsUl can be calculated based on Willhite’s equation [51] which considers a multilayered wellbore structure(3)1Ul=rtortlhl+rtoln⁡rto/rtlktub.+rtoln⁡rins./rtokins+rtorins.hc′+hr′+rtoln⁡rco/rcikcas.+rtoln⁡rh/rcokcem.,where l=1 or 2, for two well parts. The cross-section geometry can be found in Figure 1(b) and the subscripts provide the layers.k denote the thermal conductivities for different layers. hc′ is the natural convection and conduction HTC through the annulus, and hr′ is the radiation HTC through the annulus;hℓ is the HTC between the fluid and tubing. If the effects of the annulus, casing, and cements are not taken into account, the first two terms on the right side of (3) are used to calculateUℓ. If the fluid contacts directly with the rock formation, rto=rtl andUℓ =hℓ. From (3) we also know that the OHTC is mainly determined by the minimum value of the denominators of the terms on the right side. This indicates that if the thermal conductivity of some layer is very small, this layer will work as a thermal insulator.When fluid flows in a tubing, the HTC between the fluid and tubing, that is,hl, is obtained by using the relationship Nul=hlD/kf, where Nu denotes the Nusselt number, kf is the thermal conductivity of the fluid, andD is the hydraulic diameter of the tubing. For fully developed laminar flow in a pipe with circular cross section, the Nusselt number Nu=3.66, while, for transitional and turbulent flows, the Nusselt number is obtained by using the well-known Gnielinski correlation [52](4)Nu=ξ/8Re-1000Prf1+12.7ξ/8Prf2/3-1PrfPrr0.11,ξ=0.79ln⁡Re-1.64-2,when 0.5<Prf<2000 and 3000<Re<5×106. In the present model, the Reynolds and Prandtl numbers related to fluid flowing in the tubing or wellbore and the rock formation are expressed as (5)Rel=ρfvl2rtlμ,Prf=μcfkf,Prrl=μcrlkrl,where ρf and cf denote the mass density and heat capacity, respectively, of the fluid and krl and crl denote the thermal conductivity and heat capacity, respectively, of the rock formation.


#### 2.2.2. Heat Conduction in the Rock Formation Surrounding the Wellbore

As the radial heat conduction dominates the thermal diffusion process around the wellbore, the equations are written as follows in their corresponding LCCS:(6)∂Trl∂t=κrl1rl∂∂rlrl∂Trl∂rl,l=1,2,3,whereTr1,Tr2, andTr3 denote the temperature of the formation around the injection, horizontal, and production wells, respectively. The thermal diffusivities of the rock formation are denoted by κrl=krl/(ρrlcrl), with ρrl being the mass density of the rock formation. It should be noted that the radial spatial variables for the injection, horizontal, and production wells arer1,r2, andr3, respectively.The heat transfer conditions between the rock and the whole well system are given as(7)2πrtlUlTbl-Tfl=2πrhlλrl∂Trl∂rl,atrl=rhl.Based on the above equation, if the OHTCUl=0, ∂Trl/∂rl=0, indicating no heat flowing into or out of the rock formation. From the point of view of heat extraction, a thermal insulation layer along the production well can enhance the heat production by preventing the heat loss along the well during fluid flowing upwards.


### 2.3. Boundary and Initial Conditions

The injection rateQ1 is prescribed and the injection temperature is(8)Tf1=Tin,at0,0,0inGCS.The bottomhole temperature (BHT) of the vertical well is used as the input conditions for the horizontal well and the output temperature of the horizontal well is used as the input conditions for the production well, that is, at the intersectionsA andB(9)Tf1=Tf2,at0,0,HinGCS,Tf2=Tf3,atL,0,HinGCS.The ground temperature at the surface is a constant(10)Tr1=B0,onz=0,and the initial temperature of the whole system is a function of depth(11)T0=A0z+B0.The heat extraction rate or thermal power output by the fluid with output temperatureTout is expressed as(12)W=Q3ρwcwTout-Q1ρwcwTin,where Q1=Q3=Qin and from which the total heat extracted is(13)Φ=∫0tWtdt=∫0tQinρwcwTout-Tindt.


## 2.1. Problem Description and Assumptions

The model geometry studied is shown in Figure1. It contains one simplified well containing two vertical parts and one horizontal part. The initial temperature of the system is a function of depth; that is,T0 =A0z +B0, whereA0 is the geothermal gradient andB0 is the surface ground temperature. At time t>0, a fluid (water) with a constant temperature Tin is injected into the system at a constant volumetric rateQ1. The geometrical sizes of the wellbore shown in Figure 1 are defined as follows: the lengths of the injection, horizontal, and production well sections areL1,L2, andL3 (L3 =L1), respectively; the radius of the wellbore is rhi (i=1,2,3), possibly including all layers such as steel casing and cementation if they exist as shown in Figure 1. The origin of the general coordinate system (GCS), denoted by the coordinates (x,y,z), is located at the injection point of the vertical well and coincides with the origin of the local cylindrical coordinate system (LCCS), denoted by the coordinate (r1,z), for the injection well. The origin of the LCCS, denoted by the coordinate (r2,x), for the horizontal well is at the bottom of the injection well, that is, intersection A. In a similar way, the origin of the local cylindrical coordinate system (LCCS), denoted by the coordinate (r3,z1), for the production well is located at the bottom of the production well, that is, (L2, 0,L1) in the GCS.Figure 1
Model geometry. (a) Injection well system and (b) cross section of the well system.
(a)
(b)In order to make the problem tractable, some assumptions are made:(a)
The fluid is single phase, incompressible, and Newtonian, and the rock is impermeable.(b)
The material properties of the fluid and the rock are constants, independent of temperature.(c)
When there are multiple layers around the wellbore, the heat transfer through the layers is characterized by one overall heat transfer coefficient.(d)
The heat diffusion around the wellbore can be regarded to be axis-symmetric.Assumption (d) is the most important one based on which the present model is simplified. According to JiJi [50, page 242–246], for a semi-infinite domain with a homogenous initial temperature and one time-dependent surface heat flux imposed at the boundary at time t>0, the penetration depth of the thermal layer is in the order of 6κrt1/2, where κr=λr/(ρrcr) is the thermal diffusivity with λr, ρr, and cr being the thermal conductivity, mass density, and specific heat capacity, respectively, of the media. Take shale, for example, when λr=1.42 W/(m·K), ρr=2057 Kg/m3, and cr=2151 J/(Kg·K) (pp 106, [8]), the thermal diffusivity κr=3.209×10-7 m2/s. This means that, after 60 days and 30 years, the thermal layer is about 3.16 m and 42.68 m, respectively, illustrating the slow movement of the thermal front in rock. Because the reservoir size is large enough (say in the order of 1 km), the interaction of heat transfer in the horizontal and vertical directions around the bottom of the vertical wells can be neglected.Based on the above assumptions, the current model can be simplified greatly without causing large errors in predicting the temperature profiles.


## 2.2. Governing Equations

The governing equations can be written in the corresponding local cylindrical coordinate system (LCCS) for the vertical and horizontal parts of the well. Besides, temperature continuity is enforced at the intersection point at the heel of the horizontal well.


### 2.2.1. Heat Exchange along the Wellbore

As the heat flow is uniform along the well, the transient energy equations for the fluid flow are written as follows, according to their LCCSs [28]:(1)ρfcfA1v1∂Tf1∂z+2πrt1U1Tf1-Tb1=-ρfcfA1∂Tf1∂t,(injection),ρlclA3v3∂Tf3∂z1+2πrt3U3Tb3-Tf3=ρfcfA3∂Tf3∂t,(production),for the vertical wells, and (2)ρfcfA2v2∂Tf2∂x+2πrt2U2Tf2-Tb2=-ρfcfA2∂Tf2∂t,for the horizontal well, where ρf and cf are the mass density and specific heat capacity, respectively, of the fluid and rtl, Al, Ul, vl, Tfl, and Tbl (l=1 for injection, l=2 for horizontal, and l=3 for production) denote the radius of fluid flow channel, the areas of the wellbore cross section, overall heat transfer coefficients (OHTCs), fluid velocity in the wellbore, fluid temperature, and the temperature at the wellbore wall, respectively. The value of the radius of fluid flow channel will be determined based on the condition as follows. If fluid flows in a pipe, it is equal to the pipe radius; if the fluid contacts directly with the rock formation, it is equal to the wellbore radius; that is, rtl=rhl.The OHTCsUl can be calculated based on Willhite’s equation [51] which considers a multilayered wellbore structure(3)1Ul=rtortlhl+rtoln⁡rto/rtlktub.+rtoln⁡rins./rtokins+rtorins.hc′+hr′+rtoln⁡rco/rcikcas.+rtoln⁡rh/rcokcem.,where l=1 or 2, for two well parts. The cross-section geometry can be found in Figure 1(b) and the subscripts provide the layers.k denote the thermal conductivities for different layers. hc′ is the natural convection and conduction HTC through the annulus, and hr′ is the radiation HTC through the annulus;hℓ is the HTC between the fluid and tubing. If the effects of the annulus, casing, and cements are not taken into account, the first two terms on the right side of (3) are used to calculateUℓ. If the fluid contacts directly with the rock formation, rto=rtl andUℓ =hℓ. From (3) we also know that the OHTC is mainly determined by the minimum value of the denominators of the terms on the right side. This indicates that if the thermal conductivity of some layer is very small, this layer will work as a thermal insulator.When fluid flows in a tubing, the HTC between the fluid and tubing, that is,hl, is obtained by using the relationship Nul=hlD/kf, where Nu denotes the Nusselt number, kf is the thermal conductivity of the fluid, andD is the hydraulic diameter of the tubing. For fully developed laminar flow in a pipe with circular cross section, the Nusselt number Nu=3.66, while, for transitional and turbulent flows, the Nusselt number is obtained by using the well-known Gnielinski correlation [52](4)Nu=ξ/8Re-1000Prf1+12.7ξ/8Prf2/3-1PrfPrr0.11,ξ=0.79ln⁡Re-1.64-2,when 0.5<Prf<2000 and 3000<Re<5×106. In the present model, the Reynolds and Prandtl numbers related to fluid flowing in the tubing or wellbore and the rock formation are expressed as (5)Rel=ρfvl2rtlμ,Prf=μcfkf,Prrl=μcrlkrl,where ρf and cf denote the mass density and heat capacity, respectively, of the fluid and krl and crl denote the thermal conductivity and heat capacity, respectively, of the rock formation.


### 2.2.2. Heat Conduction in the Rock Formation Surrounding the Wellbore

As the radial heat conduction dominates the thermal diffusion process around the wellbore, the equations are written as follows in their corresponding LCCS:(6)∂Trl∂t=κrl1rl∂∂rlrl∂Trl∂rl,l=1,2,3,whereTr1,Tr2, andTr3 denote the temperature of the formation around the injection, horizontal, and production wells, respectively. The thermal diffusivities of the rock formation are denoted by κrl=krl/(ρrlcrl), with ρrl being the mass density of the rock formation. It should be noted that the radial spatial variables for the injection, horizontal, and production wells arer1,r2, andr3, respectively.The heat transfer conditions between the rock and the whole well system are given as(7)2πrtlUlTbl-Tfl=2πrhlλrl∂Trl∂rl,atrl=rhl.Based on the above equation, if the OHTCUl=0, ∂Trl/∂rl=0, indicating no heat flowing into or out of the rock formation. From the point of view of heat extraction, a thermal insulation layer along the production well can enhance the heat production by preventing the heat loss along the well during fluid flowing upwards.


## 2.2.1. Heat Exchange along the Wellbore

As the heat flow is uniform along the well, the transient energy equations for the fluid flow are written as follows, according to their LCCSs [28]:(1)ρfcfA1v1∂Tf1∂z+2πrt1U1Tf1-Tb1=-ρfcfA1∂Tf1∂t,(injection),ρlclA3v3∂Tf3∂z1+2πrt3U3Tb3-Tf3=ρfcfA3∂Tf3∂t,(production),for the vertical wells, and (2)ρfcfA2v2∂Tf2∂x+2πrt2U2Tf2-Tb2=-ρfcfA2∂Tf2∂t,for the horizontal well, where ρf and cf are the mass density and specific heat capacity, respectively, of the fluid and rtl, Al, Ul, vl, Tfl, and Tbl (l=1 for injection, l=2 for horizontal, and l=3 for production) denote the radius of fluid flow channel, the areas of the wellbore cross section, overall heat transfer coefficients (OHTCs), fluid velocity in the wellbore, fluid temperature, and the temperature at the wellbore wall, respectively. The value of the radius of fluid flow channel will be determined based on the condition as follows. If fluid flows in a pipe, it is equal to the pipe radius; if the fluid contacts directly with the rock formation, it is equal to the wellbore radius; that is, rtl=rhl.The OHTCsUl can be calculated based on Willhite’s equation [51] which considers a multilayered wellbore structure(3)1Ul=rtortlhl+rtoln⁡rto/rtlktub.+rtoln⁡rins./rtokins+rtorins.hc′+hr′+rtoln⁡rco/rcikcas.+rtoln⁡rh/rcokcem.,where l=1 or 2, for two well parts. The cross-section geometry can be found in Figure 1(b) and the subscripts provide the layers.k denote the thermal conductivities for different layers. hc′ is the natural convection and conduction HTC through the annulus, and hr′ is the radiation HTC through the annulus;hℓ is the HTC between the fluid and tubing. If the effects of the annulus, casing, and cements are not taken into account, the first two terms on the right side of (3) are used to calculateUℓ. If the fluid contacts directly with the rock formation, rto=rtl andUℓ =hℓ. From (3) we also know that the OHTC is mainly determined by the minimum value of the denominators of the terms on the right side. This indicates that if the thermal conductivity of some layer is very small, this layer will work as a thermal insulator.When fluid flows in a tubing, the HTC between the fluid and tubing, that is,hl, is obtained by using the relationship Nul=hlD/kf, where Nu denotes the Nusselt number, kf is the thermal conductivity of the fluid, andD is the hydraulic diameter of the tubing. For fully developed laminar flow in a pipe with circular cross section, the Nusselt number Nu=3.66, while, for transitional and turbulent flows, the Nusselt number is obtained by using the well-known Gnielinski correlation [52](4)Nu=ξ/8Re-1000Prf1+12.7ξ/8Prf2/3-1PrfPrr0.11,ξ=0.79ln⁡Re-1.64-2,when 0.5<Prf<2000 and 3000<Re<5×106. In the present model, the Reynolds and Prandtl numbers related to fluid flowing in the tubing or wellbore and the rock formation are expressed as (5)Rel=ρfvl2rtlμ,Prf=μcfkf,Prrl=μcrlkrl,where ρf and cf denote the mass density and heat capacity, respectively, of the fluid and krl and crl denote the thermal conductivity and heat capacity, respectively, of the rock formation.


## 2.2.2. Heat Conduction in the Rock Formation Surrounding the Wellbore

As the radial heat conduction dominates the thermal diffusion process around the wellbore, the equations are written as follows in their corresponding LCCS:(6)∂Trl∂t=κrl1rl∂∂rlrl∂Trl∂rl,l=1,2,3,whereTr1,Tr2, andTr3 denote the temperature of the formation around the injection, horizontal, and production wells, respectively. The thermal diffusivities of the rock formation are denoted by κrl=krl/(ρrlcrl), with ρrl being the mass density of the rock formation. It should be noted that the radial spatial variables for the injection, horizontal, and production wells arer1,r2, andr3, respectively.The heat transfer conditions between the rock and the whole well system are given as(7)2πrtlUlTbl-Tfl=2πrhlλrl∂Trl∂rl,atrl=rhl.Based on the above equation, if the OHTCUl=0, ∂Trl/∂rl=0, indicating no heat flowing into or out of the rock formation. From the point of view of heat extraction, a thermal insulation layer along the production well can enhance the heat production by preventing the heat loss along the well during fluid flowing upwards.


## 2.3. Boundary and Initial Conditions

The injection rateQ1 is prescribed and the injection temperature is(8)Tf1=Tin,at0,0,0inGCS.The bottomhole temperature (BHT) of the vertical well is used as the input conditions for the horizontal well and the output temperature of the horizontal well is used as the input conditions for the production well, that is, at the intersectionsA andB(9)Tf1=Tf2,at0,0,HinGCS,Tf2=Tf3,atL,0,HinGCS.The ground temperature at the surface is a constant(10)Tr1=B0,onz=0,and the initial temperature of the whole system is a function of depth(11)T0=A0z+B0.The heat extraction rate or thermal power output by the fluid with output temperatureTout is expressed as(12)W=Q3ρwcwTout-Q1ρwcwTin,where Q1=Q3=Qin and from which the total heat extracted is(13)Φ=∫0tWtdt=∫0tQinρwcwTout-Tindt.


## 3. Dimensionless Formulation

The governing equations, boundary and initial conditions, are simplified with the following transformation:(14)Z=zLz,X=xLx,Rl=rlrhl,τ=κr1trh12,Θfl=TflT∗,Θbl=TblT∗,Θfl=TflT∗,εl=κrlrh12κr1rhl2,Bil=rtlUlλrl,χl=Llκrlrh2vl,αl=2πLlrtlUlρfcfQl,γ1=A0HT∗,Θr10=Θf10=γ1Z-1+1,T∗=A0H+B0,Lz=L1,Lx=L2,where Θfl, Θbl, and Θrl (l=1 for injection, l=2 for horizontal, and l=3 for production) denote dimensionless temperature in the fluid, at the wall and in the formation, respectively.By using the above variables, the simplified formulation for the three subproblems of the original model is listed as follows.(a) Injection Well. The governing and heat balance equations for the injection well are(15)∂Θr1∂τ=1R1∂∂R1R1∂Θr1∂R1,Bi1Θb1-Θf1=∂Θr1∂R1atR1=1,χ1∂Θf1∂τ+∂Θf1∂Z=α1Θb1-Θf1,with the initial and boundary conditions(16)Θr1=Θf1=γ1Z-1+1,onτ=0,Θf1=Θin,at1,0inLCCS.(b) Horizontal Well. The governing and heat balance equations for the horizontal well are(17)∂Θr2∂τ=ε2R2∂∂R2R2∂Θr2∂R2,Bi2Θb2-Θf2=∂Θr2∂R2atR2=1,χ2∂Θf2∂τ+∂Θf2∂X=α2Θb2-Θf2,with the initial and boundary conditions becoming(18)Θf2=Θf1,at1,0inLCCS,Θr2=Θf2=1,onτ=0.(c) Production Well. The governing and heat balance equations for the production well are(19)∂Θr3∂τ=ε3R3∂∂R3R3∂Θr3∂R3,Bi3Θb3-Θf3=∂Θr3∂R3atR3=1,χ3∂Θf3∂τ+∂Θf3∂Z=α3Θb3-Θf3,with the initial and boundary conditions(20)Θr3=Θf3=γ3Z-+1,onτ=0,Θf3=Θf2,at0,0inLCCS,where Z-=1-Z denotes the distance of the point from the bottom of production well, γ3=-γ1. The subproblem for production well is similar to that for injection well as it can be regarded as a model with the origin of the CS at the bottom of production well and with a negative geothermal gradient γ3=-γ1.


## 4. Solution Method

By using the Laplace transformation, the analytical solutions for the three stages are obtained. In the following equations, the symbol∧ denotes the variables which are Laplace transformed and s is a complex number as the Laplace symbol. For example, f^s denotes the Laplace transform of the function ft.


### 4.1. Injection Well

From the Laplace transform of the governing equation, that is, the first of (15), for the heat diffusion in the rock around the injection well, we obtain the rock temperature(21)Θ^r1=F1Z,sK0sR1+F2Z,sI0sR1+γ1Z-1+1s,where Θ^r1 is the Laplace transform of the dimensionless temperature Θr1; In and Kn are the modified Bessel functions of the first and second kind of ordern, respectively; andF1(s) andF2(s) are unknowns to be determined by the boundary conditions. As Θr1 is finite when R1→+∞, it is easy to know that F2Z=0. Therefore, the temperature at the wellbore wall is obtained(22)Θ^b1=Θ^r1R1=1=F1Z,sK0s+γZ-1+1s.By substituting (21) and (22) into the second of (15) we have (23)Θ^f1=F1Z,sΔ+γZ-1+1s,Δ=sK1sBi1+K0s,which, after being used in the third of (15), produces an ordinary differential equation with respect to the unknown F1(Z,s); that is,(24)dF1dz+β1Δ1F1+γ1sΔ1=0,whereβ1=χ1sΔ1+α1Δ1-K0s.The solution of the above equation is found to be(25)F1Z,s=-γ1sβ+e-βZ/ΔC1s,where the unknown functionC1(s) can be determined by using (23), (25), and the injection boundary condition (16)(26)C1=Θ^insβ+βγ1-β+Δγ1sβΔ=Θ^ins+γ1-1sΔ+γ1sβ.Therefore, the solutions for the temperatures of the rock,Θ^r1, and fluid, Θ^f1, are obtained and rewritten as(27)Θ^r1=Θ^ins+γ1-1sA1+γ1sA2e-Z/A3-γ1sA2+γ1Z-1+1s,Θ^f1=Θ^ins+γ1-1s+γ1sA3e-Z/A3-γ1sA3+γ1Z-1+1s,where the functionsA1,A2, andA3 are defined as(28)A1=Bi1K0sR1/K1ss+Bi1K0s/K1s,A2=Bi1K0sR1/K1sχ1ss+Bi1K0s/K1s+α1s,A3=s+Bi1K0s/K1sχ1ss+Bi1K0s/K1s+α1s,where Θ^in is the Laplace transform of the injection temperature.The Laplace transformation of the bottomhole temperature (BHT) of the injection well is(29)BHT1=Gs=Θ^ins+γ1-1s+γ1sA3e-1/A3-γ1sA3+1s,which will be used as the input condition for the horizontal well.


### 4.2. Horizontal Well

Initially, the temperature along horizontal well is identical due to the same depth. Based on the above calculations, the thermal front only penetrates less than 50 meters into the rock within 30 years. In the region between the upper and lower thermal front, the largest temperature is less than 4°C. Therefore, it is reasonable to assume a constant temperature in the rock formation around the horizontal well.In a similar way, the solutions for the horizontal well are obtained(30)Θ^r2=sGs-1se-X/B3B1+1s,Θ^f2=sGs-1se-X/B3+1s,where the functionsB1 andB3 are defined as(31)B1=Bi2K0s/ε2R2/K1s/ε2Bi2K0s/ε2/K1s/ε2+s/ε2,B3=s/ε2+Bi2K0s/ε2/K1s/ε2χ2ss/ε2+Bi2K0s/ε2/K1s/ε2+α2s/ε2.The Laplace transformation of the temperature at the end of the horizontal well, also BHT of the production well, is obtained by usingX=1; that is,(32)BHT2=Hs=sGs-1se-1/B3+1s,which will be used as the input condition for the production well.


### 4.3. Production Well

The temperatures for the fluid in the production well and neighboring rock formation can be obtained in the same way as those for injection and horizontal well and the calculation details are omitted here(33)Θ^r3=Hss-1sC1+γ3sC2e-Z-/C3-γ3sC2+γ3Z-+1s,Θ^f3=Hss-1s+γ3sC3e-Z-/C3-γ3sC3+γ3Z-+1s,where the functionsC1,C2, andC3 are defined as(34)C1=Bi3K0s/ε3R1/K1s/ε3Bi3K0s/ε3/K1s/ε3+s/ε3,C2=Bi3K0s/ε3R1/K1s/ε3χ3sBi3K0s/ε3/K1s/ε3+s/ε3+α3s/ε3,C3=s/ε3+Bi3K0s/ε3/K1s/ε3χ3ss/ε3+Bi3K0s/ε3/K1s/ε3+α3s/ε3.The Laplace transformation of the temperature at the end of the horizontal well, also BHT of the production well, is obtained by usingX=1; that is,(35)BHT2=Hs=sGs-1se-1/B3+1s,which will be used as the input condition for the production well.


## 4.1. Injection Well

From the Laplace transform of the governing equation, that is, the first of (15), for the heat diffusion in the rock around the injection well, we obtain the rock temperature(21)Θ^r1=F1Z,sK0sR1+F2Z,sI0sR1+γ1Z-1+1s,where Θ^r1 is the Laplace transform of the dimensionless temperature Θr1; In and Kn are the modified Bessel functions of the first and second kind of ordern, respectively; andF1(s) andF2(s) are unknowns to be determined by the boundary conditions. As Θr1 is finite when R1→+∞, it is easy to know that F2Z=0. Therefore, the temperature at the wellbore wall is obtained(22)Θ^b1=Θ^r1R1=1=F1Z,sK0s+γZ-1+1s.By substituting (21) and (22) into the second of (15) we have (23)Θ^f1=F1Z,sΔ+γZ-1+1s,Δ=sK1sBi1+K0s,which, after being used in the third of (15), produces an ordinary differential equation with respect to the unknown F1(Z,s); that is,(24)dF1dz+β1Δ1F1+γ1sΔ1=0,whereβ1=χ1sΔ1+α1Δ1-K0s.The solution of the above equation is found to be(25)F1Z,s=-γ1sβ+e-βZ/ΔC1s,where the unknown functionC1(s) can be determined by using (23), (25), and the injection boundary condition (16)(26)C1=Θ^insβ+βγ1-β+Δγ1sβΔ=Θ^ins+γ1-1sΔ+γ1sβ.Therefore, the solutions for the temperatures of the rock,Θ^r1, and fluid, Θ^f1, are obtained and rewritten as(27)Θ^r1=Θ^ins+γ1-1sA1+γ1sA2e-Z/A3-γ1sA2+γ1Z-1+1s,Θ^f1=Θ^ins+γ1-1s+γ1sA3e-Z/A3-γ1sA3+γ1Z-1+1s,where the functionsA1,A2, andA3 are defined as(28)A1=Bi1K0sR1/K1ss+Bi1K0s/K1s,A2=Bi1K0sR1/K1sχ1ss+Bi1K0s/K1s+α1s,A3=s+Bi1K0s/K1sχ1ss+Bi1K0s/K1s+α1s,where Θ^in is the Laplace transform of the injection temperature.The Laplace transformation of the bottomhole temperature (BHT) of the injection well is(29)BHT1=Gs=Θ^ins+γ1-1s+γ1sA3e-1/A3-γ1sA3+1s,which will be used as the input condition for the horizontal well.


## 4.2. Horizontal Well

Initially, the temperature along horizontal well is identical due to the same depth. Based on the above calculations, the thermal front only penetrates less than 50 meters into the rock within 30 years. In the region between the upper and lower thermal front, the largest temperature is less than 4°C. Therefore, it is reasonable to assume a constant temperature in the rock formation around the horizontal well.In a similar way, the solutions for the horizontal well are obtained(30)Θ^r2=sGs-1se-X/B3B1+1s,Θ^f2=sGs-1se-X/B3+1s,where the functionsB1 andB3 are defined as(31)B1=Bi2K0s/ε2R2/K1s/ε2Bi2K0s/ε2/K1s/ε2+s/ε2,B3=s/ε2+Bi2K0s/ε2/K1s/ε2χ2ss/ε2+Bi2K0s/ε2/K1s/ε2+α2s/ε2.The Laplace transformation of the temperature at the end of the horizontal well, also BHT of the production well, is obtained by usingX=1; that is,(32)BHT2=Hs=sGs-1se-1/B3+1s,which will be used as the input condition for the production well.


## 4.3. Production Well

The temperatures for the fluid in the production well and neighboring rock formation can be obtained in the same way as those for injection and horizontal well and the calculation details are omitted here(33)Θ^r3=Hss-1sC1+γ3sC2e-Z-/C3-γ3sC2+γ3Z-+1s,Θ^f3=Hss-1s+γ3sC3e-Z-/C3-γ3sC3+γ3Z-+1s,where the functionsC1,C2, andC3 are defined as(34)C1=Bi3K0s/ε3R1/K1s/ε3Bi3K0s/ε3/K1s/ε3+s/ε3,C2=Bi3K0s/ε3R1/K1s/ε3χ3sBi3K0s/ε3/K1s/ε3+s/ε3+α3s/ε3,C3=s/ε3+Bi3K0s/ε3/K1s/ε3χ3ss/ε3+Bi3K0s/ε3/K1s/ε3+α3s/ε3.The Laplace transformation of the temperature at the end of the horizontal well, also BHT of the production well, is obtained by usingX=1; that is,(35)BHT2=Hs=sGs-1se-1/B3+1s,which will be used as the input condition for the production well.


## 5. Validation and Numerical Results

### 5.1. Method Validation

In the present model, in order to obtain analytical solutions for temperature prediction in a closed-loop geothermal reservoir, the wellbore for fluid flow and heat transfer is divided into three parts with the output of one well part as the input conditions of the next well part. In order to show that these assumptions are reasonable and to show the error resulting from decoupling of the whole process, the results from the present model are compared with those from numerical simulator TOUGH 2. An automatic Laplace inversion technique developed by D’Amore et al. [53] based on Fourier series is used in the present model to obtain the values in the time space.Figure2 shows spatial discretization of the rock formation for the numerical calculation by using TOUGH2. As the temperature of the fluid in the radial direction in the pipe is assumed to be the same, it is denoted by blue node. The node for the rock is located in the middle of the each element. Figure 3 compares the temperature variations at pointsA,B, andC, which corresponds to the output points of the injection, horizontal, and production well parts, respectively, between TOUGH2 [37] and the present model. The solid curves denote results from TOUGH2 and the dashed ones are results from the present model. Figures 3(a) and 3(b) are for temperatures at short (20 days) and long (2.74 years) times of circulation, respectively, for the case Q=0.05 m3/s, while Figures 3(c) and 3(d) are for small Q=0.01 m3/s and large Q=0.20 m3/s, respectively. Other parameters used here are L1=2000 m, L2=3000 m,B0 = 20°C, and Tin=30°C.Figure 2
Grid for numerical calculation by using TOUGH2. The blue node is for the fluid and red node for the rock formation.Figure 3
Comparison of the temperatures at pointsA, B, andC. (a) For small time and (b) for large time when Q=0.05 m3/s; (c) for small Q=0.01 m3/s and (d) for large Q=0.2 m3/s.
(a)
(b)
(c)
(d)It can be seen that the temperature results at pointsA,B, andC predicted from both methods do not match at short time period, say 10 days, as shown in Figure 3(a) for the case when Q=0.05 m3/s. The temperature difference is mainly caused by the neglect of thermal interaction at the heel of the injection and production wells (i.e., around pointsA andB) in the present model. However, for large time of circulation, the difference in the temperature results from both methods is very small, generally less than 4°C, as shown in Figure 3(b). The temperature differences after 365 days and 1000 days between the present model and TOUGH2 are −3.87°C and −3.65°C, respectively. In order to further show the accuracy of the present model, the temperatures at these three points in another two cases, that is, Q=0.01 m3/s and Q=0.20 m3/s, are also compared with those obtained from TOUGH2 and show good agreement. This means that the present model can be used to predict approximately the output temperature of the closed-loop system without causing large errors.


### 5.2. Numerical Results

In this section, the case with fluid fully contacted (without any tubing, casing, and cement) with the wellbore is first investigated. Under this condition, the radius of the tubing for fluid flow,rti, is equal to the wellbore radius, rhi. By doing so, two objectives are achieved. First, the important factors affecting the fluid and heat flow behaviour are identified; second, the critical location determining whether the heat is flowing into or out of the rock formation is identified for thermal isolation design, especially for production well. Water is chosen for the fluid; sandstone and shale are chosen for the rock around the vertical and horizontal wells, respectively, for the following calculations. The parameters used for the examples are listed in Table 1 unless otherwise specified.Table 1
Physical parameters for the present calculation (the parameter values for thermal properties of shale are from Eppelbaum et al. [8] and the value of geothermal gradient is from Quick et al. [9]).
Parameter
Value
Wellbore radius,rh (m)
0.1
Length of vertical well,H (m)
2000
Length of horizontal well (m)
3000
Geothermal gradient,A0 (K/m)
0.047
Surface temperature,B0 (°C)
20
Fluid density (Kg/(m3))
900
Injection rate,Q (m3/s)
0.05
Injection temperature (°C)
30.0
Fluid specific heat (J/(Kg·K))
4200
Fluid thermal conductivity (W/(m·K))
0.68
Fluid viscosity (Pa·s)
0.0004
Shale thermal conduct (W/(m·K))
1.42
Shale specific heat (J/(Kg·K))
2151
Shale density (Kg/(m3))
2057


#### 5.2.1. Temperature along the Vertical Wells

Figure4 displays the dimensionless temperature changes, which are defined by physical temperature changes with respect to the initial state divided by (A0H+B0), along the injection and production wells, respectively. The red curves refer to the initial state. For the injection well, we can see from Figure 4(a) that (1) before the cold fluid reaches the bottom of the wellbore, the temperature change of the fluid on the upper part of the wellbore decreases linearly with depth, while that on the lower part of the wellbore is almost a constant; (2) the temperature along the injection well decreases very quickly, with a temperature of 46.3°C at the bottom after 3 hours of circulation; and (3) as the injection temperature is lower than the initial ground temperature for most of the well length, most of the fluid flowing downwards absorbs heat from the surrounding rock except a small portion of fluid on the upper part of the wellbore.Figure 4
Temperature change profiles along the (a) injection and (b) production wells whenQ=0.05 m3/s,L1 = 3500 m, andL2 = 6000 m.
(a)
(b)Before the cold fluid reaches the production well, the temperature along the whole production well increases because of hot fluid flowing upwards, as shown by solid and dashed curves in Figure4(b). The production temperature increases from 57°C to 152°C when the time varies from 10 mins to 1.0 hours and then decreases to 126°C after 3 hours. The small difference in temperature change between times t=10 days and t=30 years means that the heat transfer in the fluid approaches pseudo-steady state after 10 days. Moreover, the critical position which indicates if the fluid absorbs heat from or releases heat to the rock formation is around Zc=0.2 or zc=700 m.


#### 5.2.2. Effect of Flow Rates on Output Performance

Figure5 presents the output temperature for different injection rates ranging from Q=0.01 m3/s to 0.14 m3/s. It can be seen that the injection rate plays a major influence on the output temperature. First, for each injection rate, the output temperature increases quickly to a maximum value, as shown in Figure 5(b), and then decreases steadily to a constant value, as shown in Figure 5(a). Second, with increasing the injection rate, the production temperature is decreased. This is due to shortened time as a result of larger fluid speed for heat exchange between fluid and the surrounding rock. In addition, the temperature difference after 30 years between the cases Q=0.01 m3/s and Q=0.14 m3/s is around 40°C. However, when the injection rate is increased to some value (0.05 m3/s for the present case), there is no large difference in the final production temperature, as shown in Figure 5(a).Figure 5
Variation with time of output temperature of the production well under different injection rates whenL1 = 3500 m andL2 = 6000 m: (a) for general view and (b) for local view.
(a)
(b)The thermal power and total heat extracted for the above cases are plotted in Figure6. According to the definition of thermal power in (12), thermal power is linearly proportional to the injection rate and temperature. However, as lower production temperature is obtained when the injection rate is increased, the curves for thermal power shown in Figure 6(b) do not exhibit similar trends to the production temperature in Figure 5(b). It has to be mentioned that the long-time thermal power ranges from 1.71 MWs to 2.21 MWs when Q is changed from 0.01 m3/s to 0.14 m3/s.Figure 6
Thermal power and total heat extracted for different injection rates: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)Figure6(c) displays the total heat extracted from the reservoir up to 30 years under different injection rates. The total heat extracted is expected to increase when the injection rate is increased. However, the total heat extracted after 30 years when Q=0.14 m3/s is only 32.75% larger than that when Q=0.01 m3/s. This means that the cost related to higher pressure required for larger flow rates has to be balanced with the total energy extracted.


#### 5.2.3. Effect of Injection Temperatures on Output Performance

The effect of injection temperatures on output performance is displayed in Figures7 and 8. As for output temperature, it increases evenly with increasing the injection temperature, as shown in Figure 7(a), and the maximum value of output temperature is not sensitive to the injection temperature, as shown in Figure 7(b).Figure 7
Variation with time of output temperature of the production well under different injection temperatures: (a) for large time period and (b) for short time period.
(a)
(b)Figure 8
Thermal power and total heat extracted for different injection temperatures: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)Figure8 presents the thermal power and total heat extracted up to 30 years. As the thermal power is linearly proportional to the temperature and the injection rate is kept constant, the thermal power exhibits a similar trend to the output temperature, as shown in Figures 8 and 7. In addition, it is interesting to find that although higher output temperature is obtained under larger injection temperature, the total heat extracted from the reservoir is reduced, as shown in Figure 8(c).


#### 5.2.4. Effect of Fluid Viscosity on Output Performance

Figures9 and 10 display the final output performance for different fluid viscosities. The water viscosity ranges from about 2.8 × 10−4 Pa·s to 1.3 × 10−3 Pa·s when the temperature changes from 100°C to 10°C, which can be found in the website link (https://en.wikipedia.org/wiki/viscosity). Here μ = 10−4, 5 × 10−4, 2.5 × 10−3, and 1.2 × 10−2 Pa·s are chosen to show the effect of the viscosity on the output temperature. From Figure 9 the output temperature follows almost the same curve under these four cases and the decreasing trend for rapid temperature drop in less than 5 hrs is insensitive to the fluid viscosity.Figure 9
Variation with time of output temperature of the production well under different viscositiesμ: (a) for large time period and (b) for short time period.
(a)
(b)Figure 10
Thermal power and total heat extracted for different fluid viscosities: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)


#### 5.2.5. Effect of Geothermal Gradient on the Output Performance

The output temperatures under different gradient geothermal gradients,A0, are displayed in Figure 11 when other parameters are kept the same. It can be found that (1) as the depth of the vertical well is the same, larger geothermal gradient produces larger maximum output temperature (MOT), which is increased by around 42.2°C if A0 has an increase of 1.5°C. For example, the MOTs for the cases withA0 = 0.020, 0.035, 0.050, 0.065, 0.080, and 0.095 K/m are 76.2, 118.4, 160.5, 202.7, 244.9, and 287.0°C, respectively. However, the output temperature after a long time (say 30 years) is found to be increased by around 4°C whenA0 has an increase of 1.5°C, as shown in Figure 11(b), where the output temperature is 34.2, 38.2, 42.1, 46.0, 49.9, and 53.8°C, respectively, for the above six cases. This means that, compared to the injection conditions, the geothermal gradient plays much less important role in the output performance.Figure 11
Variation with time of output temperature of the production well under different geothermal gradientA0: (a) for large time period and (b) for short time period.
(a)
(b)


#### 5.2.6. Fully Insulated Case

Figure12 compares the production temperatures under fully contacted and fully insulated conditions. The parameters used for the fully insulated case are the same as those in Figure 4 except that the overall heat transfer coefficient along the production well is set to zero.Figure 12
Production temperature under fully contacted and fully insulated conditions along the production well whenQ=0.05 m3/s, L1=3500 m, and L2=6000 m.It is found from Figure12 that there is no large difference in production temperature between these two cases and the production temperature predicted from fully insulated case is a little smaller than that predicted from fully contacted case. From the curves for t=10 days and t=30 years in Figure 4(b) without any casing and cementing, it can be found that (1) the final output temperature decreases very quickly and approaches some pseudo-steady value which is dominated by the inject conditions and (2) the fluid along a larger part (around 80 percent) of wellbore length absorbs heat from the surrounding rock formation and the fluid along the rest of the wellbore releases heat into the formation. Compared with this openhole case, if an insulated layer is put along the whole wellbore, the final output temperature will become smaller although the heat loss into the formation, which occurs along a small part of wellbore length in the openhole case, is prevented in the insulated case, as confirmed in Figure 12.


## 5.1. Method Validation

In the present model, in order to obtain analytical solutions for temperature prediction in a closed-loop geothermal reservoir, the wellbore for fluid flow and heat transfer is divided into three parts with the output of one well part as the input conditions of the next well part. In order to show that these assumptions are reasonable and to show the error resulting from decoupling of the whole process, the results from the present model are compared with those from numerical simulator TOUGH 2. An automatic Laplace inversion technique developed by D’Amore et al. [53] based on Fourier series is used in the present model to obtain the values in the time space.Figure2 shows spatial discretization of the rock formation for the numerical calculation by using TOUGH2. As the temperature of the fluid in the radial direction in the pipe is assumed to be the same, it is denoted by blue node. The node for the rock is located in the middle of the each element. Figure 3 compares the temperature variations at pointsA,B, andC, which corresponds to the output points of the injection, horizontal, and production well parts, respectively, between TOUGH2 [37] and the present model. The solid curves denote results from TOUGH2 and the dashed ones are results from the present model. Figures 3(a) and 3(b) are for temperatures at short (20 days) and long (2.74 years) times of circulation, respectively, for the case Q=0.05 m3/s, while Figures 3(c) and 3(d) are for small Q=0.01 m3/s and large Q=0.20 m3/s, respectively. Other parameters used here are L1=2000 m, L2=3000 m,B0 = 20°C, and Tin=30°C.Figure 2
Grid for numerical calculation by using TOUGH2. The blue node is for the fluid and red node for the rock formation.Figure 3
Comparison of the temperatures at pointsA, B, andC. (a) For small time and (b) for large time when Q=0.05 m3/s; (c) for small Q=0.01 m3/s and (d) for large Q=0.2 m3/s.
(a)
(b)
(c)
(d)It can be seen that the temperature results at pointsA,B, andC predicted from both methods do not match at short time period, say 10 days, as shown in Figure 3(a) for the case when Q=0.05 m3/s. The temperature difference is mainly caused by the neglect of thermal interaction at the heel of the injection and production wells (i.e., around pointsA andB) in the present model. However, for large time of circulation, the difference in the temperature results from both methods is very small, generally less than 4°C, as shown in Figure 3(b). The temperature differences after 365 days and 1000 days between the present model and TOUGH2 are −3.87°C and −3.65°C, respectively. In order to further show the accuracy of the present model, the temperatures at these three points in another two cases, that is, Q=0.01 m3/s and Q=0.20 m3/s, are also compared with those obtained from TOUGH2 and show good agreement. This means that the present model can be used to predict approximately the output temperature of the closed-loop system without causing large errors.


## 5.2. Numerical Results

In this section, the case with fluid fully contacted (without any tubing, casing, and cement) with the wellbore is first investigated. Under this condition, the radius of the tubing for fluid flow,rti, is equal to the wellbore radius, rhi. By doing so, two objectives are achieved. First, the important factors affecting the fluid and heat flow behaviour are identified; second, the critical location determining whether the heat is flowing into or out of the rock formation is identified for thermal isolation design, especially for production well. Water is chosen for the fluid; sandstone and shale are chosen for the rock around the vertical and horizontal wells, respectively, for the following calculations. The parameters used for the examples are listed in Table 1 unless otherwise specified.Table 1
Physical parameters for the present calculation (the parameter values for thermal properties of shale are from Eppelbaum et al. [8] and the value of geothermal gradient is from Quick et al. [9]).
Parameter
Value
Wellbore radius,rh (m)
0.1
Length of vertical well,H (m)
2000
Length of horizontal well (m)
3000
Geothermal gradient,A0 (K/m)
0.047
Surface temperature,B0 (°C)
20
Fluid density (Kg/(m3))
900
Injection rate,Q (m3/s)
0.05
Injection temperature (°C)
30.0
Fluid specific heat (J/(Kg·K))
4200
Fluid thermal conductivity (W/(m·K))
0.68
Fluid viscosity (Pa·s)
0.0004
Shale thermal conduct (W/(m·K))
1.42
Shale specific heat (J/(Kg·K))
2151
Shale density (Kg/(m3))
2057


### 5.2.1. Temperature along the Vertical Wells

Figure4 displays the dimensionless temperature changes, which are defined by physical temperature changes with respect to the initial state divided by (A0H+B0), along the injection and production wells, respectively. The red curves refer to the initial state. For the injection well, we can see from Figure 4(a) that (1) before the cold fluid reaches the bottom of the wellbore, the temperature change of the fluid on the upper part of the wellbore decreases linearly with depth, while that on the lower part of the wellbore is almost a constant; (2) the temperature along the injection well decreases very quickly, with a temperature of 46.3°C at the bottom after 3 hours of circulation; and (3) as the injection temperature is lower than the initial ground temperature for most of the well length, most of the fluid flowing downwards absorbs heat from the surrounding rock except a small portion of fluid on the upper part of the wellbore.Figure 4
Temperature change profiles along the (a) injection and (b) production wells whenQ=0.05 m3/s,L1 = 3500 m, andL2 = 6000 m.
(a)
(b)Before the cold fluid reaches the production well, the temperature along the whole production well increases because of hot fluid flowing upwards, as shown by solid and dashed curves in Figure4(b). The production temperature increases from 57°C to 152°C when the time varies from 10 mins to 1.0 hours and then decreases to 126°C after 3 hours. The small difference in temperature change between times t=10 days and t=30 years means that the heat transfer in the fluid approaches pseudo-steady state after 10 days. Moreover, the critical position which indicates if the fluid absorbs heat from or releases heat to the rock formation is around Zc=0.2 or zc=700 m.


### 5.2.2. Effect of Flow Rates on Output Performance

Figure5 presents the output temperature for different injection rates ranging from Q=0.01 m3/s to 0.14 m3/s. It can be seen that the injection rate plays a major influence on the output temperature. First, for each injection rate, the output temperature increases quickly to a maximum value, as shown in Figure 5(b), and then decreases steadily to a constant value, as shown in Figure 5(a). Second, with increasing the injection rate, the production temperature is decreased. This is due to shortened time as a result of larger fluid speed for heat exchange between fluid and the surrounding rock. In addition, the temperature difference after 30 years between the cases Q=0.01 m3/s and Q=0.14 m3/s is around 40°C. However, when the injection rate is increased to some value (0.05 m3/s for the present case), there is no large difference in the final production temperature, as shown in Figure 5(a).Figure 5
Variation with time of output temperature of the production well under different injection rates whenL1 = 3500 m andL2 = 6000 m: (a) for general view and (b) for local view.
(a)
(b)The thermal power and total heat extracted for the above cases are plotted in Figure6. According to the definition of thermal power in (12), thermal power is linearly proportional to the injection rate and temperature. However, as lower production temperature is obtained when the injection rate is increased, the curves for thermal power shown in Figure 6(b) do not exhibit similar trends to the production temperature in Figure 5(b). It has to be mentioned that the long-time thermal power ranges from 1.71 MWs to 2.21 MWs when Q is changed from 0.01 m3/s to 0.14 m3/s.Figure 6
Thermal power and total heat extracted for different injection rates: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)Figure6(c) displays the total heat extracted from the reservoir up to 30 years under different injection rates. The total heat extracted is expected to increase when the injection rate is increased. However, the total heat extracted after 30 years when Q=0.14 m3/s is only 32.75% larger than that when Q=0.01 m3/s. This means that the cost related to higher pressure required for larger flow rates has to be balanced with the total energy extracted.


### 5.2.3. Effect of Injection Temperatures on Output Performance

The effect of injection temperatures on output performance is displayed in Figures7 and 8. As for output temperature, it increases evenly with increasing the injection temperature, as shown in Figure 7(a), and the maximum value of output temperature is not sensitive to the injection temperature, as shown in Figure 7(b).Figure 7
Variation with time of output temperature of the production well under different injection temperatures: (a) for large time period and (b) for short time period.
(a)
(b)Figure 8
Thermal power and total heat extracted for different injection temperatures: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)Figure8 presents the thermal power and total heat extracted up to 30 years. As the thermal power is linearly proportional to the temperature and the injection rate is kept constant, the thermal power exhibits a similar trend to the output temperature, as shown in Figures 8 and 7. In addition, it is interesting to find that although higher output temperature is obtained under larger injection temperature, the total heat extracted from the reservoir is reduced, as shown in Figure 8(c).


### 5.2.4. Effect of Fluid Viscosity on Output Performance

Figures9 and 10 display the final output performance for different fluid viscosities. The water viscosity ranges from about 2.8 × 10−4 Pa·s to 1.3 × 10−3 Pa·s when the temperature changes from 100°C to 10°C, which can be found in the website link (https://en.wikipedia.org/wiki/viscosity). Here μ = 10−4, 5 × 10−4, 2.5 × 10−3, and 1.2 × 10−2 Pa·s are chosen to show the effect of the viscosity on the output temperature. From Figure 9 the output temperature follows almost the same curve under these four cases and the decreasing trend for rapid temperature drop in less than 5 hrs is insensitive to the fluid viscosity.Figure 9
Variation with time of output temperature of the production well under different viscositiesμ: (a) for large time period and (b) for short time period.
(a)
(b)Figure 10
Thermal power and total heat extracted for different fluid viscosities: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)


### 5.2.5. Effect of Geothermal Gradient on the Output Performance

The output temperatures under different gradient geothermal gradients,A0, are displayed in Figure 11 when other parameters are kept the same. It can be found that (1) as the depth of the vertical well is the same, larger geothermal gradient produces larger maximum output temperature (MOT), which is increased by around 42.2°C if A0 has an increase of 1.5°C. For example, the MOTs for the cases withA0 = 0.020, 0.035, 0.050, 0.065, 0.080, and 0.095 K/m are 76.2, 118.4, 160.5, 202.7, 244.9, and 287.0°C, respectively. However, the output temperature after a long time (say 30 years) is found to be increased by around 4°C whenA0 has an increase of 1.5°C, as shown in Figure 11(b), where the output temperature is 34.2, 38.2, 42.1, 46.0, 49.9, and 53.8°C, respectively, for the above six cases. This means that, compared to the injection conditions, the geothermal gradient plays much less important role in the output performance.Figure 11
Variation with time of output temperature of the production well under different geothermal gradientA0: (a) for large time period and (b) for short time period.
(a)
(b)


### 5.2.6. Fully Insulated Case

Figure12 compares the production temperatures under fully contacted and fully insulated conditions. The parameters used for the fully insulated case are the same as those in Figure 4 except that the overall heat transfer coefficient along the production well is set to zero.Figure 12
Production temperature under fully contacted and fully insulated conditions along the production well whenQ=0.05 m3/s, L1=3500 m, and L2=6000 m.It is found from Figure12 that there is no large difference in production temperature between these two cases and the production temperature predicted from fully insulated case is a little smaller than that predicted from fully contacted case. From the curves for t=10 days and t=30 years in Figure 4(b) without any casing and cementing, it can be found that (1) the final output temperature decreases very quickly and approaches some pseudo-steady value which is dominated by the inject conditions and (2) the fluid along a larger part (around 80 percent) of wellbore length absorbs heat from the surrounding rock formation and the fluid along the rest of the wellbore releases heat into the formation. Compared with this openhole case, if an insulated layer is put along the whole wellbore, the final output temperature will become smaller although the heat loss into the formation, which occurs along a small part of wellbore length in the openhole case, is prevented in the insulated case, as confirmed in Figure 12.


## 5.2.1. Temperature along the Vertical Wells

Figure4 displays the dimensionless temperature changes, which are defined by physical temperature changes with respect to the initial state divided by (A0H+B0), along the injection and production wells, respectively. The red curves refer to the initial state. For the injection well, we can see from Figure 4(a) that (1) before the cold fluid reaches the bottom of the wellbore, the temperature change of the fluid on the upper part of the wellbore decreases linearly with depth, while that on the lower part of the wellbore is almost a constant; (2) the temperature along the injection well decreases very quickly, with a temperature of 46.3°C at the bottom after 3 hours of circulation; and (3) as the injection temperature is lower than the initial ground temperature for most of the well length, most of the fluid flowing downwards absorbs heat from the surrounding rock except a small portion of fluid on the upper part of the wellbore.Figure 4
Temperature change profiles along the (a) injection and (b) production wells whenQ=0.05 m3/s,L1 = 3500 m, andL2 = 6000 m.
(a)
(b)Before the cold fluid reaches the production well, the temperature along the whole production well increases because of hot fluid flowing upwards, as shown by solid and dashed curves in Figure4(b). The production temperature increases from 57°C to 152°C when the time varies from 10 mins to 1.0 hours and then decreases to 126°C after 3 hours. The small difference in temperature change between times t=10 days and t=30 years means that the heat transfer in the fluid approaches pseudo-steady state after 10 days. Moreover, the critical position which indicates if the fluid absorbs heat from or releases heat to the rock formation is around Zc=0.2 or zc=700 m.


## 5.2.2. Effect of Flow Rates on Output Performance

Figure5 presents the output temperature for different injection rates ranging from Q=0.01 m3/s to 0.14 m3/s. It can be seen that the injection rate plays a major influence on the output temperature. First, for each injection rate, the output temperature increases quickly to a maximum value, as shown in Figure 5(b), and then decreases steadily to a constant value, as shown in Figure 5(a). Second, with increasing the injection rate, the production temperature is decreased. This is due to shortened time as a result of larger fluid speed for heat exchange between fluid and the surrounding rock. In addition, the temperature difference after 30 years between the cases Q=0.01 m3/s and Q=0.14 m3/s is around 40°C. However, when the injection rate is increased to some value (0.05 m3/s for the present case), there is no large difference in the final production temperature, as shown in Figure 5(a).Figure 5
Variation with time of output temperature of the production well under different injection rates whenL1 = 3500 m andL2 = 6000 m: (a) for general view and (b) for local view.
(a)
(b)The thermal power and total heat extracted for the above cases are plotted in Figure6. According to the definition of thermal power in (12), thermal power is linearly proportional to the injection rate and temperature. However, as lower production temperature is obtained when the injection rate is increased, the curves for thermal power shown in Figure 6(b) do not exhibit similar trends to the production temperature in Figure 5(b). It has to be mentioned that the long-time thermal power ranges from 1.71 MWs to 2.21 MWs when Q is changed from 0.01 m3/s to 0.14 m3/s.Figure 6
Thermal power and total heat extracted for different injection rates: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)Figure6(c) displays the total heat extracted from the reservoir up to 30 years under different injection rates. The total heat extracted is expected to increase when the injection rate is increased. However, the total heat extracted after 30 years when Q=0.14 m3/s is only 32.75% larger than that when Q=0.01 m3/s. This means that the cost related to higher pressure required for larger flow rates has to be balanced with the total energy extracted.


## 5.2.3. Effect of Injection Temperatures on Output Performance

The effect of injection temperatures on output performance is displayed in Figures7 and 8. As for output temperature, it increases evenly with increasing the injection temperature, as shown in Figure 7(a), and the maximum value of output temperature is not sensitive to the injection temperature, as shown in Figure 7(b).Figure 7
Variation with time of output temperature of the production well under different injection temperatures: (a) for large time period and (b) for short time period.
(a)
(b)Figure 8
Thermal power and total heat extracted for different injection temperatures: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)Figure8 presents the thermal power and total heat extracted up to 30 years. As the thermal power is linearly proportional to the temperature and the injection rate is kept constant, the thermal power exhibits a similar trend to the output temperature, as shown in Figures 8 and 7. In addition, it is interesting to find that although higher output temperature is obtained under larger injection temperature, the total heat extracted from the reservoir is reduced, as shown in Figure 8(c).


## 5.2.4. Effect of Fluid Viscosity on Output Performance

Figures9 and 10 display the final output performance for different fluid viscosities. The water viscosity ranges from about 2.8 × 10−4 Pa·s to 1.3 × 10−3 Pa·s when the temperature changes from 100°C to 10°C, which can be found in the website link (https://en.wikipedia.org/wiki/viscosity). Here μ = 10−4, 5 × 10−4, 2.5 × 10−3, and 1.2 × 10−2 Pa·s are chosen to show the effect of the viscosity on the output temperature. From Figure 9 the output temperature follows almost the same curve under these four cases and the decreasing trend for rapid temperature drop in less than 5 hrs is insensitive to the fluid viscosity.Figure 9
Variation with time of output temperature of the production well under different viscositiesμ: (a) for large time period and (b) for short time period.
(a)
(b)Figure 10
Thermal power and total heat extracted for different fluid viscosities: thermal power (a) for large time and (b) for short time and (c) total heat extracted.
(a)
(b)
(c)


## 5.2.5. Effect of Geothermal Gradient on the Output Performance

The output temperatures under different gradient geothermal gradients,A0, are displayed in Figure 11 when other parameters are kept the same. It can be found that (1) as the depth of the vertical well is the same, larger geothermal gradient produces larger maximum output temperature (MOT), which is increased by around 42.2°C if A0 has an increase of 1.5°C. For example, the MOTs for the cases withA0 = 0.020, 0.035, 0.050, 0.065, 0.080, and 0.095 K/m are 76.2, 118.4, 160.5, 202.7, 244.9, and 287.0°C, respectively. However, the output temperature after a long time (say 30 years) is found to be increased by around 4°C whenA0 has an increase of 1.5°C, as shown in Figure 11(b), where the output temperature is 34.2, 38.2, 42.1, 46.0, 49.9, and 53.8°C, respectively, for the above six cases. This means that, compared to the injection conditions, the geothermal gradient plays much less important role in the output performance.Figure 11
Variation with time of output temperature of the production well under different geothermal gradientA0: (a) for large time period and (b) for short time period.
(a)
(b)


## 5.2.6. Fully Insulated Case

Figure12 compares the production temperatures under fully contacted and fully insulated conditions. The parameters used for the fully insulated case are the same as those in Figure 4 except that the overall heat transfer coefficient along the production well is set to zero.Figure 12
Production temperature under fully contacted and fully insulated conditions along the production well whenQ=0.05 m3/s, L1=3500 m, and L2=6000 m.It is found from Figure12 that there is no large difference in production temperature between these two cases and the production temperature predicted from fully insulated case is a little smaller than that predicted from fully contacted case. From the curves for t=10 days and t=30 years in Figure 4(b) without any casing and cementing, it can be found that (1) the final output temperature decreases very quickly and approaches some pseudo-steady value which is dominated by the inject conditions and (2) the fluid along a larger part (around 80 percent) of wellbore length absorbs heat from the surrounding rock formation and the fluid along the rest of the wellbore releases heat into the formation. Compared with this openhole case, if an insulated layer is put along the whole wellbore, the final output temperature will become smaller although the heat loss into the formation, which occurs along a small part of wellbore length in the openhole case, is prevented in the insulated case, as confirmed in Figure 12.


## 6. Conclusion

This paper deals with the heat extraction from a closed-loop geothermal system. Through reasonable assumptions and the temperature continuity conditions at the intersection, the whole well is divided into three portions, two vertical wellbores and one horizontal wellbore, which can be solved independently based on the input conditions. Conclusions that can be made based on our numerical results are as follows:(1)
The difference in the output temperature predicted from the proposed model and TOUGH2 is very small, less than 4°C in the present numerical study, and thus can be used to approximately predict the heat extraction from the closed-loop geothermal system.(2)
The injection rate plays the dominant role in affecting the output performance; although it can increase the thermal power, the reduced temperature as a result of high flow rate compromises its production performance.(3)
Higher injection temperature produces larger output temperature but decreases the total heat extracted given a specific time.(4)
The output performance of geothermal reservoir is insensitive to fluid viscosity.(5)
There exists a critical point that indicates if the fluid above and below this point releases heat into or absorbs heat from surrounding formation.(6)
This approximation model proposed in this work runs in terms of seconds on a personal notebook computer and thus provides an efficient tool for reservoir optimization.



---

*Source: 2041072-2017-08-16.xml*