Journal of Hydrology (2006) 329, 647– 660

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/jhydrol

Two-dimensional modelling of preferential
water flow and pesticide transport froma tile-drained field
Annemieke I. Gärdenäs a,*, Jirka Šimůnek b, Nicholas Jarvis a,
M.Th. van Genuchten caDepartment of Soil Sciences, Swedish University of Agricultural Sciences, P.O. Box 7014, 750 07 Uppsala, Sweden
Department of Environmental Sciences, University of California, Riverside, 900 University Avenue, A135 Bourns Hall,Riverside, CA 92521, USAc
George E. Brown, Jr. Salinity Laboratory, USDA-ARS, 450 West Big Springs Road, Riverside, CA 92507, USAb

Received 6 April 2005; received in revised form 16 March 2006; accepted 20 March 2006

KEYWORDS

Summary Preferential flow through soil macropores in tile drained soils can significantly
increase the risk of pollution of surface water bodies by agricultural chemicals such as pesticides. While many field studies have shown the importance of preferential flow in tile-drained
fields, few have included detailed numerical modelling of the processes involved. The objectiveof this study was to compare four conceptually different preferential flow and/or transport
approaches for their ability to simulate drainage and pesticide leaching to tile drains. The different approaches included an equilibrium approach using modified hydraulic properties nearsaturation, and three non-equilibrium approaches: a mobile–immobile solute transport model,
a dual-porosity approach, and a dual-permeability formulation. They were implemented into the
HYDRUS-2D software package. The model predictions were compared against measurements ofdrainage and pesticide concentrations made at an undulating, tile-drained field in southern Sweden (Näsbygård) during a period of 6 weeks following spray application of the herbicide MCPA.
The dual-permeability approach most accurately simulated preferential drainage flow, eventhough this approach somewhat overestimated the drainage rates. The equilibrium andmobile–immobile approaches largely failed to capture the preferential flow process. The dualporosity approach predicted much more distinct and higher drainage flow events as comparedto the dual-permeability approach. Differences between measured and simulated tile drainage

Preferential flow
and transport;
Dual-permeability;
Pesticide leaching;
Tile drainage;HYDRUS-2D* Corresponding author. Tel.: +46 18 672294; fax: +46 18 672795.
E-mail addresses: Annemieke.Gardenas@mv.slu.se (A.I. Gärdenäs), Jiri.Simunek@ucr.edu (J. Šimůnek), Nicholas.Jarvis@mv.slu.se(N. Jarvis), rvang@ussl.ars.usda.gov (M.Th. van Genuchten).

0022-1694/$ - see front matter c 2006 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2006.03.021

648

A.I. Gärdenäs et al.
rates using the dual-permeability approach could be partly explained by water by-passing thetile drains and recharging the deeper aquifer. The dual-permeability and dual-porosity
approaches closely captured the dynamics in measured pesticide concentrations. Both the equilibrium and mobile–immobile approaches completely failed to match measured MCPA leachingby underestimating the peak concentrations by more than two orders of magnitude. We conclude that two-dimensional models are suitable tools for studying pesticide leaching from undulating fields with large spatial variability in soil properties.
c 2006 Elsevier B.V. All rights reserved.

Introduction

Rapid non-equilibrium flow through soil macropores is a
widespread phenomenon that can profoundly affect solutetransport into and through the unsaturated zone. The
resulting preferential flow process can dramatically increase the risk of groundwater pollution by surface-appliedchemicals (e.g., see reviews by Flury, 1996; and Jarvis,
2002). The importance of macropore flow for solute transport has been demonstrated repeatedly by means of breakthrough experiments on small, undisturbed soil columns.
Results of such experiments have often been interpreted
in terms of one-dimensional dual-porosity and dual-permeability models appropriate to that scale (e.g. Saxena
et al., 1994; Brown et al., 1999; Castiglione et al., 2003;
Köhne et al., 2004; Pot et al., 2005). However, quantitative
information about the significance of non-equilibrium processes for solute transport at the field-scale is still verymuch lacking. Replicated breakthrough experiments haveshown large variations in transport behaviour between individual columns, ranging from rapid breakthrough due tomacropore flow, to essentially zero leaching (e.g. Jarvis
et al., 1994; Lennartz, 1999). This may cause difficulties
in upscaling from the column to the field scale in a onedimensional framework, since effective macropore flow
parameters are often poorly identifiable (e.g., Czapar
et al., 1992; Kamra et al., 2001; Akhtar et al., 2003).
Two- and three-dimensional modelling approaches
should be more appropriate for field situations, where large
within-field spatial variations in soil properties occur andwhere differences in topography can lead to significant lateral flow along hillslopes. For example, Christiansen et al.
(2004) showed in a recent study of flow and transport at
the catchment scale, that the degree of preferential flow
significantly varied with topography and depth to groundwater within the catchment.
Preferential flow and transport can be especially important in macroporous field soils that are also tile drained,
thus significantly increasing the risk of pollution of nearby
surface water by agricultural chemicals such as pesticides
or nutrients (Stamm, 1997; Kohler et al., 2001). While many
field studies have shown the importance of preferential flowin tile-drained fields, few have included detailed numericalmodelling of the processes involved (Abbaspour et al., 2001;
Kohler et al., 2003; Köhne and Gerke, 2005).
The objective of this study was to compare various preferential flow and/or transport models in terms of their ability to simulate the measured dynamics of water drainage
and pesticide leaching from an undulating, tile-drained agricultural field in Southern Sweden. Four different and com-monly-used conceptual approaches were implemented in
the HYDRUS-2D two-dimensional transport model (Šimůnek
et al., 1999, 2003): (a) an equilibrium approach using soil
hydraulic properties modified close to saturation to accountfor accelerated flow (Vogel and Cı́slerová, 1988), and three
non-equilibrium approaches involving (b) the mobile–immobile water concept for non-equilibrium transport combinedwith standard equilibrium water flow (van Genuchten and
Wierenga, 1976), (c) a dual-porosity approach for bothwater flow and solute transport (Šimůnek et al., 2003),
and (d) a dual-permeability formulation assuming non-equilibrium flow and transport in both domains (Gerke and van
Genuchten, 1993a, 1996).

TheoryThe modified version of HYDRUS-2D (Šimůnek et al., 1999,
2003) includes one equilibrium approach for using modified
hydraulic properties near saturation, and three non-equilibrium approaches i.e., the mobile–immobile water,
dual-porosity and dual-permeability approaches. They havein common that water flow is simulated using the Richards
equation and solute transport using the convection-dispersion equation. The modelling approaches differ in terms of(a) the division of the pore space between matrix and macropore domains, (b) in which part of the porous medium(i.e., the matrix and/or macropore domain) water flow
and solute transport take place, and (c) the definition ofmass transfer, if present, between the two pore domains(Table 1). For example, since the equilibrium approach assumes a single-pore system, this model is strictly not a
preferential flow model in that the approach still predicts
uniform flow and transport. The dual-porosity approach assumes that water flow and transport are restricted only tothe macropore domain, but with possible exchange ofwater and solute between the macropore and matrix

Table 1 Pore regions that are active for water flow andsolute transport using the different modeling approaches

Approach

Water flow

Solute transport

Matrix Macropore Matrix Macropore
Equilibrium+
Mobile–immobile +Dual-porosity
Dual-permeability +

+
+
+

+

+
+
+

Two-dimensional modelling of preferential water flow and pesticide transport from a tile-drained fielddomains. Detailed descriptions of the various modelling approaches are given below separately for water flow andsolute transport.

