Published online June 21, 2006

Comparison of Pesticide Transport Processes in Three Tile-Drained Field Soils
Using HYDRUS-2D
Arnaud Boivin,* Jirka Šimůnek, Michel Schiavon, and Martinus Th. van Genuchten
performed using repacked soil columns in the laboratory(Davidson and Chang, 1972; van Genuchten et al., 1977;
Gamerdinger et al., 1990; Kookana et al., 1993; Guo andWagenet, 1999). While useful, such experiments are
generally not representative of field conditions andoften underestimate the leaching potential of chemicals
in the field because of preferential flow or other processes (Flury, 1996).
Field experiments provide a more realistic assessmentof pesticide transport processes. Tile-drained agricultural fields in particular have been shown to be attractive for studying field-scale flow and solute transport
processes. Many studies have demonstrated the merit ofsuch systems for studying the transport of tracers(Vanderborght et al., 2002; Fox et al., 2004; Gerke andKöhne, 2004), N (Mohanty et al., 1998; de Vos et al.,
2000; de Vos, 2001), and pesticides (Abbaspour et al.,
2001; Kohler et al., 2001; Zehe and Flühler, 2001; Jaynes
et al., 2001; Novak et al., 2003; Larsbo and Jarvis, 2005).
Tile drainage systems, in essence, function as large undisturbed field-scale lysimeters (Gerke and Köhne,
2004; Kung et al., 2000) since they integrate the sitespecific flow and transport processes within the fulltransport domain that contributes to drainage outflow.
Optimally interpreting tile drainage response (e.g.,drain discharge rates and solute concentrations) is still achallenge because of field-scale spatial and temporalvariability and the related problems of preferential flow
(Southwick et al., 1995; van Genuchten et al., 1999;
Larsson et al., 1999; Hendrickx and Flury, 2001; Kohler
et al., 2003). During the past two decades, several models have been developed to simulate water flow andsolute transport in the vadose zone. These models range
from relatively simplistic approaches to more complex,
physically based, dual-porosity, dual-permeability, andmultiregion models (Jarvis, 1994; Banton et al., 1995;
Abbaspour et al., 2001; Kohler et al., 2003; Pruess, 2004;
Larsbo et al., 2005). A review of various approaches tosimulate preferential, nonequilibrium flow and transport
was recently given by Šimůnek et al. (2003). Although
the more sophisticated of these models seem to be able
to describe the outflow patterns from tile-drained fields,
this may be partly due to the large number of parametersthat are often fitted to the observed field data (e.g., Haws
et al., 2005). The resulting parameter estimation process
by itself can be quite problematic because of problems ofparameter uniqueness, identifiability, and stability (e.g.,
Hopmans and Šimůnek, 1997).
One challenge when simulating the outflow response
from tile drains at a given site is that different conceptual
formulations sometimes will lead to equally acceptable

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

ABSTRACTThe purpose of this study was to assess the transport of the herbicide bentazone [3-(1-methylethyl)-1H-2,1,3-benzothiadiazin-4(3H)one 2,2-dioxide] in three contrasting tile-drained cultivated fieldsoils subject to otherwise similar experimental conditions. Observed
drain discharge rates and chemical concentrations in the drainage
water reflected the different transport processes at the three sites
in the same area of northeastern France: a sandy loam site (Villey),a silt loam site (Bouzule-1), and a silty clay site (Bouzule-2). The
sandy loam site showed very little tile drainage (240 m3 ha21) during
the 100-d study in the spring of 2002, as well as low chemical losses
in the drainage water (0.16% [v/v] of the applied amount). While
little drainage was observed also for the silty clay soil (175 m3 ha21),
observed pesticide losses were considerably larger (1.25% of the applied amount). The silt loam soil, in comparison, showed much moredrainage (521 m3 ha21) and the highest chemical loads in the drainage water (2.7% of the applied amount). Numerical simulations ofdrain discharge with the HYDRUS-2D variably saturated flow andsolute transport model compared well with the observed data forthe relatively homogeneous sandy loam (Villey) and the silt loam
(Bouzule-1) soils. The saturated hydraulic conductivity of the bottom
layer in both cases was key to correctly predicting the drainage fluxes.
Accurate predictions of the silty clay field data (Bouzule-2) could
be obtained only when the soil hydraulic functions were modified toaccount for preferential flow through drying cracks near the soil
surface. Chemical concentrations could be better described using adual-porosity (mobile–immobile water type) transport model for all
three soils, including the sandy loam. Results indicate that water andpesticide transport in soils is governed by site-specific processes. Optimal use of the HYDRUS-2D flow and transport model allowed areasonable description of the field-scale pesticide processes using only
a limited number of adjustable parameters.

A

PPLICATION of pesticides to cultivated fields often
leads to their unintentional leaching through thevadose zone toward underlying and down-gradient water resources. Since excessive leaching compromises soil
and water quality, and indirectly human health (Garcia,
2003; Wong et al., 2003; Miersma et al., 2003), manystudies have been performed to assess the fate andtransport of a variety of pesticides in different soils andfor different soil management practices. Experiments focusing on the governing transport processes are often

Boivin and M.Th. van Genuchten, USDA-ARS, George E. Brown,Jr., Salinity Lab., 450 West Big Springs Rd., Riverside, CA 92507.
Jirka Šimůnek, Dep. of Environmental Sciences, Univ. of California,
Riverside, CA 92521. A. Boivin and M. Schiavon, Laboratoire Sols etEnvironnement, UMR INPL-ENSAIA/INRA, 2 Avenue de la Forêt deHaye, B.P. 172, 54505 Vandoeuvre-lès-Nancy Cedex, France. Received
22 July 2005. *Corresponding author (arnaud.boivin@lyon.cemagref.fr).
Published in Vadose Zone Journal 5:838–849 (2006).
Original Researchdoi:10.2136/vzj2005.0089ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: ADE, advection–dispersion equation; MIM, mobile–immobile solute transfer model.

838

839

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

www.vadosezonejournal.org

predictions of observed field data (Blazkova et al.,2002). Problems related to identification of the actual
physical processes at a site often become more apparent
when solute transport is simulated simultaneously withwater flow (e.g., Abbaspour et al., 2001; de Vos et al.,1999, 2002). A good example is given by Pang et al.
(2000), who studied the field-scale transport of Br andthree pesticides using the HYDRUS-2D variably
saturated flow and transport model of Šimůnek et al.
(1999). They found that the model was able to accurately describe soil water contents, with only relatively
minor discrepancies between simulated and observed
data; however, predicted concentrations of Br andpesticide showed major differences with the observed
data. The HYDRUS-2D model was unable to simulate
preferential flow through drying cracks near the soil
surface, which caused higher Br and pesticide concentrations at greater depths. Results of the field study byPang et al. (2000) suggest the need for a more detailed
dual-permeability approach for simulating preferential
flow. Unfortunately, the many parameters needed forsuch detailed models are generally not available, unlessfitted to comprehensive field data sets, with resultant
problems of parameter uniqueness and parameter equifinality (Beven, 1993) in that different conceptualizations and multiple combinations of input parameters
may lead to equally acceptable simulations of the fielddata. The importance of estimating the additional parameters required for the MACRO dual-permeability
model, and the significance of data availability in general for simulating pesticide transport in a tile-drained
structured field soil, was prominently demonstrated in a
recent study by Larsbo and Jarvis (2005).
This study concerns the transport of the herbicide
bentazone, sprayed in the spring of 2002 onto threedifferent tile-drained agricultural fields: a sandy loam, asilt loam, and a silty clay soil. The soils were selected
because of their contrasting flow and transport properties, while most other experimental conditions were very
similar (i.e., weather conditions, pesticide applications,and soil management). Also, several data sets were already available for the sites (Novak et al., 2001, 2003),
thus providing additional information helpful for running numerical models and comparing observed and simulated data.
The use of bentazone is a major issue since this pesticide is relatively mobile and frequently detected insurface waters and groundwater (Halfon et al., 1996;
Lagana et al., 2002; Institut Français de l’Environnement, 2003, 2004). Previous laboratory experiments
conducted with the same soils showed that bentazone
degradation and sorption are relatively minor compared with other pesticides (Boivin et al., 2004, 2005).
Carrizosa et al. (2000), for this same reason, used bentazone as a model chemical for very mobile herbicides.
Observed tile drain discharge and chemical concentrations in this study were simulated with HYDRUS-2D.
One major challenge in using HYDRUS-2D, or any
other model, is that for most field-scale applications,only limited information is available about the key parameters and processes involved. In this study, we

