The applicability of incorporating the heat-pulse and sensible heat balance methods into a fully coupled heat and mass transfer numerical model used to predict soil water evaporation rates is verified using precision experimental laboratory data. Predicted evaporation rates and cumulative water loss were in good agreement with experimental measurements during Stage 1 and Stage 2 evaporation.


A combined heat-pulse and sensible heat balance method can be used to determine evaporation using temperature measurements and thermal property estimations. The objective of this study was to investigate the applicability of the combined heat-pulse and sensible heat balance method by comparing laboratory experimental data to both an analytical and a multiphase heat and mass transfer model. A bench-scale laboratory experiment was performed to measure soil thermal and hydraulic properties at fine spatial and temporal resolutions. Comparisons of experimental and numerical results confirmed the applicability of the heat-pulse and sensible heat balance methods to determine evaporation rates. Results showed close agreement with experimental water loss measurements. This study demonstrated the ability and versatility of using the heat-pulse and sensible heat balance methods in numerical heat and mass transfer models to determine evaporation rates. Calculated soil thermal properties were in 93.4 and 97.5% agreement with experimental results for water content values >0.05. Deviations were observed at low water contents due to sensor sensitivity. The calculated evaporation rates yielded cumulative water losses that were 96.8 and 97.7% in agreement with experimentally measured weight loss data. Late Stage 1 evaporation was overestimated due to observed temperature rises by two of the heat-pulse probes. Despite this, the combined heat-pulse and sensible heat balance method is a powerful tool that can be used to determine evaporation rates in situ.

At a fundamental level, bare soil-water evaporation (evapotranspiration in the absence of vegetation) and condensation are the governing processes responsible for controlling mass and energy fluxes between the land surface and atmosphere in arid and semiarid climates. Evaporation is a vapor flux phenomenon that couples heat and water transfer in the shallow subsurface and is a key component in water balance analyses involving water exchange between the land and atmosphere (Berge, 1990). Understanding soil-water evaporation has direct application to many areas of engineering and scientific research including: earthen landfill cover design (Yanful et al., 2003; Scanlon et al., 2005), air sparging and vapor extraction (Gierke et al., 1992; Armstrong et al., 1994), agricultural water management (Ventura et al., 2001; Fereres et al., 2003), and climate modeling with respect to land–atmosphere interactions (Judge et al., 2003; D’Odorico and Porporato, 2004; Mosthaf et al., 2011). Despite the importance of soil-water evaporation in these contexts, the phenomenon is still not fully understood. Soil-water evaporation is affected by complex interactions between atmospheric demand (humidity, temperature, radiation, and air flow) (van de Griend and Owe, 1994) and internal transport mechanisms in the soil in addition to the soil properties (e.g., vapor diffusivity, thermal and hydraulic properties), all of which are strongly coupled (Prat, 2002; Sakai et al., 2011). Due to the inherent complexities and difficulties of measuring evaporation rates in addition to a scarcity of existing field or laboratory data capable of testing and refining the existing energy and mass transfer theories at high spatial and temporal resolutions, the mutual interaction of energy and mass transfer is rarely considered in most models or practical applications (Smits et al., 2011).

Evaporation from bare soil involves two separate processes, diffusion of water vapor from the evaporative or drying front (the location marking the transition from partially saturated to completely dry soil conditions) within the soil profile and its removal by turbulent air flow along the land–atmosphere interface. This phenomenon occurs in three distinct stages, referred to as Stages 1 to 3, that describe the observed location of the drying front and the dominant transport mechanisms (Fisher, 1923; Lemon, 1956). Stage 1 is defined by an initially high and near-constant evaporation rate controlled primarily by atmospheric demand (Shokri and Or, 2011), which continues as long as the evaporative front remains in hydraulic connection to the water source via capillary action (Lehmann et al., 2008). Once the hydraulic connection is severed, evaporation becomes vapor-diffusion limited (Stage 2), marked by a sudden decrease in the evaporation rate (Lehmann et al., 2008). At this point, evaporation is no longer under the sole influence of atmospheric demand; instead, internal transport mechanisms and soil properties control the evaporation process (Shokri and Or, 2011). Stage 3, also vapor-diffusion limited and therefore often combined with Stage 2, occurs when the evaporative front is located far below the soil surface.

Many of the currently used physical measurement and analytical techniques assume that soil-water evaporation occurs at the soil or land surface (Heitman et al., 2008b, 2010). This permits the use of micrometeorological methods like the Bowen ratio energy balance or eddy covariance to indirectly quantify evaporation rates above the ground (Warrick, 2002). These approaches, however, do not take into consideration the inherent complexities that are associated with soil-water evaporation, such as the dependence on atmospheric conditions (van de Griend and Owe, 1994) and soil properties (Prat, 2002), and the constantly changing location of the evaporative front (Yamanaka et al., 1998). Microlysimetry is a commonly used in situ soil-water evaporation measurement technique that helps reduce the number of assumptions made in these methods (Evett et al., 1995; Daamen and Simmonds, 1996). The drawback of this technique is significant soil disturbance during installation of the microlysimeters, which affects how representative measurements are of undisturbed samples.

