Published online May 17, 2018

Special Section: Lysimeters in
Vadose Zone Research

Water Distribution in an Arid Zone
Soil: Numerical Analysis of Data
from a Large Weighing Lysimeter
J. Dijkema, J.E. Koonce, R.M. Shillito, T.A. Ghezzehei,
M. Berli,* M.J. van der Ploeg, and M.Th. van Genuchten

Core Ideas• We aimed to advance our mechanistic understanding of desert soilhydrodynamics.
• We calibrated and validated a
HYDRUS-1D software version ondesert soil data from a weighing
lysimeter.
• HYDRUS-1D is an attractive tool for
virtual experiments relevant for arid
ecosystems.

J. Dijkema, J.E. Koonce, R.M. Shillito,
and M. Berli, Division of Hydrologic
Sciences, Desert Research Institute,
Las Vegas, NV 89119; J. Dijkema
and M.J. van der Ploeg, Dep. ofEnvironmental Sciences, Soil Physics
and Land Management Group,
Wageningen
Univ.,
Wageningen,the Netherlands; T.A. Ghezzehei,
Dep. of Life and Environmental
Sciences, Univ. of California, Merced,
CA 95343; M.Th. van Genuchten,
Center for Environmental Studies,
CEA, São Paulo State Univ., UNESP,Rio Claro, SP, Brazil, and Dep. of Earth
Sciences, Utrecht Univ., Utrecht, theNetherlands. *Corresponding author(Markus.Berli@dri.edu).
Received 31 Jan. 2017.
Accepted 24 May 2017.
Citation: Dijkema, J., J.E. Koonce,
R.M. Shillito, T.A. Ghezzehei, M. Berli,
M.J. van der Ploeg, and M.Th. van
Genuchten. 2018. Water distribution
in an arid zone soil: Numerical analysis of data from a large weighing
lysimeter. Vadose Zone J. 17:170035.
doi:10.2136/vzj2017.01.0035© Soil Science Society of America. This is
an open access article distributed under
the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Although desert soils cover approximately one third of the Earth’s land surface, surprisingly little is known about their physical properties and how
those properties affect the ecology and hydrology of arid environments.
The main goal of this study was to advance our understanding of desert
soil hydrodynamics. For this purpose, we developed a process-based component within HYDRUS-1D to describe the moisture dynamics of an aridzone soil as a function of water fluxes through the soil surface. A modified
van Genuchten model for the dry end of the soil water retention curve wasdeveloped to better capture the basic flow processes for very dry conditions. A scaling method was further used to account for variabilities in water
retention because of changes in the bulk density vs. depth. The model wascalibrated and validated using hourly soil moisture, temperature, and mass
data from a 3-m-deep weighing lysimeter of the Scaling Environmental
Processes in Heterogeneous Arid Soils facility at the Desert Research Institute(Las Vegas, NV). Measurements and simulations during a 1-yr period agreed
better under precipitation (wetting) than under evaporation (drying) conditions. Evaporation was better simulated for wet than for dry soil surface
conditions. This was probably caused by vapor-phase exchange processes
with the atmosphere, which were unaccounted for and need to be further
explored. Overall, the model provides a promising first step toward developing a more realistic numerical tool to quantify the moisture dynamics of arid
ecosystems and their role in climate change, plant growth, erosion, andrecharge patterns.
Abbreviations: HDU, heat dissipation unit; SEPHAS, Scaling Environmental Processes in
Heterogeneous Arid Soils; TDR, time-domain reflectometry; TPHP, triple-point heat dissipation probes; WRC, water retention curve.

Desert soils cover approximately one third of the Earth’s land surface (Hare, 1985),
with 15 to 25% of these soils being located in arid climates that have an annual precipitation of <250 mm. Desert soils are known to affect cloud and dust formation, as well
as CO2 fluxes, which affect climate on a global scale (Austin et al., 2004; Huxman etal., 2004; Morgan et al., 2004; Bhattachan et al., 2012). Desert soils also support fragile
ecosystems with flora and fauna that have adapted to limited water resources (Devitt etal., 2011; Kidron and Gutschick, 2013; Moreno-de las Heras et al., 2016). Demands ondesert soils continue to increase because of an increase in the global population, and as aresult of human activities such as agricultural and horticultural operations, construction,
waste disposal, and recreation (Dregne et al., 1991; D’Odorico et al., 2013). Some of thefastest growing cities in the world are now in arid areas; their rapid development is often
coupled with above-average water use per capita (Kenny and Juracek, 2012).
Despite the considerable extent of desert soils, relatively little is known about their moisture dynamics in space and time. This may be because deserts were previously considered
unproductive for forestry and agriculture unless irrigated. Most hydrologic studies ofdesert soils have focused on deep infiltration (i.e., tens to hundreds of meters deep) to assessgroundwater recharge or contaminant transport processes (Gee et al., 1994; Andraski,

Vadose Zone Journal | Advancing Critical Zone Science1997; Constantz et al., 2003; Gee et al., 2005; Scanlon et al., 2005;
Andraski et al., 2014) and have overlooked the near-surface soil
environment (i.e., the top centimeters to meters). However, thenear-surface environment controls the infiltration, runoff, storage, and evapotranspiration of sparse precipitation in arid climates(Koonce, 2016).
The overarching goal of this study was to advance our mechanistic understanding of the moisture dynamics of near-surface desert
soils. We were particularly interested in moisture dynamics asa function of precipitation and evaporation within the top few
meters of an arid soil. To evaluate these dynamics, we comparedmeasured soil moisture and temperature data from a 3-m-deepweighing lysimeter located in the Mojave Desert of southern
Nevada with numerical simulations using a process-based soil
physical model implemented in the HYDRUS-1D software package (Simunek et al., 2016). The weighing lysimeter is part of the
Scaling Environmental Processes in Heterogeneous Arid Soils(SEPHAS) facility in Boulder City, NV, constructed and managed by the Desert Research Institute (Chief et al., 2009). The
lysimeter has been operational since 2008 and has collected >8 yrof data on precipitation, evaporation, soil moisture and temperature, and other relevant environmental parameters. The controlled
conditions of the SEPHAS lysimeter experiment allow us to studyhydrological processes of arid soils under well-defined initial andboundary conditions. The resulting database provides a unique
opportunity to test the validity of conceptual and mathematicalmodels simulating vadose zone water flow and heat transfer processes that are typical of arid and semiarid regions. We used thedata to test the capability of the HYDRUS-1D software package(Šimůnek et al., 2016) to capture the strong coupling of water andheat transport processes in the lysimeter, as well as unsaturated
flow near the lowest extreme of the soil moisture range.

Typic Torriorthent of the Arizo series (Soil Survey Division Staff,
1993) and is primarily vegetated with creosote bush [Larrea tridentate (DC.) Coville] and white bursage [Ambrosia dumosa (A.
Gray) Payne].

Lysimeter Design and Instrumentation
The SEPHAS lysimeters are cylindrical, stainless steel tankswith a diameter of 226 cm and depth of 300 cm (Type 304,
Moore’s Blacksmith Shop) that are placed on truck scales(Model Z-100, Cardinal Scale Manufacturing Company) andhoused underground in individual rooms that are accessible by atunnel (Chief et al., 2009). Figure 1 shows a schematic diagramof the dimensions and design of the lysimeter. The scales havea resolution of 450 g and can detect changes in the lysimeter
mass equivalent to a water column of 0.1 mm. The lysimetertanks were filled with soil from bottom to top in individual lifts(between 2- and 22-cm lift thickness), with each individual lift
being repacked as closely as possible to the bulk densities measured at the Eldorado Valley study site. The lysimeters have been
kept free of vegetation since their installation in 2008. Chiefet al. (2009) provides a comprehensive description of the three
SEPHAS lysimeters.
Only the data from Lysimeter 1 are described here, since thatlysimeter was the focus of our study. Lysimeter 1 has a homogenous texture profile with a fine fraction (<2 mm) consisting of
93.0% sand, 5.5% silt, and 1.5% clay and an average gravel content
of 18.9% by mass. The lysimeter was repacked to the bulk density

66Materials and Methods
Study Site

The study was conducted at the SEPHAS weighing lysimeter
facility in Boulder City (35.96° N, 114.85° W), which is ~40 km
southeast of Las Vegas, NV. The facility is located in the Mojave
Desert at an elevation of 768 m. The average annual precipitation
is 141 mm, whereas average daily temperatures range from 8.1°C in
January to 31.7°C in July (1931–2005) (Western Regional Climate
Center, 2017). The core of the SEPHAS facility consists of threeweighing lysimeter tanks filled with desert soil from a study site in
nearby Eldorado Valley, which is ~8 km southwest of the facility.
Our study focused on only one of the three lysimeters. Soils at the
Eldorado Valley study site developed from volcanic parent material(mainly andesite and rhyolite) that was deposited on a south-facing,shallow-sloped alluvial fan (0–15% slope angle) from the nearby
McCullough and Highland Ranges (Chief et al., 2009). The soil
we used has a been classified as a sandy-skeletal, mixed, thermic

VZJ | Advancing Critical Zone ScienceFig. 1. Schematic diagram of lysimeter dimensions and the distributionof triple-probe heat probe (TPHP) sensors for temperature (T), timedomain reflectometry (TDR) sensors for water content (q), and heatdissipation unit (HDU) sensors for pressure head (h) and temperature(T). Numbers adjacent to the sensor icons indicate the number of
replicates installed at each depth.

p. 2 of 17

precipitation rates were monitored using a tipping bucket rain
gauge (Model TE525WS-L, Campbell Scientific) located ~30 mwest of Lysimeter 1. The rain gauge provided precipitation data at
30-min intervals. This study covered the time period from 1 Oct.
2011 at 0000 h to 30 Sept. 2012 at 2300 h, which resulted in 8784
hourly measurements.

Laboratory Characterization of the Soil
Water Retention Curve

Fig. 2. Profile of the porosity (left panel) and corresponding linear
scaling coefficients (right panel) of the water content (D), pressure
head (A), and hydraulic conductivity (B).

profile shown in Fig. 2a, which reflects the bulk density profile ofthe undisturbed soil at the Eldorado Valley study site. Bulk densities ranged from 1442 to 1798 kg m−3, yielding a porosity range
from 0.24 to 0.31 (as calculated using an average measured particle
density of 2479 kg m−3).
Volumetric moisture contents were measured at depths of 10, 25,
50, 75, 100, 150, 200, and 250 cm using time-domain reflectometry (TDR) probes (Model CS605, Campbell Scientific). Fourprobes per depth were installed at 10, 25, 50, and 100 cm, andtwo probes per depth were installed at 75, 150, 200, and 250 cm.
Soil water pressure heads (matric potentials) were monitored at
depths of 5, 10, 25, 50, 75, 100, 150, 200, and 250 cm using
heat dissipation units (HDUs; Model 229, Campbell Scientific).
Four HDUs per depth were installed at 10, 25, 50, and 100 cm,
and two HDUs per depth at 75, 150, 200, and 250 cm. Soil
temperatures at depths of 5, 10, 25, 50, 75, 100, 150, 200, and250 cm were determined using temperature sensors installed in
the HDUs at the respective depths. Additionally, near-surface
soil temperatures were measured at depths of 12, 24, 36, and48 mm using triple-point heat dissipation probes (TPHP, East 30
Sensors). The arithmetic means of the soil water contents, pressure heads, and temperatures were calculated for each depth, andthe averages were used for the various calculations.
Hourly precipitation and evaporation rates were calculated fromlysimeter mass readings taken at 15-min intervals. Concurrently,

VZJ | Advancing Critical Zone ScienceThe calibration range of the in situ pressure head sensors did
not span the range that was experienced by the lysimeter soil.
Additional laboratory measurements were performed using a dewpoint potentiometer (WP4, Decagon Devices) to provide data forcalibrating the flow model in the dry range. Three disturbed andhomogenized samples of the bulk soil were moistened to ~0.1 kgkg−1 and dried in consecutive steps by exposing them to the ambient laboratory environment for short intervals. After each drying
interval, the samples were equilibrated overnight and their water
potentials determined using the potentiometer, which was set in its
manual mode to provide the most accurate readings. Gravimetricwater contents were determined at each step using the sample mass,
whereas the final mass was determined after oven drying at 105°C
for 24 h. The gravimetric water contents were converted to volumetric fractions using the bulk density values of the pertinent soil
layers in the lysimeter.

Model Description
One-dimensional unsaturated water flow and heat transport processes in the lysimeter were analyzed using the comprehensive
HYDRUS-1D code, which has been successfully used to simulate
vadose zone hydrologic processes of a wide range of soil and environmental problems (Šimůnek et al., 2016), including arid soil (Dixon,
1999). One major advantage of HYDRUS-1D for our study is its
capability to simulate coupled heat and nonisothermal water flow inboth the liquid and vapor phases. Nonisothermal liquid and vapor
flow of water is described below as a sum of isothermal and non-isothermal components (Philip and de Vries, 1957; Saito et al., 2006):¶q T ( h )
=
¶t
¶ é
¶h
¶T ù
ê( K + K vh ) + ( K LT + K vT )
ú - S ( h)¶ z ëê¶z¶ z ûú

[1]

where q T = q + q v [L3 L−3] is the total volumetric water content(i.e., the sum of the liquid and vapor components, respectively),with t being the time [T] and z being the soil depth [L]. The firstpart of Eq. [1] within brackets denotes the isothermal flux driven
by pressure head (h [L]) gradients and includes liquid (K [L T−1])
and vapor (Kvh [L T−1]) phase fluxes. The second part represents
water flow driven by thermal gradients (T [K]) and includes both
liquid (K LT [L2 K−1 T−1]) and vapor (KvT [L2 K−1 T−1]) phase
fluxes. The last term (S) of Eq. [1] denotes local sources and sinks.
The accompanying heat transport equation is given byp. 3 of 17¶q¶T+ L0 v =
¶z
¶z
¶q T
¶q ö
¶ æç ¶ T¶T-C w q-C w v - L0 v ÷÷÷
çç lq
è
¶z
¶z
¶z
¶z
¶z ø

Cp q

[2]

where Cp and Cw [M L−1 T−2 K−1] are the volumetric heat capacities of the porous medium and the pore water, respectively; L0[M L−1 T−2] is the heat of vaporization; l is the effective thermal
conductivity [M L T−3 K−1]; and q and qv (L T−1) are the fluxes
in liquid and vapor form. The flux terms are derived from thebrackets in Eq. [1]. The thermal conductivity of the composite
medium was calculated using an empirical relationship proposedby Chung and Horton (1987):[3]

lq = b1 + b2 q + b3 q 0.5

in which we used the parameter values b1 = 0.228, b2 = −2.406,
and b3 = 4.909, which are typical of a sandy texture, to yield l in
SI units [M L T−3 K−1]. The water retention curve (WRC) of theporous medium was described using the dual-porosity model ofDurner (1994):
S e = (1 - w ) éê1 + (a 1 h) 1 ùú
n

ë

-m1

û

+ w éê1 + (a 2 h) 2 ùú
n

ë

- m2

û

[4]

where S e = (q − q r)/(q s − q r) is effective saturation (with q s
and q r representing the saturated and residual water contents,
respectively), m = 1 − 1/n, and a [L−1] and n are shape factors.
The coefficients 1 and 2 in Eq. [4] represent distinct coarse andfine pore-size distributions, with the parameter w representing
the fine pore-size distribution fraction. The isothermal liquidphase hydraulic conductivity of the medium is then describedusing Mualem’s (1976) approach (Durner, 1994; Priesack andDurner, 2006):l
K = K S éë(1- w ) S e1 + wS e2 ùûìïé
1/ m1 m1 ùïï (1- w )a1 ê1- 1- S e1
ú
ï
ëê
ûú
´ïí
ïï
1
w
w
a
+
a
(
) 1
2
ïï
ïî
m ï
1/ m2 ù 2 ü
wa 2 éê1- 1- S e2
ú ïïï
ë
û
+
ý
(1- w )a1 + wa 2 ïï
ïï
þ

(

(

)

[5]

)

where KS is the saturated hydraulic conductivity [L T−1], and l is apore-connectivity parameter (dimensionless). The thermal liquidphase hydraulic conductivity is calculated by incorporating thetemperature dependence of surface tension in Eq. [5]. This physically based formulation does not contain soil-specific parameters.
Similarly, thermal and isothermal vapor conductivities are calculated by rewriting the thermal and isothermal vapor diffusivity
equations (Saito et al., 2006).

VZJ | Advancing Critical Zone ScienceThe unsaturated hydraulic conductivity formulation given by Eq.
[5] assumes that the soil pore system is made up of a bundle ofcapillaries. Although this assumption may be valid in the wet
and intermediate range of the soil moisture spectrum (0.0 to−0.3 MPa), soil water dynamics in the lower extreme is dominatedby adsorbed film flow, as well as vapor phase flow (Tuller andOr, 2001; Peters et al., 2015; Rudiyanto et al., 2015). Because
water films sustain liquid flow at a higher rate than capillaries
when the water content is very low, the unsaturated hydraulicconductivity at the lower end is likely always somewhat higher
than what is predicted using Eq. [5]. Unfortunately, the contributions of film flow have not been widely adopted yet and are
not included directly in the current version of HYDRUS-1D,
unless they are accounted for indirectly by using dual-porosity
hydraulic properties.

Parameterization of Soil Hydraulic Functions
Figure 2 shows that although the soil packed in the lysimeter has auniform texture, the bulk density varies considerably. Representing
this heterogeneity accurately in the model required partitioning
the model domain into multiple layers, each with a distinct setof hydraulic parameters. This setup would result in an exorbitant
number of free parameters and therefore cause the model to havevery little interpretive value. To limit this problem, we treatedthe model domain as an equivalent homogenous medium with asingle set of reference hydraulic parameters. Variability in the soil
bulk density was accounted for by using a set of scaling parameters (Vogel et al., 1991) included in the HYDRUS-1D software.
This allowed the water content, pressure head, and hydraulicconductivity at each node to be scaled independently by a set of
depth-dependent coefficients. For this approach, it was also important to relate these coefficients to deviations of the bulk densityfrom its reference value. In our analysis we used the depth of theshallowest set of sensors (0.1 m) as the reference depth. The derivations of the scaling relations are detailed in the Appendix. The
final scaling coefficients for the water content, pressure head, andhydraulic conductivity were, respectively:
r -r sq
D XY = o = obqr b -r b[6]ah çæ r b -r s ÷÷öç=÷h o ççè r ob -r b ÷ø

[7]bK çæ r b -r s ÷÷ö
B XY = o = çç o÷çè r b -r b ÷øK

[8]

A XY =

where r b is the soil bulk density (M/L3), r s is the particle density(M/L3), and the superscript o is used to indicate values calculatedusing parameters corresponding to the reference depth (0.1 m).
The coefficients a = −2.82 and b = −6.97 were determined byfitting Eq. [6–8] to water retention and hydraulic conductivity
p. 4 of 17

data of granular porous media packed to different bulk densities(Laliberte et al., 1966).
To further reduce the number of free soil hydraulic parameters, theparameters of the WRC of the reference depth were determinedby fitting the Durner model (Eq. [4]) to the measured water content and pressure head data. For the coarse pore-size fraction, weused data from the in situ water content and pressure head sensors
at 0.1 m. For the fine pore-size fraction, which was assumed to be
independent of packing density, we used the laboratory-measured
water content and pressure head data. The lysimeter and laboratory
data were used to determine the shape coefficients (a and n) of thetwo pore-size fractions and the weight of the finer pore-size fraction(w). The saturated water content (q s) was set to porosity, whereas theresidual water content (q r) was assumed to be zero. These assumptions were made based on laboratory tests with lysimeter soil thatwas exposed to air having different relative humidity values between
0 (oven dried and then stored in a desiccator) and 80%. The corresponding volumetric water contents ranged from 0 to <1%.
In summary, the porosity and WRCs of the entire lysimeter wereinferred from both in situ and laboratory measurements prior
to the HYDRUS-1D simulations. We did not have any directly
measured data that could be used to infer the saturated hydraulic conductivity and pore-connectivity coefficients. These were
therefore determined by inverse modeling using HYDRUS-1D.
The vapor hydraulic conductivity, as well as the thermal liquidphase hydraulic conductivity coefficients, were calculated in
HYDRUS-1D using physically based relations that did not requiresoil-specific parameters. The longitudinal thermal dispersivity wasset at 0.00221 m (Hopmans et al., 2002). The thermal capacitiesof the mineral fraction and water were set to default values of Cm= 2.49 105 J m−3 K−1 and Cw = 5.42 105 J m−3 K−1, respectively(Šimůnek et al., 2016).

Model Setup
In our study, we used HYDRUS-1D (V4.16.0110; Šimůnek
et al., 2016) for all of the simulations. Our objective in using
this model was twofold: to determine the saturated hydraulicconductivity using an inverse simulation, and to gain valuable
insight into the main mechanisms governing surface fluxes andhydraulic redistribution processes within the lysimeter using awell-calibrated model.
The numerical model consisted of a 3-m-deep, single-material soilcolumn that was uniformly discretized into 600 finite element
layers of 0.005-m thickness. The variability in hydraulic properties due to variations in bulk density was mapped onto the flow
domain using the scaling coefficients DXZ , A XZ , and BXZ for thewater content, pressure head, and hydraulic conductivity, respectively. The coefficients were first calculated separately for eachlift of 0.1 m used during packing and subsequently interpolated
linearly to obtain values for each layer. The coefficients of the first

VZJ | Advancing Critical Zone Scienceand last layers were extended to the top and bottom nodes of thenumerical grid.

Initial Conditions
Initial conditions of the soil moisture content were obtained byinterpolating moisture contents measured at the various depths
on 1 Oct. 2011 at 0000 h. Because there were no sensors at thesurface and bottom of the lysimeter, the values obtained from thehighest and lowest sensors were extended to the top and bottomof the model domain, respectively. Initial soil temperature values
were defined similarly, even though we had far more temperature
sensors near the soil surface.

Boundary Conditions
The hydraulic top boundary condition was defined in terms of atime-dependent flux, which was calculated using the measured
hourly changes in the lysimeter mass. This assumes that gains orlosses in the lysimeter mass were caused only by precipitation or
surface evaporation, respectively. The bottom boundary was set asa seepage face, which does not permit outflow when the boundary remains unsaturated. During the timeframe considered in thisstudy, the wetting front had not yet arrived at the bottom of thelysimeter, which justified using these conditions. Lateral flow andlysimeter wall effects were not included in the one-dimensionalmodel. Time-dependent temperature boundary conditions were
imposed at the upper and lower boundaries as derived from thetemperature measurements by TPHP sensors at 2.5 cm and a
HDU sensor at 250 cm, respectively.

Calibration and Testing
The data used in this study started on 1 Oct. 2011 at 0000 h andended a full year later on 30 Sept. 2012 at 0300 h, which resulted
in 8784 hourly measurements. Two subsets of the data were used
to separately calibrate and validate the model. An inverse simulation was conducted using a relatively wet subset of the data fromDay 323 through Day 367 of the simulated time period (hereinreferred to as the water year). Hourly soil water contents from alleight depths were used for the inversion, which generated 4803
data points. Two versions of the inverse model were tested: one
with both the saturated hydraulic conductivity (KS) and the poreconnectivity factor (l) as unknowns, and one with only KS as
unknown. The temperature data were not used for calibration. For
the validation simulation, data from Day 0 to Day 125 were used.
The model’s goodness of fit was expressed in terms of the weighted
RMSE and the coefficient of determination, R2 .

66Results and Discussion
Water Retention Curves

The best-fit calibrated hydraulic conductivity (KS) and poreconnectivity parameter (l) are given in Table 1. Figure 3 shows
measured and fitted WRCs for all depths. Water content andpressure head data pairs from the different in situ quadrants, as
p. 5 of 17

Table 1. Best-fit calibrated parameters.
Description

Value

Fitted saturated hydraulic conductivity (KS), cm h−1

2.00

Pore connectivity parameter (l)

2.0

R2

0.886

RMSE

6.223 ´ 10−5

well as the laboratory samples, are represented by distinct symbols.
Measurements outside the calibration ranges of the in situ pressure
head sensors (−10 4 cm £ h £ −102 cm) and the laboratory WP4
potentiometer data (h £ −5 ´ 10 4 cm) were not included. The
apparent hysteresis observed at 10 cm was thought to be due to a
response time lag of the pressure head sensors. Additionally, thepoor sensitivity of the in situ sensors near the upper limit of thecalibration range did lead to an apparent asymptotic rise in measured water contents at pressure heads near −100 cm. However,
because the exact thresholds of these biases were unknown, theDurner dual-porosity model was fitted to all measured data within
the calibration range at a depth of 10 cm.
A nonlinear optimization code written in R (R Core Team, 2016)
was used to fit the Durner model to the water retention data of
the reference layer. Table 2 shows the best-fit Durner parametersof the lysimeter soil at a depth of 10 cm. Water retention parameters of the remaining layers were determined using Eq. [6–8] as
described in the Appendix. Table 3 lists the corresponding scaling
parameters. The dashed lines in Fig. 3 represent the fitted retention
curves, and the solid lines represent the depth-adjusted WRCs.
The differences in water retention between the layers can be attributed mostly to differences in bulk density (and therefore porosity),
as shown in Fig. 2. Saturated water contents were set equal to themeasured depth-specific porosities (i.e., q s = f). We visually corroborated the scaled WRCs at all depths, but the goodness of fit
was most obvious at depths of 25 to 75 cm.

Lysimeter Mass
Figure 4 shows the mass of Lysimeter 1 and the near-surface(6–12 mm) soil temperature data recorded during water year 2011.
The temperature data clearly exhibited an annual cycle, with therange of the data points reflecting diurnal cycles. The lysimeter mass
exhibited a continuous decline, which indicated evaporative losses,
although there were occasional abrupt spikes during precipitation
events. A change in lysimeter mass of 1 kg corresponded to a change
in moisture depth of 0.25 mm. The observed range of changes in thewater depth during the period of interest was ~63 mm.
Figure 5 shows the daily mass, vapor density, humidity, andtemperature dynamics at four selected soil moisture periods. Whenthe time series data were differentiated numerically to calculate
lysimeter mass fluxes, small fluctuations (below the calibrated
accuracy of the scale) were amplified. These fluctuations in the

VZJ | Advancing Critical Zone Sciencedata were filtered using smoothed spline functions (R Core Team,
2016). We compared the original lysimeter mass data with the mass
reconstructed by integration of the calculated flux to ensure that
this smoothing did not alter the integrity of the data. Each periodof the plots shown in Fig. 5 lasted 5 d. Figure 5 also shows the airtemperature, relative humidity, and vapor density (measured 2 mabove ground) values, as well as the near-surface soil temperatureand lysimeter mass data. Two of the events represented relatively
wet near-surface conditions (starting within 2 d after a rain event),
whereas the other two events covered relatively dry conditions longafter a rain event. The last two events were only 4 wk apart andseparated by a period that had the largest rain event of the water
year. The overall drying rate was more intense when the near-surface
soil was still moist (1.2–1.5 kg d−1) compared with the dry periods(~0.5 kg d−1, approximately equal to the load-cell resolution).
A striking feature of Fig. 5 is the considerable diurnal mass fluctuation shown in the lower panel. One possible explanation forthis periodicity is a thermal artifact of the weighing mechanismof the lysimeter. The diurnal temperature fluctuations of thenear-surface soil (also plotted) were nearly identical at all times. A
careful examination of the four intervals showed that the nighttime mass gains were more definitive (larger than the instrument
resolution and less noisy) only when the surface soil was still wet.
During these wet times, approximately half of the daily losses
were regained during the night and early morning. In contrast,
the mass gains (or fluctuations) during the dry times were very
noisy. Additionally, the gain in mass started at or before midnightfor wet soils when the soil surface was still cooling. This suggeststhat thermal artifacts on the performance of the scale do not fully
explain the mass gain.
An alternative and more consequential explanation for the fluctuations is vapor condensation (dew) or direct vapor adsorption by thesoil caused by large daily fluctuations in the surface-soil temperature(~30°C at 6–10 cm below the surface). The observed nighttime
mass gains during the wet periods (~0.5 mm) are consistent withobservations of nightly gains of soil water in the Negev Desert, which
has a climate comparable with that of the Las Vegas area (Agam andBerliner, 2004; Hill et al., 2015). The above considerations justify
including the mass gain as a surficial deposition of moisture.

Model Calibration
The subset of data from water year 2011–2012 used for model
calibration is indicated by the shaded region in Fig. 4, marked as“Cal.” This subset commenced when the lysimeter was the lightest (driest) during the year and contains data from three major
storms, including the largest storm of the year (which added >50
mm in precipitation). Figure 6 compares the measured and bestfit modeled cumulative mass fluxes across the lysimeter surface.
Because mass losses through all surfaces of the lysimeter otherthan the soil surface were considered negligible, the cumulative
mass flux should be equivalent to the net mass change from thep. 6 of 17

Fig. 3. Water retention curves of the lysimeter soil at various depths. The symbols correspond to water contents and pressure heads measured at thecorresponding depths and quadrants (NE, northeast; NW, northwest; SW, southwest; SE, southeast) or laboratory data derived from dew-point potentiometer (WP4) measurements. The dashed lines represent the Durner’s dual-porosity model fitted to the 10-cm data, and the solid lines are curvestransformed using linear scaling coefficients.

VZJ | Advancing Critical Zone Sciencep. 7 of 17

Table 2. Directly determined hydraulic parameters of thelysimeter soil layer at the 10-cm depth.
Parameter

Value

Saturated water content (q s)

0.292

Residual water content (q r)0.000

Durner parameter a1 (coarse subcurve), cm−1

1.57 ´ 10−2

Durner parameter n1 (coarse subcurve)

2.571

Durner parameter a 2

10−4(fine subcurve), cm−1

Durner parameter n2 (fine subcurve)

1.35

Fraction of fine subcurve (w)

0.245

start of the calibration period. The model generally captured theinfiltration and evaporation dynamics well. However, the simulated
rate of drying was slower than the observations. This discrepancy
increased over time, which led to the widening gap shown in Fig. 6.
Dixon (1999) found similar results comparing measured with
HYDRUS-1D-simulated soil water storage in the bare soil of aweighing lysimeter at the Nevada Test Site in Nye County, Nevada.
Figure 7 shows a comparison of measured and simulated watercontents of all eight observation layers. Only the water content
data were used in the calibration. The reported error bands are
standard measurement errors from the four lysimeter quadrants,
except for the data at depths of 100, 200, and 250 cm, which each
had only two sensors. In general, the model accurately captured theobserved water content dynamics in the top 1 m of the lysimeter,
except for two distinct discrepancies. One discrepancy occurred
during the late stages of drying at 10 cm, where the simulated rateof drying lagged behind the measured rate. This lag is consistent
with the discrepancy in the surface mass flux noted earlier. The
other discrepancy occurred when the wetting front reached 75 cm.
The standard error of the water content at 75 cm (according to the
TDR sensors at the four quadrants) showed an abrupt increaseat the same time. A closer inspection of the data from individual

Fig. 4. Temperatures of the near-surface lysimeter soil and validated(Val.) and calibrated (Cal.) total mass of the lysimeter during water
year 2011–2012.

quadrants (data not shown) revealed that the northwest quadrant
was slightly drier than the other locations at and below 75 cm. Thissuggests that a preferential flow path may have developed at that
depth. At 100 cm, which had sensors only in the northeast andsouthwest quadrants, the arrival of the wetting front was muchslower. Incomplete contact between the sensor and the soil may
have been another reason for the differences in moisture contents
measured at 75-cm depth. Although the sensors were carefully
installed in the disturbed soil, incomplete sensor-to-soil contact
cannot be ruled out. Below 1 m, no changes in the moisture content were observed in the water year 2011–2012.
Figure 8 shows a comparison of measured and simulated temperatures of the eight observation layers. The temperature datasets
were not used to calibrate the model, which confirms that theTable 3. Fitted scaling parameters for the remaining depths.
Fitted scaling parameter†
Porosity(r b ¢ − r s)/(r b − r s)

DXY

A XY

BXY

1.70

0.292

1.000

1.000

1.000

1.000

1.74

0.275

0.941

0.941

1.165

0.693

50

1.84

0.232

0.793

0.793

1.785

0.249

75

1.64

0.303

1.037

1.037

0.913

1.244

100

1.91

0.215

0.735

0.735

2.159

0.158

150

1.76

0.268

0.916

0.916

1.245

0.591

200

1.74

0.253

0.864

0.864

1.440

0.417

250

1.74

0.251

0.857

0.857

1.471

0.396

Depth

Bulk densitycm

g cm−31025† r b ¢, bulk density; r s, particle density; DXY, linear scaling coefficient of water content; A XY, linear scaling coefficient of pressure head; BXY, linear scaling coefficient ofhydraulic conductivity.

VZJ | Advancing Critical Zone Sciencep. 8 of 17

Fig. 5. Diurnal dynamics of near-surface soil temperature, air temperature, air relative humidity, air vapor density, and lysimeter mass data during four
periods when the surface soil was wet and dry.

VZJ | Advancing Critical Zone Sciencep. 9 of 17thermal properties of the soil were set a priori. The model captured the temperature dynamics in the top 25 cm accurately, but
the results deviated from the measurements starting at 50 cm. The
deviations were probably caused by variations in soil bulk density (and therefore thermal properties) with depth, with the mostnotable variability occurring in the 40- to 120-cm region (Fig. 2).
However, the solid fraction in the model was kept uniform (setto the value corresponding to the porosity at 10-cm depth), since
HYDRUS-1D does not permit any scaling of thermal properties.
The discrepancy thus appears to be caused by inconsistent thermal
properties. The discrepancies at 50 and 100 cm are propagated by
the model to the underlying layers.

Discrepancy and Diurnal Fluctuations
Fig. 6. Measured cumulative mass fluxes, assuming that lysimeter
mass changes were due to water fluxes through the upper boundary only, and corresponding best-fit modeled mass fluxes during thecalibration period.

To better understand the main causes of the discrepancy between
the simulated and observed drying rates, we examined thesimulated surface pressure heads, as well as the simulated andmeasured surface mass fluxes at two 5-d intervals during thecalibration period (Fig. 9). The first interval started immediately before the second rain event (Day 333) andincluded the wettest period of water year 2011–2012. The second interval occurred during a dry
spell at the end of the calibration period. During
the first 3 d of the wet period (Fig. 9a), the surface
soil remained fairly wet, and the simulated andmeasured surface mass fluxes were essentially identical. However, when the soil surface began to dry
significantly, reaching a pressure head minimumof −10 6 cm (air-dry conditions), the model was
not able to reproduce the measured flux (see theafternoon of Day 338). During the dry period (Fig.
9b), the surface soil remained fairly dry (especiallyin the afternoons), and the model consistently
failed to match the observed upward fluxes. The
accumulation of these daily discrepancies resulted
in deviations between the observed and modeled
cumulative mass fluxes as shown in Fig. 6.

Fig. 7. Measured and modeled water contents at eight layers during the calibration period.
The shaded regions represent standard errors (SE) of four replicate sensors except at 100,
200, and 250 cm, where they denote the range of two sensors.

VZJ | Advancing Critical Zone ScienceThe above discrepancies suggest that the nonisothermal f low processes as described by theextended Richards equation (Eq. [1]) do not fully
represent the high drying rates when the surface
soil was at or near air-dry conditions. To better
understand the flow dynamics during times whenthe model failed to capture the observed drying
rates, we examined the vapor and liquid fluxes inthe top 100 cm on Days 336 and 357, respectively.
These dates are in the middle of the wet and dry
periods. The dashed lines in the top plot of Fig.
10 represent the points in time at which the fluxes
were calculated. Overall, vapor flow played a minor
role compared with the liquid flux. A significant
vapor flux occurred only during the second half
p. 10 of 17

of the wet day and was limited to the
top few centimeters. During the relatively dry Day 357, a downward flux was
observed during the first half of the day,with the total flux being divided almost
equally between the vapor and liquid
f luxes. However, in the afternoon of
Day 357, just when the model started to
fail in capturing the observed flux, theupward flux was driven only by liquidwater flow that mostly originated fromthe top few centimeters. No upward
flux occurred below 5 cm. The model
thus failed to fully capture the drying
process, since there was no sufficientliquid f lux. Additionally, the surface
soil reached a pressure head of −10−6 cm
at this time (Fig. 9, top). These results
are consistent with the basic premisesof the van Genuchten–Mualem model(Mualem, 1976; van Genuchten, 1980)of unsaturated hydraulic conductivity.
When the soil is in such a dry state, avery strong gradient in the pressure headcannot overcome the steep drop in thehydraulic conductivity.

Validation Period
Figure 11 shows a comparison of measured and best-fit simulated cumulative
mass f luxes at the lysimeter surface.
The validation period was marked byfour low-intensity rain events, which
added a total of ~5 m of precipitation
to the lysimeter. However, the lysimeter
Fig. 8. Measured and simulated temperatures at eight layers during the calibration period. The shaded
regions represent the standard errors of four replicate sensors except at 200 and 250 cm, where theyrecorded a net loss of ~0.6 cm throughdenote the range of two sensors.
the validation period. The wetting front
from the rain events clearly registered atModel Limitations
10 cm and showed very small and slow increases at 25 cm (Fig.
andthe Way Forward12), which indicates that most or all of the rainwater had been
lost by evaporation and not deep drainage. The model accurately
Evaporation from initially wet, bare porous media is typically
predicted the wetting dynamics during all four rain events but
divided into three distinct stages (Idso et al., 1974; Or et al.,failed to fully capture the drying dynamics. This discrepancy2013). The first two stages—commonly referred to as Stage I, or
between the model predictions and the measurements increased
the constant rate period, and Stage II, or the transition period—areduring the slow but steady drying process that followed the raincontrolled by the evaporative demand of the ambient atmospherefall events.
and last as long as there is a continuous liquid water supply tothe surface facilitated by upward capillary flow (Or et al., 2013).
In comparison, HYDRUS-1D predicted the temperature dynamThe transition from Stage I to Stage II evaporation represents theics at all depths reasonably well, as shown by the results in Fig.
transition from saturated or near-saturated flow to increasingly
13. The predicted temperature dynamics exhibited wider diurnal
unsaturated upward flow at rates that are less that the potential
fluctuations than the sensor data at 10 and 25 cm. However, there
evaporation rate. The phase transition of water from liquid to
was no lateral variability in the recorded temperatures, which isvapor during Stages I and II mostly occurs at the soil surface, andevident by the nearly invisible error margins.
this process controls the actual evaporation rate. In this regard,

VZJ | Advancing Critical Zone Sciencep. 11 of 17

HYDRUS-1D and the underlying Richards equation
satisfactorily capture subsurface liquid water flow. Thisis evident in the comparisons between the simulated andmeasured mass losses over the few days immediately after
each rain event, during both the calibration and validation
periods (Fig. 6 and 11).

Fig. 9. Modeled surface pressure head and prescribed and actual modeled surface
fluxes of water during (a) wet and (b) dry periods. HYDRUS-1D could not model
the prescribed flux when the surface soil was very dry.

During Stage II, liquid water flow to the surface occurs
under increasingly unsaturated conditions to meet the evaporative demand at the surface. When the rate of water supply
is no longer adequate to match the evaporative demand, evenwith increasingly steeper pressure head gradients near thesoil surface, the continuous stream of liquid water that supplies the evaporating surface breaks and a dry soil layer forms
at the surface. This results in a further drop in the evaporation rate and marks the transition from Stage II to StageIII evaporation. Stage III evaporation is controlled by water
vapor diffusion through the dry layer above the secondary
drying front (Shokri et al., 2008), where phase changes occur
at and below the soil surface. As the secondary drying front
recedes deeper into the soil, the length of the diffusive pathincreases, which leads to a continuous decrease in the overall
evaporation rate. The drop in the evaporation rate during

Fig. 10. Prescribed and actual simulated surface fluxes of water during wet (left) and dry (right) days and liquid- and vapor-phase fluxes in the top 100
cm at selected times. The dashed lines represent the points in time at which the fluxes were calculated.

VZJ | Advancing Critical Zone Science

p. 12 of 17

the transition between Stages II and III depends on the depth ofthe drying front when the liquid-phase discontinuity first emerges
and a secondary drying front forms (Shokri et al., 2008). Lehmann
et al. (2008) suggested that this depth is a characteristic property ofthe porous medium that depends on the pore-size distribution of
the medium. Shokri and Or (2011) estimated values between 3 and
14 mm (depending on the soil texture, ranging from coarse to veryfine) as the initial thickness of the dry layer and the initial depth of
the secondary drying front.
The results of our study highlighted two inherent limitations of
the governing equations (Eq. [1–5]) underlying HYDRUS-1Dand similar simulators. One limitation is that although the modelaccounts for nonisothermal vapor diffusion (see Eq. [1]), it does
not permit discontinuity in the liquid phase. Therefore, the vapor
pressure is assumed to be in equilibrium with a continuous but
increasingly drier matrix during Stage III. This assumption limits
the vapor concentration gradient that can beestablished between the drying front and thesurface. Model limitations associated withthis shortcoming can only be alleviated byenabling liquid-phase discontinuities with amoving boundary within the model domain,such as using a secondary drying front as
proposed by Shokri et al. (2008). This forced
continuity leads to a pressure head that asymptotically tends toward infinity near a drying
soil surface. To avoid numerical instability
and thermodynamically unrealistic pressure
head values, HYDRUS-1D imposes a minimum pressure head. In this study, we used avalue of hcrit = −106 cm, which corresponds
to typical air-dry conditions (the model results
were not sensitive to variations in the range−10 8 cm < hcrit < −10 6 cm). The bimodal
Durner WRC afforded adequate flexibility
that matched the measured behavior in thehyper-dry and wet regimes.
The other limitation is that the underlying
assumption behind Mualem’s model of theunsaturated hydraulic conductivity does notadequately describe flow in the hyper-dry
regime. The model was derived assuming thatthe porous medium can be approximated bya bundle of capillary tubes. This conceptualization adequately describes the retention
and dynamics of capillary menisci in wet andmoderately wet systems but does not represent the soil water dynamics well whensoil water is distributed in the form of thinfilms. The discrepancy is particularly noticeable when simulating liquid f low. The

VZJ | Advancing Critical Zone ScienceFig. 11. Measured cumulative mass fluxes, assuming that lysimeter masschanges were caused by water fluxes across the upper boundary only, andcorresponding predicted mass fluxes during the validation period.

Fig. 12. Measured and predicted water contents at eight layers during the validation period. The
shaded regions represent standard errors of four replicate sensors except at 200 and 250 cm, where
they denote the range of two sensors.

p. 13 of 17

soil moisture dynamics and evaporation invery dry soils. A model that can capture allof the infiltration, redistribution, and evaporation processes of water in arid soils will bea valuable tool for quantifying the moisturedynamics of arid ecosystems and their effects
on climate change, plant growth, erosion, andrecharge patterns.

66Conclusions

The SEPHAS lysimeter facility provides aunique dataset for desert hydrologic research.
The lysimeters are instrumented with reliable
and accurate sensors with robust quality control and curating mechanisms in place. The
well-defined upper boundary condition andthe accurate scale allowed us to close the waterbalance. Additionally, the high density of dataacquisition allowed for a reliable representation of water flow in unsaturated desert soil.

In this study, we tested the reliability of the conventional approaches of modeling unsaturated
flow in the dry end of the soil water spectrum.
We showed that infiltration and redistribution could be represented fairly accurately
with the Richards equation. The model also
performed well in describing the early stagesof evaporation (i.e., Stage I and parts of Stage II
evaporation), when the evaporative demand at
the soil surface is the primary driver. However,
the model consistently underestimated Stage
Fig. 13. Measured and modeled temperatures at eight layers during the validation period. The
III evaporation. We identified two potential
shaded regions represent the standard error of four replicate sensors except at 200 and 250 cm,
where they denote the range of two sensors.
causes that contributed to this discrepancy.
First, the classical Richards equation does not
permit hydraulic discontinuity and vapor-difbundle-of-capillaries assumption requires the unsaturated hydraulicfusion-limited water transfer from a subsurface drying (evaporation)
conductivity to drop indefinitely at a steep rate with lower negative
front upward to the soil surface. The forced continuity requirespressure heads. However, conceptual models of unsaturated flow
unrealistic equilibration between an extremely dry but continuousthat integrate both capillary and film flow (Tuller and Or, 2001)
thread of liquid water and vapor. Vapor flow that occurs under suchshow that the decline in hydraulic conductivity with decreasing presconditions is likely to underestimate Stage III evaporation. Second,
sure heads is more tapered. This problem can be addressed by using
the hydraulic conductivity model used in the HYDRUS-1D code
unsaturated hydraulic conductivity functions that integrate capillary
likely underestimates water flow rates for very dry soil conditions.
and film flow and possibly vapor flow (Tuller and Or, 2001; Peters etMore work is clearly needed to delineate the relative contributions ofal., 2015; Rudiyanto et al., 2015).
these two shortcomings and to correct for their effects.
The discrepancy between simulated and measured evaporation rates
in this study is likely related to both of the above limitations (i.e.,
properly accounting for evaporation from moving boundaries awayfrom the soil surface and accounting for film flow). The high-quality,
long-term data available from the SEPHAS lysimeter facility provide
a compelling opportunity to address this critical gap in modeling

VZJ | Advancing Critical Zone Science

66Appendix: Linear Scalingof Hydraulic Properties

The variability in bulk density in the lysimeter was accounted forin the simulations by using scaling coefficients. Water contents,
p. 14 of 17

pressure heads, and hydraulic conductivities at each node canbe independently scaled by coefficients in HYDRUS-1D. The
approach permits one to consider the profile as a homogeneous
medium using parameters of the reference state, but with adjustments for each separate soil layer using the three coefficients.
Consider the WRC of a soil at a reference bulk density (vanGenuchten, 1980):
n -m

q = q r + (q s -q r ) éê1 + (ah) ùú
ë
û

[A1]

where q s = f max. Similarly, the WRC of the same soil but at different bulk density can be written as:

-m
n¢ ù
é
q ¢ = q ¢r + (q ¢s -q ¢r ) ê1 + (a ¢h ¢) ú
ë
û

¢

[A2]

Equation [A1] represents the reference state and Eq. [A2] thealtered state. Following Vogel et al. (1991), the dimensionless scaling rules that are permitted in HYDRUS-1D areq ¢ -q ¢r= D XZ
q-q r

[A3]

h¢
= A XZh[A4]Fig. A1. Water retention curves of (a,d) Touchet silt loam, (b,e) Columbia sandy loam, and (c,f ) unconsolidated sand packed to varying degrees ofbulk density. Original data derived from Laliberte et al. (1966) are in the top row. Water contents and pressure heads were rescaled by the bulk density
according to the proposed formulation (bottom row). The solid lines represent van Genuchten water retention fits to the lowest bulk density sampleof each soil.

VZJ | Advancing Critical Zone Sciencep. 15 of 17

Haploxerolls), Columbia sandy loam (coarse-loamy, mixed, superactive, nonacid, thermic Oxyaquic Xerofluvents), and unconsolidated
sand (Laliberte et al., 1966). The best-fit value of a was found to be−2.83. Similarly, Fig. A2 shows scaled KS values vs. porosity data
for all three soils. The data exhibit consistent scaling with a best-fit
value of 6.97 for b. The two scaling coefficients a and b consistently
worked well for the three granular media and were assumed to be
applicable to the SEPHAS lysimeter as well.

Acknowledgments

This paper is based on work supported by the US National Science Foundation under Grants
no. IIA-1301726 and EPS-0447416, as well as grants from the Desert Research Institute’s Division of Hydrologic Sciences International Student Program and the Royal Dutch Institutefor Engineers (KIVI NIRIA). We are grateful to Ms. Nicole Damon for carefully editing themanuscript, as well as Dr. Markus Flury and three anonymous reviewers for their valuable comments and suggestions.

References

Agam, N., and P.R. Berliner. 2004. Diurnal water content changes
in the bare soil of a coastal desert. J. Hydrometeorol. 5:922–933.
doi:10.1175/1525-7541(2004)005<0922:DWCCIT>2.0.CO;2Fig. A2. Saturated hydraulic conductivity (KS) of (a) Touchet silt
loam, (b) Columbia sandy loam, and (c) unconsolidated sand packed
to varying degrees of bulk density (rb) derived from Laliberte et al.
(1966). The solid line represents the proposed scaling relation. The
superscript o is used to indicate values calculated using parameters
corresponding to the reference depth (0.1 m).

K S¢= B XZKS[A5]

Substituting Eq. [A3] and [A4] in Eq. [A1] and comparing theresult with Eq. [A2] leads tor ¢ -r ¢s
D XY = br b -r s[A6]

This shows that we can then use a power-law relation to describe
the other two scale factors:
æ r b¢ -r s ö÷a÷A XY = ççççè r -r ÷÷ø[A7]

æ r ¢ -r s ö÷b÷÷
B XY = ççç bèç r b -r s ø÷

[A8]

b

s

Therefore, what remains to be found are the appropriate values ofa and b. A similar approach was used by Assouline et al. (1997).
Figure A1 shows how Eq. [A6] and [A7] can be used to rescale theWRC of compacted soils (top row) to match the uncompacted WRC(bottom row). The data correspond to variably packed columns ofTouchet silt loam (coarse-silty, mixed, superactive, mesic Cumulic

VZJ | Advancing Critical Zone Science

Andraski, B.J. 1997. Soil-water movement under natural-site and wastesite conditions: A multiple-year field study in the Mojave Desert, Nevada. Water Resour. Res. 33:1901–1916. doi:10.1029/97WR01502
Andraski, B.J., W.A. Jackson, T.L. Welborn, J.K. Bohlke, R. Sevanthi, and
D.A. Stonestrom. 2014. Soil, plant, and terrain effects on natural perchlorate distribution in a desert landscape. J. Environ. Qual. 43:980–
994. doi:10.2134/jeq2013.11.0453
Assouline, S., J. Tavares, and D. Tessier. 1997. Effect of compaction on soil physical and hydraulic properties: Experimental results and modeling. Soil Sci. Soc. Am. J. 61:390–398.
doi:10.2136/sssaj1997.03615995006100020005x
Austin, A.T., L. Yahdjian, J.M. Stark, J. Belnap, A. Porporato, U. Norton, et al.
2004. Water pulses and biogeochemical cycles in arid and semiarid
ecosystems. Oecologia 141:221–235. doi:10.1007/s00442-004-1519-1
Bhattachan, A., P. D’Odorico, M.C. Baddock, T.M. Zobeck, G.S. Okin,
and N. Cassar. 2012. The Southern Kalahari: A potential new dust
source in the Southern Hemisphere? Environ. Res. Lett. 7:024001.
doi:10.1088/1748-9326/7/2/024001
Chief, K., M.H. Young, B.F. Lyles, J. Healey, J. Koonce, E. Knight et al. 2009.
Scaling environmental processes in heterogeneous arid soils: Construction of large weighing lysimeter facility. Desert Res. Inst., Las Vegas, NV.
Chung, S.O., and R. Horton. 1987. Soil heat and water-flow with a partial
surface mulch. Water Resour. Res. 23:2175–2186.
Constantz, J., S.W. Tyler, and E. Kwicklis. 2003. Temperature-profile methods for estimating percolation rates in arid environments. Vadose
Zone J. 2:12–24. doi:10.2136/vzj2003.1200
D’Odorico, P., A. Bhattachan, K.F. Davis, S. Ravi, and C.W. Runyan. 2013.
Global desertification: Drivers and feedbacks. Adv. Water Resour.
51:326–344. doi:10.1016/j.advwatres.2012.01.013
Devitt, D.A., L.F. Fenstermaker, M.H. Young, B. Conrad, M. Baghzouz, andB.M. Bird. 2011. Evapotranspiration of mixed shrub communities inphreatophytic zones of the Great Basin region of Nevada (USA). Ecohydrology 4:807–822. doi:10.1002/eco.169
Dixon, J. 1999. An evaluation of unsaturated flow models in an arid climate. M.S. thesis. Univ. of Nevada, Las Vegas. doi:10.2172/759261
Dregne, H., M. Kassas, and B. Rozanov. 1991. A new assessment of theworld status of desertification. Desertification Control Bull. 20:6–18.
Durner, W. 1994. Hydraulic conductivity estimation for soil with heterogeneous pore structure. Water Resour. Res. 30:211–223.
doi:10.1029/93WR02676
Gee, G.W., J.M. Keller, and A.L. Ward. 2005. Measurement and prediction of deep drainage from bare sediments at a semiarid site. Vadose
Zone J. 4:32–40. doi:10.2113/4.1.32
Gee, G.W., P.J. Wierenga, B.J. Andraski, M.H. Young, M.J. Fayer, and
M.L. Rockhold. 1994. Variations in water-balance and rechargepotential at 3 western desert sites. Soil Sci. Soc. Am. J. 58:63–72.
doi:10.2136/sssaj1994.03615995005800010009x

p. 16 of 17

Hare, F.K. 1985. Climate variations, drought and desertification. World
Meteorol. Organ., Geneva, Switzerland.
Hill, A.J., T.E. Dawson, O. Shelef, and S. Rachmilevitch. 2015. The
role of dew in Negev Desert plants. Oecologia 178:317–327.
doi:10.1007/s00442-015-3287-5
Hopmans, J.W., J. Šimůnek, and K.L. Bristow. 2002. Indirect estimation ofsoil thermal properties and water flux using heat pulse probe measurements: Geometry and dispersion effects. Water Resour. Res. 38(1).
doi:10.1029/2000WR000071
Huxman, T.E., K.A. Snyder, D. Tissue, A.J. Leffler, K. Ogle, W.T. Pockman etal. 2004. Precipitation pulses and carbon fluxes in semiarid and arid
ecosystems. Oecologia 141:254–268. doi:10.1007/s00442-004-1682-4
Idso, S.B., R.J. Reginato, R.D. Jackson, B.A. Kimball, and F.S. Nakayama.
1974. The three stages of drying of a field soil. Soil Sci. Soc. Am. J.
38:831–837. doi:10.2136/sssaj1974.03615995003800050037x
Kenny, J.F., and K.E. Juracek. 2012. Description of 2005–10 domestic water use for selecting U.S. cities and guidance for estimating domestic
water use. USGS, Washington, DC.
Kidron, G.J., and V.P. Gutschick. 2013. Soil moisture correlates withshrub-grass association in the Chihuahuan Desert. Catena 107:71–79.
doi:10.1016/j.catena.2013.02.001
Koonce, J.E. 2016. Water balance and moisture dynamics of an arid andsemi-arid soil: A weighing lysimeter and field study. Ph.D. diss., Univ. of
Nevada, Las Vegas.
Laliberte, G.E., A.T. Corey, and R.H. Brooks. 1966. Properties of unsaturated porous media. Hydrol. Pap. 17. Colorado State Univ., Fort Collins.
Lehmann, P., S. Assouline, and D. Or. 2008. Characteristic lengths affecting evaporative drying of porous media. Phys. Rev. E 77:056309.
doi:10.1103/PhysRevE.77.056309
Moreno-de las Heras, M., L. Turnbull, and J. Wainwright. 2016. Seed-bankstructure and plant-recruitment conditions regulate the dynamics ofa grassland-shrubland Chihuahuan ecotone. Ecology 97:2303–2318.
doi:10.1002/ecy.1446
Morgan, J.A., D.E. Pataki, C. Korner, H. Clark, S.J. Del Grosso, J.M. Grunzweig et al. 2004. Water relations in grassland and desert ecosystems exposed to elevated atmospheric CO2. Oecologia 140:11–25.
doi:10.1007/s00442-004-1550-2

Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous materials under temperature gradients. Am. Geophys. Union Trans. 38:222–
232. doi:10.1029/TR038i002p00222
Priesack, E., and W. Durner. 2006. Closed-form expression for the multimodal unsaturated conductivity function. Vadose Zone J. 5:121–124.
doi:10.2136/vzj2005.0066
R Core Team. 2016. R: A language and environment for statistical computing. R Found. Stat. Comput., Vienna.
Rudiyanto, M. Sakai, M.Th. van Genuchten, A.A. Alazba, B.I. Setiawan,
and B. Minasny. 2015. A complete soil hydraulic model accounting for capillary and adsorptive water retention, capillary andfilm conductivity, and hysteresis. Water Resour. Res. 51:8757–8772.
doi:10.1002/2015WR017703
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. doi:10.2136/vzj2006.0007
Scanlon, B.R., D.G. Levitt, R.C. Reedy, K.E. Keese, and M.J. Sully. 2005.
Ecological controls on water-cycle response to climate variability indeserts. Proc. Natl. Acad. Sci. 102:6033–6038.
Shokri, N., P. Lehmann, P. Vontobel, and D. Or. 2008. Drying front andwater content dynamics during evaporation from sand delineated by neutron radiography. Water Resour. Res. 44:W06418.
doi:10.1029/2007WR006385
Shokri, N., and D. Or. 2011. What determines drying rates at the onset ofdiffusion controlled stage-2 evaporation from porous media? Water
Resour. Res. 47:W09513. doi:10.1029/2010WR010284
Šimůnek, J., M.T. van Genuchten, and M. Sejna. 2016. Recent developments and applications of the HYDRUS computer software packages.
Vadose Zone J. 15(7). doi:10.2136/vzj2016.04.0033
Soil Survey Division Staff. 1993. Soil survey manual. Agric. Handb. 18. USDA
Soil Conserv. Serv., Washington, DC.
Tuller, M., and D. Or. 2001. Hydraulic conductivity of variably saturated
porous media: Film and corner flow in angular pore space. Water Resour. Res. 37:1257–1276. doi:10.1029/2000WR900328
van Genuchten, M.T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–
898. doi:10.2136/sssaj1980.03615995004400050002x

Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
doi:10.1029/WR012i003p00513

Vogel, T., M. Cislerova, and J.W. Hopmans. 1991. Porous media with linearly variable hydraulic properties. Water Resour. Res. 27:2735–2741.
doi:10.1029/91WR01676

Or, D., P. Lehmann, E. Shahraeeni, and N. Shokri. 2013. Advances in soil evaporation physics: A review. Vadose Zone J. 12(4).
doi:10.2136/vzj2012.0163

Western Regional Climate Center. 2017. Cooperative climatological
data summaries. Western Regional Clim. Ctr., Reno, NV. http://www.
wrcc.dri.edu/cgi-bin/cliMAIN.pl?nv1071 (accessed 7 Apr. 2017

Peters, A., S.C. Iden, and W. Durner. 2015. Revisiting the simplified evaporation method: Identification of hydraulic functions considering vapor, film
and corner flow. J. Hydrol. 527:531–542. doi:10.1016/j.jhydrol.2015.05.020

VZJ | Advancing Critical Zone Science

p. 17 of 17