Water flow
Equilibrium approach with modified hydraulic properties
Variably-saturated flow is calculated using the Richards
equation:


oh
o
oh
ð1Þ
K ij
þ K iz  S¼
ot ox iox jwhere h is the volumetric water content (L3L3), h is thepressure head (L), S is a sink term (T1), xi (i = 1, 2) are
the spatial coordinates (L), t is time (T), and Kij are components of the unsaturated hydraulic conductivity tensor(LT1). We assumed the soil in each layer to be isotropic,
with both entries, Kxx and Kzz, equal to the unsaturated
hydraulic conductivity function, K(h). The van Genuchten–Mualem model (van Genuchten, 1980) was used tocharacterize the unsaturated soil hydraulic properties
hs  hrhðhÞ ¼ hr þ m1 þ jahjnhm i 2KðhÞ ¼ K s K r ðhÞ ¼ K s Sle 1  1  S1=mewhile water and/or solutes are allowed to exchange between the two regions. We will use here the subscript ffor the fracture or macropore region and m for the soil
matrix.
The dual-porosity formulation as used here for variablysaturated flow employs the Richards Eq. (1) to describe flowin the macropores, and a mass balance equation to describe
moisture dynamics in the matrix as follows (Šimůnek et al.,2003):
 

ohf
o
oh
¼
 Sf  C w
K ij
þ K iz
ox i
ox j
ot
ð6Þ
ohm
¼ Sm þ Cw
ot
where Sf and Sm are sink terms for both regions (T1), andCw is the transfer rate for water from the inter- to the intra-aggregate pores (T1). The water mass transfer rate in(6) is assumed to be proportional to the difference in thepressure head between the macropore and matrix regions(Gerke and van Genuchten, 1993b):
Cw ¼ aw ðhf  hm Þ

ð7Þ

ð2Þ

in which aw is a first-order mass transfer coefficient(L1T1):

ð3Þ

aw ¼

where hr and hs denote the residual and saturated watercontents (L3L3), respectively; Ks(LT1) and Kr() are thesaturated and relative hydraulic conductivities, respectively; a is related to the inverse of a characteristic pore radius (L1), n is a pore-size distribution index (), l is a poreconnectivity parameter (), and m = 1  1/n().
The effect of macropores on variably-saturated waterflow is often described using composite functions such as
those suggested by Durner (1994) and Mohanty et al.
(1997), among others. In this study, we use the composite
model suggested by Vogel and Cı́slerová (1988). The approach assumes that capillary-dominated flow can be described using Eq. (3) for the unsaturated hydraulicconductivity at pressure heads less than a critical pressure
head, hk, while a linear increase in conductivity is assumed
between hk and saturation to account for non-capillary,
macropore dominated flow, i.e.:8>< K s K r ðhÞ h 6 hkhh
ð4ÞKðhÞ ¼ K k  hk k ðK s  K k Þ hk < h < 0>:
Ks h P 0
where Kk is the hydraulic conductivity at hk.
Dual-porosity approach
Dual-porosity models assume that water flow is restricted tothe macropore or fracture domain (the inter-aggregate poredomain) and that water in the matrix (the intra-aggregate
pore domain) does not move at all (van Genuchten and
Wierenga, 1976; Šimůnek et al., 2003). The approach assumes that the liquid phase can be partitioned into mobile(flowing), hf, and immobile (stagnant), hm, regions such thatthe total water content is given byh ¼ hf þ hm

649

ð5Þ

b
d2

K a ðhÞcw

ð8Þ

where d is an effective diffusion path-length (i.e. half theaggregate width) (L), b is a shape factor that depends onthe geometry of the soil aggregates (), cw is a scaling factor () (Gerke and van Genuchten, 1993b), and Ka is theeffective hydraulic conductivity of the fracture–matrix
interface (LT1) determined as a simple arithmetic average
involving both hf and hm as follows:
K a ðhÞ ¼ 0:5½K a ðhm Þ þ K a ðhf Þð9ÞDual-permeability approach
Dual-permeability approaches assume that water flow andsolute transport can take place in both the macroporesand the matrix. Some models invoke similar equations forflow in both regions, while others use different formulations(e.g. as in Larsbo and Jarvis, 2003). In this study we use theapproach of Gerke and van Genuchten (1993a, 1996) who
applied the Richards equation to both pore regions. The flow
equations for the macropore (subscript f) and matrix (subscript m) pore systems are then given by:


ohf
o
ohf
Cw
¼
ð10Þ
K fij
þ K fiz  Sf 
ot
ox j
w
ox i


ohm
o
ohm
Cw
¼
ð11Þ
 Sm þ
Km
þ Km
ij
iz
oxi
ot
ox j
1w
where w is the ratio of the volume of the macropore domain
to the total soil volume. Following Gerke and van Genuchten (1993a), the mass transfer term in the dual-permeability
approach is given by Eqs. (7)–(9). Note that the water contents hf and hm in (10) and (11) have different meanings thanin (6) where they represented water contents of the total
pore space (i.e., h = hf + hm), while they refer here to watercontents of the two separate (fracture or matrix) pore domains such that h = whf + (1  w)hm.

650

A.I. Gärdenäs et al.

Solute transport
Equilibrium approach
Solute transport in the physical equilibrium approach is
described using the convection-dispersion equation givenby:
ohc oqsooc
oq c i /
hDij
ð12Þ
þ
¼
ot
ot
ox i
ox j
oxi
where c (ML3) and s (MM1) are solute concentrations in
the liquid and solid phases, respectively, qi is the vector
for the volumetric water flux density (LT1), q is the soil
bulk density (ML3), and Dij is the dispersion tensor(L2T1). In this study, the decay term, / (ML3T1), represents the effects of first-order degradation of MCPA as controlled by the rate constants lw and ls for the liquid andsolid phases, respectively.
Mobile–immobile water
The HYDRUS -2D code (Šimůnek et al., 1999) allows the concepts of mobile and immobile water to be used for solutetransport in combination with equilibrium water flow. As
for the dual-porosity formulation, this approach assumes
that the liquid phase can be partitioned into mobile, hf,
and immobile, hm, regions, with convective-dispersive
transport being restricted to only the mobile region:
ohf cf ofqsf
o
ocf
oq cf
þ
¼
 i  /f  Cs
hf Dfij
ox i
ot
ot
ox j
oxi
ð13Þ
ohm cm oð1  fÞqsm
þ¼ /m þ Cs
otot
where cf and cm are the liquid concentrations in the macropore (mobile water) and matrix (immobile water) regions(ML3), respectively; sf and sm are the adsorbed concentrations in the mobile and immobile regions (MM1), respectively, f is the dimensionless fraction of sorption sites incontact with mobile water (), Dfij is the dispersion tensor
in the mobile region (L2T1), /f and /m are reactions in
the mobile and immobile regions (ML3T1), respectively,
and Cs is the solute transfer rate between the two regions(ML3T1). This transfer rate is assumed to be proportional
to the difference in concentration between the tworegions:
Cs ¼ xs ðcf  cm Þ

ð14Þ

where xs( = him as) is a rate coefficient containing both theimmobile water content (L3L3) and the mass-transfer coefficient as (T1). The above mobile–immobile water solutetransport model assumes that the immobile water content
is constant during the simulations.
Dual-porosity approach
Solute transport in a dual-porosity system can be describedusing the same governing equations as for the mobile–immobile water approach, except that the mass transfer
term now also considers convective solute transfer withwater moving between the two regions, in addition to diffusive mass transfer resulting from concentration gradients:
Cs ¼ xs ðcf  Cm Þ þ Cw c
where c* is equal to cf for Cw > 0 and cm for Cw < 0.