assessed the utility of HYDRUS-2D to simulate draindischarge and solute concentrations with limited input
data and for different field settings. Our aim was to compare observed and simulated water fluxes and soluteconcentrations using readily available information.
Different options in HYDRUS-2D were used in attempts to better simulate the preferential flow of waterand solutes. These options included the use of modified
soil hydraulic properties (Vogel and Cı́slerová, 1988) toaccount for increased flow near saturation, and the useof a dual-porosity (mobile–immobile water) type model
for solute transport (Šimůnek et al., 2003).
MATERIALS AND METHODSThe Three Field Sites
Three field sites in northeastern France (Lorraine) wereselected for the experiments, one site in Villey-St-Etienne(designated as Villey), and two sites in La Bouzule (designatedas Bouzule-1 and Bouzule-2), approximately 25 km from theVilley site. Winter wheat (Tritcum aestivum L.) was grown on
the three sites using the same standard management practices
in the area. Because of their close proximity, the three sites
were also subject to approximately the same weather conditions. Figure 1 shows a layout of the tile drainage systems at thethree sites. Cross-sections of the layered profiles, including thelocations of the drains, are shown in Fig. 2. The Villey site
(Table 1) consisted of sandy loam derived from alluvial deposits (Moselle River) on a marl layer (a very soft, calcareousrock with a chalky clay texture usually formed from a mixtureof clay and limestone). The PVC tile drains (0.05-m diameter)at this site were located at an average depth of 0.9 m, with aspacing of 12 m (Fig. 1 and 2).
The Bouzule-1 and Bouzule-2 sites were located on adjacent
fields. The 2.83-ha Bouzule-1 site consisted of a mediumtextured silt loam (Table 1), while the 1.85-ha Bouzule-2 site
contained a finer textured silty clay. The latter site had a 5 toVilleySandy loam soil
(2.5 ha)Diversion tube
Auto-sampler

Bouzule-1
Silt loam soil
(2.83 ha)Tile drain system(Ø = 0.05 m)

Bouzule-2
Silty clay soil
(1.85 ha)Collector(Ø = 0.1 m)

Auto-samplers
Fig. 1. Schematics of the layouts of the tile drainage systems at theVilley, Bouzule-1, and Bouzule-2 experimental sites.

840

VADOSE ZONE J., VOL. 5, AUGUST 2006

Precipitation and potential evapotranspiration rates andother meteorological variables were obtained from the MétéoFrance weather station in Toul (www.meteofrance.fr), about
8 km from the Villey site, and the INRA weather station inChampenoux, about 2 km from the two Bouzule sites (Fig. 3).
Measured precipitation during February 2002 in Toul was
143.0 mm and in Champenoux 123.5 mm, more than twice themonthly average (52 mm) for February in this area of France
during the previous three decades (1967–1997). The high rainfall rates caused the silty clay and silt loam at the two Bouzulesites to be very wet at that time. Similar conditions occurred at
the Villey site. Compared with February, far less precipitation
was recorded in March 2002 (a total of 70.2 mm at Villey and49.0 mm at the Bouzule sites), with most rainfall (94% at the
Bouzule sites and 93% at Villey site) occurring during theweek immediately after application of the bentazone. Little
or no precipitation was recorded thereafter until the end of thespring season.

Bentazone Application

25
20

10
5

Layers

30

Loam

Sand

Bulk density
23

C content%

g cm
Villey (sandy loam)%

Bouzule

11
6.2
6.2
54.2

19
3.7
3.7
40.5

70
1.52†
90.1
1.52†
90.1
1.50‡
5.3
1.40‡
Bouzule-1 (Silt loam)

0.87
0.46
-§
–

5.8
6.5
–
–

0–25
25–40
40–60
60–80
80–100

26
31.2
37.2
40.2
39.8

52
52.1
48.8
49
52.6

22
1.28¶
16.7
1.37¶
14
1.55¶
10.8
1.53¶
7.6
1.58¶
Bouzule-2 (Silty clay)

1.36
1.36
0.41
0.41
0.32

6.6
6.6
5.8
5.8
5.2

0–30
30–60
60–100

54.2
54.1
53

40.5
42.8
43

5.3
3.1
4

1.9
0.88
0.51

8.2
8.2
8.3

† from Bartoli et al. (2006).
‡ Estimated.
§ Not determined.
¶ From Waleed (1983).

2/1 2/12 2/23 3/6 3/17 3/28 4/8 4/19 4/30

pH in water

0–30
30–60
60–200
.200

1.17¶
1.30¶
1.40¶

Bentazone
application

15

0

Clay

20
18
16
14
12
10
8
6
4
2
0

VilleyTable 1. Selected physical and chemical properties of the VilleySt-Etienne (Villey) and La Bouzule (Bouzule-1 and Bouzule-2)
tile-drained field sites in the Lorraine area of northeasternFrance.
cm

20
18
16
14
12
10
8
6
4
2
0

Temperature (°C)

6% slope toward the Bouzule-1 silt loam site, which had aslope of only 1 to 2% (Fig. 1). The two Bouzule sites were
equipped with individual subsurface PVC tile drain systems(0.05-m diameter), with tile drains at an average depth of 0.9 m
and having spacings of 12 m in the silt loam (Bouzule-1) and
8 m in the silty clay (Bouzule-2).

30
Precipitation (mm)Fig. 2. Distribution of different soil horizons in the transport domain
representing half the drain spacings at the Villey sandy loam site,
the Bouzule-1 silt loam site, and the Bouzule-2 silty clay site.

Temperature (°C)

When wheat at the sites reached the second- or third-leaf
stages, the three fields each received an application of bentazone, equivalent to a dose of 1.2 kg ha21. The spray application
occurred on 4 Mar. 2002 at the Villey site but not until 11
March at the two Bouzule sites to allow the soils there to dry
out enough for the tractor to be able to go onto the fields. No
drainage was observed at any of the three sites between thePrecipitation (mm)Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Meteorological Conditions

25
20

Bentazone
application

15
10
5
0
2/1 2/12 2/23 3/6 3/17 3/28 4/8 4/19 4/30

Fig. 3. Daily precipitation rates (solid bars) and temperatures (lines)
recorded by the Météo France meteorological station in Toul(adjacent to the Villey site) and the INRA meteorological station atChampenoux (adjacent to the two Bouzule sites) during the experiments. The two vertical arrows indicate the times of bentazone
application.

841

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

www.vadosezonejournal.org

