Repurposing a Geothermal Exploration Well as a Deep Borehole Heat Exchanger: Understanding Long-Term Effects of Lithological Layering, Flow Direction, and Circulation Flow Rate

Abstract

In the drive to achieve net-zero carbon emissions, decarbonisation of heating is essential. This can be facilitated by geothermal energy, but drilling geothermal wells is associated with high risks and costs. The use of preexisting wells (e.g., exhausted hydrocarbon wells or failed geothermal exploration boreholes) offsets this cost while potentially turning liabilities into assets. The Newcastle Science Central Deep Geothermal Borehole (NSCDGB) is a geothermal exploration well that was drilled to target the Carboniferous Fell Sandstone Formation at 1418.5 to 1795 m depth. However, low hydraulic conductivities prevented the development as a conventional “wet” geothermal abstraction well; therefore, new alternative methods of development are being explored. This work investigates the repurposing of the NSCDGB as a deep borehole heat exchanger (DBHE), focusing on the sustainable operation of the system in the long term by employing a constant heat load designed to contribute to local buildings or a heat network. Numerical modelling was undertaken by using OpenGeoSys software to analyse the thermal and hydraulic performance of the system. Both homogeneous and heterogeneous models were developed to compare the influence of lithological layering in contrast to a homogeneous (nonstratified) subsurface geological model. Results from homogeneous simulations modelling the DBHE to a depth of 922 m show that a 50-kW heat load can be supported for a lifetime of 25 years. This corresponds to a 65-kW building load when coupled to a heat pump with a coefficient of performance of 4.33. Thus, the DBHE could meet up to 72% of the heat demand of the adjacent urban sciences building. Rather than being a purely hypothetical case study, this work considers a real existing borehole, adjacent to a building cluster which could make use of the geothermal heat. Heterogeneity, which has been considered for the first time at the NSCDGB site, exhibits a minor impact in comparison to homogeneous simulation results. Flow direction and mass flow rate also exhibited small effects on the system performance, whereas if the exploration well could be repurposed to increased depths, the heat load could be increased. This is the first study of a coaxial DBHE at the NSCDGB site considering long-term effects of mass flow rate, heterogeneity, and flow direction. The study evaluates the feasibility of repurposing an exploratory geothermal well in the UK as a DBHE that can be used as a low-carbon heat source for space heating, thus converting liabilities into potential “green energy” assets.

1. Introduction

Building heating is the third major source of greenhouse gas emissions in the UK after surface transport and industry, accounting for 17% of the greenhouse gas emissions [1]. Thus, building heating, which currently relies on fossil fuels, must be decarbonised to achieve net-zero emissions as agreed in the Glasgow Climate Pact [2,3,4]. The use of renewable energy technologies, such as ground source heat pumps, has the potential to replace fossil fuel-based heating [5,6]. These heat pumps are often used in combination with an array of underground borehole heat exchangers up to 300 m in depth [7]. Going deeper (>500 m) makes more geothermal heat available with fewer spatial requirements [8,9]. This is particularly relevant in urban cities with limited space and high heating demand [10]. However, deep borehole heat exchangers are still not commonly used due to high initial drilling costs, issues with integrity of materials at depth, lack of experimental systems, hydraulic issues, and poor construction technologies [11,12,13].

One way to offset capital expenditure during drilling is to use an existing borehole, such as a failed exploration well, or decommissioned oil and gas well [14]. By doing so, an existing asset is leveraged while also aiding in climate change mitigation. Several studies [15,16,17,18,19,20,21,22,23] have investigated the possibility of repurposing an existing borehole for geothermal heat extraction by using a closed-loop deep borehole heat exchanger (DBHE). According to [24], flow rate and geothermal gradient are key parameters in repurposing oil and gas wells for geothermal energy extraction. This is because heat conduction influences the heat transfer fluid at lower temperatures, and increased geothermal gradients result in higher net power. Building on measurements from the existing 2.5-km borehole in Aachen, Dijkshoorn et al. [25] studied a deep coaxial borehole heat exchanger installed to meet the cooling/heating demands of a university building at RWTH Aachen. Results indicated that the heating requirements during winter were met, but during summer, the building used a climatic control adsorption unit with a minimum temperature of 55 ${}^{\circ }$C. This was impossible to achieve in the long term (20 years). Furthermore, better insulation of the inner piping was recommended, but due to the high additional cost implications, the operation of the system was declared to be economically unfeasible [25,26].

DBHEs commonly use a coaxial tube configuration (instead of a U-tube) due to better hydraulic efficiency and ease of installation [27]. This comes as a result of higher surface area and volume of fluid flow compared to U-tube heat exchangers. Moreover, thermal resistance between the borehole and the circulating fluid is also reduced by using the coaxial design [16,27]. The experimental studies on deep coaxial borehole heat exchangers are limited, with recent studies focusing on an installed plant in Xi’an, China (e.g., [8,28]). Several numerical approaches have been used to model subsurface heat transfer in deep coaxial borehole heat exchangers such as finite element method (FEM) [29], finite volume method (FVM) [30], and finite difference method (FDM) [10,21]. Deng et al. [28] studied a 2.5-km-deep borehole heat exchanger in Xi’an, China by using experiments and a finite volume heat transfer model. They grouped factors affecting system performance into external factors (geothermal gradient and thermal conductivity of the ground), internal factors (DBHE depth and thermal conductivities of inner and outer pipes) and synergic adjustments (flow rate and inlet water temperature). They showed that increasing the thermal conductivity of the outer pipe above 10 W/(m·K) has no significant effect because thermal resistance of the ground dominates. However, increasing the heat transfer area of the outer pipe will reduce the ground thermal resistance, which enhances heat transfer. Other studies have shown that insulating the inner pipe and increasing its velocity (and flow rate) by reducing the inner pipe diameter reduces heat loss. In the latter case, however, the pumping power requirement also increases, thereby necessitating a tradeoff [10,31,32].