Another less known and relatively new method for in situ estimation of soil-water evaporation uses a novel sensor technology called a heat-pulse (HP) probe. These sensors are used to measure changes in soil temperature in response to a small measured heat input. These terms are, in turn, used to calculate the soil thermal properties (thermal diffusivity, volumetric heat capacity, and thermal conductivity) according to an approach referred to as the HP method (Campbell et al., 1991; Bristow et al., 1993, 1994; Heitman et al., 2008a). The thermal property data are then used to calculate the sensible heat flux (Cobos and Baker, 2003; Ochsner et al., 2006) and sensible heat storage (Ochsner et al., 2007) components of the sensible heat balance (SHB) equation devised by Gardner and Hanks (1966) and improved on by Heitman et al. (2008a, 2008b) and Sakai et al. (2011), ultimately yielding 
where E (m s−1) is the evaporation rate, L (J m−3) is the volumetric latent heat of vaporization, λ (W m−1 °C−1) is the soil thermal conductivity, T (°C) is the soil temperature, z (m) is depth (z = 0 at the soil surface), and ΔS (W m−2) is the change in sensible heat storage. The subscript i indicates the soil layer depth increment as shown in Fig. 1. Equation [1] shows the usefulness of this approach in terms of its independence from the soil’s hydraulic properties or coupled heat and water transport coefficients (Heitman et al., 2008b).

The HP probe in its most basic form consists of two needles (Campbell et al., 1991), with one needle containing a heating element (commonly generating 8-s, 750 J m−1 pulses of heat; Kamai et al., 2008) and the other needle containing a thermocouple or thermistor for measuring temperature. This sensor design has since been modified by others (e.g., Ren et al., 2000; Xiao et al., 2012; Deol et al., 2012). Ren et al. (2000) created a three-needle HP probe that contained a thermocouple in the outer needles and a heating element in the center needle. The Heitman et al. (2008a) three-needle model adopted in the present study contains an additional thermocouple located in the center needle with the heating element to help increase the spatial resolution by a factor of two.

A large body of literature exists with respect to the application of the HP method to help determine soil thermal properties (e.g., Campbell et al., 1991; Kluitenberg et al., 1993; Bristow et al., 1994; Bristow, 1998; Hopmans et al., 2002; Mori et al., 2003) and soil moisture (e.g., Campbell et al., 1991; Bristow et al., 1993; Basinger et al., 2003; Heitman et al., 2003). Thermal and soil moisture data collected using HP probes have been applied in numerous studies to analyze soil-water flow (e.g., Ren et al., 2000; Ochsner et al., 2005; Mortensen et al., 2006; Saito et al., 2007; Kamai et al., 2008). To date, however, there have been few studies that used the combined HP and SHB method, whether through numerical simulation (e.g., Sakai et al., 2011) or field or laboratory experimentation (e.g., Heitman et al., 2008a, 2008b; Xiao et al., 2011; Deol et al., 2012). Sakai et al. (2011) performed a numerical experiment, applying generated parameters (temperature, relative humidity, water content, and soil properties) from a numerical simulation as inputs for the SHB. The results demonstrated that the SHB method can capture the two-stage (Stages 1 and 2) evaporation behavior of soil water. Sakai et al. (2011) also found that smaller observation grids (sensor spacing) led to more accurate prediction of evaporation rates. Heitman et al. (2008a, 2008b) performed a field experiment where 10 HP probes were installed within the top 7.2 cm of an agricultural topsoil. Throughout the duration of the experiment, the site was repeatedly exposed to natural rainfall events, diurnal temperature cycles, and varying levels of solar radiation. The results showed that evaporation rates, calculated as a function of depth and time, were in good agreement with other traditional measurement techniques, i.e., microlysimetry and the Bowen ratio. The sensitivity of these sensors was also demonstrated in terms of their ability to accurately capture transient precipitation events; however, the lack of control over environmental forcings in field settings results in complicated data. More recently, Deol et al. (2012) performed a laboratory column experiment to evaluate the use of an 11-needle heat pulse probe to measure temperature and thermal property distributions.