ð15Þ

Dual-permeability approach
The dual-permeability approach assumes that convectivedispersive transport can occur in both the macropore andmatrix domains:
ohf cf oqsf
oocf
oqf cfCsþ
¼
ð16Þ i  /f 
hf Dfij
ot
ot
ox j
ox i
w
ox i


ohm cm oqsm
oocm
oqm cm
Cs
þ
¼
ð17Þ
 i
hm Dm
 /m 
ij
ox i
ot
ot
oxj
ox i
1w
2 1
where Dfij and Dm
ij (L T ) are the dispersion tensors in themacropore and matrix regions, respectively, and qfi and qmi(LT1) are the Darcian flux densities in the macropore andmatrix regions, respectively. The transfer rate, as, for solutes between the macropore and matrix regions (ML3T1)
is the sum of diffusive and convective transport as follows(Gerke and van Genuchten, 1996):

Cs ¼ as ð1  w f Þhm ðcf  cm Þ þ Cw cð18Þ

where c* is equal to cf for Cw > 0 and cm for Cw < 0, and as is a
first-order solute mass transfer coefficient (T1) of the form
as ¼

b
d2

Da

ð19Þ

in which Da is an effective diffusion coefficient (L2T1), representing the diffusion properties of the macropore–matrix
interface.

Boundary conditions
Atmospheric boundary conditions were applied to the soil
surface using (a) prescribed fluxes of precipitation, potential transpiration and evaporation, (b) a prescribed zerohead (saturation) during ponding, and/or (c) equilibrium between the soil surface water content and atmospheric water
vapour when atmospheric evaporative demand could not be
met by the soil. Zero fluxes were ascribed to the lateral andbottom boundaries. Drains were represented by equivalent
nodal sinks in a regular (structured) finite element meshusing an adjusted hydraulic conductivity, Kdrain(LT1) following Vimoke et al. (1963) and Fipps et al. (1986):
K drain ¼ KðhÞCd

ð20Þ

where Cd is a correction factor () determined from the ratio of the effective drain diameter de (L) and the side lengthof the square formed by finite elements with the adjusted
hydraulic conductivity (Fipps et al., 1986).

Material and methods
Study area and soil sampling
The Näsbygård site is located in the southern-most part of
Sweden in the vicinity of the city of Ystad (55°26 0 N;
13°27 0 E, Fig. 1). The climate is warm-temperate (Köppenclassification; Liljequist, 1970) with a mean annual temperature of 7.2 °C and an average annual precipitation rate of
662 mm. The 32-ha field is part of the Vemmenhög catchment (9 km2), formed in a glacial till with Eutric Cambisols(FAO, 1988) being the dominant soil type. The topography
is strongly undulating with height differences of up to
3.5 m within 50 m distances. Previous studies at this fieldTwo-dimensional modelling of preferential water flow and pesticide transport from a tile-drained field

Soil texture (<2 mm) was analyzed by wet-sieving and sedimentation using the pipette method (Ljung, 1987).
The transect was placed perpendicular to the drainage
system to minimize disturbance of the drains, which werelocated 15 m apart and at 1 m depth, but adjusted to the
topography and natural drainage pathways (Gustafson,1987). Wells at the site were used to collect overland flow
(open circles, Fig. 1). Tile drainage discharge of the entire

Precipitation / Drainage (mm)site (Jarvis et al., 2001; Roulier and Jarvis, 2003) showedthat the various landscape elements (i.e., hilltop, mid-slope
and hollows) differ greatly in texture, structure and soil organic matter content, and hence in terms of the risk for pesticide leaching, including the potential for preferential flow
and transport. Jarvis et al. (2001) found considerable spatial
variability in soil texture and organic C content, mostly as afunction of landscape position, while Roulier and Jarvis(2003) found similar dependencies of various macropore
characteristics on landscape position in column tracer
breakthrough experiments calibrated against the MACRO
dual-permeability model. The organic C content was foundto be as high as 5% in the topsoil of the hollows, as compared to 1.2% and 1.4% in the topsoil of the hilltops andmid-slopes, respectively. Pesticide sorption hence should
be more extensive in the hollows, thus locally reducing
the risk of pesticide transport. The hilltops have somewhat
higher clay contents (24%, 17.3%, and 19.8% on average forthe topsoils of hilltop, slope and hollow positions, respectively), leading to lower saturated hydraulic conductivitiesof the matrix, a better developed structure, and hence
more pronounced macropore flow (Roulier and Jarvis,
2003).
A 50-m transect, from a hollow to a hilltop, in the NW–SE
direction, was selected to represent the entire field and toinvestigate short-distance spatial variations in soil texture
and organic matter content (Fig. 1). Soil texture and organic
matter were sampled every two meters at 10–15, 35–40,
60–65 and 85–90 cm depth. The soil organic matter content
was analyzed by loss on ignition at 550–600 °C during 3 h.

651

20 May30 May9 June19 June29 JuneFigure 2 Hourly sums of precipitation and drainage during
the simulation period 16 May–30 June. MCPA spraying occurred
on 21 May 1999.

Figure 1 (a) Location of the field site within Sweden and (b) map of field site with the tile drainage system (tiles are placed at 15m distance intervals and at 1-m depth), and positions of the transect, the tile outlet monitoring station, piezometers and wells atNäsbygård. Hilltops are shown in light grey and hollows in dark grey (source: Gustafson, 1987).

652

A.I. Gärdenäs et al.

field had been monitored at the south-west corner of the
field since 1973 (Fig. 1). ‘Grab’ water samples for pesticide
analyses were taken on 13 occasions during the period May
19–June 19 1999 (Jarvis et al., 2001). The field was cropped
with spring wheat (Triticum aestivum L.), sowed on 4 April
1999. It was conventionally tilled, with autumn ploughing
and spring harrowing prior to planting. The field was sprayedwith the herbicide MCPA on 21 May (40 mg m2) by the farmer using his normal equipment. Total rainfall and tile drainage amounts during the simulation period from 16 May to 30
June were 91 mm and 11 mm, respectively. After a long dry
period with recessing flow, a 29 mm rain storm occurred
within a few hours during the night of 3–4 June. An additional 29 mm of rain fell during the following four days,
which caused peak discharges on 5–7 June 1999 (Fig. 2).

Model applications
Simulated domain, boundary and initial conditions
The simulated domain had a width of 50 m, a depth of 3 m
and a height difference between the hollow and hilltop of
2.5 m (Fig. 3), since the typical length of a sequence of hilltop, mid-slope and hollow was 50 m for this field (Jarvis
et al., 2001), The flow domain contained three drains at
1 m depth and 15 m apart (10, 25 and 40 m from the hollowside). The drains were simulated with the correction factor
Cd and effective drain diameter de in Vimoke’s model set to
4 and 0.04 m, respectively (Šimůnek et al., 1999).
Based on soil texture sampling, the soil profile was divided into 11 different materials (Fig. 3), with the top threelayers each being 25 cm thick, and the bottom layer 2.25 m
thick. The hilltop and mid-slope parts of the profiles shared
the same bottom material, as well as a large part of thethird layer, since the differences in texture here were insignificant (Fig. 3). The flow domain contained 23,128 unstructured triangular finite elements, with higher nodal densities
in the upper 140 cm of the soil profile as well in three 1-m
wide vertical columns around the drains.
The initial level of the groundwater table was set at 1 m
along the hollow side, with a slope of 0.05° in the Eastern
direction, based on groundwater level measurements using
piezometers at four different locations in the field (Fig. 1,
three piezometers along the northern border and one close
to the monitoring station). The piezometric heads were recorded only once a month, which caused some uncertainty
in the initial water level position since the groundwater level decreased considerably during May 1999.

