Investigating the Influence of Groundwater Flow and Charge Cycle Duration on Deep Borehole Heat Exchangers for Heat Extraction and Borehole Thermal Energy Storage

Abstract

Decarbonisation of heat is essential to meeting net zero carbon targets; however, fluctuating renewable resources, such as wind or solar, may not meet peak periods of demand. Therefore, methods of underground thermal energy storage can aid in storing heat in low demand periods to be exploited when required. Borehole thermal energy storage (BTES) is an important technology in storing surplus heat and the efficiency of such systems can be strongly influenced by groundwater flow. In this paper, the effect of groundwater flow on a single deep borehole heat exchanger (DBHEs) was modelled using OpenGeoSys (OGS) software to test the impact of varying regional Darcy velocities on the performance of heat extraction and BTES. It is anticipated that infrastructure such as ex-geothermal exploration or oil and gas development wells approaching the end of life could be repurposed. These systems may encounter fluid flow in the subsurface and the impact of this on single well deep BTES has not previously been investigated. Higher groundwater velocities can increase the performance of a DBHE operating to extract heat only for a heating season of 6 months. This is due to the reduced cooling of rocks in proximity to the DBHE as groundwater flow replenishes heat which has been removed from the rock volume around the borehole (this can also be equivalently thought of as “coolth” being transported away from the DBHE in a thermal plume). When testing varying Darcy velocities with other parameters for a DBHE of 920 m length in rock of thermal conductivity 2.55 W/(m·K), it was observed that rocks with larger Darcy velocity (1e-6 m/s) can increase the thermal output by up to 28 kW in comparison to when there is no groundwater flow. In contrast, groundwater flow inhibits single well deep BTES as it depletes the thermal store, reducing storage efficiency by up to 13% in comparison to models with no advective heat transfer in the subsurface. The highest Darcy velocity of 1e-6 m/s was shown to most influence heat extraction and BTES; however, the likelihood of this occurring regionally, and at depth of around or over 1 km is unlikely. This study also tested varying temporal resolutions of charge and cyclicity. Shorter charge periods allow a greater recovery of heat (c. 34% heat injected recovered for 1 month charge, as opposed to <17% for 6 months charge).

1. Introduction

In order to meet net zero carbon emission targets and curb global warming, it is essential to decarbonise heat. Renewable energy can provide alternative sources of heat; however, supply can fluctuate due to seasonality or irregular weather patterns, which often do not match periods of high demand. Therefore, methods of thermal storage must be considered to store energy in low-demand periods (summer) that can be utilised during high-demand periods (winter). While many underground thermal energy storage technologies exist, this paper considers the potential use of BTES in repurposed deep wells, either from ex-geothermal exploration wells or oil and gas wells approaching end of life.

Many have considered the use of DBHEs exploiting energy in closed-loop geothermal systems (e.g., [1,2,3,4,5,6]), with some looking at the potential of repurposed oil and gas wells (e.g., [7,8,9]), or ex-geothermal exploration wells (e.g., [10,11,12,13]). Yet, few have investigated the potential to repurpose DBHEs for BTES. Those that do have investigated systems which rely on conductive heat transfer only in the subsurface, neglecting the influence of groundwater flow (see Figure 1), which can be common in the shallow subsurface or in deeper reservoirs. In shallow to medium depth BTES arrays, this has been proven to be an important factor influencing storage efficiency (e.g., [14,15]). In particular for the UK, it has been shown that there is a strong potential for groundwater flow to occur in deep Permo-Triassic aquifers [16,17].

Xie et al. [18] tested 5 parameters for a single well BTES system in a repurposed oil and gas well. They concluded the most influential parameters, in order, were inlet temperature during extraction, inlet temperature during charge, flow rate and geothermal gradient. Further work from Brown et al. [11] built on this idea, but for the repurposing of an ex-geothermal exploration well, testing 10 parameters and defining a new metric for storage efficiency. They concluded that the most influential parameters based on rankings of the global sensitivity analysis were inlet temperature during charge, inlet temperature during extraction, flow rate and geothermal gradient. When using the new metric defined in the study by Brown et al. [11], the authors demonstrated that deep borehole heat exchangers struggle to achieve a storage efficiency of greater than 17% (where storage efficiency is defined as increase in output from a DBHE as a fraction of heat injected) for 6 month charge/discharge cycles. In light of this, DBHEs appear to be not optimum candidates for BTES, unless there are no other viable uses for surplus heat and the construction of other forms of thermal energy storage are economically unviable. Both the aforementioned studies of Xie et al. [18] and Brown et al. [11] provide useful analysis when determining the potential of such systems; however, they neglect the influence of groundwater flow and do not consider the longevity or varying periods of charge on such systems. Therefore, this study aimed to bridge this gap by considering the most influential parameters (flow rate, inlet temperature during extraction and inlet temperature during charge) with varying levels of groundwater flow. Long term simulations were then undertaken to consider the efficiency of BTES over several cycles of operation of a single well DBHE.

