This is an open access article distributed under the CC BY-NC-ND license (


In this study, we numerically and experimentally evaluated heat transfer in soils under unsaturated conditions in the context of simulating a laboratory-scale, three-dimensional soil-borehole thermal energy storage (SBTES) system. Most previous studies assumed that soil thermal and hydraulic properties are constant and that heat transfer in soil occurs only in the form of conduction, neglecting convective and latent heat transfer. In addition, there is a lack of data from controlled experiments to validate multiphase numerical models that can be used to better study SBTES systems installed in the vadose zone. The goal of this study was to evaluate the significance and impact of variable soil thermal and hydraulic properties, as well as different heat transfer mechanisms, in unsaturated soils. Four laboratory experiments were performed using a three-dimensional laboratory-scale SBTES system, incorporating sensors to collect soil temperature and moisture data at high spatial and temporal resolutions. Experimental data were then used to validate a numerical model that solves for water and vapor flow and considers nonequilibrium phase change. Results revealed that for the test conditions studied, convective heat transfer was higher than conductive heat transfer in the middle of the borehole array. Moreover, for the experiments on unsaturated sand, about 10% of the total heat transfer was in the form of latent heat. Simulation results demonstrated the importance of including both convection and latent heat in SBTES system modeling. Results also revealed a need for using saturation-dependent effective thermal conductivity in modeling SBTES systems in unsaturated soils rather than using constant values such as those obtained from system thermal response tests.

Soil-borehole thermal energy storage (SBTES) systems are used to store heat generated from renewable resources (e.g., solar energy) in the subsurface for later extraction and use in the heating of buildings (Sibbitt et al., 2007; Pinel et al., 2011; McCartney et al., 2013; Başer et al., 2015; Catolico et al., 2016). Seasonal storage of thermal energy in geothermal borehole arrays has been proposed as an alternative to energy storage in shallow aquifers due to the scarcity of such aquifers in arid and semiarid regions and the management difficulties associated with the primary use of aquifers as water supply sources and potential recharge challenges (Bear et al., 1991). These SBTES systems have attracted growing interest owing to their numerous advantages over other energy storage systems (e.g., batteries, phase change materials, etc.). For example, they permit the storage of energy from renewable sources (solar thermal panels) and thus have lower environmental impacts, are space efficient (i.e., underground), and can be used in both rural and populated areas, and because they can store locally generated energy, they do not require long-distance energy transmission.

Despite the benefits of SBTES systems, an important issue limiting their cost-effective implementation is their low efficiency in terms of the extraction of energy injected into the system due to heat loss to the surrounding environment, which depends on the setting because most SBTES systems are insulated only at the ground surface. For instance, Zhang et al. (2012) predicted an energy recovery of only 27% after 10 yr of operation at the Drake Landing Solar Community in Alberta, Canada, which was later confirmed by Catolico et al. (2016) with a more in-depth analysis. The main parameters affecting the overall heat loss from a SBTES system include the system size and shape (number of boreholes and spacing), the heat injection–withdrawal scheme (i.e., heat transfer rate and duration), and soil thermal properties (Nordell and Hellström, 2000; Başer and McCartney, 2015).

There have been several experimental and numerical studies intended to better understand heat transfer in seasonal energy storage systems and increase the efficiency of these systems through analyzing the effects of the parameters listed above. Consequently, many numerical models have been developed to consider heat transfer at different geometric scales and dimensions (e.g., Eskilson, 1987; De Carli et al., 2010; Ozudogru et al., 2014; Lazzarotto, 2014; Lanini et al., 2014). Inspection of the parameterization of these models shows that most only consider subsurface heat transfer by conduction and that the soil thermal properties are constant in time and space and therefore do not depend on changes in the degree of saturation and/or temperature in the soil (e.g., Rees and He, 2013; Ozudogru et al., 2014; Diao et al., 2004). For SBTES modeling and design purposes, the thermal properties of the soil are often obtained using point measurements or system-level measurements such as a thermal response test, and in both cases are assumed constant for the entire heat storage system (Zhang et al., 2015). This assumption decouples heat transfer from soil moisture (water in both liquid and vapor forms) transport and air flow, whereas, in reality, the coupling is very strong, especially in the presence of large temperature gradients. The coupling in unsaturated soils may lead to a significant amount of heat being transferred by convection and latent heat transfer in addition to conduction.

Although convective heat transfer due to buoyancy-driven and groundwater flows have been shown to be significant in fully saturated soils at high temperatures, large injection rates, and high groundwater velocities (van den Brink and Hoogendoorn, 1983; Nordell and Hellström, 2000; Angelotti et al., 2014; Catolico et al., 2016), only a few studies have investigated the effect of convective heat transfer, latent heat transfer, and variable thermal and hydraulic properties either independently or together in seasonal storage systems installed in unsaturated soils (e.g., Bear et al., 1991; Reuss et al., 1997; Leong et al., 1998; Moradi et al., 2015a; Başer et al., 2016). For instance, Bear et al. (1991) developed a continuum-scale one-dimensional model with variable thermal properties that takes into account the effect of convective heat flux and latent heat transfer. Using analytical solutions of their one-dimensional model, they showed that heat storage systems in unsaturated soils could be unstable due to moisture redistribution in the system. Their study was among the first attempts to provide guidelines for a suitable choice of soil and initial conditions for heat storage systems. Reuss et al. (1997) performed a series of experiments on different soil types using a single vertical heat exchanger. These experiments were conducted to verify the theory of coupled heat and mass transfer developed by Philip and de Vries (1957) for temperatures up to 90°C. They showed that the thermal properties and overall thermal performance of seasonal heat storage systems change considerably with the degree of soil saturation. Moradi et al. (2015a) performed a set of two-dimensional laboratory experiments on different sandy soils under unsaturated conditions to evaluate the physics of coupled heat and mass transfer in unsaturated soils that are likely to be encountered in SBTES systems. In their experiments, no-mass flux boundary conditions were assumed for all boundaries of the two-dimensional bench-scale setup. Constant-temperature boundary conditions were considered for the left and right boundaries, and the top and bottom boundaries were thermally insulated. In each experiment, the measured ambient temperature was considered as the initial temperature for the entire domain. For the initial and boundary conditions assumed in their study, convective heat flux was shown to be considerably larger than conductive heat flux. Furthermore, they indicated that the initial degree of saturation, soil thermal conductivity, soil water retention curve (SWRC; e.g., Brooks and Corey, 1964; van Genuchten, 1980; Lu and Khorshidi, 2015) properties and porosity are important to properly model heat transfer under unsaturated conditions.

Although there are several field-scale seasonal heat storage systems in operation (e.g., Argiriou, 1997; Sibbitt et al., 2007; Bauer et al., 2010; Nussbicker-Lux, 2012; Bjoern, 2013; Başer and McCartney, 2015), the data necessary to properly develop and validate multiphase numerical models at these sites are difficult and costly to obtain. Furthermore, application of a fully discretized three-dimensional field-scale model is computationally expensive or, in most cases, impractical, albeit necessary to properly understand, for example, processes such as transient heat and mass transfer within heat exchanger arrays. Laboratory intermediate-scale (∼1–10 m) experimental and numerical studies offer the opportunity to investigate heat transfer mechanisms and fundamental processes (e.g., convective and latent heat transfer) under well-controlled initial and boundary conditions, providing precise data sets that can be used for model validation. In fact, these experiments will be beneficial for the comparison of numerical models and experimental data because the parameter estimation uncertainty will be reduced under controlled conditions (Oostrom et al., 2007). Intermediate scale, as seen here, is defined as the scale at which small-scale processes (i.e., evaporation and condensation) occur so that their relative significance in the overall system behavior (i.e., heat transfer here) can be determined (Lenhard et al., 1995). This requires that the environmental conditions are well controlled and that the experimental setup is compatible with the measurement techniques used. Among the previous seasonal energy storage studies, none to date includes a three-dimensional, intermediate-scale laboratory study under unsaturated conditions.

In this study, the influence of convective and latent heat transfer and variable soil thermal and hydraulic properties on heat transport in a laboratory, intermediate-scale SBTES system was studied. Numerical simulations were performed using a three-dimensional, fully coupled heat and mass transfer model that solves for multiphase (liquid and gas) flow and allows for nonequilibrium liquid water–water vapor phase change, as discussed below. To account for variable heat flux from U-shaped geothermal heat exchangers through which hot fluid is circulated to transfer heat into the soil, the fluid velocity inside heat exchangers was also modeled by solving the laminar Navier–Stokes equation. To validate the numerical model, a series of three-dimensional laboratory experiments with different initial conditions, injection temperatures, and flow rates was performed, simulating SBTES system operation. The setup was instrumented with a series of sensors to monitor the temperature and degree of saturation at high spatial and temporal resolutions. The validated numerical model was then used to evaluate the system behavior and quantify the significance and impact of variable soil thermal and hydraulic properties as well as different heat transfer mechanisms.

Materials and Methods

Experimental Apparatus

Laboratory-scale experiments were conducted in a three-dimensional tank with dimensions of 2.1 m long, 1.2 m wide and 0.6 m high (Fig. 1). Plywood panels (30 mm thick) were used to construct the tank structure, with aluminum t-slot beams providing external structural support. The inner part of the tank was lined with thin polyvinyl chloride (PVC) sheets (3 mm thick) sealed together with marine glue to assure a watertight seal. Polystyrene rigid foam was used to insulate the top surface of the tank. The foam was selected to represent an insulation layer, similar to those found at field-scale SBTES systems, which helps to restrict moisture and vapor exchange with or losses to the atmosphere. The sides and bottom of the tank were not insulated to mimic an unrestricted heat-flow boundary condition likely to be present in a field setting. The no-flow fluid boundary conditions do not represent the field setting, but no groundwater flow was applied in these experiments, and the zone of thermally induced water flow was expected to only be close to the heat exchangers and well away from the boundaries. Two constant-head devices were used to adjust and maintain the water table in the soil layer at a predetermined level during the experiments. Schematics of the three-dimensional tank, heat exchangers, and sensor bars are depicted in Fig. 1.

Four U-shaped heat exchangers (i.e., model boreholes) were constructed using stainless-steel tubing with inner and outer diameters of 3.8 and 6.2 mm, respectively. The heat exchangers were connected to a heated circulating bath (Polyscience Model AD07R-20) that supplied hot water to the U-shaped heat exchangers at a constant temperature and flow rate. All external tubing (∼1 m long) used to connect the boreholes to the heat bath circulator was wrapped with 5-mm-thick insulation tape to minimize heat loss before the inlets. A total of 24 dielectric sensors (Model 5TM, Decagon Devices) were used to measure the temperature with an accuracy of ±0.1°C and the volumetric water content with an accuracy of ±3%. These sensors were installed throughout the tank to measure the vertical and horizontal distributions of the degree of saturation and temperature with time. Each sensor was secured to a PVC rod (referred to here as a sensor bar) to allow for ease of installation and ensure that each sensor maintained its position while packing the tank with sand, as discussed below. The method developed by Sakaki et al. (2008) that accounts for sensor-to-sensor variability in analog-to-digital converter (ADC) counts (i.e., raw sensor outputs) was used to calibrate the dielectric sensors for volumetric water content readings. In this method, ADC readings under dry and water-saturated conditions are used in a two-point α-mixing model to calculate the degree of saturation for each sensor location. The location of each sensor can be seen in Fig. 1. As shown in Fig. 1, the inlet and outlet ports were instrumented with ST-100 thermistors (70 mm long, 6-mm-diameter rubber housing enclosed; Apogee Instruments) to monitor the temperature and consequently determine the heat transfer rate in each experiment. To monitor heat loss through the tank walls, a total of 50 ST-100 and ST-200 thermistors (Apogee Instruments) having an accuracy of ±0.1°C were installed on the inner and outer surfaces of the vertical walls. An AM16/32B multiplexer (Campbell Scientific) interfaced with a CP1000 datalogger (Campbell Scientific) was programmed to supply excitation voltage to measure the temperature. Room temperature was continuously monitored with an ST-100 thermistor located in the vicinity of the tank. The data from all the different sensors were collected every 15 min.

After installation of the sensors in the tank, sand was wet packed to the desired porosity (0.35) in 50-mm lifts using a procedure similar to that outlined by Sakaki and Illangasekare (2007) to achieve an initially saturated condition with minimal occluded air inside the sand. Deionized water was used as the pore fluid in all experiments.


Well-characterized uniform specialty silica sand (no. 12/20) under the trade name Accusand (Unimin Corp.) was used in all experiments. This sand was selected based on its hydraulic properties (i.e., air-entry suction) so as to allow the development of an unsaturated zone within the available height of the soil tank. The grain density of the sand is 2650 kg m−3, with a dominant mineral composition of quartz (99.8%). The uniformity coefficient for the sand is approximately 1.2. Hydraulic and thermal properties of this soil are: median grain diameter (d50), 1.040 mm; porosity (ϕ), 0.35 m3m−3; residual volumetric water content (θr), 0.017 m3m−3; saturated hydraulic conductivity (Ks), 3.76 × 10−3 m s−1; saturated thermal diffusivity, 11 × 10−7 m s−2; λdry = 0.04 W m−1 K−1; λsat = 3.05 W m−1 K−1; and van Genuchten (1980) water retention curve model parameters α = 0.00816 kPa−1 and n = 12.69. The hydraulic and effective thermal properties were determined using a small Tempe cell apparatus as outlined by Sakaki and Illangasekare (2007).

Experimental Procedure

This study consisted of four experiments involving soil layers having different degrees of saturation and different heat transfer rates controlled using different inlet heat exchange fluid temperatures and flow rates, as summarized in Table 1. The heat exchange fluid (water) temperatures (i.e., 70 and 80°C) were selected based on heat transfer rates that can be obtained by SBTES systems as well as the upper temperature operating limits of some of the sensors used in this study (Başer et al., 2015). The first experiment was conducted with a saturated soil layer (degree of water saturation Sw = 1) as a baseline case for comparison with the unsaturated cases. To perform the experiments on unsaturated soils, the water level in the tank was lowered to a predetermined level (i.e., 50 mm above the bottom of the tank) by adjusting the constant-head devices. For each experiment, after the desired saturation condition was achieved (either saturated or unsaturated), circulation of the heat transfer fluid was started through the soil borehole heat exchangers. Although not shown, Experiments EX1 and EX2 were replicated to demonstrate experimental repeatability.

Numerical Model Development

Although the experiments were designed to provide insight into the spatial and temporal distribution of temperature and degree of soil saturation, they do not permit conclusions to be drawn regarding heat and mass transfer processes and their impacts on SBTES systems. Therefore, numerical models are needed to obtain a more complete representation of heat transfer in these systems. By utilizing two- and three-dimensional models, more complex phenomena and additional transport mechanisms such as convective and latent heat transfer, rather than only conduction, can be considered (Zhang et al., 2015). In this study, a numerical model was used to investigate the impacts of different physical processes and assumptions on heat transfer inside the soil. For this modeling effort, a multiphase flow model previously developed by Moradi et al. (2015a) in COMSOL Multiphysics software was used in this study. In this model, conduction, convection and latent heat transfer are taken into account to simulate coupled heat and mass transfer in a two-dimensional bench-scale experimental setup. The model was modified to include a three-dimensional domain and different boundary and initial conditions. To better simulate the spatially and temporally variable heat flux from the boreholes, the fluid velocity inside the boreholes was also modeled. For an in-depth discussion on model assumptions and formulations, see Moradi et al. (2015a) and Smits et al. (2011). However, we include a brief description here of the main governing equations (Table 2) as well as new model developments.

Nonisothermal Multiphase Flow

Darcy’s law is used to describe the nonisothermal multiphase flow in the soil. Two separate partial differential equations (Eq. [1–2]) for both the liquid and gas phases are coupled through the capillary pressure (Pc = PgPw, where Pg is gas pressure and Pw is water pressure). The capillary pressure is related to the degree of saturation through the SWRC. The SWRC and relative permeability function are defined using the models of van Genuchten (1980) and van Genuchten–Mualem (van Genuchten, 1980), respectively. To account for the effect of temperature on the capillary pressure, the SWRC measured previously at room temperature, is modified using Eq. [3]. Moreover, the She and Sleep (1998) formulation (Eq. [4]) is used to describe the change in residual water content at elevated temperatures. Phase change occurs due to a difference in chemical potential between the liquid and gas phases (e.g., Chammari et al., 2003). Previous studies at the pore and representative elementary volume scales have shown that the equilibrium phase change assumption, in which the phase change is considered to occur instantaneously, is not suitable under certain conditions such as high temperatures (e.g., >50°C) and high air flow due to the hygroscopic properties of the porous medium (Cherblanc et al., 2007; Bénet et al., 2009; Smits et al., 2011; Ouedraogo et al., 2013; Trautz et al., 2015). Nonequilibrium phase change accounts for the finite time associated with phase change and consequently latent heat transfer (Trautz et al., 2015) and can better represent the dynamics of evaporation and condensation. For this reason, in this study, a nonequilibrium phase change formulation (Eq. [5]) was used. The formulation selected was used based on its stable model convergence and because it was successfully used in other tests (e.g., Zhang and Datta, 2004; Smits et al., 2011).

Water Vapor Transport

The mass balance for water vapor in the gas phase is given by Eq. [6]. The total diffusion coefficient of water vapor in the soil (Eq. [7]) is related to the binary diffusion coefficient in air (Eq. [8]) and the tortuosity (Eq. [9]). The Millington and Quirk (1961) formulation is used to express the tortuosity, and the vapor enhancement factor (Eq. [10]) of Cass et al. (1984) is introduced to address the underestimation of diffusive vapor fluxes by Fick’s law (Gurr et al., 1952) under nonisothermal conditions.

Heat Transfer

Conduction, convection, and latent heat transfer due to phase change are the three main heat transfer mechanisms in soils. Heat conduction is described by Fourier’s Law, which relates the time rate of heat conduction to the negative gradient of temperature through a proportionality constant known as the thermal conductivity (λ). Heat conduction in unsaturated soils is complicated due to the multicomponent nature of soil, consisting of soil grains, air, and water, thus the λt of soil is defined as an effective property that depends on the state of the system in both space and time. For instance, the λ of the soil components varies across three orders of magnitude (6–11 W m−1 K−1 for quartz minerals, 0.58 W m−1 K−1 at 25°C for water, and 0.026 W m−1 K−1 at 20°C for air). Because the λ of water is two orders of magnitude greater than that of air, an increase in the amount of water in the pore space leads to a higher λt. Therefore, in this study, λt was defined as a function of the degree of saturation. Convective heat flux, on the other hand, occurs when kinetic energy is transported by the bulk movement of the pore fluid. In our model, both wetting and nonwetting phases (i.e., water and air) are considered as pore fluids. Therefore, the velocities of both phases are used in the convective component of the heat transfer equation. In the experiments on unsaturated soils performed in this study, the contribution of liquid-phase velocity in the overall convective heat transfer was minimal due to the relatively small volume of the saturated region. Conductive, convective, and latent heat transfer are described through the energy balance equation (Eq. [11]). This equation assumes local thermal equilibrium among phases (i.e., Tgas = Twater = Tsoil = T), allowing effective thermal properties to be used. In general, the λt of the soil changes with the degree of saturation, and several mathematical models have been proposed to estimate this relationship based on easily measurable soil parameters (e.g., Côté and Konrad, 2005; Dong et al., 2015; Lu and Dong, 2015). In this study, λt as a function of the degree of saturation was determined separately using a small Tempe cell apparatus as outlined by Smits et al. (2010) and was adjusted based on the porosity of the domain using the Campbell model (Campbell et al., 1994). The effect of temperature on the effective thermal conductivity was neglected to simplify the numerical simulation, as discussed below. The effective heat capacity is expressed as a weighted sum of the heat capacity values of each phase, as given by Eq. [12]. Furthermore, the temperature dependency of the density and viscosity of water and air for nonisothermal conditions were defined and incorporated into the model through temperature-dependent nonlinear functions (Hillel, 1980; Lide, 2001).

Fluid Flow inside Borehole Heat Exchangers

To simulate the physics of heat transfer in closed-loop geothermal borehole heat exchangers, the numerical model should be solved for two separate domains in a coupled manner: (i) the time-dependent heat transfer in the soil outside of the boreholes; and (ii) the transient fluid flow and forced convective heat transfer inside the boreholes (Ozudogru et al., 2014). According to the experimental data in this study, the heat conduction from the heat exchangers’ walls toward the soil results in a decrease in water temperature flowing inside the heat exchangers as monitored by ST-100 sensors installed at the inlet and outlet ports of the boreholes. This temperature difference alters the outward heat flux from the surface of the heat exchangers. The heat inside the borehole heat exchangers is transferred by a laminar (Reynold number <1600), pressure-driven flow of fluid, which is coupled with heat transfer within the surrounding soil domain in the experiments. Coupling these two physics in the numerical model can lead to numerical complexities due to different time scales associated with transient heat transfer inside and outside of the heat exchangers. For instance, short time steps are needed to capture rapid changes in temperature inside the borehole, while for the surrounding soil, larger time steps may be sufficient. A domain decomposing approach might be used to address this problem, as discussed by Kim et al. (2011), but it is complex and outside the scope of this study. Instead, a simplified, decoupled approach was utilized in this study to decrease the computational time.

In our simplified approach, an initial simulation was first performed to obtain the steady-state velocity field inside the heat exchangers. It should be mentioned that the water properties (i.e., density and dynamic viscosity) were determined using the average inlet temperature (at steady state). Furthermore, the pressure drop due to viscous shear, friction heat, thermal resistance of the borehole wall, and the effect of temperature variation inside the heat exchangers on water properties were neglected. The velocity field was then implemented in the convective heat transfer component of the general heat transfer equation that solves for both domains inside and outside of the borehole heat exchangers (Eq. [11]). In fact, in the transient simulation of heat transfer for the coupled domains, the velocity field of the circulating fluid inside the boreholes was introduced as an initial condition for convective heat transfer. The Navier–Stokes and continuity equations for incompressible Newtonian fluid built into COMSOL were used to calculate the velocity field inside the boreholes (Eq. [15–16]). Because the heat exchange fluid circulation rate remained constant during the experiments, the fluid flow was assumed to be at steady state.

Initial and Boundary Conditions and Numerical Solver

Initial and boundary conditions were defined based on the characteristics of the soil tank experimental measurements. Neumann boundary conditions (no-mass flux) for liquid water and water vapor flow were assigned for all boundaries of the tank. The top boundary was thermally insulated (zero heat flux), and convective heat flux boundary conditions were assigned to the bottom and side boundaries (Eq. [14]). In this type of boundary condition, a heat transfer coefficient should be assigned for each convective surface (h in Eq. [14]). The heat transfer coefficient can be related to the geometry of the cooling surface (i.e., length scale), the physical properties of the cooling fluid (i.e., density and viscosity of air around the tank), and the temperature difference between the surface and external fluid bulk across a fictitious thermal boundary layer. To estimate these coefficients under natural convection, a set of correlations for vertical and horizontal surfaces (Incropera and DeWitt, 2002, Eq. 9.26.– 9.27 and 9.30–9.32) were incorporated. These correlations are defined based on three dimensionless numbers (Grashof number, Rayleigh number, and Nusselt number). At the inlet and outlet of boreholes, experimentally measured time series of temperature were used to prescribe Dirichlet boundary conditions. To create the initial saturation within the soil domain, the fully saturated domain was drained by applying a constant-pressure boundary condition at the bottom of the tank to match the initial experimental condition. Therefore, the initial water pressure was hydrostatic and the total gas pressure was set to atmospheric pressure. The initial vapor concentration was set to zero, and the average temperature measured inside the soil tank at the beginning of the experiment was considered the initial temperature for the entire domain including the fluid inside the boreholes.

The system of differential equations was solved using COMSOL Multiphysics software. The domain was discretized into 587,527 elements (1,415,590 degrees of freedom), with finer elements being used around the boreholes. Due to the large number of elements in the domain, the simulations were performed in a single node of a cluster with 64 GB of memory and 16 processors using a batch script in UNIX. Simulations were run for 8 d correlating with the experiment durations after creating the initial saturation condition based on the experimental measurements.

Results and Discussion

Temperature Behavior inside the Soil

To demonstrate the dynamics of temperature variations in all of the experiments, two sample points (A2 and D2 in Fig. 1) in the domain were selected. Figure 2 shows the temperature trends of these two points for all four experiments. Evaluation of the temperature inside the soil indicates that steady-state conditions were established within the first 3 d of the test for all unsaturated experiments and 5 d for the saturated case. The longer period for the saturated case is mainly due to the smaller thermal diffusivity of saturated soil at the location of A2 compared with unsaturated soil, which is confirmed by the numerical model results (11 × 10−7 vs. 13.5 × 10−7 m2 s−1). Because thermal diffusivity describes the rate at which heat is conducted through soil, in the saturated experiment the rate of temperature increase was smaller, as seen in Fig. 2a and 2b. It should be mentioned that we refer to the volume of the domain within the array of borehole heat exchangers as the impact zone. In a full-scale SBTES system, the impact zone would be the primary location for storing heat, although heat is also stored outside of the borehole array. The cross-section of the impact zone at the soil surface is highlighted in Fig. 1. Figure 2a shows the temperature change with time for all four experiments at the location of A2, located outside of the heat exchangers’ impact zone. As seen from the figure, temperature trends for the unsaturated experiments (EX2, EX3, and EX4) were similar due to similar thermal diffusivity. Note that the steady-state temperature at Point A2 for the saturated case was slightly higher than that of the two unsaturated experiments (EX3 and EX4), whereas they had a 10°C higher inlet temperature. Figure 2b shows the temperature trends for Sensor D2, located inside the boreholes’ impact zone in the middle of the domain. The temperature trends were slightly different inside the impact zone and outside. For instance, while the temperature outside of the impact zone (A2) in the fully saturated experiment (EX1) was almost the same as in EX3 and EX4 (with 10°C higher inlet temperature) except in the initial stages, a lower temperature was observed for the same experiment inside the impact zone. This implies that heat loss from the system in the experiment on the saturated sand was higher for two reasons. First, the saturated sand had a higher λt within the domain. Second, convective heat fluxes (buoyancy-driven water flow) around the boreholes, as shown below, resulted in lower temperatures inside the borehole array and consequently higher soil temperatures outside the impact zone. Similar behavior was observed in all of the sensors throughout the tank.

The measured ambient temperature varied slightly for each experiment (Fig. 2a) due to climatic variation in the laboratory where the experiments were conducted. The mean ambient temperatures were 23.7, 23.8, 24.2, and 23.9°C with variances of 3.26, 3.57, 3.39, and 2.97 for EX1, EX2, EX3, and EX4, respectively. As Fig. 2a shows, the ambient temperature did not have a direct impact on the transient behavior of the soil temperature because no diurnal fluctuation was observed in the temperature trends inside the soil. However, the ambient temperature affected the steady-state temperature through the natural convective heat flux from the tank walls (Eq. [14]).

Steady-state temperature profiles for sensors located along Transects A, B, C, and D (see Fig. 1) along with the temperature values for the ST-100 and ST-200 sensors installed on the inner and outer surfaces of the vertical walls are plotted vs. distance from the center of the tank in Fig. 3. Note that the ST-100 and ST-200 sensors were not exactly aligned with Transects A, B, C, and D. Transects A and C were located 120 mm below the soil surface (480 mm from the bottom of the tank), while Transects B and D were at a 480-mm depth (120 mm from the bottom of the tank). As shown in Fig. 3, higher temperatures were observed closer to the boreholes. This was mainly due to radial heat conduction and, in part, the heat loss through the uninsulated walls of the tank. It should be mentioned that the temperature did not change linearly between measuring points mainly due to radial heat conduction; the dashed lines are used only to represent the average temperature gradient that correlates with the slope of the connecting lines. Figure 3a shows the temperature profile along Transect A, which was close to the insulated surface of the tank. As seen in this figure, the temperature profile of EX3 and EX4 with 80°C injection temperature are almost the same. However, in EX1 and EX2 with an equal injection temperature of 70°C, different average temperature gradients were observed. Compared with the unsaturated experiment (EX2), in the experiment on the saturated sand (EX1), the average temperature gradient was slightly higher in the middle parts of the domain but decreased in close proximity to the boreholes. Specifically, outside of the boreholes’ impact zone, the temperature gradient for the saturated experiment was smaller than that of the unsaturated experiment (EX2). This was due to the λt of the saturated sand in EX1 being higher than the unsaturated sand in the other experiments and the effect of buoyancy-driven water flow, as mentioned above. No such deviation in temperature gradients were observed along Transect B, located 480 mm below the soil surface (Fig. 3b). This was because the sensors along Transect B were located in the capillary fringe (i.e., ∼100 mm of capillary fringe above ∼50 mm of a fully saturated region), where the soil had the same saturation and effective thermal properties in all experiments. Similar trends were observed along two other transects, as seen in Fig. 3c and 3d. Although not shown here, evaluation of steady-state temperatures for the rest of the sensors revealed that temperature trends in the saturated experiment started to deviate from the trends in the unsaturated experiment with the same injection temperature (EX2) toward the soil surface. This was due to the change in the degree of saturation and consequently the effective thermal properties of the soil. However, the temperature trends at all depths were almost the same for the last two unsaturated experiments (EX3 and EX4), which had equal injection temperatures (80°C) but different injection flow rates.

Total Heat Injected into the Soil

For each experiment, temperatures of the circulating fluid at the inlet and outlet of the boreholes were measured to estimate the total heat injected into the system. Although not shown here, the highest average temperature difference between the inlet and outlet (highest heat injection) was observed in Experiment EX1 (3.5°C) under fully saturated conditions, correlating with the higher effective heat capacity of the saturated soil in the domain. The lowest temperature difference (2°C) was seen in EX2, in which the inlet temperature of 70°C was applied to the unsaturated soil. Although Experiments EX3 and EX4 were performed under unsaturated conditions, the inlet temperature of both experiments was equally increased by 10°C compared with the first two experiments. In addition, the flow rate in EX4 was decreased by 50% with respect to the base flow rate used in EX1 to EX3, as highlighted in Table 1. A comparison of temperature differences between the inlet and outlet for the last two experiments showed that the decrease in flow rate resulted in a slight increase in temperature difference (2.4–2.6°C). This difference was not large enough to significantly impact the temperature distributions within the domain as discussed above. By numerically integrating over time, the total average heat injection can be determined using Eq. [17] in Table 2. In this equation, Δt (min) is the measurement time interval, n is the total number of time steps in the experimental data set, (kg min−1) is the mass flow rate of circulated water through the boreholes for the time interval, and c (J kg−1 °C−1) is the specific heat capacity of water. Figure 4 represents the results of these calculations. As indicated by Eq. [17], the total injected heat correlates with the flow rate and the temperature difference between the inlet and outlet of the heat exchanger. Therefore, similar behavior is expected.

Model Validation


The experimental results were used to validate the three-dimensional numerical model by comparing the predicted and measured temperatures. These comparisons at a single point, located in the middle of the domain (D2), for all experiments are shown in Fig. 5. Although not shown, similar results were obtained for other locations throughout the tank. The numerical results are generally in good agreement with the measured values regarding both magnitude and trends. The model provides a better prediction of temperature in the saturated experiment, as evident in Fig. 5a, while it slightly overestimated the steady-state temperature in the unsaturated cases (Fig. 5b, 5c, and 5d). This is in part due to the additional physics involved in the heat transfer processes in the experiments involving unsaturated sand (i.e., vapor transport, latent heat transfer, etc.) compared with the experiments involving saturated sand. These issues increase the model output uncertainty based on additional model inputs. Although the numerical model predicted the experimental data well, some discrepancies existed as statistically quantified by R2 values of 0.994, 0.975, 0.983, and 0.979 for EX1, EX2, EX3 and EX4, respectively. These discrepancies are, in part, due to model uncertainty and natural variability associated with mathematical modeling (Isukapalli and Georgopoulos, 2001; Moradi et al., 2015b). The model uncertainty pertains to simplifications made in the model structure and the constitutive relationships incorporated in the model. For instance, continuum-scale formulations of heat and mass transfer cannot fully capture the physics of multiphase flow caused by various pore-scale driving mechanisms such as capillarity and latent heat transfer. Furthermore, the decoupled approach that was used to model variable heat flux from the heat exchangers can introduce some error in model predictions. Natural variability is introduced by the spatial and temporal variability of environmental parameters and variables such as porosity, soil bulk density, thermal properties, and the vapor diffusion coefficient. For instance, soil bulk density increases during compaction and can influence the SWRC as a result (Assouline, 2006). In fact, all the model inputs in our numerical model were deterministic, with most of them being subject to environmental randomness.

Degree of Saturation

The simulated and observed degrees of saturation along Sensor Bars D and B at 0 and 8 d for EX2 are shown in Fig. 6. Comparison between predicted and measured values at different locations within the test tank shows that the model captured the drying behavior around the heat sources (Fig. 6b). However, some discrepancies exist, especially for Sensor Bar D, located in the middle of the tank. Although the simulated results demonstrate a decrease in the degree of saturation along Sensor Bar D, no considerable drying was observed in the experimental data. This is, in part, due to the accuracy of the sensors. The differences between measured and predicted values are within the accuracy of the sensors (1–2% in volumetric water content or degree of saturation) using medium-specific calibration. The range of variation for each sensor reading is graphically represented by error bars in Fig. 6. It should be mentioned that elevated temperatures could introduce another source of error by affecting the sensor electronics as well as changing the electrical conductivity of both water and air. However, as shown by Kizito et al. (2008), temperature sensitivity to soil water content for the family of ECH2O sensors used in this study is sufficiently small. Temperature effects are expected to be even less significant under residual saturation conditions due to the temperature dependency of air electrical conductivity and ADC counts being smaller than that of water. Similar drying patterns were observed in all simulations. It is important to note that, unlike the two-dimensional bench-scale results of Moradi et al. (2015a), there was no detectable thermally induced moisture increase in the experiments. This was probably due to the effect of the different boundary conditions and dimensionality in this work. In the three-dimensional experiments presented here, the amount of evaporated or condensed water was considerably less than the total amount of pore water available. Phase change occurred in a relatively larger volume and water vapor moved in three dimensions, thus decreasing the amount of water that condensed at a given location.

Convective and Conductive Heat Transfer

The magnitude of conductive and convective heat fluxes for EX2 are shown on a logarithmic scale in Fig. 7a and 7b at time t = 8 d. It should be noted that the magnitudes of both convective and conductive heat flux were significantly higher in the vicinity of the heat exchangers, and for this reason, a logarithmic-scale plot can better reflect the variability of heat flux within the domain. As seen from the results in Fig. 7a, the conductive heat transfer magnitude was greater near the heat exchangers due to larger thermal gradients. In addition, higher conductive heat was observed in the bottom of the tank where the soil had a higher degree of saturation and consequently higher λt. The regions with higher conductive heat flux did not extend to the middle part of the domain as the temperature gradient decreased due to symmetric placement of the heat exchangers. Similar to conductive flux, convective flux was greatest close to the heat exchangers. This was due to the higher density-driven flow of both phases resulting from the presence of higher temperature gradients.

Although both convective and conductive heat fluxes were small in the middle part of the domain, the convective flux was relatively larger (up to 80 times) than the conductive flux, as shown in Fig. 7c. The volume in which convective flux was dominant did not extend all the way through the impact zone, suggesting that SBTES systems can be engineered (i.e., through borehole heat exchanger configuration and spacing, injection temperature, and initial saturation condition) in a way that convective heat transfer could potentially play an important role in the overall heat transfer. It should be noted that the convective flux of the gas phase can also affect the latent heat transfer because they are strongly coupled under nonisothermal unsaturated conditions.

To represent the transient behavior of convective and conductive heat transfer under both saturated and unsaturated conditions, two points that are located at the same distance from a heat exchanger (BH1) outside and inside of the impact zone were selected for comparison. Figure 8 represents the change of these fluxes at Points D3 and G3 (shown in Fig. 1) as a function of time for first two experiments. In the saturated case (EX1) (Fig. 8a and 8b), the conductive flux is higher than the convective flux, showing that conduction was the main form of energy transfer under fully saturated conditions. Because the temperature gradients outside of the impact zone were higher than those inside the domain, higher conductive heat flux was observed outside of the boreholes’ impact zone, as shown in Fig. 8a and 8b. Buoyancy-driven water flow drives the convective heat transfer in the saturated case. This is due to temperature gradients impacting the density of the pore water and therefore creating convective water currents in close proximity to the heat sources. In addition to the temperature gradients, another factor affecting the convective currents is the soil permeability. Van den Brink and Hoogendoorn (1983) showed that for sandy soils with permeabilities >10−11 m2, these convective currents could adversely affect the efficiency of seasonal heat storage systems because they can carry heat upward and away from the storage volume. The sand used in the experiments has a permeability of 1.08 × 10−11 m2. Therefore, the high permeability of the sand along with the higher temperature gradients in the outer parts of the impact zone can result in greater buoyancy-driven water flow, as confirmed by the numerical simulation results (not shown), and consequently higher convective heat flux in the outer domain of EX1 than the impact zone.

Unlike the observations from the experiment on saturated sand, convective flux appears to be the dominant mode of heat transfer inside the impact zone, especially in the upper, less saturated parts of the domain in the experiment on unsaturated sand (EX2) (also shown in Fig. 7c). Conduction was higher outside of the impact zone, similar to that of the saturated case, due to higher temperature gradients. It should be noted that both of the selected points for EX2 were located in the unsaturated region. Therefore, convective flux occurred only due to the bulk movement of the gas through the pore space and not due to liquid movement as in the saturated case.

Phase Change and Latent Heat Transfer

In a dynamic, unsaturated, nonisothermal system, thermally induced vapor flow occurs due to two temperature-dependent mechanisms: diffusive transport that occurs because of a vapor concentration gradient and advective transport due to the bulk movement of the gas phase in the system. To understand the dynamics of vapor transport in this study, normalized vapor concentrations as well as total vapor flow directions at initial and final stages of EX2 are plotted in Fig. 9. The white regions in Fig. 9 represent the saturated zone and capillary fringe in which the air-filled porosity and consequently the vapor concentration is zero. As the system started to warm up, the residual water close to the heat exchangers began to volatilize, and a region with lower vapor concentration developed. During this time, the middle part of the impact zone as well as the outer part of the domain were cooler than the regions close to the heat exchangers. Because vapor generally moves from hotter to colder regions of the domain, it accumulated in the central region of the tank, leading to an increase in vapor concentration, as is highlighted in red in Fig. 9a. Similar behavior led to higher concentrations outside of the impact zone. As the temperature of the soil at the center of the tank approached steady state (after ∼2 d), a region with lower vapor concentration developed throughout the central domain, resulting in a uniform distribution of vapor, with a lower concentration in the impact zone (highlighted in green in Fig. 9b). However, due to the availability of pore water at the bottom of the tank, a region with high vapor concentration was observed above the saturated region. By the end of the simulation, the concentration gradient between the relatively dry zone in the middle and outside of the impact zone changed the direction of vapor flux, as illustrated by the white arrows in Fig. 9b. It is important to note that, in field-scale systems, it might take a long time to reach equilibrium, depending on the borehole configuration and soil thermal properties. In addition, this process might not occur in the same form in an open system (i.e., a field-scale SBTES surrounded by an infinite soil domain in which vapor can potentially escape the impact zone) as opposed to an isolated tank system with defined soil boundaries.

Latent heat transfer due to evaporation and condensation is one of the main heat transfer mechanisms in unsaturated soils, especially at elevated temperatures. Water vapor moves from hot boundaries and condenses as it reaches colder parts, thus releasing the absorbed energy. At high temperatures, latent heat transfer can significantly contribute to the overall heat flux (Cary, 1965). However, the significance of latent heat in the efficiency of SBTES systems installed in the vadose zone is yet to be understood. Obviously, phase change and consequently latent heat transfer depend on several state variables such as temperature, water content, gas pressure, and soil properties. To evaluate the dynamics of latent heat transfer in the system, the integral of the latent heat rate over the soil volume (J s−1) was calculated for EX2. Figure 10 shows the change in the total latent heat transfer rate within the soil domain as a function of time for the simulation period. The increase in the latent heat transfer rate represents energy absorption from the system, while the decrease correlates with the release of energy due to condensation. Lozano et al. (2009) showed that evaporation and condensation do not occur at identical sites of pore spaces, thus resulting in different water distributions and transport characteristics with time, which is referred to as the asymmetric behavior of evaporation and condensation. We do not have experimental evidence of such behavior in this study. In addition, it is not possible to confirm the effect of pore geometry using a numerical model because in a continuum-scale formulation such as the one used in current study, the effect of pore space geometry is not considered. However, investigating the dynamic of phase change (Eq. [5]) revealed that soil temperature has a minimal effect on the asymmetric behavior of the latent heat transfer rate, and it occurs mainly due to the combined effect of the Swveq − ρv) term, where ρveq is equilibrium vapor density and ρv is vapor density, in Eq. [5]. Apparently, as the system moves toward an equilibrium state, the gradient between the actual vapor concentration and equilibrium vapor concentration decreases and the latent heat transfer rate approaches zero. Similar behavior was observed in two other unsaturated experiments. It is important to note that the total mass of evaporated and condensed water (about 0.15 kg during 8 d of experiment) was significantly less (<1%) than total amount of water available in the unsaturated regions of the tank as residual water (about 18 kg). Therefore, for the experimental conditions of this study, the condensation of water vapor cannot significantly contribute to increasing the degree of saturation as discussed above. This was confirmed by both experimental observations and numerical simulations.

