Published online May 26, 2006

Numerical Analysis of Coupled Water, Vapor, and Heat Transport in the Vadose Zone
Hirotaka Saito,* Jiri Šimůnek, and Binayak P. Mohanty
for long-term landfill management, requires estimates ofwater and heat movement through engineered landfill
covers (e.g., Khire et al., 1997). Understanding surface
energy and water balances as well as liquid water, water
vapor, and heat movement in soils is critical for the performance evaluation of engineered surface covers for
waste containment in landfills in arid or semiarid regions(Scanlon et al., 2005). Traditional engineered covers
typically consist of multilayered resistive barriers withlow hydraulic conductivities to minimize water movement into the underlying waste. Recently, evapotranspiration (ET) covers have been developed (e.g., Hauser
and Gimon, 2004) to increase the water storage capacityof landfill covers and their protective function by removing water through evapotranspiration. These ET
cover designs thus need to be assessed mainly in terms ofthe near surface water and energy balance.
Another engineering application requiring an assessment of coupled liquid water, water vapor, and heattransport involves nuclear waste repositories. Radioactive decay in repositories may create steep temperature
gradients in adjacent soils and rocks, which can lead tounwanted consequences. Generated heat may cause
evaporation of soil water, and subsequent migration andcondensation of water vapor in cooler areas. This process may significantly change the physical characteristics
of surrounding materials due to the precipitation and dissolution of various minerals (e.g., Spycher et al., 2003).
Finally, in military applications, an accurate prediction ofwater contents and temperatures near the soil surface
may allow improved detection of buried land mines
since their presence can significantly change these twovariables (Šimůnek et al., 2001).
Early pioneering studies on interactions between liquid water, water vapor, and heat movement were reported by Philip and de Vries (1957), who provided amathematical description of liquid water and water
vapor fluxes in soils driven by both pressure head(isothermal) and soil temperature (thermal) gradients.
They derived the governing flow equation for nonisothermal flow as an extension of the Richards equation, which originally considered only the pressure headgradient. The theory of Philip and de Vries (1957) waslater extended by Nassar and Horton (1989), who additionally considered the effect of an osmotic potential
gradient on the simultaneous movement of water, solute, and heat in soils.
Heat transport and water flow are coupled by the
movement of water vapor, which can account for significant transfer of latent energy of vaporization. Soil
temperatures may be significantly underestimated whenthe movement of energy associated with vapor transport
is not considered. For example, Cahill and Parlange

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT
Vapor movement is often an important part in the total water flux inthe vadose zone of arid or semiarid regions because the soil moisture isrelatively low. The two major objectives of this study were to develop a
numerical model in the HYDRUS-1D code that (i) solves the coupledequations governing liquid water, water vapor, and heat transport,
together with the surface water and energy balance, and (ii) provides
flexibility in accommodating various types of meteorological information to solve the surface energy balance. The code considers the
movement of liquid water and water vapor in the subsurface to be
driven by both pressure head and temperature gradients. The heattransport module considers movement of soil heat by conduction,
convection of sensible heat by liquid water flow, transfer of latent heatby diffusion of water vapor, and transfer of sensible heat by diffusionof water vapor. The modifications allow a very flexible way of using
various types of meteorological information at the soil–atmosphereinterface for evaluating the surface water and energy balance. The
coupled model was evaluated using field soil temperature and watercontent data collected at a field site. We demonstrate the use of standard daily meteorological variables in generating diurnal changes inthese variables and their subsequent use for calculating continuouschanges in water contents and temperatures in the soil profile. Simulated temperatures and water contents were in good agreement withmeasured values. Analyses of the distributions of the liquid and vapor
fluxes vs. depth showed that soil water dynamics are strongly associated with the soil temperature regime.

T

HE simultaneous movement of liquid water, water

vapor, and heat in the vadose zone plays a criticalrole in the overall water and energy balance of the nearsurface environment of arid or semiarid regions in many
agricultural and engineering applications. Moisture near
the soil surface is influenced by evaporation, precipitation, liquid water flow, and water vapor flow, most ofwhich are strongly coupled. In arid regions, vapor movement is often an important part of the total water flux
since soil moisture contents near the soil surface usuallyare very low (Milly, 1984).
In agricultural applications, spatial and temporalchanges in surface soil moisture need to be well understood to achieve efficient and optimum water management. Vapor transport is also important since the actual
contact area between liquid water and seeds is oftenvery small such that seeds need to imbibe water fromvapor to germinate (Wuest et al., 1999).
In engineering applications, accurate assessment of
the volume of leachate from landfills, which is essential

H. Saito and J. Šimůnek, Dep. of Environmental Sciences, Univ. of
California, Riverside, CA 92521; B.P. Mohanty, Biological and Agricultural Engineering, Texas A&M Univ., College Station, TX 778432117. Received 12 Jan. 2006. *Corresponding author (hirotaka.saito@
ucr.edu).
Published in Vadose Zone Journal 5:784–800 (2006).
Original Researchdoi:10.2136/vzj2006.0007ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: DOY, day of the year; TDR, time domain reflectometry.

784

785

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org(1998) reported that 40 to 60% of the heat flux in the top
2 cm of a bare field soil of Yolo silt loam was due towater vapor flow. Fourier’s law describing heat transport
due to conduction (e.g., Campbell, 1985) thus needs tobe extended to include heat transport by liquid water
and water vapor flow. The general heat transport model
then considers movement of soil heat by (i) conduction,(ii) convection of sensible heat by liquid water flow, andtransfer of (iii) latent and (iv) sensible heat by diffusionof water vapor (e.g., Nassar and Horton, 1992).
The soil–atmosphere interface is an important boundary condition affecting subsurface movement of liquidwater, water vapor, and heat under field conditions.
Direct measurement of all components needed to fully
evaluate the water and energy balance at the soil surface, such as precipitation, evapotranspiration, heat flux,
and radiation, are rarely available at short time intervals.
In many situations, only standard daily meteorological
data from weather stations are available. If detailed
simulations of water and heat fluxes are required, thecomponents of the surface mass and energy balancesneed to be calculated from standard daily information,such as daily maximum and minimum air temperatures,
daily maximum and minimum relative humidities, theshortwave incoming solar radiation, and daily averagewind speed. One common approach to address this issue
is to first calculate approximate diurnal changes in themeteorological variables using relatively simple models(Ephrath et al., 1996) and to then use estimated values
for further investigation. Very few studies have demonstrated and evaluated such approaches for simulating
soil moisture and temperature in field soils.
Although it is widely recognized that the movementof liquid water, water vapor, and heat are closely coupled and strongly affected by each other, their mutualinteractions are rarely considered in practical applications. The effect of heat transport on water flow is often
neglected, arguably because of model complexity and alack of data to fully parameterize the model, but partly
also because of a lack of user-friendly simulation codes.
The two major objectives of this study thus were todevelop a numerical model that (i) solves the coupledequations governing liquid water, water vapor, and heattransport, together with the surface water and energybalance, and (ii) provides considerable flexibility inaccommodating various types of meteorological information that can be collected at selected time intervals(daily, hourly, or other time intervals). Complete evaluation of the movement of liquid water, water vapor, andheat in the subsurface can be accomplished by simultaneously solving the system of equations for the surface
water and energy balances, subsurface heat transport,
and variably saturated flow. These equations need tobe solved numerically since they are highly nonlinear
and strongly coupled. This coupled movement of liquidwater, water vapor, and heat in the subsurface, as wellas interactions of these subsurface processes with theenergy and water balances at the soil surface, was implemented in the HYDRUS-1D code (Šimůnek et al.,1998a). The HYDRUS-1D is a widely used, welldocumented, and tested public domain code for simu-

lating water and solute transport in soils (e.g., Scanlon
et al., 2002).
Model performance of the modified HYDRUS-1D
code was tested in this study by simulating continuouschanges in water contents, temperatures, and fluxes of abare field soil. Results were compared against field soil
temperature and water content data collected at different depths at a field site near the University of California Agricultural Experimental Station in Riverside, CA,
during the fall of 1995. We demonstrate how standard
daily meteorological variables can be used for simulating diurnal changes in soil water contents, temperatures,and liquid water, water vapor, and heat fluxes. In addition, the impact of vapor flow on predictions of soil
water contents and soil temperatures was investigated.
MATERIALS AND METHODS
ModelsLiquid Water and Water Vapor Flow
The governing equation for one-dimensional flow of liquidwater and water vapor in a variably saturated rigid porous
medium is given by the following mass conservation equation:

]u
]qL
]qv
52
2
2S
]t
]z
]z

[1]

where u is the total volumetric water content (m3 m23), qL andqv are the flux densities of liquid water and water vapor (ms21), respectively, t is time (s), z is the spatial coordinate positive upward (m), and S is a sink term usually accounting forroot water uptake (s21). The total volumetric water content isdefined as[2]u 5 ul 1 uv

where ul and uv are volumetric liquid water and water vapor
contents (expressed as an equivalent water content, m3 m23).
The flux density of liquid water, qL, is described using amodified Darcy law as given by Philip and de Vries (1957):

qL 5 qLh 1 qLT 5 2KLh

ð

Þ

]h
]T
1 1 2 KLT
]z
]z

[3]