In this study, the Newcastle Science Central Deep Geothermal Borehole (NSCDGB) was considered as a case study as plans are in place to repurpose this for testing as a pilot DBHE in the UK (study area highlighted in Figure 2). It is a former geothermal exploration well that encountered low-permeability strata within the Carboniferous Fell Sandstone Formation at depths in excess of 1400 m [19]. It was initially identified as a geothermal prospect due to the geothermal heat flux and elevated geothermal gradient associated with the deeply concealed radiogenic North Pennine Batholith [20,21]. Groundwater conditions in the target reservoir (the Fell Sandstone) were not extensively tested, although very low permeability (~8.1 × 10⁻¹⁰ m s⁻¹) was inferred from the lack of appreciable fluid yield from the borehole [19]. Nevertheless, groundwater was identified during well drilling at higher stratigraphic levels in the well [19]. Therefore, the influence of groundwater movement on extraction and BTES was tested in this study using a modelling approach. OGS software uses the finite-element method for spatial discretisation and has been verified/validated against a series of solutions for shallow and DBHE examples (e.g., [10,12,22,23,24,25]). Further benchmarking was also undertaken in this study in comparison to T2Well-EOS1/TOUGH2. The ‘Dual Continuum’ method is used to model the wellbore, with components of the DBHE (grout, fluid, casing) treated as a 1D line source, integrated within a 3D medium for the subsurface rocks.

This work has further universal implications for DBHEs operating to extract heat only, as previous work has suggested the influence of groundwater flow is minimal when thin aquifers are present [23]; however, this work investigates the presence of large-scale groundwater flow over the length of the DBHE. Groundwater flow was modelled as a regional Darcy velocity (groundwater flow per unit perpendicular cross-sectional area of aquifer), which is defined as the product of hydraulic conductivity and hydraulic gradient. Locally to the UK, low-enthalpy Mesozoic Basins, such as the Cheshire and Wessex Basins, have large thicknesses of sandstones up to 2 km thick (Knutsford-1 well in the Cheshire Basin—see [26]) and a high geothermal potential (e.g., [27,28,29,30,31]). Future geothermal exploration wells or oil and gas wells approaching end of life (e.g., [32,33]) may be considered for repurposing as a DBHE for heat extraction or BTES. These will penetrate large quantities of aquifers with advective regimes which must be considered. Results from this study will therefore be beneficial in any decision making in repurposing of wells in areas of thick aquifer concentrations.

To summarize, the objectives of this paper were (i) to model regional groundwater flow for deep borehole heat exchangers for both heat extraction and BTES, (ii) to test the impact of charge-discharge cycle duration on storage efficiency and (iii) to understand if operational duration has an impact on the performance of single well deep BTES.

2. Methods

2.1. Governing Equations

OGS software was used to simulate BTES in a single DBHE. The mesh was spatially discretised using the finite-element method which incorporates a ‘dual-continuum’ approach which treats the DBHE as a 1D line source and the surrounding rock as a 3D mesh (Figure 3). In this study, the model used a CXA configuration (coaxial pipe with annular inlet) (Figure 1) as previous work shows the influence of flow direction (i.e., inlet through the annulus versus the central pipe) has minimal impact on DBHEs [11]. In the CXA configuration, fluid is circulated down the annular space, exchanging heat or coolth with the subsurface via conduction through the borehole wall, before being circulated back to the surface through the central pipe. Groundwater flow in the subsurface was modelled using a constant regional Darcy velocity ($v_f$) (Equation (1)) (see [34] for methods of modelling groundwater flow), which is a common approach in similar studies for shallow arrays (e.g., [35,36]). It is worth noting that this is unlikely to replicate the true subsurface hydrogeological conditions at the NSCDGB, but will help to improve the understanding of the impact of groundwater flow for DBHEs and deep BTES.

The DBHE was modelled with a simplified single grout layer and casing. Four governing equations model heat transfer within (1) the rock formation, (2) the grout, (3) the borehole casing and (4) the central coaxial pipe. The governing equation in the rock is given by the energy balance (e.g., [22,23]):