The objective of this study was to investigate the applicability of the combined HP and SHB method using precision laboratory data in both an analytical and a fully coupled heat and mass transfer model. We constructed and outfitted a soil tank with an array of state-of-the-art sensor technologies for the continuous and autonomous collection of soil water content, temperature, and relative humidity data at high spatial and temporal resolutions. The tank was packed with a uniform test sand for which the hydraulic and thermal properties have been characterized (Sakaki and Illangasekare, 2007; Smits et al., 2010). We incorporated the experimental temperature and heat data collected from the HP probes into a numerical MATLAB model containing the HP and SHB methods. The applicability of integrating the SHB method into fully coupled numerical models (see Sakai et al., 2011) was tested using the aforementioned experimental data set for validation. This was performed in a modified, fully coupled, nonisothermal, nonequilibrium multiphase flow (heat, liquid water, and water vapor) COMSOL model first developed by Smits et al. (2011). The same atmospheric boundary conditions and soil properties of the experimental portion of the present study were used as inputs for this numerical model.

Mathematical Theory

Heat-Pulse and Sensible Heat Balance Methodology

Soil-water evaporation studies commonly make the simplifying assumption that latent heat flux (LE) and, as a result, evaporation occur exclusively at the soil surface (Mayocchi and Bristow, 1995). Gardner and Hanks (1966) proposed that the sensible heat balance could be modified so that soil-water vapor flow (evaporation) can be determined at any depth below the soil surface: 
where L (J m−3) is the latent heat of vaporization; E (m s−1) is evaporation; H0 and H1 (W m−2) are sensible heat fluxes measured at two different depths; and ΔS (W m−2) is the change in sensible heat storage between the two depths. The latent heat of vaporization is calculated according to (Forsythe, 1964): 
where Tmean (°C) is the mean temperature for a given soil depth layer as a function of time.
The sensible heat fluxes (H0 and H1) are calculated using Fourier’s law of conduction, which is defined as the heat flow (dQ/dt) per unit cross-sectional area or the product of the negative temperature gradient and the thermal conductivity (λ, W m−1 °C−1) (Parker, 2003). The incorporation of Fourier’s law of conduction can be seen in Eq. [1]. The change in sensible heat storage is dependent on the change in temperature, depth, volumetric heat capacity (C, J m−3 °C−1), and time (t, s) (Ochsner et al., 2007): 
where the subscript j denotes the time step and the subscript i again indicates the soil depth layer.
Equations [1] and [4] are both dependent on two inherent soil thermal properties, λ and C; λ is calculated as the product of C and thermal diffusivity (κ, m2 s−1) (Bristow et al., 1994). To determine C, there are many different formulations available (e.g., Campbell et al., 1991; Kluitenberg et al., 1993; Knight and Kluitenberg, 2004). Campbell et al. (1991) first applied a temperature distribution around an instantaneously heated line source modeled after Carslaw and Jaeger (1959) to simulate the heating element in a heat-pulse needle. In reality, however, the heating elements are actually finite heated line sources that are not instantaneously heated. For this reason, Kluitenberg et al. (1993) applied the temperature solution for a finite heating period based on the work of de Vries (1952) to derive the equation for the volumetric heat capacity: 
where q′ (W m−1) is the energy per unit length of heater in the probe, t0 (s) is the duration of the heat pulse, r (m) is the radial distance from the line source, κ (m2 s−1) is the thermal diffusivity, and tmax (s) is the time required to reach the maximum temperature change ΔTmax (°C). The term −Ei(−x) represents the exponential integral function with input x. Bristow et al. (1994) recognized that Eq. [5] is a global maximum; when its derivative is taken and set equal to zero, κ can be determined as: 
Errors associated with Eq. [6] that result from uncertainties in ΔTmax and tmax data are compounded in the calculation of C as shown by Eq. [5]. A second formulation developed by Knight and Kluitenberg (2004) that uses a Taylor series approximation and is independent of κ is used to address this problem of compounding error: 
where e is the natural logarithm base and ε = t0/tmax.
In addition to indirectly measuring soil-water evaporation (Eq. [1]), HP probes are used to calculate the volumetric water content (θw, m3 m−3). The volumetric water content can be approximated as a function of the volumetric heat capacity of water (Cw, J m−3 °C−1), θw, the soil bulk density (ρb, kg m−3), and the specific heat of the soil minerals (cs, J kg−1 °C−1) (Campbell et al., 1991): 

The greatest obstacle to using this approach to determine water content in the field is the required foreknowledge of the bulk density and specific heat of the soil minerals—this is not a problem for laboratory experimentation because the soil type and packing conditions can be controlled. Bristow et al. (1993) showed that the applicability of this approach is moisture dependent. Error increases significantly with decreasing water content as the soil dries, limiting the use of this approach at low water contents (Bristow et al., 1993).

Multiphase Heat and Mass Transfer Model

In this work, we used a fully coupled numerical model developed by Smits et al. (2011) to solve the governing equations for heat, liquid water, and water vapor transport in the soil. The code implemented a nonisothermal solution that accounted for nonequilibrium liquid–gas phase change with gas-phase vapor diffusion. Smits et al. (2011) demonstrated the importance of accounting for nonequilibrium phase change in simulation efforts to capture gas-phase transport under transient field conditions. Below follows a summary of the Smits et al. (2011) model used in the present study; see Smits et al. (2011) for an in-depth discussion.