where qLh and qLT are isothermal and thermal liquid water flux
densities (m s21), respectively, h is the pressure head (m), T is
the temperature (K), and KLh (m s21) and KLT (m2 K21 s21) are
the isothermal and thermal hydraulic conductivities for liquidphase fluxes due to gradients in h and T, respectively.
The flux density of water vapor, qv, can also be separated
into isothermal, qvh, and thermal, qvT, vapor flux densities(m s21) as follows (Philip and de Vries, 1957):

qv 5 qvh 1 qvT 5 2Kvh

]h
]T
2 KvT
]z
]z

[4]

where Kvh (m s21) and KvT (m2 K21 s21) are the isothermal andthermal vapor hydraulic conductivities, respectively. Combining Eq. [1], [3], and [4], we obtain the governing liquid water
and water vapor flow equation:



]u
]
]h
]T
]h
]T
5
KLh
1 KLh 1 KLT
1 Kvh
1 KvT
]t ]z
]z
]z
]z
]z


]
]h
]T
KTh
1 KLh 1 KTT
2S
[5]
2S5
]z
]z
]z

786

VADOSE ZONE J., VOL. 5, MAY 2006

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

where KTh (m s21) and KTT (m2 K21 s21) are the isothermal andthermal total hydraulic conductivities, respectively, and where:

KTh 5 KLh 1 Kvh[6]

KTT 5 KLT 1 KvT[7]Soil Hydraulic Properties
The pore-size distribution model of Mualem (1976) wasused to predict the isothermal unsaturated hydraulic conductivity function, KLh(h), from the saturated hydraulic conductivity and van Genuchten’s (1980) model of the soil water
retention curve, Eq. [9]:hi2mKLh (h) 5 Ks Sle 1 2 (1 2 S1/m
e )

[8]

where Ks is the saturated hydraulic conductivity (m s21), Se is
the effective saturation (unitless), and l and m are empiricalparameters (unitless). The parameter l was given a value of
0.5 as suggested by Mualem (1976). The parameter m was determined by fitting van Genuchten’s analytical model (vanGenuchten, 1980)

ul (h) 55

ur 1 hus 2 ur
1 1 |ah|n

im h , 0

1 dg
g0 dT

Þ

[10]

where GwT is the gain factor (unitless), which quantifies thetemperature dependence of the soil water retention curve(Nimmo and Miller, 1986), g is the surface tension of soil water(J m22), and g0 is the surface tension at 258C (571.89 g s22).
The temperature dependence of g is given byg 5 75:6 2 0:1425T 2 2:38 3 1024 T 2

u7/3
a
u2s

[11]

where g is in g s22 and T in 8C.
The isothermal, Kvh, and thermal, KvT, vapor hydraulicconductivities are described as (e.g., Nassar and Horton, 1989;
Noborio et al., 1996b; Fayer, 2000)[15]

and Da is the diffusivity of water vapor in air (m2 s21) at temperature T (K):Da 5 2:12 3 10

25

ð

T
273:15Þ

2[16]

The saturated vapor density, rsv (kg m23), as a function oftemperature is expressed as

rsv 5 1023

ð

exp 31:3716 2

Þ6014:79
2 7:92495 3 1023 T
T
T
[17]

while the relative humidity, Hr, can be calculated from thepressure head, h, using a thermodynamic relationship between
liquid water and water vapor in soil pores as (Philip and deVries, 1957):[9]

to water retention data using the RETC code (van Genuchten
et al., 1991). In Eq. [9], us and ur are the saturated and residual water contents (m3 m23), respectively, and a (m21), n(unitless), and m (51 2 1/n) are empirical shape parameters.
The thermal hydraulic conductivity function, KLT, in Eq. [5]is defined as follows (e.g., Noborio et al., 1996b):ð

H 5

Hr 5 exp

us h $ 0

KLT 5 KLh hGwT

where ua is the air-filled porosity (m3 m23), H is the tortuosity
factor as defined by Millington and Quirk (1961):

ð Þ
hMg
RT

[18]

When the liquid and vapor phases of water in soil pores arein equilibrium, the vapor density of the soil can be expressed as
the product of the saturated vapor density and the relative
humidity. The HYDRUS-1D code uses an enhancement factor, h, to describe the increase in the thermal vapor flux as aresult of liquid islands and increased temperature gradients inthe air phase (Philip and de Vries, 1957). An equation for theenhancement factor was derived by Cass et al. (1984) and isexpressed as

"
#4
u
2:6 u
h 5 9:5 1 3 2 8:5exp 2 1 1 pﬃﬃﬃﬃ
us
f c us

5 ð

Þ 6

[19]

where fc is the mass fraction of clay in the soil (unitless).
Heat Transport
The governing equation for the movement of energy in avariably saturated rigid porous medium is given by the following energy conservation equation:

]Sh
]qh
52
2Q
]t
]z

[20]

D
Mg
Kvh 5
r
Hr
rw sv RT

[12]

where Sh is the storage of heat in the soil (J m23), qh is the totalheat flux density (J m22 s21), and Q accounts for sources andsinks of energy (J m23 s21). The storage of heat is given byD
drhHr svdT
rw[13]

Sh 5 Cn Tun 1 Cw Tul 1 Cv Tuv 1 L0 uv 5 (Cn un 1 Cw ul
1 Cv uv )T 1 L0 uv 5 Cp T 1 L0 uv
[21]

KvT 5

where D is the vapor diffusivity in soil (m2 s21), rw is thedensity of liquid water (kg m23), rsv is the saturated vapor
density (kg m23), M is the molecular weight of water (Mmol21, 50.018015 kg mol21), g is the gravitational acceleration(m s22, 59.81 m s22), R is the universal gas constant (J mol21
K21, 58.314 J mol21K21), h is the enhancement factor(unitless, Cass et al., 1984), and Hr is the relative humidity(unitless). The vapor diffusivity, D, of the soil is defined as

D 5 H ua D a[14]

where T is a given temperature (K), un is the volumetric fraction of solid phase (m3 m23), Cn, Cw, Cv, and Cp are volumetricheat capacities (J m23 K21) of the solid phase, liquid water,water vapor, and moist soil (de Vries, 1963), respectively, and
L0 is the volumetric latent heat of vaporization of liquid water(J m23) given by

L0 5 Lw rw[22]

where Lw is the latent heat of vaporization of liquid water(J kg21; Lw [J kg21] 5 2.501 3 106 2 2369.2T [8C]).

787

www.vadosezonejournal.org

The total heat flux density, qh, 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, andlatent heat by vapor flow (de Vries, 1958):

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

]T
qh 5 2l(u)1 Cw TqL 1 Cv TqE 1 L0 qv
]z

[23]

where l(u) is the apparent thermal conductivity of soil
(J m21s21K21), and qL and qv are the flux densities of liquidwater and water vapor (m s21), respectively.
Combining Eq. [20], [21], and [23] results in the governing
equation for heat movement (e.g., Nassar and Horton, 1992;
Fayer, 2000):



]Cp T
]uv
]
]T
]qL T
]qv
1 L0
5
l(u)
2 Cw
2 L0
]t
]t
]z
]z
]z
]z
]qv T
[24]
2 Cw ST
2 Cv]zwhere the last term on the right side represents a sink of energy associated with root water uptake.
Soil Thermal Properties
The apparent thermal conductivity of soil, l(u), in Eq. [23]combines the thermal conductivity of the porous medium inthe absence of flow and the macrodispersivity, which is assumed to be a linear function of velocity (de Marsily, 1986;
Hopmans et al., 2002). The apparent thermal conductivity, l(u),may then be expressed as (e.g., Šimůnek and Suarez, 1993):l(u) 5 l0 (u) 1 bCw |qL |[25]

where b is the thermal dispersivity (m). Since the thermal
dispersivity plays an important role only when the liquid water
flux is very large, only a few studies have actually determined its
values (e.g., Hopmans et al., 2002). The thermal conductivity,l0(u), accounts for the tortuosity of the porous medium, andcan be described with a simple equation given by Chung andHorton (1987):

l0 (u) 5 b1 1 b2 u 1 b3 u0:5

[26]

where b1, b2, and b3 are empirical regression parameters (Wm21 K21). Chung and Horton (1987) provided average values
for the b coefficients for three textural classes (i.e., clay, loam,
and sand), which are implemented in the HYDRUS-1D code.
Soil thermal properties are much more influenced by watercontent than by textural differences since thermal propertiesof the different mineral components of the solid phase are
all approximately of the same order of magnitude (Jury andHorton, 2003).
Initial inspection of the fully coupled set of equations forliquid water, water vapor, and heat transport may suggest thatsignificantly more parameters are needed compared with more
approximate decoupled water flow and heat transport models.
This actually is not the case. The soil hydraulic properties are
fully described by six soil-specific hydraulic parameters thatappear in the van Genuchten–Mualem model (i.e., ur, us, a, n,
Ks, and l); these parameters are needed to simulate variably
saturated liquid water flow based on the Richards equation.
The soil thermal properties are described by means of three b
coefficients that appear in the thermal conductivity function.
Considering additionally thermal liquid flow, and both thermal
and isothermal vapor flow, does not increase the demand for
input parameters since KLT, Kvh, KvT, and Sh (Eq. [10], [12],[13], and [21]) can be fully described using literature values ofvarious physical properties, such as surface tension, the dif-fusivity of water vapor, and the heat capacities of the different
soil components.
Surface Energy Balance
Surface precipitation, irrigation, evaporation, and heat
fluxes are used as boundary conditions for liquid water andwater vapor flow and heat transport in field soils. Evaporationand heat fluxes can be calculated from the surface energy balance (e.g., van Bavel and Hillel, 1976; Noborio et al., 1996a;
Boulet et al., 1997):