two applications since little or no precipitation was recorded
during this period. The herbicide was applied using the commercial formulation Bazamais (BASF Agro, Genay, France)by spraying onto the fields as a standard spring treatment by
farmers using their own equipment (a 2 by 6 m tractor draw
ramp sprayer). Bentazone is a selective contact herbicide
locally used on many different crops. The chemical is weakly
acidic (pKa 3.3), has a relative molecular mass of 240.3 g mol21
and a solubility of 570 mg L21 (in water at pH 7, 228C).

Bentazone Sorption and Degradation
Previous laboratory experiments using surface samples
from the three sites showed very low mineralization rates forbentazone. Boivin et al. (2004) observed that, after 7 d of incubation, the mineralized amounts (as a percentage of the applied dose) ranged from 0.11 to only 0.3% for the three soils.
Sorption was also found to be very minor, with 97 to 98% of
the applied dose being available immediately after applicationof bentazone to the different soil samples. After 7 d of incubation, 91 to 96% of the applied herbicide was still available.
Boivin et al. (2005) found that bentazone adsorption was notcorrelated with the soil organic matter contents of 13 soils,
which included soil samples from the three sites in this study.
Bentazone adsorption data could be described using linear
isotherms with distribution coefficients (kD) of 1.4 6 0.1, 1.4 6
0.1, and 1.2 6 0.1 L kg21 for the Villey, Bouzule-1, andBouzule-2 sites, respectively.

Tile Drain Sampling
The sampling system at each site consisted of a Sigma 900PMAX automatic sampler (Sigma-Aldrich, St. Louis, MO)
coupled to a container with a V-notch weir. This system madeit possible to link drainage-water herbicide sampling with thewater flux recordings (taken at 10-min intervals). Water samples were collected with the automatic sampler proportionallyto the measured flow rate and stored into 10-L glass bottles.
Bentazone in the tile drainage water was analyzed using
high performance liquid chromatography coupled with a C18
Merck Lichro CART 250–4 (250 mm long, 4 mm i.d.) column
and a diode-array detector (Varian 9065, Varian BV, Middelburg, the Netherlands). Concentrations were determined bytwice extracting the drain-water samples with dichloromethane (5:1 v/v water/dichloromethane). Since bentazone isweakly acidic, drainage water samples were acidified to pH 2
by adding HCl (0.1 M). The dichloromethane phase was recovered and evaporated to dryness under reduced pressureusing a Heidolph 94200 evaporator (Heidolph Instruments,
Schwabach, Germany) and then redissolved in high-quality
methanol (2 mL). The extracted samples were stored in
the dark at 2208C. The recovery rate of this method was 92 6
4%. Analyses were performed using the following conditions: a wavelength of 220 nm, an injection volume of 20 mL,
30-min run times, 0.8 mL min21 flow rates, and elution with anacetonitrile plus phosphate buffer solution (70:30 by volume).
The detection limit for these experimental conditions was
0.1 mg L21.

THEORETICAL ANALYSIS
Governing Flow and Transport Equations
We used the HYDRUS-2D software package (Šimůnek
et al., 1999) to simulate the flow and pesticide transport
processes at the three field sites. Since the various processes in HYDRUS-2D are described in detail in theuser manual, we provide here only a brief synopsis of
those features most relevant to our experiments. Twodimensional variably saturated flow was calculatedusing the standard Richards equation:
"
#
]u
]
A ]h
A
5
K Kij
1 Kiz 2 S
[1]
]t
]xi
]xj

ð

Þ

where u is the volumetric water content (L3 L23), h is thepressure head (L), S is a sink term for root water uptake(T21), xi (i 5 1,2) are the spatial coordinates (L), t is time(T), and Kij are components of the unsaturated hydraulic conductivity tensor (L T21). In this study, we assumed
isotropic media such that the conductivity tensor is diagonal, with both entries, Kxx and Kzz, equal to the unsaturated hydraulic conductivity function, K(h). The soil
water retention, u(h), and hydraulic conductivity, K(h),properties needed in Eq. [1] were described using thefunctions of van Genuchten (1980) as follows:
us 2uru 1h,0nu(h) 5 r[2][1 1 |ah| ]m
us
h$0

5

m 2
K(h) 5 Ks Sle [1 2 (1 2 S1/m
e ) ]

[3]

where ur and us are the residual and saturated watercontents (L3 L23), respectively; Ks (L T21) is the saturated hydraulic conductivity, a (L21) and n (unitless) areempirical shape parameters, m 5 1 2 1/n (unitless), andSe is effective saturation given byu 2 ur[4]
Se ¼
us 2 ur
The effects of macroporosity on variably saturated flow
are often described using composite hydraulic properties (e.g., Durner et al., 1999; Mohanty et al., 1997). In
this study, we used the composite hydraulic conductivity
model of Vogel and Cı́slerová (1988), which assumes
that Eq. [3] can be used for pressure heads less thansome critical value, hk, but with the conductivity increasing linearly from Kk to Ks between hk and saturation to account for noncapillary, macropore-dominatedflow as follows:
Ks Kr (h)h # hk(h 2 hk )(Ks 2 Kk )K(h) 5 Kk 1
hk , h , hs
hs 2 hkKs
h $ hs[5]5where Kk is the hydraulic conductivity at hk. Thisequation assumes that the predicted hydraulic conductivity function is matched to a measured value of thehydraulic conductivity, Kk 5 K(hk), at some water content u(hk) less than us (Vogel and Cı́slerová, 1988). Mostof the hydraulic parameters in Eq. [2] and [3] for thedifferent soil layers at the three sites (Table 2) wereobtained from direct measurements, while others were
estimated from available soil textural data using the
Rosetta pedotransfer functions of Schaap et al. (2001).
Equation [5] was applied only to the more macroporous
Bouzule-2 silty clay site.

842

VADOSE ZONE J., VOL. 5, AUGUST 2006

Table 2. Soil hydraulic and solute transport parameters used in the simulations of the three tile-drained field sites: ur is the residual watercontent, us is the water content at saturation, Ks is the saturated hydraulic conductivity, a and n are empirical parameters, uk is the watercontent at the transition point in hydraulic conductivity function close to saturation, Kk is the hydraulic conductivity at uk, uim is theimmobile water content (set equal to ur), and as is the first-order mass transfer coefficient for exchange between mobile and immobile
liquid regions.
Layer

ur

us

uk

n

a

Ks

21

cm
Villey (sandy loam)

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

Kk

uim

21

as
21

cm d

d

1
2
3
4

0.258†
0.258†
0.258§
0.104§

0.477†
0.477†
0.477§
0.510§–‡–––

0.029§
2.000§
0.031§
2.341§
0.031§
2.341§
0.014§
1.334§
Bouzule-1 (silt loam)

140†
100†
70§
0.1§

–‡
–
–
–

0.258¶
0.258¶
0.258¶
0.104¶

0.05#
0.05#
0.05#
0.05#

1
2
3
4
5

0.179††
0.246††
0.229††
0.235††
0.235††

0.528††
0.504††
0.450††
0.459††
0.459††–––––

0.007§
1.588§
0.008§
1.527§
0.010§
1.419§
0.011§
1.401§
0.011§
1.397§
Bouzule-2 (silty clay)

225‡‡
162‡‡
100‡‡
49‡‡
0.1‡‡

–
–
–
–
–

0.179¶
0.246¶
0.229¶
0.235¶
0.235¶

0.08#
0.08#
0.08#
0.08#
0.08#

1
2
3