In this work, we take the Newcastle Science Central Deep Geothermal Borehole (NSCDGB) as a case study [33] (Figure 1). The NSCDGB is currently located within the Newcastle Helix area [34] which is a hybrid state-of-the-art quarter at the heart of Newcastle. The borehole was finished in 2014 and intended to target a conventional geothermal system in which fluid from the formation is extracted to add to the energy mix of the area. A conventional system was explored as it can often contribute more energy in comparison to DBHEs and provide greater financial reward (e.g., [35,36]). After drilling, a high temperature was observed at the underlying Fell Sandstone Formation, but the hydraulic conductivity was low, making it impossible to obtain reasonable flow rates (less than 0.1 L/s [33]). Therefore, the EPSRC-funded project “NetZero GeoRDIE–Net Zero Geothermal Research for District Infrastructure Engineering” (Grant number EP/T022825/1) embarked on exploring the potential of repurposing the borehole as a closed-loop coaxial DBHE [37]. This work is carried out as part of the project. The nearby demand from local buildings and high heat flow associated with the geothermal resource constrained to the radiogenic North Pennine Batholith make this an ideal candidate for repurposing (Figure 2) [38,39]. Because the NSCDGB was found to be unsuitable as a conventional geothermal well, it could potentially be repurposed as a DBHE for space heating in the adjacent urban sciences building. This research models the operation of the borehole for space heating by using a numerical model.

Although there are studies in the literature investigating DBHEs, few of these studies look at long-term performance for using constant heat loads as a contribution to heating for the lifetime of a project and usually focus on intermittent operation (e.g., [21,29,40,41,42]). Many of these studies focus on homogeneous geology and do not consider the effects of heterogeneous geological layers in comparison. The study by [21] considered long-term operation but employed constant inlet temperature for the simulations. Studies by [29,41,42] considered constant heat-load operation but lithological layering of the rocks was not considered. Moreover, modelling studies evaluating the repurposing of geothermal exploration wells are few, whereas in studies that consider repurposing oil and gas wells (e.g., [43]), lithological layering is not considered. This work is of particular importance as studies looking at repurposing a geothermal exploration well are rare and this provides an opportunity to reduce drilling risk by offering a method of exploiting a resource if failure occurs during exploration. Rather than being a purely hypothetical case study, this work considers a real existing borehole, adjacent to a building cluster which could make use of the geothermal heat. The aim of this work is to therefore study the long-term operation of the NSCDGB when repurposed as a DBHE. Numerical modelling is undertaken by using the open-source, finite element code OpenGeoSys (OGS). A range of simulations are performed to test the most important parameters as identified in the literature; these include the heat load, depth, and mass flow rate [21,28]. With this modelling study, a predictive hypothesis regarding the borehole’s thermal behaviour during heat extraction is generated. This hypothesis has the possibility of being tested subsequently against empirical data from a real planned thermal response test (TRT) [37]. This study thus provides a case study of how well predictive tools stand up against future thermal response tests.

2. Methods

In the “dual-continuum” FEM modelling approach, which is adopted in this work, the DBHE is discretised by using one-dimensional (1D) finite elements, while the surrounding rock is discretised by using three-dimensional (3D) finite elements [44,45]. A thermal capacitor resistor network is used to estimate the heat transfer rate [45]. This modelling approach has been implemented in OpenGeoSys (OGS) [29,46,47], which is an open-source finite element multiphysics simulator. OGS version 6 is employed in this work.

2.1. Governing Equations

Heat can be extracted with a coaxial DBHE by using two configurations: (a) the CXA configuration (shown in Figure 3 left) where cold circulating fluid flows in down the annulus and leaves through the central coaxial pipe after receiving heat through the borehole walls; and (b) the CXC configuration (shown in Figure 3 right), wherein the fluid flows in down the central pipe and leaves through the annulus. For the governing equations, we consider the CXA configuration. There is usually a grout material between the outer or annular pipe and the surrounding rock for better stability of the borehole and to prevent uncontrolled fluid flow behind the casing. There are four governing equations for the heat transfer problem between the DBHE and the surrounding rock: one governing equation each for the rock, the grout, the inner pipe (central coaxial pipe), and the outer pipe (annular pipe). For the rock and grout, conduction is taken to be the dominant heat transfer mechanism, while convection dominates the annular and central pipes.

The governing equation for the rock is given by [29,48]

$$\frac{\partial }{\partial t}\left[\varphi {\rho }_f{c}_f+(1-\varphi ){\rho }_s{c}_s\right]T_s-\nabla ·({\mathsf{\Lambda }}_s·\nabla T_s)=H_s$$

(1)

with $\varphi$, the solid rock porosity, ${\rho }_f$, the fluid density, $c_f$, the fluid specific heat capacity, ${\rho }_s$, the solid rock density, $c_s$, the solid rock specific heat capacity, $T_s$, the solid rock temperature, ${\mathsf{\Lambda }}_s$, the bulk saturated effective thermal conductivity, and $H_s$, the source or sink term. The heat flux between the solid rock and the DBHE is expressed as

$$q_{n{T}_s}=-({\mathsf{\Lambda }}_s·\nabla T_s).$$

(2)

The governing equation for the grout is [49]

$$(1-{\varphi }_g){\rho }_g{c}_g\frac{\partial T_g}{\partial t}-\nabla ·\left[(1-{\varphi }_g){\lambda }_g·\nabla T_g)\right]=H_g\mathrm{in}{\mathsf{\Omega }}_g,$$

(3)

