Masaru Sakai*Dep. of Environmental Sciences
Univ. of California
Riverside, CA 92521

Nobuo Toride
Graduate School of Bioresources
Mie Univ.
Tsu City
Japan 514-8507

Jiří Šimůnek

Dep. of Environmental Sciences
Univ. of California
Riverside, CA 92521

The diffusion of warm, humid air into an initially cold, dry, sandy column was analyzed to
study the movement of water vapor and liquid water under nonisothermal and low watercontent conditions. The analysis was performed using the HYDRUS-1D code. While thewater retention curve of sand was measured experimentally, the unsaturated hydraulicconductivity function was inversely estimated from the observed water content proﬁles in
the column. The estimated unsaturated hydraulic conductivity function displayed a shapethat reﬂected distinct processes of capillary pore water ﬂow and ﬁlm ﬂow at high and low
water contents, respectively. Four components of the total water ﬂux, including thermal
and isothermal liquid water and water vapor ﬂuxes, were evaluated using the calibrated soil
hydraulic properties. Evaporation and condensation rates were calculated based on watermass balance. Water vapor entered the soil column at the hot surface and condensed at thecold bottom. Subsequently, liquid water moved upward and evaporated at the moisture front
in the middle of the column where the relative humidity decreased below unity. Liquid water
and water vapor then circulated between the bottom and the moisture front, accompaniedby condensation and evaporation processes. The impact of the enhancement factor in thethermal vapor diffusion term could not be clearly identiﬁed from available experimental
water content proﬁles. Increases in liquid water ﬂow and the evaporation rate could becompensated for by increases in vapor ﬂow and the condensation rate. Additional data
would be needed to fully evaluate the effect of the enhancement factor.

C

omplex processes involved in the surface energy balance
produce diurnal ﬂuctuations of the air temperature at thesoil–atmosphere interface, inducing either water vapor evaporation or condensation at different times of the day. Since
water vapor ﬂow in semiarid and arid regions can represent a
major part of the overall water ﬂow, it is important to take itinto account together with liquid water ﬂow when evaluatinghydrologic ﬂuxes and solute transport (Milly, 1984; Scanlon
et al., 2003). Additionally, water vapor ﬂow signiﬁcantly affects the movement of heat, since it transports a substantial
amount of energy as the latent heat of vaporization (Cahill andParlange, 1998; Saito et al., 2006). Hence, the simultaneous
evaluation of liquid water, water vapor, and heat movement insoils is essential for many applications in research as well as inwater management.
Penman (1940) ﬁrst developed a theory of vapor ﬂow
in soils based on Fick’s law. Later, Philip and de Vries (1957)
mathematically extended Penman’s theory to describe coupled
movement of liquid water and water vapor under nonisothermal conditions. They divided the total water ﬂux into four
Soil Sci. Soc. Am. J. 73:707-717
doi:10.2136/sssaj2008.0094
Received 14 Mar. 2008.
*Corresponding author (sakai@ucr.edu).
© Soil Science Society of America
677 S. Segoe Rd. Madison WI 53711 USA
All rights reserved. No part of this periodical may be reproduced ortransmitted in any form or by any means, electronic or mechanical,
including photocopying, recording, or any information storage andretrieval system, without permission in writing from the publisher.
Permission for printing and for reprinting the material contained
herein has been obtained by the publisher.
SSSAJ: Volume 73: Number 3 • May–June 2009

SOIL PHYSICS

Water and Vapor Movement with Condensationand Evaporation in a Sandy Column

components, consisting of liquid water and water vapor ﬂuxes
driven by either water content or temperature gradients. Philipand de Vries (1957) also introduced an enhancement factor forvapor ﬂux due to the temperature gradient. Water vapor could
move through a “liquid island” in a soil pore, condensing at
one end and evaporating at the other. Another reason for theenhanced vapor movement is that the actual local temperature
gradient in air-ﬁlled pores may be signiﬁcantly greater than
the average temperature gradient for the bulk soil. Since theenhancement factor is difﬁcult to measure directly, Cass et al.
(1984) estimated this parameter indirectly based on the pressure head dependence of the thermal conductivity. Nassar andHorton (1989) and Nassar et al. (1992b) further extended thePhilip and de Vries (1957) model by including osmotic effects
on liquid water and water vapor movement.
It is essential to precisely determine soil hydraulic properties at low water saturations when describing liquid water ﬂow
in dry soils. Although many methods exist to determine soil
hydraulic properties (Klute and Dirksen, 1986), unsaturated
hydraulic conductivity at low water saturations is very difﬁcult
to measure. Only methods based on ﬂow under accelerated
gravitational (centrifugal) conditions seem to be practical for
low water saturations (e.g., Nimmo, 1990). In recent years,
several inverse methods, such as the multistep outﬂow method(Eching and Hopmans, 1993; Inoue et al., 1998) and the evaporation method (Šimůnek et al., 1998; Fujimaki and Inoue,2003) have been developed to estimate unsaturated hydraulicproperties from experimentally obtained data involving transient water ﬂow. Although these methods usually cover only
a relatively high water content range (i.e., above the measurement range of a tensiometer), it should be possible to extend
the application of these methods to drier soils, providing that a707

soil hydraulic property model correctly describes retention andconductivity properties in the dry water content range.
One of the most widely used functions for describing unsaturated hydraulic properties is van Genuchten’s (1980) set
of closed-form equations for the soil water retention curve,
coupled with Mualem’s (1976) pore-size distribution modelfor the unsaturated hydraulic conductivity function. The van
Genuchten (1980) model sometimes fails to describe soil water
retention properties, however, and underestimates the unsaturated hydraulic conductivity at low water contents. Tuller andOr (2001) showed a relation between the unsaturated hydraulic conductivity and the pressure head based not only on the
pore-size distribution, but also ﬁlm ﬂow on the surface of soil
particles, a process that dominated at low saturations. Fayer
and Simmons (1995) modiﬁed the van Genuchten soil water
retention model to better describe soil water retention at low
water contents and described the unsaturated hydraulic conductivity function by coupling their model with the Mualem(1976) model. The Mualem model contains the pore-connectivity coefﬁcient that accounts for the effects of pore connectivity and tortuosity, which can have a signiﬁcant impact on
the gradient of the unsaturated hydraulic conductivity withrespect to the pressure head. Although Mualem (1976), basedon an analysis of 45 soil samples, suggested the average value ofthe pore-connectivity coefﬁcient to be 0.5, some recent studies had utilized other values (e.g., Romano and Santini, 1999;
Hopmans et al., 2002; Schaap et al., 2001).
Although numerical simulations of the simultaneous
movement of liquid water, water vapor, and heat in soils withmedium or high water contents have been performed in thepast (e.g., Nassar et al., 1992b), few simulations have been
performed for relatively dry soils. This is largely because of
uncertainty about unsaturated hydraulic conductivity. The objective of this study was to evaluate the vapor condensation
experiments performed by Miyazaki (1976) using the Philipand de Vries (1957) theory and the Fayer and Simmons (1995)
hydraulic property model implemented in the HYDRUS1D code (Šimůnek et al., 2008a,b; Saito et al., 2006). Theselaboratory experiments involved coupled movement of liquidwater, water vapor, and heat in sandy columns under an imposed temperature gradient. During the experiments, water
vapor diffused from the hot and humid end of the columns
and subsequently condensed at the cold dry end. In our study,
the water retention curve parameters were ﬁrst ﬁtted to independently measured water retention curve data using the Fayer
and Simmons (1995) model. The numerical model was thencalibrated against experimental data of volumetric water content proﬁles by inversely estimating the pore-connectivity coefﬁcient to determine the unsaturated hydraulic conductivity
function for low water contents while using independently ﬁtted water retention curve parameters. Further analyses of various processes were then performed using the calibrated model.
Second, four components of the total water ﬂux, including
liquid and vapor ﬂuxes due to pressure head and temperature
gradients, were evaluated using calculated pressure head andtemperature proﬁles. Third, condensation and evaporation
rates were calculated from the mass balance and liquid water
and water vapor ﬂuxes. Fourth, the impact of the enhancement
factor on the coupled transport involving condensation and708