$$\frac{\partial }{\partial t}\left[\varphi {\rho }_f{c}_f+\left(1-\varphi \right){\rho }_r{c}_r\right]T_r+\nabla \cdot \left({\rho }_f{c}_f{v}_f{T}_r\right)-\nabla \cdot \left({\Lambda }_r\cdot \nabla T_r\right)=H_r$$

(1)

where $\varphi$ is the effective porosity, ${\rho }_r$ is the matrix density of the rock, ${\rho }_f$ is the density of the circulating fluid, $T_r$ is the temperature of the saturated rock, $c_r$ is the matrix specific heat capacity of the rock and $c_f$ is the specific heat capacity of the fluid. $H_r$ is the source term and ${\Lambda }_r$ is the thermal hydrodynamic dispersion tensor, which depends on the bulk thermal conductivity of the rock ${\lambda }_r$. A heat flux ($q_{n{T}_r}$) boundary condition was adopted between the DBHE and surrounding rock:

$$q_{n{T}_r}=-({\Lambda }_r\cdot \nabla T_r)$$

(2)

Conductive heat transfer dominates the thermal regime within the grout:

$$(1-{\varphi }_g){\rho }_g{c}_g\frac{\partial T_g}{\partial t}-\nabla \cdot \left[(1-{\varphi }_g){\lambda }_g\cdot \nabla T_g\right]=H_g$$

(3)

where the subscript $g$ represents the grout. In the annular space (i.e., inlet, subscript $i$) and central coaxial pipe (i.e., outlet, subscript $o$), heat transfer is governed by advection:

$${\rho }_f{c}_f\frac{\partial T_i}{\partial t}+{\rho }_f{c}_f{\mathit{v}}_i\cdot \nabla T_i-\nabla \cdot ({\Lambda }_f\cdot T_i)=H_i$$

(4)

$${\rho }_f{c}_f\frac{\partial T_o}{\partial t}+{\rho }_f{c}_f{\mathit{v}}_o\cdot \nabla T_o-\nabla \cdot ({\Lambda }_f\cdot T_o)=H_o$$

(5)

where ${\mathit{v}}_i$ and ${\mathit{v}}_o$ are the inlet and outlet fluid velocity vectors, respectively. ${\Lambda }_f$ represents the hydrodynamic thermo-dispersion tensor which in this example can be simplified to equal the fluid thermal conductivity (${\lambda }_f$).

The horizontal thermal resistance to heat flow within the DBHE was modelled as an analogue to a resistor network (e.g., see [37]). Thermal resistance for heat flow was modelled between (a) the rock and grout ($R_{gr}$), (b) the grout and borehole casing ($R_{fig}$) and (c) the borehole casing and the central coaxial pipe ($R_{ff}$). Using the outer surface area at the relevant interface, the thermal resistances can be expressed as heat transfer coefficients ($\Phi$) which appear in the boundary conditions for the grout, borehole casing, and central coaxial pipe. The boundary condition for Equation (3) can be expressed as:

$$q_{n{T}_g}=-{\Phi }_{gr}(T_r-T_g)-{\Phi }_{fig}(T_i-T_g)$$

(6)

In a similar procedure, the boundary conditions for Equations (4) and (5) can be expressed respectively as:

$$q_{n{T}_i}=-{\Phi }_{fig}(T_r-T_i)-{\Phi }_{ff}(T_o-T_i)$$

(7)

$$q_{n{T}_o}=-{\Phi }_{ff}(T_i-T_o)$$

(8)

The heat transfer coefficients in Equations (6)–(8) are a function of the borehole casing diameter ($d_{ca\sin g}$), the central coaxial pipe diameter ($d_{central}$), and the borehole diameter ($D_b$). The heat transfer coefficients can be given by: ${\Phi }_{gr}={1/R_{gr}\pi D}_b$, ${\Phi }_{fig}=1/R_{fig}\pi d_{ca\sin g}$, and ${\Phi }_{ff}=1/R_{ff}\pi d_{central}$. See, for example [38], on how to compute the thermal resistances.

2.2. Model Set Up, Initial Conditions and Parameterisation

