Journal of Hydrology 321 (2006) 390–408
www.elsevier.com/locate/jhydrol

Bromide transport at a tile-drained field site: experiment,
and one- and two-dimensional equilibrium andnon-equilibrium numerical modeling
Sigrid Köhnea,b,*, Bernd Lennartzb, J. Maximilian Köhneb, Jirka ŠimůnekaaDepartment of Environmental Sciences, University of California, Riverside, 900 University Avenue,
A135 Bourns Hall, Riverside, CA 92521, USAb
Institute for Land Use, Faculty for Land use and Agricultural Sciences, University Rostock,
Justus-von-Liebig-Weg 6,18044 Rostock, Germany
Received 7 October 2004; revised 6 August 2005; accepted 11 August 2005

Abstract
Systematically tile drained field sites have been recognized as one major source for surface water contamination with
agrochemicals. To study the effects of tile drainage and physical non-equilibrium on solute transport in structured soil, bromide(BrK) transport experiments were carried out on three plots (A, B, C) with different tile drain depths (128, 101, 96 cm) andspacings (16, 18 and 12 m) at the Infeld experimental field site (North–West Germany). Tile drain outflow along with BrKconcentrations were monitored over a half-year period. For all three plots, fast BrK breakthrough was observed with BrKconcentrations fluctuating around levels below 6 mg/L during the experiment without a distinct concentration maximum.
Experimental observations of plot B were analyzed using one-dimensional (1D) and two-dimensional (2D) single-porosity
model (SPM) and mobile-immobile model (MIM) approaches. All SPM and MIM parameters were obtained from independentmeasurements, except for the calibrated MIM water and solute transfer coefficients. Water flow and BrK transport were thenpredicted for the A and C plots using the model parameters as obtained for plot B. Measured plateau-like BrK concentrations
could only be consistently calibrated (plot B) and predicted (A and C) using the 2D-MIM approach, while the 1D-MIM, 2DSPM, and 1D-SPM approaches (in this order) increasingly deviated from the experimental data. The MIM simulations
suggested that solute transfer into the immobile region represented more than 60% of the surface applied BrK. Non-equilibriumtransport with advective and diffusive mass transfer increased early mass loss and entailed extended, slower leaching
as compared to equilibrium transport. Furthermore, model simulations suggested that the two-dimensional flow field as
induced by tile drains enhanced BrK dispersion and accelerated BrK appearance in the drain. This study showed that both

* Corresponding author. Address: Institute for Land Use, Faculty for Land use and Agricultural Sciences, University Rostock, Justus-vonLiebig-Weg 6,18044 Rostock, Germany. Tel.: C1 979 779 3430.
E-mail addresses: sigrid.koehne@uni-rostock.de (S. Köhne), bernd.lennartz@uni-rostock.de (B. Lennartz), max.koehne@uni-rostock.de(J.M. Köhne), jiri.simunek@ucr.edu (J. Šimůnek).

0022-1694/$ - see front matter q 2005 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2005.08.010

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

391

the variably-saturated 2D flow field and physical non-equilibrium transport should be explicitly accounted for in physically
based model simulations of solute transport in tile-drained structured field soils.
q 2005 Elsevier B.V. All rights reserved.
Keywords: Field-scale experiment; Subsurface-drainage; Physical non-equilibrium; Preferential flow and transport; Mobile-immobile model;
two-dimensional modeling

1. Introduction
The loss of nutrients and pesticides from agricultural fields to subsurface drainage systems impairs
ground and surface water quality (Richard andSteenhuis, 1988; Kladivko et al., 1991; Southwick
et al., 1992). Physical non-equilibrium transportconditions initiated by preferential flow resulted into
faster solute transport towards drains, and caused
higher peak concentrations, than was expected underequilibrium conditions (e.g. Kladivko et al., 1991;
Mohanty et al., 1998a,b; Lennartz et al., 1999; Gishet al., 2004). In this context, the term physical nonequilibrium is defined as local differences of watercontents, pressure heads, and solute concentrations
between bypasses (also denoted fractures, mobile ormacropore regions/domains) and the soil matrix(Jarvis, 1998; Gerke and Van Genuchten, 1993).
The extent of bypassing and physical non-equilibrium
during flow and transport varies depending on various
factors, such as soil water content, soil texture andstructure, tillage, precipitation amount and intensity,and lag-time between solute application and rainfall(e.g. White, 1986; Bouma, 1990; Jaynes et al., 2001;
Fortin et al., 2002).
Physically based numerical models are frequently
used to analyze and predict solute transport in tile
drained fields. Some models, for instance DRAINMOD (Skaggs, 1980) or SWMS_2D (Šimůnek et al.,1994), consider the drain-induced two-dimensional
flow field, but are based on the single-porosity
assumption, where water contents and pressure
heads are always at local equilibrium. Such models
were used by, among many others, Mohanty et al.
(1998a,b) or Šimůnek and de Vos (1999) to simulate
nitrate leaching from tile-drained fields. Other models
relied on the dual-permeability concept to account forphysical non-equilibrium transport, but were constrained to only one dimension (1D). For example,
MACRO (Jarvis, 1994) was applied by Larsson andJarvis (1998), Larsson et al. (1999), and Andreu et al.

(1994) to simulate physical non-equilibrium solutetransport data from subsurface drained field plots. The
1D MACRO model approximates the effect of 2D
saturated groundwater flow using a specific lower
boundary condition. Recently, Kohler et al. (2001)
developed and applied a 2D numerical dual-permeability model for analysis of preferential solutetransport in tile-drained fields assuming 1D verticalflow in the macropore domain. There is still a lack ofstudies comparing 1D and 2D physical equilibriumand non-equilibrium models to systematically separate the spatial and preferential flow effects on solutetransport at tile-drained field sites.
A common model for describing physical nonequilibrium solute transport in soils is the mobileimmobile model (MIM), also named dual-porosity
model or two-region model (Van Genuchten and
Wierenga, 1976). The MIM assumes that the soils’
liquid phase is composed of a mobile region, in whichtransport is governed by convection and dispersion,and an immobile (stagnant) region. The two regions
can exchange solutes by means of molecular
diffusion. The MIM concept has recently been
extended to consider transient, variably saturated
water flow in the mobile region, and exchange of bothwater and solutes with the immobile region (Šimůnek
et al., 2003). The MIM has the advantage of requiringconsiderably fewer parameters than the dual-permeability models with two interacting mobile regions(Gerke and Van Genuchten, 1993; Jarvis, 1994).
A variably saturated flow, simplified 1D MIM
model with drain-boundary condition (an extension of
the MACRO model) was used by Larsson et al. (1999)
to analyze bromide transport data from a tile-drained
field. Abbasi et al. (2004) applied a 2D MIM to
analyze data from furrow irrigation experiments.
However, the immobile water content was constant,preventing water transfer between mobile andimmobile regions and allowing only for diffusive
solute transfer between the regions (Abbasi et al.,2004). In this study, the MIM approach for variably

392

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

saturated flow and transfer of water and solute
between the regions (Šimůnek et al., 2003) is
extended to two dimensions and is implemented into
the HYDRUS-2D model (Šimůnek et al., 1999).
The objective of this study was to investigate theeffects of the tile-drain induced, two-dimensional flow
field and of physical non-equilibrium on solutetransport through a structured tile-drained soil. We
approached our objective by first performing BrK
tracer tests on three plots with different tile drain
spacing and depths at the Infeld tile-drained field site(Germany). These tests were then followed by thecomparison of the experimental observations withsimulation results obtained with one- and twodimensional, single- and dual-porosity models.