in which ${\varphi }_g$ is the grout porosity, ${\rho }_g$ is the grout density, $c_g$ is the specific heat capacity of the grout, ${\lambda }_g$ is the thermal conductivity of the grout, $T_g$ is the grout temperature, $H_g$ is the source/sink term and ${\mathsf{\Omega }}_g$ is the grout domain. The following heat flux boundary condition is used,

$$q_{n{T}_g}=-{\mathsf{\Phi }}_{gs}(T_s-T_g)-{\mathsf{\Phi }}_{fig}(T_{in}-T_g)\mathrm{on}{\mathsf{\Gamma }}_g,$$

(4)

in which the heat transfer coefficient between the grout and rock, ${\mathsf{\Phi }}_{gs}$, is a function of borehole diameter D and the thermal resistance between the grout and solid rock $R_{gs}$ [48]:

$${\mathsf{\Phi }}_{gs}=\frac{1}{R_{gs}\pi D}.$$

(5)

The annular inflow is governed by

$${\rho }_f{c}_f\frac{\partial T_{in}}{\partial t}+{\rho }_f{c}_f{\mathit{v}}_i·\nabla T_{in}-\nabla ·({\mathsf{\Lambda }}_f·T_{in})=H_i\mathrm{in}{\mathsf{\Omega }}_i,$$

(6)

with $T_{in}$, the fluid inflow temperature, ${\mathit{v}}_i$, the inflow fluid velocity, $H_i$, the source/sink term, ${\mathsf{\Omega }}_i$, the inflow domain, and ${\Lambda }_f$, the hydrodynamic thermodispersion tensor expressed as

$${\mathsf{\Lambda }}_f=({\lambda }_f+{\rho }_f{c}_f{\alpha }_L\parallel v_f\parallel )\mathit{I}.$$

(7)

The following boundary condition is adopted,

$$q_{n{T}_{in}}=-{\mathsf{\Phi }}_{fig}(T_s-T_{in})-{\mathsf{\Lambda }}_{ff}(T_{out}-T_{in})\mathrm{on}{\mathsf{\Gamma }}_i,$$

(8)

in which $T_{out}$ is the outflow temperature, ${\mathsf{\Phi }}_{fig}$ is the heat transfer coefficient between the inflow (annular) pipe and grout, and ${\mathsf{\Phi }}_{ff}$ is the heat transfer coefficient between the inflow (annular) pipe and the outflow (central coaxial) pipe. These heat-transfer coefficients are related to their respective thermal resistances as [48]

$${\mathsf{\Phi }}_{fig}=\frac{1}{R_{fig}\pi d_{casing}}$$

(9)

$${\mathsf{\Phi }}_{ff}=\frac{1}{R_{ff}\pi d_{central}},$$

(10)

where $d_{casing}$ is the borehole (annular) casing diameter and $d_{central}$ is the central coaxial pipe diameter (see [45] for how to compute the thermal resistances). The fluid outflow governing equation is

$${\rho }_f{c}_f\frac{\partial T_{out}}{\partial t}+{\rho }_f{c}_f{\mathit{v}}_f·\nabla T_{out}-\nabla ·({\mathsf{\Lambda }}_f·T_{out})=H_o\mathrm{in}{\mathsf{\Omega }}_o,$$

(11)

with boundary condition

$$q_{n{T}_{out}}=-{\mathsf{\Phi }}_{ff}(T_{in}-T_{out})\mathrm{on}{\mathsf{\Gamma }}_o.$$

(12)

The finite element method (FEM) is employed to solve the four governing Equations: (1), (3), (6) and (11). The backward Euler finite difference is used for time discretisation [45]. In the finite element solution procedure, a weak form of the governing equations is obtained by using a test function. Next, the obtained result is integrated by parts. The discretised primary variables at the nodes are then used to obtain the solution to the problem by using shape or basis functions [45]. OGS adopts this FEM implementation.

2.2. Incorporating a Heat Pump

The heat load on a DBHE (${\dot{Q}}_{\mathrm{DBHE}}$) is the same as the building heat load (${\dot{Q}}_{\mathrm{building}}$) it feeds when there is no heat pump. However, a heat pump is often used to elevate the fluid temperature from the DBHE before it is used for heating a building space. Hence, with a heat pump, ${\dot{Q}}_{\mathrm{building}}$ is not the same as ${\dot{Q}}_{\mathrm{DBHE}}$. The performance of the heat pump is estimated by using the coefficient of performance (COP) which is given by [29]

$$\mathrm{COP}=\frac{{\dot{Q}}_{\mathrm{building}}}{W_{\mathrm{hp}}},$$

(13)

in which $W_{\mathrm{hp}}$ is the electricity consumed by the heat pump. The heat extracted from the DBHE or heat load on the DBHE is the building load less the electricity consumption of the heat pump:

$${\dot{Q}}_{\mathrm{DBHE}}={\dot{Q}}_{\mathrm{building}}-W_{\mathrm{hp}}.$$

(14)

By using Equation (13), a relationship between ${\dot{Q}}_{\mathrm{DBHE}}$ and ${\dot{Q}}_{\mathrm{building}}$ is obtained as follows:

$${\dot{Q}}_{\mathrm{building}}=\frac{\mathrm{COP}}{\mathrm{COP}-1}{\dot{Q}}_{\mathrm{DBHE}}.$$

(15)

To obtain initial estimates for the NSCDGB, a linear relationship for COP which depends on the outflow temperature ($T_{out}$) of the DBHE is adopted here,

$$\mathrm{COP}=a{T}_{out}+b,$$

(16)

where $a=0.083$ and $b=3.925$ assume a heating temperature of 35 ${}^{\circ }$C [29,50], which is considered to be a suitable outlet temperature for the local demand in the Newcastle Helix area.