Under initial conditions, a linear geothermal gradient (33.4 °C/km) was established, increasing with depth. This is in line with those values recorded at the wellsite and modelled in the literature (e.g., see [19,20,21,30,39,40]). The surface temperature was set at a fixed constant Dirchlet boundary corresponding to 9 °C and the remaining lateral and basal boundaries were fixed as Neumann no-flow boundaries with heat flux set equal to zero. Under initialisation of the simulation, all wellbore components (grout, pipe, fluid) were set to equal initial conditions (i.e., increase with depth under the natural geothermal gradient). Lateral and basal boundaries of the model were set to extend to a distance where thermal interaction with the DBHE would be minimal (and not occur). For the setup of groundwater flow, a constant Darcy velocity was established in the x direction throughout the model. Finally, the domain size was set at 500 × 500 × 1300 m (x, y, z) (Figure 3).

The base case model parameters are listed in

Table 1. Thermo-physical parameters of the model. Model parameters are either taken from the literature, assumed unpublished values (compiled by [42,43]), calculated values or given as the most likely value. Note the inner pipe is the coaxial pipe and the outer pipe is the casing. The real nature of the casing situation is notably more complex than that modelled. The thermal properties of the rock in the subsurface are taken as the weighted average from Kolo et al. [10].

Parameter	Value	Units	Symbol
Borehole Depth [19]	920	m	$L$
Borehole Diameter [19]	0.216	m	$D_b$
Outer Diameter of Inner Pipe	0.1005	m	-
Thickness of Inner Pipe	0.00688	m	-
Thickness of Outer Pipe	0.0081	m	-
Thickness of Grout	0.01905	m	-
Thermal Conductivity of Polyethylene Inner Pipe	0.45	W/(m·K)	-
Thermal Conductivity of Steel Outer Pipe	52.7	W/(m·K)	-
Density of Rock [10,44]	2480	kg/m3	${\rho }_r$
Thermal Conductivity of Rock [10,19,45,46,47]	2.55	W/(mK)	${\lambda }_r$
Specific Heat Capacity of Rock [10,48,49]	950	J/(kg·K)	$C_r$
Volumetric heat capacity of rock	2.356	MJ/(m3·K)	-
Density of Grout	995	kg/m3	${\rho }_g$
Thermal Conductivity of Grout	1.05	W/(m·K)	${\lambda }_g$
Specific Heat Capacity of Grout	1200	J/kgK	$C_g$
Density of Fluid [1]	998	kg/m3	${\rho }_f$
Thermal Conductivity of Fluid	0.59	W/(m·K)	${\lambda }_f$
Specific Heat Capacity of Fluid	4179	J/kg·K	$C_f$
Surface Temperature [45]	9	°C	-
Geothermal Gradient [19,45]	33.4	°C/km	-
Porosity	20	%	$\varphi$
Volumetric Flow Rate	0.005	m3/s	$Q$

and have been collected from in-situ data or literature. The geology was modelled as a homogenous rock and the subsurface parameters (i.e., density, specific heat capacity and thermal conductivity) listed are the bulk saturated properties for the rock taken as a weighted average through the vertical dimension of the model [10]. The initial inlet temperature during charge was fixed at 95 °C, whilst during extraction it was fixed at 5 °C. An initial extraction and charge period of 6 months was used, which corresponds to the typical heating season in the UK.

Further parameters tested throughout the study are highlighted in

Table 2. Parameters tested in this study.

Parameter	Minimum	Maximum	Units
Groundwater Velocity (Darcy velocity)	None (conduction only)	1e-6	m/s
Inlet Temperature (charge)	65	95	°C
Inlet Temperature (extraction)	5	20	°C
Internal Fluid Flow Rate	1	7	L/s

. The groundwater flow velocity was considered typical in shallow aquifers and is consistent with aquifer values reported in the literature (e.g., [35,41]). The regional groundwater flux was applied to the Darcy velocity in Equation (1). For context, a Darcy velocity of 1e-7 m/s results from a hydraulic conductivity of 1e-5 m/s (or 0.9 m/d, corresponding to a reasonably clean poorly lithified sand or sandstone aquifer) and a hydraulic gradient of 1% [36]. It is worth noting that the higher Darcy velocities modelled are more likely to be encountered in shallower settings rather than deep, but were considered in this study as a maximum upper limit of groundwater flow for more universal applications.

2.3. Evaluation Metrics

The (time variant) thermal power ($P$) at any given time during both charge and extraction can be calculated as (e.g., [50,51]):

$$P={\rho }_f{c}_{f}Q\left(T_{out}-T_{in}\right)$$

(9)