2. Material and methods
2.1. The Infeld study site
The Infeld experimental station is located near
Nordenham (Germany) in a Marsh-region that wassystematically drained by open ditches and tile drains
several decades ago. At the pasture cropped Infeld
site, nutrient losses via tile drains have beenmonitored since 1975 by the county agricultural
authority Landwirtschaftskammer Weser-Ems. The
Infeld site (4.16 ha) is subdivided into 10 plots (0.1–
0.24 ha) with different tile drain depths (79–130 cm)
and spacings (10–18 m). At each plot, three tile-lines
are routed to a measurement station. Plots arehydrologically isolated by surrounding discard drains,which are connected to the ditch system (see Fig. 1).
Since 1975, drainage from the individual plots has
been monitored continuously using tipping bucket
samplers.
The soil at Infeld is a Fluvi-Dystric Gleysol
according to FAO-classification. It was described indetail by Neuhaus (1983) and Diez and Weigelt(1987). Four main horizons were identified: H1 (0–25 cm depth), H2 (25–70 cm), H3 (70–135 cm), andH4 (135–200 cm). The H1 surface layer wascharacterized by subangular to crumb structure,
whereas the subsoil layers showed angular toprismatic (H2) or dense prismatic (H3, H4) structure.
The soil characteristics do not vary spatially over theexperimental site. Soil texture as determined with a

combined sieve and sedimentation method (German
Industrial Standard DIN 19683 T2) was very fine withmore than 40% clay content in each horizon (Table 1).
Additional soil properties are listed in Table 1.
In 1999, six undisturbed 250-cm3 samples weretaken from each soil horizon by driving stainless steel
cylinders vertically into the soil. Additionally, bulk
soil material was collected. Water retention curves
were measured by stepwise desaturation at 30, 50,
100, 150, 200, 300 and 600 cm suctions and in apressure chamber at 15000 cm. The soil hydraulic van
Genuchten’s (1980) parameters a, n, and the saturated
water content, qs, were estimated from the retention
data using the RETC program (van Genuchten et al.,1991). The residual water content, qr, was fixed at
60% of the average water content of a particular soil
horizon as measured in the pressure chamber (at15000 cm). Additionally, the bulk density of each soil
horizon was calculated as the arithmetic average ofdry soil weights divided by sample volumes (Table 1).
Saturated hydraulic conductivities, Ks, of H1 andH2 were determined using the volume samples and afalling head permeameter. To determine the depthaveraged saturated hydraulic conductivity of the(combined H3 and H4) subsurface layers, the auger
hole method according to Hooghout and Ernst (Twitty
and Rice, 2001) was applied at seven boreholes (230cm depth, 8-cm diameter). During the measurement
period, the hydrostatic groundwater level in sevenholes was located in the third (H3) layer and during
the auger hole tests was lowered into the fourth (H4)layer.
2.2. Bromide transport experiments
A bromide tracer experiment was carried out in thewinter of 1999/2000 on the plots N2, N4 and N8 of theInfeld experimental station. The plots N2, N4, and N8
are renamed A, B, and C for simplification from here
onwards. A solution carrying 65.5 g/L bromide (BrK)as potassium bromide (KBr) was sprayed on
November 22th 1999 onto each plot with a volumeof 0.066 L/m2. At the time of application, tile drainage
had not started yet, and the water table was below the
deepest drain level at Infeld of 1.30 m. Table 2
summarizes the most important characteristics of theexperiment. Rainfall (Hellmann rain gauge at 1 mheight), relative humidity, temperature, wind speed,

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

393

Fig. 1. Location of the Infeld experimental station in Germany and schematic view of the site.

and sunshine duration were measured daily. The
rainfall data were corrected according to Richter(1995) to account for measurement errors caused bywind, evaporation and wetting of the Hellmann rain
gauge. The potential evapotranspiration was calculated using the Turc-Wendling method (Wendlinget al., 1991), while the actual evapotranspiration(ETa) was determined employing the soil waterbalance model of Wendling (DVWK, 1996). Fig. 2
shows the atmospheric conditions at the site during
the experimental period. We note that the experiment
was conducted in a relatively mild winter, for whichdaily mean temperatures below 0 8C were recorded
only on 3 days (12/15/99, K0.3 8C; 1/10/00, K0.2 8C;
1/23/00, K1.8 8C).
Drainage outflow (cm/d) was measured daily forplot A and every 6 h for plots B and C. Water samplesof 90 cm3 each were collected in 15-minute timeintervals from the drainage effluent and wereaccumulated to 6-hour samples. The 6-hour samples
were analyzed for BrK concentration at the B and C
plots, whereas for the A plot, four (6-hour) samples ofa particular day were mixed to yield daily samples of
BrK concentrations.

394

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

Table 1
Selected soil properties at Infeld experimental site.
Horizon

Clay
(%)

Silt
(%)

Sand
(%)

Corg
(%)

Bulk density(g cm-1)

pH
(K)

H1
H2
H3
H4

42.9
57.4
50.8
43.7

47.4
38.7
45.2
53.1

10.0
4.3
3.6
1.1

4.6
0.9
0.9
n.d.

1.35
1.31
1.09
0.91

5.9
6.8
6.9
n.d.

Equivalent diameter of Clay, 0–0.002 mm; Silt, O0.002–0.063 mm; Sand, 0.063–2 mm. Corg, organic carbon. n.d., not determined.

Bromide was analyzed using a Metrohm Ion
Chromatography system with a conductivity detector
and chemical suppression. Anions were separated on a
Metrosep Dual 2 column (75!4.6 mm; Metrohm,
Switzerland). The flow rate of the mobile phase—a
1.3 mmol/L Na2CO3 and 2 mmol/L NaHCO3 solution—was adjusted at 0.8 mL/min. Sample volume
was 20 mL. The detection limit for Bromide was
0.15 mg/L.
2.3. Model description
Four different numerical model approaches werecompared to describe and interpret the experimental
observations: (i) one-dimensional single-porosity
model (1D-SPM), (ii) one-dimensional dual-porosity
mobile-immobile model (1D-MIM), (iii) two-dimensional single-porosity model (2D-SPM), and (iv) twodimensional dual-porosity mobile-immobile model(2D-MIM). An extended version of HYDRUS-1D(Šimůnek et al., 1998, 2003) was used for 1D-SPM
and 1D-MIM simulations, while similarly modified
HYDRUS-2D (Šimůnek et al., 1999) was used for 2DSPM and 2D-MIM simulations.
In the one-dimensional single-porosity model (1DSPM), vertical variably-saturated water flow and nonreactive solute transport in the soil was describedusing the Richards’ and convection dispersion

equations, respectively:



vq
v
vh
Z
KðhÞ
C1
vt
vz
vz

(1)



vqc
v
vc
vqc
Z
qD
K
vt
vz
vz
vz

(2)

where q is the volumetric water content [L3LK3], his the pressure head [L], z the vertical coordinate(positive upward) [L], t is time [T], K(h) is thehydraulic conductivity function [LT-1], c is thesolute concentration [MLK3], q is the flux density[LT-1], and D is the dispersion coefficient [L2TK1]
defined as:
qD Z lq C D0 qt
tZ

(3a)

q7=3
q2s

(3b)

where l is the dispersivity [L], D0 is the molecular
diffusion coefficient [L2T-1], t is the tortuosity
coefficient [K], evaluated according to Millington
und Quirk (1961), and qs is the saturated watercontent.
The one-dimensional mobile-immobile model(1D-MIM) describes vertical variably-saturated
water flow and non-reactive solute transport in a soil
characterized by mainly two pore systems. The MIM

Table 2
Characteristics of the bromide tracer experiments (22 Nov 1999 to 29 Apr 2000) at plots A, B, and C of the Infeld experimental station
Plot

Drain spacing(m)

Drain depth(m)Plot size(m2)

