Numerical SimulaƟon of Flow Dynamics
during Macropore–Subsurface Drain
InteracƟons Using HYDRUS

OR

Onur Akay, Garey A. Fox,* and Jirka Šimůnek
Macropores, such as those created by deep-burrowing earthworms, have the poten al to be hydraulically connected not
only to the soil surface but also to subsurface drains. This hydraulic connec on may lead to rapid movement of surfaceapplied chemicals to receiving waters as they bypass the bulk of the soil matrix. In this study, a numerical model (HYDRUS)
that solves the three-dimensional Richards equa on for both matrix and macropore domains was used to analyze previously conducted experiments that contained a single, surface-connected or buried, ar ficial macropore and a subsurface
drain installed in a laboratory soil column. Both matrix and macropore domains were parameterized using con nuous soil
hydraulic func ons. Simula ons confirmed that surface-connected macropores were highly eﬃcient preferen al flow
paths that substan ally reduced arrival mes to the subsurface drainage outlet, with this reduc on being directly related
to the length of the macropore. Surface-connected macropores need to extend at least halfway to the drain to have ano ceable eﬀect (>50% reduc on) on the arrival me. No significant changes were observed in total drain ou lows forcolumns with laterally shi ed macropores (away from a drain) compared with centered macropores unless the macropore
depth extended significantly (>75%) into the profile. The model predicted that buried macropores became ac ve and contributed to the total ou low only when pressure heads in the soil profile became posi ve. The eﬀect of buried macroporeson drain flow was inves gated for a case where an ini ally par ally saturated profile was drained. Under these condi ons,
the numerical simula ons suggested that buried macropores could contribute up to 40% of the total ou low, which confirms laboratory observa ons with subsurface-drained soil columns with macropores.
A: EF, express frac on; NS, Nash–Sutcliﬀe.

C

about the rapid transport of contaminants, such as pesticides (Fox et al., 2004), pathogens(Joy et al., 1998; Geohring et al., 1999; Jamieson et al., 2002;
Shipitalo and Gibbs, 2005), and nutrients from the soil surface togroundwater through macropores (Magesan et al., 1995; Kladivko
et al., 1999; Shipitalo and Gibbs, 2000). Commonly observed
short-circuiting in many subsurface drainage ﬁeld studies has
been hypothesized to be due to direct hydrologic connectivity
between macropores and subsurface drains (Fox et al., 2004, 2007;
Shipitalo and Gibbs, 2000, 2005). Directly connected macropores can result in the rapid transport of contaminants from the soil
surface into the subsurface drains and then into adjacent receiving streams and channels, bypassing the soil’s ﬁltering capacity.
The ability to model the interrelationship between macropore-

facilitated contaminant transport and subsurface drainage systems,
where soil is consistently near saturation, is important for evaluating potential environmental contamination (Hoorman et al.,2005). In fact, Dorner et al. (2006) reported model simulations
that predicted that most microorganisms entering streams enterfrom tile drainage rather than overland transport.
Column studies have been performed to answer questions
regarding conditions when macropore ﬂow occurs (Phillips et al.,
1989; Trojan and Linden, 1992). Phillips et al. (1989) showedthat water under negative pressures can enter simulated macropores after ﬁrst establishing a continuous water ﬁlm on the fulllength of the macropore walls. Trojan and Linden (1992) used
10 19-L, 29-cm-diameter buckets each packed with topsoil. After
packing the buckets, 14 worms of species Aporrectodea tuberculata
Eisen were placed on the soil surface. They concluded that thequantity of water that potentially can be delivered to a burrow
is a function of initial soil water content, rainfall intensity andamount, hydraulic conductivity, and surface contributing area.
Joschko et al. (1989) allowed earthworms (Lumbricus terrestrisL.) to burrow into a soil column (19 cm wide and 40 cm long),
artiﬁcially compacted to a pore volume of 40%. They demonstrated an exponential relationship between burrow depth andsaturated hydraulic conductivity for the mainly vertically oriented
L. terrestris burrows.
Several studies have included laboratory tests using undisturbed soil samples taken from macroporous ﬁelds (Shipitalo etal., 1990) or in situ tests using pan lysimeters or isolated blocks(Shipitalo and Edwards, 1996; Shipitalo et al., 1994). Thisapproach permits an intact soil section to be analyzed underO. Akay, Water Resources Group, PBS&J, Henderson, NV, formerly Dep.
of Biosystems and Agricultural Engineering, Oklahoma State Univ.; G.A.
Fox, Dep. of Biosystems and Agricultural Engineering, Oklahoma State
Univ., 120 Agricultural Hall, S llwater, OK 74078; J. Šimůnek, Dep. ofEnvironmental Sciences, A135 Bourns Hall, Univ. of California, Riverside,
CA 92521. Received 29 Aug. 2007. *Corresponding author (garey.fox@
okstate.edu).
Vadose Zone J. 7:909–918
doi:10.2136/vzj2007.0148© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmi ed in any form or by any means, electronic or mechanical,
including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng from the publisher.

www.vadosezonejournal.org · Vol. 7, No. 3, August 2008

909