where $Q$ is the volumetric flow rate (e.g., Table 1). The total thermal energy injected or extracted during a charge or discharge period can be calculated by integrating the power over the time of the charge or discharge cycle to obtain a value in TJ (10¹² joules). The storage efficiency or recovery of heat is essential to understanding the performance of BTES. While different methods are used for calculating this value, in this study the method by Brown et al. [11] was used. Storage efficiency ($S{E}_{new})$ was defined as the difference in energy extracted with and without charge with respect to the total energy injected:

$$S{E}_{new}=\frac{Total energy extracted post charge-Total energy extracted no charge}{Total energy injected during charge}$$

(10)

Therefore, the increase in energy recovered after charge can be determined. For this study, the constant inlet temperature boundary condition was applied for modelling purposes, but in realistic operational conditions a constant thermal power would be used to reflect the surface supply of heat for charge and consumer demand for heat.

2.4. Benchmarking

OGS has been tested against a range of solutions generally providing little discrepancy; these include an analytical solution by Beier [52] (see [4,12,23]) and a numerical solution by Brown et al. [1] (see [10,11]), among many others. At present, data is unavailable for the NSCDGB site. Therefore, further benchmarking was undertaken in this study against a solution in T2Well-EOS1/TOUGH2 (see [53] for software description or [3] for DBHE setup). T2Well-EOS1/TOUGH2 is a wellbore-reservoir simulator that uses the integrated finite-volume method for the surrounding rocks and the drift flux method for solving fluid flow within the wellbore. It also models density, specific heat capacity and thermal conductivity of fluid in the wellbore as a function of pressure and temperature, which is in contrast to OGS. For T2Well-EOS1/TOUGH2, the pressure was assumed to increase with depth from 1 atm at the surface to a maximum pressure of 9 MPa at the base of the model. It also uses a 2D axis-symmetrical layout with the spatial domain set to 1423 and 1500 m (x, y). In this scenario, the parameters defined in Section 2.2 (Table 1) were used for comparison between the two software with no groundwater flow to test the DBHE model (i.e., conductive heat transfer in the surrounding rocks only). Other boundary and initial conditions were set up identical to those outlined in Section 2.2, other than for the base of the model where an infinite heat source term boundary was adopted.

It can be observed in Figure 4 that the models have a good fit with time (Figure 4a), with minimal variations in outlet temperature during extraction and charge. At the end of the charge period the difference in outlet temperature was ~0.6 °C and at the end of the extraction period it was ~0.2 °C. Similarly, there were minor differences in temperature corresponding to depth recorded at the end of the charge and extraction periods, respectively (Figure 4b,c). This could be due to either imposing an initial pressure boundary condition or the non-constant thermo-physical properties of the working fluid on T2Well-EOS1/TOUGH2, which will influence pressure and temperature with time within the DBHE.

3. Results

In this section the results of the numerical models are provided. Initially, the results without charge were discussed as they were required to evaluate the storage efficiency highlighted in Section 2.3 (Equation (10)). Then the impact of groundwater flow on varying thermal energy storage scenarios was evaluated. Four charge-discharge cycle types were investigated: (1) Continuous extraction only for a period of 6 months, (2) a baseline annual case of 6 months charge, followed by 6 months extraction, (3) an annual case of 1 month charge followed by 11 months extraction, (4) an annual simulation of 3 months charge followed by 9 months extraction and (5) an extended simulation of case 2 but for 5 years.

3.1. Extraction Only

Due to the method of calculating storage efficiency for BTES (i.e., Equation (10)), a period of “extraction only” must be evaluated for each case of thermal energy storage. As highlighted in Figure 5, the influence of groundwater has minimal impact on DBHEs when the Darcy velocity is less than or equal to 1e-7 m/s. When the Darcy velocity significantly exceeds this value (i.e., for 1e-6 m/s), it can be observed that there is a reduction in cooling of the borehole due to asymmetrical thermal drawdown and the cool thermal plume being transported away from the DBHE (Figure 5c,d). At values less than this, the advective transport of heat away from the DBHE by groundwater is minor. With a Darcy velocity of 1e-6 m/s, a thermal plume extends downstream from the DBHE, reaching ~40 m after 6 months, with a width of <10 m. In the absence of groundwater flow, the conductive radius of thermal influence extends to ~9 m from the borehole after 6 months.