Rn 2 H 2 LE 2 G 5 0

[27]

22

where Rn is net radiation (W m ), H is the sensible heat flux
density (W m22), LE is the latent heat flux density (W m22), L
is the latent heat (J kg21), E is the evaporation rate (kg m22s21),
and G is the surface heat flux density (W m22). While Rn and G
are positive downward, H and LE are positive upward.
Net radiation is defined as (e.g., Brutsaert, 1982)

Rn 5 Rns 1 Rnl 5 (1 2 a)St 1 (es Rldﬂ 2Rlu› )

[28]

22

where Rns is the net shortwave radiation (W m ), Rnl is the net
longwave radiation (W m22), a is the surface albedo (unitless),
St is the incoming (global) shortwave solar radiation (W m22),
es is the soil surface emissivity (unitless) representing thereflection of longwave radiation at the soil surface, Rldﬂ is theincoming (thermal) longwave radiation at the soil surface(downward flux) (W m22) as emitted by the atmosphere andcloud cover, and Rlu› is the sum of the outgoing (thermal)
longwave radiation emitted from the surface (vegetation andsoil) into the atmosphere (W m22).
The surface albedo can be described with a simple linear
expression relating the albedo with the surface water content
since it depends, especially for bare soils, on the soil surface
wetness. Van Bavel and Hillel (1976) proposed the following
simple formulae to calculate the surface albedo from thesurface wetness:a 5 0:25
a 5 0:10
a 5 0:35 2 utop

utop , 0:1
utop $ 0:25
0:1 # utop , 0:25

[29]

where utop is the water content at the surface.
Using the Stefan–Boltzmann law (e.g., Monteith andUnsworth, 1990), the net longwave radiation can be rewritten
as (e.g., Brutsaert, 1982):

Rnl 5 es Rldﬂ 2 Rlu› 5 es ea sTa4 2 es sTs4[30]

where the subscripts a and s are used for variables of the atmosphere and the soil, respectively. The emissivity of the atmosphere depends both on the air temperature and humidity,
while the emissivity of the soil surface depends on the watercontent and the vegetation at the soil surface. The emissivity ofbare soil can be expressed as a function of the volumetric watercontent (van Bavel and Hillel, 1976) as

es 5 min(0:90 1 0:18utop ; 1:0)

[31]

leading to an emissivity of 0.9 for a dry surface and |0.98 when
the soil is saturated (u | 0.45). Following Idso (1981), the atmospheric emissivity can be calculated as follows:ea 5 0:70 1 5:95 3 1025 ea exp(1500/Ta )[32]

where ea is the atmospheric vapor pressure (kPa), which can be
expressed as a function of Ta asea 5 0:611exp
17:27ðTa 2 273:15Þ
Hr
Ta 2 35:85

[33]

788

VADOSE ZONE J., VOL. 5, MAY 2006

For a sky partially covered with clouds, the longwave radiation is calculated with (Monteith and Unsworth, 1990)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

hi
Rldﬂ 5 (1 2 0:84c)ea 1 0:84c sTa4

[34]

where c (unitless) is the fraction of cloud cover (or the cloudiness factor) of the sky. This equation was derived from the average temperature difference between the atmosphere and theclouds in England (Monteith and Unsworth, 1990). The cloudiness factor is readily estimated from the atmospheric transmission coefficient for solar radiation, Tt, as (Campbell, 1985)

0 # c 5 2:33 2 3:33Tt # 1[35]

A value of the incoming shortwave solar radiation St (Wm22) at any given time and location can be calculated by taking
into account the position of the sun (e.g., van Bavel and Hillel,1976) using the following equation (Campbell, 1985):

St (t) 5 max(Gsc Tt sin e, 0)

[36]

22

where Gsc is the solar constant (1360 W m ), and Tt (unitless)is defined as the ratio of the measured daily global solar radiation Stm (W m22) and the daily potential global (extraterrestrial) radiation Ra (Wm22):

Tt 5

Stm
Ra

[37]

The last term of Eq. [36], e (rad), is the solar elevation angle
given by (Monteith and Unsworth, 1990)

sin e 5 sin B sin d 1 cos B cos d cos

2p(t 2 t0 )24[38]

Ts 2 TarH[39]

where Ta is the air temperature (K), Ts is the soil surface
temperature (K), Ca is the volumetric heat capacity of air(J m23 K21), and rH is the aerodynamic resistance to heattransfer (s m21).
Since surface evaporation is controlled by atmospheric conditions, surface moisture, and moisture transport in the soil,a model is needed that accounts for all of these factors.
Evaporation, E, is often calculated asr 2 rva
E 5 vsrv 1 rs

rs 5 2805 1 4140(us 2 utop )[40]

where rvs is the water vapor density at the soil surface (kg m23),
rva is the atmospheric vapor density (kg m23), rv is the aerodynamic resistance to water vapor flow (s m21), and rs is thesoil surface resistance to water vapor flow that acts as anadditional resistance along with the aerodynamic resistance(s m21) (Camillo and Gurney, 1986). Evaporation models commonly use only the aerodynamic resistance in the denominatorof Eq. [40] (e.g., Milly [1984] and some existing codes, such as
SHAW [Flerchinger et al., 1996] and UNSATH [Fayer, 2000]).
This approach may not be correct since it overestimates theevaporation rate for dry soils because of the assumption of Eq.
[18] that equilibrium exists in the soil (Kondo et al., 1990).
When the soil surface is dry, water vapor in larger pores isdynamically transported toward the atmosphere and its
density in pores is not in equilibrium with the average watercontent of a given depth. The soil surface resistance is thus

[41]

Calculation of the sensible heat flux, H, or evaporation rate,
E, using Eq. [39] or [40] requires knowledge of the aerodynamic resistances to heat transfer, rH, and to water vapor
flow, rv, respectively. These may be calculated from the surface
wind speed as follows (Campbell, 1985):


1
zref 2 d 1 zH
lnrH 5 rv 5
1 cH
zH
U*k

ð

Þ

[42]

where U* is the friction velocity (m s21), k is the von Karman
constant, which has a value of 0.41, zref is the reference heightof temperature measurements (m), zH is the surface roughness
for the heat flux (m), d is the zero-plane displacement (m), andcH is the atmospheric stability correction factor for the heat
flux (unitless). The assumption that the resistance to heat
flux is equal to that to vapor flow has long been used (e.g., vanBavel and Hillel, 1976; Flerchinger et al., 1996). The friction
velocity in Eq. [42] is defined as (Campbell, 1985)"

ð

U* 5 uk ln

where t is the local time within a day and t0 is the time of solar noon.
Among a number of available accepted models (e.g., Wang
and Bras, 1998), the sensible heat flux, H, can be simply defined as (e.g., van Bavel and Hillel, 1976)

H 5 Caused to account for an extra resistance to water vapor flow insoil pores. This resistance depends strongly on both soil structure and texture. The relationship between the soil surface
resistance, rs, and the surface water content has been empirically formulated in many studies and typically has an
exponential form (e.g., van de Griend and Owe, 1994), indicating that the surface resistance increases dramatically as thesoil dries out. The soil surface resistance in this study is determined using the following formula (Camillo and Gurney, 1986):

Þ

zref 2 d 1 zm
1 cm
zm

21
[43]

where u is the mean wind speed (m s21) at height zref, zm is thesurface roughness for the momentum flux (m), and cm is theatmospheric stability correction factor for the momentum flux
(unitless). For bare soils, d is equal to zero, while typicalsurface roughness values of 0.001 m are used for both zH andzm (Oke, 1978). In general, calculations of rH and rv require an
iterative procedure since stability correction factors must be
evaluated from variables that themselves depend on the stability correction factors. Camillo and Gurney (1986) proposed
an approximate approach to calculate the stability correction
factor, which is used also in this study so no need exists for theiterative procedure.
Meteorological Variables Generation
Solving the energy balance equation (Eq. [27]) at a timeinterval of interest requires values of meteorological variablesat the same or similar time intervals. Weather stations, however, do not always provide standard meteorological data attime intervals of interest. Although several models exist forcalculating diurnal changes from daily average values (e.g.,
Ephrath et al., 1996), finding the best model was not an objective of this study. We used relatively simple approaches togenerate continuous values of the meteorological variablesfrom available daily information.
When the daily maximum and minimum air temperatures
are the only information provided by a weather station, continuous estimates of the air temperature, Ta, can be obtainedusing a trigonometric function with a period of 24 h as follows(Kirkham and Powers, 1972):
t 2 13
Ta 5 T 1 AT cos 2p
24

ð

Þ

[44]

789

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

where T is an average daily temperature (8C), AT is the amplitude of the cosine wave (8C) that can be calculated from thedifference between the daily maximum and minimum temperatures, and t is a local time within the day. The argument of
the cosine function shows that the highest temperature isassumed to occur at 1300 h and the lowest at 0100 h.
Since daily temperature variations typically show a cyclic
behavior throughout the day, it is reasonable to assume thatthe relative humidity also shows such a cyclic pattern. Similarly
as for air temperature (Eq. [44]), a trigonometric function canbe used to calculate continuous values of the relative humidityfrom daily information (e.g., Gregory et al., 1994). When both
the daily maximum and minimum relative humidity values areavailable, continuous values of the relative humidity can begenerated as well using a cosine function with a period of 24 h
as follows:
t 2 tmax
Hr 5 H 1 AH cos 2p
24