Figure 3

Local climate conditions and root water uptake
Precipitation was recorded each day at 7 a.m. at 500 m distance from the field. Precipitation amounts measured atthat time were attributed to the previous day. We assumed
that all precipitation events occurred at the end of the daywith a maximum rate of 5 mm per hour. Following FOCUS(2000), we assumed that 50% of the pesticide application
was intercepted by the growing crop, and that 10% of theintercepted amount was washed off by each subsequent
mm of rain.
Potential transpiration and evaporation were estimatedwith the Penman–Monteith equation (Monteith, 1965) using
daily values of air temperature, radiation, wind speed, precipitation and the vapour pressure deficit (VPD) obtained
from nearby meteorological stations. Radiation was partitioned between transpiration and soil evaporation according
to Beer’s law as a function of Leaf Area Index (LAI). Canopy
resistance was estimated as a function of VPD, radiation and
LAI according to the Lohammar equation (Lohammar et al.,
1980; Gärdenäs and Jansson, 1995), while aerodynamic
resistance above the canopy was estimated as a functionof crop height and within the canopy as a function of LAI.
Interception evaporation was calculated from LAI, maximum storage capacity per LAI and intercepted water (Jansson and Karlberg, 2004). Potential transpiration wascorrected for interception on a daily basis. Evaporation
equalled the sum of interception and soil evaporation.
Parameterizations of crop characteristics and development,such as LAI (0.1–5), height (0.1–0.8 m), extinction factor
for net radiation (0.5), albedo (25%) and canopy resistancewere taken from Lewan (1993), Blombäck et al. (1995)
and Myrbäck (1998). Soil surface resistance was set to 200
sm1 (Lewan, 1993).
Root water uptake was simulated using the model of Feddes et al. (1978). The maximum rooting depth was assumed
to be 1 m, with the highest root density in the upper 70 cm
of soil. The critical pressure heads in the water-stressresponse function of Feddes et al. (1978) were based on a
study by Wesseling (1991), but adjusted for Swedish conditions according to Kätterer and Andrén (1995), leading toa value of 0.1 m for the critical wet-end pressure headhopt, 15 m for the critical pressure head h3, and 160 m
for the wilting point h4.
Soil hydraulic properties
Hydraulic properties of the 11 soil materials were estimatedusing the soil texture data from the 50-m transect and pedo-

Simulated soil domain with the 11 different soil materials, and locations of the three drains (indicated by bold numbers).

Two-dimensional modelling of preferential water flow and pesticide transport from a tile-drained fieldtransfer functions (see Table 2 for details) for the van
Genuchten–Mualem parameters based on local data by
Svensson (1999) and Jarvis et al. (2001). The saturated
water content, hs, was estimated as a function of bulk density q. The bulk density and the van Genuchten shapeparameter a were calculated as a function of the organiccarbon content and depth, and the pore-size distribution index, n, as a function of texture, organic carbon content anddepth. The pedotransfer function for the saturated hydraulic conductivity, Ks, included both porosity and the 50%
water content, h50, as variables. The residual water content, hr, was assumed to be negligible for all materials
and in all approaches, while the pore-connectivity, l, wasset to 0.5 in all cases.
For the equilibrium approach, the threshold water content, hk, in Eq. (4) was set equal to the water content at
a pressure head, hk, of 10 cm. The hydraulic conductivity at this threshold, Kk, was set equal to the saturated
matrix conductivity, Ksm, in the dual-permeability model,
estimated with a pedotransfer function (see Tables 2 and3).
For the dual-porosity approach, the total saturated
water content, hs, was divided into the mobile saturated
water content, hsf and the immobile saturated water content, hsm, with the assumption that the mobile water content was equal to 0.0225 m3m3 for all soil materials (i.e.
the same saturated macropore water content as for thedual-permeability approach). The water retention coefficients af, and nf for the mobile zone were set to typicalvalues for a coarse-textured soil (10 m1 and 1.8, respec-

Table 2653

tively) so as to represent macropores; these values wereassumed to be the same for all soil materials. The masstransfer coefficient, aw, was based on calibrations carriedout on 20-cm long soil columns sampled from the hollow,
mid-slope and hilltop locations in the same field (Roulier
and Jarvis, 2003).
For the dual-permeability approach, the ratio between
total macropore and total soil volume, w, for all soil layers
was set to 0.03, being the mean value estimated for soilmaterials 5–11 belonging to the slope and hill-top landscape elements. The macropore saturated water content,
hsf, (0.75) was estimated in a similar way (multiplying w
with hsf gives the same saturated macropore water content
per soil volume as for the dual-porosity approach). Following Gerke and van Genuchten (1993a,b), the saturated
hydraulic conductivity of the fracture–matrix interface,Kas, was taken to be 1% of the saturated conductivity ofthe matrix, the shape factor, b, was set equal to 3 ()and the scaling factor, c, to 0.04 (). Values for the masstransfer coefficient, aw, of the various layers were basedon results of Roulier and Jarvis (2003).
Solute transport and reaction properties
The bulk density, q, the dispersion tensor, D, the first-orderdegradation rate constants in liquid and solid phases, lw andls, respectively, and the transfer coefficients, as and xs,
were allowed to vary with landscape position, based on studies at the same field by Jarvis et al. (2001), Roulier and Jarvis(2003), and Lindahl et al. (2005) (Table 4). The transverse
dispersivity, DT , was assumed to be 10% of the longitudinal

Pedotransfer functions for Näsbygård (data from Jarvis et al. (2001) and Svensson (1999))

Parameter

FunctionFit and significanceBulk density, q (g cm3) and porosity, e (m3 m3)

q = 1.464 + 0.05 z  0.08foc
e = 1  q/qs
qs = 2.7  1.7foc

R2 = 0.76, p < 0.0001

Pore-size distribution index, n ()n = l(0.973  0.345q) + 1

R2 = 0.44, p = 0.005

van Genuchten am (m1)

Log a = 1.321  1.838q

R2 = 0.65, p = 0.0001

Saturated matrix water content, hsm (m3 m3)

h10 ¼ hs ð1 þ ð10aÞn Þð1=nÞ1hs ¼ 0:96–0:351q
hsm = h10/(1  w)
w = 0.03

Saturated matrix and macropore conductivity,
Ksm and Ksf (m d1)

K10 = 25dg

R2 = 0.40, p < 0.0001

Log Ks =  1.19 + 2.88 log(100(e  h50))
h50 ¼ hs ð1 þ ð50aÞn Þð1=nÞ1
Ksf = ((Ks  K10)/w).(24/1000)
Ksm = (K10/(1  w)). (24/1000)

R2 = 0.65, p < 0.0001

qs = particle size density (g cm3).
foc = fractional organic carbon content ().
z = horizon number (1 = 0–25 cm, 2 = 25–50 cm, 3 = 50–75 cm, 4 = 75 + cm).
l = particle size distribution index of van Genuchten type equation (, based on Lindahl et al. (2005)).
h10 and h50 = Water content at w = 10 and 50 cm respectively (m3 m3).
hs = Saturated water content (m3m3).
K10 = Hydraulic conductivity at w =  10 cm (cm h1).
dg = characteristic particle size of van Genuchten type equation (mm, based on Lindahl et al. (2005)).

R2 = 0.84, p < 0.0001

654