evaporation processes was investigated. Finally, the uncertaintyof the inverse parameter estimation procedure using the available data set was determined and the need for additional data
to guarantee a unique solution was recognized.

MATERIALS AND METHODS
Condensation Experiments
Miyazaki (1976) conducted one-dimensional laboratory column
experiments to observe the diffusion of water vapor from hot and humid air into cold and dry soil and to evaluate the various mechanisms
involved in water vapor ﬂow. The experiments were performed using
Hamaoka dune sand, which had the typical particle size distributionof dune sand. It contained 3% of soil particles in the range of 0.002
to 0.02 mm, 7% of soil particles in the range of 0.02 to 0.1 mm, and
90% of soil particles >0.1 mm. The saturated hydraulic conductivityof 34.6 m d−1 was measured using the falling-head method. Hamaoka
dune sand with an initial volumetric water content of 0.0045 m3 m−3
was packed uniformly at a dry bulk density of 1.6 g cm−3 in six acrylic
columns with 10-cm height and 10-cm diameter. The sand columns
were then placed vertically in a chamber with a constant air temperature of 37°C and relative humidity of 85 to 90%. The tops of thecolumns were exposed to the moist air, while the closed bottoms ofthe columns were maintained at 20°C by circulating constant-temperature water along the base of the columns using a water pump.
The walls of the columns were insulated using foam polystyrene toachieve one-dimensional heat ﬂow. Since differences between the observed temperatures at the center and close to the wall of the column
were <0.6°C (or 12°C m−1), the experiment was assumed to have
one-dimensional heat ﬂow.
Because of the difference in the vapor density between the twoends of the columns, water vapor diffused from the hot and humid
end of the column toward the cold and dry end, where it condensed.
Subsequently, the condensed liquid water moved upward due to thepressure head gradients. Temperature proﬁles in the columns weremonitored with thermocouples located at depths of 0, 2, 5, 8, and10 cm at the center of the columns. The amount of cumulative water
vapor diffusion into the column was measured based on the column
weights. Identical experiments were repeated on all six columns to obtain water content proﬁles at different times (2, 5, 10, 16, 21, and 30d). Each sand column at those times was sectioned into seven 1- and2-cm-thick layers to gravimetrically determine the water contents.

Numerical Model
Unsaturated Water Flow
Darcy’s law for liquid water ﬂux, qL, under time-variable temperature conditions has to consider both the ﬂow driven by the temperature gradient, accounting for the temperature dependence of thesurface tension, as well as the term driven by the pressure head gradient (Philip and de Vries, 1957). Since the vapor density is also a function of the pressure head and the temperature, Fick’s law for the vapor
ﬂux, qv, also consists of two (thermal and isothermal) components.
Thus, the total water ﬂux, qTotal, in a variably saturated soil is the sumof the liquid water ﬂux and the water vapor ﬂux:

qTotal = qL + qv = qLh + qLT + qvh + qvT∂T∂h∂T [1]⎛ ∂h ⎞= − K Lh ⎜ + 1⎟ − K LT
− K vh− K vT∂z∂z∂z⎝ ∂z ⎠SSSAJ: Volume 73: Number 3 • May–June 2009

where qLh and qLT are the isothermal and thermal liquid water ﬂuxes[L T−1], respectively; qvh and qvT are the isothermal and thermal water
vapor ﬂuxes [L T−1], respectively; h is the pressure head [L]; T is thetemperature [K]; z is the spatial coordinate positive upward [L]; KLh[L T−1] and KLT [L2 K−1 T−1] are the isothermal and thermal hydraulic conductivities for liquid water ﬂuxes, respectively; and Kvh [L T−1]
and KvT [L2 K−1 T−1] are the isothermal and thermal vapor hydraulicconductivities, respectively. Inserting Eq. [1] into the mass conservation equation leads to the Richards equation for nonisothermal conditions (e.g., Nassar and Horton, 1997; Noborio et al., 1996):∂q∂θ= − Total
∂t
∂z
∂ ⎡
∂T
∂h
∂T ⎤
⎛ ∂h ⎞= ⎢ K Lh ⎜ + 1⎟ + K LT
+ K vh+ K vT
∂z ⎣
∂z
∂z
∂z ⎥⎦
⎝ ∂z ⎠[2]

where θ is the total volumetric water content [L3 L−3] and t is time [T].
This mass conservation equation can be divided into two equations that account for liquid water and water vapor contents (Nassar
and Horton, 1992):∂θL∂q=− L −E
∂t
∂z

[3]

∂θv
∂q
=− v +E∂t∂z[4]where θL is the volumetric liquid water content [L3 L−3], θv is the volumetric vapor content expressed as an equivalent amount of liquid water
per unit volume of the soil [L3 L−3], and E is the evaporation or condensation rate [T−1] (positive for evaporation and negative for condensation).

Soil Hydraulic Properties
A model proposed by Fayer and Simmons (1995) (denoted below
as the Fayer model) was used for the soil water retention curve. They
modiﬁed the residual water content in the van Genuchten model (vanGenuchten, 1980) to better represent the soil water retention curve at
low water contents. The model accounts for the adsorption of water
on soil under dry conditions:nθL = χθa + (θs − χθa ) ⎡1 + (−αh ) ⎤
⎣
⎦

−m

[5]

where θs is the saturated water content [L3 L−3], and α [L−1], n (dimensionless), m ( = 1 − 1/n), and θa [L3 L−3] are empirical shapeparameters. The ﬁrst term, χθa, denotes adsorption of water on soiland χ is described asχ (h ) = 1 −ln (−h )ln (−hm )[6]where hm is the pressure head at the water content equal to 0, and isgenerally taken to be −107 cm (Rossi and Nimmo, 1994).
The closed form of the isothermal hydraulic conductivity function
for the liquid water ﬂux, KLh, is obtained for the Fayer model by substituting Eq. [5] into the pore-size distribution model of Mualem (1976):
2