Start drainage

Total rainfall(cm)

Total drainage(cm)

Applied BrK mass
(kg/Plot; g/m2)

A
B
C

16
18
12

1.28
1.01
0.96

2208
2430
1620

3 Dec 1999
6 Dec 1999
2 Dec 1999

47.7
47.7
47.7

34.2
30.6
31.5

9.50; 4.3
10.43; 4.3
6.95; 4.3

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

395

5Precipitation intensity
Cumulative precipitaion

40

4

Cumulative ETa

30

3

20

2

10

1

0

0
11.21.99 12.11.99 12.31.99

Precipitaion intensity (cm/d)

Cumulative precipitation and ETa (cm)50

1.20.00

2.9.002.29.00

3.20.00

4.9.004.29.00

Date
Fig. 2. Rain intensity, cumulative precipitation and actual evapotranspiration (ETa) at the Infeld experimental site during the tracer experiment(Nov. 21st 1999 to April 29th 2000).

assumes that the soil water content, q, is partitioned
into a mobile, qm, and an immobile, qim, fraction as
follows:q Z qm C qim(4)

This concept was incorporated into the Richards’ (1)and convection-dispersion (2) equations as follows(e.g. Simunek et al., 2003):



vqm
v
vhm
Km ðhÞZ
C 1 KGw
(5a)
vz
vt
vz
vqim
Z Gw
vt

Water exchange, Gw, was assumed to be proportional
to the difference in relative saturations of the mobile
and immobile regions:
Gw Z

(6a)

vqim cim
Z Gs
vt

(6b)

where indices m and im refer to the mobile andimmobile regions, respectively. The rates of water andsolute exchange between mobile and immobile
regions, Gw [TK1] and Gs [ML-3TK1], respectively,
are described by the first-order terms (7) and (9).

(7)

where aw is a first-order water transfer coefficient [T1im], and Sm
e and Se are the dimensionless effective
saturations of the mobile and immobile regions given
by
Sm
e Z

qm Kqm;rqm;s Kqm;r

(8a)

Sim
e Z

qim Kqim;rqim;s Kqim;r

(8b)

(5b)



vqm cm
v
vcm
vqc
q D
Z
K m KGsvz m m vz
vt
vz



vqim
im
Z a w Sme KSevt

where qm,s, qm,r (qim,s, qim,r) are the saturated andresidual mobile (immobile) water contents [L3LK3],
respectively.
The solute exchange term, Gs, includes diffusive(proportional to the concentration gradient between
the mobile and immobile regions) and convective(proportional to water transfer and, depending on
the direction of water transfer, the concentration ineither mobile or immobile region) components as396

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

follows:

2.4. Modeling strategy(

Gs Z as ½cm Kcim  C

Gw cm ; if Gw O 0
Gw cim ; if Gw ! 0(9)where as represents a first-order solute transfer
coefficient [T-1].
In the two-dimensional single-porosity model (2DSPM), variably-saturated water flow and non-reactive
solute transport in a two-dimensional vertical crosssection of the soil is described with the following
forms of Richards’ and convection dispersion
equations:vqvKðhÞ
Z V½KðhÞVh C
vt
vz

(10)



vqc
Z V qDij Vc KVqc
vt

(11)

where P is the vector differential operator (v/vx,v/vz) and Dij is the dispersion coefficient tensor[L2TK1] evaluated according to Bear (1972) as
qDij Z lT jqjdij C ðlL KlT Þ

qj qi
C qD0 tdij
jqj

(12)

where jqj is the absolute value of the water flux [LT-1],
dij is the Kronecker delta function (dijZ1 if iZj, anddijZ0 if isj), lL and lT are the longitudinal andtransverse dispersivities, respectively [L]. Eq. (10)assumes isotropic porous medium.
In the two-dimensional mobile-immobile model(2D-MIM), two-dimensional variably saturated waterflow and physical non-equilibrium tracer transport in
the soil with a secondary porosity was described byextending Eqs. (10) and (11) to yield: vK ðh Þvqm
Z V, Km ðhm ÞVhm C m m KGw
vz
vt

(13)

vqm cmZ V½qm Dij Vcm KVqcm KGs
vt

(14)

We note that in (14) the dispersion coefficient tensor,
Dij, is assigned to the mobile region. Furthermore, therelationships between q, h, and K, or qm, hm and Km
were defined by the van Genuchten-Mualem’sfunctions (Van Genuchten, 1980) for both the singleand dual-porosity (i.e. mobile-immobile) models.

The basic idea was to first evaluate various oneand two-dimensional models with respect to theirability to describe the observations of tile drain
outflow fluxes, cumulative tile drainage water flow,
and BrK concentrations in drainage water at the Bplot, and then to use these model parameters to predict
respective observations collected at the two otherplots (A and C) with different drainage systems. Using
the model approaches as described in Section 2.3, fivedifferent model cases were distinguished as follows:
The 1D-SPM case and the 2D-SPM case usedparameters that were either independently derived or
fixed at estimated (not calibrated) values. The 1D- and2D-MIM cases relied on inverse (1D) or ‘semiinverse’ (2D) estimation of the transfer parameters aw
and as based on the experimental data collected at plot
B. Finally, the 1D-MIMforward case used the transfer
parameters of the 2D-MIM case. Procedures for
estimating the input model parameters are describedin the sections below.
Water flow and BrK transport observations of theplots A and C were then predicted using the 2D-MIM,
the 1D-MIM, and the 1D-MIMforward cases.
2.5. Model parameters
Input parameters for various models are presented
in Table 3. All 1D- and 2D-SPM parameters wereindependently derived or were fixed at estimated (notcalibrated) values. Water retention and hydraulicconductivity parameters were estimated frommeasured laboratory and field data as describedabove in Section 2.1. The dispersivity was set to
0.5 cm, a value in the range of laboratory observations
obtained with undisturbed columns (Meyer-Windel,
1998; Maraqa et al., 1997).
Based on the range of values for the fraction ofmobile water in field soils or large undisturbed soilcolumns as reported in Vanderborght et al. (1997),
and considering the Infeld soil profile as well as using
Eq. (4), mobile and immobile water contents, qm,s andqim,s, for the MIM model were roughly estimated. For
the top (H1) layer with a dense network of subangular
to crumb structure, half of the total saturated watercontent was assigned to both mobile (qm,sZ0.5qs) andimmobile (qim,sZ0.5qs) regions. To represent more

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

distinct interaggregate flow path regions of the (H2,
H3, and H4) subsoil layers with angular to prismatic
structure, a smaller fraction of qs was assumed to be
mobile (i.e. qm,sZ0.2qs), while the remaining fraction
was assumed to be immobile (i.e. qim,sZ0.8qs). The
residual water content of the mobile region, qm,r, wasassumed to be zero.
The three different modeling cases for the mobile
immobile model differed with respect to the treatmentof the solute and water transfer parameters: (i) 2DMIM case: The transfer parameters aw and as of theH1–H4 layers for the 2D-MIM approach were
estimated from the observations of plot B using a
combined inverse analysis and subsequent manual‘fine-tuning’ calibration. To reduce the number offitted parameters from 8 (4 layers) to 6, aw and as
were constrained to similar values in the lower twolayers (i.e. aw(H3)Zaw(H4) and as(H3)Zas(H4)),
which was consistent with the assumption of a single
Ks value and with the observation of relatively similar
structure in the H3 and H4 layers. (ii) 1D-MIM case:
We attempted to fit the 1D-MIM water and solutetransfer parameters, aw and as, similarly as done forthe 2D-MIM, in order to evaluate the capacity of the
1D-MIM for the best possible representation of thedata, and to compare corresponding fitted 1D-MIM
and 2D-MIM parameters. However, the simultaneous
optimization of aw and as in the 1D-MIM failed (i.e.
led to program abortion). Alternatively, we inversely
estimated as while keeping the previously (2D-MIM)
obtained aw values. (iii) 1D-MIMforward: both aw
and as were used as obtained with the 2D-MIM case.
For the inverse parameter identification in the 1Dand 2D-MIM approaches, the Levenberg-Marquardt
optimization algorithm was used to minimize theweighted squared deviations between observations
and simulated values (Simunek et al., 1999):FðbÞ Z

m
X
jZ1

vj

n
X

wi;j ½Oj ðz; ti ÞKEj ðz; ti ; bÞ2

(15)

iZ1

where m is the number of different sets ofmeasurements (cumulative drainage, drainage rates,
and BrK concentrations of the tile-drain effluent), n
represents the number of observations in a particularmeasurement set, Oj(z,ti) are observations at time ti at
location z, Ej(z,ti,b) are the corresponding estimatedspace-time variables for the vector b of optimized

397

parameters and vj and wi,j are weighting factors
associated with a particular measurement set or point,respectively. In the 1D-MIM (2D-MIM) optimization,
the vector b consisted of as (aw and as) for the three(4) soil horizons. Furthermore, we assumed that theweighting coefficients wi,j in (15) were equal to one,which is similar to assuming equal variances of theerrors within a particular measurement set.
2.6. Initial and boundary conditions
Initial soil moistures and pressure heads were notknown, but were fixed to represent a field capacity
type moisture profile based on following considerations. Since drainage flow stopped on July 10, 1999,
and did not start since, the groundwater table wascertainly well below the tile drain level. On the other
hand, repeated rainfall (15.2 cm since mid September1999) prevented the soil from drying out. The initialcondition for water flow was set to a linear pressure
head increase between K600 h Pa at the upperboundary (q(K600 h Pa)Z0.468 for qsZ0.509) andK100 h Pa at a 200-cm depth (q(K100 h Pa)Z0.634
for qsZ0.640). Initial BrK concentrations in the soil
were assumed to be zero.
Daily values of rainfall and actual evapotranspiration were used as the upper boundary condition forwater flow. The flux (Neumann) condition automatically changed into a prescribed ponding head
(Dirichlet) boundary condition when a ponding
depth of 0.2 cm was reached. When the calculated
flux became less than the prescribed flux, theprescribed head (ponding of 0.2 cm) boundary
condition switched back to the prescribed flux. During
ponding, the difference between calculated andprescribed flux was considered to be surface runoff(overland flow). Zero flux boundary conditions were
imposed for both vertical sides and the lower
boundary of the 2D domain. A third-type (Cauchy)boundary condition was used for the solute condition.
The lower boundary condition for 1D simulations
represented groundwater flow towards drains ascalculated according to the Hooghoudt equation(e.g. Simunek et al., 1998). The drain parameters,
i.e. the entrance resistance and the wet perimeter were
set to 5 and 25 cm, respectively; the drainage
geometry (drain spacing and depth) was specified
according to plot characteristics (Table 2), and the398

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

saturated hydraulic conductivity, Ks, in the Hooghoudt equation was calculated from the borehole
method (see above).
In two-dimensional simulations the drain pipe wasrepresented by a system-dependent water flowboundary condition at a single node that switchedbetween zero flux for an unsaturated drain node to azero pressure head upon saturation. To represent theentrance resistance of the drain, the hydraulicconductivities at the nodes adjacent to the drainnode were reduced according to the method of Fipps
and Skaggs (1986). An effective drain diameter of
4 cm was assumed and the standard correction factor(Fipps and Skaggs, 1986) was applied without
additional reduction factor. A zero gradient (Neumann) boundary condition was assumed for solutetransport out of the tile drain.
2.6.1. Numerical Grids
The size of the transport domain for onedimensional simulations was set equal to thecorresponding tile depth (127 cm for plot A, 101 cm
for B, and 96 cm for C). Each domain was subdivided
into 101 nodes at regular vertical increments.
The depths of the transport domain for twodimensional simulations was set equal to the verticaltile location plus 100 cm, in order to minimize theeffect of the lower impermeable boundary. The
domain width was equal to the half-distance between
drains. Resulting transport domains were subdivided
into 46 (A), 51 (B), and 35 (C) vertical columns, and
49 (A) or 44 (B, C) horizontal rows at variable spacingwith higher spatial resolutions at layer interfaces andnear the drain node. The model drain node was locatedat one side of the transport domain at 128 (101, 96) cm
depth for plot A (B, C). The domain was subdivided
into four horizontal layers representing H1–H4 asdescribed in Table 2.

3. Results and discussion
3.1. Experimental data and simulation resultsfor plot B
Outflow from the plot B tile lines started on
December 6th 1999 and the cumulative drainage was
equal to 30.6 cm (64% of the rainfall) over theduration of the experiment (Table 2). Fig. 3 compares
the observed cumulative drainage (Fig. 3a) anddrainage rates (Fig. 3b, c) from plot B with the oneand two-dimensional single-porosity model (SPM)
and mobile-immobile model (MIM) cases. The
measured cumulative drainage (Fig. 3a) was matchedclosely by the 2D-SPM and 2D-MIM cases. Cumulative drainage tended to be overestimated by both
1D-cases in the beginning but was underestimated andmatched by the 1D-SPM and the 1D-MIM, respectively, in the end. Measured drainage rates reachedmaxima of about 1.9 cm/d, but frequently dropped
below 0.5 cm/d. Both 2D models, i.e. 2D-SPM and
2D-MIM, matched the drainage rate peaks above
1 cm/d considerably better than the corresponding 1Dmodels (compare Fig. 3b and c).
Bromide concentration data and model simulation
results are compared in Fig. 4. Measured BrK
concentrations in the tile drainage effluent fluctuatedaround a value of 5 mg/L throughout the experiment(Fig. 4a), without obvious reactions to variations intile-drain discharge. Bromide was found in theeffluent from the start of drainage onward. By
contrast, the 2D-SPM simulated the first appearanceof BrK in the drain to occur 1 month (or 10 cm ofdrainage) later, and the 1D-SPM simulation even gave
a 2.5-month (or about 20 cm of drainage) time-lag
between application and Br K arrival in draindischarge (Fig. 4a). Accordingly, both SPM models
underestimated the BrK mass recovery in drain
effluent at the end of the drainage period: thesimulations suggested recoveries of 2.5% (1D-SPM)
and 8.6% (2D-SPM) as compared to the 25% actually
retrieved in the field experiment (Table 4). Similar to
the measurements, the 2D-SPM simulation showed
short term fluctuations of the Br concentrations
whereas the 1D-SPM simulated BTC was smooth(Fig. 4a).
Fig. 4b compares measured BrK BTCs with thosesimulated using 1D and 2D dual-porosity models. The
2D-MIM matched the measured BrK concentrations
very well in terms of the first appearance of BrK, theconcentration level and the short-term fluctuation
pattern. The 1D-MIMforward case (using 2D-MIM
fitted transfer parameters) predicted a 1 week (or 5 cmof drainage) delayed arrival and slow increase of BrK
concentrations, while it overestimated BrK concentrations in the last third of the experiment. The 1D-

(a)

35

Cumulative drainage (cm)S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

30

399

25
Data

20

1D-SPM

15

2D-SPM

10

1D-MIM

5

2D-MIM0
11/21/99 12/11/99 12/31/99

1/20/00

2/9/002/29/00

3/20/00

4/9/004/29/00

Date

Drainage (cm/d)(b) 2.0
Data

1.5

1D-SPM
2D-SPM

1.0
0.5
0.0
11/21/99 12/11/99 12/31/99

1/20/00

2/9/002/29/00

3/20/00

4/9/004/29/00(c)

2.0

Drainage (cm/d)

Date

1.5

Data
1D-MIM
2D-MIM

1.0

0.5

0.0
11/21/99 12/11/99 12/31/99

1/20/00

2/9/002/29/00

3/20/00

4/9/004/29/00

Date
Fig. 3. Time series of drainage data, and one- and two -dimensional single-porosity model (SPM) and mobile-immobile model (MIM) results forplot B (18 m drain spacing, 1.01 cm drainage depth): a) cumulative drainage, b) drainage data and SPM results, and c) drainage data and MIM
results.

