Numerical Modeling of a Soil-Borehole Thermal Energy Storage System

Growing concerns about greenhouse gas emissions and fossil fuel consumption have motivated the increased development of renewable energy systems including solar thermal energy harvesting technologies for the heating and cooling of buildings. In recent decades, borehole thermal energy storage (BTES) systems with heat derived from solar technology are rapidly gaining attention and use worldwide (Claesson and Hellstrom 1981; Dalenbäck and Jilar, 1985; Nordell and Hellström, 2000; Sanner and Knoblich, 1999; Sanner et al., 2003; Morofsky, 2007; Sibbitt et al., 2007, 2011, 2012; Wang et al., 2010; Dehkordi and Schincariol, 2014a; Acuna and Palm, 2013; Başer and McCartney, 2015). In BTES systems, a series of U-tube pipes placed in closely spaced (1.5–2.5 m) vertical boreholes are connected to form a closed-loop heat exchanger (Fig. 1). Heat exchange is achieved by circulation of a heat carrier fluid through the closed-loop U-tube pipes. Many BTES systems store heat collected from solar thermal panels in the summer months until it can be extracted for use during the winter months. These thermal energy storage systems present a potentially economical and environmentally sustainable alternative to traditional heating systems because they permit the storage of renewable energy in a space-efficient manner underground.

A major challenge facing BTES systems is their relatively low heat extraction efficiency. Annual efficiency is a measure of a thermal energy storage system’s performance, defined as the ratio of the total energy recovered from the subsurface storage to the total energy injected during a yearly cycle (Dincer and Rosen, 2007). Efficiencies for the first 6 yr of operation of the BTES system at the DLSC range from 6 to 54% (Sibbitt et al., 2011). Nonetheless, the BTES system at DLSC is still able to supply >90% of the thermal energy for heating 52 homes (Sibbitt et al., 2012).

Few studies have investigated the impact of hydrologic parameters on the efficiency of BTES. Common guidelines require a subsurface that is drillable, of high heat capacity, and with low hydraulic conductivity (Schmidt et al., 2003; Pavlov and Olesen, 2012). Several factors may affect BTES system performance, including the volume and geometry of storage space, the number and dimensions of boreholes, the injection–extraction scheme, and the soil and rock mechanical, hydrological, and thermal properties (Ohga and Mikoda, 2001; Dehkordi and Schincariol, 2014b; Başer and McCartney 2015). The influences of these factors are poorly understood, and no specific parameter guidelines have been suggested for enhancing the efficiency of BTES systems. Advancements in the understanding of the specific effects of these parameters are needed because improved efficiency would greatly increase the economic viability of BTES technology.

In the past few decades, numerical modeling has become the standard in the analysis of geothermal energy systems (Breger et al., 1996; O’Sullivan et al., 2001; Diersch et al., 2011); however, models that study the parameters affecting the performance of BTES systems are lacking. Existing models of thermal energy storage systems are limited to aquifer thermal energy storage systems (Molson et al., 1992; Tenma et al., 2003; Lee, 2008, 2010; Lee and Sang, 2008; Sommer et al., 2013). These studies focused on the configuration of borehole design and operational parameters in such systems. Recent studies have been primarily concerned with engineering layout (Lanini et al., 2014; Başer and McCartney, 2015) and modeling methodology (Lazzarotto, 2014). Dehkordi and Schincariol (2014b) examined the sensitivity of a single borehole heat exchanger system to thermal–hydrogeological properties; however, they did not consider the role of such properties in a system that incorporates thermal energy additions to the subsurface. Existing models have not explored the role of soil thermal and hydrological properties such as thermal conductivity, permeability, and degree of saturation on BTES system performance, nor have they looked at the effect of groundwater flow on a multi-borehole system. There exists a critical need for numerical modeling to be applied to these BTES scenarios to better understand heat transport in BTES systems and provide guidance for site selection and system design.

The objective of this study was to examine how soil hydrological and thermal parameters affect the BTES system efficiency. Little is known about the role of site geology in BTES systems despite the subsurface acting as its primary means of thermal storage. Several questions were raised: How do thermal properties of the soil affect BTES system efficiency? How will background groundwater flow affect the BTES system’s ability to store thermal energy? How will the soil permeability affect its ability to store and transfer heat? Will heat transfer processes differ with the degree of saturation in a way that will significantly affect BTES efficiency? To address these questions, a three-dimensional transient coupled fluid and heat transfer model was developed to simulate heat transport in soils for an existing BTES system at DLSC. The model was calibrated using field data from the existing BTES system. The calibrated model was then used as a basis to analyze the BTES system’s sensitivity to soil thermal conductivity, hydraulic gradient, saturated soil permeability, and unsaturated soil permeability.