A.I. Gärdenäs et al.

Table 3 Hydraulic parameters for the matrix (subscript m) and macropore (subscript f) regions of the different soil materials as used in thevarious modeling approaches: (a) equilibrium, (b) dual-porosity and (c) dual-permeability(a) Equilibrium
Landscapeelement

Soil
depth
(cm)

Bulk characteristics
hr (m3m3)hs (m3m3)a (m1)n ()Ks(md1)l ()Threshold characteristics
hk (m3m3)

Kk (md1)

Hollow

0–25
25–50
50–75
75–300

0
0
0
0

0.489
0.473
0.466
0.517

6.5
5.5
5.1
10.0

1.165
1.153
1.153
1.150

1.747
1.107
0.964
0.400

0.5
0.5
0.5
0.5

0.439
0.430
0.425
0.453

0.025
0.022
0.023
0.019

Slope

0–25
25–50
50–75

0
0
0

0.475
0.442
0.417

5.5
3.7
2.7

1.149
1.134
1.128

1.056
0.411
0.209

0.5
0.5
0.5

0.433
0.410
0.391

0.025
0.024
0.025

Slope and
Hilltop75–300

0

0.395

2.1

1.131

0.025

0.5

0.372

0.019

Hilltop

0–25
25–50
50–75

0
0
0

0.476
0.447
0.423

5.6
3.9
2.9

1.129
1.120
1.123

0.815
0.376
0.228

0.5
0.5
0.5

0.436
0.415
0.396

0.024
0.018
0.006

Landscapeelement

Soil
depth
(cm)

Matrixhsm (m m )am (m )nm ()

Ksf (md )hsf (m m )af (m )

nf ()

aw (m1d1)

Hollow

0–25
25–50
50–75
75–300

0.464
0.448
0.441
0.492

6.5
5.5
5.1
10.0

1.165
1.153
1.153
1.15

22.28817.31316.86632.321

0.02250.02250.02250.0225

10
10
10
10

1.8
1.8
1.8
1.8

0.02
0.02
0.02
0.02

Slope

0–25
25–50
50–75

0.450
0.417
0.392

5.5
3.7
2.7

1.149
1.134
1.128

15.7118.142
4.820

0.02250.02250.0225

10
10
10

1.8
1.8
1.8

0.01
0.01
0.01

Slope and
Hilltop75–300

0.370

2.1

1.131

3.305

0.0225

10

1.8

0.01

Hilltop

0–25
25–50
50–75

0.451
0.422
0.398

5.6
3.9
2.9

1.129
1.12
1.123

13.041
7.784
5.813

0.02250.02250.0225

10
10
10

1.8
1.8
1.8

0.001
0.001
0.001

(b) Dual-porosity
Macropore
3

3

1

1

Transfer3

3

1(c) Dual-permeability
Landscapeelement

Soilmaterial(cm)

Hollow

MatrixMacropore

Transfer

hsm(m3m3)

am
(m1)

nm
()

Ksm
(md1)

hsf(m3m3)

af
(m1)

nf
()

Ksf
(md1)w()aw(m1d1)

b
()

c
()

Ksa
(md1)

0–25
25–50
50–75
75–300

0.480
0.480
0.464
0.457

6.5
5.5
5.1
10.0

1.165
1.153
1.153
1.150

0.025
0.022
0.023
0.019

0.75
0.75
0.75
0.75

10
10
10
10

1.8
1.8
1.8
1.8

22.28817.31316.86632.321

0.03
0.03
0.03
0.03

0.01
0.01
0.01
0.01

3
3
3
3

0.4
0.4
0.4
0.4

0.000250.000250.000250.00025

Slope

0–25
25–50
50–75

0.509
0.466
0.432

5.5
3.7
2.7

1.149
1.134
1.128

0.025
0.024
0.025

0.75
0.75
0.75

10
10
10

1.8
1.8
1.8

15.711
8.142
4.820

0.03
0.03
0.03

0.01
0.01
0.01

3
3
3

0.4
0.4
0.4

0.000250.000250.00025

Slope and
Hilltop75–300

0.406

2.1

1.131

0.019

0.75

10

1.8

3.305

0.03

0.01

3

0.4

0.00025

Hilltop

0–25
25–50
50–75

0.384
0.468
0.438

5.6
3.9
2.9

1.129
1.120
1.123

0.024
0.018
0.006

0.75
0.75
0.75

10
10
10

1.8
1.8
1.8

13.041
7.784
5.813

0.03
0.03
0.03

0.01
0.01
0.025

3
3
3

0.4
0.4
0.4

0.000250.000250.00025

dispersivity, DL, consistent with limited literature data (e.g.,
Fetter, 1979). The fraction of the sorption sites in the macropore volume f was taken to be 1% of the total sorption sites(Roulier and Jarvis, 2003). The distribution coefficient Kd for
MCPA was calculated as a function of the organic carbon content and the organic carbon partition coefficient, Koc, of

75 cm3g1. A linear MCPA adsorption isotherm was assumedin accordance with Helweg (1987) and Riise et al. (1994). The
MCPA first-order degradation rate in the solid phase, ls, wasassumed to be equal to that in the liquid phase, lw.
In the mobile–immobile approach, the immobile watercontent, him, was set equal to 0.15. Since him in theTwo-dimensional modelling of preferential water flow and pesticide transport from a tile-drained field

655

Table 4 Values of solute transport and reaction parameters of the different soil materials used in the various modelingapproaches

Landscapeelement

Soil depth All approaches(cm)

Dual-permeability

Mobile–immobileand dual-porosityq (kg m3) DL (m) DT (m) Kd (l kg1) lw (d1)as (d1)xs (d1)

Hollow

0–25
25–50
50–75
75–300

1200
1200
1200
1200

0.1
0.1
0.1
0.1

0.01
0.01
0.01
0.01

1.35
1.45
1.75
3.61

0.33
0.33
0.33
0.33

0.096
0.096
0.096
0.096

0.0336
0.0336
0.0336
0.0336

Slope

0–25
25–50
50–75

1400
1400
1400

0.1
0.1
0.1

0.01
0.01
0.01

0.102
0.734
0.576

0.38
0.38
0.38

0.024
0.024
0.024

0.00840.00840.0084

Slope and
Hilltop75–300

1400

0.1

0.01

0.331

0.38

0.024

0.0084

Hilltop

0–25
25–50
50–75

1500
1500
1500

0.2
0.2
0.2

0.02
0.02
0.02

0.979
0.598
0.410

0.42
0.42
0.42

0.001760.001760.001760.000616
0.000616
0.000616

mobile–immobile approach is constant during the simulations, the total water content can never decrease below thisvalue. This in contrast with the dual-porosity approach,where the immobile water content may vary with time.
The maximum immobile water content in the dual porosity
approach was assumed to be 0.35. The solute mass transfer
coefficient of the dual-permeability approach, as, was estimated from the results of Roulier and Jarvis (2003). Finally,
we note that the rate coefficient, xs, in Eqs. (14) and (15)for the mobile–immobile and dual-porosity approaches,
lumps the effects of the traditional mass transfer coefficient, as, and the immobile (matrix) water content.