MIM case (with fitted as, solid, thin line in Fig. 4b)gave slightly better simulation results than the 1DMIMforward case.
The effect of using a 2D model vs. a 1D model on
the simulation results is discussed in the following.
The BTC simulated with the two-dimensional models(both SPM and MIM) started earlier and better

represented the observed concentration fluctuations
than the corresponding 1D cases. This can be
explained as follows: in the two-dimensional flow
field, solute is transported along different pathways interms of path lengths, orientations and water
saturations. Directly above the drain, BrK travel
distance is shortest and is oriented vertically towards

400

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

Concentration (mg/l)(a) 10Data

1D-SPM

2D-SPM8
6
4
2
0
11.21.99 12.11.99 12.31.99

1.20.00

2.9.002.29.00

3.20.00

4.9.004.29.00

2.29.00

3.20.00

4.9.004.29.00

Date

Concentration (mg/l)(b) 108

Data

2D-MIM

1D-MIM

1D-MIMforward

64
2
0
11.21.99 12.11.99 12.31.99

1.20.00

2.9.00Date
Fig. 4. Measured bromide concentrations, and one- and two -dimenstion single-porosity model (SPM) and mobile-immobile model (MIM)
results for plot B (18 m drain spacing, 1.01 cm drainage depth): a) data and SPM results, b) data and MIM results.

the drain. At half drain-spacing, the total transportdistance was largest, with mainly vertical orientation
but shorter path length through the unsaturated zone,
followed by a predominantly lateral (curved) flow
paths in groundwater towards the drain. This variationin transport path lengths, directions and velocities inunsaturated zone and groundwater enhanced mixingof BrK with a larger soil water volume and caused
early arrival, dispersion, and fluctuation of simulated
BrK concentrations at the tile drain, similar to the fieldobservations.
The additional effect of considering dual-porosity
vs. equilibrium models becomes obvious by the factthat the simulation results of all MIM cases werecloser to the measured BrK concentrations than both
1D- and 2D-SPM cases (Fig. 4a,b). Accordingly,
MIM-simulated BrK mass losses were much closer tothe observed loss (Table 4). The improvement yieldedwith dual-porosity over single-porosity models to
reproduce the observed rapid solute breakthrough andhigher total solute losses was found earlier and was
explained with preferential flow in structured soils(e.g. Larsson and Jarvis, 1999; Haws et al., 2005).
However, for surface-applied solutes, preferential
flow conditions were typically associated with thehighest concentration peak occurring after the firstmajor rainfall event following solute application (e.g.
Kladivko et al., 1991; Jardine et al., 1993; Lennartz
et al., 1999; Gerke and Köhne, 2004; Gish et al.,2004), a close correspondence between concentration
peaks and outflow peaks (Lennartz et al., 1999; Kunget al., 2000; Jaynes et al., 2001; Kohler et al., 2003)
and high tracer recoveries in tile drain discharge
amounting to at least 40–60% of the applied mass at
the end of the drainage season (e.g. Lennartz et al.,
1999; Köhne and Gerke, 2005). The BTC as observed

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

Relative immobile mass (–)

at the Infeld site lacked most of these preferential flow
characteristics, except for the early solute arrival. In
this context, Skaggs et al. (1998) and Kohler et al.
(2003) stated that at tile drained fields, transport
through the saturated zone may effectively remove
any concentration pattern of preferential transport.
However, this cannot be true for all tile drained fields,
where preferential solute breakthrough in tile drainage
effluent was often observed, as found in many of thestudies mentioned above. Utermann et al. (1990)
showed for simulations of solute transport towards adrain using a bimodal probability density function
approach for the unsaturated zone and a streamline
model for the saturated zone, that the bimodal
concentration distribution simulated in the vadosezone developed lower and wider peaks after passage
through the saturated zone. Yet similar to our results,
in the simulations of Utermann et al. (1990) the early
first appearance of a tracer as caused by preferential
flow was not retarded by transport through thesaturated zone.
Based on our 2D-MIM simulation results, the main
mechanism leading to the weakly preferential BTC
characteristics at the Infeld site is large mobileimmobile solute transfer, which will be discussed in
the following. Fig. 5 shows the dual-porosity
simulated (1D-MIM, 1D-MIMforward, 2D-MIM)
relative BrK mass in the immobile region (as afraction of the applied mass) over time. All MIM
approaches simulated an initially increasing immobilization of surface-applied BrK, showed a slightinflection at the onset of tile outflow, and a continued
rise until December 15th 1999, when approximately

401

67% of the applied BrK mass was immobilized. Only
after heavy rainfall and increased drainage rates in thefirst week of March 2000 did the BrK net mass in theimmobile region gradually decrease, but it remained
above 60% throughout the drainage season. For
comparison, Meyer-Windel (1999) found in BrKtransport experiments with structured soil columns
subject to transient water flow conditions that 20 to
60% of the applied BrK mass was immobilized in thesoil after the BrK concentrations in the effluent
decreased to values below the detection limit.
Obviously, solute release from immobile into mobile
regions is very slow as compared to the solute transferfrom mobile into immobile regions. For a surface
applied tracer, fast net immobilization may be causedby high-rate advective transfer during initial infiltration, where differences in saturation between mobile
and immobile regions are higher. Relatively large
fitted aw values suggest that this was a majormechanism for the two upper layers. Likewise,
diffusive mass transfer is intensified due to highconcentration gradients since the solute is entering thesoil via the mobile region only, whereas the immobile
regions are initially solute free. The backward process
is taking place at a smaller rate, since on one handconcentration gradients are less steep (both regions
carry solute), and on the other hand, with each new
infiltration cycle following redistribution, advective
transport from mobile to immobile regions isreiterated. Our model results suggest that relatively
large solute transfer into the immobile regions
restrained preferential solute transport by reducing
the initial BrK peak and overall BrK loss. Inter-region

0.8
0.6
0.4
0.2
2D-MIM0
11.21.99 12.11.99 12.31.99 1.20.00

1D-MIM

2.9.002.29.00

1D-MIMforward

3.20.00

4.9.004.29.00

Date
Fig. 5. Time series’ of the relative bromide mass in the immobile region (fraction of the applied mass) of the soil profile as calculated using oneand two-dimensional mobile-immobile models (ID-MIM and 2D-MIM, respectively).

Following transport parameters were assumed to be identical in all model simulation: dispersivity: lZ0.5 cm; the effective bromide diffusion coefficient: D0 Z1.797 cm2 dK1aFitted with 2D-MIM and B data.
b
Fitted with 1D-MIM and B data.
c
135–230 for simulations of plot A.

0.020a 0.020b
0.021a 0.023b
0.075a 0.016b
0.075a0.12aO
0.10aO
0.07aO
0.07a
0.190
0.200
0.180
0.210
0.254
0.433
0.460
0.512
0
0
0
0
0.254
0.108
0.115
0.128
4.8
14.9
35.0
35.0
1.59
1.39
1.53
1.67
0.00100.00080.00100.0014
0.19
0.20
0.18
0.21
0–25
25–70
70–135135–200c
H1
H2
H3
H4