Drake Landing Solar Community (DLSC), located in Okotoks, AB, Canada, consists of 52 houses, an 800-panel garage-mounted array of solar collectors, short-term thermal storage tanks, and approximately 34,000 m³ of subsurface BTES space (Drake Landing Solar Community, 2015). The BTES system is comprised of 24 series of six connected U-tube pipes that extend from the center to the outer boundary in a radial direction (Fig. 2). During charging, heated fluid is circulated into the boreholes at the center of the array first, then progressively through boreholes toward the boundary. In this case, the center of the array reaches the greatest temperatures. During the winter discharging months, cold water is injected into the outer boundary pipes, is heated by the thermal energy stored in the soil, and then is extracted from the center pipes. The warm water is pumped to the homes and heat is transferred to the homes through forced-air-fan coil units. The system is bounded above by an insulation layer and hydraulic barrier that prevent heat and fluid exchange with the surface (Fig. 2).

No report on soil properties is available for the DLSC; however, a geotechnical analysis was conducted in 2004 at the nearby Twin Valley Resort (G Tech Earth Sciences Corp, 2004). In this study, 15 test boreholes were drilled to study the geology of the subsurface. The holes showed a near-surface layer averaging 4.07 m thick that varied among silty, clayey, and gravely sands. In all cases, this layer overlaid glacial clay till or bedrock. The water table ranged between 2.64 and 5.90 m below the surface. Surface temperatures annually average 5°C but can range from 37°C in the summer to −35°C in the winter. The system has been operational since 12 Apr. 2007. Recorded data are available since 1 July 2007 including: (i) temperature and flow rates of injection and extraction; (ii) temperature for a lateral array of sensors near and below the ground surface between the vertical boreholes; and (iii) temperature at varied depths representing core BTES temperatures near the center of the system.

To examine how soil hydrological and thermal parameters affect BTES system efficiency, a three-dimensional numerical model of the BTES system at DLSC was developed. A three-dimensional transient fluid and heat flow code (TOUGH2; Pruess, 1991), based on an integrated finite difference method, along with the PetraSim graphical interface Version 5 (Thunderhead Engineering) was used. The model was calibrated using the measured temperatures and validated using the calculated efficiencies. The governing mass and energy balance equations for fluid flow and heat transfer in TOUGH2 are described in the supplement. Many of the initial values for the current model were taken from a previous model of the BTES system at DLSC (Zhang et al., 2012). The current model expands on the previous model through the inclusion of a larger data set for temperature comparison, modified soil parameters, and an improved injection temperature scheme. Whereas the previous model utilized an annually repeating charging scheme derived from 1 yr of injection temperature data, the current model incorporates time-varying injection temperatures for 6 yr of operation that reflect the actual injection parameters at the site.

The subsurface storage area with piping has dimensions of 30 m long by 30 m wide by 36 m deep (Fig. 3). By symmetry of the borehole array, one quarter of the storage domain was modeled. Thus, 36 boreholes, or six series each containing six connected U-tubes, were modeled. The borehole region is 15 m long by 15 m wide by 36 m deep. The domain extended an extra 35 m past the outer boreholes. The three-dimensional domain and the imposed boundary conditions are shown in Fig. 3a, with dimensions of 50 m long by 50 m wide by 80 m deep to minimize boundary effects. The model was discretized by 62 varying-space grids in the x (horizontal) and y (horizontal) directions and 26 varying-space grids in the z (vertical) direction. This grid was chosen through systematic mesh refinement such that further grid refinement achieved the same model result. The placement of the U-tube series in relation to the domain is shown in Fig. 3b. All pipes are simulated with square cross-section areas equal to the real dimensions of the circular boreholes. The pipes are anisotropic in their permeability, confining fluid flow to only along the pipe direction while permitting heat conduction in all directions. The connection between boreholes is as shown in Fig. 2b for each series of six U-tubes. There is one center injection–extraction site and one outer injection–extraction site for each series, and thus two points where the series intersects the model surface. The injection site is set to a constant injection temperature and flow rate into the series. Because there is zero permeability between the piping grid cells and the surrounding medium, fluid flow is confined to one direction in the pipe away from the injection site. The fluid flow moves around the “curves” of the pipe due to this anisotropic permeability of the grid cells, thus the “downward” and “upward” U-tubes are all connected and receive the flow from the previous grid cell. An equivalent volume and flow rate of liquid is removed from the model at the extraction site, allowing a proper mass balance. The model contains six of these series of U-tubes and therefore 12 injection–extraction sites.