controlled conditions, but it can become diﬃcult to distinguishmacropore eﬀects from those introduced by the inherent complexity of natural soil systems.
In addition to the strategies discussed above to study macropore ﬂow, there have been laboratory studies that utilized artiﬁcial
macropores (Köhne and Mohanty, 2005; Ghodrati et al., 1999;
Li and Ghodrati, 1997; Castiglione et al., 2003), which involves
the creation of a macropore of known dimensions in a standardized soil. Although highly artiﬁcial, a standardized artiﬁcial soil
eliminates the inﬂuence of textural and structural variations on
the soil’s physical behavior and permits study of the eﬀect ofindividual macropore properties on soil physicochemical behavior(Buttle and Leigh, 1997). For example, Li and Ghodrati (1997)
showed that at given normalized ﬂuxes, preferential ﬂow tends
to occur to a greater extent in ﬁne-textured soils than coarser
textured soils, since coarser soils have much greater matrix conductivity than ﬁne soils. Preferential ﬂow also occurred underunsaturated water flow conditions. Laboratory studies with
artiﬁcial macropores can also give an important insight into theeﬀects of macropore tortuosity on breakthrough curves and solute
distribution because of their well-deﬁned macropore geometries(Allaire-Leung et al., 2000a,b).
The number of studies focused on the interaction between
macropores and subsurface drains, however, are limited (Shipitalo
and Gibbs, 2000; Shipitalo et al., 2004). Shipitalo and Gibbs(2000) observed macropores created by deep burrowing speciesof earthworms (L. terrestris) that allowed water to move directly
to subsurface drains in a silt loam soil. Their research included
the use of smoke injected into drain lines to observe its transmission to the soil surface. Smoke-emitting macropores werelocated within 50 cm of the drain line and the distance from thesubsurface drains correlated with the inﬁltration rate. The rate at
which water entered earthworm burrows declined with the logarithm of distance from the drain tile. Shipitalo and Gibbs (2005)
investigated the importance of this connectivity by observing thestructure of macropores and their surface connectivity (Fig. 1).
This “direct connectivity” phenomenon was veriﬁed by
Akay and Fox (2007) by conducting inﬁltration experiments
in a laboratory soil column with an artiﬁcial macropore placed
directly above or shifted away from a subsurface drain. The
experimental setup allowed the length of surface-connected or
buried macropores to be varied without unpacking or disturbing
the soil column between experiments. It was observed that thelonger the buried macropore (i.e., as the macropore approached
the soil surface), the more rapid the response at the drain outlet
and the higher the fraction (35–40%) of the total drain ﬂow
that came from the macropore. Arrival times of drain ﬂow forcolumns with surface-connected macropores (approximately
eight times earlier than for columns with only matrix ﬂow)
decreased signiﬁcantly compared with columns with buried
macropores (only about two times earlier than for columns with
only matrix ﬂow). For experiments with 55-cm-deep surfaceconnected macropores, the average ratio of steady-state ﬂow
rates through macropores and matrix decreased as the lateral
distance of the macropore from the drain increased: 2.4, 2.1,
and 1.6 for distances of 0, 6.25, and 12.5 cm, respectively.
An extrapolation of this trend suggests that only macropores
located within 20 to 25 cm of the drain can be considered tobe hydrologically directly connected.
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008

F . 1. A photograph of a natural, directly connected macropore
in the field preserved with polymer resin, demonstra ng the connec on between a macropore and a subsurface drain (courtesy of
Mar n Shipitalo, USDA-ARS, Coshocton, OH).

Steenhuis et al. (1997) modeled preferential transport tosubsurface drains using a conceptual model based on two linear
reservoirs located near the soil surface and closer to the drain.
Their results suggested that methods based on ﬁeld-averaged
preferential ﬂow characteristics are needed to model preferential
ﬂow to tile drains. Fox et al. (2004) modiﬁed a pesticide transport model, the Root Zone Water Quality Model (RZWQM),
to include a direct connectivity using an express fraction (EF),which routes a user-speciﬁed fraction (estimated to be 2%) ofwater and contaminants through macropore ﬂow directly fromthe soil surface to the drains. The estimated 2% EF was based on
the assumption that directly connected macropores were locatedwithin 20 to 25 cm of the drain (a 10-m distance between subsurface drain lines) at the modeled ﬁeld site in Allen County,
Indiana, with silty clay soils. The modiﬁed model more appropriately captured the immediate breakthrough of pesticides during
rainfall events shortly after pesticide application (Fox et al., 2004).
Fox et al. (2007) evaluated the EF modiﬁcation at an additional
ﬁeld site and suggested that EF must be calibrated to site-speciﬁc
conditions, and ranged between 2 and 5% for the simulated ﬁeld
conditions. It is hypothesized that the larger site-speciﬁc EF could
be a result of macropore–drain connectivity leading to shorter
contact times for reactive solutes.
Because of the relatively limited modeling of macropore–subsurface drainage interactions using process-based models, we
do not currently fully understand the water ﬂow dynamics ofdrain systems inﬂuenced by surface-connected or buried macropores. Most laboratory column studies are limited to seepage faceas opposed to subsurface drained bottom boundary conditions
and most ﬁeld studies are limited to measurements of breakthrough curves in tile drain ﬂow. Therefore, the objective of thisresearch was to improve our understanding of water ﬂow dynamics in soils containing macropores in the vicinity of subsurface
drains using a three-dimensional, ﬁnite-element, variably saturated ﬂow numerical model. The modeling was veriﬁed using910

the laboratory soil column experiments reported in Akay andFox (2007). Speciﬁc questions to be addressed in this research
included the following:• Can a numerical model that uses diﬀerent water retention andsoil hydraulic functions to describe soil hydraulic properties ofboth macropore and matrix domains properly simulate ﬂow
in drained systems with either surface connected or buried
macropores?
• Can the model conﬁrm the commonly observed relationshipbetween an increase in drain ﬂow (and decrease in ﬂow arrivaltime) with the macropore distance from the drain?
• Under what conditions do buried macropores become active
contributors to subsurface drain ﬂow?

Soil pressure heads were monitored every 10 s at 10, 40,
and 70 cm above the drain using four tensiometers (bubblingpressure = 100 cm H2O, Soil Measurement Systems, Tucson,
AZ) equipped with pressure transducers (ASDX001, Invensys,
Milpitas, CA) spaced along a horizontal plane at 12.5, 25.0,
31.5, and 38.0 cm from the left side of the column. The tensiometer–transducer system had a range of 5 to −70 cm H2O. A
datalogger (CR10X, Campbell Scientiﬁc, Logan, UT) continuously received and transmitted information to a computer for
automated pressure monitoring. Digital scales (EK-12Ki, A&D,
Milpitas, CA) recorded the outﬂow from the subsurface drain
for shifted and centered surface-connected macropore experiments. For centered, buried macropore experiments, a tube was
run from the top of the drain where the macropore entered
to a separate scale. The scales recorded outﬂow from both themacropore (discharge from the tube connected to the bottomof the macropore) and the soil matrix (collected from the subsurface drain tube) for centered buried macropore experiments(see Akay and Fox, 2007, Fig. 2).
For all experiments (i.e., centered and shifted, surface-connected and buried macropores), a Marriott-type inﬁltrometercreated a 1-cm ponded water layer at the soil surface for either
two or 3 h (referred to as the inﬁltration component of theexperiments). Data collection of pressure heads and drain ﬂow
continued for 24 h after the water supply was turned oﬀ (referred
to as the drainage component of the experiments).