0.509
0.541
0.576
0.640

qm,s
(m3 mK3)
Ks(cm dK1)
na(cmK1)
qr(m3 mK3)
qs(m3 mK3)
Depth(cm)
HorizonTable 3
Soil hydraulic parameters for the one- and two-dimensional single- and dual-porosity model simulations

qm,r(m3 mK3)qim,s(m3 mK3)qim,r(m3 mK3)

as(dK1)

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

aw
(dK1)

402

mass transfer and lateral mixing have been identified
previously as a major factor controlling preferential
solute transport (e.g. Flühler et al. 1996). Jaynes et al.
(2001) and Fortin et al. (2002) suggested that mobileimmobile solute transfer may considerably reduceleaching losses.
Here, relatively large values for the first-order rate
coefficient for water transfer, aw, controlling advective solute transfer, were obtained between 0.12 and0.07 dK1 for top and subsoil layers, respectively(Table 3). For simulations of field-scale preferential
flow and solute transport, aw values were found whichwere within similar range (Köhne and Gerke, 2005),
or smaller, i.e. around 0.02 d-1 (Gerke and Köhne,2004) or even below 0.004 dK1 (Haws et al., 2005)
when pronounced preferential flow conditions prevailed. On the other hand, values for the solutetransfer rate coefficient, as, between 0.02 and 0.07 d-1
for top and subsoil layers were within the range of
10K3 to 100 d-1 reported by others for physical nonequilibrium conditions (e.g. Jaynes et al., 1995;
Lennartz and Meyer-Windel, 1995, Köhne andGerke, 2005; Haws et al., 2005). Abbasi et al.
(2004) obtained much higher as values in theirfurrow-irrigation field experiment with Bromide.
Although variably saturated flow conditions prevailed
in the experiment, the 2D MIM model used by Abbasi
et al. (2004) did not consider advective solute transfer,
such that their high as values accounted for advective
and diffusive transfer in a lumped fashion.
The above model simulation results indicate thatthe plateau-like, preferential BrK BTC at the Infeld
tile drained site is a product of both, the twodimensional flow field and non-equilibrium transportcharacteristics with pronounced transfer between
mobile and immobile regions. Additionally to the
spreading of the BTC caused by the 2D flow field, thesolute mass transfer from mobile to immobile regions
has buffered concentrations at the tile drain, since
solute has been taken out of the mobile flow regions.
Although pronounced solute mass transfer governed the tracer transport at Infeld, the simulatedleaching loss during the 1999/2000 drainage season
was still higher under the physical non-equilibriumtransport conditions involving mass transfer (17%,2D-MIM in Table 4) than under equilibrium conditions (8.6%, 2D-SPM in Table 4). To evaluate theeffect of physical non-equilibrium transport and mass

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

403

Table 4
Observed and simulated bromide recovery (kg and (% of applied mass)) in drain outflow
Plot

Data

1D-SPM

2D-SPM

1D-MIM

2D-MIM

A

1.96 (20.6)

n.s.

n.s.

1.19 (12.5)

B

2.59 (24.8)

0.26 (2.5)

0.90 (8.6)

C

1.88 (27.0)

n.s.

n.s.

1.35 (14.3)a0.81 (14.3)b2.22 (21.3)a1.86 (17.8)b2.32 (33.4)a2.71 (31.2)b

1.83 (17.5)
1.40 (20.1)

n.s., not simulated.
a
1D-MIM case.
b
1D-MIMforward case.

transfer on long term solute loss, we compared 2DSPM and 2D-MIM model predictions of solute losses
until the end of the year 2000 (no atmospheric data
were available thereafter). These two model cases
equally simulated a drainage depth of 6.6 cm for
November/December 2000 (i.e. less than 25% of thedrainage of the entire 1999/2000 drainage season).
Fig. 6 shows the simulated relative BrK mass losses
until the end of 2000. About the same amount (8.4%of the applied mass) of BrK as in the entire 1999/2000
drainage season (8.6%) is leached in the following
half-season in the 2D-SPM equilibrium case.
Although the total loss is higher in the 2D-MIM
case assuming physical non-equilibrium and masstransfer, the leaching loss in the 2nd (half) drainage
season is with 6% of the applied mass only one thirdof the 1999/2000 loss and is 2.5% less than that of the
2D-SPM case. This means that conservative soluteleaching under physical non-equilibrium conditions

may in the long run proceed slower than forequilibrium or no-mass transfer conditions.
3.2. Predictions for different drainage systems(A and C plots)
Simulation results for the B plot showed that thedual-porosity concept represented BrK transportobservations at the Infeld site better than the singleporosity concept. As a next step, we applied the 1Dand 2D-MIMs with their plot B parameters to predictwater flow and BrK transport observations at the A
and C plots, which differed from B with regard todrainage depths and spacing. The purpose of carryingout the MIM predictions was to further evaluate theimportance of considering two dimensions in MIM bycomparing how well the 1D- and 2D-MIM approaches
can predict solute transport data at other plots withdifferent drainage systems.

0.30

Relative mass loss (-)

2D-MIM
2D-SPM

0.20

0.10

0.00
11/21/99

02/19/00

05/19/00

08/17/00

11/15/00

Date
Fig. 6. Time series’ of relative cumulative bromide mass loss (fraction of applied mass) as simulated with two-dimensional mobile-immobile(2D-MIM) and single porosity (2D-SPM) models.

404

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

Cumulative drainage (cm)(a) 40

30

20

10

0
11.21.99 12.11.99 12.31.99

1.20.00

2.9.00Data

2D-MIM

1D-MIM

1D-MIMforward

2.29.00

3.20.00

4.9.004.29.00

2.29.00

3.20.00

4.9.004.29.00

Date

Concentration (mg/l)(b) 108

Data

2D-MIM

1D-MIM

1D-MIMforward

64
2
0
11.21.99 12.11.99

12.31.99

1.20.00

2.9.00Date

Fig. 7. Measured cumulative drainage (a) and bromide concentrations (b) compared with one- and two-dimensional mobile-immobile model(MIM) predictions for the plot A (16 m drain spacing, 128 cm drainage depth).

The A plot had a tile drainage system characterized
by a considerably larger depth (128 cm) and somewhat narrower spacing (16 m) than the B plot (101 cm
depth, 18 m spacing) (see Table 2). Fig. 7 compares
the observed tile drain outflow and BrK concentrations with one- and two-dimensional MIM predictions for the A plot. Compared with B, the drainage at
A started 3 days earlier and yielded a 3.6 cm higher
total value at the end of the experiment, which isconsistent with the deeper tile-drain level (Fig. 7a,
Table 2). The cumulative drainage was matched wellby the 2D-MIM and was most of the time slightly
overestimated by the 1D-MIM cases, similarly as for
the B plot (compare Figs. 3a and 7a). As the soil
hydraulic parameters were the same for both 1D-MIM
and 1D-MIMforward, their simulated drain outflow
curves were identical.

Measured BrK concentrations quickly reached an
initial level of 5–6 mg/L during December 1999
(Fig. 7b). After a month of missing data, observed
BrK concentrations fluctuated around a lower level of
3–4 mg/L, with a dip in concentrations in March
2000. Neither of the MIM predictions matched theBrK concentration data very well (Fig. 7b). Simulated BrK mass recoveries in drain effluent underestimated the measurement of the A experiment(Table 4). Still, the 2D-MIM was the best in
predicting the onset of the observed BrK concentrations and the later plateau-like shape of the BrK
BTC. The 1D-MIMforward (using water and solutetransfer coefficients obtained for 2D-MIM for plot B)
initially predicted a late appearance of BrK indrainage and significantly underestimated BrK concentrations. During the last third of the experiment,

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