At the DLSC, no known natural flux boundaries or barriers exist adjacent to or below the BTES system. The following boundary conditions were assigned to represent the in situ conditions. A calibrated temperature of 8.1°C and hydrostatic pressure were imposed at the two far-field boundaries and the bottom boundary. The far-field and bottom boundaries are sufficiently distant from the borehole area so that the boundary conditions do not affect the resulting temperatures. The top boundary was set to the annual average surface temperature of 5.1°C and a pressure of 1.0 × 10⁵ Pa. By symmetry, a zero temperature gradient and zero pressure gradient boundary condition were imposed on the vertical planes that cross the center of the BTES system.

An intrinsic permeability of 1.5 × 10⁻¹⁴ m² is representative for the glacial till (Fetter, 2001) present at the site. The sensitivity of BTES system performance to this value is examined below. Thermal conduction was allowed between the pipe cells and the surrounding medium. Within the pipe cells, we set the thermal conductivity to that of the carrier fluid. Between the pipe cells and the surrounding medium, the thermal conductivity is equal to that of the piping material. A thermal conductivity of 0.51 W m⁻¹ °C⁻¹ was used to represent the polyethylene U-tube piping. Between the piping and ground surface is an insulation layer with a thermal conductivity of 0.23 W m⁻¹ °C⁻¹ (Zhang et al., 2012). The low thermal conductivity of the insulating layer allows some minimal heat loss, as would be expected in the real-life scenario but is sufficient to support a constant-temperature boundary condition for the time and spatial scales considered in our model. For simplicity, this layer extends across the entire top of the model, although at the site this insulation layer only overlaps the BTES system by 5 m. In a real-life scenario, it could be conceived that once the thermal plume extended radially beyond the bounds of the insulated area, greater heat loss to the surface would occur. The relevant properties of the different materials are summarized in Table 1.

The total thermal energy in joules injected during charging stress periods (or total thermal energy in joules discharged during discharging stress periods) during the course of a year (J_t) were divided into each stress period and weighted by the length of the stress period. Although the total energy calculation weighted by the length of the stress period may seem rough, it serves to maintain the total energy balance when incorporated into Eq. [5] because it accounts for the overestimated average charging temperatures and underestimated discharging temperatures in Eq. [4]. In Eq. [6], t_i is the length of the stress period and t_t is the total charging (or discharging) time for the year. Injection enthalpies are calculated by TOUGH2 from the steam table using average injection temperatures calculated in Eq. [4] and atmospheric pressure. Details of the injection scheme are listed in Supplemental Table S1.

The evolution of thermal storage plumes with time is shown in Fig. 5 for the borehole region at x = 6.5 m and z = −17.5 m. Modeled temperatures are shown for September corresponding to the end of the longest charging period and for March corresponding to the end of the longest discharging period. The radius of thermal contours at the depth of 17.5 m are shown in Supplemental Fig. S1.

Several observations can be made from these simulations. The borehole region shows a trend of increasing temperature with time (Fig. 5). By Year 6, average temperatures at all sensor locations were about 20°C higher than in Year 1. In the charging period of the first year (Fig. 5a and 5b), the elevated temperature region was restricted within the borehole region (x < 15 m, z < 40 m), reaching a maximum weighted core temperature of about 52°C. By the end of the discharging period (Fig. 5c and 5d), center temperatures decreased as some of the heat was extracted and some was conducted outward into the surrounding medium. As time elapsed, the elevated temperature region expanded both outward and downward. The maximum simulated weighted core temperature in the sixth year was about 70°C, 18°C higher than that in the first year. Correspondingly, the minimum weighted core temperature in Year 6 was about 48°C, about 19°C warmer than that at the end of the discharging period in the first year.