At the end of the 6-month period of heat extraction, for the varying groundwater flows the minimum outlet temperature was 7.84 °C (no groundwater flow, conduction only) and the maximum was 9.04 °C (with 1e-6 m/s Darcy velocity). With no groundwater flow, the thermal power extracted was 66.4 kW (average over the 6 months) and 59.1 kW as recorded at the end of the 6 months. When groundwater Darcy velocity was 1e-6 m/s, the thermal power extracted was 85.8 kW (average) and 84.3 kW recorded at the end of the 6 months. This represents an increase of 25.2 kW in final thermal output. A variety of constant inlet temperatures (Figure 6) (from 5 °C to 10 °C, 15 °C and 20 °C) were investigated; it was found that outlet temperatures increase and thermal output decreases approximately linearly with increasing inlet temperature. The effect of flow rate within the borehole was analysed, and it was found that outlet temperature decreases and thermal power output increases (both non-linearly) with increasing flow rate (Figure 7). Increasing the Darcy velocity positively impacts both parameters by removing the coolth in proximity to the borehole, but the effect only becomes significant above 1e-7 m/s. The maximum increase in thermal power was observed for the highest circulation flow rate within the DBHE of 7 L/s, where thermal power recorded at the end of the 6 month simulation increased by 27.94 kW for the case with ($v_f$= 1e-6 m/s) relative to the case without groundwater flow (Figure 7b). Similarly, for the lowest inlet temperature during extraction (5 °C) the thermal power increased by 25.1 kW for with ($v_f$= 1e-6 m/s) and without groundwater flow.

3.2. Borehole Thermal Energy Storage

3.2.1. Influence of Inlet Temperature during Charge

The effects of DBHE inlet temperature during charge cycles was investigated for values between 65 and 95 °C as these might be typical inlet temperature for BTES systems around the world (e.g., [54]). An increase in charge inlet temperature corresponds to more thermal energy stored in the subsurface across all scenarios (Figure 8a). For lower Darcy velocities (1e-7 m/s or less), this corresponds to more energy recovered from the subsurface during the extraction period (Figure 8b). The increase in energy extracted as charge inlet temperature increases is relatively small when compared to the case without charge (Figure 8c); therefore, the storage efficiency reduces with greater charge inlet temperatures (Figure 8d).

When the Darcy velocity in the subsurface was set to 1e-6 m/s, far more energy was injected during charge for all temperatures as compared with lower Darcy velocities (Figure 8). This can be observed by the reduced outlet temperature during all charge periods (Figure 9); however, the outlet temperature during extraction was almost identical to the other lower velocity scenarios. This corresponds to minimal changes in the energy recovered after charge, regardless of the inlet temperature (Figure 8). The increase in energy stored was due to heat being transported away from the DBHE in a thermal plume by the movement of groundwater which was not recoverable. This aids in maintaining a larger difference in temperature between the DBHE and surrounding rock such that higher heat transfer occurs.

Inlet temperature during charge had a minor negative correlation to the storage efficiency and recovery of heat (Figure 8d). This was due to more heat stored in the formation which cannot be extracted efficiently, particularly with increased groundwater flow. The highest storage efficiency was recorded as 16% when there was no groundwater flow and charge inlet temperature was set at 65 °C, while the poorest storage efficiency was recorded at ~1% for the highest inlet temperatures (85 and 95 °C) and maximum Darcy velocity imposed.

3.2.2. Influence of Inlet Temperature during Extraction

Inlet extraction temperatures of 5 to 20 °C were investigated, all of which are lower than the average undisturbed ground temperature along the length of the DBHE (i.e., 24.36 °C). If the value exceeds this average, then there is a risk that heat will tend to enter the geological environment, rather than being extracted from it. Furthermore, the lowest temperature was chosen to ensure no freezing within the DBHE (or potential heat pump) (e.g., [23]); however, if an antifreeze solution was used within the DBHE then lower temperatures could be used for the inlet. It was observed that the total energy recovered during extraction, with and without charge, significantly increases approximately linearly with declining inlet temperatures (Figure 10b,c). The impact of varying the extraction temperature on storage efficiency was negligible, as highlighted in Figure 10d. On the one hand, therefore, operating with a lower input temperature maximises energy extraction, but it has little effect on storage efficiency and lower inlet and outlet temperatures would, in turn, mean a lower operational efficiency for any heat pump being employed for heat delivery to a consumer.

3.2.3. Influence of the Deep Borehole Heat Exchanger Internal Flow Rate