2⎡ dS h ⎤e∫l ⎡ Γ s (h )⎤0⎢⎥= K s Se ⎢K Lh = K s Se⎥
⎢ 1 dS h ⎥⎣ Γ max ⎦⎢⎣ ∫0 e⎥⎦Sel

SSSAJ: Volume 73: Number 3 • May–June 2009[7]

where Ks is the saturated hydraulic conductivity [L T−1], Se (= θL/θs)is the effective liquid saturation (dimensionless), and l is the poreconnectivity coefﬁcient that was assumed to be 0.5 in the original
study. The deﬁnitions of Γs(h) and Γmax are presented in Fayer andSimmons (1995).
The thermal hydraulic conductivity for the liquid water ﬂux,
KLT, is deﬁned as (e.g., Noborio et al., 1996)⎛
1 dγ ⎞K LT = K Lh ⎜ hGwT⎟
γ 0 dT ⎠⎝[8]

where γ is the surface tension of soil water [M T−2] and γ0 is thesurface tension at 25°C (= 71.89 g s−2). The gain factor, GwT, whichcorrects the temperature dependence of the surface tension, was ﬁxed
at 7 in this study (Nimmo and Miller, 1986).
Hydraulic conductivities for vapor phase ﬂuxes due to pressure
head, Kvh, and temperature, KvT, gradients are deﬁned as

K vh =

D
Mg
ρvs
Hr
ρw
RTabs

[9]

K vT =dρD
ηH r vs
ρw
dT

[10]

respectively, where D is the vapor diffusivity in soil [L2 T−1], ρw
is the density of liquid water [M L−3], ρvs is the saturated vapor
density [M L−3], M is the molecular weight of water [M mol−1](= 0.018015 kg mol−1), g is the gravitational acceleration [L T−2] (=9.81 m s−2), R is the universal gas content [M L2 T−2 mol−1 K−1] (=8.341 J mol−1 K−1), Tabs is the absolute temperature [K], Hr is therelative humidity, and η is the enhancement factor. A formulation
derived by Cass et al. (1984) and Campbell (1985) was used for η:
3
⎧ ⎡⎛
⎫
θL
2.6 ⎞ θL ⎤ ⎪⎪η = a + 3 − (a − 1)exp ⎨− ⎢⎜1 +
⎥
⎟
⎬
⎜
θs
f c ⎟⎠ θs ⎥⎦ ⎪⎪⎩ ⎢⎣⎝⎭[11]

where fc is the mass fraction of clay in the soil (dimensionless) (ﬁxed at
0.02 in this study) and a is a constant. As discussed below in relation to
Fig. 11, we assumed that a = 8 for the Hamaoka dune sand, based on
the study of Cass et al. (1984) on lysimeter sand. The vapor diffusivity
in soil can be derived from the diffusivity of water vapor in air, Da [L2 T−1],
multiplied by the tortuosity factor in the gaseous phase, τa (Millington and
Quirk, 1961), and the volumetric air content, θair [L3 L−3]:

θ7/3
D = τa θair Da = air2 θair Da
θs

[12]

Selected parameters used in the numerical simulations are shown in
Table 1.

Heat Transport
The principle of heat conservation is described as (de Vries, 1958)
∂CpT
∂t

+ L0

∂θv ∂ ⎡
∂q T
∂q
∂q T
∂T ⎤
= ⎢λ (θL ) ⎥ − Cw L − Ca v − L0 v [13]
∂t
∂z ⎣
∂z ⎦
∂z
∂z
∂z

where L0 (= Lwρw) is the volumetric latent heat of vaporization ofwater (J m−3), Lw (= 2.501 × 106 − 2369.2T) is the latent heat of vaporization of water (J kg−1), and λ(θL) is the soil thermal conductivity[M L T−3 K−1] described as (Chung and Horton, 1987)

709

Table 1. Parameters used in numerical simulations.
Parameter
Surface tension of soil water (γ), g s−2
Density of liquid water (ρw), kg m−3
Saturated vapor density (ρvs), kg m−3
Diffusivity of vapor in air (Da), m2 s−1
Relative humidity (Hr)

Equation†

Reference

Hillel (1971)
γ = 75.6 − 0.1425T − 2.38 × 10−4T2
Hillel (1971)ρw = 1 − 7.37 × 10−6(T − 4)2+3.79 × 10−8(T − 4)3
ρvs = exp(31.37 − 6014.79Tabs−1 − 7.92 × 10−3Tabs)Tabs−1 × 10−3 Campbell (1985)
Campbell (1985)
Da = 2.12 × 10−5(Tabs/273.15)2
Philip and de Vries (1957)
Hr = exp(hMg/RTabs)

Latent heat of vaporization of water (Lw), J kg−1 Lw = 2.501 × 106 − 2369.2T

Monteith and Unworth (1990)† T, temperature; Tabs, absolute temperature; h, pressure head; M, molecular weight of water; g, gravitational acceleration; R, universal gas constant.

λ (θL ) = b1 + b2 θL + b3θ0.5
L[14]

where b1, b2, and b3 are constants (for sands, b1 = 0.228 W m−1 K−1,
b2 = −2.406 W m−1 K−1, and b3 = 4.909 W m−1 K−1, as implemented
in the HYDRUS-1D code). The two terms on the left-hand side of
Eq. [13] represent changes in the energy content and the latent heat ofthe vapor phase. The terms on the right-hand side of Eq. [13] represent soil heat ﬂow by conduction, the convection of sensible heat withﬂowing water, the transfer of sensible heat by diffusion of vapor, andthe transfer of latent heat by diffusion of vapor. The volumetric heat
capacity of the soil, Cp [M L−1 T−2 K−1] is deﬁned as the sum of thevolumetric heat capacities of the solid Cn (= 1.92 MJ m−3 K−1), liquid
Cw (= 4.18 MJ m−3 K−1), and air Ca (= 6.3 kJ m−3 K−1) phases multiplied by their respective volumetric fractions θ (de Vries, 1958):

Cp = Cn θn + Cw θL + Ca θa[15]

where θn is the volumetric fraction of the solid phase [L3 L−3]. The
third term is usually neglected because it is signiﬁcantly smaller thanthe other two terms.

Initial and Boundary Conditions
Initial and boundary conditions were based on the condensation
experiment. The initial volumetric water content was assumed to beuniform throughout the soil column (i.e., θL = 0.0045). The corresponding initial pressure head, calculated from the water retention
curve (Eq. [5]), was approximately −7 × 105 cm. The upper boundary
condition for water ﬂow was given by the total incoming water ﬂux:

qTotal (0, t ) = (qLh + qLT + qvh + qvT ) z =0(t > 0 ) [16]The total surface ﬂux, qTotal(0,t), was assumed to be equal to thechange in weight of the entire column, thus reﬂecting the condensation of water vapor in the column. We ﬁtted a quadratic curve to theobserved increase in weight of the column (cumulative condensation),and determined qTotal(0,t) from the time derivative of this curve.
It was assumed that there was no water ﬂux at the lower boundary describing the closed end of the column:

qTotal (−10, t ) = (qLh + qLT + qvh + qvT ) z =−10 = 0(t > 0 ) [17]The sum of the liquid water ﬂux (the sum of qLh and qLT) and thevapor ﬂux (the sum of qvh and qvT) must be zero to satisfy the no-ﬂux
boundary condition.
Since only measured temperature proﬁles for the 6th and 16thdays were available from Miyazaki (1976), the initial temperature wasassumed to be a constant room temperature of 23.5°C throughout
the soil proﬁle. Dirichlet boundary conditions with speciﬁed temperatures were considered for heat transport calculations. As it is unlikelythat the boundary temperatures of the soil are equal to imposed temperatures (Farlow, 1993, p. 19–26), the upper and lower boundary
temperatures were set equal to the average temperatures observed atthose locations nearest to both boundaries, i.e., 36.3°C at the soil
surface and 23.5°C at the soil bottom.
Governing Eq. [2] and [13] describing coupled movement ofliquid water, heat, and water vapor were solved numerically using the
HYDRUS-1D code that included vapor transport (Šimůnek et al.,2008a,b). The soil column was discretized uniformly into ﬁnite elements of 0.1 cm. The time step was allowed to vary between the initialand maximum allowed time steps of 1 × 10−10 and 0.01 d, respectively.

Hydraulic Parameters

Figure 1a shows water retention
curves for the Hamaoka dune sandmeasured using both drainage andimbibition processes. Since Miyazaki(1976) did not provide the water retention curve, both processes were observed using a hanging water column.
A time domain reﬂectometry (TDR)
probe (three 6.5-cm-long metallic
rods of 0.15-cm diameter with 0.5-cm
separations) and a tensiometer (having
a 1.5-cm-long and 0.6-cm-diameterporous cup) with a pressure transducer were installed horizontally in asoil sample with a height of 2 cm andaninternal diameter of 10 cm. After
Fig. 1. (a) The soil water retention curve of the Hamaoka dune sand, with the Fayer and Simmons (1995)
model fitted to the imbibition data; and (b) unsaturated hydraulic conductivity functions evaluated using draining from saturation to a pressure
pore-connectivity coefficient (l ) values of 0.5, 4.91, 6.00, and 6.17.
head of −60 cm, the imbibition pro710

SSSAJ: Volume 73: Number 3 • May–June 2009

cess was subsequently performed by connecting the soil sample witha Mariotte tank and returning it to saturation. Finally, the drainage
curve was measured and the ﬁnal water content at −60 cm was determined gravimetrically. The TDR measurements were then calibrated
based on the ﬁnal water content and the cumulative amount of drainage. Using this calibration, dielectric constants measured during theimbibition process were then converted to volumetric water contents.
We used the water content values estimated by TDR in Fig. 1a. The
water retention curve for soils with lower water contents was thenmeasured using the pressure-plate apparatus for pressure heads lower
than −100 cm (see Fig. 1a) and assuming that hysteresis could be
ignored for these low pressure heads.
Since the Hamaoka dune sand displayed a relatively importanthysteretic behavior and because the condensation experiment represented a wetting process during the entire period, the imbibitionwater retention curve was used in the calculations. The Fayer model(Eq. [5]) was ﬁtted to the observed data, resulting in values of θa =0.027, θs = 0.325, α = 0.0656 cm−1, and n = 4.71. Note that θs
was about 0.06 smaller than the porosity value of 0.385 calculated
from bulk and particle densities. This was probably partly due to theinevitably entrapped air during the drainage and imbibition processes.
Furthermore, since we estimated θs based on the ﬁnal water content at
h = −60 cm and the cumulative drainage amount, θs might be slightly
underestimated. Since the condensation experiment was performed in
the pressure head range below −20 cm, we believe that these possible
experimental errors near saturation did not affect our calculations.
Figure 1b shows the relation between the pressure head, h, andthe unsaturated hydraulic conductivity, KLh, derived by substituting
Fayer’s water retention curve model (Eq. [5]) into Mualem’s model(Eq. [7]). The observed value of the saturated hydraulic conductivity(Ks = 34.6 m d−1) and the value of the pore-connectivity coefﬁcient
suggested by Mualem (l = 0.5) were used to obtain the unsaturated
hydraulic conductivity function. The gradient, dlogKLh/dlogh, is larger for the higher pressure head range (−10 ≥ h ≥ −100 cm) than for the
lower pressure head range (h < −100 cm).
When the a parameter in the enhancement factor is assumed tobe known (i.e., equal to 8 (Eq. [11]), the pore-connectivity coefﬁcientl is the only remaining unknown parameter (when not assumed to be0.5). Therefore, the value of l was optimized using an inverse analysis(e.g., Šimůnek et al., 1998), during which the numerical solution was
ﬁtted to volumetric water contents observed in the vapor condensation experiment. Parameters θa, θs, α, n, and Ks were kept constant atthe values discussed above during the inverse analyses. The objectivefunction, Φ, was deﬁned asnjniΦ = ∑∑ ⎡⎣θL* (z j , ti )− θL (z j , ti )⎤⎦2

Fig. 2. Enhancement factors for different values of the a parameter(Eq. [11]).

water content proﬁles and water ﬂuxes change when different valuesof the enhancement factor (with a = 5, 8, or 15 or η = 1) were used.

RESULTS AND DISCUSSION
Unsaturated Hydraulic Conductivity
The liquid water content proﬁles were ﬁrst calculated withthe original values of the pore-connectivity coefﬁcient (l = 0.5)
suggested by Mualem (1976) and the enhancement factor (a= 8) suggested by Cass et al. (1984). Figure 3 compares thecalculated and observed water content proﬁles. Similar to theobserved water contents, the calculated water contents also increased from the bottom due to vapor condensation; however,the shape of the simulated water content proﬁles was quite different from those observed. It was found that the simulated water contents were much smaller than the observed values near[18]

j =1 i =1

where θL*(zj, ti) and θL(zj, ti) were the observed and calculated watercontents at time ti for the spatial coordinate zj, respectively, ni was thenumber of measurement times, and nj was the number of measurement locations. Minimization of the objective function, Φ, was accomplished using the Levenberg–Marquardt nonlinear minimization
method (Marquardt, 1963).
The impact of the enhancement factor, η, on the optimized l parameter in the unsaturated hydraulic conductivity, KLh, was evaluatedby performing the model calibration with different values of the a parameter (= 5, 8, or 15) in Eq. [11] (Fig. 2). We also demonstrated how

SSSAJ: Volume 73: Number 3 • May–June 2009

Fig. 3. Comparison of observed water content profiles with those
calculated assuming a pore-connectivity coefficient l of 0.5 and an aparameter in the enhancement factor of 8.
711

relation between the pressure head,h, and the unsaturated hydraulic conductivity, KLh, is shown in
Fig. 1b. Higher values of l lead toa faster decrease in KLh with decreasing h. Figure 4a shows watercontent proﬁles calculated with an
optimized value of l = 6.17 (for a =8). In this case, the calculated watercontents agreed very well with theobserved water contents. Figure 4b
shows the corresponding pressure
head proﬁles. While pressure heads
increased signiﬁcantly from the initialvalue (h = −7 × 105 cm) in the lower
part of the soil column (−30 ≤ h ≤ −20cm), the condensation process caused
Fig. 4. (a) Comparison of observed water content profiles with those calculated assuming a porethem to remain very low (−1 × 105 ≤connectivity coefficient l of 6.17 and an a parameter in the enhancement factor of 8; and (b) simulatedh ≤ −3 × 104 cm) in the upper part.
pressure head profiles.
Since similar results were obthe bottom of the soil column. For example, while the calculattained with the other two combinaed θL was 0.03, the observed value was 0.096 at a depth of 9.5
tions of parameters a and l, these two parameters were mucm after 10 d. On the other hand, calculated water contents
tually correlated and could not be estimated simultaneously
were overestimated in the middle of the column (e.g., calculatfrom available data, i.e., water content proﬁles only. For a given
ed θL = 0.025 vs. observed θL = 0.0075 at a depth of 5 cm after
experimental setup, only one of these two parameters can be10 d). This disagreement between observations and calculations
optimized, while the other parameter has to be estimated indemight be the result of either overestimation of the upward liquid
pendently and ﬁxed during calibration. Notice, however, thatwater ﬂow or underestimation of the downward vapor ﬂow.
differences between the optimized hydraulic conductivity funcSince our focus was mainly on the uncertainty in the untions for different enhancement factors were relatively small
saturated hydraulic conductivity, KLh, for lower water contents,
and probably below the precision of the available measurement
we ﬁxed the water vapor ﬂow parameters (i.e., a = 5, 8, andtechniques for low water contents, for which hydraulic con15 in the enhancement factor) and attempted to decrease theductivity measurements are extremely difﬁcult to do.
upward liquid water ﬂow by adjusting KLh through changing
The agreement between the measured and calculated
the pore-connectivity coefﬁcient, l. Optimized values of l uswater contents indicates that the relation of KLh and h in theing the inverse analysis of the observed water contents were
Hamaoka dune sand at lower water contents can be described
6.00, 6.17, and 4.91 for a = 5, 8, and 15, respectively. The
well using Mualem’s (1976) pore-size distribution model with
the Fayer water retention curve model. The KLh(h) function
shows different gradients, dlogKLh/dlogh, for higher and lower
pressure heads. This shape for the KLh(h) function is similar
to that proposed by Tuller and Or (2001), which considers
ﬁlm ﬂow for lower water contents as well as pore water ﬂow
for higher water contents. While the sand retains water in soil
pores by capillary forces at higher water contents, at lower water contents water is retained on soil particles as water ﬁlms.
These different retention mechanisms may alter the dependency of the hydraulic conductivity on the pressure head. A
similar argument was used for solute mixing in terms of the hydrodynamic dispersion in the dune sand (Toride et al., 2003).
Solutes mix well in larger pores with higher pore-water velocities and higher water contents. On the contrary, slow solute
mixing, primarily by transverse diffusion, becomes dominant
for lower water contents.
Figure 5 shows the unsaturated hydraulic conductivity functions, KLh(h), calculated using the van Genuchten–
Fig. 5. Relations between the unsaturated hydraulic conductivity andMualem (VG) model (van Genuchten, 1980) with various
the pressure head described using the van Genuchten–Mualem model
valuesof l (−1 ≤ l ≤ 5). Although the VG model shows differ(VG), KLh = Ks Sel [1 − (1 − Se1/m)m]2, for various values of the poreent values of dlogKLh/dlogh for different values of l, it neverconnectivity factor l (Ks is the saturated hydraulic conductivity, Se is
resembles the shape of the Fayer model, which has different
the effective liquid saturation, and m is an empirical shape parameter),values of dlogKLh/dlogh for different ranges of pressure heads.
and the Fayer model with l equal to 6.17.
712

SSSAJ: Volume 73: Number 3 • May–June 2009

This indicates that the Fayer model not only modiﬁes the VG
model for the soil water retention curve, but it also has a profound effect on the unsaturated hydraulic conductivity function for lower water contents.
Since the numerical model that uses the unsaturated hydraulic conductivity function KLh(h) with the pore connectivity parameter l = 6.17 that was optimized for the enhancement
factor with a = 8 (Fig. 1b and 5) describes the experimental data
(Fig. 4) well, the following discussion will be based on calculations with this unsaturated hydraulic conductivity function. It
needs to be emphasized here, however, that the optimized unsaturated hydraulic conductivity function may not be uniqueand could not be independently validated against directly measured hydraulic conductivities at low water contents.

Temperature Profiles
Figure 6 presents the simulated and observed temperature
proﬁles obtained during the experiment. A linear decrease intemperature between the upper and lower boundary values wasreached after only 0.1 d. Observed temperatures after 6 and16 d closely followed this linear distribution, except for small
convex deviations at a depth of 9 cm after 6 d and 7 cm after
16 d. This small deviation from the linear proﬁle can be explained by the difference in thermal conductivities between thelower and upper part of the column as a result of their different
water contents. For example, the convex deviation occurred at
a depth of 7 cm after 16 d, which agreed well with the depthof the moisture front (Fig. 4a). The thermal conductivity calculated using Eq. [14] for the lower part of the column (θL >0.075) was >1.4 W m−1 K−1, while it was <0.7 W m−1 K−1 forthe upper part of the column (θL < 0.01).

Liquid Water and Water Vapor Fluxes

Fig. 6. Comparison of observed and calculated soil temperature profiles.

por ﬂux, qvT, due to the large downward temperature gradient(Fig. 6). The isothermal vapor ﬂux, qvh, was small except near
the moisture front (in the 9-cm depth after 6 d and the 7-cm
depth after 16 d), where upward vapor ﬂux was calculated dueto the large positive pressure head gradient (Fig. 4b). Since qvh
was much smaller than qvT, the actual vapor ﬂux, qv, was downward after both 6 and 16 d.
Upward water movement occurred mainly as isothermal
liquid water ﬂux, qLh, due to the upward pressure head gradient, which developed at the bottom of the soil column as aresult of vapor condensation. The isothermal liquid water ﬂux
had a peak at both times near the moisture front due to theextremely large pressure head gradients. The isothermal liquidwater ﬂux was also large at the bottom of the column after 16 d
(0.067 cm d−1) because of the large unsaturated hydraulic conductivity, KLh, for high water contents. As the water content
gradually increased at the column bottom and the moisture
front moved upward, so did the maxima and absolute valuesof qLh and qvT. The thermal liquid water ﬂux, qLT, showed

Once the pressure head and temperature proﬁles were obtained, it was possible to evaluate the liquid and vapor ﬂuxes
due to pressure head (isothermal) and temperature (thermal)
gradients using Eq. [1]. Figure 7 shows the calculated total waterﬂux proﬁles, as well as their four components, after 6 and 16 d.
Positive values represent upward movements, while negative values represent downward movements.
The total water ﬂux, qTotal,
was always directed downward(negative) throughout the column from the warm end to thecold end, similarly as in previous studies of nonisothermal ﬂow
in column experiments (Jones
and Kohnke, 1952; Nassar et al.,1992a). While the total water ﬂux,
qTotal, was large in the upper partof the soil column above the moisture front (−0.017 cm d−1 after 6 dand −0.02 cm d−1 after 16 d), it
was much smaller in the lower partof the column with higher watercontents (Fig. 4a), and it was zeroat the bottom of the column. The Fig. 7. Calculated profiles of the thermal liquid and vapor (q and q , respectively) and isothermalLTvT
downward water movement oc- liquid and vapor (qLh and qvh, respectively) fluxes and the total water flux (qTotal) at 6 d (left) and 16 dcurred mainly as thermal water va- (right) (positive fluxes are directed upward).
SSSAJ: Volume 73: Number 3 • May–June 2009713

observed diffusion ﬂux of water vapor
into the column during calculations.
Figure 8a presents the calculated surface pressure heads during the experiment. Although not
clearly visible in Fig. 8a, at the beginning of the experiment the surface pressure head ﬁrst decreased
quickly from the initial pressure
head of −7 × 105 cm, and then increased because of the rapid change
in the surface temperature from theinitial value of 23.5°C to the imposed boundary value of 36.3°C.
After that, the surface pressure
Fig. 8. Calculated (a) surface pressure heads and (b) surface thermal liquid and vapor (qLT and qvT, head started increasing steadily unrespectively) and isothermal liquid and vapor (qLh and qvh, respectively) fluxes and the total water flux til it reached a value of around −1.5
(qTotal) as a function of time.
× 105 cm. The increasing pressure
head reﬂects vapor condensation at
downward water movement near the bottom of the column.
the soil surface. As the amount of condensed water was relaAlthough qLT was small after 6 d, it increased after 16 d
tively small, water content increases were barely noticeable.
(−0.06 cm d−1 at the bottom) due to a larger thermal hydraulicFigure 8b shows the four components of the surface ﬂux
conductivity, KLT (Eq. [8] with the large KLh). Because qLh
qTotal(0,t) as a function of time. Liquid water ﬂuxes (qLh andwas larger than qLT at every depth, the actual liquid water ﬂux,
qLT) were signiﬁcantly smaller (and in fact negligible) thanqL, was always upward. The value of qL was negligibly small atvapor ﬂuxes (qvh and qvT) because the hydraulic conductividepths above the moisture front after both 6 and 16 d.
ties for liquid water were several orders of magnitude smaller
Miyazaki (1976) calculated the increase in water contents(KLh ≈ 10−19 cm d−1 and KLT ≈ 10−15 cm2 K−1 d−1) thanat the bottom of the columns by evaluating the water vapor ﬂux
vapor hydraulic conductivities (Kvh ≈ 10−6 cm d−1 anddue to the temperature gradient using the Philip and de Vries
KvT ≈ 10−1 cm2 K−1 d−1). Although the direction of qvh wastheory and the Penman theory (Penman, 1940), which was thechanging at the beginning of the experiment, due to oscillaoriginal theory of water vapor ﬂow without the enhancement
tions of the surface pressure head in time (Fig. 8a), downward
factor. The water content calculations signiﬁcantly overestiqvT and upward qvh dominated the overall water ﬂux qTotal(0,mated the observed values, however, because of the difﬁcultyt), and deﬁned the water vapor diffusion into the column.
of calculating the simultaneous upward liquid ﬂux resulting
from the water content increase. Our numerical simulations ofCondensation and Evaporation Rates
the coupled movement of water vapor, liquid water, and heat
The evaporation rate (E) at each internal discretization
could consider simultaneously all four ﬂux components andnode can be calculated with a discretized form of Eq. [3]:thus lead to better agreement between simulated and observedj −1/2
θLj i − θLj −i 1 qLji−−1/2water contents.
j
1/2 − qLi +1/2
[19]

Ei = −

Surface Boundary Condition
Both Milly (1984) and Saito et al. (2006) evaluated surface
water ﬂux. The former used meteorological variables, while thelatter utilized the surface energy balance. Since evaluating water movement inside the soil column was the primary purposeof this study, the observed amount of water vapor diffused into
the column, qTotal(0,t), was determined from weight changes
in the entire soil column and used as the surface boundary condition for water ﬂow. We assumed that qTotal(0,t) was equal to
the sum of liquid and vapor ﬂuxes passing through the surface
boundary, as described by Eq. [16]. The four components ofthe total ﬂux at the soil surface were then determined using thepressure head and temperature gradients, as described by Eq.
[1]. Each ﬂux may have a different magnitude and direction
depending on the pressure head and temperature gradientsat the soil surface. Since the temperature in the soil column
quickly established steady-state conditions (during about 0.1 d
as shown in Fig. 6), the pressure head at the soil surface was the
only other variable that could adjust qTotal to be equal to the
714

Δt

+

Δz

where the subscript i represents the position of a lattice point
and the superscripts j and j − 1 denote the current and previous
time points, respectively.
The evaporation rate at locations i = 1 + 1/4 and i = N −1/4 can be evaluated using a method similar to Eq. [19], wherei = 1 and N represent the lower and upper boundary nodal
points, respectively:

E1j+1/4 = −

j
j −1
θL1
q j −1/2 − qL1j −+1/2+1/4 − θ L1+1/4
1/2
+ L1
Δt
Δz 2

[20]

ENj −1 4 = −θLj N −1/4 − θLj −N1−1/4 qLjN−1/2 − qLjN−1/2
−1/2
+
Δt
Δz 2

[21]

Values at i = 1 + 1/4 and i = N − 1/4 (e.g., θL1+1/4) can beobtained using a linear interpolation of adjacent node values.
Note that qL1 does not have to be zero because, as describedby Eq. [17], only the total boundary ﬂux is assumed to be zero.
SSSAJ: Volume 73: Number 3 • May–June 2009

The evaporation rate at the surface
boundary, E1, is evaluated by linearly extrapolating E1+1/4 and E2
toward the surface. A similar linear extrapolation from EN−1/4 andEN−1 is made to estimate EN.
Figure 9a displays evaporation
rate proﬁles calculated using Eq.
[19–21]. Note that negative evaporation represents vapor condensation. The maximum evaporation
rate was located at a depth of 9 cm
after 6 d, and moved upward to adepth of 7 cm after 16 d. This localized evaporation resulted in the
maximum downward thermal vapor ﬂux, qvT, behind the moisture Fig. 9. Calculated profiles of (a) the evaporation rate and (b) the relative humidity.
front as shown in Fig. 7. Figure 9b
shows relative humidity proﬁles
the observed water content proﬁles as well as those calculated(Table 1) at 6 and 16 d. The relative humidity decreased fromwith different enhancement factors after 30 d with an optialmost unity near the bottom of the column to 0.98 at themized value of l = 6.17. In calculations without the enhancemoisture front. The location of the maximum evaporation rate
ment factor (η = 1), calculated water contents were underesticorresponded exactly with the position of the abrupt decrease
mated in the lower part of the column and overestimated inin the relative humidity. As shown in Fig. 4b, a large pressure
the upper part. More vapor condensed in the upper part of thehead drop occurred at the moisture front (h = −2 × 104 cm),
column due to underestimation of the vapor ﬂow. This simularesulting in a decrease in the relative humidity.
tion clearly demonstrates the importance of considering the
Vapor condensation started just below the maximum
enhancement factor. Calculated water content proﬁles were
evaporation peak, where the relative humidity increased to alsimilar when different values of the a coefﬁcient (a = 5, 8, andmost unity. Vapor had to condensate at this position because15) were used to calculate the enhancement factor (Fig. 11).
the air could not hold all the water vapor coming from the upFigures 12a and 12b show proﬁles of the liquid water and water
per location. As shown in Fig. 6, vapor moved farther down andvapor ﬂuxes and the evaporation rate, calculated with a = 8 andeventually condensed at the cold end. The maximum condensa15. Both the downward vapor ﬂux and the upward liquid watertion (minimum evaporation) rate occurred at the bottom of theﬂux increased as the enhancement factor increased, and so didcolumn: −0.123 d−1 after 6 d and −0.085 d−1 after 16 d (Fig. 9a).
both the maximum evaporation and condensation rates. The
The total amount of vapor condensation in an entire soil
amount of water circulating in the soil column increased ascolumn can be obtained by integrating the condensation rate
the enhancement factor increased. These two opposing effects
proﬁle (negative evaporation rate in Fig. 9a) with depth and(evaporation and condensation) compensated for each other,time. Figure 10 shows the cumulative amount of condensation
as a function of time. The total quantity of condensation in thesoil proﬁle includes not only the condensation at the bottom,
but also that in the rest of the soil proﬁle. The condensation at
the bottom dominated initially and accounted for 40% of thetotal condensation during the ﬁrst 5 d. Condensation inside
the column subsequently increased and reached about 75% of
the total condensation (25% at the bottom) after 20 d. The
observed cumulative amount of water vapor diffused into thecolumn was also plotted in Fig. 10. Note that the total amountof condensation in the column was greater than the cumulative amount of water vapor diffused into the column due tointernal evaporation (Fig. 9a). In other words, the difference
between the total amount of condensation and the cumulative
amount of diffusion into the column was equal to the amountof evaporation inside the column. The amount of internal
evaporation accounted for around 25% of the total condensation, regardless of the time.

Enhancement Factor

The impact of the enhancement factor, η, can also bedemonstrated using water content proﬁles. Figure 11 showsSSSAJ: Volume 73: Number 3 • May–June 2009

Fig. 10. Cumulative condensation amounts in the whole soil column andat the bottom of the soil column compared with the cumulative amountof water vapor diffused through the surface as a function of time.
715

ment of water vapor and liquid water in the sand column were
evaluated using the Philip and de Vries (1957) model. The coupled transport of water vapor, liquid water, and heat were calculated using the HYDRUS-1D code that included vapor ﬂow.
The numerical model was calibrated against experimental data by optimizing the pore-connectivity coefﬁcient, l, ofthe unsaturated hydraulic conductivity function for different
values of the enhancement factor η. The Fayer and Simmons(1995) model, describing soil hydraulic properties with anestimated l parameter of 6.17 (for the enhancement factor
with a = 8), provided a very good agreement between simulated and observed water contents. The rate of decrease in theunsaturated hydraulic conductivity with decreasing pressure
head, described by the Fayer model, differed for high and low
water contents, reﬂecting capillary pore water ﬂow and ﬁlm
ﬂow, respectively. It was concluded that the Fayer model withthe calibrated l parameter of Mualem’s pore-size distribution
model worked well for simulating water ﬂow in sandy soils
at low water contents. Due to the lack of directly measured
Fig. 11. Water content profiles after 30 d calculated using various
values of the enhancement factor (a = 5, 8, 15 and η = 1).
unsaturated hydraulic conductivities, however, the calibrated
hydraulic conductivity function could not be independently
making the impact of the enhancement factor less important
validated. Mutual correlation between the pore-connectivity
for the water content proﬁle in the condensation process.
coefﬁcient and the enhancement factor also raises a question of
We have also evaluated the impact of η on the estimation ofuniqueness of the optimized unsaturated hydraulic conductivthe pore-connectivity coefﬁcient, l. Optimized l values were 6.00,
ity function.
6.19, and 4.91 for the enhancement factor a parameters of 5, 8,
Additionally, we quantitatively evaluated the four compoand 15, respectively. As shown in Fig. 1b, there is only a relatively
nents of the total water ﬂux: liquid and vapor ﬂuxes driven
small range of uncertainty in the unsaturated hydraulic conductivby pressure head and temperature gradients. The evaporation
ity function due to differences in η for lower pressure heads. As
and condensation rates inside the soil column were describeddiscussed above, however, these two parameters are mutually highfurther by analyzing these water ﬂuxes. Water vapor that enly correlated for given experimental data and only one of them cantered the soil column from the hot surface moved downward
be optimized independently. Additional data need to be measured
due to the temperature gradient and condensed at the coldso that both the pore-connectivity parameter l and the enhancebottom end of the column. Liquid water subsequently moved
ment factor η can be estimated simultaneously.
upward due to the pressure head gradient and evaporated at
the moisture front where the relative humidity decreased fromCONCLUSIONS
1 to 0.98. Vapor generated by evaporation then moved downWe have numerically analyzed the vapor condensation
ward together with vapor coming from the soil surface. A cerexperiments reported by Miyazaki (1976). These experiments
tain amount of vapor condensed just below the location withinvolved vapor diffusion and condensation in laboratory sandmaximum evaporation at the moisture fronts, and the maxicolumns due to an imposed temperature gradient. The movemum condensation occurred at the cold bottom of the column(25–40% of the total condensationin the column). Liquid water andwater vapor circulated between
the bottom of the column and the
moisture front, accompanied bycondensation and evaporation processes in the sand column.
The impact of the enhancement factor on vapor ﬂuxes inside
the soil column was also demonstrated. Downward vapor ﬂux was
signiﬁcantly underestimated whenthe enhancement factor was neglected. Calculated water content
proﬁles were very similar when thea parameter from the function of
the enhancement factor was >5.
Higherenhancement factors proFig. 12. A comparison of profiles of (a) liquid water and water vapor fluxes (qL and qv, respectively) andduced larger upward liquid water(b) evaporation rates simulated using a parameter values of 8 and 15 in the enhancement factor.
716

SSSAJ: Volume 73: Number 3 • May–June 2009ﬂuxes and evaporation rates, as well as downward vapor ﬂuxes
and condensation rates. These two opposing effects compensated for each other, making the water content proﬁles in thecondensation experiment less sensitive to the enhancement
factor. This indicates that water content proﬁles do not provide enough information to simultaneously estimate both theunsaturated hydraulic conductivity function and the enhancement factor. Additional data are needed so that both the poreconnectivity parameter l and the enhancement factor η can be
estimated simultaneously.

ACKNOWLEDGMENTS
We would like to thank Dr. Tsuyoshi Miyazaki from the Universityof Tokyo for his kindness in offering us his experimental data as wellas the Hamaoka dune sand that was used in his column experiments
for additional measurements of the water retention curve. This study
was supported in part by the Terrestrial Sciences Program of the Army
Research Ofﬁce (Terrestrial Processes and Landscape Dynamics andTerrestrial System Modeling and Model Integration), by the National
Science Foundation (NSF) Biocomplexity programs no. 04-10055
and NSF DEB 04-21530.

REFERENCES
Cahill, A.T., and M.B. Parlange. 1998. On water vapor transport in ﬁeld soils.
Water Resour. Res. 34:731–739.
Campbell, G.S. 1985. Soil physics with BASIC. 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 ﬂow with a partial
surface mulch. Water Resour. Res. 23:2175–2186.
de Vries, D.A. 1958. Simultaneous transfer of heat and moisture in porous
media. Trans. Am. Geophys. Union 39:909–916.
Eching, S.O., and J.W. Hopmans. 1993. Optimization of hydraulic functions
from transient outﬂow and soil water pressure data. Soil Sci. Soc. Am. J.
57:1167–1175.
Farlow, S.J. 1993. Partial differential equations for scientists and engineers.
Dover, New York.
Fayer, M.H., and C.S. Simmons. 1995. Modiﬁed soil water retention functions
for all matric suctions. Water Resour. Res. 31:1233–1238.
Fujimaki, H., and M. Inoue. 2003. A transient evaporation method for determining
soil hydraulic properties at low pressure. Vadose Zone J. 2:400–408.
Hillel, D. 1971. Soil and water: Physical principles and processes. Academic
Press, New York.
Hopmans, J.W., J. Šimůnek, and K.L. Bristow. 2002. Indirect estimation of soil thermal
properties and water ﬂux using heat pulse probe measurements: Geometry anddispersion effects. Water Resour. Res. 38(1):1006, doi:10.1029/2000WR000071.
Inoue, M., J. Šimůnek, J.W. Hopmans, and V. Clausnitzer. 1998. In situ
estimation of soil hydraulic functions using a multistep soil-water
extraction technique. Water Resour. Res. 34:1035–1050.
Jones, H.E., and H. Kohnke. 1952. The inﬂuence of soil moisture tension onvapor movement of soil water. Soil Sci. Soc. Am. Proc. 16:245–248.
Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity:
Laboratory methods. p. 687–734. In A. Klute (ed.) Methods of soil
analysis. Part 1. Physical and mineralogical methods. 2nd ed. Agron.
Monogr. 9. ASA and SSSA, Madison, WI.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of nonlinear
parameters. SIAM J. Appl. Math. 11:431–441.
Millington, R.J., and J.M. Quirk. 1961. Permeability of porous solids. Trans.
Faraday Soc. 57:1200–1207.
Milly, P.C.D. 1984. A simulation analysis of thermal effects on evaporation.

SSSAJ: Volume 73: Number 3 • May–June 2009

Water Resour. Res. 20:1087–1098.
Miyazaki, T. 1976. Condensation and movement of water vapor in sand undertemperature gradient. (In Japanese.) Trans. Jpn. Soc. Irrig. Drain. Rural
Eng. 61:1–8.
Monteith, J.L., and M.H. Unworth. 1990. Principles of environmental physics.
2nd ed. Edward Arnold, London.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12:513–522.
Nassar, I.N., A.M. Globus, and R. Horton. 1992a. Simultaneous soil heat andwater transfer. Soil Sci. 154:465–472.
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 inporous media: I. Theoretical development. Soil Sci. Soc. Am. J. 56:1350–1356.
Nassar, I.N., and R. Horton. 1997. Heat, water, and solute transfer inunsaturated porous media: I. Theory development and transport
coefﬁcient evaluation. Transp. Porous Media 27:17–38.
Nassar, I.N., R. Horton, and A.M. Globus. 1992b. Simultaneous transfer ofheat, water, and solute in porous media: II. Experiment and analysis. Soil
Sci. Soc. Am. J. 56:1357–1365.
Nimmo, R.J. 1990. Experimental testing of transient unsaturated ﬂow theory at
low water content in a centrifugal ﬁeld. Water Resour. Res. 26:1951–1960.
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. 1996. Two-dimensional model
for water, heat, and solute transport in furrow-irrigated soil: I. Theory.
Soil Sci. Soc. Am. J. 60:1001–1009.
Penman, H.L. 1940. Gas and vapor movement in soil: I. The diffusion ofvapors in porous solids. J. Agric. Sci. 30:437–462.
Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous materials
under temperature gradients. Trans. Am. Geophys. Union 38:222–232.
Romano, N., and A. Santini. 1999. Determining soil hydraulic functions from
evaporation experiments by a parameter estimation approach: Experimental
veriﬁcations and numerical studies. Water Resour. Res. 35:3343–3359.
Rossi, C., and J.R. Nimmo. 1994. Modeling of soil water retention fromsaturation to oven dryness. Water Resour. Res. 30:701–708.
Saito, H., J. Šimůnek, and B.P. Mohanty. 2006. Numerical analysis of coupled water,
vapor, and heat transport in the vadose zone. Vadose Zone J. 5:784–800.
Scanlon, B., K. Keese, R.C. Reedy, J. Šimůnek, and B.J. Andraski. 2003. Variations
in ﬂow and transport in thick desert vadose zones in response to paleoclimatic
forcing (0–90 kyr): Field measurements, modeling, and uncertainties.
Water Resour. Res. 39(7):1179, doi:10.1029/2002WR001604.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA:
A computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. J. Hydrol. 251:163–176.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M.Th. van Genuchten. 2008a. The
HYDRUS-1D software package for simulating the movement of water, heat,and multiple solutes in variably saturated media. Version 4.0. HYDRUS
Softw. Ser. 3. Dep. of Environ. Sci., Univ. of California, Riverside.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008b. Development andapplications of the HYDRUS and STANMOD software packages, andrelated codes. Vadose Zone J. 7:587–600.
Šimůnek, J., O. Wendroth, and M.Th. van Genuchten. 1998. Parameter
estimation analysis of the evaporation method for determining soil
hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.
Toride, N., M. Inoue, and F.J. Leij. 2003. Hydrodynamic dispersion in an
unsaturated dune sand. Soil Sci. Soc. Am. J. 67:703–712.
Tuller, M., and D. Or. 2001. Hydraulic conductivity of variably saturated
porous media. Film and corner ﬂow in angular pore space. Water Resour.
Res. 37:1257–1276.
van Genuchten, M.Th. 1980. A closed-form equation for predicting hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.

717