Measured and simulated temperature profiles are shown in Fig. 6 for several locations including the center and outer boreholes, lateral sensors near the surface, and depth sensors near the center of the system. The simulated temperatures at the center (Fig. 6a) and outer boreholes (Fig. 6b) represent the fluid injection and extraction temperatures. Approximated injection temperatures based on the measured data were imposed at the center boreholes during charging and at the outer boreholes during discharging periods. This allows comparison between the measured and simulated extraction temperatures at the outer boreholes during charging periods and the center boreholes during discharging periods. The simulated annual average extraction temperatures at the center boreholes ranged from 47.6 to 62.6°C. The simulated annual average extraction temperatures at the outer boreholes ranged from 30.0 to 41.2°C. The average simulated temperatures at the center boreholes were slightly above the average measured temperatures, with a maximum difference of 5.7°C. The average simulated temperatures at the outer boreholes were slightly below the average measured temperatures, with a maximum difference of 2.9°C.

Measured soil temperatures from several near-surface locations were compared with the simulated temperatures. Soil temperature profiles for two of these locations, one near the center boreholes (TS-24-2) and one near the outer edge of the borehole region (TS-24-7), are shown in Fig. 6c and 6d. The difference in simulated annual average temperature and measured temperature at the lateral temperature sensors ranged from 0.1 to 3.8°C. The largest discrepancy in each year was at Sensor TS-24-4, which lies in the middle of the borehole region, farthest from the injection–extraction sites. Simulated average annual simulated temperatures at the lateral sensors were slightly higher than those measured for all years.

Temperature sensors are used to continuously measure the subsurface temperature at different depths in the center of the BTES system. The temperatures at two of these depths as a function of time are shown in Fig. 6e and 6f (temperatures at additional depths are detailed in Supplemental Fig. S2 and S3). These sensors represent depths near the top of the borehole region (TS-22-1) and at three quarters of the system’s depth (TS-22-5). Several of the sensors failed before this simulation and were not available for comparison. While temperatures overall increased with time, the oscillating signal of the injection temperatures was dampened with depth. The difference in simulated and measured annual average temperatures at the sensor locations ranged from 0.7 to 4.0°C. Simulated temperatures were slightly higher than measured temperatures near the surface and slightly lower than measured at depth. The simulated temperature time series at depth also shows smaller oscillating amplitude than the measured oscillations, as well as a slight delay in peak temperatures. These discrepancies may be a result of the assumed homogeneity of soil parameters with depth while the actual subsurface is probably heterogeneous. For example, the soil thermal conductivity may differ with geological layering. Although the homogeneous thermal conductivity model appropriately replicated the average injection, extraction, and core temperatures, actual conductivity may vary with depth, thus affecting the distribution of the thermal plume.

The largest discrepancies between the measured and simulated temperatures occurred when there were frequent alternations between charging and discharging. These simulation periods correspond to the time periods when the flow direction changed frequently, sometimes multiple times in a day. This frequent change of flow direction is smoothed out in the model (Fig. 4) and could lead to the discrepancies during these periods.

Overall, the simulated temperatures at the sensor locations closely resemble the daily average measured temperatures at the corresponding locations. This suggests that the method of averaging flow rate and temperature while preserving total energy input properly replicates the system behavior. The model sufficiently approximated heat transport and storage in the BTES system at the DLSC and provides a solid basis for future study of factors affecting BTES system performance.

The annual heat extraction efficiency for the 6 yr of simulation were calculated according to Eq. [1] and plotted in Fig. 7. The simulated efficiencies from Years 1 to 6 are 5.2, 20.6, 35.3, 54.5, 31.1, and 48.7%, respectively. These values closely mirror the rise and fall of the measured efficiencies. The annual heat extraction efficiency generally increases with time, with the exception of Year 5. It is believed that efficiency in Year 5 was lower than the previous year due to the high amount of energy extracted and lower amount of solar thermal energy injected into the system in Year 4. This resulted in a smaller thermal plume and less available energy for extraction at the start of Year 5. It should also be noted that the maximum efficiencies are much higher than those reported by Zhang et al. (2012). This is because their estimates used an annually repeated injection scheme derived from the first year of injection data, whereas this model considered the actual energy inputs from 2007 through 2013.

Heat transport and storage processes in the BTES system can be further dissected through examining the thermal energy distribution in terms of heat injected, heat recovered, heat stored in the borehole region, and heat lost from the borehole region. The heat loss decreased with time as the system became able to more efficiently recover thermal energy, but there is always heat loss into the surrounding porous media because there is heat flow directed away from the borehole region even during the discharging months (Supplemental Fig. S4). After 6 yr, 31.5% of the heat injected into the system was recovered in total, 21.9% remained in the borehole region, and 46.7% escaped the borehole region (Supplemental Fig. S5). After the first year, the cumulative percentage of heat leaving the borehole region remained relatively constant. The fraction of energy stored in the region decreased with time while the energy extracted increased with time, approaching a steady state and reflecting the increase in efficiency.