Flow rates for the water circulating internally within the co-axial DBHE of 1 to 7 L/s were investigated. With increasing flow rates, more energy was stored and extracted; this is because increasing fluid flow rate reduces the temperature difference between inlet and outlet, which in turn raises or lowers the average internal DBHE fluid temperature (during charge and extraction, respectively) and increases the temperature differential between the DBHE and the formation. Low groundwater velocities (<1e-8 m/s) result in minimal impact during extraction or charge. When the groundwater flow is 1e-7 m/,s there is a minor impact on the recovery of heat (Figure 11b), which consequently reduces the storage efficiency across all flow rates. For flow rates of 1 L/s in the wellbore the storage efficiency was reduced in comparison to the no groundwater flow simulation by 0.5%, whilst for 7 L/s it 1% lower.

When the Darcy velocity was set at 1e-6 m/s, flow rate had minimal impact on the recovery of heat (Figure 8b) in contrast to heat extraction without storage (Figure 8c). Furthermore, storage was very inefficient for higher Darcy velocities and most of the energy injected was transported away from the DBHE and the recovery of heat was always less than 3% for different fluid flow rates in the wellbore.

3.2.4. Influence of Varying Charge Periods

The charge period was reduced to (1) 3 months, followed by a discharge period of 9 months, and (2) 1 month, followed by a discharge period of 11 months. This enabled the testing of different charge periods to evaluate if the storage efficiency improves with varying modes of operation. With shorter charge periods there is less time for heat to propagate away, increasing the storage efficiency across all groundwater velocities (Figure 12). When the groundwater velocity was less than or equal to 1e-7 m/s the storage efficiency was always greater than 21 and 32% for the 3 months and 1 month charge periods, respectively (Figure 12d). Maximum storage efficiencies were recorded for the conductive heat flux only scenario with storage efficiencies reaching 23 and 34%, respectively, for 3- and 1-month charge periods.

This suggests that DBHEs are most effectively used for underground thermal energy storage, when used with short periods of intense (high power) heat charge, followed by longer periods of discharge. This is because, during the short charge period, heat does not have time to disseminate far from the DBHE, before discharge commences, the radial thermal gradient is reversed and the heat is induced back towards the DBHE.

3.2.5. Long Term Simulations

The simulation has been extended to cover a 5-year period, comprising annual charge and discharge cycles each of 6 months duration. Inlet temperatures during charge have been constant at 95 °C and during discharge at 5 °C, while the internal DBHE fluid flow rate was 5 L/s. In this study, the recovery of heat improves slightly with time, with outlet temperature increasing from the end of the first to last (year 5) charge periods from 8.4 to 8.8 °C for the no groundwater flow case (Figure 13). This is simply because, for each year, more heat is injected (~3.6 TJ per cycle) than is extracted (~1.6 TJ per cycle), so the ground has an inevitable tendency to warm. For the highest groundwater Darcy velocity of 1e-6 m/s, the outlet temperature at the end of the extraction period remains constant with time at 9.04 °C. This emphasises the observation that the importance of groundwater flow for enhancing heat extraction from a DBHE (irrespective of any charge cycle) is more than any thermal storage effect (which declines with increasing ground water flow). For the highest groundwater Darcy flow of 1e-6 m/s, the length of the corresponding thermal plume increases with time, in excess of 250 m after 5 years as the plume reaches the outer boundaries of the model (i.e., in comparison for the first cycle v last in Figure 14).

4. Discussion

4.1. Implications of an Active Groundwater Flow on Heat Extraction Only

Groundwater flow has a positive impact when extracting heat from DBHEs, compared to cases where there is no advective transport of heat in the ground. The advective transport of heat through the DBHE (equivalent to the advective transport of “coolth” away from it) only becomes significant when the groundwater Darcy velocity exceeds 1e-7 m/s (Figure 5c). Such a Darcy velocity is only achieved when hydraulic conductivity exceeds 1e-5 m/s (or 0.9 m/d) for a hydraulic gradient of 0.01. Such hydraulic conductivities and gradients are the exception rather than the rule, and unlikely in most systems (including the NSCDGB), at depths of 1–2 km, explaining Chen et al.’s [23] conclusion that groundwater only has a minor impact on DBHEs. However, this study concludes that, in cases where DBHEs penetrate thick aquifers with a Darcy velocity in the region of 1e-6 m/s or greater, available thermal outputs will be positively impacted by groundwater flow.

4.2. Implications of an Active Groundwater Flow on Borehole Thermal Energy Storage

