Estimating the effect of shallow groundwater on diurnal heat transport in a vadose zone

Abstract

The influence of shallow groundwater on the diurnal heat transport of the soil profile was analyzed using a soil sensor automatic monitoring system that continuously measures temperature and water content of soil profiles to simulate heat transport based on the Philip and de Vries (PDV) model. Three experiments were conducted to measure soil properties at depths of 5 cm, 10 cm, 20 cm, and 30 cm when groundwater tables reached 10 cm, 30 cm, and 60 cm (Experiments I, II, and III). Results show that both the soil temperature near shallow groundwater and the soil water content were effectively simulated by the PDV model. The root mean square errors of the temperature at depths of 5 cm, 10 cm, and 20 cm were 1.018°C, 0.909°C, and 0.255°C, respectively. The total heat flux generated the convergent and divergent planes in space-time fields with valley values of –161.5 W·m⁻² at 7:30 and –234.6 W·m⁻² at 11:10 in Experiments II and III, respectively. The diurnal heat transport of the saturated soil occurred in five stages, while that of saturated-unsaturated and unsaturated soil profiles occurred in four stages because high moisture content led to high thermal conductivity, which hastened the heat transport.

Keywords

heat transport / heat flux / soil temperature / water content / shallow groundwater table

Introduction

The heat transport in a vadose zone plays a significant role in the water and energy balance of a near surface environment in many agricultural applications. Philip and de Vries ( 1957) investigated thermal influence on unsaturated soils. Milly ( 1982) later extended the work to heterogeneous and hysteretic media. Theories that soil water content strongly affects heat transfer dynamics by changing soil thermal properties were later developed ( Jabro, 2009; Deb et al., 2011). Nassar and Horton ( 1997) indicated that latent heat transport is positive for vapor and negative for water evaporation. Moreover, the thermally induced vapor movement absorbs or releases energy through evaporation and condensation processes, resulting in rapid heat transfer in evaporation conditions ( Kane et al., 2001; Grifoll et al., 2005). Heat transport research has mainly focused on the climate-controlled energy balance of the soil–atmosphere interface and has thus ignored the lower boundary that is affected by groundwater ( Bittelli et al., 2008; Zeng et al., 2011; Schaefli et al., 2012). Fan et al. ( 2011) found that the atmospheric warming leads to a notable increase in deep soil temperature, with the increase being higher in wet sites than in dry sites because heat transport in the former results mainly from conduction. Saito and Šimunek ( 2009) assumed that heat transfer across the lower boundary occurs only through the conversion of liquid to vapor when the groundwater table is far below the domain of interest. However, the effect of shallow groundwater table on heat transfer across soil profiles has rarely been analyzed, although the link between groundwater table and surface temperature has been discussed with remote sensing data ( Corbari et al., 2015).

Heat transport in a vadose zone is typically studied under non-isothermal conditions based on numerical models. However, choosing the best model for estimating soil heat transport in highly complex conditions is very difficult. The movement of water and heat in soil are coupled. The Philip and de Vries (PDV) model ( Philip and de Vries, 1957) remains the main basis of soil–atmosphere continuum modeling ( Zeng et al., 2009). The contemporary version of the PDV model has been slightly modified by numerous studies ( Nassar and Horton, 1997; Cahill and Parlange, 1998; Saito et al., 2006). Cass et al. ( 1984) measured apparent soil thermal conductivity (λ) under various air pressures and developed a method for separating the contributions of conduction and the latent heat transfer of vapor to overall soil heat flow. Lu et al. ( 2011) treated conduction heat transfer as the apparent soil thermal conductivity associated with infinite atmospheric pressure. The goal of this paper is to justify if the PDV model can be properly applied to the soil heat transport above shallow groundwater tables.

Very few studies have simultaneously and continuously demonstrated and evaluated soil heat dynamics in time and space, although the theory of heat transport in boreal soil is widely developed. The most common way to address this issue is either to assess the time-series information at specific depths ( Wang, 2002; Novak, 2010; Laredj et al., 2012) or to analyze profile temperature and heat flux information at specific times ( Gran et al., 2011). In this study, both the spatial and temporal variations of heat flux are presented. The space-time information represents a two dimension field while a dependent heat flux represents a third dimension. These factors are directly established through an interpolation and smoothing procedure. The primary objective is to measure soil water content and temperature changes and to subsequently produce temporal and spatial information on the heat transport based on the PDV model. The heat transfer dynamics of soil profiles above shallow groundwater tables are analyzed in this study as well.

Material and methods

Experiments

Three experiments (Experiments I, II, and III) involving three groundwater table depths (10 cm, 30 cm, and 60 cm) were conducted under non-isothermal evaporation conditions. The soil columns were divided into three representative types; namely, saturated, saturated-unsaturated, and absolutely unsaturated zones, according to the capillary water rise height.

Soil samples were obtained from the saline wasteland of the China–Singapore eco-city, Tianjin, China (116°43′E−118°194′E, 38°34′N−40°15′N). The samples were air-dried, crushed, and separated by a 10 mm sieve. Before packing the soil samples into the soil box measuring 100 cm ×50 cm×100 cm (Fig. 1), fine sand layers with a thickness of 20 cm were packed into the bottom of the soil box to represent groundwater aquifers. A mesh was placed in the sand-layer boundary to avoid flow erosion. Two side water boxes measuring 25 cm×50 cm×100 cm were used as an inlet and an outlet to adjust the groundwater table. Subsequently, the soil samples were packed into 10 cm layers to achieve a natural bulk density of approximately 1.53 g·cm⁻³. Each soil layer was lightly raked prior to packing the next layer to minimize the discontinuities between layers. The total heights of the soil columns were 10 cm, 30 cm, and 60 cm in Experiments I, II, and III, respectively.