Because available solar energy and input temperatures vary from year to year, it is difficult to predict the long-term heat extraction efficiency of this particular system; however, general trends in long-term efficiency can be examined using a simplified model with an annual repeating injection cycle. The simplified model can be used to quantitatively examine the controlling factors on long-term BTES system efficiency such as the soil and rock hydrological and thermal properties. Sensitivity analyses are used to understand how different geological settings may influence BTES heat extraction efficiency. Through sensitivity study, the role of soil thermal conductivity, background hydraulic gradient, soil permeability, and soil saturation in BTES systems can be assessed. Most previous models have assumed saturated conditions, but BTES systems are probably located in unsaturated zones. Thus it is important to examine the role of saturation in the thermal energy storage and recovery processes.

The sensitivity of BTES system heat extraction efficiency to soil thermal and hydrological parameters was examined using a simplified model for computational efficiency. The simplified model contained two boreholes and an annually repeating injection scheme. The background soil temperature was set to 15°C. All other parameters were taken from the DLSC model. Two boreholes were simulated with a plane of symmetry next to the center borehole. The domain had dimensions of 100 m long by 50 m wide by 50 m deep. Each U-tube pipe was 30 m deep, with a separation of 5 m between the two U-tubes. The simulated domain was discretized by 38 varying-space grids in the x direction, 43 in the y direction, and 41 in the z direction. The discretization of the domain and the imposed boundary conditions are shown in Fig. 8a. Placement of the U-tube pipes in relation to the domain are shown in Fig. 8b.

A simple annually repeating injection scheme was imposed, with 8 mo of charging at 60°C followed by 4 mo of discharging at 30°C. These temperatures are close to the average injection temperatures at the DLSC site. The length of the charging period was exaggerated to obtain a sizable heat storage plume in the simplified system containing only two boreholes. Thus, the sensitivity results are not intended to portray realistic BTES system numbers but rather to demonstrate how a range in soil hydrologic and thermal properties affects system performance.

A base scenario assumed that the subsurface has homogeneous and isotropic soil similar to the DLSC site and no background groundwater flow. In the sensitivity analyses, the following properties were systematically altered to observe their effect on BTES system performance: soil thermal conductivity, groundwater flow, and soil permeability. All material properties and boundary conditions were maintained with the exception of the variable of interest in each analysis. All sensitivity analyses were conducted under saturated conditions. An additional set of soil permeability sensitivity analyses was conducted under unsaturated conditions as described below. All analyzed variables and their ranges are presented in Table 2. Heat flux into the subsurface was calculated according to Eq. [8]. The average heat flux into the system for the simplified model ranged from 20 to 50 W m⁻² depending on the parameter scenarios.

The baseline scenario for the simplified model used the same parameters from the DLSC model except the new simplified annually repeating injection scheme. The base scenario model was run for 20 yr to observe long-term trends in heat extraction efficiency. The simulation results show efficiency starting at 13.3% for the first year and increasing rapidly in the following few years. The efficiency asymptotically approached a steady-state value of 22.9% by about Year 6 (Fig. 9). The efficiency represents the proportional amount of energy recovered by the system each year. The amount of energy into the system, calculated as in Eq. [3], decreased with time toward a steady state. This can be explained by the decrease in the thermal gradient between the U-tube piping and soil as the thermal storage plume developed. The amount that can be extracted, calculated as in Eq. [2], increased with time toward a steady state. This can similarly be explained by an increase in the thermal gradient in the discharging months as the thermal plume developed and temperatures in the borehole area increased compared with the circulating fluid.

These results suggest that with a given set of base scenario parameters, there exists an upper limit steady-state heat extraction efficiency. The steady-state efficiency of the simplified system was lower than that of the BTES system at the DLSC due to the inclusion of fewer boreholes and a different injection scheme. Efficiencies obtained from sensitivity analyses are compared with the base scenario.

It has been observed that thermal conductivity is influenced by mineralogy, rock density, degree of saturation, and temperature (Farouki 1981; Brandon and Mitchell, 1989; Smits et al., 2013). Heat conduction is more robust at higher thermal conductivities, but the rate for heat conduction away from a reservoir also increases. Thus, it has been suggested that medium thermal conductivities may be desirable for thermal energy storage systems (Socaciu, 2011; Başer and McCartney, 2015).