Results
Measured and simulated drainage discharge
dynamics
Measured and simulated drainage rates using the different
approaches are shown in Fig. 4. Notice that the vertical axis
in each plot was scaled to the peak discharge rate. No separate results for the water balance of the mobile–immobile
approach are given since the water balance is the same as
for the equilibrium approach. Table 5 lists cumulative drainage amounts for selected periods during the experiment.
Whereas the equilibrium approach simulated the correct
amount of cumulative drainage for the period 3–9 June (Table 5), the largest peak in the simulated drainage flow was
only one third of the measured peak rate (Fig. 4a and b).
More importantly, the equilibrium approach simulated theflow peaks several days too late, thus failing to capture
the preferential flow character of the drainage response.
The dual-permeability and dual-porosity approaches
both simulated the timing of the first drainage flow peakfairly accurately (Fig. 4a, c and d). Simulated results
showed a lag of about 10 h behind the measured data. We
attribute this due to the uncertainty in the precipitation
data, which had been recorded only once a day at 7 a.m.

the following morning. However, the dual-porosity anddual-permeability models both overestimated the measured
peaks in the drainage rates, especially the first peak. Total
cumulative drainage as calculated with the dual-porosityand dual-permeability overestimated the measurements byfactors of 2.5 and 5.6, respectively (Table 5). Notice thatthe dual-porosity approach predicted very high peaks inthe drainage outflow, but of much shorter duration, as compared to the dual-permeability approach. In the dual-porosity approach, all water flow occurs in the fast flow region(see Table 1). For the dual-permeability, some of water volume is exchanged with the matrix which reduces peak flow
rates.
Better agreement between measured and simulated tile
drainage rates could probably have been obtained by adjusting the soil hydraulic properties for each modelling approach separately. This was not done here in view of ouraim to compare the effects of the different conceptual approaches to modelling preferential flow and transport. Separate adjustment of the soil hydraulic properties for eachapproach would have biased the comparisons. Instead, thehydraulic properties for the different approaches wereparameterized as consistently as possible with independentmeasurements and available data. We found that the dualpermeability model best captured the measured dynamics
in the drainage flow rate at our field site, although overestimating total tile drainage flow. Macropore flow dominatedtotal drainage (matrix flow is the difference between total
and macropore flow, see Fig. 4d), even though the macropores were estimated to comprise only 2.25% of the total soilvolume.
One possible explanation for overestimation of the cumulative drainage and peak outflows is that in reality somewater may have bypassed the drains to the deeper subsurface and underlying aquifer, while we used a zero-flux condition at the bottom boundary of the simulated flow
domain. Although the subsoils of these clayey till soils aredense and quite impermeable, Gustafson (1987) made waterbalance calculations showing annual groundwater recharge

656

A.I. Gärdenäs et al.

Figure 4 Measured and simulated tile drainage rates for the different approaches (a) measured, (b) equilibrium model withmodified hydraulic properties, (c) dual-porosity model, and (d) dual-permeability model with total (bold solid line) and macropore
flow (dashed line).

rates of between 20 and 70 mm under this field. If assumed
to be constant over the year, this would have amounted to
between 3 and 8 mm for the 44-day time period of our experiments, which is not enough to explain the relatively large
differences in Table 5. Measurements with, some water
bypassing the drains is very likely. Still, measurements with

the isotope O18 on samples from the piezometer closest to
the drainage monitoring station (see Fig. 1 for location) indicated significant bypass of infiltrating ‘new’ water to shallow
groundwater (Barbro Ulén, Dept. of Soil Science, SLU, Uppsala, pers. comm., 9 Feb., 2005). Another possible reason
for the overestimation is uncertainty in the initial groundwa-

Two-dimensional modelling of preferential water flow and pesticide transport from a tile-drained field
Table 5657

Measured and simulated cumulative drainage amounts for different periods

Approach

Cumulative drainage (mm)Measured

Equilibrium, mobile–immobileDual-porosity

Dual-permeability

16 May–2 June
3–9 June
10–29 June

5.6
1.9
2.6

1.4
2.2
0.1

14.0
31.8
11.2

7.9
7.3
10.2

Total

10.1

3.7

57.0

25.4

ter level. We re-ran the simulations using initial groundwaterlevels 20 cm higher and 10 cm lower than those assumed in
the earlier calculations (1 m). Initial groundwater levels at

depths of 0.8 and 1.10 m resulted in cumulative drainage
rates of 39 and 19 mm, respectively, for the dual-permeability approach (as compared to 25.4 mm for the initial water

Figure 5 Measured and simulated pesticide concentrations of the drainage water for the different approaches (a) measured(crosses) and dual-permeability approach (solid line), (b) equilibrium (solid line) and mobile–immobile (dashed line) and (c) dualporosity.

658
table being at 1 m, and a measured cumulative rate of
10.1 mm). Cumulative drainage hence was found to be quitesensitive to the initial groundwater level.

Measured and simulated pesticide concentrations
Pesticide concentrations (Fig. 5) in the drainage water simulated with the different approaches showed similar patterns versus time as the simulated tile drainage rates. The
equilibrium and mobile–immobile approaches significantly
underestimated the measured pesticide concentrations,and simulated the increase in concentrations much too late
and too slow compared to the measured data. The mobile–immobile approach was able to increase the pesticide concentration by one magnitude compared to equilibriumtransport, but produced only a slightly earlier appearanceof the MCPA in the tile drains (Fig. 5b). By comparison,
the dual-permeability and dual-porosity approaches captured the dynamics of the measured pesticide concentrations quite well. Since the automatic water sampling
equipment failed and samples had to be taken manually instead, some uncertainty existed as to whether or not theevents with the largest concentrations were sampled correctly. For these reasons we did not further calculate
whether and by how much the simulations overestimatedthe measured concentrations. Irrespective of this uncertainty, Fig. 5a show very good agreement between the
dual-permeability modelling results and the measured data.
Results probably could have been further improved by treating the pesticide mass transfer coefficient as a calibration
parameter. Similarly as for the soil hydraulic parameters,
the mass transfer coefficient was taken from Roulier andJarvis (2003).
The sampling for MCPA concentration measurements
started on 19 May 1999 (i.e. two days before MCPA spraying)
and continued until 19 June 1999. During periods of drainage recession it is not unusual to find traces of pesticides
in the tile drainage water from earlier applications. One
measurement just above the detection limit of 0.02 lg/l
was made 2 days before spraying on 19 May. Another measurement on 25 May was below the detection limit, while
several measurements during 15–19 June 1999 were justabove the detection limit. The dual-permeability anddual-porosity models were found to simulate this pattern
well (see close-up in Fig. 5a for the dual-permeability
results).
Similarly as for the drainage water, we expected someoverestimation of the total pesticide load in the tile drainage since groundwater recharge between the drains was neglected in the simulations. Unfortunately, the quality ofpesticide data, together with the uncertainties in the waterbalance, did not allow us to obtain a reliable mass balance.

Discussion and conclusions
The dual-permeability approach involving two mobile flow
regions most accurately simulated the measured dynamics
of pesticide leaching from this undulating, tile-drained field
in southern Sweden, even though it overestimated totaldrainage. Gerke and Köhne (2004) also found that thedual-permeability approach could capture the presumed