Soil probes, including Hydra probes (Steven, USA) and SWT4R tension meters (Delta-T, England), were inserted into the soil samples at depths of 5 cm, 10 cm, 20 cm, and 30 cm to measure temperature and water content changes. The inlet and outlet of the water valves were opened to maintain a constant level of groundwater table. The evaporation simulator was operated from 8:30 to 17:30 daily to control fluctuation of the air temperature from (27±1)°C to (35±1)°C, similar to summer temperatures at mid-latitudes ( Zeng et al., 2009). Temperature, water content, and water potential data were recorded every 10 min by a data logger (CR1000, Campell, USA) and stored in a disk. The relative humidity and air temperature were measured hourly by a psychrometer (wi270740, China). Each experiment was conducted for five days. This procedure was duplicated to test the repeatability of Experiments I, II, and III.

Model description

Theory

Heat transport

The HYDRUS-1D code, which refers to coupled water, water vapor, and heat transport in soil, was applied in the simulation of soil heat fluxes. Total heat flux density, q_h (J·m⁻²·s⁻¹), is defined as the sum of the conduction of sensible heat as described by Fourier’s law—sensible heat by convection of liquid water and water vapor, and latent heat by vapor flow ( de Vries, 1958), given as follows:

(1) q_h=-{\lambda }_c(\theta ){\displaystyle \frac{\partial T}{\partial z}}+C_{W}T{q}_L+C_{v}T{q}_v+L{q}_v,

where {\lambda }_c(\theta ) is the thermal conductivity of soil (W·m⁻¹·K⁻¹), T is a given temperature (K), z is the soil depth (m), C_W and C_v are the volumetric heat capacities (J·m⁻³·K⁻¹) of the liquid water and water vapor, respectively ( de Vries, 1963), q_L and q_v are the flux densities of liquid water and water vapor (m·s⁻¹), respectively, and L is the latent heat of vaporization (~2.45×10⁶ J·kg⁻¹).

Liquid water and water vapor transport

When root water uptake is not considered in a soil water flux simulation, the governing equation for the one-dimensional flow of liquid water and water vapor in a variably saturated rigid porous medium is given by the following mass conservation equation ( Saito et al., 2006):

(2) {\displaystyle \frac{\partial \theta }{\partial t}}=-{\displaystyle \frac{\partial q_L}{\partial z}}-{\displaystyle \frac{\partial q_v}{\partial z}},

where \theta is the volumetric water content (m³·m⁻³), q_L and q_v are the flux densities of liquid water and water vapor (m·s⁻¹), respectively, t is time (s), and z is the spatial coordinate positive upward (m).

The flux density of liquid water, q_L, is described using the following modified Darcy law ( Philip and de Vries, 1957):

(3) q_L=q_{Lh}+q_{LT}=-K_{Lh}\left({\displaystyle \frac{\partial h}{\partial z}}+\mathrm{1}\right)-K_{LT}{\displaystyle \frac{\partial T}{\partial z}},

where q_{Lh} and q_{LT} are the isothermal and thermal liquid water flux densities (m·s⁻¹), respectively, and K_{\mathrm{Lh}} (m·s⁻¹) and K_{LT} (m²·K⁻¹·s⁻¹) are the isothermal and thermal hydraulic conductivities, respectively, for liquid phase fluxes resulting from gradients in h and T.

The flux density of water vapor, q_v , can also be classified into the isothermal, q_{vh} , and the thermal, q_{vT} , vapor flux densities (m·s⁻¹) as follows:

(4) q_v=q_{vh}+q_{vT}=-K_{vh}{\displaystyle \frac{\partial h}{\partial z}}-K_{vT}{\displaystyle \frac{\partial T}{\partial z}},

where K_{vh} (m·s⁻¹) and K_{vT} (m²·K⁻¹·s⁻¹) are the isothermal and thermal vapor hydraulic conductivities, respectively.

The water retention curve (WRC) is one of the most fundamental hydraulic characteristics of solving the flow equation of water in soil. The soil water retention equation is as follows ( van Genuchten, 1980):