A range of soil thermal conductivity values from 0.5 to 3 W m⁻¹ °C⁻¹ was considered. The results show that BTES heat extraction efficiency decreases with increasing soil thermal conductivity (Fig. 10). A soil with a low thermal conductivity of 0.5 W m⁻¹ °C⁻¹ showed 47% higher efficiency in the first year than the base scenario, with a thermal conductivity of 2.0 W m⁻¹ °C⁻¹. Energy injection into the soil decreases for lower soil thermal conductivity values, but the ability to extract energy showed slight increases. The development of the thermal plume for low and high thermal conductivity soils is shown in Fig. 11. In both cases, the thermal plume grew outward despite the system being in a heat discharging state. Heat transport away from the piping is always at play as long as there is a thermal gradient between the plume and the surrounding medium to drive conduction. In the case with higher soil thermal conductivity (3.0 W m⁻¹ °C⁻¹), a much larger thermal plume resulted from the increased ability to receive injected energy and conduct the heat into the surrounding soil in the charging phase. This increased ability of the heat to move outward away from the pipes resulted in lower temperatures near the pipes (about 52.5°C) and thus a lower gradient for heat to flow back toward the piping in the discharge phase. At the end of charging in the low-conductivity case, higher temperatures remained near the pipes (about 55.5°C), creating a higher gradient to drive heat flow back toward the piping in the discharging phase. The higher gradient in the low-conductivity case allowed a greater portion of the heat to be recovered, whereas in the high-conductivity case, the thermal plume continued to expand outward.

These results suggest that the thermal properties of the subsurface play a significant role in the heat extraction efficiency of BTES systems. Low thermal conductivity soils allow higher temperature gradients in the discharging phase and thus higher system efficiency. Presently, only large-scale geothermal systems conduct regular thermal response tests to determine the effective thermal conductivity of the subsurface (Stauffer et al., 2013). Proper knowledge of soil thermal conductivity is essential for suitable emplacement and design of BTES systems.

In saturated systems, a background hydraulic gradient and therefore a groundwater flow field may exist. A test of this background groundwater flow influence on BTES heat extraction efficiency was conducted by incorporating a hydraulic gradient of 0.05, corresponding to a groundwater velocity of 3.6 × 10⁻⁸ m s⁻¹. The result showed significant alteration of the thermal plume. A cross-section of the thermal plume after 1 yr is shown in Fig. 12 for the base scenario and the groundwater-flow scenario. The base scenario produced a thermal plume that symmetrically surrounded the piping area. When groundwater flow was incorporated, the thermal plume migrated in the direction of groundwater movement. This relatively small addition of groundwater movement resulted in a decrease in efficiency of 0.6% for the first year of operation. This is because groundwater flow promotes heat advection and carries heat away from the BTES area, whereas in the base scenario, conduction is the primary means of heat transfer. Because greater hydraulic gradients would cause increased advection away from the BTES area, one can infer that BTES system areas of appreciable background hydraulic gradient would have lower system efficiencies. These results indicate that groundwater flow has a limited influence on BTES performance for the given small hydraulic gradient in this case. The significance of groundwater flow influence under a higher hydraulic gradient warrants further study through examining the relative magnitude of heat conduction and heat convection.

When a Rayleigh number exceeds the critical value, the onset of convection in a porous medium occurs. The thermal plume development and heat flow in a cross-section through the piping area for permeability values corresponding to Rayleigh numbers that fall above and below the critical Rayleigh value are shown in Fig. 14. Conduction is the dominant mode of heat transfer for the lower permeability values. The vectors show the flow of heat (W m⁻²) and support the presence of convection at Rayleigh numbers surpassing the critical value. For this system, the Rayleigh number exceeded the critical value, leading to convection, for soil permeability >1.5 × 10⁻¹² m². At Rayleigh numbers slightly higher than the critical value, heat convection starts outside the borehole region (Fig. 14b). At even higher Rayleigh numbers, significant convection occurs outside the borehole region (Fig. 14c). This convective process is accompanied by lateral heat loss and a sharp decrease in system efficiency. A synthesis plot of Rayleigh numbers for differing permeability values is shown in the Supplemental Fig. S6.