ð

Þ

[45]

where H is the average daily relative humidity (unitless), AH is
the amplitude of the cosine wave that can be calculated fromthe difference between the maximum and minimum relative
humidity values, t is a local time within the day, and tmax is thehour of the day when the relative humidity is at its maximum.
The relative humidity is generally at its minimum during daytime when the air temperature and the wind speed are at theirmaxima. Gregory et al. (1994) showed that the relative humidity was highest at around 0500 and 0600 h in Lubbock, TX.
In this study, we used 0500 h for the peak time of the relative
humidity. The diurnal change in the saturation vapor density
may be calculated from generated temperature values using
Eq. [17], while the corresponding atmosphere vapor density,rva, can be calculated by multiplying the generated relative
humidity (Eq. [45]) and the saturation vapor density.
It is well known that wind transports heat and effectively
mixes the soil–atmospheric boundary layer. Wind is generally
highly variable in space and time since it involves mostly turbulent flow and as such is characterized by random fluctuations in its speed and direction (Campbell, 1977). Because
many meteorological models require continuous inputs of thewind speed, continuous values of the wind speed must somehow be calculated from daily information. In this study, weused the following simple approach based on the maximum-tominimum wind speed ratio Ur, which is defined as

Umax
Ur 5
Umin

[46]

where Umin and Umax (m s21) are the unknown minimum andmaximum wind speeds of the day, respectively. This ratio may
be determined from prior knowledge or calibrated on available data. The maximum and minimum wind speeds can thenbe calculated from the daily average wind speed as follows:

Umax 52UrU
1 1 Ur[47]

Umin 52U
1 1 Ur[48]

where U (m s21) is a daily average wind speed. Cyclic behaviorof the wind speed U during the day is then modeled using thedaily average wind speed U as follows (Gregory, 1989):


t 2 tmax
U 5 U 1 (Umax 2 U)cos 2p
24

ð

Þ

[49]

where t is the time of the day, and tmax is the time when the
maximum wind speed occurs. In Gregory’s (1989) study, themaximum wind speed occurred at 1530 h.
Calculation of continuous values of the incoming shortwave
radiation using Eq. [36] requires continuous values of the transmission coefficient, Tt (Eq. [37]). When only the daily incoming
shortwave solar radiation, Stm, is available, a daily average transmission coefficient must be used since no reliable generation isavailable. Calculation of the daily average transmission coefficient using Eq. [37] requires values of the daily potential global
solar radiation, which depends on the latitude of the location ofinterest and the day of year as follows (Campbell, 1985):

Ra 524(60)
Gsc dr (vs sin B sin d 1 cos B cos d sin vs )p[50]

where dr is the relative distance from the earth to the sun, vs is
the sunset hour angle (rad), d is the angular distance of the sun(i.e., solar declination, rad), and B is the latitude of the location(rad).

Field Data
The coupled numerical model and other meteorologicalmodules implemented in HYDRUS-1D were tested using asoil temperature and moisture data set collected at the fieldsite near the University of California Agricultural Experimental Station in Riverside, CA. Temporal soil temperature andwater content variations were measured using horizontally inserted thermocouples and time domain reflectometry (TDR)
probes, respectively, near the soil–air interface at depths of 2,
7, and 12 cm during the fall of 1995 (24 November (DOY 328)–
5 December (DOY 339); see Mohanty et al., 1998). The TDR
probes consisted of three 20-cm-long rods spaced 2.5 cm apart(total width 5 cm). The experimental field was irrigated twice
during the measurement period on DOY 334 and 335 with 0.55
and 0.20 cm of water, respectively, using a sprinkler system
containing several laterals and outlet ports. Irrigation rates
were measured using catch cans near the measuring location.
Soil temperatures and water contents at each depth were
measured at 20- and 40-min intervals, respectively.
Field Measurements
Figure 1 shows measured water contents and temperaturesat three depths at the measurement location. Soil water contents were not recorded from noon of DOY 330 to noon of
DOY 331 due experimental problems. The soil temperature
data show a typical sinusoidal diurnal behavior at all threedepths, while the soil water content data do not reveal a distinct diurnal pattern, only noisy fluctuations. When TDR isused to measure the soil water content, noisy fluctuations areoften inherent in the interpretation of TDR waveforms (e.g.,
Cahill and Parlange, 1998). Irrigations on DOY 334 and335 increased the water contents at depths of 2 and 7 cm, while
the water content remained almost constant at the depth of
12 cm during the experiments. As a result of irrigation, averagewater contents before and after irrigation were significantly
different (P , 0.01), even at the depth of 7 cm. More details on
temperature and water content variations at the site can befound in Mohanty et al. (1998).
Soil Hydraulic Data
The soil at the experimental site was characterized as an
Arlington fine sandy loam (coarse-loamy, mixed, thermic
Haplic Durixeralf). Soil water retention curves were measured

790

VADOSE ZONE J., VOL. 5, MAY 200630

Temperature [ oC]

20

15

10

5

+

2cm
7cm
12cm

++++
++++
+
+
+
+++
++ +
+
+
+
+
+
+
+
+
+
++
+++
++++
+
+++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+++
+
+
+
+
+
+
+
+
+
+
+
+++
++
+
++ + ++ + + +
+
+
+
+++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
++++
+
+
+
+
+
+ ++ + +
+
+
+ +++
+
++ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ + +
+++
+
+
+
+
+
++++
+
+++
+
+
+
++
+
+
+
+
+
+ +++ + + +++ ++ +++++++++++++ +++ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
++ + ++ + ++ + ++ +
++++
+
+
+++
+ ++ + +++++ ++ +
+
+++
+++
++++
+

328

329

330

331

332

333

334

335

336

337

338

339

DOY
0.3
0.25

Water Content [m3m-3]

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

25

0.2
0.15

+
+ ++++++ ++++ ++++++++ + ++ + + ++ ++ ++++ ++ +++ ++++ +++++++++ +++ +++++++++++++++++++++++ + +++
++ ++
+
+
++++++++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
++++ ++++++ ++++
+
+
+
+
+
+
+
+ +++++++++++++++++ ++
+
+
+
+ ++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
++
+++ ++++ ++++++ + ++ ++++ + +++ + + ++ + +++++ +++++++++ + + ++++++
+ +
+++

+
+++
++
+++++++++++++++++++++++++++++++++++++++++++++++
++
+
+
+
+
+
++ ++ ++++ +

0.1
0.05
0

328

329

330

331

332

333

334

335

336

337

338

339

DOY
Fig. 1. Soil temperatures and water contents measured at three depths (2, 7, and 12 cm) from 24 November (Day of the Year [DOY] 328) to 5
December (DOY 339) at the experimental site in Riverside, CA.

in the laboratory using the Tempe Cell (Soil Moisture Equipment Corp.) method for pressure heads down to 28 m and apressure chamber (Soil Moisture Equipment Corp.) for lower
pressure heads (down to 2150 m). Figure 2 shows a plot of thewater retention curve for a soil core sample collected 4 m eastof the measurement location. The figure includes both measured drainage (drying) and imbibition (wetting) data. The van
Genuchten (1980) analytical model (Eq. [9]) was fitted to theretention data using the RETC code (van Genuchten et al.,1991), leading to ur 5 0.011 m3 m23, us 5 0.445 m3 m23, a 5
0.0277 cm21, and n 5 1.38.
Several studies previously measured the saturated hydraulicconductivity, Ks, of the Arlington fine sandy loam soil using
different methodologies (e.g., Šimůnek et al., 1998b; Wanget al., 1998a, 1998b). Table 1 lists soil hydraulic conductivities
obtained by the different researchers. Measured values varydepending on the methodology used. In this study, we usedthe value of 34.2 cm d21 obtained by Wang et al. (1998a) using
a tension infiltration experiment for two reasons. First, theirmeasurements were conducted at the same site as used inour study. Also, their Ks value was close to the median value
obtained with the three different methods listed in Table 1.
The clay fraction of the Arlington fine sandy loam at Riverside needed for the enhancement factor (Eq. [19]) was 8.8%
(Šimůnek et al., 1998b).

Meteorological Data
Standard daily meteorological data in California can beobtained from approximately 400 weather stations that are apart of the California Irrigation Management Information
System (CIMIS) Network operated by the California Department of Water Resources (DWR). Data collected daily by
DWR were taken from the University of California Integrated
Management Program website (http://www.ipm.ucdavis.edu/
WEATHER/about_weather.html, accessed November 2005,
verified 26 Mar. 2006). Data were obtained for the automatic
CIMIS weather station located at UC Riverside, Agricultural
Operations Department Area, CA (338589N, 1178209W, elevation 310.9 m). Daily meteorological variables used in this study
are summarized in Table 2. Data were also available from the
CIMIS website (http://wwwcimis.water.ca.gov/cimis/welcome.
jsp, accessed November 2005, verified 26 Mar. 2006), fromwhich we downloaded daily and hourly values of standard
meteorological variables. We used hourly values from this website to validate the various generations of continuous weather
parameters (see also Table 2 for variables used in this study).

Numerical Implementation
The total variably saturated water flow and heat transport
equations (Eq. [5] and Eq. [24]) were both solved numerically

791

www.vadosezonejournal.org

Table 2. Weather variables recorded at the University of California-Riverside Agricultural Operations Department Area.