A deep borehole heat exchanger (BHE) has higher hydraulic losses compared to a shallow BHE. Thus, an electrically powered circulation pump is needed to overcome the pressure drop ($\Delta P$) calculated as [9]

$$\Delta P=\frac{L{\rho }_f{v}_f^2}{2D_h{[0.79\ln \left(Re\right)-1.64]}^2},$$

(17)

where L is the length of the DBHE, $v_f$ is the fluid velocity, $D_h$ is the hydraulic diameter of pipe and $Re$ is the turbulent flow Reynolds number dependent on the hydraulic diameter. The pressure drop can be calculated for inflow or outflow and the corresponding pipe (annular or central, respectively) is used to obtain the hydraulic diameter and Reynolds number. The pressure drop in the central (outflow) pipe plus the pressure drop in the inflow (annular) pipe has been used in this work, and all flows have been confirmed to be turbulent.

The electricity consumed by the circulating pump ($W_{cp}$) is computed by using [29]:

$$W_{\mathrm{cp}}=\frac{\dot{V}\Delta P}{\eta },$$

(18)

in which $\dot{V}$ is the volumetric flow rate of the circulating fluid, and $\eta$ is the efficiency of the circulating pump, assumed here to be 70% [29]. The system performance can be evaluated by using a coefficient of system performance (CSP), which takes into account the electricity consumption of both the heat pump and the circulating pump [29]:

$$\mathrm{CSP}=\frac{{\dot{Q}}_{\mathrm{building}}}{W_{\mathrm{hp}}+W_{\mathrm{cp}}}.$$

(19)

2.3. Parameter Outline for Newcastle Science Central Deep Geothermal Borehole

The NSCDGB is located to the northeast of England in Newcastle where heat flows are elevated due to the presence of the North Pennine Batholith (see Figure 2). The 1821-m-deep borehole cuts a clastic and carbonate Carboniferous succession of strata [33]. The borehole has a 7-inch casing stopping at 941 m and a hole of 6-inch nominal diameter extending to the total depth. There are other, outer concentrically grouted casing strings associated with the real borehole, but, for simplicity, these have been ignored in the model. From 922 m to 1651 m, a 4.5-inch Techniseal liner has been inserted in 2014 for pump tests. This lower portion of the well is considered too narrow to realistically install a DBHE. Therefore, the initial scenario only considers the DBHE being installed in the well to a depth of 922 m.

To develop a heterogeneous model of the geological strata, relevant layer data are required. These include thickness, thermal conductivity, density and specific heat capacity for each layer. Starting with the Pennine Middle Coal Measures Formation (161 m thickness) at the ground surface, the borehole penetrates the Pennine Lower Coal Measures Formation (157 m), Millstone Grit (58 m), Yordale Group (1042.5 m) and Border Group (401.5 m) as listed in [33]. Within the Yordale Group, it cuts through the Stainmore Formation (293 m), Alston Formation (191 m), a fault zone (68 m) and Tyne Limestone Formation (490.5 m) [33]. The Border Group contains the Fell Sandstone Formation (376.5 m) and ends with the Lyne Formation (25 m). Figure 4 shows the geological layers from the ground surface to a depth of 1651 m at the Fell Sandstone Formation. Below 1651 m, there are abandoned rock bits in the borehole. Within the Newcastle Helix area [34], the NSCDGB is located adjacent to Newcastle University buildings such as the urban sciences building. These buildings could provide a proximal demand for the extracted heat or alternatively, the extracted heat could be fed into a newly developed heat network.

The work by Younger et al. [33] outlined formation parameters such as thermal conductivity particularly from a depth of 500 m to the end of the borehole but has little data on other parameters needed for modelling such as density and specific heat capacity. Thus, additional sources have been used to compile the needed data for the complete set of parameters of the borehole and are listed in

Table 1. Thermo-physical modelling parameters for the geological layers of the NSCDGB. Sources: [33,51,52,53,54,55,56,57,58,59,60]. Thermal conductivities, densities and specific heat capacities in this table are believed to be representative of bulk, water-saturated values. For input to the model, porosities and mineral matrix densities/specific heat capacities have been selected to result in the appropriate bulk values.

Formation	End Depth [m]	Thickness [m]	Thermal Conductivity [W/m/${}^{\circ }$C]	Density [kg/m3]	Specific Heat [J/kg/${}^{\circ }$C]
Pennine Middle Coal Measures	161	161	2.35	2451.09	857.40
Pennine Lower Coal Measures	318	157	2.35	2451.09	857.40
Millstone Grit Group	376	58	2.90	2460	930
Stainmore Formation	669	293	2.61	2500	1000
Alston Formation	860	191	2.63	2500	1000
Fault Zone	928	68	2.75	2500	1000
Tyne Limestone	1418.50	490.50	2.86	2460	930
Fell Sandstone	1651	232.50	2.91	2460	930
Homogenised average (1651 m)	1651	1651	2.69	2472	939
Homogenised average (922 m)	922	922	2.55	2480	946
Homogenised average (500 m)	500	500	2.48	2464	901

[51,52]. The corresponding homogeneous model properties based on averaging over all layers up to 1651 m and other depths considered (500 m and 922 m) are also listed. These are applied for homogeneous simulations in Section 3.1.

To model the NSCDGB, a lateral rock domain of 200 m × 200 m is used, extending to a depth of 1807 m. The boundary distances were tested prior to the study to ensure no detrimental impacts were caused by boundary influence. The 922-m DBHE is centralised in the domain. A borehole diameter of 0.216 m is considered, and a mass flow rate of 8.33 kg/s is used [51]. The heterogeneous finite element mesh is shown in Figure 5. Other parameters are listed in