405

Cumulative drainage (cm)(a) 4030
20
10
011.21.99 12.11.99 12.31.99 1.20.00

2.9.00Data

2D-MIM

1D-MIM

1D-MIMforward

2.29.00

3.20.00

4.9.004.29.00

3.20.00

4.9.004.29.00

Date

Concentration (mg/l)(b) 15
Data

2D-MIM

1D-MIM

1D-MIMforward

10

5

0
11.21.99 12.11.99 12.31.99 1.20.00

2.9.002.29.00

Date
Fig. 8. Measured cumulative drainage (a) and bromide concentrations (b) compared with one- and two -dimensional mobile-immobile mode(MIM) predictions for the plot C (12 m drain spacing, 96 cm drainage depth).

the 1D-MIMforward predicted BrK concentrations
exceeding the measured values by about 2 mg/L
(Fig. 7b). The BrK BTC predicted with the 1D-MIM(with the solute transfer coefficients fitted to the Bdata) somewhat better represented the initial BrK
concentration increase, but predicted a peak of 6 mg/
L in early March 2000 that considerably exceeded thedata (Fig. 6b).
The C plot had the narrowest tile drain spacing(12 m) among the three plots (Table 2). C also had theshallowest tile drains (96 cm), but with tile drain
depth only slightly smaller than at B (101 cm). Fromthe beginning, the measured drainage curve wasslightly overestimated by the 2D-MIM, and more soby the 1D-MIM (Fig. 8a). Plot C, similarly as A, had aplateau-like measured BrK BTC with gradually
decreasing concentrations (Fig. 8b). The observed
BrK concentrations were closely approximated by the
2D-MIM prediction. Although the 1D-MIMforward
predicted a late start of the BrK BTC, it matched very

well the concentration level of 7 mg/L measured at the
beginning of the BrK BTC. Starting from January
2000, the 1D-MIMforward overestimated themeasured BrK concentrations. The 1D-MIM predicted an early peak of almost 15 mg/L and asecondary peak that again overestimated Br K
measurements.
Neither the 1D-MIM nor the 1D-MIMforward case(using transfer parameters obtained from plot B data)
adequately predicted observations of plots A and C
having different drainage geometries than B. The
failure of the 1D-MIM prediction became particularly
obvious by the unrealistically high Br K peaksimulated for the C plot. This failure revealed that acalibrated 1D-MIM model cannot be used to predict
solute transport in tile-drained field plots, since it does
not consider the tile-drain induced flow field. When
we fitted as on the B data (1D-MIM case), theresulting as values for two upper layers with mainly
vertical flow hardly changed (H1: asZ0.020 dK1, H2:

406

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

asZ0.023 d-1) and hence appeared to be close to the
optimal values; only as of the H3 horizon was reduced
from 0.075 to 0.016 dK1 (Table 3) to compensate forthe 2D-flow field effects on the BrK BTC. Recently,Köhne and Gerke (2005) stated the relevance of the2D transport scheme at tile drained fields with nonequilibrium flow. Moreover, Abbaspour et al. (2001)
found that multiple parameter fitting does not yieldsatisfactory transport predictions unless all relevant
processes are described. In our study, the solutetransfer coefficients in the 1D-MIM approach thatwere fitted to B data must be viewed as effective
parameters that are specific only to the B plot. The
2D-MIM gave a more consistent description ofobserved BrK transport (and water flow) at the 3plots with different tile-drainage geometries. Hence,parameters of the 2D-MIM that were calibrated
against observations from one plot could be used toadequately predict physical non-equilibrium solutetransport in plots with different tile drain
configurations.

4. Conclusions
The observed rapid bromide breakthrough at the
tile-drained Infeld site without high initial peak wascharacterized as weakly preferential. The systematiccomparison between 1D and 2D single- and dualporosity model simulations suggested two major
underlying transport processes. First, the tile-drain
induced flow field with variable lengths of unsaturated
and saturated transport pathways from the surface
towards the tile drain dispersed the BrK tracer
concentrations in drain effluent. Second, physicalnon-equilibrium gave rise to rapid solute leaching, butlarge solute mass transfer into immobile regions
prevented the high early solute losses typical for
preferential flow.
The fact that 2D-MIM simulations were moreaccurate for different tile drain configurations than 1DMIM simulations imply that the two-dimensional flow
should be considered when aiming at a mechanistic
description of physical non-equilibrium solute transport in tile-drained fields. The 2D-MIM simulations
explained the observed BrK concentrations level withover 60% of the applied BrK mass being initiallytransferred into immobile soil regions, and stored
there throughout the drainage season.
The MIM (compared with SPM) analyses showedthat even when physical non-equilibrium transport
increases early mass loss, in the long run mass transfer
may slow down solute leaching. In case of degradable
chemicals, however, immobilization of a solute as
caused by mass transfer and slow release may leavetime for efficient decay of a chemical.
Acknowledgements
This work was financed by Deutsche Forschungsgemeinschaft grant no. LE 945 4-1 and promoted by
the Dept. of Environmental Sciences, University ofCalifornia in Riverside. The authors greatly thank the
Landwirtschaftskammer Weser-Ems for providing the
instrumented experimental field site, preliminary data,
and technical support. Furthermore, we acknowledge
Stefan Arndt for his contribution in the experimentation and data compilation. The authors sincerely
appreciate the support of Dr. Schäfer from Bodentechnisches Institut Bremen des Niedersächsisches
Landesamt für Bodenforschung (BTI/NLFB) for thisresearch project.
References
Abbasi, F., Feyen, J., van Genuchten, M.Th., 2004. Two dimensional
simulations of water flow and solute transport below furrows:
model calibration and validation. J. Hydrol. 290, 63–79.
Abbaspour, K.C., Kohler, A., Simunek, J., Fritsch, M., Schulin, R.,
2001. Application of a two-dimensional model to simulate flow
and transport in a macroporous agricultural soil with tile drains.
Eur. J. Soil Sci. 52, 433–447.
Andreu, L., Moreno, F., Jarvis, N.J., Vachaud, G., 1994.
Application of the model MACRO to water movement andsalt leaching in drained and irrigated marsh soils, Marismas,
Spain. Agric. Water Manage. 25, 71–88.
Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, NewYork.
Bouma, J., 1990. Using morphometric expressions for macropores toimprove soil physical analyses of field soils. Geoderma 46, 3–11.
Diez, T., Weigelt, H., 1987. Böden unter landwirtschaftlicher
Nutzung. BLV Verlagsgesellschaft, München pp. 108–109.
DVWK, 1996. Ermittlung der Verdunstung von Land- undWasserflächen. Merkblätter zur Wasserwirtschaft, Heft 238.
Wirtschafts- und Verlagsgesellschaft Gas und Wasser mbH,
Bonn.
Fipps, G., Skaggs, R.W., 1986. Drains as a boundary condition in
finite elements. Water Resour. Res. 22, 1613–1621.

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408
Fortin, J., Gagnon-Bertrand, E., Vézina, L., Rompé, M., 2002.
Preferential bromide and pesticide movement to tile drains under
different cropping practices. J. Environ. Qual. 31, 1940–1952.
Gerke, H.H., Köhne, J.M., 2004. Dual-permeability modeling ofpreferential bromide leaching from a tile drained glacial tillagricultural field. J. Hydrol. 289, 239–257.
Gerke, H.H., Van Genuchten, M.Th., 1993. A dual-porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29, 305–319.
Gish, T.J., Kung, K.-J.S., Perry, D.C., Posner, J., Bubenzer, G.,
Helling, C.S., Kladivko, E.T., Steenhuis, T.S., 2004. Impact ofpreferential flow at varying irrigation rates by quantifying massflux. J. Environ. Qual. 33, 1033–1040.
Haws, N.W., Rao P.S.C., Šimůnek, J., Poyer, I.C., 2005. Single-porosityand dual-porosity modeling of water flow and solute transport in
subsurface-drained field using effective field-scale parameters. J.
Hydrol. (in press) doi: 10.10161j.hydrol.2005.03.035.
Jardine, P.M., Jacobs, K.G., Wilson, G.V., 1993. Unsaturated
transport processes in undisturbed heterogeneous porous media:
I. Inorganic contaminants. Soil Sci. Soc. Am. J. 57, 945–953.
Jarvis, N.J., 1994. The MACRO Model. (Vers. 3.1.) Technicaldescription and sample simulations. Reports and Dissertations
19. Department of Soil Science. Swedish Univerity ofAgricultural Science, Uppsala, Sweden, p. 51.
Jarvis, N.J., 1998. Modeling the impact of of preferential flow on
nonpoint source pollution. In: Selim, H.M., Ma, L. (Eds.),
Physical Non-equilibrium in Soils: Modeling and Application.
Ann Arbor Press, Chelsea, MI, pp. 195–221.
Jaynes, D.B., Logsdon, S.D., Horton, R., 1995. Field method formeasuring mobile/immobile water content and solute transfer
rate coefficient. Soil Sci. Soc. Am. J. 59, 352–356.
Jaynes, D.B., Ahmed, S.I., Kung, K.-J.S., Kanwar, R.S., 2001.
Temporal Dynamics of preferential flow to subsurface drain.
Soil Sci. Soc. Am. J. 65, 1368–1376.
Kladivko, E.T., van Scoyoc, G.E., Monke, E.J., Oats, K.M., Pask,
W., 1991. Pesticide and Nutrient movement into subsurface tile
drains on a silt loam soil in Indiana. J. Environ. Qual. 20, 264–
270.
Kohler, A., Abbaspour, K.C., Fritsch, M., van Genuchten, M.Th.,
Schulin, R., 2001. Simulating unsaturated flow and transport ina 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 export to subsurface
drains by preferential flow. Vadoze Zone J. 2, 68–75.
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.
Kung, K.-J.S., Steenhuis, T.S., Kladivko, E.T., Gish, T.J., Bubenzer,G., Helling, C.S., 2000. Impact of preferential flow on thetransport of adsorbing and non-adsorbing tracers. Soil Sci. Soc.
Am. J. 64, 1290–1296.
Kung, K.-J.S., Kladivko, E.T., Gish, T.J., Steenhuis, T.S., Bubenzer,G., Helling, C.S., 2000. Quantifying preferential flow bybreakthrough of sequentially-applied tracers: silt loam soil.
Soil Sci. Soc. Am. J. 64, 1296–1304.