Water Content [m3m-3]

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

0.5

Variable

0.4

Reporting interval

Condition
1.5 m above the surface

Precipitation, inches

Daily average, dailymax. & min., hourly
Daily total, hourly

Relative humidity, %Vapor pressure, mbars

Daily max/min, hourly
HourlySolar radiation,Langleys

Daily total, hourly

Air temperature, °F

0.3

0.2

0.1

0
0
10

drainage
imbibition
VG ModelWind speed, miles h& direction
1

10

2

10

3

10

Pressure Head [cm]Fig. 2. Water retention curve measured on a soil core sample collected
at the experimental site in Riverside, CA. Both drainage and imbibition curves are plotted (open symbols: drainage, closed symbols:
imbibition). The van Genuchten (1980) model (VG) is fit to themeasured data.

using the finite element method for spatial discretizationand finite differences for the temporal discretization. The
HYDRUS-1D software package (Šimůnek et al., 1998a) was
modified to consider the simultaneous movement of liquidwater, water vapor, and heat, as well as the surface energy andwater balances. The modified version of HYDRUS-1D either
can consider continuous values of net radiation and other meteorological variables measured at any time interval, which arethen used directly in the energy balance equation (Eq. [27]), or
can calculate continuous values from daily meteorological data
as described above.
At each time step, the net radiation (Rn), the sensible heat
flux (H), and the latent heat flux (LE) were calculated to obtain the surface heat flux (G) by solving the surface energybalance equation (Eq. [27]), which is subsequently used as aknown heat flux boundary condition. Calculation of Rn, H
and LE requires knowledge of the temperature and pressure
head of the soil surface. These were obtained by sequentiallysolving the governing equations of the surface energy balance, variably saturated water flow, and heat transport twice,while updating temperatures and water contents between thetwo runs.
Numerical solutions of transport equations often exhibit
undesired numerical oscillations. An appropriate combinationof space and time discretization can often prevent such oscillations. In the original HYDRUS-1D, undesired oscillations
Table 1. Soil hydraulic conductivities of Arlington fine sandy loam
estimated by different researchers using different methodologies. Listed values are all geometric means of measured values.
Saturated hydraulicconductivity
21

Šimůnek et al.
(1998b)Wang et al.
(1998a)Wang et al.
(1998b)Tension disk infiltrometer with
Wooding’s method
Tension disk infiltrometer with
Wooding’s method
Tension disk infiltrometer with
the Darcy–Buckingham method
Tension disk infiltrometer withthe sorptivity method
Tension disk infiltrometer with
Wooding’s method
Guelph permeameter

Daily average, hourly

4

10

Methodology

21

In a 20 cm diameter
gauge
1.5 m above the surface
Calculated frommeasured relative
humidity and airtemperature data
By pyranometer at
2.0 m above thesurface
2.0 m above the surface

cm d32.0, 20.4

30.6 6 3.0
34.2
53.2 6 13.2
19.9 6 7.81
22.0 6 17.0

were minimized or eliminated by using two dimensionless numbers that characterize the space and time discretization, i.e.,
the grid Peclet and Courant numbers (Šimůnek et al., 1998a).
A dimensionless Courant number is used also in the modified
HYDRUS-1D code to stabilize the numerical solution of theheat transport equation. Neglecting terms associated with thetransfer of latent heat by water vapor in Eq. [24], the Courant
number, Cu, may be defined as

Cu 5(Cw qL 1 Cv qv )Dt
Cp Dz

[51]

where Dt and Dz are temporal and spatial steps, respectively.
The maximum permitted time step, Dtmax, is then calculated
assuming a value of 1 for the Courant number.
The numerical solution can also be stabilized by specifying
the maximum nodal temperature difference, DTmax, allowed
during a given time step (e.g., 0.258C). The maximum allowed
time step is calculated as

Dtmax 5

DTmax
Dt
DTi

[52]

where DTi is the maximum nodal temperature change at theprevious time step Dt. The maximum permitted time step forthe heat transport calculations is then the smaller value of Eq.
[51] or [52].
Initial and Boundary Conditions
The soil profile was considered to be 50 cm deep, with observation nodes located at depths of 2, 7, and 12 cm for comparing calculated temperatures and volumetric water contents
with measured values. A constant nodal spacing of 2 mm wasused, leading to 251 discretization nodes across the problem
domain. Calculations were performed for a period of 12 d from24 Nov. (DOY 328) to 5 Dec. (DOY 339) 1995. No naturalprecipitation was recorded at the weather station during thistime period.
Initial water contents (a constant value of 0.13 m3 m23) andsoil temperatures (variable with depth) were determined frommeasured values on DOY 328. Boundary conditions at the soil
surface for liquid water, water vapor, and heat transport were
determined from the surface water and energy balance equations. When the soil was irrigated on DOY 334 and 335 (0.55
and 0.20 cm of water, respectively), irrigation rates were usedas the flux boundary condition. For heat transport, a heat flux
at the soil surface, G, was calculated from the surface energybalance equation (Eq. [27]). The temperature of irrigationwater was assumed to be equal to the air temperature. Zeropressure and temperature gradients were used as the bottom

792

VADOSE ZONE J., VOL. 5, MAY 2006

boundary conditions at the depth of 50 cm. These conditions
assume that the water table is located far below the domain ofinterest and that heat transfer across the lower boundary occurs only by convection of liquid water and water vapor.

Generated continuous changes in the air temperature, relative humidity, and wind speed are comparedwith hourly measured values obtained from the nearby
weather station in Fig. 3 for the entire simulation period.
A conventional value of 3 was used for the maximum-tominimum wind speed ratio, Ur (Eq. [46]). Most calcu-

Air temperature [oC]35
Measured
Approximated

30
25
20
15
10
5
0
327

329

331

333

335

337

339

335

373

339

DOY
Relative humidity [%]

80
70
60
50
40
30
20
10
0
327

329

331

333

DOY
8

Wind speed [m/s]Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

RESULTS AND DISCUSSION
Meteorological Model Verification

6

4

lated values compared well with hourly measured values, even for the wind speed, which fluctuated morethan the other variables due to inherent randomness.
Although the general trends of the relative humidity
were well described, some maximum values were slightly
overestimated. Using the generated continuous air temperatures and relative humidities, diurnal variations in
the atmosphere vapor pressure were calculated using
Eq. [17] and compared with hourly vapor pressuresobtained from the weather station (Fig. 4). As shown in
Fig. 4, diurnal changes in the vapor pressure calculated
from the generated meteorological variables closely fitted the hourly vapor pressure values at the weather station. Results presented in Fig. 3 and 4 demonstrate thevalidity of using the generated continuous diurnal variations in meteorological variables.
Meteorological variables, including net radiation, airtemperature, and relative air humidity, were measured
at the experimental site every 5 min from 21 Oct. (DOY294) to 1 Nov. (DOY 305) and every 10 min from 2 Nov.
(DOY 306) to 25 Nov. (DOY 329) 1995 before the soil
temperature and water content were measured. To investigate whether or not generated diurnal changes in
the meteorological variables can be used to reproduce
continuous changes in the net radiation, calculations
were performed from DOY 294 to 329 with the same
domain conditions as described above. Daily meteorological data during this period were again obtained
from the weather station in Riverside (Table 2). Diurnalchanges in the meteorological variables needed to calculate the net radiation at any given time were calculated
from the daily data.
Figure 5 shows temporal changes in both the measured and simulated net radiations for a 10-d period.
Positive values indicate an incoming flux to the land surface (downward), while negative radiation values imply
an outgoing flux from the land surface (upward). Calculated and measured daytime net radiations matchedreasonably well, whereas nighttime agreement was lessgood, especially around 0000 h of DOY 310, 311, 315,
and 319. This may have been caused by clouds or fog onthese days, which reduced the net radiation toward zero.
The current model takes cloudiness into account only
during the daytime (Eq. [35]); modeling such effects
during nights is not easy. Figure 5 also depicts simulated
net shortwave and longwave radiations, which were used
to calculate net radiation values. Simulated net radiation
values before and after the displayed period in Fig. 5matched measured values similarly well. Results presented in Fig. 5 demonstrate that our choice of meteorological models to calculate the net radiation is reasonable.

2

Simulated Soil Temperatures and Water Contents
0
327

329

331

333

335

337

339

DOY
Fig. 3. Diurnal changes of three meteorological variables—air temperature, relative humidity, and wind speed—generated from daily
information during the simulation period along with hourly values
measured at the weather station near the study site from 23 November (Day of the Year [DOY] 327) to 6 December (DOY 340).

Soil water contents and soil temperatures at the experimental site were numerically simulated from DOY 328
to 339 using the coupled liquid water, water vapor, andheat transport module implemented in the HYDRUS1D code. Figure 6 shows measured and simulated watercontents at depths of 2, 7, and 12 cm. Predicted watercontents follow fairly well the measured values at all

793

www.vadosezonejournal.org

6
5
0.6

4

0.4

3
2

0.2

0
327

Vapor pressure
Vapor density

Vapor density [gm-3]

Vapor pressure [kPa]

7

0.8

1
0

329

331

333

335

337

339

DOY
Fig. 4. Vapor densities calculated from generated air temperatures (Eq. [44]) and relative humidities (Eq. [45]) along with vapor pressurescalculated hourly at the weather station from 23 November (Day of the Year [DOY] 327) to 6 December (DOY 340).