Table 2. Parameters for the NSCDGB as a DBHE [51,52].

Parameter	Symbol	Value	Unit
Model domain		200 × 200 × 1807	m3
Borehole depth		1651	m
Borehole diameter	D	0.216	m
Annular pipe outer diameter	$d_{casing}$	0.1779	m
Central pipe outer diameter	$d_{central}$	0.1005	m
Annular pipe wall thickness		0.0081	m
Central pipe wall thickness		0.00688	m
Surface temperature		9	∘C
Reference temperature		9	∘C
Fluid flow rate	$\dot{V}$	0.00833	m3/s
Fluid dynamic viscosity		0.0008	Pa·s
Geothermal gradient		33.4	∘C/km
Fluid volumetric heat capacity	${\rho }_f{c}_f$	$4.1581\times {10}^6$	J/m3/∘C
Grout volumetric heat capacity	${\rho }_g{c}_g$	$1.824\times {10}^6$	J/m${}^3$/${}^{\circ }$C
Fluid thermal conductivity	${\lambda }_f$	0.59	W/m/∘C
Grout thermal conductivity	${\lambda }_g$	1.05	W/m/${}^{\circ }$C
Thermal conductivity of inner pipe (HDPE)		0.45	W/m/${}^{\circ }$C
Thermal conductivity of outer pipe (steel)		52.7	W/m/${}^{\circ }$C

. Freshwater is assumed as the circulating fluid. An annular steel pipe is adopted, and the central coaxial pipe is assumed to be made of high-density polyethylene [51]. An API Grade G water-saturated cement is used as grout material. In OGS, the reference temperature is the initial temperature of the fluid in the central coaxial pipe when using a boundary condition with specified power (heat load).

2.4. Boundary and Initial Conditions

The initial temperature ($T_0$) of the formation (in ${}^{\circ }$C), depends on the surface temperature and the geothermal gradient

$$T_0\left(z\right)=9-0.0334z,$$

(20)

where z (negative downwards) represents depth. As sourced from [56], a surface temperature of 9 ${}^{\circ }$C and a linear geothermal gradient of 33.4 ${}^{\circ }$C/km have been used similar to [61]. At all the domain boundaries (top, bottom and four lateral), a constant temperature Dirichlet boundary condition has also been used following equation (20). The initial conditions of the DBHE are also assumed to follow equation (20).

2.5. Benchmarking

Several studies have verified OGS as a reliable tool for modelling DBHEs [29,61]. Herein, a case of the NSCDGB with a 922-m DBHE is considered. This is because the current installation, which has a Techniseal liner, will only allow a DBHE to this depth; it is hydraulically unattractive to operate at depths exceeding this due to the narrow diameter [51]. Extraction of 50 kW is simulated for 25 years, which is representative of a long-term operation [29,41,62]. Results are compared with the analytical model by Beier [63] which was developed for DBHEs fully accounting for a geothermal gradient. The inlet and outlet temperatures at the DBHE top are shown on Figure 6. A very good comparison is observed between the analytical solution of Beier and the finite element solution of OGS. At 25 years, the outlet temperatures are 3.83 ${}^{\circ }$C and 3.46 ${}^{\circ }$C for Beier and OGS, respectively, a temperature difference of 0.37 ${}^{\circ }$C. This small discrepancy may be traced to the different approaches between the models for estimating Reynolds and Nusselt numbers, or due to the fact that Beier’s solution is designed to simulate short-term thermal response tests. The results show the capability of OGS to model a DBHE system operation to a reasonable level of accuracy and justifies its adoption for investigating long-term performance of the NSCDGB as a DBHE.

3. Numerical Results

In this section, analysis of the Newcastle Science Central Deep Geothermal Borehole (NSCDGB) as a DBHE is first performed by using a homogeneous model with averaged rock properties as listed in Table 1. The results are then analysed in comparison with a vertically heterogeneous setting taking the properties of different layers into account. From the analyses, an estimate is made of the amount of heat that can be extracted from the NSCDGB as a heat load applied directly to the DBHE. Depth and mass flow rate have significant effect on the operation of a DBHE [21,28,29]. Thus, a parametric study of these parameters is conducted; the effects of coaxial flow configuration and heterogeneity are also explored. Next, incorporation of a heat pump is considered in the operation of the DBHE. This affects the amount of load that can be supplied to a building.

3.1. Homogeneous Model

3.1.1. Heat Load and Temporal Analysis

As mentioned earlier, we have assumed that a DBHE can be installed only to a depth of 922 m, due to the 4.5-inch Techniseal tubing in the borehole beyond this depth. Results showing the fluid inflow ($T_{in}$) and outflow temperatures ($T_{out}$) for different applied heat loads for the DBHE of 922 m depth are presented in Figure 7. For all applied heat loads, there is a large temperature drop in the initial years of extraction. Although this period of large temperature drop is just two years for 50-kW extraction, it is over 10 years for a load of 150 kW. After this initial period of large temperature drop, a quasisteady state is reached. The figure also shows that 50 kW can be extracted for 25 years without the fluid inlet temperature going below 0 ${}^{\circ }$C. This corresponds to 10.95 GWh of energy extracted over 25 years. At higher heat loads, the fluid inlet temperature drops well below 0 ${}^{\circ }$C, which is the freezing point of water. For design considerations in this work, freshwater without any antifreeze constituents is being assumed. Hence, a load which causes the fluid inlet temperature to fall below 0 ${}^{\circ }$C cannot be supported. This shows that the DBHE can only support a load of 50 kW in the long term. Thus, subsequent simulations will focus on a 50-kW heat load.