A.I. Gärdenäs et al.
non-equilibrium transport processes of a field tracer experiment at a tile-drained field in Germany. An attractive feature of the dual-permeability model is the ability tosimulate both the peak flow resulting from macropore flow,
and the base flow reflecting matrix characteristics. The
equilibrium and mobile–immobile approaches with modified hydraulic properties close to saturation failed to simulate the preferential nature of the flow process, while thedual-porosity approach assuming only mobile water in themacropore volume overestimated the degree of preferential flow.
Differences between measured and simulated tile drainage rates using the dual-permeability approach could be
partly explained by water bypassing the tile drains andrecharging the deeper aquifer. Better agreement between
measured and simulated tile drainage rates could probably
be obtained by adjusting the soil hydraulic properties foreach modelling approach separately. This was not done here
in view of our aim to compare the effects of the different
conceptual approaches to modelling preferential flow andtransport. An alternative approach for model comparison
would be to allow separate parameterizations for each approach using inverse modelling (e.g., Köhne et al., 2004; Coquet et al., 2005). While this is possible at least in theory,
practical problems with ill-posedness are to be expected
for studies of the type considered here, involving a relatively large spatial scale with many different soil materials,
as well as alternative conceptualizations and numbers of unknown parameters (see review by Šimůnek et al., 2003).The
dual-permeability and dual-porosity approaches in thisstudy captured the dynamics in measured pesticide concentrations quite well. Unfortunately, the limited number ofavailable pesticide samples hampered accurate evaluationof the simulated pesticide concentrations, and did not permit accurate mass balance calculations. In hindsight, simultaneously collected data on conservative tracer transport
would have helped us to better identify the major transportmechanisms and pathways as compared to the current study
for which only flow and reactive chemical transport data
were available.
Two-dimensional dual-permeability modelling for applications like the current field experiment can require considerable CPU time. The simulations in our case involved some
23,000 finite elements and 11 soil materials, requiring almost 130,000 iterations for the 1.5-month long simulation.
For these reasons the use of one-dimensional models may
be justified if they can capture the most important processes involved. In our case, the choice between one- andtwo-dimensional modelling depended very much on the possible importance of lateral flow. We knew from earlier studies that the hollows in our experimental field have relatively
high organic matter contents (Jarvis et al., 2001). If transported laterally from the hilltops to the hollows, the pesticides may be adsorbed in the hollows and subsequently
decomposed. We examined the importance of this lateralwater flow using the dual-permeability approach. Simulations showed that the vertical component of macropore
flow dominated in the vadose zone, while in the saturated
zone the lateral macropore flow component was much more
important. The lateral flow velocities in the macropore domain were found to be relatively constant in time during thesimulations. Drain no. 3 at the boundary of the hilltop andTwo-dimensional modelling of preferential water flow and pesticide transport from a tile-drained fieldmid-slope area remained in the unsaturated zone during theentire simulation. Macropore flow from the hilltop area
hence did not end up in the nearest tile drain, but instead
moved mostly vertically into the saturated zone and thenlaterally to drains 1 and 2. Clearly, these lateral flow processes cannot be captured using a one-dimensional model.
Another suitable process to study with a two-dimensionalmodel in a future project would be to evaluate the importance of the position of the drains in order to minimize
leaching of agricultural chemicals. For example, shifting
the drains in our study several meters towards the hollow
may reduce pesticide leaching since there the flow velocities are much smaller, while the potential for pesticides
sorption is much higher. We conclude that two-dimensionalmodels are suitable tools for studying pesticide leachingfrom undulating fields with large spatial variability in soil
properties.

Acknowledgements
We wish to acknowledge several colleagues at the Dept. ofSoils Sciences, Swedish University of Agricultural Sciences:
Christina Öhman for her valuable help with soil physicalanalyses, Erarso Etana and Laurant Brodeau for assistancein the field, and Göran Johanson and Barbro Ulén for supplying data and valuable input in the discussions. We also liketo acknowledge Näsbygård (Kjell Thuresson) for giving us access to the site and supplying crop and management data.
Swedish Research Council (VR, formerly NFR) supported
the study financially. This study was also partly supportedby the SAHRA Science Technology Center as part of NSF
grant EAR-9876800, and by Carl Von Horn’s Foundation
and Formas through travel grants. Two anonymous reviewers are acknowledged for many useful suggestions that improved the manuscript.

References
Abbaspour, K.C., Kohler, A., Šimůnek, J., Fritsch, M., Schulin, R.,
2001. Application of a two-dimensional model to simulate flow
and transport in a macroporous agricultural field with tile drains.
Eur. J. Soil Sci. 52, 433–447.
Akhtar, M.S., Steenhuis, T.S., Richards, B.K., McBride, M.B., 2003.
Chloride and lithium transport in large arrays of undisturbed silt
loam and sandy loam soil columns. Vadose Zone J. 2, 715–727.
Blombäck, K., Stähli, M., Eckersten, H., 1995. Simulation of waterand nitrogen flows and plant growth for a winter wheat stand incentral Germany. Ecol. Model. 81, 157–167.
Brown, C.D., Marshall, V., Deas, A., Carter, A.D., Arnold, D., Jones,
R.L., 1999. Investigation into the effect of tillage on solute
movement to drains through a heavy clay soil. II. Interpretation
using a radio-scanning technique, dye tracing and modelling. Soil
Use and Management 15, 94–100.
Castiglione, P., Mohanty, B.P., Shouse, P.J., Šimůnek, J., vanGenuchten, M.Th., Santini, A., 2003. Lateral water diffusion in
an artificial macroporous system: Modeling and experimental
evidence. Vadose Zone J. 2, 212–221.
Christiansen, J.S., Thorsen, M., Clausen, T., Hansen, S., Refsgaard,
J.C., 2004. Modelling of macropore flow and transport atcatchment scale. J. Hydrol. 299, 136–158.
Coquet, Y., Šimůnek, J., Coutadeur, C., van Genuchten, M.Th., Pot,
V., Roger-Estrade, J., 2005. Water and solute transport in a659

cultivated silt loam soil: 2. Numerical analyses. Vadose Zone J.
4, 587–601.
Czapar, G.F., Horton, R., Fawcett, R.S., 1992. Herbicide and tracer
movement in soil columns containing an artificial macropore. J.
Environ. Qual. 21, 110–115.
Durner, 1994. Hydraulic conductivity estimation for soils withheterogenous pore structure. Water Resour. Res. 30, 211–224.
FAO, 1988. FAO/Unesco Soil map of the world. Revised legend.
World Resources Report 60, FAO, Rome, Reprinted as TechnicalPaper 20, ISRIC, Wageningen, 1989.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of fieldwater use and crop yieldSimulation Monographs. Centre for
Agricultural Publishing and Documentation, p. 189.
Fetter, C.W., 1979. Contaminant Hydrology. MacMillan Publ. Co.,New York.
Fipps, G., Skaggs, R.W., Nieber, J.L., 1986. Drains as a boundary
condition in finite elements. Water Resour. Res. 22, 1613–1621.
Flury, M., 1996. Experimental evidence of transport of pesticides
through field soils – a review. J. Environ. Qual. 25, 25–45.
FOCUS, 2000. ‘‘FOCUS groundwater scenarios in the EU review ofactive substances’’ Report of the FOCUS Groundwater Scenarios
Workgroup, EC Document Reference SANCO/321/2000 rev. 2, p.
202.
Gärdenäs, A., Jansson, P.-E., 1995. Simulated water balance of
Scots Pine stands in Sweden for different climate change
scenarios. J. Hydrol. 166, 107–125.
Gerke, H.H., van Genuchten, M.T., 1993a. A dual-porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29, 305–319.
Gerke, H.H., van Genuchten, M.T., 1993b. Evaluation of a firstorder water transfer term for variably saturated dual-porosity
flow models. Water Resour. Res. 29, 1225–1238.
Gerke, H.H., van Genuchten, M.Th., 1996. Macroscopic representation of structural geometry for simulating water and solute
movement in dual-porosity media. Adv. Water Resour. 19, 343–
357.
Gerke, H.H., Köhne, J.M., 2004. Dual-permeability modeling of
prefrential bromide leaching from a tile-drained glacial tillagricultural field. J. Hydrol. 289, 239–257.
Gustafson, A., 1987. Simulation of water discharge rates from aclay-till soil over a ten year period. J. Hydrol. 92, 263–274.
Helweg, A., 1987. Degradation and adsorption of 14C-MCPA in soil –Influence of concentration, temperature and moisture content
on degradation. Weed Res. 27, 287–296.
Jansson, P.-E., Karlberg, L., 2004. Coupled heat and mass transfer
model for soil–plant–atmosphere systems. TRITA-LWR Report
3087, p. 427.
Jarvis, N.J., 2002. Macropore and preferential flow. In: Plimmer, J.
(Ed.), The Encyclopedia of Agrochemicals, vol. 3. J. Wiley andSons, Inc., pp. 1005–1013.
Jarvis, N.J., Stähli, M., Bergström, L., Johnsson, H., 1994. Simulation of dichlorprop and bentazon leaching in soils of contrasting texture using the MACRO model. J. Env. Sci. Health A29,
1255–1277.
Jarvis, N., Gärdenäs, A., Alavi, G., Kreuger, J., 2001. In: Holman, I.,
Hollis, J. (Eds.), CAMSCALE: upscaling, predictive models andcatchment water quality. Final Report ENV4-CT97-0439, Soil
Survey and Land Research Centre, Silsoe, United Kingdom. pp.
12–18.
Kamra, S.K., Lennartz, B., van Genuchten, M.Th., Widmoser, P.,
2001. Evaluating non-equilibrium solute transport in small
columns. J. Contam. Hydrol. 48, 189–212.
Kätterer, T., Andrén, O., 1995. Measurements and simulations ofheat and water balance components in a clay soil cropped withwinter wheat under drought stress or daily irrigation andfertilization. Irrig. Sci. 16, 65–73.
Kohler, A., Abbaspour, K.C., Fritsch, M., van Genuchten, M.Th.,
Schulin, R., 2001. Simulating unsaturated flow and transport in a