After the onset of convection at soil permeability >1.5 × 10⁻¹² m², the amount of energy into the system increases and the amount of energy extracted decreases (Fig. 13). The rate of heat transfer into the soil increases because the convective processes maintain a high temperature gradient away from the pipe. The convective processes effectively carry heat upward and away from the piping, creating an asymmetric thermal plume. The system is unable to extract the heat as efficiently as in cases of little convection.

These model results suggest that low-permeability soils are preferable for enhancing the efficiency of BTES systems in saturated soil layers. Heat losses due to convection away from the piping area in high-permeability soils outweigh the potential benefits of convection occurring within the piping region. Future studies should consider different engineering layouts that may utilize convective processes in a more productive manner. For example, a greater distance between the boreholes may allow increased extraction temperatures associated with higher convection-induced temperature gradients within the borehole array itself. Another potential scenario might be to incorporate vertical subsurface thermal barriers surrounding the storage region to reduce lateral convective heat losses.

Borehole thermal energy storage systems are probably located in unsaturated zones, in part to take advantage of the lower thermal conductivity with degree of saturation (Smits et al., 2013). It has been suggested that heat convection within the piping area may affect the performance of BTES systems in terms of both heat storage and the rate of heat transfer (McCartney et al., 2013). Lu (2001) found that the rate of heat transfer in a convective cell in an unsaturated soil layer may be up to 10 times faster than heat transfer by conduction. Previously, the influence of the soil intrinsic permeability on convective processes within a BTES system in saturated soil was shown. Here, the sensitivity of the BTES model to soil permeability and thermal conductivity in an unsaturated soil was examined.

To model the effect of an unsaturated soil setting, the simplified model was initialized with 25% liquid saturation in the upper 30 m and run to hydrostatic equilibrium before the transient simulation. The relative permeability of liquid water in the unsaturated case varies with the liquid degree of saturation. A van Genuchten–Mualem model (Mualem, 1976; van Genuchten, 1980) for relative permeability and soil water retention was used (details are shown in the supplement, Supplemental Fig. S7 and S8). We chose water retention properties representative of a sandy soil (air-entry capillary pressure of 0.1 Pa), although future studies may look at varying this parameter for different soil types. As the degree of saturation increases, the liquid-phase permeability dominates, while at low saturation the gas-phase permeability dominates.

The unsaturated soil modeling intended to address the following specific questions: How does the heat extraction efficiency of the unsaturated system compare with the saturated case for similar intrinsic permeability values? How does the onset of convection differ in the unsaturated and saturated cases? How does convection affect the BTES system heat extraction efficiency?

As in the saturated soil case, conduction appears to be the dominant mode of heat transfer for the lower permeability values, while convection occurs in higher permeability soils. Figure 16 shows the thermal plume development and heat flow in a cross-section through the piping area for permeability values corresponding to Rayleigh numbers that fall above and below the critical Rayleigh value. The heat flow vectors demonstrate the presence of convectional processes in the unsaturated zone at higher Rayleigh values. Here, the Rayleigh number exceeded the critical value for soils with intrinsic permeability >1.5 × 10⁻⁹ m². Convection contributes to lateral heat loss and correlates to the decrease in efficiency shown in Fig. 15. A synthesis plot of Rayleigh numbers for differing permeability values is shown in Supplemental Fig. S6.

After the onset of convection in soils with intrinsic permeability >1.5 × 10⁻⁹ m², the amount of energy into the system increased and the amount of energy recovered decreased (Fig. 15). As in the saturated case, the rate of heat transfer into the soil increased because the convective processes maintained a high temperature gradient away from the pipe. Convective processes carry heat upward and away from the piping, creating an asymmetric thermal plume and an associated decrease in efficiency.

The unsaturated soil produced higher energy outputs than the saturated case for an intrinsic soil permeability between 1.5 × 10⁻¹² and 1.5 × 10⁻⁹ m² due to the delayed onset of convection in the unsaturated scenario. An unsaturated soil will have a lower effective permeability than its intrinsic permeability because the connection between pores is decreased when less water is available. The low permeability hinders potential convective processes that would otherwise occur in the saturated case. Thus, if the subsurface has a permeability in the range of 1.5 × 10⁻¹² and 1.5 × 10⁻⁹ m², an unsaturated soil is desirable because it will allow heat transfer by conduction rather than convection and the system will run efficiently. Energy output is maximized in low-permeability soils. For soil permeability <1.5 × 10⁻¹² m², either a saturated or unsaturated subsurface will produce similar energy output. Finally, if the soil permeability is >1.5 × 10⁻⁹ m², neither a saturated nor unsaturated subsurface will produce a strong energy output due to the dominance of convectional heat transfer.