The mass balance equation for liquid water in porous media is given by (Bear, 1972): 
where ρw (kg m−3) is the density of water as a function of temperature, θw (m3 m−3) is volumetric water content, pc (Pa) is capillary pressure, uw (m s−1) is the mean pore velocity of liquid water calculated according to Darcy’s law, and Rgw (kg m−3 s−1) is the rate of phase change between liquid water and water vapor through the process of either evaporation or condensation and is given by (Bixler, 1985; Zhang and Datta, 2004): 
where the term in square brackets is the mass transfer coefficient in which b (s m−2) is an empirical fitting parameter, θrw (m3 m−3) is the residual water content, R (N m mol−1 K−1) is the universal gas constant, T (K) is temperature, and Mw (kg mol−1) is the molecular weight of water. The coefficient ρg (kg m−3) denotes the density of the gas phase expressed according to the ideal gas law, wv is the mass fraction of water vapor in the gas phase, cvs (kg kg−1) is the temperature-dependent saturated vapor concentration in the gas phase, and Hre (Pa Pa−1) is the relative humidity calculated according to Kelvin’s equation; the product of cvs and Hre is the equilibrium water vapor concentration.
Similar to the mass balance for liquid water, a mass balance is established for the vapor in the gas phase: 
where θg (m3 m−3) is the volumetric gas constant, ug (m s−1) is the mean pore velocity of the gas phase (calculated according to Darcy’s law [Bear, 1972]), and Dv (m2 s−1) is the effective vapor diffusion coefficient, defined as the product of the vapor diffusion coefficient in air, tortuosity, and an enhancement factor (η). The Penman (1940) tortuosity model (τ = 0.66 θg) and the vapor enhancement factor described by Cass et al. (1984) and Campbell (1985) were applied in the present study.
Energy conservation is written at a macroscopic scale for the entire system containing three phases (solid, liquid, and gas): 
where Cg (J kg−1 K−1) is the specific heat capacity of gas, λT is the effective thermal conductivity of the soil (Campbell et al., 1994), and Qs (J) is the heat loss. The energy conservation Eq. [12] is coupled to the combined mass balance equations for liquid water and water vapor phases Eq. [9] and [11] to solve for the total water-phase mass transport, yielding: 
where φ (m3 m−3) is the porosity of the soil.

Materials and Methods

Soil Material

In the present study, we used a uniform sand, Accusand no. 30/40 (Unimin Corp.), for which the soil thermal and hydraulic properties are well characterized (Smits et al., 2010, 2011, 2012). According to the manufacturer, the sand is 99.8% quartz, has a uniformity coefficient of 1.2, a mean grain density of 2.66 g cm−3, and is rounded in shape.

Experimental Apparatus

For the experimental portion of this study, a rectangular tank (9.2 cm wide, 58.4 cm long, and 45.7 cm tall) was constructed using 1.27-cm-thick polycarbonate. The tank was outfitted with an array of sensors to continuously and automatically monitor the θ and T distributions from the soil surface to a depth of 27.94 cm below ground surface (bgs); the sensors were installed horizontally through the tank walls at 5.08-cm increments starting at a depth of 2.54 cm bgs. The HP probes, described in more detail below, were installed across the depth interval of 0 to 4.8 cm bgs, which allowed very shallow temperature gradients to be captured at high spatial and temporal resolutions.

Thirteen ECH2O EC-5 dielectric soil moisture sensors (Decagon Devices) and 11 EC-T temperature sensors (Decagon Devices) were installed to measure the water content and temperature, respectively. The EC-5 sensors were calibrated according to the method presented by Sakaki et al. (2008). Temperature and relative humidity were measured (EHT temperature/relative humidity sensor, Decagon Devices) at the soil surface and at a height of 40.6 cm above the soil surface. Two EHT temperature/relative humidity sensors located at the soil surface were placed in close contact with the sand to reflect the temperature and humidity conditions directly at the soil surface, whereas the sensor located at a height of 40.6 cm captured the ambient atmospheric conditions. The data collected by the three relative humidity sensors were recorded using ECH2O (Decagon Devices) dataloggers.

Four HP probes installed in the alternating pattern shown in Fig. 2, constructed and calibrated by East 30 Sensors, were used in the present study. These probes, modeled after those of Heitman et al. (2008a, 2008b), consisted of three stainless steel needles (1.2-mm diameter, 30 mm long) arranged in parallel within an epoxy body as shown in Fig. 1. A thermistor, used for measuring temperature, was installed in each needle; the center needle also contained a resistance heater for generating a heat pulse. The 8-s, 750 J m−1 heat pulses were generated by controlling a direct-current power supply via a CR1000 datalogger (Campbell Scientific).