0.305††
0.324††
0.308††

0.566††
0.519††
0.505††

0.556
0.509
0.495

0.018§
0.016§
0.015§

25#
10#
1.7#

10‡‡
1.5‡‡
0.1‡‡

0.305¶
0.324¶
0.308¶

0.027#
0.027#
0.027#

1.311§
1.318§
1.322§† Measured data from Bartoli et al. (2006).
‡ Not applicable.
§ Estimated with ROSETTA (Schaap et al., 2001).
¶ Same as residual water content.
# Adjusted to the modeled system.
†† Measured data from Waleed (1983).
‡‡ Measured data from Ailliot (1972).

Solute transport was described using the advection–
dispersion equation (ADE) given by:

ð

Þ

]uc
]rs
]
]c
]qi c
1
5
uDij
2
]t
]t
]xi
]xj
]xi

(i,j 5 1,2)

[6]

where c (M L23) and s (M M21) are solute concentrations in the liquid and solid phases, respectively, qi isrepresents water flux density (L T21), r is the bulk density (M L23), and Dij is the dispersion tensor (L2 T21),
which was described using standard expressions (e.g.,Bear, 1972) in terms of the longitudinal and transverse
dispersivities and contributions arising from molecular
diffusion in the liquid phase (Šimůnek et al., 1999). Degradation of bentazone in Eq. [6] was neglected based onlaboratory experiments that showed little or no degradation within time periods such as those used for thefield experiments (Boivin et al., 2004).
Equation [6] represents the standard ADE, which is
known to have limitations for transport in many fieldsoils (Šimůnek et al., 2003). We additionally considereda dual-porosity (mobile–immobile water) type formulation, which assumes that the liquid phase can be partitioned into mobile, um, and immobile, uim, regions, withadvective–dispersive transport being restricted to themobile region, and with first-order mass transfer between the mobile and immobile liquid regions. The
transport equations for this model, further referred to as
the MIM model, are given byðÞ]um cm]]cm]qi cm5
um Dij2
2 as (cm 2 cim )
]t
]xj
]xi
]xi
]uim cim]rsim
1
5 as (cm 2 cim )
]t
]t

[7a]

[7b]

where cm and cim are the liquid concentrations in themobile (macropore) and immobile (matrix) regions(M L23), respectively, sim is the adsorbed concentrationin the immobile region (M M21), Dij is the dispersion
tensor in the mobile region (L2 T21), and as is the solute
mass transfer coefficient between the two regions (T21).
In this study, we assumed that the immobile water content, uim, in Eq. [7a] and [7b] was constant during thesimulations and that sorption occurred only in the immobile region.
We used a longitudinal dispersivity of 15 cm and atransversal dispersivity of 1 cm for the solute transport
simulations. These values are well within the range ofdispersivities (5–20 cm) generally found for field-scaletransport (Anderson, 1984; Jury and Sposito, 1985;
Beven et al., 1993), and are the same as those used by deVos et al. (1999) in their tile-drained field flow andtransport simulations. Also, following Boivin et al.
(2005), we assumed linear equilibrium partitioning between the adsorbed and solution concentrations in bothEq. [6] and [7] using laboratory-measured kD values.

Numerical Implementation
For each soil profile, a finite element grid was generated using the MeshGen-2D module of HYDRUS-2D(Šimůnek et al., 1999), stratified in accordance with observed soil profile layering (Fig. 2). As an example, theBouzule-2 silty clay loam soil profile was taken to be 4 m
wide (half the spacing between the tile drains) and 1.2 m
deep, with the flow domain being divided into threelayers and assigned soil hydraulic parameters as described in Table 2. A total of 1237 nodes and 2298 triangular elements was used to discretize the Bouzule-2
transport domain, using relatively fine discretizations

843

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

www.vadosezonejournal.org

near the soil surface and around the tile drain to accommodate higher pressure gradients in those parts of thedomain (Fig. 4), and coarser grids elsewhere. The tile
drain was specified at a depth of 90 cm along the side of
the modeled system.
Discretizations of the Bouzule-1 and the Villey sites
were analogous to those in Fig. 4 for Bouzule-2 (grids
not further shown here). The transport domain for the
Bouzule-1 silt loam system was 6 m wide and 1.2 m deep
(Fig. 2). The soil profile at this site was divided into fivelayers with soil hydraulic properties as listed in Table 2.
The flow domain itself was discretized using 1874 nodes
and 3494 triangular elements. The sandy loam soil profile of the Villey site was modeled to be 6 m wide and2.2 m deep, with the soil profile being divided into fivelayers (Fig. 2) using a total of 1833 nodes and 3465 triangular elements.
Boundary conditions for the simulations included noflow boundaries on the sides of the modeled systems,
atmospheric conditions at the soil surface (daily precipitation and evaporation rates), and free drainage (zeropressure gradients) at the bottom boundary. The freedrainage condition reduces to a zero-flux boundary condition when the saturated conductivity of the bottom
layer is made very small or zero. The tile drains were implemented using seepage face boundary conditions, thus
implying that the drains functioned only when the surrounding soil was saturated. Root water uptake was calculated from potential transpiration values using the
reduction model of Feddes et al. (1978), with model uptake parameters for wheat as compiled previously by
Wesseling (1991) in the SWATRE model. Sensitivity
analyses with HYDRUS-2D showed that transpiration
and root water uptake affected mostly water contents
near the soil surface and had little effect on calculated
drain discharge rates and bentazone concentrations in
the drainage water during the 100-d experimental period.
Preliminary flow simulations starting on 10 Jan. 2002(|2 mo before the pesticide application) were performed with an initial equilibrium pressure head profile
so as to obtain realistic pressure head and water content
distributions in the profile at the time of the pesticide
application (4 March at Villey and 11 March at the
Bouzule sites). The initial solute concentration was assumed zero throughout the profile, except for a 1-cm
layer near the soil surface in which the concentration
was defined such that the total solute mass in that layer
was equal to the applied dose (corresponding to theFig. 4. Finite element grid for the transport domain representing half
the drain spacing of the modeled silty clay (Bouzule-2) site. The
grid consisted of 2298 triangular elements and 1237 nodes.

same 1.2 kg ha21 for all three sites). Simulations in allcases continued for 100 d.

RESULTS AND DISCUSSION
Measured and simulated instantaneous and cumulative drain discharge rates are compared in Fig. 5, 6, and 7
for the Villey, Bouzule-1, and Bouzule-2 sites, respectively. These figures also show measured and calculated
bentazone concentrations in the drainage water, as well
as daily and cumulative amounts of bentazone exportedwith the drainage water as a percentage of the appliedamounts. The same scales were used for similar plots inthe three figures to better visualize differences between
the sites. Notice that all three sites produced drainage
during the experiments, with the Bouzule-1 site (Fig. 6b)
producing more drainage than the other two sites. From
10 Mar. to 20 Apr. 2002, the cumulative volume of drainwater measured at the Bouzule-1 silt loam site (Fig. 6b)
was 521 m3 ha21, compared with volumes of 240 and
175 m3 ha21 recorded at the Villey sandy loam (Fig. 5b)and Bouzule-2 silty clay (Fig. 7b) sites, respectively.
Bentazone concentrations in the drainage water at
Bouzule-1 and Bouzule-2 were relatively high, leading
to considerable amounts of pesticide being exported
from those two sites, especially from Bouzule-1 (Fig. 6e)because of the higher drainage volumes at this site thanat Bouzule-2. By comparison, total pesticide export
from the Villey site was relatively small, mostly becauseof lower pesticide concentrations in the drainage water(Fig. 5c and 5e).