The overall heat extraction efficiency is higher in the unsaturated system than that of the saturated system for all permeability values. Specifically, in the unsaturated case, a higher proportion of the injected heat can be recovered. To explain this, one can look to thermal conductivity, which varies with saturation. Lower saturation corresponds to lower thermal conductivity because air makes up a great portion of the void space. The rate of radial heat loss may be lower due to air’s insulating effect. As was shown in the thermal conductivity analysis, a lower thermal conductivity leads to greater energy recovery due to the distribution of the thermal gradients during the charging and discharging phases. The total amount of energy stored in the unsaturated system is also higher than that of the saturated system.

Saturation and permeability were shown to play a large role in BTES system heat extraction efficiency. In both the saturated and unsaturated models, low-permeability soils allowed a higher BTES system efficiency. In high-permeability soils, both saturated and unsaturated systems showed decreased efficiency from convective heat loss. Finally, the overall heat extraction efficiency of unsaturated systems at all permeability values was shown to be higher than for saturated systems due to lower effective thermal conductivity.

A three-dimensional transient fluid flow and heat transfer model for the subsurface BTES system at the DLSC was established using TOUGH2. The model used the time-dependent injection temperatures and energy input measured at the site. A total of 6 yr of annual cycles were simulated from July 2007 through July 2013. The modeled time-dependent temperatures within the borehole region agreed well with the measured values. The model sufficiently replicated the BTES system behavior at the DLSC and provides the basis for future study of factors affecting BTES system performance. The model confirmed a trend of increasing BTES system heat extraction efficiency with time toward a steady state. Heat flow out from the borehole region decreased due to a decreasing thermal gradient, thus more thermal energy was recovered as a result. However, as long as there is a thermal gradient between the thermal storage plume and its surroundings, there will always be some heat loss into the surroundings.

Sensitivity analyses were performed to quantitatively examine how soil hydrological and thermal properties affect BTES system performance. The examined parameters included soil thermal conductivity, hydraulic gradient, and soil permeability under both saturated and unsaturated conditions. Sensitivity analyses led to the following conclusions:

1. The BTES system heat extraction efficiency increases with decreasing soil thermal conductivity. Low-conductivity soils allow more concentrated thermal plumes and thus higher temperature gradients near the pipes during the discharging phase. The higher gradients allow a greater portion of thermal energy to be recovered in the winter months and therefore yield higher system efficiency.

2. The BTES system heat extraction efficiency decreases with the presence of regional groundwater flow. When a hydraulic gradient is incorporated, advective processes can cause the thermal plume to migrate in the direction of groundwater movement. This leads to less thermal energy availability near the piping during the recovery phase and an associated decrease in system efficiency.

3. The BTES system heat extraction efficiency decreases with convective heat losses associated with high soil permeability for both saturated and unsaturated soils. The high permeability allows convective processes to carry heat upward and away from the piping.

4. Saturation was shown to affect the BTES system heat extraction efficiency in two primary ways. First, convection onset occurs at higher intrinsic permeability values in unsaturated soils than saturated soils. This is a product of the influence of water content on the relative permeability of the soil. Lower relative permeability in unsaturated soils hinders the convective processes below a critical threshold intrinsic permeability, leading to higher efficiencies than for saturated soils. Second, the overall heat extraction efficiency of unsaturated systems at all permeability values was shown to be higher than for saturated systems due to the insulating effect of air and a lower effective thermal conductivity.

Insights gained from this model study suggest that site-specific characteristics of the subsurface deserve greater consideration in the planning and design of underground thermal energy storage systems. Based on the sensitivity analyses, a desirable geological setting for BTES systems should consist of low thermal conductivity soils, a low hydraulic gradient, low intrinsic permeability, and low saturation vadose zones. The role of convective processes may be further explored to increase the heat extraction efficiency if they can be confined within the piping region.

We would like to thank Bill Wong of Leidos Canada and Bruce Sibbitt from NRC Canada for sharing the data from DLSC. We also thank Ning Lu and Rong Lei Zhong for initial discussions on TOUGH2 modeling. Funding from NSF Grant CMMI 1230237 is greatly appreciated. The views in this paper are those of the authors alone.