407

Larsson, M.H., Jarvis, N.J., 1999. A dual-porosity model to quantify
macropore flow effects on nitrate leaching. J. Environ. Qual. 28,
1298–1307.
Larsson, M.H., Jarvis, N.J., Torstensson, G., Kasteel, R., 1999.
Quantifying the impact of preferential flow on solute transport
to tile drains in a sandy field soil. J. Hydrol. 215, 116–134.
Lennartz, B., Meyer-Windel, S., 1995. The role of immobile waterin unsaturated substrates. Hydrogéologie 4, 75–83.
Lennartz, B., Michaelsen, J., Wichtmann, W., Widmoser, P., 1999.
Time variance analysis of preferential solute movement in a tile
drained field site. Soil Sci. Soc. Am. J. 63, 39–47.
Maraqua, M.A., Wallace, R.B., Voice, T.C., 1997. Effects of degreeof water saturation on dispersivity and immobile water in sandy
soil columns. J. Contam. Hydrol. 25, 199–218.
Meyer-Windel, S., 1998. Transport processes of bromide andherbicides under steady-state and transient flow conditions,
Schriftenreihe Institut für Wasserwirtschaft u. Landschaftsökologie der Christian-Albrechts Universität Kiel, Heft 28 p. 113.
Mohanty, B.P., Bowman, J.M., Hendrickx, J.M.H., Šimůnek, J., vanGenuchten, M.Th., 1998a. 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., Skaggs, T.H., van Genuchten, M.Th., 1998b. Impactof saturated hydraulic conductivity on the prediction of tile flow.
Soil Sci. Soc. Am. J. 62, 1522–1529.
Neuhaus, H., 1983. Nitratmobilität bei Grünlandnutzung auf
gedränten tonreichen Marschböden. Zeitschrift für Kulturtechnik und Flurbereinigung 24, 347–357.
Richard, T.L., Steenhuis, T.S., 1988. Tile drain sampling ofpreferential flow on a field scale. J. Contam. Hydrol. 3, 307–325.
Richter, D., 1995. Ergebnisse methodischer Untersuchungen zur
Korrektur des systematischen Meßfehlers des Hellmann-Niederschlagsmessers, Berichte des Deutschen Wetterdienstes, p. 194.
Šimůnek, J., de Vos, J.A., 1999. Inverse optimization, calibrationand validation of simulation models at the field scale. In: Feyen,J., Wiyo, K. (Eds.), Modelling Transport Processes in Soils at
Various Scales in Time and Space Proc. of Intern. Workshop ofEurop. Agric. Eng. Field of Interest on Soil and Water, Leuven,
Belgium, 24.-26. Wageningen Pers, Wageningen, The Netherlands, pp. 431–445.
Šimůnek, J., Šejna, M., van Genuchten, M. Th., 1998. The
HYDRUS-1D software package for simulating one-dimensional
movement of water, heat, and multiply solutes in variably
saturated media. Version 2.0, IGWMC-TPS-70, International
Ground Water Modelling Center, Colorado School of Mines,
Golden, Colorado, p. 251.
Šimůnek, J., Šejna, M., van Genuchten, M. Th., 1999. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat, and multiply solutes in variably
saturated media. Version 2.0, IGWMC-TPS-53, International
Ground Water Modelling Center, Colorado School of Mines,
Golden, Colorado, p.202.
Šimůnek, J., Jarvis, N.J., Van Genuchten, M.T.h., 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.

408

S. Köhne et al. / Journal of Hydrology 321 (2006) 390–408

Skaggs, R.W., Skaggs, R.W., 1980. Drainmod reference report,
Methods For Design And Evaluation Of Drainage-water Management Systems For Soils With High Water Tables. USDA, Soil
Conservation Service, Fort Worth, TX.
Skaggs, T.H., Kabala, Z.J., Jury, W.A., 1998. Deconvolution of a
nonparametric transfer function for solute transport in soils.
J. Hydrol. 207, 170–178.
Southwick, L.M., Willis, G.H., Selim, H.M., 1992. Leaching ofatrazine from sugarcane in southern Louisiana. J. Agric. Food
Chem. 40, 1264–1268.
Twitty, K., Rice, J., 2001. National Engineering Handbook Part 624,
Chaper 10, Water Table Control, USDA, NRCS, 104 p.
Utermann, J., Kladivko, E.J., Jury, W.A., 1990. Evaluating pesticidemigration in tile-drained soils with a transfer function model.
J. Environ. Qual. 19, 707–714.

Vanderborght, J., Mallants, D., Vanclooster, M., Feyen, J., 1997.
Parameter uncertainty in the mobile-immobile solute transportmodel. J. Hydrol. 190, 75–101.
Van Genuchten, M.Th., 1980. A closed-form equation forpredicting 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.
Wendling, U., Schellin, H.G., Thomä, G., 1991. Bereitstellung vontäglichen Informationen zum Wasserhaushalt des Bodens für
die Zwecke der agrar-meteorologischen Beratung. Z. Meteorologie 41, 468–475.