(5) \theta (h)=\{\begin{array}{ll}{\theta }_r+{\displaystyle \frac{{\theta }_s-{\theta }_r}{{[\mathrm{1}+{\left|\alpha h\right|}^n]}^m}}\hfill & h\le \mathrm{0}\hfill \\ {\theta }_s\hfill & h>\mathrm{0}\hfill \end{array}

where {\theta }_r and {\theta }_s are the residual and saturated water contents, respectively (cm³·cm⁻³), a (>0, in cm⁻¹) is related to the inverse of the air-entry pressure, and n (>1) is a measure of the pore size distribution affecting the slope of the retention function (m= 1‒1/n).

The pore size distribution model of Mualem ( 1976) was used to predict the isothermal unsaturated hydraulic conductivity function, K_Lh(h), from the saturated hydraulic conductivity and van Genuchten’s ( 1980) model of soil WRC.

(6) K_{Lh}(h)=K_s{S}_e^l{[\mathrm{1}-{(\mathrm{1}-S_e^{\mathrm{1}/m})}^m]}^{\mathrm{2}},

where K_s is the saturated hydraulic conductivity (m·s⁻¹), S_e is the effective saturation (unitless), and l and m are empirical parameters (unitless). Parameter l was set to 0.5 as suggested by Mualem ( 1976). Parameter m was determined by fitting van Genuchten’s analytical model.

(7) S_e={[\mathrm{1}+{(\alpha h)}^n]}^m.

The thermal hydraulic conductivity function, K_LT, in Eq. (3) is defined as follows ( Noborio et al., 1996):

(8) K_{LT}=K_{Lh}\left(h{G}_{WT}{\displaystyle \frac{\mathrm{1}}{{\gamma }_{\mathrm{0}}}}{\displaystyle \frac{d\gamma }{dT}}\right),

where G_WT is the gain factor (unitless) that quantifies the temperature dependence of the soil WRC ( Nimmo and Miller, 1986), \gamma is the surface tension of soil water (J·m⁻²), and {\gamma }_0 is the surface tension at 25°C (71.89 g·s⁻²). The temperature dependence of \gamma is as follows:

(9) \gamma =75.6-0.1425T-2.38\times {10}^{-4}{T}^2.

The isothermal, K_{vh} , and thermal, K_{vT} , vapor hydraulic conductivities are described as follows:

(10) K_{vh}={\displaystyle \frac{\alpha aD}{{\rho }_w}}{\rho }^{*}{\displaystyle \frac{Mg}{RT}}{H}_r,

(11) K_{vT}={\displaystyle \frac{\alpha aD}{{\rho }_w}}\eta H_r{\displaystyle \frac{d{\rho }^{*}}{dT}},

where {\rho }_w is the pure water density (kg·m⁻³).

Cass et al. ( 1984) developed an empirical model to estimate η from soil water content and porosity:

(12) \eta =9.5+3{\displaystyle \frac{\theta }{{\theta }_s}}-\mathrm{8.5}\exp \left\{-{\left[\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2.6}}{\sqrt{f_c}}}\right){\displaystyle \frac{\theta }{{\theta }_s}}\right]}^{\mathrm{4}}\right\},

where f_c is the mass fraction of clay in soil (unitless).

Meteorological variable generation

The air temperature and relative humidity over the soil surface did not always match the standard meteorological data at the time intervals of interest. Thus, the diurnal changes of both factors must be calculated from the available daily average values using meteorological models ( Philip and de Vries, 1957). The continuous values for air temperature, T_a, can be obtained from the daily maximum and minimum air temperatures using a trigonometric function with a period of 24 h.

(13) T_a=\overline{T}+A_t\cos \left[\mathrm{2}\pi {\displaystyle \frac{(t-t_{\max })}{\mathrm{24}}}\right],

where \bar{T} is the average daily temperature (°C), A_t is the daily amplitude of the cosine wave (°C) calculated from the difference between the daily maximum and minimum temperatures, t is a local time during the day (h), and t_max is the time of highest temperature.

Relative humidity can be assumed to demonstrate the same cyclic behavior throughout the day due to the variance in typical diurnal air temperature ( Gregory et al., 1994). Therefore, a trigonometric function can be used to calculate the continuous values of relative humidity based on the daily maximum and minimum values within a 24 h period.

(14) H_r=\bar{H}+A_r\cos \left[\mathrm{2}\pi {\displaystyle \frac{(t-t_{\max })}{\mathrm{24}}}\right],

where \bar{H} is the average daily relative humidity (unitless), A_r is the daily amplitude of the cosine wave calculated from the difference between the maximum and minimum relative humidity values, t is the local time within the day, and t_max is the hour of the day when the relative humidity is at its maximum.

Surface Energy Balance

Surface precipitation, irrigation, evaporation, and heat fluxes were used as boundary conditions for liquid water and water vapor flow and for heat transport through unsaturated soils. Evaporation and heat fluxes can be calculated from the major components of the vertical energy balance at the soil surface ( Saito et al., 2006; Deb et al., 2011).

(15) R_n-H-LE-G=\mathrm{0},

where R_n is the net radiation (positive downward) (M·J·m⁻²·d⁻¹); H is the sensible heat flux density (positive upward) (M·J·m⁻²·d⁻¹); LE is the latent flux density (positive upward) (M·J·m⁻²·d⁻¹), where L and E are the latent heat (M·J·kg⁻¹) and evaporation rate (kg·m⁻²·d⁻¹), respectively, and G is the surface heat flux density (positive upward) (M·J·m⁻²·d⁻¹).

Net radiation is defined as follows ( Brutsaert, 1982):

(16) R_n=R_{ns}+R_{nl}=(\mathrm{1}-a)S_t+({\epsilon }_s{R}_{ld\downarrow }-R_{lu\uparrow }),

where R_{ns} is the net shortwave radiation (W·m⁻²), R_{nl} is the net longwave radiation (W·m⁻²), a is the surface albedo (unitless), S_t is the incoming (global) shortwave solar radiation (W·m⁻²), {\epsilon }_s is the soil surface emissivity (unitless) representing the reflection of the longwave radiation at the soil surface, R_{ld\downarrow } is the incoming (thermal) longwave radiation at the soil surface (downward flux) (W·m⁻²) as emitted by the atmosphere and cloud cover, and R_{lu\uparrow } is the sum of the outgoing (thermal) longwave radiation emitted from the surface (vegetation and soil) into the atmosphere (W·m⁻²).