Villey Sandy Loam Site
Recent excavations at this site indicated the presenceof a low-permeability marl layer at a depth between 1.8
and 2 m, but no intermediate (mixed) layer between the
sandy deposit and the bottom layer. The excavations
also showed evidence of oxidation–reduction processes
between 1.8- and 2-m depth (L. Florentin, personal communication, 2004), consistent with perched water tableconditions. Following Jacquin and Florentin (1980), theKs of the marl material was set equal to 0.1 cm d21 in the
numerical calculations. The HYDRUS-2D simulations
for the entire soil profile with the marl layer included
correctly produced perched water above the lowconductive marl layer, with the tile drain calculations
closely matching the observed field data (Fig. 5a and 5b).
Simulations with the standard ADE model did notlead to any pesticide in the tile drainage water at theVilley site during the experimental time period (Fig. 5cand 5d). Although the Villey sandy loam soil contained
a large sand fraction (Table 1), significant amounts ofloam, clay, and organic matter were also present, which
promoted good soil structure. Hence, we simulated theexperiment also with the MIM transport model, assuming uim to be the same as ur. We used a relatively highvalue of 0.258 for uim, consistent with results of recent infiltrometer experiments by Bartoli et al. (2006),
which indicated that the residual water contents at thissite were much higher than values predicted with the
Rosetta pedotransfer functions (Schaap et al., 2001).

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

844

VADOSE ZONE J., VOL. 5, AUGUST 2006

Fig. 5. Observed and simulated (a) instantaneous and (b) cumulative tile drain discharge rates, (c) bentazone concentrations in the drainage water,(d) daily amounts of bentazone exported with the drainage water, and (e) the cumulative amount of bentazone exported for the Villey sandy loamsite. Solute transport simulations were performed with both the advection–dispersion equation (ADE) and the mobile–immobile nonequilibriumtransport model (MIM).

The HYDRUS-2D simulation using the MIM transport model with uim 5 0.258 and a single value of
0.05 d21 for the mass transfer coefficient as (visually
fitted to the data) produced excellent agreement withthe observed concentrations (Fig. 5c), as well as with theamount of pesticide exported with the tile drainage
water (Fig. 5d and 5e). The use of the MIM model isqualitatively consistent with previous studies by Boivin
et al. (2004) on the time-dependent availability of bentazone. They showed that, immediately after bentazone
application, 97% of the applied dose was available using
water extraction (0.01 M CaCl2), while an additional
0.1% could be removed only with additional methanol
extractions. On the other hand, after 7 d only 77% wasavailable with water extraction (0.01 M CaCl2) and 19%
using subsequent methanol extraction. Boivin et al.
(2004) concluded that slow solute diffusion into the soil
microporosity may have resulted in entrapment of asignificant part of the chemical, which subsequently
could be removed only by means of a more effective
solvent (methanol).

Bouzule-1 Silt Loam Site
Simulations for the Bouzule-1 site again showed the
key role of the Ks of the bottom layer in determining theamount of drain discharge. A sensitivity analysis (Fig. 6aand 6b) showed that small changes in Ks can cause large
variations in both the instantaneous and cumulative tile
drainage rates. Excellent results were obtained with avalue of 0.1 cm d21 for Ks. The only noticeable difference between the simulation and observed data was aslightly larger calculated drain discharge rate immediately after the rainfall event (Fig. 6a). The optimal valueof 0.1 cm d21 was similar to a laboratory-measured value
reported by Ailliot (1972). The low permeability of the845

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

www.vadosezonejournal.org

Fig. 6. Observed and simulated (a) instantaneous and (b) cumulative tile drain discharge rates, (c) bentazone concentrations in the drainage water,(d) daily amounts of bentazone exported with the drainage water, and (e) the cumulative amount of bentazone exported for the Bouzule-1 silt
loam site. Drain discharge calculations were carried using three values of the saturated hydraulic conductivity (Ks) of the bottom layer 5 (Fig. 2b).
Solute transport simulations were performed with the mobile–immobile nonequilibrium transport model (MIM) using three values of the masstransfer coefficient (as).

bottom layer explains why this site becomes quickly
water saturated during large rain events, and why
artificial tile drains are needed in this area, where theaverage annual precipitation is only |800 mm (mostly
restricted to the fall and winter periods). Previous fieldstudies by Novak et al. (2001) also indicated the presence of a perched water table during part of the year at
this site, leading to considerable tile drainage, especiallyduring the winter season. In our case, having no draindischarge in the beginning of the simulation indicatedthat little or no perched water must have been present
at the time of the pesticide application, with the major
rainfall event after the bentazone application causing
considerable tile drainage flow.
Similarly as for the Villey site, ADE solute transport
simulations that considered only equilibrium transport
did not produce any pesticide transport into the tile
drains at Bouzule-1 (Fig. 6c and 6d). Although we wereable to predict some pesticide leaching by increasing thelongitudinal dispersivity DL (results not further shown),
the total amount of leached bentazone remained severely underestimated. Good agreement between simulated and observed concentrations could be obtained
again with the MIM. Calculated outflow concentrations
were found to be very sensitive to the solute mass transfer coefficient as in Eq. [6]. A sensitivity analysis (Fig. 6c,
6d, and 6e) revealed good agreement between simulated
and observed concentrations and amounts of pesticides
exported when as was set equal to 0.08 d21, somewhat
higher than the optimal value of 0.05 d21 found for theVilley site.

Bouzule-2 Silty Clay Site
The Bouzule-2 site showed less drainage (175 m3 ha21)than the adjacent Bouzule-1 silt loam site (521 m3 ha21),
but higher bentazone concentrations in the drainage
water, probably due to preferential flow. Evidence ofReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

846

VADOSE ZONE J., VOL. 5, AUGUST 2006

Fig. 7. Observed and simulated (a) instantaneous and (b) cumulative tile drain discharge rates, using the standard van Genuchten–Mualem (VGM)
and the modified van Genuchten–Mualem (MVGM) hydraulic functions; (c) bentazone concentrations in the drainage water; (d) daily amountsof bentazone exported with the drainage water; and (e) the cumulative amount of bentazone exported for the Bouzule-2 silty clay site. Solutetransport simulations were performed with both the advection–dispersion equation (ADE) and the mobile–immobile nonequilibrium transportmodel (MIM).

preferential flow at this same site was previously provided by Novak et al. (2001, 2003), who observed large
amounts of both Br2 and the pesticide metolachlor in
the drainage water, especially during the spring andsummer seasons. Novak et al. (2003) concluded from theexperiments that tile drainage water at Bouzule-2 probably came directly from the soil surface through shortcircuiting of much of the soil matrix and rapid leachingof large amounts of pesticides. The Ks values obtained
from our laboratory experiments did not account for
the influence of small drying cracks (up to 2 mm wide)in especially the surface horizons at this site, and thus
substantially underestimated the actual Ks values in thefield. The drying cracks at this site were particularly
noticeable in the late spring and summer months. Direct
field experiments of Ks of the first soil layer (10-cm
depth) using a ring infiltrometer showed the importanceof these drying cracks. When the infiltration ring wasplaced on top of surface cracks, the measured infiltration
rates increased by factors of 2 to 10 relative to areas
without cracks (Bartoli et al., 2006).
Preferential flow at Bouzule-2 could not be simulatedusing the standard van Genuchten–Mualem hydraulicfunctions. Predictions in this case severely underestimated observed drain discharge rates (Fig. 7a and 7b).
To account for flow through the drying cracks, thecomposite soil hydraulic functions of Vogel et al. (2000)
were used to mimic preferential flow of water and solute(modified van Genuchten–Mualem hydraulic functions). Parameters for these functions (Eq. [4] and [5])are listed in Table 2. By using higher values for the hydraulic conductivity near saturation (Ks . Kk), good
agreement was obtained between the simulated andobserved drain discharges (Fig. 7a and 7b). Minor differences between the observed and simulated data still
occurred at the end of the drainage event, with theReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