Figure 8 shows the contours of rock temperature. At the initial time (0 days), the DBHE temperature is the same as the rock temperature. The extraction of heat from the rock is visible at 92 days within a limited radius of influence (see also Figure 17). The radius of influence develops further with time, and at 25 years the region of low temperature has covered the full depth of the DBHE. The corresponding inflow and outflow fluid temperatures are shown in Figure 9 along the DBHE depth. The temperature of the flowing fluids decrease with time as the available heat energy reduces due to extraction. The temperature drop is rapid in the early stage of extraction. The temperature at the DBHE bottom is 20 ${}^{\circ }$C after two days; this rapidly reduces to ∼13 ${}^{\circ }$C at 92 days, a 7 ${}^{\circ }$C drop. However, the drop from 92 days to 480 days is only ∼3 ${}^{\circ }$C. As explained earlier, after two years of extraction, the temperature modestly drops with time. Approximately 5 ${}^{\circ }$C drop is observed at the DBHE bottom between 480 days and 25 years of extraction.

3.1.2. Impact of Depth

Figure 10 shows the effect of different DBHE depths on $T_{in}$ and $T_{out}$. A 500-m DBHE cannot support the 50-kW heat load, whereas for 922 m and 1651 m, the load can be supported. The fluid inflow temperature at 25 years is 3.46 ${}^{\circ }$C and 24.84 ${}^{\circ }$C for DBHE depths of 922 m and 1651 m, respectively. For a depth of 1651 m, a much higher heat load can be supported. This is further investigated in Figure 11. With a 1651-m DBHE, a heat load of 100 kW can be supported without the fluid inlet temperature going below 12 ${}^{\circ }$C at 25 years. However, when a load of 150 kW is applied, the fluid inlet temperature reaches 0 ${}^{\circ }$C at 20 years, indicating that 150 kW can only be supported for 20 years. These results indicate that depth has a significant effect on the performance and sustainability of the DBHE system.

3.1.3. Impact of Flow Configuration

Using an annular inlet and central coaxial outlet (CXA) gives a slightly better performance compared to the reverse configuration, which uses a central inlet and annular outlet (CXC); however, the difference between the two configurations lessens with time. This is demonstrated in Figure 12. At 5 years, the inlet fluid temperature is 6 ${}^{\circ }$C for CXA and 5.58 ${}^{\circ }$C for CXC, a difference of 0.42 ${}^{\circ }$C. The difference continues to decrease until it reaches 0.04 ${}^{\circ }$C at 24 years and only 0.03 ${}^{\circ }$C at 25 years.

A plot of the inlet and outlet temperatures along the depth of the DBHE is presented in Figure 13 at 480 days and 25 years. The reason for the slightly improved performance in the CXA configuration is associated with the direction of fluid flow in comparison to the thermal gradient of the surrounding rock. As described by Banks [64], the flow of the heat-gaining fluid in the annular space is moving in the direction of the increasing temperature gradient in the surrounding formation, leading to more efficient heat transfer. The improved performance of CXA configuration is more clearly observed at 480 days whereas the difference is marginal at 25 years. The difference in temperatures between the different configurations at the DBHE bottom is more at 480 days compared to 25 years when a quasisteady state is approached. This explains why the difference between CXA and CXC is more pronounced in the initial years of extraction.

3.1.4. Impact of Heat Transfer Fluid Flow Rate

For the same heat load and the same fluid specific heat capacity, changing the mass flow rate will result in a change in the temperature differential between inflow and outflow fluid ($\Delta T$). As shown in Figure 14, a decrease in mass flow rate causes an increase in temperature drop for 50 kW. Thus, a very low mass flow rate can result in a large temperature drop that causes the fluid inlet temperature to go below 0 ${}^{\circ }$C. For all the mass flow rates considered, a fluid inlet temperature above 0 ${}^{\circ }$C is maintained, indicating that all values can technically be used. However, going below 3 kg/s might lead to fluid inlet temperatures falling below 0 ${}^{\circ }$C. A higher mass flow rate results in more sustainable use of the DBHE due to higher fluid inlet temperature. However, the cost of pumping as a function of flow rate is also an important design consideration.

Figure 15 shows the fluid temperatures with depth for the different mass flow rates at 25 years. The figure shows that with an increase in mass flow rate, the fluid temperature at the bottom of the DBHE lowers. The large convexity in the up-flow curve for a flow rate of 3 kg/s, means that at low flow rates, a lot of heat is being lost across the inner pipe to the down-flowing cold fluid (i.e., borehole thermal short-circuiting). Higher flow rates result in lower $\Delta T$. A lower $\Delta T$ means less thermal short circuiting between annular fluid and coaxial fluid across the inner pipe. This effectively means lower borehole thermal resistance. Both 5 kg/s and 8.33 kg/s show very close values for $T_{out}$ but different values for $T_{in}$. A high change in temperature usually results in lower efficiency. This, as well as the cost of pumping, is discussed further in Section 3.3.

3.2. Vertically Homogeneous vs. Heterogeneous Model

In this section, the stratigraphy of the rock layers is fully considered as a heterogeneous model. Heterogeneous here implies stratified, i.e., different strata are assigned different properties, and within each stratum, properties are still homogeneous. Results are compared to the previous section where a homogeneous rock layer has been used.

3.2.1. Impact of Depth

Figure 16 shows a comparison of the homogeneous and heterogeneous models for different depths. For all models, the heterogeneous model results in slightly higher temperatures for the fluid inlet ($T_{in}$) and outlet temperatures ($T_{out}$). Although the 922-m homogeneous model shows a $T_{in}$ value of 3.46 ${}^{\circ }$C at 25 years, the heterogeneous model shows a $T_{in}$ value of 3.54 ${}^{\circ }$C—just a 0.08 ${}^{\circ }$C difference. The difference is 0.29 ${}^{\circ }$C for the 500-m DBHE model and 0.92 ${}^{\circ }$C for the 1651-m DBHE model. The heterogeneous model records a higher temperature because the thermal conductivity of deep layers beyond 928 m reaches 2.9 W/m/${}^{\circ }C$ whereas the homogeneous model has an effective thermal conductivity of only 2.69 W/m/${}^{\circ }C$. At shallower depths, the layers have thermal conductivities very close to the average, which explains why the difference in fluid temperatures is significantly reduced.