Van Bavel and Hillel ( 1976) proposed the following simple formulae to calculate the surface albedo from surface wetness:

(17)

where {\theta }_{\mathrm{top}} is the water content at the surface.

The incoming shortwave solar radiation S_t (W·m⁻²) is calculated by considering the position of the sun as follows ( Campbell, 1985):

(18) S_t(t)=\max (G_{sc}{T}_t\sin e,\mathrm{0}),

where G_sc is the solar constant (1,360 W·m⁻²), and T_t (unitless) is defined as follows:

(19) T_t={\displaystyle \frac{S_{tm}}{R_a}},

where S_tm is the daily global solar radiation (W·m⁻²), and R_a is the daily potential global (extraterrestrial) radiation (W·m⁻²).

The last term of Eq. (18), \sin e , is the solar elevation angle given by ( Monteith and Unsworth, 1990)

(20) \sin e=\sin \phi \sin \delta +\cos \phi \cos \delta \cos {\displaystyle \frac{\mathrm{2}\pi }{\mathrm{24}}}(t-t_{\mathrm{0}}),

where \phi is the local geographic latitude, and \delta is the sun declination of the day.

The net longwave radiation can be rewritten as follows based on the Stefan–Boltzmann law:

(21) R_{nl}={\epsilon }_{s}R{}_{ld\downarrow }-R_{lu\uparrow }={\epsilon }_s{\epsilon }_a\sigma T_a^{\mathrm{4}}-{\epsilon }_s\sigma T_s^{\mathrm{4}},

where T_a is the air temperature (K), T_s is the soil surface temperature (K), {\epsilon }_a is the atmospheric emissivity (unitless), and σ is the Stefan-Boltzman constant (5.67×10⁻⁸ J·s⁻¹·m⁻²·K⁻⁴).

The emissivity of bare soil can be expressed as a function of the volumetric water content.

(22) {\epsilon }_s=\min (\mathrm{0.90}+\mathrm{0.18}{\theta }_{\mathrm{top}};\mathrm{1.0}).

Following Idso (1981), the atmospheric emissivity can be calculated as follows:

(23) {\epsilon }_a=\mathrm{0.70}+\mathrm{5.95}\times {\mathrm{10}}^{-\mathrm{5}}{e}_a\exp (\mathrm{1500}/T_a),

where e_a is the atmospheric vapor pressure (kPa) expressed as a function of T_a.

(24) e_a=\mathrm{0.611}\exp \left[{\displaystyle \frac{\mathrm{17.27}(T_a-\mathrm{273.15})}{T_a-\mathrm{35.85}}}\right]H_r.

The sensible heat flux, H, can be simply defined as follows ( van Bavel and Hillel, 1976):

(25) H=C_a{\displaystyle \frac{T_s-T_a}{r_H}},

where C_a is the volumetric heat capacity of air (J·m⁻³·K⁻¹), and r_H is the aerodynamic resistance to heat transfer (s·m⁻¹).

The computation of evaporation, E, is defined as follows ( Camillo and Gurney, 1986):

(26) E={\displaystyle \frac{{\rho }_{vs}-{\rho }_{va}}{r_v+r_s}},

where r_v is the aerodynamic resistance for water vapor transfer (s·m⁻¹), r_s is the soil surface resistance for water vapor transfer (s·m⁻¹), {\rho }_{va} is the atmosphere vapor concentration (mol·mol⁻¹), and {\rho }_{vs} is the vapor concentration of the soil surface (mol·mol⁻¹).

The aerodynamic resistance for heat (r_H) and vapor (r_v) transport depends on surface roughness properties and wind speed ( Campbell, 1977).

(27) r_H=r_v={\displaystyle \frac{\mathrm{1}}{u^{*}k}}\left[\ln \left({\displaystyle \frac{z_{ref}-d+z_H}{z_H}}\right)+{\phi }_H\right],

where u^* is the friction velocity (m·s⁻¹), k is the von Karman’s constant equal to 0.41, z_ref is the reference height of the temperature measurement (m), z_H is the surface roughness parameter for the heat flux (m), d is the zero plane displacement (m), and {\phi }_H is the atmospheric stability correction factor for the heat flux. Friction velocity is defined as follows:

(28) u^{*}=uk\left[\ln \left({\displaystyle \frac{z_{ref}-d-z_M}{z_M}}\right)+{\phi }_M\right],

where u is wind speed (m·s⁻¹) at height z_ref, z_M is the surface roughness for the momentum flux, and {\phi }_M is the atmospheric stability correction factor for momentum flux. The typical value of soil surface roughness is 0.01 m, which was employed when calculating r_v ( Campbell and Norman, 1998). The stability function was used to calculate {\phi }_M and {\phi }_H , ( Webb, 1970).

(29) {\phi }_H={\phi }_M=-\xi f,

where \xi =(z-z_M-d)/\mathrm{\Lambda } with \mathrm{\Lambda } being the Monin-Obukhov’s stability parameter, and f, a constant equal to 4.7. The stability parameter is computed as follows:

(30) \mathrm{\Lambda }={\displaystyle \frac{u^{*\mathrm{3}}{C}_h{T}_a}{kgH}},

where C_h is the volumetric heat of air (1,200 J·m⁻³·C⁻¹ at 20°C, at sea level), and g is the gravitational constant (9.81 m·s⁻¹).

Initial and boundary conditions