three depths during the entire simulation period, including the gradual decrease in the water content at the 2-cm
depth before irrigation. The observed noise in measured
values is inherent to the TDR measurements (e.g., Cahill
and Parlange, 1998). Rapid increases in the water content at the 2-cm depth after 0.55 and 0.20 cm of irrigationwater was applied on DOY 334 and 335, respectively,
were predicted reasonably well, with a slight overestimation after the first irrigation and some underestimation after the second irrigation. The small increase in u atthe 7-cm depth after the two irrigations was also well
predicted. The simulated and measured water contents
at the 12-cm depth were almost constant during the simulation period.
Figure 7 depicts simulated and measured soil temperatures at three depths. The simulated and measured
temperatures both show typical sinusoidal diurnal behavior. The amplitude of both the simulated and measured daily temperature variations decreased with depth
due to attenuation of the transported heat energy. Be-fore the irrigation on DOY 334, the simulated temperature amplitude was larger than the observed one. Juryand Horton (2003) presented an analytical solution for
the temperature amplitude at a particular depth z:rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
rﬃﬃﬃﬃﬃﬃﬃﬃﬃ
vv Cp
A 5 A0 exp z
5 A0 exp z
[53]
2KT
2 l

ð

Þ

ð

Þ

where A0 and A are temperature amplitudes at the soil
surface and the specified depth (8C), respectively, v is theangular frequency (s21), and KT is the thermal diffusivity(m2 s21), being the ratio of the thermal conductivity andthe heat capacity. The temperature amplitude can be reduced by either increasing the soil heat capacity (Eq. [21])or decreasing the soil thermal conductivity, both of which
would lead to a lower thermal diffusivity. For example,
reducing the thermal diffusivity by 50% would lead to amuch better fit of the measured temperatures at all threedepths. After the irrigation, simulated values almost perfectly reproduced the measured ones. Apart from the50
Measured
Net Radiation
Short Rad
Long.Rad.

40

Radiation [MJ/m2/day]

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1

30

20

10

0

-10
310

311

312

313

314

315

316

317

318

319

320

DOY
Fig. 5. Net radiation calculated from 6 Nov. 1995 (Day of the Year [DOY] 310) to 16 Nov. 1995 (DOY 320) using the standard daily meteorological
data obtained from the weather station at the UC-Riverside Agricultural Operations Department Area in Riverside, CA. Measured hourly net
radiation values were obtained at the experimental site near the University of California Agricultural Experimental Station in Riverside. Net
radiation was calculated as the sum of the net incoming shortwave solar radiation and net outgoing longwave radiation.

VADOSE ZONE J., VOL. 5, MAY 2006

depth = 2cm

Temperature [oC]

30

0.2

0.1
Observed
Simulated

0
327

329

331

25
20
15
10
5

333

335

337

0
327

339

depth = 2cm

Observed
Simulated

329

331

30

Temperature [oC]

depth = 7cm

3-3

Volumetric Water Content [m m ]

335

337

339

335

337

339

335

337

339

35

0.3

0.2

0.1

depth = 7cm

25
20
15
10
5

0
327

0
327
329

331

333

335

337

331

333

DOY
35

0.3

3

329

339

DOY

30

depth = 12cm

Temperature [oC]

-3

333

DOY

DOY

Volumetric Water Content [m m ]Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

35

0.3

3

-3

Volumetric Water Content [m m ]

794

0.2

0.1

depth = 12cm

25
20
15
10
5

0
327

0
327
329

331

333

335

337

339

DOY
Fig. 6. Simulated and measured volumetric water contents during thesimulation time period (Day of the Year [DOY] 328–DOY 339) atthree depths (2, 7, and 12 cm) at the experimental site.

temperature amplitude, simulated and measured temperatures generally agreed at all three depths, with most simulated values being within a few degrees of the measured
values. A large number of formulas can be used to calculate various meteorological parameters, such as air
emissivity or boundary layer resistances. For example,
Mahfouf and Noilhan (1991) showed that, dependingon the choice of an evaporation model, simulated soil
surface temperatures can be easily altered by 38C or more.
Differences between simulated and measured values
in our study decreased with depth and were ,28C at the
12-cm depth.
Calculated surface energy fluxes during the simulation period are depicted in Fig. 8. Although the net radiation values were fairly constant during the simulation

329

331

333

DOY
Fig. 7. Simulated and measured soil temperatures during the simulation time period (Day of the Year [DOY] 328–DOY 339) at threedepths (2, 7, and 12 cm) at the experimental site.

period, variations in daily surface heat and sensible heat
fluxes were more pronounced. This confirms that the dynamics of the surface energy components and the partitioning of the energy are influenced by many different
land and atmospheric attributes and as such cannot be
predicted solely from the net radiation. While surface
heat fluxes are needed as boundary conditions for subsurface heat transport, their variations are directly related to simulated soil temperature variations (Fig. 7).
Latent heat fluxes are in general very small during dry
periods because of a lack of soil moisture for evaporation (i.e., before irrigation); however, after the soil was
irrigated on DOY 334 and 335, the latent heat fluxes
increased significantly during the next few days as soil
moisture near the surface increased. The peak of thelatent heat flux was predicted to be approximately at795

www.vadosezonejournal.org

Net Radiation
Latent Heat Flux
Sensible Heat Flux
Surface Heat Flux

Heat Flux [MJ/m2/day]

20

10

0

-10

-20
327

329

331

333

335

337

339

DOY
Fig. 8. Simulated surface energy balance components at the experimental site: net radiation, sensible heat flux, latent heat flux, and surface heat flux
from 24 November (Day of the Year [DOY] 328) to 5 December (DOY 339). All fluxes are positive upward, except for the net radiation where
positive value indicates incoming radiation to the soil surface (downward) and negative value outgoing radiation from the soil surface (upward).

the pressure head and relative humidity during the night
are discussed below.
Figure 10 shows temporal variations in the calculatedtransmission coefficient, the surface albedo, the surface
emissivity, and the air emissivity during the simulation
period. All of these variables are necessary to calculate
the net radiation. The transmission coefficient and air
emissivity show large temporal variations. Variationsin the air emissivity were due to changes in the vapor
pressure (Eq. [32]), temporal variations of which arealso distinct (Fig. 4). The soil surface albedo and soil
emissivity, on the other hand, were fairly constant during
the simulation period, except for the albedo immediately after irrigation on DOY 334 and 335. Wet surfaces
do not reflect radiation as much as dry surfaces, resulting
in higher net radiation values during and after the irri-

1200 h every day as expected. The sensible heat flux decreases after irrigation since a large part of the incomingenergy is then used for evaporation (Fig. 8). The direction of the sensible heat flux varies depending on thetime of day; the soil surface is heated during daytime,leading to upward sensible heat fluxes (positive values),
while the sensible heat flux in general is downward during
night (negative values) when the soil surface cools down.
The temporal variation in the latent heat flux may beinvestigated also by analyzing temporal variations in thepressure head and relative humidity at the soil surface(Fig. 9). As soil water evaporates into the atmosphere
during the daytime, the pressure head of the soil surface
is reduced to as low as 210 000 cm (pF 5 4). The pressure head (or relative humidity) then increases sharply
during the nighttime. The mechanisms of increases in
0

80

-4000
60
-6000
40
-8000
Pressure head

-1000020

Relative Humidity [%]

100

-2000

Pressure Head [m]Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

30

Relative humidity-12000
327

0
329

331

333

335

337

339

DOY
Fig. 9. Temporal changes of the pressure head and relative humidity at the soil surface during the simulation period from 24 November (Day of theYear [DOY] 328) to 5 December (DOY 339). The corresponding surface relative humidity was calculated from the surface pressure head using
Philip and de Vries’ equilibrium model (Eq. [18]).

796

VADOSE ZONE J., VOL. 5, MAY 2006
0.3
0.25

Surface albedo [-]

0.6
0.5
0.4
0.3
0.2

0.2
0.15
0.1
0.05

0.1
329

331

333
DOY

335

337

0
328

339

1

1

0.95

0.95

Air emissivity [-]Transmission coefficient [-]

0.7

0
327

Surface emissivity [-]Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

0.8

0.9

0.85

0.8

0.75
327

330

332

334
DOY

336

338

0.9

0.85

0.8

329

331

333

335

337

339

DOY

0.75
327

329

331

333

335

337

339

DOY

Fig. 10. Temporal changes of transmission coefficient, surface albedo, surface emissivity, and atmosphere emissivity required for calculations of the
net radiation at a given time of day for the simulation time period (Day of the Year [DOY] 328–DOY 339). These variables are calculated directlyby the program.

gations (Fig. 7). Figure 10 demonstrates that temporalvariations in most variables are not simple and need to
be well estimated before they should be used in calculations of the net radiation.