simulations producing no more tile drainage while thefield data still indicated some drainage.
Using the modified van Genuchten soil hydraulic functions did produce much more realistic simulations of
the drain discharge rate at the Bouzule-2 silty clay site.
While adjustment of the Kk values of the different soil
horizons resulted in a close match with the observed
data at this site, it is unlikely that the selected combination of Kk values in our simulations represented a single
unique description of the flow process. Please note that
we used a trial-and-error approach to obtain the best setof Kk values, rather than the full parameter estimation
capabilities of HYDRUS-2D. While we calibrated thisparameter in our study, it possibly could be measured
also independently using tension infiltrometers, or evenestimated from texture using pedotransfer functions(Jarvis et al., 2002).
Similarly as for the Villey and Bouzule-1 sites, theamount of leached pesticide was severely underestimated unless the MIM model was used. As for the othertwo sites, uim was again assumed to be the same as ur;
however, the mass transfer coefficient as was again fitted
to the data. A value of 0.027 d21 for as produced excellent agreement between measured and simulated
results, both for the bentazone concentrations in the tile
drainage water (Fig. 7c) and the total amount ofbentazone exported with the drainage water. The fitted
value of 0.027 d21 for as is somewhat lower than thevalues of 0.08 and 0.05 d21 found for the other two sites.
These values are all well within the range of values
reported for intermediate-scale and field studies (e.g.,
Pang and Close, 1999; Maraqa, 2001; Jørgenson et al.,2004), and mostly lower than laboratory studies (vanGenuchten and Wierenga, 1976; Maraqa, 2001; Leeet al., 2000). We note that in our study we adopted asingle value for as for each site as a whole. This approach
ignores any dependency of as on the water content or
the flow velocity, which have been shown to occur in several previous studies (e.g., Nkedi-Kizza et al., 1983; De
Smedt and Wierenga, 1984; Maraqa, 2001; Jaynes, 2002).

CONCLUSIONSThe purpose of this study was to compare the transport of the herbicide bentazone in three contrasting tiledrained cultivated field soils subject to otherwise similar
experimental conditions. Observed drain discharge rates
and amounts of pesticide exported from the sites were
mostly a reflection of the soil textural properties andsoil layering at the three sites. Numerical simulations ofdrain discharge with HYDRUS-2D compared well withthe observed drainage data for the relatively homogeneous sandy loam (Villey) and the silt loam (Bouzule-1)
soils. The saturated hydraulic conductivity of the bottom
layer in both cases was key to correctly predicting thedrainage fluxes.
Accurate results could be obtained with the originalvan Genuchten equations for the Villey and Bouzule-1
sites, but not for the more fine-textured Bouzule-2 silty
clay site. Acceptable tile drainage predictions at this sitewere possible only when the soil hydraulic functions

847

were modified to account for preferential flow throughdrying cracks near the soil surface. While the selectedcombination of Kk values in the modified hydraulicproperties (Eq. [6]) provided a good description of thedrain discharge rates, the values probably were notunique. Hence, we acknowledge a need for further studies that more precisely delineate the exact and possibly
time-dependent (seasonal) nature of the preferential
flow processes at especially the Bouzule-2 site. Morecomprehensive long-term data sets may be needed forthis purpose, including in situ measured water contents
and concentrations.
Simulations of pesticide concentrations in the drainage water were unsuccessful when the equilibrium ADE
transport model was used. Much better results were obtained with the nonequilibrium MIM approach, both forthe concentrations of the tile drainage water and total
amounts of solute being exported with the drainage
water. In this study we assumed the immobile watercontent to be the same as the residual water content, anassumption that appeared very much adequate for ourdata; however, as was found to be a critical parameterthat needed to be fitted to the data. Numerical results were very sensitive to the adopted value of as (e.g.,
Fig. 6e), even if we limited this parameter to values wellwithin the range of values previously documented in theliterature. This points to the importance of methods todirectly measure as in the field (e.g., Jaynes, 2002) at spatial and time scales typical of tile drainage experiments.
Our study shows that the use of a limited amountof input data can still lead to acceptable simulations ofdrain discharge rates and chemical concentrations of thedrainage water. Most of the data we used were derived
from previous studies at the sites (Table 2), or using
reasonable assumptions (e.g., for the dispersivities andthe immobile water content). Notable exceptions, again,
were the modified hydraulic parameters (Kk) to accountfor preferential flow at the fine-textured Bouzule-2 site,
and the mass transfer coefficient (as) at all three sites.
These parameters had to be fitted to the field data. Thissuggests that future studies should concentrate especially on these two parameters, including their possible
time dependency during the year as a result of agricultural management practices and shrink–swell phenomena.
Finally, we note that maximum bentazone concentrations in the drainage water ranged between 13.6 mg L21
at Villey to 223 mg L21 at Bouzule-2, while intermediateconcentrations were observed at Bouzule-1 (120 mg
L21). These concentrations are significantly higher thatthe limit of 0.1 mg L21 set by the European Union
for maximum pesticide concentrations in groundwater(Focus Leaching Modelling Workgroup, 1995), andhigher than 3 mg L21 used by the U.S. Environmental
Protection Agency for the maximum contaminant levelfor bentazone (USEPA, 2004). Equally important, thetotal amount of bentazone exported with the drainage
water during our study was very large for two of thethree sites. This was especially the case for the Bouzule-1
silt loam site, where 2.7% of the applied amount ofbentazone ended up in the tile drains (large volumes ofdrainage water with relatively high solute concentrations).

848

VADOSE ZONE J., VOL. 5, AUGUST 2006

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

ACKNOWLEDGMENTS
We wish to thank L. Florentin (LSE) for his valuable advice;
B. Colin and S. Colin (LSE) for theirs skills in building andmaintaining the experimental devices; F. Bartoli (CNRS-LSE)
and J.C. Begin (CNRS-LSE) for the work done characterizing
the hydraulic properties; and R. Cherrier (CRAL, Nancy) for
his help in improving the analytical methods. We are also
grateful to A. Granier (INRA, Champenoux) for providing
the meteorological data, and H.P. Guimont and G. Bertrand(ENSAIA experimental farm, La Bouzule) as well as M.
Grojean (Villey-St-Etienne) for their collaboration. This paper
is based on work supported in part by the SAHRA Scienceand Technology Center through NSF grant EAR-9876800.