The soil profile was 60 cm deep, and a spatial discretization of 0.5 cm led to 121 nodes across the flow domain. The observation nodes for the temperature and matric potential were located at depths of 5 cm, 10 cm, 20 cm, and 30 cm. Calculation was performed in a 120 h period in August 2013. Soil profile was assumed to be vertically homogenous. The initial conditions for temperature and pressure head were based on a linear interpolation of the measured values at the site.

The upper boundary conditions for heat transport were determined by the surface energy balance equation. HYDRUS-1D (version 4.16) uses continuous meteorological data (Eq. (16)) to obtain surface heat flux, G, which is a known heat flux boundary condition. Surface evaporation is calculated as a surface boundary condition for soil moisture transport ( Deb et al., 2011). In this study, the lower boundary condition for heat transport was temperature BC when the groundwater table was 10 cm, while that for Experiments II and III was a zero temperature gradient. The water flow at the lower boundary maintained the constant pressure head because of the shallow water tables.

HYDRUS-1D uses measured soil matric potential and temperature data at all depths to inversely estimate water and heat transport parameters. In this study, some of the hydraulic and thermal parameters were also optimized to consider the inherent soil variability within the experiment conditions, although soil water retention, hydraulic properties, and soil thermal properties were obtained. The parameter values used for the numerical applications are shown in Table 1.

Results and discussion

Meteorological measurements and calculations

The diurnal changes of the air temperature and relative humidity over the soil surface were calculated through trigonometric functions (Eqs. (13) and (14)). Figure 2 shows the good agreement between the simulation and the observation of the temperature and relative humidity in a time series.

According to the monitored data (Fig. 2(a)), the peak temperature of (36±1)°C occurred at 14:00, and the valley value of (27±1)°C was recorded at 2:00. These values are similar to the summer temperature changes in Tianjin, China. The diurnal cycles of atmospheric relative humidity can be calculated by using a cosine function (Eq. (14)). Gregory et al. ( 1994) pointed out that relative humidity in Lubbock, Texas peaks between 5:00 and 6:00. In the present study, the peak relative humidity was recorded at 5:00 (Fig. 2(b)).

Soil temperature and water content

The observed soil temperatures and water content were compared with those calculated by the HYDRUS-1D. The simulation’s goodness of fit was quantified with the following root mean square error (RMSE) measure:

\mathrm{RMSE}=\sqrt{{\displaystyle \sum _{i=\mathrm{1}}^{N_w}{(M_i-C_i)}^{\mathrm{2}}}/N_w,}

where N_w is the number of measurements, and M_i and C_i are the measurements and calculations, respectively.

The diurnal changes of soil temperature are mainly influenced by evaporation and focus on the upper soil depth. Thus, the observed results were compared with the simulation at the upper layers (depths of 5 cm, 10 cm, and 20 cm) when the groundwater table was 60 cm. The diurnal temperature cycles were reasonably predicted by the PDV model (Fig. 3). The RMSEs of the temperatures at depths of 5 cm, 10 cm, and 20 cm were 1.018°C, 0.909°C, and 0.255°C, respectively. The results indicate that the soil temperature near the water table was effectively simulated by the PDV model. The experimental and simulated temperatures displayed a typical diurnal sinusoidal behavior. The amplitude decreased and the valley and peak times were delayed with depth. These conditions are similar to those in field trials conducted in other studies ( Saito et al., 2006; Bittelli et al., 2008; Deb et al., 2011).

Groundwater recharge significantly affects soil moisture redistribution through capillary rise, especially for moisture at deeper depths. Therefore, water content changes observed in the deep layers (depths of 5 cm, 20 cm, and 30 cm) were measured in this study. A comparison of the measurement and simulation of water content is provided in Fig. 4.

The general trends of the simulation agreed with the measured values at soil depths below 20 cm. The RMSEs of the water content at depths of 20 cm and 30 cm were 0.004 m³·m⁻³ and 0.002 m³·m⁻³, respectively. For the observed and simulated water content at a depth of 5 cm, a discrepancy was observed in the maximum deviation of 1.4% at approximately the 50^th hour. Other studies have also reported such discrepancy levels, especially under low moisture and high temperature conditions. For example, Grifoll et al. ( 2005) used a general model formulation to simulate water and energy transport, and found that the average discrepancy between measured and simulated moisture content is 27.5% near the dry surface (0 cm to 1.27 cm). Nassar and Horton ( 1999) pointed out the good agreement between predicted and measured values given mean soil temperatures of 10°C and 27°C. This finding reveals that high soil water content equates to good simulation results obtained from the PDV model with non-isothermal evaporation.

Surface energy flux

Surface heat fluxes are directly influenced by soil temperature variations and other factors ( Chen et al., 2014). The diurnal surface energy fluxes calculated based on Eqs. (15) and (26) are depicted in Fig. 5. The calculation period was set to the last 24 hours in the experiment to minimize the anthropogenic effect of the initial conditions. The surface heat flux at daytime increased and resulted in upward sensible heat fluxes (positive values). At night, the surface heat flux generally decreased (negative values) as the soil surface cooled (Fig. 5(a)). As explained by Bittelli et al. ( 2008), latent heat flux is negative during the day because of water evaporation from the soil surface. When the groundwater table increased from 10 cm to 60 cm in the present work, the peak heat flux decreased from 148.05 W·m⁻² to 109.11 W·m⁻², and the soil heat contribution for the vapor flux decreased from 2.59 W·m⁻² to 0.01 W·m⁻² (Fig. 5(b)). These results, which match that of Saito et al. ( 2006), confirm that the heat flux and latent heat flux increased significantly as the soil moisture increased.