Liquid Water and Water Vapor Fluxes
Calculated distributions of the liquid water and water
vapor fluxes vs. depth for both isothermal and thermal
components before and after irrigation are depicted inFig. 11 and 12, respectively. While the plots in Fig. 11
show fluxes on a typical dry day (from 1200 h of DOY
329 to 1200 h of DOY 330), the impact of irrigation on
DOY 335 on these fluxes is shown in Fig. 12. The thermal liquid flux was almost always negligible comparedwith the other three fluxes on the dry day.
Distributions of the water content and soil temperature vs. depth before and after irrigation are shown in
Fig. 13. While a large downward thermal water vapor
flux occurred at noon of DOY 329 (Fig. 11) due to alarge downward temperature gradient (Fig. 13), the isothermal liquid water and water vapor fluxes are upward(Fig. 11) because of an upward pressure head gradient.
As the downward thermal vapor flux component
decreased in the top 1 cm of the profile, the net water
vapor flux in this layer was upward because of a large
upward isothermal water vapor flux. This upward water
vapor flux served as a source of water for evaporationfrom the soil surface to the atmosphere during the day.
The source for this upward water vapor flux was clearlyliquid water transported by the isothermal liquid water
flux from the depths of |6 cm and above (Fig. 11).
Liquid water arriving at |1-cm depth changed to water
vapor, which then moved toward the soil surface. Thisimplies that actual vaporization (phase change) of liquidwater did not occur at the soil surface, but at a depth of|1 cm. The depth where water vapor flow starts dominating the overall water flux will be referred to as thedrying front in the remainder of our discussion.
Most fluxes were small and upward at 0000 h of thedry day (Fig. 11, DOY 330) when the temperature gradient changed from downward to upward (Fig. 13, top
right). Analysis of isothermal fluxes reveals again thatliquid water is moving from the deeper layers (|10-cm
depth or above) toward the drying front (|1-cm depth),
where it is converted to water vapor. Since the latent
heat flux during the night is small, not much water vapor
is transported to the atmosphere. Instead, liquid water
moving toward the soil surface accumulates near the
drying front, resulting in a small increase in the watercontent (Fig. 13, top left). This process, which can also be
observed in the surface pressure head or relative humidity (Fig. 9), substantially increases during the night.
Soil moisture is eventually transported back to deeper
layers by downward thermal water vapor flow.
Around 1200 h of DOY 330, the four fluxes were very
similar to those of the previous day (1200 h of DOY329), except that the drying depth was slightly lower at|1.5 cm below the soil surface (Fig. 11). Since there wasno precipitation during this time period, the drying front
moved downward, leading to a thicker dry surface layer.
The overall water dynamics in the soil during a typicaldry day can be summarized as follows. A continuous
upward isothermal liquid water flux occurs below thedrying front, accompanied by a large downward thermal
water vapor flux during the daytime and an upward but
small thermal water vapor flux during the night. This
leads to a daily circulation pattern of soil moisture below

797

www.vadosezonejournal.org

0

DOY=329.50-5Depth [cm]Depth [cm]-10-15Isothermal liquid
Thermal liquid
Isothermal vapor
Thermal vapor

-20

-0.1

0

0.1

-10
-15
-20
-25
-0.2

0.2

-0.1

Flux [cm/day]

0

0.1

0.2

Flux [cm/day]
DOY=330.50Depth [cm]

-5
-10
-15
-20
-25
-0.2

-0.1

0

0.1

0.2

Flux [cm/day]Fig. 11. Calculated vertical distributions of the thermal and isothermal fluxes of liquid water and water vapor at 1200 and 0000 h of a typical dry 24-h
period during Days of the year [DOY] 329 and 330. Positive values indicate upward and negative values downward fluxes.

the drying front. The processes are quite different above
the drying front where the upward isothermal vapor fluxis the most dominant flux component. Liquid water
moving upward by the isothermal liquid water flux (andwater vapor by the thermal water vapor flux during thenight) evaporates at the drying front and is further

0

DOY=335.0-5Depth [cm]Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

-5

-25
-0.2

DOY=330.0-10-15Isothermal liquid
Thermal liquid
Isothermal vapor
Thermal vapor

-20
-25
-0.2

-0.1

0

0.1

0.2

Flux [cm/day]Fig. 12. Vertical profiles of thermal and isothermal fluxes of liquid water
and water vapor at 0000 h of Day of the Year (DOY) 335. Lightirrigation was applied at 1200 h of DOY 334. Positive values indicate
upward fluxes, while negative values indicate downward fluxes.

transported as water vapor, mainly by the pressure headgradient, toward the surface. During the daytime, water
vapor escapes into the atmosphere from the soil surface(i.e., evaporation), while during the night, some moisture accumulates near the surface.
Flux profiles are significantly different after irrigation
when the isothermal liquid water flux dominates(Fig. 12). As a result of irrigation, the water content
near the surface increases, leading to a large downward
isothermal liquid water flux. Close to the soil surface, alarge upward isothermal liquid water flux exists, which
supplies water for evaporation. A small upward thermal
liquid water flux, which was almost always negligible forthe dry day, was observed near the soil surface. There
was also an upward thermal vapor flux due to the temperature gradient as the surface cooled down due tohigh latent heat fluxes after irrigation (Fig. 8). The upward isothermal and thermal liquid water and thermal
water vapor fluxes all contribute to high latent heat
fluxes, even during the night.
Minor sinusoidal diurnal behavior observed in thesimulated water content at 2-cm depth before the irrigation (Fig. 6, top) is attributed to upward and downward
water vapor flow resulting from temperature gradientvariations. Even though it is difficult to see daily sinusoidal behavior in our water content data because of

0

0

-5

-10

Depth [cm]Depth [cm]VADOSE ZONE J., VOL. 5, MAY 2006

-10

-15
DOY 329.5
DOY 329.75
DOY 330
DOY 330.25
DOY 330.5

-20

-25

0

0.05

-20

-30

-40

0.1

0.15
3

-50

0.2

0

-5

-10

Depth [cm]

0

-10

-15

-25
0

DOY 334.5
DOY 334.75
DOY 335
DOY 335.25
DOY 335.4

0.05

20

30

Temperature [ C]

0

-20

10

o

-3

Water Content [m m ]Depth [cm]Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

798

-20

-30

-40

0.1

0.15

0.2

Water Content [m3m-3]

-50

0

10

20

30

o

Temperature [ C]

Fig. 13. Simulated profiles of water contents (left) and soil temperatures (right) at 1200 h of Day of the Year (DOY) 329 (top row) and DOY 334(bottom row), and 6, 12, 18, 24 h after. Light irrigation was applied at 1200 h of DOY 334.

the strong background noise (Fig. 6), other researchers
also observed similar diurnal variations (e.g., Cahill andParlange, 1998). There is an ongoing debate whether
the sinusoidal variations are real or simply a reflection oftemperature effects on the measurement apparatus,
such as TDR (Or and Wraith, 2000). Although this debate is still unresolved, it is clear that the observed sinusoidal variations of the water content can be produced
only when vapor flow is considered. Since vapor fluxes
decrease significantly after irrigation (liquid water fluxes
become dominant), the diurnal variations of water content are much smaller in wet soils.

SUMMARY AND CONCLUSIONS
We modified the HYDRUS-1D software package toincrease its flexibility in accommodating various typesof meteorological information and to solve the coupledequations describing liquid water, water vapor, and heattransport in soils. The liquid water, water vapor, andheat transport equations were solved simultaneously sothat interactions between different processes involved
can be taken into account. Boundary conditions at
the soil surface were obtained by solving the water andenergy balance equations. When only daily meteorological information is available, continuous values of energybalance components can be calculated using generated
values of meteorological variables.
The model was evaluated using soil temperature andwater content data collected at three depths (2, 7, and12 cm) at an experimental site near the University ofCalifornia Agricultural Experimental Station in Riverside, CA, during the fall of 1995. Using measured soil
hydraulic properties and daily meteorological data, reasonable agreement between simulated and measured
soil temperatures and water contents was obtained. Numerical analysis showed the dynamics of soil moisturein field soils, with both water flow and heat transport
being strongly affected by vapor transport. Consideration of vapor transport may result in small daily sinusoidal variations in water contents at shallow depths.
Vapor flow, often neglected in previous numericalsimulations of water flow and heat transport, was fully
considered in this study. Temporal dynamics of soil waterand energy fluxes at the soil surface cannot be describedwithout considering water vapor transport. Since water
vapor transport is mainly driven by temperature gradients, heat transport must be taken into account as well.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

This implies that it is important to consider coupledliquid water, water vapor, and heat transport in predictions of soil water and heat dynamics in field soils.
This study shows that limited daily meteorological information can still be used to obtain reasonable predictions of soil temperatures and water contents when the
coupled liquid water, water vapor, and heat transportmodel is used. Detailed measurements of meteorological data, such as net radiation variations, may not benecessary even when knowledge of continuous (diurnal)
soil temperatures and water contents is required.
ACKNOWLEDGMENTS
This work was supported in part by SAHRA (Sustainabilityof Semi-Arid Hydrology and Riparian Areas) under the STC
Program of the National Science Foundation, Agreement no.
EAR-9876800 and the Terrestrial Sciences Program of the
Army Research Office (Terrestrial Processes and Landscape
Dynamics and Terrestrial System Modeling and Model Integration). B.P. Mohanty acknowledges the support by NASA–
Global Water and Energy Cycle grant no. NAG5-11702. We
would like to thank Rien van Genuchten for his constructive
comments on the manuscript.