The needle spacing (the distance between two adjacent needles) of the HP probes used in the present study was 3 mm. This spacing was chosen in response to the numerical study performed by Sakai et al. (2011), who showed that a smaller needle spacing increases the resolution at which evaporation rates can be captured by decreasing the undetectable zone (the layer of soil located between the top needle and the midpoint between the top two needles of any given HP probe); this has also been addressed using a different sensor structure (see Deol et al., 2012; Xiao et al., 2012). This needle spacing is smaller than the industry standard (6 mm) and therefore may introduce error as a result of the physical volume of the needles (Kluitenberg et al., 1995; Knight et al., 2012). The expected error, based on the theoretical error analysis of Kluitenberg et al. (1995), which assumes that the HP probes behave as infinite line sources, ranges from 0.01 to 12% for thermal diffusivity and 0.01 to 5% for volumetric heat capacity. This compares with the expected error ranges for the 6-mm spacing of 0.01 to 8 and 0.01 to 2% for thermal diffusivity and volumetric heat capacity, respectively. The significance of this error is discussed below.

The HP probes were calibrated in saturated soil according to the method devised by Campbell et al. (1991) and discussed in detail by Ham and Benson (2004). The saturated volumetric heat capacity value used in the calibration process to determine the apparent needle spacing was determined using two KD-2 Pro thermal property analyzer 30-mm dual-needle heat pulse sensors (SH-1, Decagon Devices). In addition to measuring the temperature and heat input data collected during the experiment using the heat pulse probes, λ, C, and k as a function of θw were determined independently in a separate small column experiment using the SH-1 and KD-2 Pro.

In a separate series of tests, the four HP probes were operated under saturated and dry conditions. The thermal properties (κ and C), as determined using the HP methodology, were compared with KD2 Pro estimates using two different SH-1 sensors. For each measurement, the same thermal sensor was installed in the sand, the sample was compacted, and the temperature heat inputs were recorded to calculate the thermal properties. The sensor was then removed and the entire procedure was repeated 10 times; all measurements were recorded under identical conditions (i.e., ambient temperature and packing porosity) to allow precise comparison. The mean values of the calculated κ and C values were found to lie within the 95% confidence intervals of the κ and C values measured using the SH-1 sensors under saturated conditions. The calculated values fell within the 90% confidence level under dry conditions, which is not unexpected given the known inaccuracy of these probes under these conditions (Bristow et al., 1993). The correlation of the HP probe and SH-1 sensor show that the 3-mm needle spacing did not lead to significant errors in thermal property estimation.

The experimental portion of the present study was performed in a climate-controlled, closed-circuit wind tunnel located at the wind tunnel/porous media test facility at the Center for Experimental Study of Subsurface Environmental Processes (CESEP) at the Colorado School of Mines in Golden, CO. The wind tunnel, capable of sustaining airflow between 0.865 and 10 m s−1, has a test section with a 1- by 1-m cross-sectional area and length of 7.6 m. The test section itself is located within a climate-controlled chamber that allowed the soil tank to be maintained at constant temperature. The soil tank was interfaced with the wind tunnel by removing the floor and raising the tank until the soil surface was flush with the test section floor (Fig. 3). The test section was outfitted with a pitot-static probe (Davis Instruments) and Laser Doppler Velocimetry system (Innova 70C Argon laser, Coherent) to measure and monitor the wind speed.

Experimental Procedure and Data Acquisition

The soil tank was wet packed using no. 30/40 sand and deionized water to a final porosity of 0.327 m3 m−3. Sand was added to the tank in 1-cm depth increments, maintaining approximately 5 cm of water above the sand surface to ensure saturation and a uniform bulk density throughout. After each sand layer was in place, the tank walls were repeatedly tapped with a rubber mallet to further compact the sand as described by Sakaki and Illangasekare (2007). This packing approach allowed greater densities to be achieved in comparison to the standard ASTM Method D 4253-00 (Youd, 1973). Vibratory compaction (American Society for Testing and Materials, 2006) was not used so as to prevent damaging the sensors. Once the tank was fully packed, it was covered with a plastic sheet to prevent water loss from occurring before the start of the experiment. The tank was placed on a weighing scale (Sartorius Model 11,209-95, 65 kg ± 1 g range) to measure mass losses of water.

The experimental portion of the present study was conducted for a total of 10 d under steady conditions of approximately 1.1 m s−1 air flow and near-constant diurnal ambient temperature. For the duration of the experiment, the mass loss of water, water content, temperature, and relative humidity were sampled every 10 min. The HP probes operated on a 15-min sampling cycle to allow ample time for the temperature of the surrounding soil to recover to ambient temperature between heating events (e.g., Smits et al., 2012; Heitman et al., 2008a, 2008b; Sakai et al., 2011). The temperature was measured using these probes at 1-s intervals for a 5-min duration starting at the beginning of each heating cycle. The initial temperature of the soil at the start of each heating cycle was subtracted from the maximum recorded temperature to yield ΔTmax. The MATLAB “max” command was used to automatically locate the maximum temperature value and its associated time (tmax) of occurrence for each probe’s heating cycles.