REFERENCES
Abbaspour, K.C., A. Kohler, J. Šimůnek, M. Fritsch, and R. Schulin.
2001. Application of a two-dimensional model to simulate flow andtransport in a macroporous agricultural soil with tile drains. Eur.
J. Soil Sci. 52:433–447.
Ailliot, B. 1972. Etude du drainage d’un pélosol avec référence à un sol
brun lessivé marmorisé. Ph.D. diss. Institut National Polytechniquede Lorraine, Nancy, France.
Anderson, M.P. 1984. Movement of contaminants in groundwater:
Groundwater transport—advection and dispersion. p. 37–45. In
Groundwater contamination: Studies in geophysics. Natl. Acad.
Press, Washington, DC.
Banton, O., M. Larocque, P. Lafrance, M. Montminy, and M.A.
Gosselin. 1995. Développement d’un outil d’évaluation des pertes
environnementales de pesticides: Intégration d’un module PestiFlux au logiciel AgriFlux. INRS-Eau, Ste Foy, Québec, Canada.
Bartoli, F., A. Boivin, and M. Schiavon. 2006. Water and herbicidetransient flow transport in field-dried topsoils during controlled
infiltration: I. Water capillary and gravity-driven flows: A preliminary examination. Geoderma (in press).
Bear, J. 1972. Dynamics of fluid in porous media. Elsevier, New York.
Beven, K.J. 1993. Prophecy, reality and uncertainty in distributedhydrological modelling. Adv. Water Resour. 16:41–51.
Beven, K.J., D.E. Henderson, and A.D. Reeves. 1993. Dispersion
parameters for undisturbed partially saturated soil. J. Hydrol. 143:
19–43.
Blazkova, S., K. Beven, P. Tacheci, and A. Kulasova. 2002. Testing the
distributed water table predictions of TOPMODEL (allowing for
uncertainty in model calibration): The death of TOPMODEL?
Water Resour. Res. 38:1257, doi:10.1029/2001WR000912.
Boivin, A., R. Cherrier, C. Perrin-Ganier, and M. Schiavon. 2004.
Time effect on bentazone sorption and degradation in soil. Pest
Manage. Sci. 60:809–814.
Boivin, A., R. Cherrier, and M. Schiavon. 2005. A comparison of five
pesticides adsorption and desorption processes in thirteen contrasting field soils. Chemosphere 61:668–676.
Carrizosa, M.J., M.J. Calderon, M.C. Hermosin, and J. Cornejo. 2000.
Organosmectites as sorbent and carrier of the herbicide bentazone.
Sci. Total Environ. 247:285–293.
Davidson, J.M., and R.K. Chang. 1972. Transport of picloram in relation to soil physical conditions and pore-water velocity. Soil Sci.
Soc. Am. Proc. 36:257–261.
De Smedt, F., and P.J. Wierenga. 1984. Solute transfer through columns of glass beads. Water Resour. Res. 20:225–232.
de Vos, J.A. 2001. Monitoring nitrate leaching from submerged drains.
J. Environ. Qual. 30:1092–1096.
de Vos, J.A., D. Hesterberg, and P.A.C. Raats. 2000. Nitrate leaching ina tile-drained silt loam soil. Soil Sci. Soc. Am. J. 64:517–527.
de Vos, J.A., P.A.C. Raats, and R.A. Feddes. 2002. Chloride transport
in a recently reclaimed Dutch polder. J. Hydrol. 257:59–77.
de Vos, J.A., J. Šimůnek, P.A.C. Raats, and R.A. Feddes. 1999. Identification of the hydraulic characteritics of a layered silt loam.
p. 783–798. In M.Th. van Genuchten et al. (ed.) Characterizationand measurement of the hydraulic properties of unsaturated porous
media. Part 1. Univ. of California, Riverside.
Durner, W., B. Schultze, and T. Zurmühl. 1999. State-of-the-art ininverse modeling of inflow/outflow experiments. p. 661–681. In

M.Th. van Genuchten et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media.
Univ. of California, Riverside.
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation offield water use and crop yield. Simulation monographs. PUDOC,
Wageningen, the Netherlands.
Flury, M. 1996. Experimental evidence of transport of pesticides
through field soils—A review. J. Environ. Qual. 25:24–45.
Focus Leaching Modelling Workgroup. 1995. Leaching models and
EU registration. Guidance document 4952/VI/95. Commission of
the European Communities, Directorate-General for Agriculture
VI B II–I, Brussels, Belgium.
Fox, G.A., R. Malone, G.J. Sabbagh, and K. Rojas. 2004. Interrelationship of macropores and subsurface drainage for conservative
tracer and pesticide transport. J. Environ. Qual. 33:2281–2289.
Gamerdinger, A.P., R.J. Wagenet, and M.Th. van Genuchten. 1990.
Application of two-site/two-region models for studying simultaneous nonequilibrium transport and degradation of pesticides. Soil
Sci. Soc. Am. J. 54:957–963.
Garcia, A.M. 2003. Pesticides exposure and women’s health. Am. J.
Ind. Med. 44:584–594.
Gerke, H.H., and J.M. Köhne. 2004. Dual-permeability modeling ofpreferential bromide leaching from a tile-drained glacial tillagricultural field. J. Hydrol. 289:239–257.
Guo, L., and R.J. Wagenet. 1999. Evaluation of alachlor degradation
under transport conditions. Soil Sci. Soc. Am. J. 63:443–449.
Halfon, E., S. Galassi, R. Bruggemann, and A. Provini. 1996. Selectionof priority properties to assess environmental hazard of pesticides.
Chemosphere 33:1543–1562.
Haws, N.W., P.S.C. Rao, J. Šimůnek, and I.C. Poyer. 2005. Singleporosity and dual-porosity modeling of water flow and solutetransport in subsurface-drained fields using effective field-scaleparameters. J. Hydrol. 313:257–273.
Hendrickx, J.M.H., and M. Flury. 2001. Uniform and preferential flow:
Mechanisms in the vadose zone. p. 149–187. In Conceptual modelsof flow and transport in the fractured vadose zone. Natl. Acad.
Press, Natl. Res. Counc., Washington, DC.
Hopmans, J.W., and J. Šimůnek. 1997. Review of inverse estimation ofsoil hydraulic properties. p. 643–660. In M.Th van Genuchten et al.
(ed.) Characterization and measurement of hydraulic properties of
unsaturated porous media. Univ. of California, Riverside.
Institut Français de l’Environnement. 2003. Pesticides in water.
Annu. Rep. 5th. 2001 data. Institut Français de l’Environnement,
Orléans, France.
Institut Français de l’Environnement. 2004. Pesticides in water.
Annu. Rep. 6th. 2002 data. Institut Français de L’Environnement,
Orléans, France.
Jacquin, F., and L. Florentin. 1988. Atlas des sols de Lorraine. Presses
universitaires de Nancy, Nancy, France.
Jarvis, N.J. 1994. The MACRO model. Version 3.1. Technicaldescription and sample simulations. Rep. and diss. 19. Dep. of
Soil Science, Swedish Univ. of Agricultural Sciences, Uppsala.
Jarvis, N.J., L. Zavattaro, K. Rajkai, W.D. Reynolds, P.-A. Olsen, M.
McGechan, M. Mecke, B. Mohanty, P.B. Leeds-Harrison, and D.
Jacques. 2002. Indirect estimation of near-saturated hydraulic conductivity from readily available soil information. Geoderma 108:
1–17.
Jaynes, D.B. 2002. Determination of uim and a using a disk permeameter. p. 1419–1423. In J.H. Dane and G.C. Topp (ed.)
Methods of soil analysis. Part 4. Physical methods. SSSA Book
Ser. 5. SSSA, Madison, WI.
Jaynes, D.B., S.I. Ahmed, K.-J.S. Kung, and R.S. Kanwar. 2001. Temporal dynamics of preferential flow to a subsurface drain. Soil Sci.
Soc. Am. J. 65:1368–1376.
Jørgenson, P.R., T. Helstrup, J. Urup, and D. Seifert. 2004. Modeling of
non-reactive solute transport in fractured clayey till during variablysaturated flow rate and time. J. Contam. Hydrol. 68:193–216.
Jury, W.A., and G. Sposito. 1985. Field calibration and validation ofsolute transport models for the unsaturated zone. Soil Sci. Soc. Am.
J. 49:1331–1341.
Kohler, A., K.C. Abbaspour, M. Fritsch, and R. Schulin. 2003. Using
simple bucket models to analyze solute export to subsurface drains
by preferential flow. Vadose Zone J. 2:68–75.
Kohler, A., K.C. Abbaspour, M. Fritsch, M.Th. van Genuchten, andReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