660
macroporous soil to tile drains subject to an entrance head:model development and preliminary evaluation. J. Hydrol. 254,
67–81.
Kohler, A., Abbaspour, K.C., Fritsch, M., Schulin, R., 2003. Using
simple bucket models to analyze solute transport to subsurface
drains by preferential flow. Vadose Zone J. 2, 68–75.
Köhne, J.M., Köhne, S., Mohanty, B.P., Šimůnek, J., 2004. Inverse
mobile–immobile modeling of transport during transient flow:
Effects of between-domain transfer and initial water content.
Vadose Zone J. 3, 1309–1321.
Köhne, J.M., Gerke, H.H., 2005. Spatial and temporal dynamics ofpreferential bromide movement towards a tile drain. Vadose
Zone J. 4, 79–88.
Larsbo, M., Jarvis, N., 2003. MACRO 5.0. A model of water flow andsolute transport in macroporous soil. Technical description.
Emergo 2003:6. Swedish. Univ. of Agric. Sci., Dept. of Soil Sci,
Div. of Env. Phys. ISBN 91-576-6592-3.
Lennartz, B., 1999. Variation of herbicide transport parameters
within a single field and its relation to water flux and soil
properties. Geoderma 91, 327–345.
Lewan, E., 1993. Evaporation and discharge from arable land withcropped or bare soils during winter. Measurements and simulations. Agric. Forest Meteorol. 64, 131–159.
Liljequist, G.H., 1970. Climatology. Generalstabens Litografiska
Anstalt, Stockholm, p. 427. (In Swedish).
Lindahl, A., Kreuger, J., Stenström, J., Gärdenäs, A., Alavi, G.,
Roulier, S., Jarvis, N., 2005. Stochastic modelling of diffuse
pesticide losses from a small agricultural catchment. J. Environ.
Qual. 34, 1174–1185.
Ljung, G., 1987. Mekanisk analys. Beskrivning av en rationell metodför jordartsbestämning. Swed. Univ. of Agric. Sci., Dept. of Soil
Sci., Div. Agric. Hydrotechnics, Communication 87:2, p. 9 (InSwedish).
Lohammar, T., Larsson, S., Linder, S., Falk, S.O., 1980. FAST –simulation models of gaseous exchange in Scots pine. In Persson,
T. (Ed.), Structure and Function of Northern Coniferous Forests– An Ecosystem Study, Stockholm, Ecol. Bull. vol. 32, pp. 505–
523.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten,
M.Th., 1997. New piecewise-continuous hydraulic functions formodeling preferential flow in an intermittent flood-irrigatedfield. Water Resour. Res. 33, 2049–2063.
Monteith, J.L., 1965. Evaporation and environment. In: Fogg, G.E.
(Ed.), The State and Movement of Water in Living Organisms,
19th Symp. Soc. Exp. Biol., Cambridge, The Company ofBiologists, pp. 205–234.
Myrbäck, Å., 1998. Swedish agricultural and horticultural crops.
Swedish National Chemical Inspectorate PM 1/98. p. 44.

A.I. Gärdenäs et al.
Pot, V., Šimůnek, J., Benoit, P., Coquet, Y., Yra, A., Martı́nezCordón, M.-J., 2005. Impact of rainfall intensity on the transportof two herbicides in undisturbed grassed filter strip soil cores. J.
of Contaminant Hydrology 81, 63–88.
Riise, G., Eklo, O.M., Pettersen, M.N., Salbu, B., 1994. Associationof MCPA, dichlorprop, tribenuron-methyl, atrazine and dimethoate with different soil types – Laboratory experiments. Norw.
J. Agric. Sci. 13 (supplement), 17–29.
Roulier, S., Jarvis, N., 2003. Modeling Macropore Flow Effects on
Pesticide Leaching: Inverse Parameter Estimation using MicroLysimeters. J.Environ. Qual 32, 2341–2353.
Saxena, R.K., Jarvis, N.J., Bergström, L., 1994. Interpreting nonsteady state tracer breakthrough experiments in sand and clay
soils using a dual-porosity model. J. Hydrol. 162, 279–298.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1999. The HYDRUS2D software package for simulating two-dimensional movementof water, heat, and multiple solutes in variably saturated media.
Version 2.0, IGWMC – TPS – 53, International Ground Water
Modeling Center, Colorado School of Mines, Golden, Colorado, p.
251.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.Th., Gärdenäs, A.,
2003. Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadosezone. J. Hydrol. 272, 14–35.
Stamm, C., 1997. Rapid transport of phosphorus in drained
grasslands soils. PhD Dissertation No. 12486, Swiss Federal
Institute of Technology, Zürich.
Svensson, O., 1999. Markkarakterisering av ett avrinningsområde isödra Skåne. Division of Water Quality Management, Swedish
University of Agricultural Sciences, Uppsala, Sweden. Seminars
and MSc thesis no. 31, p. 15. (in Swedish with English abstract).
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., Wierenga, P.J., 1976. Mass transfer studies
in sorbing porous media, I. Analytical solutions. Soil Sci. Soc.
Am. J. 40, 473–481.
Vimoke, B.S., Tura, T.D., Thiel, T.J., Taylor, G.S., 1963. Improvements in construction and use of resistance networks forstudying drainage problems. Soil Sci. Soc. Am. Proc. 26 (2),203–207.
Vogel, T., Cı́slerová, M., 1988. On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention
curve. Transp. Porous Media 3, 1–15.
Wesseling, J.G., 1991. Meerjarige simulaties van grondwateronttrekking voor verschillende bodemprofielen, grondwatertrappen
en gewassen met het model SWATRE, Report 152, WinandStaring Centre, Wageningen (In Dutch).