Soil heat flux

The study of heat flux variation is essential to gain a better understanding of the mechanism of heat transfer in soil profiles. In this study, the heat flux comprising conduction and latent heat fluxes was calculated by using a function of soil temperature and moisture content (Eq. (1)).

Figure 6 shows the changes in the diurnal heat flux in the three typical soil profiles (saturated, saturated-unsaturated, and unsaturated profiles). The heat flux variation was analyzed in the last 24 hours of the experiment not only to minimize the impact of the initial conditions, but to ensure a high goodness of fit in the model. The heat flux of the saturated soil profile at a groundwater table depth above 10 cm clearly varied with the amplitude of 99.92 W·m⁻² at the lower boundary (Fig. 6(a)). The changes in the heat flux of the saturated-unsaturated and absolutely unsaturated soil profiles were primarily detected at the upper 20 cm and 15 cm depths, respectively (Figs. 6(b) and 6(c)). The result is supported by the finding of Deb et al. ( 2011) who found that the dramatic fluctuation of heat flux is observed at the upper 10 cm depth and remains constant with depth. However, their work is inconclusive as they failed to consider the influence of groundwater, which is significant in soil heat transfer. In this paper, the result indicates that the depth of heat flux fluctuation decreases with increasing groundwater table depth due to the decreasing heat conductivity.

An interpolation and smoothing program (Surfer 11) was used to create a space-time field for interpreting the diurnal pattern of heat flux. A dependent specific variation of heat flux was employed as the third dimension. The contour charts of the total heat flux, including the conduction and latent heat fluxes (Fig. 7), show the continuous variations above the three groundwater table depths (10 cm, 30 cm, and 60 cm). When the groundwater table depth was 60 cm, the heat flux variation below 30 cm was too small to be ignored (Fig. 7(c)). Total heat flux is positive when energy is transferred into soil and is negative in the upward direction.

When the groundwater table was at a 10 cm depth before dawn (0:00−9:00), the upward heat flux of the soil surface varied from −62.16 W·m⁻² to 0 W·m⁻², with the valley time being 6:00. Heat transport first increased and then decreased from the inferior to the cooler soil surface (Fig. 7(a)). At daytime (10:00 to 19:00), the temperature of the solar-heated soil surface was higher than that of the interior. This condition resulted in a downward heat transport. The surface heat flux increased from 0 W·m⁻² to 147.9 W·m⁻², and the maximum flux was recorded at 14:00. After 20:00, the surface heat flux demonstrated an upward trend until 24:00. As shown in Fig. 7(a), two zero heat flux planes were observed in the soil profile and were classified into convergent and divergent planes based on the definitions of zero heat flux planes.

At groundwater table depths of 30 cm and 60 cm (Figs. 7(b) and 7(c)), another divergent zero heat flux plane was observed below the 4 cm depth from 0:00 to 10:00. Obvious ellipses were recorded at the 5 cm depth between the planes of divergence and convergence. The valley heat fluxes of the ellipses were −161.5 W·m⁻² at 7:30 and −234.6 W·m⁻² at 11:10 for Experiments II and III, respectively. The amplitude of the heat flux increased with the delay in the valley time and with the increase in groundwater table depth. This result indicates that the heat transport affected by the soil temperature and vapor transport was sensitive to the lower boundary condition ( Zeller and Nikolov, 2000; Prunty and Bell, 2005). Furthermore, the contours in the upper 10 cm depth were the most intensive. Therefore, the fluctuation of the heat flux was strongest near the surface. This result is consistent with the view of Kane et al. ( 2001).

The latent heat flux in soil can be represented by the last term of Eq. (1). In this study, the latent heat flux was calculated based on Eqs. (1), (4), (11), and (12). Figure 8 shows that the latent heat flux clearly varied in the soil profiles and that its contribution to the total heat flux should not be neglected. When groundwater table depths were 30 cm and 60 cm, the peak values of the latent heat flux in the soil profiles reached 32.6 W·m⁻² and 53.3 W·m⁻² respectively. The latent heat flux in soil was smaller than that on the surface at a groundwater table depth of 10 cm. This result can be attributed to the saturated soil having little water vapor and to the occurrence of main water vapor variation on the soil surface. Bittelli et al. ( 2008) reported that the average daytime soil heat contribution from vapor flux at a depth of 2 cm is in the order of 20 W·m⁻² (ranging from 2 W·m⁻² to 63 W·m⁻²).

Latent heat flux plays a significant role in the soil energy budget, and is mainly affected by soil moisture content ( Fan et al., 2011). Saito et al. ( 2006) found that latent heat fluxes significantly increased after irrigation. The reason is that the high soil moisture content increases thermal conductivity, which could increase latent heat flux. This result is very interesting and is worthy of further study.

Heat transport dynamics

As indicated by the spatial-temporal distributions of the heat flux, the diurnal heat transport occurred in more than four stages at shallow groundwater table depths above 10 cm, 30 cm, and 60 cm (Fig. 9). Given the high moisture content of the saturated soil, the thermal conductivity increased, and the heat transfer occurred more quickly than it did in the unsaturated soil.