Modeling Procedures

Thermal properties and sensible heat balance components as described in Eq. [3], [4], [5], [11], and [12] were calculated using MATLAB (The MathWorks). This software was chosen for its capability of handling large data sets quickly and efficiently and its ability to read several file formats. The temperature and heat input data generated using the HP probes were therefore easily imported into MATLAB as a single file and stored as two matrices: one for temperature and another for heat. The calculation of the HP and SHB components using these data in the MATLAB analytical code is referred to here as Approach 1.

A finite element method based software, COMSOL Multiphysics, was used for the numerical simulation portion of this study. As discussed above, we used a model first developed by Smits et al. (2011) that couples nonisothermal, nonequilibrium multiphase flow (heat, liquid water, and water vapor), herein called Approach 2. In this model, we constrained all mass and energy losses to occur from the top boundary through the careful definition of the boundary conditions. As in the experimental portion of the present study, the sides and bottom of the simulated tank were treated as no-flow boundaries. The top boundary was time invariant with respect to liquid water flow and total gas-phase pressure, whereas it was time dependent with respect to temperature and vapor concentration. A no-flow boundary was assigned to the liquid water, and the total gas-phase pressure was assumed to be equal to the atmospheric pressure. The sides of the tank could be treated as thermally insulated boundaries due to the experimental apparatus being located within the climate-controlled chamber. A heat sink term was included in the model to account for any heat loss from the experimental apparatus. The experimental soil surface temperature and relative humidity and the internal temperature were used as the time-dependent inputs for the temperature and vapor-concentration boundary conditions at the top and bottom of the tank, respectively. A total of 150, one-dimensional elements consisting of 151 grid points, resulting in a 3-mm mesh size, were used in this second approach. This corresponds with the same physical needle spacing used in the experimental portion of the present study. UMPACK, a direct matrix solver incorporated into the COMSOL software, was used to implicitly solve the model. To remain consistent with our experimental analysis, this simulation was allowed to run for a total of 10 d. The model was used to generate the necessary inputs (water content, thermal conductivity, and temperature) at a mesh size equivalent to our physical needle spacing; see Sakai et al. (2011) for an in depth description of the modeling approach.

Results and Discussion

Temperature Analysis and Heat-Pulse Probe Data

Understanding temperature variations in this work is critical. The ambient air temperature within the wind tunnel fluctuated ±0.5°C diurnally and increased a total of 4°C by the end of the experiment. This warming trend was observed by all EC-T temperature sensors and HP probes (Fig. 4) in addition to being predicted in the multiphase heat and mass transfer model (Approach 2). Figure 4 provides a comparison of temperature measurements made at depths of 2.40 and 2.54 cm bgs using a HP probe and EC-T temperature sensor, respectively. The general trends of the two measurements are very similar, the only difference being a slight difference in magnitude (~0.2°C). This difference remained constant throughout the experiment and can be explained in terms of the difference in depth of the sensors and therefore slight differences in temperature behavior. Liquid water has a greater heat absorbing capacity (specific heat ~4.187 kJ kg−1 K−1) than air (specific heat ~1.005 kJ kg−1 K−1) (Lide, 2003). Therefore, as the drying front moves deeper, the “temperature buffer” previously provided by liquid water no longer exists, in turn leading to an overall temperature increase within the unsaturated soil profile.

Calculated Soil Thermal Properties and Volumetric Water Content

To further verify the HP probe measurements, HP data were used to calculate the soil thermal properties κ, C, and λ using Approaches 1 and 2 and compared with independently measured soil thermal properties obtained by Smits et al. (2010). Smits et al. (2010) performed a small column experiment to determine the soil thermal properties as a function of soil moisture of no. 30/40 sand under boundary conditions similar to those of the present study (see Smits et al. [2010] for a detailed description of their experimental study). Figure 5 shows a comparison of the thermal properties (κ, C, and λ) of the test sand as a function of water content for Approaches 1, 2, and Smits et al. (2010). Data obtained from the HP probe located at a depth of 2.4 cm bgs were used for Approaches 1 and 2. Although not shown, the other HP probes showed consistent results.

Figure 5 shows that the general shape of the three curves are similar, although the magnitude is slightly different. Approaches 1 and 2 closely follow the curves of Smits et al. (2010) across the intermediate and saturated water contents but tend to deviate at low water contents (at and below residual saturation). In analyzing the behavior by comparing the statistical results, the modified index of agreement (MIA) values between Smits et al. (2010) and Approaches 1 and 2 ranged from 93.4 to 97.5% for water content values >0.05. The overall differences in the thermal properties were most likely due to differences in packing conditions (i.e., soil porosity) between this study and Smits et al. (2010). Although the models captured the thermal properties with time across the majority of water contents reasonably well, the lack of agreement, especially at low saturations, demonstrates that there are still some limitations in the approaches. Deviations in HP probe thermal property measurements at low water contents have been similarly observed by Bristow et al. (1993), Bristow (1998), and Basinger et al. (2003).

The linear trend of the volumetric heat capacity predicted by Approach 2 in Fig. 5b was expected; volumetric heat capacity and water content are linearly related according to Eq. [8]. The observed increase in volumetric heat capacity in Approach 1 (Fig. 5b) at water contents between 0.01 and 0.05 is the direct result of a large increase in the maximum temperature recorded between Days 7.8 and 10. This in turn led to a smaller ΔTmax and therefore larger C. The HP sensors are extremely sensitive to soil moisture and temperature at low or near residual water contents (Bristow et al., 1993). This fact explains why the temperature increase between Days 3 and 5.4 did not necessarily result in an increase in volumetric heat capacity. If a similar change in the maximum temperature of a given heating cycle had been observed, then there would have been an increase in volumetric heat capacity as well. An increase in thermal diffusivity within this water content range (Fig. 5a) was not observed because this soil property is independent of temperature (see Eq. [6]).

Figure 6 provides a comparison of water content measured independently using an EC-5 sensor at a depth of 2.48 cm bgs and Approaches 1 and 2. In Approach 1, the water content was calculated using Eq. [14] and simulated in Approach 2 using the multiphase heat and mass transfer model. Approach 1 follows the EC-5 sensor data closely until Day 3, where it begins to deviate greatly due to an increase in volumetric heat capacity. The volumetric heat capacity increased during this time period in response to the temperature rise measured by the HP probe (Fig. 4). The increase in water content from Days 3 to 5.4 and Days 7.8 to 10 are therefore the direct result of the observed increases in the maximum temperature change during these two time periods (Fig. 4). The numerical model (Approach 2) underestimated the water content compared with the EC-5 data during the entire course of the experiment. Through the adjustment of soil hydraulic properties and model-fitting parameters, we could have fit the experimental and modeling data (Approach 2). However, the goal of this work was not to fit experimental and modeling data but rather to investigate soil evaporation formulated with different approaches. Therefore, with this in mind, we did not introduce additional empirical fitting parameters, even though they are widely used in numerical modeling efforts in unsaturated soils. The underestimation of water content is consistent with that seen by Smits et al. (2011), who concluded that further work focusing on small-scale processes (experiments and models) is needed to better understand the fundamental physical processes at play so as to avoid using fitting parameters.

Sensible Heat Balance Calculation of Evaporation

The total evaporation rate calculated using Approaches 1 and 2 can be compared directly with the experimental water loss measurements made using the weighing scale (Fig. 7).

Evaporation rates were calculated in Approach 1 by applying the observed temperature gradients, heat inputs, and calculated thermal properties to the SHB (Eq. [1]). Similarly, the temperatures and thermal properties simulated using the multiphase, nonisothermal model were applied to the SHB in Approach 2. A total of four evaporation rates, one for each HP probe or soil layer, were determined using both Approach 1 and Approach 2. The sum of the evaporation rates yielded a total evaporation rate (Fig. 7) (Heitman et al., 2008a, 2008b; Sakai et al., 2011; Xiao et al., 2011, 2012). The dip visible in the experimental evaporation curve between Days 6 and 8 is the result of experimental error. The weighing scale’s alignment was temporarily compromised when the experimental setup was accidentally nudged while working in the vicinity of the experimental apparatus. The calculated rise in the evaporation rate (Approach 1) between Days 1.8 and 2.2 is the result of small temperature increases measured by two of the HP probes (not shown in Fig. 4). Therefore, the summation of the four soil layers’ evaporation rates yielded the observed bump. The slight decrease in evaporation rates between Days 3.8 and 4.1 (Fig. 7) can be similarly explained in terms of an observed temperature decrease.

Visual inspection of Fig. 8 shows that Approaches 1 and 2 are in good agreement with the measured water loss rates. In Approaches 1 and 2, cumulative evaporation was determined by integrating the respective evaporation curves with respect to time in accordance to Xiao et al. (2011, 2012). Statistical MIA analysis showed a 96.8% agreement between Approach 1 and the experimental total evaporation and 97.7% agreement between Approach 2 and the experimental total evaporation. Despite the discrepancies in the water loss measurements and the two approaches, the general behavior and magnitude of the three evaporation curves (Fig. 7) and total evaporation curves (Fig. 8) are similar. The three evaporation curves in Fig. 7 are marked by initially high evaporation rates that decrease rapidly and stabilize or plateau at later times. The curves presented in Fig. 8 display high water losses at the start of the experiment and very little water loss at later times.