R. Schulin. 2001. Simulating unsaturated flow and transport in amacroporous soil to tile drains subject to an entrance head: Modeldevelopment and preliminary evaluation. J. Hydrol. 254:67–81.
Kookana, R.S., R.D. Schuller, and L.A.G. Aylomore. 1993. Simulationof simazine transport through soil columns using time-dependent
sorption data measured under flow conditions. J. Contam. Hydrol.
14:93–115.
Kung, K.-J.S., T.S. Steenhuis, E.J. Kladivko, T.J. Gish, G. Bubenzer,
and C.S. Helling. 2000. Impact of preferential flow on the transportof adsorbing and non-adsorbing tracers. Soil Sci. Soc. Am. J. 64:
1290–1296.
Lagana, A., A. Bacaloni, I. De Leva, A. Faberi, G. Fago, and A.
Marino. 2002. Occurrence and determination of herbicides andtheir major transformation products in environmental waters. Anal.
Chim. Acta 462:187–198.
Larsbo, M., and N. Jarvis. 2005. Simulating solute transport in a
structured field soil: Uncertainty in parameter identification andpredictions. J. Environ. Qual. 34:621–634.
Larsbo, M., S. Roulier, F. Semo, R. Kasteel, and N. Jarvis. 2005. An
improved dual-permeability model of water flow and solute transport in the vadose zone. Vadose Zone J. 4:398–406.
Larsson, M.H., N.J. Jarvis, G. Torstensson, and R. Kasteel. 1999.
Quantifying the impact of preferential flow on solute transport totile drains in a sandy field soil. J. Hydrol. 215:116–134.
Lee, J., R. Horton, and D.B. Jaynes. 2000. A time domain reflectometry
method to measure immobile water content and mass transfer coefficient. Soil Sci. Soc. Am. J. 64:1911–1917.
Maraqa, M.A. 2001. Prediction of mass-transfer coefficient for solutetransport in porous media. J. Contam. Hydrol. 50:1–19.
Miersma, N.A., C.B. Pepper, and T.A. Anderson. 2003. Organochlorine pesticides in elementary school yards along the Texas–Mexico
border. Environ. Pollut. 126:65–71.
Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, J. Šimůnek, andM.Th. van Genuchten. 1998. Preferential transport of nitrate to atile drain in an intermittent-flood-irrigated field: Model development
and experimental evaluation. Water Resour. Res. 34:1061–1076.
Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van
Genuchten. 1997. New piecewise-continuous hydraulic functions formodeling preferential flow in an intermittent-flood-irrigated field.
Water Resour. Res. 33:2049–2064.
Nkedi-Kizza, P., J.W. Biggar, M.Th. van Genuchten, P.J. Wierenga,
H.M. Selim, J.M. Davidson, and D.R. Nielsen. 1983. Modeling
tritium and chloride 36 transport through an aggregated Oxisol.
Water Resour. Res. 19:691–700.
Novak, S.M., O. Banton, and M. Schiavon. 2003. Modeling metolachlor exports in subsurface drainage water from two structured soils
under maize (eastern France). J. Hydrol. 270:295–308.
Novak, S.M., J.-M. Portal, and M. Schiavon. 2001. Effects of soil typeupon metolachlor losses in subsurface drainage. Chemosphere 42:
235–244.
Pang, L., and M. Close. 1999. Field-scale physical non-equilibriumtransport in an alluvia gravel aquifer. J. Contam. Hydrol. 38:447–464.
Pang, L., M.E. Close, J.P.C. Watt, and K.W. Vincent. 2000. Simulationof picloram, atrazine, and simazine leaching through two New
Zealand soils and into groundwater using HYDRUS-2D. J. Contam.
Hydrol. 44:19–46.

849

Pruess, K. 2004. The TOUGH Codes—A family of simulation tools
for multiphase flow and transport processes in permeable media.
Vadose Zone J. 3:738–746.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA:
A computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. J. Hydrol. 251:163–176.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gardenas. 2003.
Review and comparison of models for describing non-equilibrium
and preferential flow and transport in the vadose zone. J. Hydrol.
272:14–35.
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1999. The
HYDRUS-2D software package for simulating the two-dimensional movement of water, heat, and multiple solutes in variablysaturated media. Version 2.0. Int. Ground Water Model. Ctr.
TPS-53. Colorado School of Mines, Golden.
Southwick, L.M., G.H. Willis, D.C. Johnson, and H.M. Selim. 1995.
Leaching of nitrate, atrazine, and metribuzin from sugarcane in southern Louisiana. J. Environ. Qual. 24:684–690.
USEPA. 2004. Drinking water standards and health advisories. Officeof Water, USEPA, Washington, DC.
van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:
892–898.
van Genuchten, M.Th., and P.J. Wierenga. 1976. Mass transfer studies
in sorbing porous media: I. Analytical solutions. Soil Sci. Soc. Am. J.
40:473–480.
van Genuchten, M.Th., M.G. Schaap, B.P. Mohanty, J. Šimůnek, andF.J. Leij. 1999. Modelling flow and transport processes at the local
scale. p. 23–45. In J. Feyen and K. Wiyo (ed.) Modelling of transport
process in soils at various scales in time and space. Wageningen
Pers., Wageningen, the Netherlands.
van Genuchten, M.Th., P.J. Wierenga, and G.A. O’Connor. 1977. Masstransfer studies in sorbing porous media. III. Experimental
evaluation with 2,4,5-T. Soil Sci. Soc. Am. J. 41:278–285.
Vanderborght, J., P. Gähwiller, and H. Flühler. 2002. Identification oftransport processes in soil cores using fluorescent tracers. Soil Sci.
Soc. Am. J. 66:774–787.
Vogel, T., and M. Cı́slerová. 1988. On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention curve.
Transp. Porous Media 3:1–15.
Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2000. Effect of the
shape of the soil hydraulic functions near saturation on variablysaturated flow predictions. Adv. Water Resour. 24:133–144.
Waleed, J. 1983. Comportement hydrodynamique de sols lourds
drainés en fonctions des technologies de drainage. Ph.D. diss. Inst.
Natl. Polytechnique de Lorraine, Vandoeuvre-lès-Nancy, France.
Wesseling, J.G. 1991. Meerjarige simulaties van groundwateronttrekking voor verschillende bodemprofielen, grondwatertrappen en
gewassen met het model SWATRE. Rep. 152. Winand Starring
Center, Wageningen, the Netherlands.
Wong, W.Y., G.A. Zielhuis, C.M.G. Thomas, H.M.W.M. Merkus, and
R.P.M. Steegers-Theunissen. 2003. New evidence of the influenceof exogenous and endogenous factors on sperm count in man. Eur.
J. Obstet. Gynecol. Reprod. Biol. 110:49–54.
Zehe, E., and H. Flühler. 2001. Preferential transport of isoproturon at
a plot scale and a field scale tile-drained site. J. Hydrol. 247:100–115.