Figure 9(a) shows that the heat transport in the saturated soil profile occurred in five stages. The first stage was from 0:00 to 9:00. The heat migrated upward from the bottom to the cool soil surface affected by the temperature gradient. During the second stage (9:00 to 10:00), solar radiation increased the energy of the soil surface to generate the downward heat flux, which gathered with the upward heat flux to generate the zero convergent plane for the soil profile. With the increase in the downward temperature gradient at the third stage (10:00 to 16:00), the heat transport demonstrated a downward trend throughout the entire soil profile. The fourth stage started from 16:00 to 20:00. At this point the upward heat flux occurred at 5 cm to generate the upper convergent plane and lower divergent plane as a result of the spatial heterogeneity of the soil. These conditions led to the positive energy and negative energy that were first cancelled out at the 5 cm depth. The heat transition from 20:00 to 24:00 was analyzed as the fifth stage. The depth of the divergent plane increased from 8 cm to 10 cm.

The heat transfer process in the saturated-unsaturated soil is shown in Fig. 9(b). At the first stage (0:00 to 6:00), the divergent plane range was observed at depths of 3.5 cm to 10 cm, above which the heat migrated to the soil surface. At the second stage (6:00 to 10:00), the convergent plane remained constant at a depth of 10 cm, beneath which the divergent plane extended down to a depth of 30 cm. The convergent plane appeared from the surface to an 8 cm depth during the third stage (10:00 to 19:00), thereby indicating that the downward heat transfer was dominant in the upper soil layer. At the fourth stage (19:00 to 24:00), the upward heat flux gradually increased to generate two plane types at depths of 0 cm to 3 cm and 8 cm to 12 cm.

The heat transport in the unsaturated soil profile occurred in four stages (Fig. 9(c)). In the first stage (0:00 to 9:00), the heat transport demonstrated an upward and then a downward trend from the divergent plane at depths of 3 cm to 20 cm, which is similar to that of the saturated-unsaturated soil profile. Meanwhile, convergent and divergent planes were observed in the second stage (9:00 to 12:00). In the third stage (12:00 to 19:00), heat energy gathered at the convergent plane as the upward heat transport decreased. The depth of the convergent plane increased from the surface to 2.5 cm at daytime. In the fourth stage, both convergent and divergent planes existed in the soil profile. The depth of the convergent plane remained constant at 4 cm while that of the divergent plane increased from 4 cm to 11 cm.

Conclusions

As an important factor in the overall water and energy balance of the near surface environment, soil heat transport affected by shallow groundwater under non-isothermal evaporation is seldom studied. In this work, soil temperature and water content were measured to produce temporal and spatial information on heat transport based on the PDV model. The heat transfer dynamics of the soil profiles above shallow groundwater tables were analyzed. The simulations of soil temperature and water content agreed with the measured values under non-isothermal evaporation. Soil temperature and water content near the shallow groundwater table was effectively simulated by the PDV model. Furthermore, the space-time fields of the diurnal pattern of heat fluxes were created. Convergent and divergent planes of heat flux existed in the space-time fields, and clear ellipses occurred between the two planes. The diurnal heat transport of saturated soil underwent five stages while that of saturated-unsaturated and unsaturated soil underwent four stages because the high moisture content of the saturated soil led to high thermal conductivity, which hastened heat transport. The PDV model effectively simulated the heat transport of saturated soil, especially with the shallow groundwater table as the lower boundary. However, additional descriptive and experimental work is needed to further understand the response of heat transport to other conditions, such as soil types and vegetation cover.

References

[1] Bittelli M, Ventura F, Campbell G S, Snyder R L, Gallegati F, Pisa P R (2008). Coupling of heat, water vapor, and liquid water fluxes to compute evaporation in bare soils. J Hydrol (Amst), 362(3−4): 191–205

[2] Brutsaert W (1982). Evaporation into the atmosphere: theory, history, and applications. Dordrecht: D. Reidel Publ.

[3] Cahill A T, Parlange M B (1998). On water vapor transport in field soils. Water Resour Res, 34(4): 731–739

[4] Camillo P J, Gurney R J (1986). A resistance parameter for bare-soil evaporation models. Soil Sci, 141(2): 95–105

[5] Campbell G S (1977). An Introduction to Environmental Biophysics. New York: Springer-Verlag

[6] Campbell G S (1985). Soil Physics with Basic, Transport Models for Soil-Plant Systems. Amsterdam: Elsevier

[7] Campbell G S, Norman J M (1998). An Introduction to Environmental Biophysics (2nd edition). New York: Springer-Verlag

[8] Cass A, Campbell G S, Jones T L (1984). Enhancement of thermal water vapor diffusion in soil. Soil Sci Soc Am J, 48(1): 25–32

[9] Chen S H, Su H B, Zhan J Y (2014). Estimating the impact of land use change on surface energy partition based on the Noah model. Front Earth Sci, 8(1): 18–31

[10] Corbari C, Mancini M, Li J, Su Z (2015). Can satellite land surface temperature data be used similarly to ground discharge measurements for distributed hydrological model calibration? Hydrological Sciences Journal, 60(2): 202–217

[11] de Vries D A (1958). Simultaneous transfer of heat and moisture in porous media. Trans Am Geophys Union, 39(5): 909–916

[12] de Vries D A (1963). The thermal properties of soils. In: van Wijk W R, ed. Physica of Plant Environment. Amsterdam: North-Holland Pub. Co, 210–235