The trends observed in the evaporation rate and water loss curves in Fig. 7 and 8 can be explained in terms of the two stages of evaporation discussed above. Stage 1 evaporation occurs from Days 0 to1.75 (experimental data set and Approach 2) and Days 0 to 2.2 (Approach 1); evaporation during these time periods was driven primarily by atmospheric demand. Days 1.75 and 2.2 mark the start of the transition period of evaporation to its vapor-diffusion limited Stage 2 for the experimental data set/Approach 2 and Approach 1, respectively. The rapidly falling evaporation rates during this stage are common for the uniform no. 30/40 Accusand used in the experimental portion of this study. Stage 2 evaporation was reached by Day 4 for Approach 1 and Day 6 for the experiment and Approach 2. The initially high evaporation rates observed during Stage 1 evaporation and the long transition period (2–4.25 d) are the result of the high hydraulic conductivity (0.108 cm s−1) and large displacement head (16.1 cm) of the no. 30/40 sand used in the present study (Sakai et al., 2011). The large hydraulic conductivity and displacement pressure allowed the soil surface to remain in hydraulic connection with the evaporation front.

Approach 1 underestimated the evaporation rate (Fig. 7) throughout the majority of the experiment, with the notable exception being late Stage 1 evaporation and early stage transition. This, in turn, yielded the observed underestimation of water loss shown in Fig. 8. Consistent with previous results, however, the overall underestimation was expected given the 1.5-mm undetectable zone at the soil surface in which the soil thermal properties are unknown (Heitman et al., 2008b; Sakai et al., 2011). As explained above, the deviation of Approach 1 from the experimental results from Day 1.8 to 2.2 was the result of the measured temperature increase by two of the sensors. This allowed a higher steady evaporation rate and therefore water loss to be maintained longer than expected. The short transition period (2 d) in comparison with the 4.25 d in the experimental data and Approach 2 can be explained in terms of the HP and SHB methods’ sensitivity to temperature (Sakai et al., 2011). The large temperature increase starting at approximately Day 2 therefore led to rapid soil drying and a fast transition to Stage 2 evaporation for Approach 1. As mentioned above, the goal of this work was not to fit the numerical multiphase heat and mass transfer model (Approach 2) to experimental results using fitting parameters. Rigorous parameter adjustment and incorporation of physical processes that may be of importance to the model would have allowed us to match Approach 2 and the experimental data with almost exact precision. The overestimation of Approach 2 with respect to water loss (Fig. 8) is therefore consistent with the results of Smits et al. (2011, 2012).


With the goal of verifying the combined HP and SHB method and its applicability to numerical models, we modified a fully coupled, nonisothermal, nonequilibrium, multiphase flow (heat, liquid water, and water vapor) model developed by Smits et al. (2011) to include the HP and SHB method. Results from the modified model were compared with an analytical HP and SHB model as well as the measured data. We developed and used a unique soil tank equipped with a network of different soil moisture, temperature, and relative humidity sensors that interfaced with a wind tunnel placed above. The results from the analytical and heat and mass transfer numerical simulations were compared with the experimental data.

The results demonstrate that both models can capture the soil thermal property behavior (thermal diffusivity, volumetric heat capacity, and thermal diffusivity) across a large range of water contents. However, agreement differed at water contents at and below residual saturation. The variation at low water contents was most likely due to a smaller ΔTmax value in addition to the sensor technology becoming very sensitive to slight changes under these conditions. Evaporation rates estimated using both modeling approaches were in close agreement with the experimentally measured evaporation rates. The multiphase heat and mass transfer model overestimated the evaporation rate, not capturing some of the physical processes at play. The analytical HP and SHB model underestimated the late part of Stage 1 evaporation due to an observed temperature increase between Days 1.5 and 2.2 by two of the HP probes. These findings support the use of the HP and SHB methods with a multiphase heat and mass transfer model to predict evaporation rates and cumulative water loss at fine spatial and temporal scales. The methodology tested and discussed in this study provides the basis for future experimental efforts to predict evaporation rates under heterogeneous and more complex boundary conditions conducted at larger scales. It provides an invaluable tool for calculating evaporation rates in situ without reliance on hydraulic properties and with limited reliance on fitting parameters.

This research was funded by the U.S. Army Research Office Award W911NF-04-1-0169 and the Engineering Research and Development Center (ERDC) and the National Security Science and Engineering Fellowship (NSEFF) Award no. 9559-10-1-0139. We acknowledge David Stepp from the Army Research Office and Stacy Howington, John Peters, and Matthew Farthing from the ERDC for support and technical contributions.

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Open Access