To quantify the significance of the total latent heat transfer with respect to the total heat injected into the system, depicted in Fig. 4, we calculated the area under the total latent heat rate–time curve. It appeared that in EX2, 9.3% of the total injected heat was transferred in the form of latent heat. We note that latent heat transfer is related to soil’s ability to retain water. Therefore, the contribution of latent heat to the overall heat transfer increases by an increase in the degree of saturation. On the other hand, an increase in soil saturation above a certain degree can hinder vapor transport and consequently latent heat transfer by decreasing the air-filled porosity. As shown in Fig. 9b, symmetric heating can lead to the development of a zone with a uniformly distributed low vapor concentration inside the impact zone. Under such conditions, latent heat transfer will be hindered after a certain time, which might change based on the size and configuration of the SBTES system as well as the initial and boundary conditions. Alternatively, an asymmetrical heat injection scheme might be useful in sustaining the evaporation and condensation processes within the impact zone. For instance, if heat injection from a few selected central heat exchangers is temporarily shut down, local cold regions will develop inside the impact zone. These cold regions can prevent water vapor from escaping the system and create a continuous circulation of vapor inside the impact zone. Sustaining evaporation and condensation processes within the impact zone can probably increase the time that a system needs to reach an equilibrium state, leading to an increase in the area under the total latent heat rate–time curve. Moreover, as shown by Zhang et al. (2015), a variable inlet temperature improves the thermal performance of borehole heat exchangers by alleviating the heat buildup around the borehole.

Thermal Property Variation within the Domain

Thermal properties are key parameters in analyzing heat transfer in unsaturated soils. It has been a common practice in most previous numerical studies of SBTES systems to assume constant thermal properties for the entire soil domain. However, the soil’s degree of saturation can significantly change the λt. In this study, the effect of temperature on λt was neglected according to the experimental studies by Smits et al. (2013). They found that for high and near-residual degrees of saturation at temperatures <50°C, the thermal property enhancement due to a temperature increase is insignificant. In all of our experiments, except in the immediate proximity of the boreholes, the soil temperature was <50°C.

To illustrate the change in λt during the unsaturated experiments, a volume plot of λt at times t = 0 and 8 d for EX2 is shown in Fig. 11. The variation in the initial degree of saturation due to self-weight drainage of an initially saturated soil tank to a predetermined water table (50 mm above the bottom of the tank) resulted in vertically variable λt (1–3 W m−1 K−1) as seen in Fig. 11a. By applying temperature gradients from heat exchangers, regions with a lower degree of saturation started to develop close to the hot boundaries, reducing the λt of the soil (Fig. 11b). Considerable drying adjacent to the heat exchangers can hinder heat transfer to the soil. Consistent with the modeling results, at lower degrees of saturation, λt is extremely sensitive to the degree of soil saturation while the effective heat capacity does not considerably change with saturation. In macroscopic (i.e., continuum) formulations, the effective heat capacity is defined as a weighted average of the heat capacities of all three phases (Eq. [12]). Therefore, at lower saturations, its dependency on the degree of soil saturation is minimal. However, unlike the effective heat capacity, λt is a function of the state of the pore fluid distribution and geometrical arrangement of phases (soil–water interaction regimes) at the pore scale rather than the amount of available water in the pore space (Cahill and Parlange, 1998).

Moradi et al. (2015a) evaluated the impact of a constant λt value on heat transfer under unsaturated conditions in a two-dimensional, bench-scale study. They showed that assigning a constant λt to the entire domain could result in unrealistic model predictions. In their study, it was shown that the temperature distribution in unsaturated soil (even near residual moisture contents) is sensitive to the thermal conductivity. To evaluate the impact of such assumptions at the three-dimensional intermediate scale domain, three simulations were performed. In these simulations, three values (0.5λsat, 0.75λsat, and1.5λsat) were assigned to the entire soil domain. All other hydraulic properties, initial conditions, and boundary conditions mimicked the base-case scenario (EX2). The trends in temperature for all scenarios at a single point (D2) are shown in Fig. 12 along with a comparison with the base case (black line) in which a variable λt was used. As seen in the figure, all three scenarios underestimated the temperature compared with the base case, especially in the case with higher λt (= 1.5λsat). This shows that a higher λt does not necessarily lead to higher temperatures inside the boreholes’ impact zone, which is due to an increase in heat loss from the soil. Although a higher λt increases the rate of heat transfer to the impact zone, it also simultaneously contributes to the lateral heat loss from the system. It is important to note that, in field-scale problems, the underestimation of temperature might not be as pronounced as in our intermediate-scale study due to differences in lateral convective boundary conditions and the degree of saturation profile. However, this example indicates that defining a proper λt for the domain could have a significant effect on temperature trends and that this should be a major design consideration for SBTES systems in the vadose zone.