Figure 17 depicts the rock temperature at selected depths taken radially from the 922 m DBHE. There is almost no distinguishable difference between the homogeneous and heterogeneous models. The radius of thermal influence appears to be around 60 m from the DBHE (measured to within 0.5 degree of initial conditions). The radius of influence increases with depth; although it is approximately 40 m at 161-m depth, it grows to 60 m at 908-m depth. The figure also shows that the temperature drop increases with depth. Increased heat exchange between the circulating fluid and the deeper rocks, which are hotter, leads to higher temperature drops with depth.

3.2.2. Impact of Flow Configuration

Similar to the homogeneous case, the CXA configuration has slightly higher fluid temperatures compared to the CXC configuration. Again, the heterogeneous model shows higher temperatures compared to the homogeneous model. However, the difference between the two configurations does not vanish with time in the heterogeneous model. At 25 years, the heterogeneous model estimates $T_{in}$ for the CXA and CXC configurations to be 3.54 ${}^{\circ }$C and 3.23 ${}^{\circ }$C, respectively (in Figure 18). These are higher than 3.46 ${}^{\circ }$C and 3.43 ${}^{\circ }$C estimated by the homogeneous model for the CXA and CXC configurations, respectively. This difference is demonstrated further with depth in Figure 19. The difference between the homogeneous and heterogeneous model for CXA configuration is small compared to CXC configuration. Looking at the CXC configuration, the heterogeneous model shows higher flow temperatures than the homogeneous model at an early stage of the simulation (480 days). However, at later stages (25 years), the flow temperatures in the homogeneous model are higher.

3.2.3. Impact of Heat Transfer Fluid Flow Rate

As shown in Figure 20, both the homogeneous and heterogeneous models demonstrate that for a mass flow rate of 3 kg/s, the DBHE cannot be used for more than 25 years because the fluid inlet temperature reduces to 0.4 ${}^{\circ }$C and 0.33 ${}^{\circ }$C, respectively, for the heterogeneous and homogeneous models. Generally, both models show similar results and estimates are within 1 ${}^{\circ }$C between the homogeneous and heterogeneous models for a sustainable load of 50 kW.

3.3. Incorporating a Heat Pump

For a 50 kW load, the homogeneous model shows a fluid outlet temperature of 4.90 ${}^{\circ }$C after 25 years whereas the heterogeneous model estimates the fluid outlet temperature to be 4.98 ${}^{\circ }$C. The results of the heterogeneous model are almost identical with the homogeneous model results. For the homogeneous model, this represents a heat pump COP of 4.33. Using Equation (15) results in a building load of 65.01 kW. Similarly, the heterogeneous model gives a COP of 4.34 and a building load of 64.98 kW. Although the DBHE can only support a heat extraction rate of 50 kW, together with a heat pump, it can supply ∼65 kW to a building. The difference between the homogeneous and heterogeneous models is marginal for this load. Both models predict an outlet temperature greater than 4 ${}^{\circ }$C which is required to maintain a reasonable COP [29]. Hence, the maximum DBHE load that can be supported when using a heat pump is 50 kW (equivalent to ∼65 kW building load). It is to be noted that even though heat pumps can operate below 4 ${}^{\circ }$C, there is the risk of formation of ice crystals within the evaporator below this temperature (unless an antifreeze or salt is added to the circulating fluid). With a high flow rate (low $\Delta T$), heat pumps can be used below this temperature or an antifreeze can be used. However, this is beyond the scope of this work.

Earlier, it was shown that different mass flow rates (5 kg/s and 8.33 kg/s) resulted in the same outlet temperature, $T_{out}$. This implies the same COP value but this will result in a difference in pressure drop and system performance (CSP). Hydraulic and overall system performance is demonstrated in Figure 21 for different mass flow rates. In Figure 21a, the pressure drop $\Delta P$ clearly increases with increasing mass flow rate. However, there is almost no increase in $T_{out}$ when the mass flow rate is increased beyond 5 kg/s; after 6 kg/s $T_{out}$ even starts to decrease. The increasing $\Delta P$ causes an increase in circulation pump electricity consumption as seen in Figure 21b. When the mass flow rate is increased from 5 kg/s to 6 kg/s, the COP remains almost constant (Figure 21c), but this leads to a decrease in system performance (Figure 21d) due to additional electricity consumption by the circulation pump. Increasing the mass flow rate beyond 6 kg/s causes a decrease in $T_{out}$ and COP, and hence the mass flow rate should not go beyond 6 kg/s. The optimal mass flow rate for operating the DBHE is approximately 5 kg/s and it should not exceed 6 kg/s. It is to be noted that an assumption of a constant viscosity has been maintained. However, fluid viscosity changes with temperature, and this could affect $\Delta P$.

4. Discussion

An OpenGeoSys model of the NSCDGB has been developed, incorporating the geological layers of the area. The model has been verified against Beier’s analytical solution [63], showing very good agreement with a difference below 0.4 ${}^{\circ }$C. This shows the reliability of the model developed to analyse the NSCDGB in the hypothetical case that it is repurposed as a deep borehole heat exchanger.