Materials and Methods
Laboratory Column Studies

This research used data from the column studies of Akay andFox (2007), which were constructed to simulate macropores inthe vicinity of subsurface drainage. The 50-cm-wide by 28-cmdeep by 95-cm-long Plexiglas soil column consisted of a 6-cm
perforated polyvinyl chloride subsurface drain tube installed in
the center and approximately 2 cm above the impervious bottomof the soil column (Fig. 2). A sandy loam soil (55% sand, 32%
silt, and 13% clay) was packed into the column to a height of 75
cm at a bulk density of 1.6 g cm−3.
A 1-cm-diam. macropore was formed by wrapping a wooden
rod with aluminum and nylon meshes. Macropores were either
open at the soil surface (surface connected) and extending to different depths (15, 35, 55, and 75 cm) or fully buried, starting at
a certain depth (3, 15, 35, and 55 cm) and ending at the depthof the drain. Macropores were located either at the center of thesoil column above the drain or shifted 6.25 or 12.5 cm laterally.
Matrix ﬂow conditions were simulated with the wooden rod fully
inserted through the column down to the drain. The wooden
rod could be partially extracted to simulate various lengths ofa buried macropore. For surface-connected macropores, thesoil was packed to the lower part of the macropore and aroundthe remaining length of the wooden rod, which was later fully
extracted. Ponded water at the soil surface was allowed to entersurface-connected macropores but not buried macropores, which
consisted of a soil buﬀer between the soil surface and the beginning of the macropore.

Numerical Model Simula ons

The experiments were simulated using the three-dimensional
HYDRUS code. The HYDRUS program is a finite-element
model for simulating two- and three-dimensional movementof water, heat, and multiple solutes in variably saturated media.
The HYDRUS code numerically solves the Richards equation forsaturated–unsaturated water ﬂow (Šimůnek et al., 2006). Dueto sharp interfaces between the macropore and matrix domains,
a ﬁne ﬁnite-element grid around the macropore (Δxmin = 0.15cm), as well as small time steps (Δtmin = 10−10 h), were used
to maintain numerical stability. In contrast to earlier studies inwhich dual-porosity and dual-permeability one-dimensionalmodels with either ﬁrst-order (Castiglione et al., 2003; Šimůnek
et al., 2003) or second-order (Köhne and Mohanty, 2005)
mass transfer terms were used, this study involved a single
ﬂow (Richards) equation used for the entire three-dimensional transport domain, including both the matrix and amacropore. Earlier work by Köhne and Mohanty (2005)
demonstrated that, for a laboratory column packed with soils
having two distinctive soil hydraulic functions, utilization ofa pseudo-three-dimensional axisymmetrical Richards equation provided more accurate results than dual-permeability
models with ﬁrst- and second-order mass transfer terms. The
van Genuchten–Mualem model (van Genuchten, 1980) wasused to describe the water retention, θ(h), and conductivity,
K(h), functions for both matrix and macropore domains:

F . 2. The laboratory soil column and the corresponding three-dimensional
computa onal domain in HYDRUS. The macropore is located at the centerof the rectangular cross-sec on.
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008911⎧⎪
θs − θr⎪⎪ θ r +⎡1 + α h n ⎤ m
θ ( h ) = ⎪⎨⎣⎢⎦⎥⎪⎪⎪⎪ θ s⎩⎫h < 0 ⎪⎪⎪⎪⎬⎪⎪h ≥ 0 ⎪⎪⎭

[1]

m ⎤2
⎡
K ( h ) = K s S el ⎢1− (1− S e1/ m ) ⎥⎣⎦m = 1 −1 n , n > 1

though controlled drainage has been found to reduce NO3–N
and P losses by 35 to 45% (Skaggs et al., 1994), the response of
the preferential ﬂow paths is not yet understood.
For the hypothetical scenarios with buried macropores, theboundary condition at the bottom of the macropore was set
to zero pressure head to simulate its direct connection to thesubsurface drainage system (buried macropores are open to theatmosphere through the subsurface drain), whereas the seepage
face boundary condition was used for the remainder of the drain
in contact with the matrix domain. A no-ﬂux boundary condition was used at the soil surface. Simulations included buried
macropores with various lengths (i.e., 15, 35, 55, and 80 cm) andvarious initial water table levels (i.e., 15, 35, 55, and 80 cm above
the top of the drain). The column dimensions remained the sameas in simulations with surface-connected macropores.
The van Genuchten (1980) model parameters for the matrix
domain were obtained using inverse estimation and the laboratory column inﬁltration experiment without a macropore, i.e.,
involving only matrix ﬂow. Inverse estimation methods allow themost suitable model description of experimental data (Köhne andMohanty, 2005). During inverse estimation, an objective function that represents deviations between measured and calculated
space–time variables and soil hydraulic properties is minimized.
Minimization of this objective function can be accomplishedusing the Levenberg–Marquardt nonlinear minimization method(Šimůnek et al., 2006). The space–time variable chosen for theparameterization of the matrix ﬂow model (without a macropore)
was the cumulative subsurface drain ﬂow obtained from the laboratory experiment. The soil hydraulic parameters Ks, α, and nwere optimized, whereas θs and θr were ﬁxed and hence excluded
from the inverse estimation process. Initial estimates of the soil
hydraulic parameters for Ks, α, and n and the ﬁxed values of θsand θr were obtained from the water retention analysis conducted
on a soil core sample extracted from the laboratory column of
Akay and Fox (2007). Water retention was determined on undisturbed soil cores using Tempe cells for 0- to 1000-cm pressures.
For pressures >1000 cm, water retention was determined on disturbed soil samples by the pressure plate extractor method. The
RETC code (van Genuchten et al., 1991) was used to ﬁt the van
Genuchten parameters of θs, α, and n (Table 1) based on datafrom this undisturbed soil core (Akay and Fox, 2007).
The Ks value of the macropore domain was the only van
Genuchten (1980) parameter adjusted by a trial-and-error process
to match the cumulative drain ﬂow for the laboratory experimentwith the 55-cm, centered, surface-connected macropore. The
other parameters were held constant to reﬂect a small air-entry
pressure (high α value), below which the hydraulic conductivity
reduced to negligible values (high n value) (Castiglione et al.,2003). Using the calibrated value of the macropore-domain Ks,
various simulations were performed with centered and shifted,
surface-connected macropore lengths (i.e., 15, 25, 35, 45, 55, 65,
and 70 cm) and the ratios of total to matrix steady-state discharge
relative to macropore length were compared.
In contrast to surface-connected macropores, the van
Genuchten (1980) parameters of θs, α, and n were found to signiﬁcantly inﬂuence the exchange between the matrix and buried
macropores. Two sets of ﬁxed parameters (θs = 1.0, α = 1 cm−1,
n = 2 or 6) were used in the presented simulations. Two valuesof n (2 and 6) were used to show the sensitivity of the model[2]

where Se = (θ − θr)/(θs − θr) is the eﬀective saturation; α [L−1], n,and l are empirical parameters; θs is the saturated water content [L3L−3]; θr is the residual water content [L3 L−3]; and Ks [L T−1] is thesaturated hydraulic conductivity. A three-dimensional numerical
grid comprising 31 vertical layers was constructed (Fig. 2). The
number of finite-element nodes ranged between 25,460 and
53,692, depending on the length of the simulated macropore.
The boundary conditions for the numerical model werechosen according to the experimental conditions used by Akay
and Fox (2007). The constant pressure head (1-cm) boundary
condition was used at the top of the soil proﬁle during the ﬁrst 2
or 3 h to simulate inﬁltration. The boundary condition was thenswitched to a no-ﬂux boundary condition at the end of the inﬁltration period to simulate the redistribution and drainage parts of
the laboratory experiment. The seepage face boundary condition
was speciﬁed along a cylinder representing the 6-cm-diametersubsurface drain that was located 2 cm above the no-ﬂux baseboundary. A seepage face boundary acts as a zero pressure headboundary for nodes that are saturated and a no-ﬂux boundaryfor nodes that are unsaturated. When simulating the laboratory
experiments, the initial condition in the ﬂow domain was speciﬁed based on tensiometer readings obtained from the laboratory
study, with a pressure head of 2 cm at the bottom (i.e., maintaining a zero pressure head at the base of the subsurface drain) and−77 cm at the soil surface.
Observations were made regarding the conditions necessary
for the development of ﬂow in the buried macropore regime.
Buried macropores become active only when positive pressuresdevelop in their vicinity in the matrix domain. In simulating theinﬁltration experiments with buried macropores, the model was
not able to predict macropore ﬂow due to simulating the drain
as a seepage face boundary condition (i.e., no pressure buildup
was simulated in the matrix in the vicinity of the macropore).
The subsurface drain represented as a seepage face boundary inthe model was, in fact, a perforated tube with <10% open area.
In addition, ﬁne migrating particles may become lodged in thesesmall openings, clogging the soil–drain interface and increasing
the resistance to ﬂow. Therefore, the laboratory drain did not act
as a true seepage face boundary.
Since no detailed information was available on the degreeof restriction at the drain, a hypothetical drainage scenario wassimulated with buried macropores. Numerical simulations were
performed for a scenario in which the water table was initiallyabove the subsurface drain. Such cases can be observed in subsurface drain management applications. For example, controlled
drainage systems are being applied to subsurface-drained ﬁelds toreduce NO3 losses by placing a weir in the drainage outlet so as
to raise the water level in the outlet and reduce subsurface drainage rates (Skaggs et al., 1995). This type of management practice
can raise the average groundwater level to 20 to 25 cm below theground level (Wesstrom and Messing, 2007; Skaggs et al., 1995).
As a result of the high water table, the otherwise ineﬀective buried
macropores can transmit signiﬁcant ﬂow due to their relatively
high conductance during lowering of the groundwater level. Even
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008912

to this parameter, where n = 2 was selected to represent a value
for a typical soil and n = 6 was selected to obtain an immediate release of water from saturated to residual moisture content.
Theoretically, an inﬁnite n would be required to represent theimmediate decline in θ from θs to θr with increases in the soilpressure head. The percentage of drain ﬂow entering the drainfrom the macropore and the matrix was calculated.

Results and Discussion
Calibra on to Matrix Flow

While values of the calibrated parameters for the matrix
ﬂow regime are given in Table 1, water retention and hydraulicconductivity functions are shown in Fig. 3. The water retention
curves for the macropore domain, with either n = 2 or n = 6,
demonstrate a more immediate release of water with increased
pressure head (Fig. 3). The 95% conﬁdence intervals of optimized parameters are generally small compared with the inversely
estimated values (Table 1). The saturated hydraulic conductivity,
Ks, measured on the core sample taken from the laboratory soilcolumn was approximately four times higher than the calibrated
Ks, while the α parameter ﬁt by RETC to the experimentallymeasured retention curve was approximately 2.4 times smaller
than the optimized value.
A comparison between experimental and HYDRUS-simulated
cumulative drain ﬂow as a function of time demonstrates that thenumerical model was able to capture both the water front arrival
time and the cumulative drain ﬂux (Fig. 4). The Nash–Sutcliﬀe(NS) eﬃciency was slightly greater than 0.99, suggesting a nearperfect ﬁt between the simulated and observed values. The model
also predicted reasonably well pressures heads measured in thelaboratory using tensiometers (Fig. 5). Good agreement (i.e., NS= 0.97, 0.78, and 0.96 during inﬁltration and NS = 0.88, 0.84,
and 0.89 during drainage for tensiometers at 10, 40, and 70 cm
above the subsurface drain, respectively) between measured andpredicted pressure heads 10, 40, and 70 cm above the subsurface
drain demonstrates the ability of the numerical model to simulate
well the subsurface drainage boundary condition (Fig. 5).

F . 3. Soil hydraulic func ons represen ng the matrix and macropore domains: (a) soil moisture characteris c curve and (b) hydraulicconduc vity func on. Parameters for the matrix domain wereobtained using inverse es ma on. The n values for the macropore
domain were selected as follows: n = 2 to represent a value for atypical soil and n = 6 to obtain an immediate release of water fromsaturated to residual moisture content.

of the macropore; however, an increase in Ks of the macropore
domain (i.e., from 4.0 to 60,000 cm h−1), in an attempt to match
Surface-Connected Macroporesthe cumulative drain ﬂux under these conditions, suggested an
Flow conditions in the 55-cm-long, centered, surface-conupper bound for Ks. Since there is a 20-cm-thick matrix layer
nected macropore were only slightly aﬀected by the soil hydraulicbetween the bottom of the macropore and the drain (the top
parameters θr, θs, α, or n of the macropore domain. Flow condiof the subsurface drain is 75 cm below the surface), an additions were aﬀected mainly by the saturated hydraulic conductivity
tional increase in Ks could not generate a further increase incumulative drain ﬂow. The cumulative drain
ﬂux remained constant when higher Ks values
T
1. Soil hydraulic parameters† for the matrix and macropore domains using the van
Genuchten (1980) model. The 95% confidence intervals are given for parameters op were used, as the entire macropore became
mized with either HYDRUS or RETC.
almost instantaneously fully saturated. For
2example,Ks values of 100, 1000, 10,000, andDomain
n
K
l
R
θr
θs
α
s
60,000cm
h−1 produced simulated cumula3−3−1−1—m m —
cm
cm htive drain ﬂows of 7200, 10,200, 12,200, andCore sample‡
0.03§ 0.46 ± 0.01 0.005 ± 0.001 1.63 ± 0.24 6.22¶–
1.0
Model, matrix#
0.03§ 0.46§
0.012 ± 0.001 1.56 ± 0.08 1.46 ± 0.08 0.5 1.0
12,300 cm3, respectively, compared with theModel, macropore
0.00
0.46, 1.0†† 0.1
2, 6††
60,0000.5–
observed cumulative drain flow of 15,200† θ r, residual volumetric water content; θ s, saturated volumetric water content; α, n, and l, empiricm3. This cumulative drain ﬂow was simulated
cal parameters; Ks, saturated hydraulic conduc vity.
only when a 70-cm-long, centered, surface‡ Fi ed using RETC (van Genuchten et al., 1991).
connected
macropore was considered in the§ Fixed (not fi ed).
model.
¶ Determined on undisturbed soil cores by the constant-head method.
# Es mated using HYDRUS (Šimůnek et al., 2006).
Using the macropore Ks of 60,000 cm†† Used only for simula ons with buried macropores.
h−1, the model was able to capture the more
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008913

immediate response of laboratory tensiometers (Fig. 6) than when
the ﬂow system had no macropores (Fig. 5), although with lower
NS values (i.e., NS = 0.41, 0.98, and 0.64, respectively, for thetensiometers 10, 40, and 70 cm above the drain) than the matrix
ﬂow simulation. The reason for the lower NS value for the tensiometer 10 cm above the drain was the lack of observed data dueto the pressure exceeding the limit of the tensiometer. The model
did not predict as rapid a response as observed at the tensiometer
70 cm above the drain (Fig. 5a). In general, the model performed
better in matching the drainage part (i.e., 2–24 h) of the experiment (i.e., NS above 0.75 for all tensiometers) with the 55-cm,
centered, surface-connected macropore than its inﬁltration part(i.e., 0–2 h). Figure 7 presents simulated pressure head proﬁlesthat show ﬂow between the macropore and matrix on initiation of
the inﬁltration due to hydraulic nonequilibrium conditions. Thisﬁgure indicates the maximum pressure head of 55 cm observedat the bottom of the macropore. With increased time, the matrix
wetting front progressed down the column and the pressure heads
in the soil matrix surrounding the macropore increased due tolateral exchange between the macropore and matrix domains.
Simulations, as well as experimental data, showed a signiﬁcant inﬂuence of the length of the centered, surface-connected

F . 4. A comparison of measured and model-calibrated dimensionless flow volume, V* (where Vf = flow volume, φ = porosity, and Vc= column volume) vs. dimensionless me, t* (where t = me, Ks =saturated hydraulic conduc vity, and Lc = length of the column) forthe column experiment without a macropore (i.e., matrix flow).

F . 5. A comparison of measured (averaged) and calculated pressure heads 70, 40, and 10 cm above the subsurface drain during
the (a) infiltra on and (b) drainage components of the column
experiment without a macropore (i.e., matrix flow). The averaged
values represent the mean of four tensiometers along a horizontal
plane (12.5, 25.0, 31.5, and 38.0 cm from the le side) at the indicated height.
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008

F . 6. A comparison of measured (averaged) and calculated pressure heads 70, 40, and 10 cm above the subsurface drain during
the (a) infiltra on and (b) drainage components on the soil column
with the surface-connected, 55-cm-long macropore in the center ofthe profile. The averaged values represent the mean of four tensiometers along a horizontal plane (12.5, 25.0, 31.5, and 38.0 cm
from the le side) at the indicated height.
914

Butt (1999) could not verify soil homogeneityand constant macropore properties in the ﬁeld.
Also, the ﬁeld observations of Shipitalo and Butt(1999) were for soils that did not have subsurface
drains. An open burrow length in the vicinity ofthe drain appears to be an important macropore
characteristic for determining the potential eﬀect
on the drain ﬂux and potentially could be evenmore important for contaminant transport. In
fact, Allaire-Leung et al. (2000a,b) noted thatthe importance of macropore continuity seems
to increase with an increase in the adsorption
characteristics of the contaminant.
Under current domain dimensions andhydraulic conditions, the simulated arrival timeF . 7. Simulated pressure head profiles in the cross-sec on through the center of thesoil column for the infiltra on test with the surface-connected, 55-cm-long centered
for the soil column without macropores wasmacropore at (a) 0.017 h, (b) 0.17 h, and (c) 0.5 h.
1.37 h, while it was 0.33 h for the soil column
with the surface-connected, 55-cm-long, cenmacropores on the drain ﬂux. For example, 15-, 35-, 55-, andtered macropore that extended down approximately 75% of the
70-cm-long macropores increased the drain ﬂux approximately
distance to the drain. A sigmoidal relationship (R2 = 0.998) was
1.03, 1.22, 1.60, and 2.41, respectively, compared with the caseobtained between the length of surface-connected macropores andwithout macropores (Fig. 8). The eﬀect of the macropore length
the simulated outﬂow arrival times when arrival times of all surfaceon the drain ﬂux, however, diﬀers from the ﬁeld observations ofconnected macropores were considered (Fig. 9). This conﬁrms that
Shipitalo and Butt (1999), who observed no signiﬁcant relationshipthe bulk soil with a relatively slower water movement is largely
between the burrow length and inﬁltration rates. While homogebypassed when surface-connected macropores are introduced into
nous soil, as well as constant macropore characteristics, was achieved
the system.
under the laboratory conditions, the ﬁeld research of Shipitalo and
The eﬀect of surface-connected macropores on the increaseof the steady-state discharge decreases with the lateral distanceof the macropore from the drain, especially as the macropore
penetration percentage increases. Total to matrix steady-state discharge became more uniform across all shifted distances whenthe macropore penetration of the shifted macropores decreased(Fig. 8). Since the 75-cm-long, centered, surface-connected
macropore was directly connected also with the drain (i.e., pipe
ﬂow phenomena), the data point for this case is not included inFig. 8. In this case, the cumulative drain ﬂow becomes directly
proportional to Ks.
The laboratory data for the soil proﬁle with the 55-cm-long,surface-connected, centered macropore matched the results ofF . 8. Flux ra os vs. the lateral distance of the macropore fromthe drain. The flux ra o is the ra o between the steady-state flux
for scenarios with macropores and the one without a macropore.
Experimental data are for the soil column with the surfaceconnected, 55-cm-long macropore. Calcula ons were performed
for macropores of diﬀerent lengths. The flux ra o is 1 when there
is no eﬀect of the macropore on the steady-state flux. Error bars
represent one standard devia on calculated from triplicate lab
experiments. Lines represent HYDRUS simula on results.
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008

F . 9. Simulated reduc on in the flow arrival me as a funcon of the macropore penetra on for the soil column containing
surface-connected macropores compared with condi ons without
macropores.
915

numerical simulations with the 70-cm-long,surface-connected, centered macropore; the65-cm-long, surface-connected macropore
shifted 6.25 cm away from the drain; and the
55-cm-long, surface-connected macropore
shifted 12.5 cm away from the drain (Fig. 8).
Several potential reasons exist for the discrepancies between the laboratory experiments and thenumerical model. Because of the pressure heads
generated, the most likely reason is ﬁngering or
additional preferential ﬂow created in the laboratory experiments in the soil layer below thebottom of the macropore and the drain. Thisadditional preferential ﬂow had a larger impact F . 10. Simulated pressure head profiles in the cross-sec on through the center of theon steady-state ﬂux when the macropore was soil column for the drainage test with the buried 55-cm-long centered macropore andcentered over the drain than when it was shifted the ini al water table level at 80 cm at (a) 0.017 h, (b) 0.17 h, and (c) 0.5 h.
away from the drain.
In general, the numerical simulations suggested an even
the water table was initially above the subsurface drain, demongreater contributing area with directly connected macropores
strated that the percentage of macropore ﬂow compared with thethan suggested by Akay and Fox (2007) and Fox et al. (2004). The
total drainage outﬂow (through both the matrix and macropore
need for calibration observed by Fox et al. (2007) is supported bydomains) ranged from 32 to 48% for n = 2 and from 19 to 43%
these numerical simulations that show that the contributing areafor n = 6 (Table 2). For their imperfect drain condition wheredepends on the site-speciﬁc macropore depth of penetration.
the water table built up over the drain, Akay and Fox (2007)
observed that approximately 35% of the total ﬂow came directly
Buried Macropores
from 20- and 40-cm-long buried macropores and approximately
Since buried macropores may be connected to the drain but
40% from 60- and 72-cm-long buried macropores. Distinctive
not to the soil surface, the hydraulic conditions in the domain
responses of the matrix and buried macropore domains to drainchange drastically from conditions with surface-connected
age conditions are plotted in terms of cumulative outﬂow curves
macropores (Fig. 10). The early hydraulic nonequilibriumin Fig. 11. Although outﬂow from both domains started simultaobserved during simulations with surface-connected macroponeously, initially steep lines representing outﬂow from macropores
res does not exist in simulations
2. The part of the drain ou low arriving directly through the buried macropore (n values ofwith buried macropores. The T
seepage face boundary condition the van Genuchten [1980] model represent the macropore domain) for ini al water table levels of+15, +35, +55, and +80 cm.
could not be used for the section
Drain ou low from buried macropore
Buriedof the boundary where the buried
macropore+15 cm+35 cm+55 cm+80 cmmacropore enters the subsurface
length
n=2
n=6
n=2
n=6
n=2
n=6
n=2
n=6
drain since the model calculates————————————————————————— % —————————————————————————
here negative pressure heads, which 15 cm
38.93
31.07
33.15
22.02
31.86
19.58
31.68
18.75
39.70
31.64
41.38
33.62
39.06
29.47
38.55
28.29
leads to no-ﬂux conditions. For this 35 cm
39.72
31.81
41.71
34.04
44.26
37.17
43.63
35.72
reason, the seepage face boundary 55 cm
80 cm
39.81
31.97
41.80
34.24
44.51
37.61
48.30
42.91
condition was assigned aroundthe drain, except for the buried
macropore section. This part of theF . 11. Simulated
drain was assigned a zero pressure
cumula ve macropore
head boundary condition since thisand matrix ou lows for
section was connected to the drain
the drainage experioutlet. It should be noted that inment with the buried
55-cm-long centeredtilled soils, buried macropores initimacropore. Ini al waterate at the interface between the plowtable was 80 cm above
layer and the soil matrix. These two
the drain. The n values
layers may possess unique hydraufor the macropore
lic properties. Therefore, a layered
domain were selected
as follows: n = 2 toagricultural ﬁeld soil with buried
represent a value for a
macropores could behave differtypical soil and n = 6 to
ently from the artiﬁcial soil system
obtain an immediate
under study.
release of water from
Simulations reﬂecting the pracsaturated to residual
moisture content.
tical conditions described above forcontrolled drainage practices, where
www.vadosezonejournal.org · Vol. 7, No. 3, August 2008916

suggested an immediate release of water from buried macropores.
Even though this volume was higher than the matrix cumulative
outﬂow during the ﬁrst couple of hours, the overall contribution of the matrix reached about 60% of the total outﬂow in thelong term. It should be noted here that these ratios could change
depending on the hydraulic parameters and domain properties(i.e., a macropore radius and domain dimensions).
When water tables are initially above the subsurface
drain, two variables are expected to signiﬁ cantly aﬀ ect thewater exchange between the matrix and macropore domains:
the pressure gradient and the interfacial surface area between
the macropore and matrix domains. Th e water table position determines the pressure gradient, whereas the length of
the buried macropore determines the interfacial surface area.
When the initial water table level is at or below the upper endof the buried macropore, the percentage of macropore ﬂow
is approximately the same and independent of the macropore
length. The macropore ﬂow contribution increases proportionally to the macropore length, however, when the initial watertable is greater than the macropore length due to the increasein the interfacial surface area. For example, when the initialwater level is 55 cm (n = 2), the percentage of macropore
ﬂow increases from 32 to 44% when the macropore length
increases from 15 to 55 cm, and then remains approximately
constant for longer (e.g., 80-cm) buried macropores (Table 2).
For a given length of the buried macropore, the percentage
of macropore ﬂow maximizes when the water table is at anequivalent elevation to the macropore length. When the watertable elevation further increases, the percentage of macropore
ﬂow decreases relative to the total matrix plus macropore ﬂow.
For example, for the buried 55-cm-long macropore (n = 2),
the percentage of macropore ﬂow increased from 40 to 44%
when the initial water table was increased from 15 to 55 cm,
respectively, and then decreased approximately 1% when theinitial water table was at 80 cm (Table 2). This is due to thefact that the increase in macropore ﬂow (5.8 cm3) is smaller
than the corresponding increase in matrix ﬂow (11.1 cm3), as
shown in Table 3.

dual-permeability models to describe such conditions, thenumerical solution of the three-dimensional Richards equation with two distinctive hydraulic functions for the matrix
and macropore domains proved to be useful as well. Presented
numerical simulations veriﬁed the model’s ability to handle
sharp interfaces between the two domains, with dramaticallydiﬀerent soil hydraulic parameters. Simulations also conﬁrmed
the frequently observed relationships between drain ﬂow and
its arrival time with the macropore distance from the subsurface
drain. The importance of the macropore’s penetration distanceon this relationship is highlighted. Reduction in ﬂow arrival
times for conditions with macropores compared with matrix
ﬂow conditions was directly related to the penetration depth.
Numerical simulations suggested a greater contributing area
(i.e., the region where the macropore has an eﬀect on the totaldrain ﬂow) with directly surface-connected macropores thansuggested in the literature, with this area again depending on
the macropore length. It should be researched further whether
the contributing area should vary for water, tracers, and reactive
contaminants. As the buried macropores could not be activatedby surface inﬁltration of a limited duration, an alternative scenario with the water table initially above the subsurface drain
was also simulated. Considering an immediate release of waterfrom the drain-connected buried macropores and their further
ability to transmit 30 to 40% of ﬂow, the chemical compositionof water in macropores can dictate the water quality in the outﬂow of the subsurface drains during early stages of drainage.
A
This project has been supported by the USDA Cooperative State
Research, Education, and Extension Service (CSREES) NRI Seed
Grant (2004-35102-14890) and USDA CSREES NRI Grant (200735102-18242).

References
Akay, O., and G.A. Fox. 2007. Experimental investigation of direct interconnectivity between macropores and subsurface drains during inﬁltration. Soil
Sci. Soc. Am. J. 71:1600–1606.
Allaire-Leung, S.E., S.C. Gupta, and J.F. Moncrief. 2000a. Water and solute
movement in soil as inﬂuenced by macropore characteristics: 1. Macropore
continuity. J. Contam. Hydrol. 41:283–301.
Allaire-Leung, S.E., S.C. Gupta, and J.F. Moncrief. 2000b. Water and solute
movement in soil as inﬂuenced by macropore characteristics: 2. Macropore
tortuosity. J. Contam. Hydrol. 41:303–315.
Buttle, J.M., and D.G. Leigh. 1997. The inﬂuence of artiﬁcial macropores on water and solute transport in laboratory soil columns. J. Hydrol.
191:290–314.
Castiglione, P., B.P. Mohanty, P.J. Shouse, J. Šimůnek, M.Th. van Genuchten,
and A. Santini. 2003. Lateral water diﬀusion in an artiﬁcial macroporous
system: Modeling and experimental evidence. Vadose Zone J. 2:212–231.
Dorner, S.M., W.B. Anderson, R.M. Slawson, N. Kouwen, and P.M. Huck.
2006. Hydrologic modeling of pathogen fate and transport. Environ. Sci.

Summary and Conclusions
Several inﬁltration experiments previously conducted by
Akay and Fox (2007) were simulated using HYDRUS. Model
results demonstrated significant hydraulic nonequilibrium
conditions between the matrix and macropore domains duringinﬁltration into a soil column that contained surface-connected
macropores. The hydraulic nonequilibrium was less pronounced
for simulations involving buried macropores. In addition topreviously used one- and two-dimensional dual-porosity andT
3. Components of the total drain ou low coming from the matrix and the buried 55-cm-long macropore at the early ( me t = 10 s)
and final (t = 24 h) stage of drainage (n values of the van Genuchten [1980] model represent the macropore domain).
Ini al watertable level
cm
+15
+35
+55
+80

Cumula ve flow entering the drain
Total drainage amountFrom the buried macropore
From the matrix
10 s
24 h
10 s
24 h
10 s
24 h
n=2
n=6
n=2
n=6
n=2
n=6
n=2
n=6
n=2
n=6
n=2
n=6
—————————————————————————————————— cm3 ——————————————————————————————————
13.7
20.1
34.0
27.4
2.3
2.6
51.7
58.7
16.0
22.8
85.7
86.1
27.7
43.6
75.6
62.2
4.4
5.5
105.7
120.5
32.1
49.1
181.4
182.7
39.5
64.0
110.5
95.0
6.1
7.9
139.2
160.5
45.6
71.9
249.7
255.5
6.1
7.9
150.3
175.2
45.5
71.9
266.6
272.5
39.4
64.0
116.3
97.3

www.vadosezonejournal.org · Vol. 7, No. 3, August 2008917