REFERENCES
Boulet, G., I. Braud, and M. Vauclin. 1997. Study of the mechanisms of
evaporation under arid conditions using a detailed model of thesoil–atmosphere continuum. Application to the EFEDA I experiment. J. Hydrol. 193:114–141.
Brutsaert, W. 1982. Evaporation into the atmosphere: Theory, history,
and applications. D. Reidel Publ., Dordrecht, the Netherlands.
Cahill, A.T., and M.B. Parlange. 1998. On water vapor transport infield soils. Water Resour. Res. 34:731–739.
Camillo, P.J., and R.J. Gurney. 1986. A resistance parameter for baresoil evaporation models. Soil Sci. 141:95–105.
Campbell, G.S. 1977. An introduction to environmental biophysics.
Springer-Verlag, New York.
Campbell, G.S. 1985. Soil physics with BASIC: Transport models for
soil–plant systems. Elsevier, New York.
Cass, A., G.S. Campbell, and T.L. Jones. 1984. Enhancement of thermal water vapor diffusion in soil. Soil Sci. Soc. Am. J. 48:25–32.
Chung, S.-O., and R. Horton. 1987. Soil heat and water flow with apartial surface mulch. Water Resour. Res. 23:2175–2186.
de Marsily, G. 1986. Quantitative hydrogeology, groundwater hydrology for engineers. Academic Press, Orlando, FL.
de Vries, D.A. 1958. Simultaneous transfer of heat and moisture inporous media. Trans. Am. Geophys. Union 39:909–916.
de Vries, D.A. 1963. The thermal properties of soils. p. 210–235. In
R.W. van Wijk (ed.) Physics of plant environment. North Holland,
Amsterdam.
Ephrath, J.E., J. Goudriaan, and A. Marani. 1996. Modelling diurnal
patterns of air temperature, radiation, wind speed and relative
humidity by equations from daily characteristics. Agric. Syst. 51:
377–393.
Fayer, M.J. 2000. UNSAT-H version 3.0: Unsaturated soil water andheat flow model. Theory, user manual, and examples. Pac. Northw.
Natl. Lab., Richland, WA.
Flerchinger, G.N., C.L. Hanson, and J.R. Wight. 1996. Modeling evapotranspiration and surface energy budgets across a watershed. Water
Resour. Res. 32:2539–2548.
Gregory, J.M. 1989. Wind data generation for Great Plains locations.
Pap. 892664. In Proc. Am. Soc. Agric. Eng. Int. Winter Meet., New
Orleans. 12–15 Dec. ASAE, St. Joseph, MI.
Gregory, J.M., R.E. Peterson, J.A. Lee, and G.R. Wilson. 1994. Modeling wind and relative humidity effects on air quality. p. 183–190. In
AWMA/AGU Int. Specialty Conf. on Aerosols and Atmosph.
Optics, Snowbird, UT. 26–30 Sept. 1994. Midwest Plan Service,
Ames, IA.

799

Hauser, V.L., and D.M. Gimon. 2004. Evaluating evapotranspiration(ET) landfill cover performance using hydrologic models [Online].
Available at http://www.hqafcee.brooks.af.mil/products/techtrans/
LandfillCovers/ (accessed 7 June 2005; verified 26 Mar. 2006). Air
Force Center for Environmental Excellence, Brooks City-Base, TX.
Hopmans, J.W., J. Simunek, and K.L. Bristow. 2002. Indirect estimation of soil thermal properties and water flux using heat pulse probe
measurement: Geometry and dispersion effects. Water Resour.
Res. 38(1) doi:10.1029/2000WR000071.
Idso, S.B. 1981. A set of equations for full spectrum and 8- to 14-mm
and 10.5- to 12.5-mm thermal radiation from cloudless skies. Water
Resour. Res. 17:295–304.
Jury, W.A., and R. Horton. 2003. Soil physics. 6th ed. John Wiley &Sons, New York.
Khire, M.V., C.H. Benson, and P.J. Bosscher. 1997. Water balance
modeling of earthen final covers. J. Geotech. Geoenviron. Eng.
123:744–754.
Kirkham, D., and W.L. Powers. 1972. Advanced soil physics. John
Wiley & Sons, New York.
Kondo, J., N. Saigusa, and T. Sato. 1990. A parameterization of evaporation from bare soil surface. J. Appl. Meteorol. 29:385–389.
Mahfouf, J.F., and J. Noilhan. 1991. Comparative study of various formulations of evaporation from bare soil using in situ data. J. Appl.
Meteorol. 30:1354–1365.
Millington, R.J., and J.P. Quirk. 1961. Permeability of porous media.
Trans. Faraday Soc. 57:1200–1207.
Milly, P.C.D. 1984. A simulation analysis of thermal effects on evaporation from soil. Water Resour. Res. 20:1087–1098.
Mohanty, B.P., P.J. Shouse, and M.Th. van Genuchten. 1998. Spatiotemporal dynamics of water and heat in field soil. Soil Tillage Res.
47:133–143.
Monteith, J.L., and M.H. Unsworth. 1990. Principles of environmental
physics. 2nd ed. Edward Arnold, London.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:
513–521.
Nassar, I.N., and R. Horton. 1989. Water transport in unsaturated
nonisothermal salty soil: II. Theoretical development. Soil Sci. Soc.
Am. J. 53:1330–1337.
Nassar, I.N., and R. Horton. 1992. Simultaneous transfer of heat,water, and solute in porous media: I. Theoretical development. Soil
Sci. Soc. Am. J. 56:1350–1356.
Nimmo, J.R., and E.E. Miller. 1986. The temperature dependence ofisothermal moisture vs. potential characteristics of soils. Soil Sci.
Soc. Am. J. 50:1105–1113.
Noborio, K., K.J. McInnes, and J.L. Heilman. 1996a. Two-dimensionalmodel for water, heat, and solute transport in furrow-irrigated soil.
I. Theory. Soil Sci. Soc. Am. J. 60:1001–1009.
Noborio, K., K.J. McInnes, and J.L. Heilman. 1996b. Two-dimensionalmodel for water, heat, and solute transport in furrow-irrigated soil:
II. Field evaluation. Soil Sci. Soc. Am. J. 60:1010–1021.
Oke, T.R. 1978. Boundary layer climates. Methuen & Co., London.
Or, D., and J.M. Wraith. 2000. Comment on “On water vapor transport
in field soils” by Anthony T. Cahill and Marc B. Parlange. Water
Resour. Res. 36:3103–3106.
Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous
materials under temperature gradient. Trans. Am. Geophys. Union
38:222–232.
Scanlon, B.R., M. Christman, R.C. Reedy, I. Porro, J. Šimůnek, andG.N. Flerchinger. 2002. Intercode comparisons for simulating waterbalance of surficial sediments in semiarid regions. Water Resour.
Res. 38(12):1323 doi:10.1029/2001WR001233.
Scanlon, B.R., R.C. Reedy, K.E. Keese, and S.F. Dwyer. 2005. Evaluation of evapotranspiration covers for waste containment in arid
and semiarid regions in the southwestern USA. Vadose Zone J.
4:55–71.
Šimůnek, J., and D.L. Suarez. 1993. The UNSATCHEM-2D code forsimulating two dimensional variably saturated water flow, heattransport, carbon dioxide transport, and solute transport with major
ion equilibrium and kinetic chemistry. Version 1.1. Res. Rep. no.
128. U.S. Salinity Lab., Riverside, CA.
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1998a. The
HYDRUS-1D software package for simulating the one dimensional
movement of water, heat, and multiple solutes in variably-saturated

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

800

VADOSE ZONE J., VOL. 5, MAY 2006

media. Version 2.0. IGWMC-TPS-70. Int. Ground Water ModelingCenter, Colorado School of Mines, Golden.
Šimůnek, J., D. Wang, P.J. Shouse, and M.Th. van Genuchten. 1998b.
Analysis of a field tension disc infiltrometer experiment by parameter estimation. Int. Agrophys. 12:167–180.
Šimůnek, J., J.M.H. Hendrickx, and B. Borchers. 2001. Modelingtransient temperature distributions around landmines in homogeneous bare soils. p. 387–397. In A.C. Dubey et al. (ed.) Proc. Detection and Remediation Technologies for Mines and Minelike
Targets, 5th. Vol. 4394. SPIE, Bellingham, WA.
Spycher, N.F., E.L. Sonnenthal, and J.A. Apps. 2003. Fluid flow andreactive transport around potential nuclear waste emplacement
tunnels at Yucca Mountain, Nevada. J. Contam. Hydrol. 62–63:
653–673.
van Bavel, C.H.M., and D.I. Hillel. 1976. Calculating potential andactual evaporation from a bare soil surface by simulation of concurrent flow of water and heat. Agric. For. Meteorol. 17:453–476.
van de Griend, A.A., and M. Owe. 1994. Bare soil surface resistance toevaporation by vapor diffusion under semiarid conditions. Water
Resour. Res. 30:181–188.
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44:892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC codefor quantifying the hydraulic functions of unsaturated soils. Version
1.0. EPA/600/2–91/065. U.S. Salinity Lab., Riverside, CA.
Wang, D., S.R. Yates, and F.F. Ernst. 1998a. Determining soil hydraulicproperties using tension infiltrometers, time domain reflectometry,
and tensiometers. Soil Sci. Soc. Am. J. 63:318–325.
Wang, D., S.R. Yates, B. Lowery, and M.Th. van Genuchten. 1998b.
Estimating soil hydraulic properties using tension infiltrometerswith varying disk diameters. Soil Sci. 163:356–361.
Wang, J., and R.L. Bras. 1998. A new method for estimation of sensible
heat flux from air temperature. Water Resour. Res. 34:2281–2288.
Wuest, S.B., S.L. Albrecht, and K.W. Skirvin. 1999. Vapor transport vs.
seed–soil contact in wheat germination. Agron. J. 91:783–787.