Comparisons between the homogeneous model that averages layer properties and the heterogeneous model that considers all layers show marginal differences generally in the range of 0.1–1 ${}^{\circ }$C for a load of 50 kW. This indicates that effects of heterogeneity are not very significant in the case of the NSCDGB when a conservative estimate of the heat load is applied. Results of the heterogeneous model generally show slightly better performance for all parameters studied due to slightly higher thermal conductivity in comparison with the homogeneous model. At a depth of 922 m, the thermal conductivity of the heterogeneous model is 2.63 W/m/${}^{\circ }C$ and a large portion of top layers have conductivity of 2.35 W/m/${}^{\circ }C$. This balances out and results in a homogenised thermal conductivity of 2.55 W/m/${}^{\circ }C$. Hence, the difference between the homogeneous and heterogeneous models is very minimal (0.08 ${}^{\circ }$C), and models show very similar results.

It has been shown that higher loads show longer periods of large temperature drops (see Figure 7 and Figure 11), and longer depths support more heat load. It can thus be inferred that the differences between homogeneous and heterogeneous models might generally increase with depth.

The temporal analysis over a period of 25 years shows that there is an initial period of large temperature drop before a quasisteady-state is reached. In the initial stages of the simulation, thermal resistance between the pipe and rock dominates, but with time, the thermal conductivity of the rock becomes dominant. At this later stage, the circulating fluid temperature (water) reaches equilibrium with the rock temperature. The radius of influence from the DBHE was found to be approximately 60 m. From the different DBHE loads considered, simulations showed that a load of 50 kW can be supported by a depth of 922 m for 25 years without the fluid inlet temperature going below 3.54 ${}^{\circ }$C based on the heterogeneous model. This is with respect to the direct load applied on the DBHE. When a heat pump with a COP of 4.34 is incorporated, the building load that can be supported is ∼65 kW. Sensitivity analyses were performed for a DBHE load of 50 kW to check the effect of depth, flow configuration and mass flow rate.

Depth has a significant effect on the long-term performance of the DBHE. For 922 m, the inlet temperature at 25 years of operation was ∼4 ${}^{\circ }$C compared to ∼25 ${}^{\circ }$C for 1651 m. This shows that installing a deeper DBHE will significantly increase the load supported by the DBHE (over 100 kW). However, based on the current site constraints, this depth cannot be reached. Additional work on the NSCDGB is required to go deeper and this is likely to come at significant cost.

By using the heterogeneous model, it was found that the CXA flow configuration with an annular inlet and central coaxial pipe outlet results in slightly better performance (0.31 ${}^{\circ }$C difference) than the CXC configuration for the DBHE, which is operating at extraction-only mode. This is in agreement with other studies in the literature [11]. For the homogeneous model, only a temperature difference of 0.03 ${}^{\circ }C$ was recorded after 25 years as the difference between the configurations reduced with time.

The mass flow rate is an operational parameter that can be tailored to suit operational conditions. Simulations showed that any mass flow rate above 3 kg/s will sustainably deliver 50 kW heat load over the 25-year period considered. However, with a higher mass flow rate, the temperature drop is reduced, thereby increasing the sustainable operational period of the DBHE. Increasing the mass flow rate beyond 5 kg/s did not increase the outflow temperature but increased the cost of circulating the DBHE fluid. This led to a reduction in the system performance. Hence, the optimum mass flow rate based on the CSP is around 5 kg/s. The maximum mass flow rate that should be used is around 6 kg/s because any value above this causes a reduction in outflow temperature, and hence the COP.

5. Conclusions

A heterogeneous finite element model developed in OpenGeoSys has been applied, for the first time, to study the feasibility of repurposing the NSCDGB as a coaxial deep borehole heat exchanger. The study has shown that the NSCDGB could be repurposed as a deep borehole heat exchanger, and although it cannot meet the heat demand of the building adjacent to the borehole [65], it can contribute to the Newcastle Helix heat network. A nonconservative analysis using the total depth of the borehole showed that more than 100 kW can be extracted over a 25-year period without the fluid inlet temperature going below 3.54 ${}^{\circ }$C. However, the current usable depth of 922 m significantly reduces the sustainable continuous heat extraction rate to approximately 50 kW. When a heat pump installation is assumed, a building load of about 65 kW (569.4 MWh per annum) can be supported for 25 years by the system. These load values are likely to meet 58 to 72% of the annual heat demand of the adjacent urban sciences building, which has an annual demand of between 792 and 989 MWh [65], a finding not inferred before. The borehole is best utilised to contribute to a heat network within the Newcastle Helix area because it cannot provide the required heat demand of the complete building. Another option is to consider a much smaller building with less heating demand within the Newcastle Helix site. Parametric studies have indicated that for a DBHE of 922 m depth, the following applies.

• A 50-kW direct DBHE load can be supported by the NSCDGB over a period of 25 years without the fluid inlet temperature going below 3.54 ${}^{\circ }$C. Based on the assumptions made in this work, when used in combination with a heat pump, a building load of 65 kW can be supported for 25 years.
• The effects of heterogeneity appear to be relatively modest and generally within 1 ${}^{\circ }$C for the cases considered when using a DBHE load of 50 kW.
• Depth has a significant effect on the supported load in the long term. More than 100 kW can be extracted for 25 years when considering a DBHE model depth of 1651 m.
• The CXA configuration has a slightly better performance in an extraction-only system in the long term.
• The mass flow rate needs to meet a minimum threshold but a higher mass flow rate will improve the operational period of the system. This should be considered along with the cost required to maintain the chosen mass flow rate. In this study, for the first time, an optimum mass flow rate of around 5 kg/s was determined based on the long-term COP and CSP for the specific coaxial DBHE at the NSCDGB.
• With this modelling study, a predictive hypothesis for the thermal behaviour of a DBHE under heat extraction was generated. This hypothesis has the possibility of being tested subsequently against empirical data generated by a real planned TRT.