In contrast to single DBHEs that operate to extract heat only, the storage of thermal energy in DBHEs is negatively impacted by groundwater flow. Groundwater flows with a Darcy velocity of 1e-6 m/s reduce the thermal energy storage efficiency to less than 5% across all simulations for a period of 6 months injection and 6 months extraction. Previous work has suggested that deep BTES in a single borehole is less efficient than more equidimensional (and shallower) BTES arrays, but nevertheless could be useful if there is an otherwise unusable surplus of heat in proximity to an abandoned well that could be repurposed [11]. This study is in agreement that the recovery of heat after storage is limited due to the poor efficiencies. Any potential for storage could be limited even further if groundwater is present. It does, however, appear the optimal form of underground thermal energy storage in DBHEs is with short term, high-power charges of thermal energy, followed by longer periods of heat extraction. Cycles of 1 month charge, followed by 11 months extraction appear to result in storage efficiencies in excess of 30%. Therefore, there could be potential to combine this technology with curtailed wind or waste heat discharged in shorter time periods, rather than for seasonal thermal energy storage. Future work should investigate the influence of groundwater flow in shallow aquifers or fractures within the ground which are more likely to be encountered than the significant thicknesses of aquifer modelled in this study.

4.3. Implications to Prospective Areas within the UK

In the UK, and internationally, there are a series of deep sedimentary basins with thick successions of sandstone aquifers. These are usually targets for both hydrocarbon and geothermal developments and they can encounter deep groundwater (e.g., [16,17]). Within England, the four major Mesozoic basins include the Wessex, Cheshire, East England, and Worcester; their infill is dominated by Permo-Triassic aquifers that can be up to 2 km in thickness (e.g., Knutsford-1 well in the Cheshire Basin—see [26]). Therefore, if looking to develop or repurpose a DBHE then local groundwater conditions must be considered. For context, the regional Darcy velocity value of 1e-7 m/s can be related to sandstone aquifers with a hydraulic conductivity of ~1e-5 m/s and hydraulic gradient of 1% [36]. The geometric mean of core data from Permo-Triassic rocks in the UK is within this magnitude (albeit from shallower boreholes) [55] and DBHEs could be impacted by deeper flows.

4.4. Comparison with Previous Studies for Shallow BTES

Modelling results indicated regional groundwater flow impacts a deep single BHE BTES system when Darcy velocity approaches 1e-6 m/s. This is similar to that modelled in the literature for shallow BTES arrays which typically indicate that a Darcy velocity less than c. 1e-7 m/s has little impact on borehole heat exchanger performance [35,36,56,57,58,59,60]. Therefore, the influence of groundwater is not dependent on the depth of the BTES system in the subsurface.

5. Conclusions

In this study, the influence of groundwater flows on DBHEs for both heat extraction and BTES was considered. While groundwater flows positively improve the amount of heat that can be extracted it limits the possibility of deep single well BTES. A model on OGS was constructed to test the local advective heat transport associated groundwater movement on the operation of a DBHE. Initial, preliminary testing in comparison to T2Well-EOS1/TOUGH2 subsurface showed the model to have limited discrepancy (outlet temperature within 0.6 °C), giving confidence in the results of this study. A series of simulations to test the impact of the most influential parameters on BTES (as determined by [11]) and groundwater flow were then tested for heat extraction for 6 months and thermal energy storage for 1 year periods. The influence of the charge time period was also considered, before testing cyclic operation for a number of years. The key conclusions were:

• Groundwater flow from thick aquifers with a Darcy velocity around or greater than 1e-6 m/s has a positive impact on heat extraction using DBHEs and will likely increase the longevity of such systems. The impact of this reduces with increased extraction inlet temperature whilst it significantly improves the achievable thermal power with increased flow rates.
• In contrast, increasing groundwater flow (approaching or above 1e-6 m/s) for BTES in single well DBHEs negatively impacts the storage efficiency (<5%).
• Increasing the internal DBHE fluid flow rate in lower Darcy velocity conditions did improve the performance for BTES by over 5%.
• Reducing the charge period significantly increases the recovery of heat, with charge periods of 1 and 3 months (followed by 11 months and 9 months discharge) resulting in storage efficiencies of up to 34 and 23%, respectively. Therefore, it may be more beneficial for DBHEs used for thermal energy storage to apply short, intense charge periods, followed by longer discharge periods.
• Simulation over a longer (5 year) series of charge-discharge cycles only has a minor impact on the recovery of heat, at least in the “fixed inlet temperature” mode of simulation that has been adopted in this paper.
• To maximize the storage efficiencies in single well BTES systems, specific to the modelled parameters in this paper at 920 m depth, it appears that it is best to have lower charge temperatures (of 65 °C), higher circulation flow rates (of 7 L/s), lower charge periods (of 1 month of less) and target subsurface systems with aquifer Darcy velocity of 1e-7 m/s or less.