This study examined heat transfer processes in a three-dimensional intermediate laboratory scale SBTES system under both saturated and unsaturated conditions. A nonisothermal, nonequilibrium model for coupled heat and mass transfer was developed and tested using experimental data. The model predictions agreed well with the measured values, suggesting that such a model can be used to study SBTES system behavior in unsaturated soils. Note that these results are specific to the soil type and experimental conditions used in this study and are not necessarily valid for all soil types or boundary and initial conditions. Evaluation of experimental observations along with the numerical model simulations led to the following conclusions:

  • Although limited, experimental results for the temperature distribution inside the soil showed that the steady-state conditions were established faster in unsaturated experiments than the fully saturated case due to the lower thermal diffusivity of saturated soil.

  • Simulation results showed that, in experiments on unsaturated sand, convective heat flux due to gas transport is larger (up to 80 times) than the conductive flux close to the thermally insulated soil surface and inside the boreholes’ impact zone. Increasing the convective heat flux due to gas transport increases the rate of latent heat transfer by increasing vapor transport within the domain.

  • The volume in which convective flux is dominant does not extend all the way through the impact zone. Although the convective heat flux due to buoyancy-driven water movement is higher close to the boreholes, conduction is the main form of energy transfer under saturated conditions.

  • For the test conditions considered in this study, latent heat transfer was responsible for about 10% of the total heat transfer in the experiments on unsaturated sand. Latent heat transfer was higher at the initial stages of heat injection and decreased as a relatively dry zone developed inside the borehole’s impact zone. This suggests that an asymmetric heat injection scheme (i.e., temporarily shutting down the heat injection from a few central heat exchangers) might be useful in sustaining the evaporation–condensation processes within the impact zone. By sustaining evaporation and condensation within the impact zone, we can take advantage of water as a heat exchanger. The water in the liquid phase absorbs thermal energy close to heat sources, travels within the domain, and releases that heat when it reaches colder regions. This process happens as long as water vapor does not escape the impact zone due to vapor diffusion.

  • Assigning a constant λt in modeling SBTES systems installed in unsaturated soils could introduce significant errors in the temperature distribution, and a model for λt that considers the dependency on the degree of saturation is needed for accurate SBTES simulations.