[13] Deb S K, Shukla M K, Sharma P, Mexal J G (2011). Coupled liquid water, water vapor, and heat transport simulations in an unsaturated zone of a sandy loam field. Soil Sci, 176(8): 387–398

[14] Fan Z, Nefa J C, Harden J W, Zhang T, Veldhuis H, Czimczik C I, Winston G C, O'Donnell J A (2011). Water and heat transport in boreal soils: implications for soil response to climate change. Sci Total Environ, 409(10): 1836–1842

[15] Gran M, Carrera J, Massana J, Saaltink M W, Olivella S, Ayora C, Lloret A (2011). Dynamics of water vapor flux and water separation processes during evaporation from a salty dry soil. J Hydrol (Amst), 396(3−4): 215–220

[16] Gregory J M, Peterson R E, Lee J A, Wilson G R (1994). Modeling wind and relative humidity effects on air quality. In: International Conferences on Air Pollution from Agricultural Operations. Midwest Plan Service Ames IA, 183–190

[17] Grifoll J, Gastó J M, Cohen Y (2005). Non-isothermal soil water transport and evaporation. Adv Water Resour, 28(11): 1254–1266

[18] Idso S B (1981). A set of equations for full spectrum and 8- to 14-μm and 10.5- to 12.5-μm thermal radiation from cloudless skies. Water Resour Res, 17(2): 295–304

[19] Jabro J D (2009). Water vapor diffusion through soil as affected by temperature and aggregate size. Transp Porous Media, 77(3): 417–428

[20] Kane D L, Hinkel K M, Goering D J, Hinzman L, Outcalt S I (2001). Non-conductive heat transfer associated with frozen soils. Global Planet Change, 29(3−4): 275–292

[21] Laredj N, Missoum H, Bendani K, Maliki M (2012). A coupled model for heating and hydratation in unsaturated clays. Arab J Geosci, 5(5): 935–942

[22] Lu S, Ren T, Yu Z, Horton R (2011). A method to estimate the water vapour enhancement factor in soil. Eur J Soil Sci, 62(4): 498–504

[23] Milly P C (1982). Moisture and heat transport in hysteretic, inhomogeneous porous media: a matric head-based formulation and a numerical model. Water Resour Res, 18(3): 489–498

[24] Monteith J L, Unsworth M H (1990). Principles of Environmental Physics. London: Edward Arnold

[25] Mualem Y (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour Res, 12(3): 513–522

[26] Nassar I N, Horton R (1997). Heat, water, and solute transfer in unsaturated porous media: I – theory development and transport coefficient evaluation. Transp Porous Media, 27(1): 17–38

[27] Nassar I N, Horton R (1999). Salinity and compaction effects on soil water evaporation and water and solute distributions. Soil Sci Soc Am J, 63(4): 752–758

[28] Nimmo J R, Miller E E (1986). The temperature dependence of isothermal moisture vs. potential characteristics of soils. Soil Sci Soc Am J, 50(5): 1105–1113

[29] Noborio K, McInnes K J, Heilman J L (1996). Two-dimensional model for water, heat, and solute transport in furrow-irrigated soil: I. Theory. Soil Sci Soc Am J, 60(4): 1001–1009

[30] Novak M D (2010). Dynamics of the near-surface evaporation zone and corresponding effects on the surface energy balance of a drying bare soil. Agric Meteorol, 150(10): 1358–1365

[31] Philip J R, de Vries V D (1957). Moisture movement in porous materials under temperature gradient. Trans Am Geophys Union, 38(2): 222–232

[32] Prunty L, Bell J (2005). Soil temperature change over time during infiltration. Soil Sci Soc Am J, 69(3): 766–775

[33] Saito H, Šimunek J (2009). Effects of meteorological models on the solution of the surface energy balance and soil temperature variations in bare soils. J Hydrol (Amst), 373(3−4): 545–561

[34] Saito H, Šimunek J, Mohanty B P (2006). Numerical analysis of coupled water, vapor, and heat transport in the Vadose Zone. Vadose Zone J, 5(2): 784–800

[35] Schaefli B, van der Ent R J, Woods R, Savenije H H G (2012). An analytical model for soil-atmosphere feedback. Hydrol Earth Syst Sci, 16(7): 1863–1878

[36] van Bavel C H M, Hillel D I (1976). Calculating potential and actual evaporation from a bare soil surface by simulation of concurrent flow of water and heat. Agric Meteorol, 17(6): 453–476

[37] van Genuchten M Th (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J, 44(5): 892–898

[38] Wang D (2002). Dynamics of soil water and temperature in aboveground sand cultures used for screening plant salt tolerance. Soil Sci Soc Am J, 66(5): 1484–1491

[39] Webb E K (1970). Profile relationships: the log-linear range, and extension to strong stability. Q J R Meteorol Soc, 96(407): 67–90

[40] Zeller K F, Nikolov N T (2000). Quantifying simulations fluxes of ozone, carbon dioxide and water vapor above a subalpine forest ecosystem. Environ Pollut, 107(1): 1–20

[41] Zeng Y, Su Z, Wan L, Wen J (2011). A simulation analysis of the advective effect on evaporation using a two-phase heat and mass flow model. Water Resour Res, 47(10): W10529

[42] Zeng Y, Wan L, Su Z, Saito H, Huang K, Wang X (2009). Diurnal soil water dynamics in the shallow vadose zone (field site of China University of Geosciences, China). Environmental Geology, 58(1): 11–23

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg