HYDROLOGICAL PROCESSES
Hydrol. Process. 24, 1095– 1105 (2010)
Published online 8 January 2010 in Wiley InterScience(www.interscience.wiley.com) DOI: 10.1002/hyp.7571

Modelling investigation of water partitioning at a semiarid
ponderosa pine hillslope
Huade Guan,1 * Jirka Simunek,2 Brent D. Newman3† and John L. Wilson1
1 Department of Earth and Environmental Science, New Mexico Tech, Leroy Place 801, Socorro, NM 87801, USA
2 Department of Environmental Sciences, University of California, Riverside, CA 92521, USA
3 Earth and Environmental Sciences Division, MS J495, Los Alamos National Laboratory, Los Alamos, NM 87545 USA

Abstract:
The effects of vegetation root distribution on near-surface water partitioning can be two-fold. On the one hand, the roots
facilitate deep percolation by root-induced macropore flow; on the other hand, they reduce the potential for deep percolation by
root-water-uptake processes. Whether the roots impede or facilitate deep percolation depends on various conditions, including
climate, soil, and vegetation characteristics. This paper examines the effects of root distribution on deep percolation into
the underlying permeable bedrock for a given soil profile and climate condition using HYDRUS modelling. The simulations
were based on previously field experiments on a semiarid ponderosa pine (Pinus ponderosa) hillslope. An equivalent single
continuum model for simulating root macropore flow on hillslopes is presented, with root macropore hydraulic parameterization
estimated based on observed root distribution. The sensitivity analysis results indicate that the root macropore effect dominates
saturated soil water flow in low conductivity soils (Kmatrix below 107 m/s), while it is insignificant in soils with a Kmatrix larger
than 105 m/s, consistent with observations in this and other studies. At the ponderosa pine site, the model with simple rootmacropore parameterization reasonably well reproduces soil moisture distribution and some major runoff events. The results
indicate that the clay-rich soil layer without root-induced macropores acts as an impeding layer for potential groundwater
recharge. This impeding layer results in a bedrock percolation of less than 1% of the annual precipitation. Without this
impeding layer, percolation into the underlying permeable bedrock could be as much as 20% of the annual precipitation. Thissuggests that at a surface with low-permeability soil overlying permeable bedrock, the root penetration depth in the soil iscritical condition for whether or not significant percolation occurs. Copyright  2010 John Wiley & Sons, Ltd.
KEY WORDS

hillslope; percolation; root macropore modelling; recharge; semiarid; ponderosa pine; New Mexico; root
distribution; macropore flow

Received 24 February 2009; Accepted 4 November 2009

INTRODUCTION
Groundwater recharge is an important factor for longterm water resource management, especially in arid
and semiarid regions. Due to the dry climate, thickvadose zones, and vegetation water uptake, groundwater recharge in arid and semiarid basin floors is typicallyvery small (Phillips, 1994; Izbicki et al., 2000; Walvoord
et al., 2002; Flint et al., 2004; Walvoord and Phillips,
2004; Scanlon et al., 2005). Although vegetation functions in partitioning shallow subsurface water in two
ways: (1) through root-water-uptake reducing or completely eliminating deep percolation (potential recharge),
or even causing vertical groundwater discharge where
the water table is shallow, and (2) through root-induced
macropores facilitating deep percolation, the first mechanism seems to be dominant in arid and semiarid regions(Gee et al., 1994; Stothoff et al., 1999). With increasing* Correspondence to: Huade Guan, School of Chemistry, Physics andEarth Sciences, National Centre for Groundwater Research and Training,
Flinders University, GPO Box 2100, SA 5001, Australia. E-mail:huade.guan@flinders.edu.au† Present address: Brent D. Newman, Isotope Hydrology Section, International Atomic Energy Agency, PO Box 100, Wagramer Strasse 5, Vienna,
A1400, Austria.
Copyright  2010 John Wiley & Sons, Ltd.

soil wetness, the second mechanism may become moresignificant. Rasse et al. (2000) demonstrate that Alfalfa
root-induced macropores increase saturated soil hydraulicconductivity by 57%, and facilitate deep percolation.
Whether the roots impede (by root water uptake) or
facilitate (by root-induced macropores) deep percolation
depends on various conditions including climate, soil, andvegetation characteristics (Newman et al., 2006), particularly the frequency and duration at which macroporesare saturated or nearly saturated (Dong et al., 2003).
Because of increased precipitation and reduced potential evapotranspiration in the mountains, the potential for
local water yield increases. Under a conventional assumption of impermeable underlying bedrock, this water yield
will contribute to lateral surface runoff and/or interflow
along the hillslope. Recent experiments at a humid hillslope site challenge the impermeable-bedrock assumption
with observations that bedrock percolation can be a predominant hillslope water loss, even for granite (Trompvan Meerveld et al., 2007). In most situations, bedrock
percolation depends on whether, and how quickly, soil
moisture can be transmitted from the root zone into thesoil-bedrock interface.
Root induced macropores have been recognized to beeffective soil–water conduits in hillslopes (Beven and1096

H. GUAN ET AL.

Figure 1. Location map of the study site (adapted from Wilcox et al., 1997)

Germann, 1982; Noguchi et al., 1999). Thus, a deeper
root zone can be advantageous for downward water
movement into the soil profile. However, a deeper rootzone would also amplify root-water-uptake functions,
reducing the potential for deep percolation. Thus, thefirst objective of this paper is to examine how root
distribution in a semiarid hillslope affects soil waterpercolation into the underlying permeable bedrock bynumerical modelling.
The study is based on experiments conducted at aponderosa pine (Pinus ponderosa) site in northern central New Mexico, USA. The ponderosa pine site has alow slope (6%), with highly permeable tuff bedrock (saturated hydraulic conductivity, Ks D 1Ð8 ð 104 mm/year,
or permeability k D 5Ð8 ð 1014 m2 ). Field observations
suggest only negligible percolation into the bedrock(Wilcox et al., 1997). At this same site, slope–parallel
flow through root macropores was observed in a clayrich soil horizon during the snowmelt seasons (Newman et al., 2004), and during some rain events (Wilcox
et al., 1997 and Newman et al., 1998). The sum of surface runoff and interflow through the soil was estimatedto be 10–60 mm/year (Wilcox et al., 1997). The restof the water was returned to the atmosphere throughevapotranspiration (ET) (Brandes and Wilcox, 2000). A
low-permeability barrier at or above the soil–tuff interface was hypothesized to impede downward movementof water into the highly permeable tuff (Wilcox et al.,1997). However, based on field observations of a transect across the hillslope, Newman et al. (2004) suggested
that the less permeable soil matrix of the root zone
and root-induced macropores together explained the lowquantity of bedrock percolation at this site. In contrast to
Copyright  2010 John Wiley & Sons, Ltd.

the conventional notion that the soil is saturated upwards
from the soil-bedrock interface during a wetting event,they proposed that at this specific site the soil is saturated
from the top of the soil profile by root macropores. Dueto the lack of vegetation roots in the soil horizons justabove the permeable tuff, wetting to the bedrock rarely
happens.
With this in mind, the second objective of this paper
is to examine these interpretations regarding low levelsof percolation into the highly permeable bedrock at theponderosa pine hillslope site by investigating water partitioning using numerical simulations. Both root-wateruptake and root-macropore flow processes are consideredin these simulations.

THE STUDY SITE AND DATA
The study site is located at an elevation of 2315 m, andcovered by sparse ponderosa pine trees with an understory of short grasses (Figure 1). The southeast-facing
hillslope has a steepness of about 6%. The soil cover isfairly uniform, about 100 cm thick, with A, Bw, Bt, andCB horizons (Wilcox et al., 1997). Roots occur primarily
in the top 70 cm of the soil in the A, Bw, and Bt horizons(Newman et al., 2004). The estimated soil hydraulicproperties from soil core measurements are given in
Table I (columns 2 through 6). Micrometeorological data
were measured at 15-min intervals, with instruments
installed at a height of about 2 m above the ground. The
data included precipitation, solar radiation, wind speed,
relative humidity, and temperature. Precipitation wasmeasured at 1 min intervals. Surface runoff and interflow
were monitored on-site every 15 min, and measured at a
Hydrol. Process. 24, 1095– 1105 (2010)

1097

HILLSLOPE HYDROLOGICAL SIMULATION

Table I. van Genuchten (1980) model parameters and saturated hydraulic conductivities of the four soil horizons and the underlying
bedrock
Horizon˛a(1/ cm)

na

s a

r a

Kmatrix c(m/s)

Ks – x d(m/s)

Ks – z e
(m/s)

A
Bw
Bt
CB
R (tuff)

0Ð015
0Ð017
0Ð00450Ð016
0Ð0014

1Ð33
1Ð14
1Ð15
1Ð11
1Ð42

0Ð44
0Ð39
0Ð44b
0Ð47
0Ð28

0Ð06
0Ð06
0Ð08
0Ð06
0

7Ð5 ð 1075Ð7 ð 109
4Ð0 ð 109
3Ð0 ð 1095Ð8 ð 1079Ð4 ð 106
8Ð7 ð 106
9Ð4 ð 106
3Ð5 ð 1095Ð8 ð 107

3Ð1 ð 1062Ð3 ð 1062Ð5 ð 106
3Ð2 ð 1095Ð8 ð 107

a The van Genuchten model parameters.
b The saturated volumetric water content was changed from the lab-measured value of 0Ð4–0Ð44 to match the observed saturated soil water content

in the field.
c The lab-measured saturated hydraulic conductivity of the soils at various depths, and bedrock near the study site.
d,e The bulk soil saturated hydraulic conductivity in the direction paralleld and perpendiculare to the slope surface. The root perimeter aperture wasassumed to be 2Ð7% of the root diameter and the equivalent root dip angle to be 15° from the surface (calibrated based on surface runoff).

Figure 2. The atmospheric boundary condition for three water years: 1994 (left), 1995 (middle), and 1996 (right) at the ponderosa pine site. The
months are listed along the horizontal axis. Only daily data are presented here. The temporal resolution of the data, which was used in HYDRUS
simulations, is every 15 min for 1994 and 1995, and every 60 min for 1996. Precipitation is presented here as the sum of rain and snowmelt

1 min interval when there was flow. Soil moisture wasmeasured a few tens of times at various locations anddepths during each year.
The atmospheric boundary condition was generated
for the numerical simulations described in the next
section. The precipitation (P) data (heated gauge data)
were used to generate a time-series of water (rainfall C
snowmelt) readily available for infiltration and runoff.
The total precipitation from October 1993 through August
1998 was 2590 mm, with an annual mean of 527 mm.
The recorded precipitation was classified into rainfall
and snowfall based on the mean daily temperature(Wigmosta et al., 1994). In this study, we used 0 ° C for
the rainfall and snowfall separation. The snowmelt wasestimated using an empirical temperature-index model,
in which the snowmelt is linearly related to the daily
temperature above a certain threshold (0 ° C was used).
The potential evapotranspiration (PET) was calculatedusing the Penman-Monteith equation for a hypothetic
reference grass using the collected data of solar radiation,
temperature, and relative humidity. Since wind speed was
not measured above the ponderosa trees, a mean valueof 2 m/s was used (Allen et al., 1998). The potential
ET during the night (8 pm to 6 am) was assumed to
Copyright  2010 John Wiley & Sons, Ltd.

be zero. The estimated total PET from October 1993
through August 1998 was 6486 mm, with an annual meanof 1319 mm, about 2Ð5 times the precipitation during thisperiod. Data of three water years (1994–1996, each fromthe previous October through September) are shown in
Figure 2. The seasonality of PET was similar from year toyear, while rainfall and snowmelt varied between years.
The ratios of PET/P were 2Ð5, 2, and 3 for 1994, 1995,
and 1996 water years, respectively.

NUMERICAL SIMULATIONS
Since root macropore flow was observed (Newman et al.,2004), and ET represents a dominant water balance component at this study site, both are important processes,
and need to be considered in the simulations. Various computer codes with differing conceptual models
of fracture flows of various complexity (Altman et al.,
1996; Simunek et al., 2003), can be used to simulate
root-macropore flows. However, since dynamic variations in ET also have to be considered, the choice ofcomputer codes becomes rather limited. HYDRUS is oneof a few models that considers both ET and macropore
flow processes simultaneously. In order to examine root
Hydrol. Process. 24, 1095– 1105 (2010)

1098

H. GUAN ET AL.

macropore flow in two directions (parallel and orthogonal to the slope surface), a two-dimensional version of
HYDRUS, called HYDRUS-2D, is used in this study. In
HYDRUS, variably saturated water flow in porous media
is described using the Richards equation. Evaporation
is modelled using a Darcy’s law-based extraction function, and transpiration is calculated as a root-distributionweighted sink term, depending on root-zone soil-water
potential and potential transpiration. With its capacity todynamically relate near-surface water partitioning processes (including evaporation, transpiration, soil waterflow) to soil hydraulic properties, root distribution, andthe root-water-uptake function, HYDRUS allows us to
examine soil and vegetation effects on water partitioningat the hillslope site.
Root macropore modelling
It is challenging to represent the hydrological effectsof roots on soil water flow. The difficulties include(1) how to represent the quantitative contribution of anindividual root to an increase in the bulk soil hydraulicconductivity, (2) how this contribution varies with soil
water potential, (3) how this contribution varies with thedirection (e.g. parallel vs. orthogonal to the slope), and(4) how the bulk root-induced hydraulic conductivity isquantitatively related to the observed root distribution.
Dong et al. (2003) estimated macropore contribution to
the bulk soil hydraulic conductivity based on bulk soil
water potential and macropore size distribution, which
addresses issue (2) quite well. It was developed forone-dimension vertical flow, and not appropriate forsimulating water movement on hillslope. In spite of theabove challenges, when the boundary conditions anda flow field are known, it is less difficult to estimate
macropore equivalent hydraulic properties, which can be
defined as those properties subjected to the boundary
conditions (Renard and deMarsily, 1997). We will apply
the equivalent-property approach to quantify the root
macropore effects on the studied hillslope. A similar
approach has been used to model subsurface tile-drainedsoils (Carlier et al., 2007), where equivalent hydraulicconductivity is calculated with drainage density and drain
radius.
The procedure to estimate the equivalent anisotropic
hydraulic conductivities of a soil with root macroporeson a sloping surface, which is related to the observed root
distribution, is based on several assumptions. First, thatthe root macropore contribution to the saturated hydraulicconductivity is quantitatively related to root distribution(e.g. root size and density). Second, that the anisotropy ofsaturated root macropore flow is quantitatively related toan average root dip angle. Third, that the bulk hydraulicconductivity at various saturations is evaluated using
the Mualem (1976) pore size distribution model andthe retention properties of the soil matrix. Using asingle equivalent porous medium for both macropore
and matrix flow leads to a conceptual flaw in that themacropore flow does not cease, as it should, when the
Copyright  2010 John Wiley & Sons, Ltd.

macropore water is drained below a certain degree, but
is only reduced by a factor that is related to soil matrix
properties. On the other hand, for the same reason, thesingle equivalent continuum model allows for macropore
flow to occur before the soil matrix gets fully saturated,which was commonly observed at the field site (Newman
et al., 1998, 2004). Numerically, the single continuum
model should work appropriately in simulating waterflow at saturated and near-saturated conditions, when themacropore flow is most important.
To formulate the equivalent single continuum model,
an individual root macropore is considered as an annularvoid around the root, with its aperture proportional to
the root size (Figure 3). Following Pneuli and Gutfinger(1992), the saturated hydraulic conductivity (mÐs1 ) in
the longitudinal direction of the slope due to one root, i,parallel with the slope is
Ki D
and



Ci D 1 cos2  gDi C 2bi 4
Ci
Ai 
128

Di
Di C 2bi4[1  Di /Di C 2bi 2 ]2lnDi C 2bi /Di 

1


2
where  (radian) is the angle between the root and the longitudinal slope direction x (Figure 3),  (dimensionless)is the tortuosity of the root, which is defined as the ratioof the actual root length over the linear length, Di (m),bi (m), and Ai (m2 ) are the root diameter, the aperture
width, and the cross-sectional area (including aperture)of the root, and  (kgÐm3 ), g (mÐs2 ), and  (PaÐs) are
the density of liquid water, gravitational acceleration, andthe dynamic viscosity of water, respectively.
For a given layer with an isotropic root distribution
in the plane parallel to the slope surface, the rootinduced saturated hydraulic conductivity is an summationof (1) over different angles  (from /2 to /2) for all
different root sizes:[Cj Dj C 2bj 4 nj ]grootclassj
Kroot D
3256
[1/4nj Dj C 2bj 2 ]rootclassjwhere nj is the number of roots in one root class j with amean root diameter of Dj observed in the transect, Cj isdefined in Equation (2). Equation (3) gives the saturated
hydraulic conductivity in the longitudinal direction (x) of
the slope due to root macropores, assuming that all roots
are in planes parallel to the slope surface.
However, the roots may intercept the slope at a certain
angle rather than being parallel to the slope (Figure 3).
An equivalent root dip angle is used to lump all roots.
In areas with similar surface vegetation coverage andsoil water conditions, the equivalent root dip angles
should be similar. In this situation, root macropores also
contribute to water flow in the z direction, where z is
the downward direction orthogonal to x. A contribution
Hydrol. Process. 24, 1095– 1105 (2010)

HILLSLOPE HYDROLOGICAL SIMULATION

1099

Figure 3. Conceptualization of the hydraulic effects of root macropores, where D is the root diameter, b is the annular void between the root surface
and surrounding soil matrix, x is the longitudinal direction parallel to the slope surface, and z is the downward direction perpendicular to x, ˇ isaverage root dip angle, and  is the angle between individual root projection on the slope surface and the x directionof root macropores to the hydraulic conductivity in the z
direction is as follows:[Cj Dj C 2bj 4 nj ]g sin ˇ
Kzroot D
4
256
[1/4nj Dj C 2bj 2 ]where ˇ (radian) is the equivalent root dip angle withrespect to the slope surface. The saturated hydraulicconductivity in the x direction has to be modified as well,resulting in[Cj Dj C 2bj 4 nj ]g cos ˇKxroot D
5
256
[1/4nj Dj C 2bj 2 ]The bulk saturated hydraulic conductivities (mÐs1 ) inx and z directions are then given by
Kzbulk D Kzroot nroot C Kmatrix 1  nroot 6

Kxbulk D Kxroot nroot C Kmatrix 1  nroot 

7

and

where nroot is the volumetric fraction occupied by theroots and their surrounding voids, defined as
nroot D

size
class jD1
rj Dj C 2bj 2 
4

8

where r j (m2 ) is the root density of the root-size classj and Kmatrix the soil matrix saturated hydraulic conductivity. The root density (r ) is measured in transects
orthogonal to the sloping direction and reported as a number of roots per unit area for each root-size class.
Copyright  2010 John Wiley & Sons, Ltd.

Evaporation and transpiration modelling
It is estimated that over 90% of soil water loss at the
study site was due to ET (Brandes and Wilcox, 2000).
It is therefore important to appropriately represent evaporation and transpiration processes in order to simulate
soil water partitioning correctly. Evaporation and transpiration processes were modelled separately, according
to the surface vegetation coverage. Potential evaporation(PE) and potential transpiration (PT) are input to themodel, defining the maximum rate of simulated evaporation and transpiration. In HYDRUS, evaporation is
calculated from PE dependent of the surface soil moisture availability. A critical water potential is prescribed todetermine the threshold above which evaporation occurs
at the rate of PE. When the surface soil water potential is below (drier than) this critical value, the boundary condition is set equal to this critical water potential
and the actual evaporation is determined by how quickly
the soil moisture moves toward the surface according toDarcy’s law. Apparently, simulated evaporation is sensitive to the prescribed critical value. According to Rassam
et al. (2003), for silty soils in the study area, a value of1000 m was applied.
In HYDRUS, transpiration is calculated from PT
multiplied by a scaling function dependent of root-zone
soil moisture availability. An S-shape root-water-uptake
function (the scaling function) (van Genuchten, 1987)
was used1˛D
9
1 C h/h50 p
where ˛ (dimensionless) is the ratio of the actual rootwater-uptake rate to its potential rate, in which thepotential rate is determined from PT and root distributionfunction, h (m) the soil water potentials, h50 (m) the
Hydrol. Process. 24, 1095– 1105 (2010)

1100

H. GUAN ET AL.

soil water potential at which ˛ is reduced by 50%, andp (dimensionless) an experimental constant. A wiltingpoint soil water potential, below which ˛ is set to zero, isprescribed to mimic the termination of transpiration due
to extreme soil moisture stress.
Separated PE and PT inputs are required for evaporation and transpiration modelling. Partitioning of PET into
potential evaporation (PE D 70% PET) and potential transpiration (PT D 30% PET) was roughly estimated basedon the fractional surface coverage of the ponderosa pinetrees at the study site. In addition to the ponderosa pine,
sub-story grass also contributes to transpiration. However, since the rooting depth of grass is relatively shallow
when compared to pine roots, grass transpiration can betechnically mimicked using evaporation. The root density
was calculated according to Newman et al. (2004), with
the A and Bw horizons containing about 65% of roots
and the Bt horizon the remaining 35%.
Numerical set-up
The simulations were based on data at location A on
the north hillslope (Figure 1). A two-dimension hillslope
section (Figure 4), 1 m wide and 2 m deep, was used
for the numerical simulations. The soil profile consistedof four soil horizons (A: 10 cm, Bw: 20 cm, Bt: 40 cm,
and CB: 30 cm thick) according to Wilcox et al. (1997)
and had a root zone in the top three soil horizons (Newman et al., 2004). The soil hydraulic properties (Table I)
and root distribution were assumed to be homogeneousin each soil horizon. A prescribed constant gradient,0Ð5 sin(2), where  (radian) is the slope angle, wasassigned to each side of the soil profile, mimicking a continuous hillslope condition. Runoff was modelled as aninfiltration-excess process consistent with the field resultsof Wilcox et al. (1997), without considering downslope
propagation. An atmospheric boundary condition with aninterval of 15–60 min was prescribed at the top end of
the flow domain, while a free drainage boundary condition was applied at the bottom of the domain (Figure 4).
The parameterization of the S-shape root-water-uptake
model was manually adjusted to fit the measured soil
water contents. Since root macropores strongly influence the soil infiltration capacity, simulations with bulk
soil hydraulic conductivities considering root-macropore
effects were tested against the observed runoff, to findthe appropriate parameterization of the root macropore
model.

RESULTS AND DISCUSSION
Sensitivity analysis of the root macropore modeland simulation results
A sensitivity analysis of the root macropore model
was conducted against the matrix hydraulic conductivity,
the size of the root-macropore aperture, and the root
dip angle (Figure 5). The relative contribution of root
macropores to the bulk hydraulic conductivity is related
to soil and root density, size, and direction. To examine
Copyright  2010 John Wiley & Sons, Ltd.

Figure 4. A schematic graph showing the two-dimension, 1-m-wide
hillslope section used for the simulations

the effect of the soil matrix, we used the root distribution
observed in the A-horizon, a root dip angle of 15° , andan aperture of 2Ð7% root diameter (this value is estimatedfrom calibration as described in the next paragraph). The
results show that the contribution of root macropores
dominates the bulk soil hydraulic conductivity when
Kmatrix is smaller than 1 ð 107 m/s, while it becomesinsignificant when Kmatrix is larger than 1 ð 105 m/s(Figure 5a). The actual Kmatrix of the A-horizon is 7Ð5 ð
107 m/s (Table I), indicating that both the soil matrix
and root macropores are important for infiltration at
the study site. Parameterization of the root macropore
model should thus be sensitive to the infiltration-excess
runoff. The bulk hydraulic conductivity in the z direction(Figure 3) was calculated using the laboratory-measured
soil hydraulic properties and observed root density forvarious root-macropore apertures (Figure 5b) and dip
angles (Figure 5c). The bulk hydraulic conductivity of the
CB horizon was not sensitive to root parameters, because
only a few roots were present at this depth. On the other
hand, the bulk hydraulic conductivity for the top threehorizons was sensitive to both root factors.
Consequently, calibration of root parameters (dip angle
and aperture) against the observed runoff did not provide
a unique solution. To overcome the problem of nonuniqueness, a root dip angle of 15° was fixed andonly the root aperture was optimized. An angle of
15° should be close to the actual root distribution in
the near-surface, where horizontal roots are dominant.
Although different root dip angles could be used at
different depths, only one root dip angle was used in
the simulations, due to the difficulty with the calibration
process. This may slightly affect the simulated valuesof interflow and percolation. Root-macropore apertureis assumed to be proportional to root diameter. Several
different aperture sizes were applied for the simulations
for the 1994 water year. The one of 2Ð7% root diameter,
Hydrol. Process. 24, 1095– 1105 (2010)

HILLSLOPE HYDROLOGICAL SIMULATION

1101

but misses the first one in the year (Figure 6). Similar
simulations, without considering root macropore effects,
were performed, with the result of a modelled runoff overfour times of the observation in the 1994 water year
(not shown). With the root-induced macropore effects
considered, the modelled runoff is about the same amountof the observation. The difference between the modelled
and observed first runoff event could be due to either
modelling failure or the heterogeneity effect in rainfall
and runoff generation between the small catchment andthe simulated location. One common problem leading
to failure in runoff simulations could be resulted froma low temporal-resolution precipitation input. Because
precipitation was input at 15-min resolution, this problem
is less likely. Other factors, such as rainfall size, and dry
interval, do not explain the apparent modelling failure
either, because similar situations occurs either in runoff
event (2) or (3).
After calibration with the 1994 water year data, simulated volumetric water content fairly well matched theobservations in all 3 years (Figure 7). The abrupt change
in the modelled volumetric soil water content at 1 m
depth is due to the prescribed material change from thesoil to the tuff. Although the soil thickness was observed
about 1 m thick, it varies slightly from one location toanother. Moreover, the soil-tuff interface is not flat andabrupt. Some surface relief and/or weathering could leadto a soil-moisture transition zone between the soil andunderlying tuff. Thus, it is appropriate to compare modelled volumetric moisture content below 1 m depth tothe observed moisture content below the transition zone
(Figure 7). Big deviation between the modelling results
and observations in the 180th day of 1995 water year
are probably related to a recent focused snowmelt event,which was not included in the simulation. In preparing theatmospheric boundary conditions, snow was assumed to
be uniformly distributed on the surface. In reality, snow
might have piled locally, and melt to feed a small area ofthe catchment. If this happened near location A (Figure 1)
some time before the 180th day, it might saturate the soil.
The saturated soil in depth could persist to the 180th dayof the water year (Figure 7). This possibility is supportedby over 0Ð45 volumetric water content, almost at saturation, in comparison to 0Ð4 in the other days.

Figure 5. Sensitivity analysis of the root macropore model. (a) The saturated hydraulic conductivity (Kz ) perpendicular to the slope surface asa function of the soil matrix hydraulic conductivity. Kz was calculatedusing the root distribution data observed in the A horizon. The two solid
symbols represent the measured A-horizon matrix hydraulic conductivities. (b) Kz calculated as a function of the root-macropore aperture, which
is defined as the ratio between b and D (Figure 3), for an equivalent root
dip angle of 15° . (c) Kz calculated as a function of the root dip angle forthe root aperture of 0Ð027

which best fits the observed runoff, was chosen. The
simulation reproduces two of the three major infiltrationexcess runoff events observed on the north hillslope,
Copyright  2010 John Wiley & Sons, Ltd.

Root-water-uptake model calibration and soil moisturepartitioning
The root-water-uptake model was calibrated against
the root-zone soil moisture for the 1994 water year
(Figure 7). The wilting point water potential was prescribed to be 250 m, close to the observed root vulnerability curve (percentage loss of xylem conductivity vs.
water potential) (Stout and Sala, 2003). The calibration
process produced h50 of 5 m and p of 2 (Equation 9),
however, this is quite different from the reported vulnerability curve (Stout and Sala, 2003). If we assume that theroot vulnerability curve can be used for root-water-uptake
model, the fitted value of h50 for the ponderosa pine
Hydrol. Process. 24, 1095– 1105 (2010)

1102

H. GUAN ET AL.

Figure 6. Observed and simulated major daily runoff events in the 1994 water year (Oct 1, 1993— Sept 30, 1994). Daily precipitation (rain C
snowmelt) time series are shown for comparison

would be 80 m and of p equal to 2Ð7. Using this parameterization, simulated root-zone soil moisture is muchdrier than that observed. A few reasons may explain theinconsistency between the calibrated h50 and that determined from the measured root vulnerability curve. First,
the transpiration process is affected by the hydraulic conductivity of the root, stem, and leaf conductance. As thevulnerability curve represents only the response of root
and/or stem to the moisture stress, it does not include
the vegetation regulation of transpiration by leaf stomata.
Thus, the parameterization obtained from the vulnerability curve represents only the lower bound for h50 . Whenleaf regulation is also considered, h50 may occur in muchwetter conditions. Second, as the effect of osmotic potential was not considered in the simulation, h50 calibrated
against the observed soil moisture can represent slightly
wetter conditions than the actual soil water potential,
which includes both the matrix and osmotic potentials.
Third, as the two parameters (h50 and p) were calibrated
simultaneously, the results could be nonunique and deviate from the real values. However, the latter two effects
generally do not play a big role. Therefore, the first effectis likely the main reason for the difference between thetwo sets of obtained root water uptake parameters. This
indicates that the root (or stem) xylem vulnerability curve
is probably not appropriate for parameterizing root-wateruptake models.
The 1994 water year had an annual precipitation close
to average long-term conditions. The initial condition
for the 1994 water-year simulation was obtained byrerunning the same 1994 climate forcing several times
to achieve a quasi-steady state. With root-macroporeand ET models calibrated against observations from the1994 water year, additional simulations were conducted
to examine the soil water partitioning for the 1995
and 1996 water years (Figure 7 and Table II). Initialconditions for the 1995 and 1996 water years were basedon simulations of the previous water years. The root
Copyright  2010 John Wiley & Sons, Ltd.

mean square errors of the simulated soil moisture at 5
depths for 20 observation times are 0Ð09, 0Ð08, and 0Ð07,
in comparison to mean observed soil moisture of 0Ð30,
0Ð35, and 0Ð28, respectively, for the three water years.
While annual percolation into the tuff at the ponderosasite was simulated to be less than 5 mm (or 1% of
precipitation), simulated evapotranspiration accounted for
94% of annual precipitation in 1994 (Table II, column#2). Percolation into the tuff in both 1995 and 1996
was similar to that of the 1994 water year, i.e. lessthan 1% of annual precipitation (Table II, columns #2 ¾
4). From these simulations, it appears that groundwater
recharge (percolation) at this location is not significant,which is consistent with the previous study of this site
(Wilcox et al., 1997). The modelled average ET for the
3-year period is 95% of the total precipitation, which isin agreement with previous water balance analyses ofthe field data (Brandes and Wilcox, 2000). Simulated
interflow during the 3 years is of trace amount, similar
to what was observed (Wilcox et al., 1997).
The effects of impeding layer and root macropore
on groundwater recharge
Wilcox et al. (1997) suggested that an impeding layerof low permeability may restrict downward water movement into the bedrock at the ponderosa pine site. At
this site, an abrupt change in soil matrix permeability
occurs at a depth of about 10 ¾ 20 cm in the Bw horizon(Table I). However, soil water in the field was observed
to pond at a depth of 70 cm, around the interface between
the Bt and CB horizons. Apparently, the soil hydraulicconductivity of the matrix cannot explain a hypothetical
impeding layer at this site. Wilcox et al. (1997) hypothesized three potential impeding layers: (1) CB horizon,(2) the base of the Bt horizon, and (3) a thin ‘smear’ oftranslocated clay at the soil–tuff interface. Newman et al.
(2004) suggested that the absence of root macropores at
Hydrol. Process. 24, 1095– 1105 (2010)

HILLSLOPE HYDROLOGICAL SIMULATION

1103

Figure 7. Simulated (lines) and observed (open diamonds) soil water contents (dimensionless) versus the depth (cm) at selected days (the number isgiven in each panel) for three water years. Selected days are either following the snowmelt season or in the summer monsoon season. Big deviations
between modelling results and observations in 1995(180) are probably due to a recent focused snowmelt near location A, which was not represented
in the simulations

Table II. Simulated water fluxes (in equivalent depths ofmm/year) for the hillslope
Water year

Precipitation (P)Infiltration
ET
RunoffPercolation
ET/P
Percolation/P

Observed soil profile
1994

1995

1996

Hypotheticalsoil profilea
1994

524
506
490
18
<5
93%
<1%

686
637
598
49
<5
87%
<1%

467
467
435
0
<5
93%
<1%

524
506
379
18
108
72%
21%

a It was assumed that the root zone directly contacts the underlying tuff
by replacing the CB horizon with tuff.

depths below 70 cm could explain an impeding layer atthis depth.
Because of modifications due to root-induced macropores, the clay-rich soil of the root zone has a vertical
Copyright  2010 John Wiley & Sons, Ltd.

bulk saturated hydraulic conductivity close to that of theloam soil (Table I). With fewer roots in the CB horizon(or the bottom of the Bt horizon), the soil hydraulicconductivity is about three orders of magnitude lower.
The root-induced macropore modification of the bulk
saturated hydraulic conductivity of this clay-rich soil (bya factor of 103 ) is much larger than the enhancement
factor of 1Ð5 found by Rasse et al. (2000) at an alfalfa
study site. This difference can be explained by thedifference in Kmatrix between the two cases, which isclearly shown in Figure 5a. The Kmatrix in the Rasseet al. (2000) study was high (3Ð2 ð 106 m/s), while in
the Bw and Bt horizons of this study, the Kmatrix is 5 ð109 m/s. With root-macropore modification considered,
the permeability contrast is located between Bt andCB horizons, or at a depth of 70 cm. This explains
the observed vertical distribution of soil moisture, andsupports Newman et al’s. (2004) interpretation and thefirst two hypotheses of Wilcox et al. (1997), because theHydrol. Process. 24, 1095– 1105 (2010)

1104

H. GUAN ET AL.

CB horizon and the base of the Bt horizon play the samerole provided that root macropores are absent at these
depths. Since the existence of the third potential impedinglayer at the soil-bedrock interface, located at the 100-cm
depth, was not supported by the vertical soil moisture
distribution (Figure 7), it was not further examined here.
To further test whether the absence of root macropores
was the cause for low percolation into the highly permeable bedrock, two additional simulations with a modified
CB horizon were conducted. In the first simulation, theCB horizon was replaced with tuff, leading to a 70-cm
thick soil with root macropores in direct contact with theunderlying tuff and resulting in significantly enhanced
percolation (108 mm, Table II, column #5). In the second
simulation, the CB horizon was replaced with a Bt horizon that included root macropores, resulting in a similarly
enhanced percolation (105 mm, not shown). Optionally,
we could add root macropores to the CB horizon. As
the matrix hydraulic conductivity of the Bt and CB horizons is similar, the resulted annual percolation would be
similar to the second simulation. Simulated percolation,
when the root zone is in direct contact with the tuff,
increases to approximately 20% of annual precipitation.
This suggests that the low-permeability soil that lacks
vegetation roots behaves as an impeding layer. It prevents
significant percolation into the bedrock. If the ponderosa
hillslope had a thinner soil cover, with the root zone indirect contact with the underlying permeable bedrock,the potential distributed recharge (percolation) could be
as large as 20% of annual precipitation. This result indicates that the diffuse recharge at this environment, andother similar environments with low-permeability soil
overlying permeable bedrock, can be highly variable,
depending on the combination of soil thickness and root
distribution.

CONCLUSIONSVegetation impacts soil moisture partitioning by two
processes: root-water-uptake and root-induced macropore flow. The apertures surrounding vegetation roots canbehave as macropores, or serve as preferential flow paths.
An equivalent-continuum model is presented to quantify
the root-induced macropore effect. The model-calculatedroot-macropore enhancement factors are consistent withthe observations at the semiarid ponderosa pine site ofthis study, and at the alfalfa site reported by Rasse et al.
(2000). The results indicate that the root macropore effectdominates saturated soil water flow in low permeablesoils (Kmatrix below 107 m/s), while it becomes insignificant in soils with Kmatrix larger than 105 m/s.
At the ponderosa pine site, the model with simple rootmacropore parameterization reasonably well reproducessoil moisture distribution and some major runoff events.
The results indicate that the clay-rich soil layer withoutroot-induced macropores behaves as an impeding layer
for potential groundwater recharge. This impeding layer
results in a bedrock percolation of less than 1% of
Copyright  2010 John Wiley & Sons, Ltd.

the annual precipitation. Without this impeding layer,
percolation into the underlying permeable bedrock canbe as much as 20% of annual precipitation. These
results support the first two potential impeding layers
hypothesized by Wilcox et al. (1997), as well as theinterpretation for little deep percolation at this site byNewman et al. (2004).

ACKNOWLEDGEMENT

This work was primarily supported by the SAHRA (Sustainability of semi-Arid Hydrology and Riparian Areas)
under the STC Program of the National Science Foundation, Agreement EAR-9876800, and partially supported
by Flinders University Establishment Funding.

REFERENCES
Allen RG, Pereira LS, Raes D, Smith M. 1998. Crop evapotranspirationGuidelines for computing crop water requirements. Food andAgriculture Organization of the United Nations Irrigation and Drainage
Paper 56.
Altman SJ, Arnold BW, Barnard RW, Barr GW, Ho CK, McKenna SA,
Eaton RR. 1996. Flow Calculations for Yucca Mountain Groundwater
Travel Time (GWTT-95). Sandia National Laboratory: Albuquerque.
Beven K, Germann P. 1982. Macropores and water flow in soils. Water
Resources Research 18(5): 1311– 1325.
Brandes D, Wilcox BP. 2000. Evapotranspiration and soil moisturedynamics on a semiarid ponderosa pine hillslope. Journal of theAmerican Water Resources Association 36(5): 965–974.
Carlier JP, Kao C, Ginzburg I. 2007. Field-scale modeling of subsurface
tile-drained soils using an equivalent-medium approach. Journal of
Hydrology 341(1– 2): 105–115.
Dong WQ, Yu ZB, Weber D. 2003. Simulations on soil water variationin arid regions. Journal of Hydrology 275(3– 4): 162– 181.
Flint AL, Flint LE, Hevesi JA, Blainey JM. 2004. Fundamental conceptsof recharge in the Desert Southwest: a regional modeling perspective.
In Groundwater Recharge in a Desert Environment: The Southwestern
United States, Water Science and Applications Series, Hogan JF,
Phillips FM, Scanlon BR (eds). American Geophysical Union:
Washington, DC; 159– 184.
Gee GW, Wierenga PJ, Andraski BJ, Young MH, Fayer MJ, Rockhold ML. 1994. Variations in water-balance and recharge potential at
3 Western Desert Sites. Soil Science Society of America Journal 58(1):63–72.
Izbicki JA, Radyk J, Michel RL. 2000. Water movement through st thickunsaturated zone underlying an intermittent stream in the western
Mojave Desert, southern California, USA. Journal of Hydrology
238(3– 4): 194–217.
Mualem Y. 1976. New model for predicting hydraulic conductivity ofunsaturated porous-media. Water Resources Research 12(3): 513–522.
Newman BD, Campbell AR, Wilcox BP. 1998. Lateral subsurface flow
pathways in a semiarid ponderosa pine hillslope. Water Resources
Research 34(12): 3485– 3496.
Newman BD, Wilcox BP, Archer SR, Breshears DD, Dahm CN, Duffy
CJ, McDowell NG, Phillips FM, Scanlon BR, Vivoni ER. 2006.
Ecohydrology of water-limited environments: a scientific vision. Water
Resources Research 42(6): 15.
Newman BD, Wilcox BP, Graham RC. 2004. Snowmelt-driven macropore flow and soil saturation in a semiarid forest. Hydrological Processes
18(5): 1035– 1042.
Noguchi S, Tsuboyama Y, Sidle RC, Hosoda I. 1999. Morphologicalcharacteristics of macropores and the distribution of preferential flow
pathways in a forested slope segment. Soil Science Society of America
Journal 63(5): 1413– 1423.
Phillips FM. 1994. Environmental tracers for water movement in desert
soils of the American southwest. Soil Science Society of America
Journal 58(1): 15–24.
Pneuli D, Gutfinger C. 1992. Fluid Mechanics. Cambridge University
Press: Cambridge, United Kingdom; 496.
Rassam D, Simunek J, van Genuchten MT. 2003. Modelling Variably
Saturated Flow with HYDRUS-2D. ND Consult, Brisbane, Australia.
Hydrol. Process. 24, 1095– 1105 (2010)

HILLSLOPE HYDROLOGICAL SIMULATION
Rasse DP, Smucker AJM, Santos D. 2000. Alfalfa root and shootmulching effects on soil hydraulic properties and aggregation. Soil
Science Society of America Journal 64(2): 725–731.
Renard P, deMarsily G. 1997. Calculating equivalent permeability: areview. Advances in Water Resources 20(5– 6): 253– 278.
Scanlon BR, Levitt DG, Reedy RC, Keese KE, Sully MJ. 2005.
Ecological controls on water-cycle response to climate variability indeserts. Proceedings of the National Academy of Sciences of the UnitedStates of America 102(17): 6033– 6038.
Simunek J, Jarvis NJ, van Genuchten MT, Gardenas A. 2003. Review
and comparison of models for describing non-equilibrium andperferential flow and transport in the vadose zone. Journal of Hydrology
272: 14–35.
Stothoff SA, Or D, Groeneveld DP, Jones SB. 1999. The effect ofvegetation on infiltration in shallow soils underlain by fissured bedrock.
Journal of Hydrology 218(3– 4): 169– 190.
Stout DL, Sala A. 2003. Xylem vulnerability to cavitation in Pseudotsuga
menziesii and Pinus ponderosa from contrasting habitats. Tree
Physiology 23(1): 43–50.
Tromp-van Meerveld HJ, Peters NE, McDonnell JJ. 2007. Effect of
bedrock permeability on subsurface stormflow and the water balanceof a trenched hillslope at the Panola Mountain Research Watershed,
Georgia, USA. Hydrological Processes 21(6): 750–769.

Copyright  2010 John Wiley & Sons, Ltd.

1105

van Genuchten MT. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal 44(5): 892– 898.
van Genuchten MT. 1987. A Numerical Model for Water and Solute
Movement in and Below the Root Zone. U.S. Salinity laboratory:
Riverside, CA.
Walvoord MA, Phillips FM. 2004. Identifying areas of basin-floor
recharge in the Trans-Pecos region and the link to vegetation. Journalof Hydrology 292(1– 4): 59– 74.
Walvoord MA, Plummer MA, Phillips FM, Wolfsberg AV. 2002. Deep
arid system hydrodynamics— 1. Equilibrium states and response times
in thick desert vadose zones. Water Resources Research 38(12): 1308,
DOI: 10.1029/2001wr000824.
Wigmosta MS, Vail LW, Lettenmaier DP. 1994. A distributed hydrologyvegetation model for complex terrain. Water Resources Research 30(6):1665– 1679.
Wilcox BP, Newman BD, Brandes D, Davenport DW, Reid K. 1997.
Runoff from a semiarid ponderosa pine hillslope in New Mexico. Water
Resources Research 33(10): 2301– 2314.

Hydrol. Process. 24, 1095– 1105 (2010)