a fitting parameter

b fitting parameter for nonequilibrium phase change, s m−2

c fitting parameter for She and Sleep (1998)

cp heat capacity for the phase, J kg K−1

Dt total diffusion coefficient

Dv binary diffusion coefficient of water vapor in air

g gravitational acceleration, m2 s−1

h heat transfer coefficient, W m2 K−1

kint intrinsic permeability of soil, m2

krg relative permeability of gas (dimensionless)

krw relative permeability of water (dimensionless)

LRgw latent heat due to phase change

Mw molecular weight of water, kg mol−1

Pc = PgPw, capillary pressure in soil, Pa

Qs heat loss from the system, J m−3 s−1

R universal gas constant, J mol−1 K−1

Rgw nonequilibrium phase change rate between water and its vapor, kg m−3 s−1

Sw water degree of saturation (dimensionless)

Sg gas degree of saturation (dimensionless)

T soil temperature, K

T ambient temperature, K

u liquid velocity inside heat exchangers, m s−1

ug gas velocity, m s−1

uw liquid water velocity, m s−1

wv mass fraction of water vapor in the gas phase, kg kg−1

η vapor enhancement factor

θr residual water content (dimensionless)

λdry effective thermal conductivity of dry soil, W m−1 K−1

λsat effective thermal conductivity of saturated soil,W m−1 K−1

λt effective thermal conductivity of soil, W m−1 K−1

μg dynamic viscosity of gas, Pa s

μw dynamic viscosity of water, Pa s

ρ density of water inside heat exchangers, kg m−3

ρg density of gas, kg m−3

ρv vapor density, kg m−3

ρveq equilibrium vapor density, kg m−3

ρw density of water, kg m−3

σ surface tension, N m−1

τ tortuosity

ϕ total porosity of soil (dimensionless)

This research was funded by National Science Foundation (NSF) Sustainable Energy Pathways (SEP) Collaborative Project Award no. CMMI-1230544. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NSF. We acknowledge Paul Schulte from Colorado School of Mines for his support with the experimental setup.