Technol. 40:4746–4753.
Fox, G.A., R. Malone, G.J. Sabbagh, and K. Rojas. 2004. Interrelationship ofmacropores and subsurface drainage for conservative tracer and pesticidetransport. J. Environ. Qual. 33:2281–2289.
Fox, G.A., G.J. Sabbagh, R.W. Malone, and K. Rojas. 2007. Modeling parent
and metabolite fate and transport in subsurface drained ﬁelds with directly
connected macropores. J. Am. Water Resour. Assoc. 43:1359–1372.
Geohring, L.D., P.E. Wright, T.S. Steenhuis, and M.F. Walter. 1999. Fecal coliforms in tile drainage eﬄuent. ASAE Pap. 992203. Am. Soc. Agric. Eng.,
St. Joseph, MI.
Ghodrati, M., M. Chendorain, and Y.J. Chang. 1999. Characterization of macropore ﬂow mechanisms in soil by means of a split macropore column.
Soil Sci. Soc. Am. J. 63:1093–1101.
Hoorman, J.J., J.N. Rausch, T.M. Harrigan, W.G. Bickert, M.J. Shipitalo, M.J.
Monnin, S.R. Reamer, F.E. Gibbs, M.I. Gangwar, H. Keener, and L.C.
Brown. 2005. Research, educational, and technical assistance priorities forliquid manure application in the Midwest. Pap. 052062. In ASAE Annu.
Int. Meet., Tampa, FL. 17–20 July 2005. Am. Soc. Agric. Eng., St. Joseph,
MI.
Jamieson, R.C., R.J. Gordon, K.E. Sharples, G.W. Stratton, and A. Madani.
2002. Movement and persistence of fecal bacteria in agricultural soils andsubsurface drainage water: A review. Can. Biosyst. Eng. 44: 1.1–1.9.
Joschko, M., H. Diestel, and O. Larink. 1989. Assessment of earthworm burrowing eﬃciency in compacted soil with a combination of morphologicaland soil physical measurements. Biol. Fertil. Soils 8:191–196.
Joy, D.M., H. Lee, C.M. Reaume, H.R. Whitely, and S. Zelin. 1998. Microbial
contamination of subsurface tile drainage water from ﬁeld applications ofliquid manure. Can. Agric. Eng. 40:153–160.
Kladivko, E.J., J. Grochulska, R.F. Turco, G.E. Van Scoyoc, and J.D. Eigel. 1999.
Pesticide and nitrate transport into subsurface tile drains of diﬀerent spacings. J. Environ. Qual. 28:997–1004.
Köhne, J.M., and B.P. Mohanty. 2005. Water ﬂow processes in a soil column
with a cylindrical macropore: Experiment and hierarchical modeling. Water Resour. Res. 41:1–17.
Li, Y., and M. Ghodrati. 1997. Preferential transport of solute through soilcolumns containing constructed macropores. Soil Sci. Soc. Am. J.
61:1308–1317.
Magesan, G.N., R.E. White, and D.R. Scotter. 1995. Inﬂuence of ﬂow rate on
the concentration of indigenous and applied solutes in mole-pipe drain
eﬄuent. J. Hydrol. 172:23–30.
Phillips, R.E., V.L. Quisenberry, J.M. Zeleznik, and G.H. Dunn. 1989. Mechanism of water entry into simulated macropores. Soil Sci. Soc. Am. J.
53:1629–1635.
Shipitalo, M.J., and K.R. Butt. 1999. Occupancy and geometrical properties of Lumbricus terrestris L. burrows aﬀecting inﬁltration. Pedobiologia
43:782–794.
Shipitalo, M.J., and W.M. Edwards. 1996. Eﬀects of initial water content onmacropore/matrix ﬂow and transport of surface-applied chemicals. J. Environ. Qual. 25:662–670.
Shipitalo, M.J., W.M. Edwards, W.A. Dick, and L.B. Owens. 1990. Initial storm
eﬀects on macropore transport of surface-applied chemicals in no-till soil.
Soil Sci. Soc. Am. J. 54:1530–1536.
Shipitalo, M.J., W.M. Edwards, and C.E. Redmond. 1994. Comparison of water movement and quality in earthworm burrows and pan lysimeters. J.
Environ. Qual. 23:1345–1351.
Shipitalo, M.J., and F. Gibbs. 2000. Potential of earthworm burrows to transmit
injected animal wastes to tile drains. Soil Sci. Soc. Am. J. 64:2103–2109.
Shipitalo, M.J., and F. Gibbs. 2005. Preferential ﬂow of liquid manure in macropores and cracks. Pap. 052063. In Proc. Annu. Conf., Am. Soc. Agric.
Biol. Eng., St. Joseph, MI.
Shipitalo, M.J., V. Nuutinen, and K.R. Butt. 2004. Interaction of earthworm
burrows and cracks in a clayey, subsurface-drained soil. Appl. Soil Ecol.
26:209–217.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review
and comparison of models for describing non-equilibrium and preferential
ﬂow and transport in the vadose zone. J. Hydrol. 272:14–35.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2006. The HYDRUS software package for simulating two- and three-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Tech. Manual.
Version 1.0. PC Progress, Prague, Czech Republic.

www.vadosezonejournal.org · Vol. 7, No. 3, August 2008

Skaggs, R.W., M.A. Breve, and J.W. Gilliam. 1994. Hydrologic and water
quality impacts of agricultural drainage. Crit. Rev. Environ. Sci. Technol.
24:1–32.
Skaggs, R.W., M.A. Breve, and J.W. Gilliam. 1995. Predicting eﬀects of watertable management on loss of nitrogen from poorly drained soils. Eur. J.
Agron. 4:441–451.
Steenhuis, T.S., M. Bodnar, L.D. Geohring, S.A.E. Aburime, and R. Wallach.
1997. A simple model for predicting solute concentration in agricultural
tile lines shortly after application. Hydrol. Earth Syst. Sci. 4:823–833.
Trojan, M.D., and D.R. Linden. 1992. Microrelief and rainfall eﬀects on water and solute movement in earthworm burrows. Soil Sci. Soc. Am. J.
56:727–733.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC code forquantifying the hydraulic functions of unsaturated soils. Version 1.0. EPA
Rep. 600/2-91/065. U.S. Salinity Lab., Riverside, CA.
Wesstrom, I., and I. Messing. 2007. Eﬀects of controlled drainage on N and P
losses and N dynamics in a loamy sand with spring crops. Agric. Water
Manage. 87:229–240.

918