Journal of Contaminant Hydrology 85 (2006) 1 – 32
www.elsevier.com/locate/jconhyd

Multi-process herbicide transport in structured soilcolumns: Experiments and model analysis
J. Maximilian Köhne a,*, Sigrid Köhne a, Jirka Šimlnek baInstitute for Land Use, Faculty for Agricultural and Environmental Sciences, University of Rostock,
Justus-von-Liebig-Weg 6, 18044 Rostock, Germanyb
University of California, Riverside, Department of Environmental Sciences, 900 University Avenue, A135 Bourns Hall,Riverside, CA 92521, United States
Received 16 June 2005; received in revised form 12 December 2005; accepted 3 January 2006
Available online 21 February 2006

Abstract
Model predictions of pesticide transport in structured soils are complicated by multiple processes
acting concurrently. In this study, the hydraulic, physical, and chemical nonequilibrium (HNE, PNE, and
CNE, respectively) processes governing herbicide transport under variably saturated flow conditions werestudied. Bromide (Br), isoproturon (IPU, 3-(4-isoprpylphenyl)-1,1-dimethylurea) and terbuthylazine(TER, N 2-tert-butyl-6-chloro-N 4-ethyl-1,3,5-triazine-2,4-diamine) were applied to two soil columns. An
aggregated Ap soil column and a macroporous, aggregated Ah soil column were irrigated at a rate of
1 cm h 1 for 3 h. Two more irrigations at the same rate and duration followed in weekly intervals.
Nonlinear (Freundlich) equilibrium and two-site kinetic sorption parameters were determined for IPU and
TER using batch experiments. The observed water flow and Br transport were inversely simulated using
mobile–immobile (MIM), dual-permeability (DPM), and combined triple-porosity (DP-MIM) numericalmodels implemented in HYDRUS-1D, with improving correspondence between empirical data andmodel results. Using the estimated HNE and PNE parameters together with batch-test derived equilibrium
sorption parameters, the preferential breakthrough of the weakly adsorbed IPU in the Ah soil could bereasonably well predicted with the DPM approach, whereas leaching of the strongly adsorbed TER was
predicted less well. The transport of IPU and TER through the aggregated Ap soil could be described
consistently only when HNE, PNE, and CNE were simultaneously accounted for using the DPM. Inverse
parameter estimation suggested that two-site kinetic sorption in inter-aggregate flow paths was reduced ascompared to within aggregates, and that large values for the first-order degradation rate were an artifact
caused by irreversible sorption. Overall, our results should be helpful to enhance the understanding and* Corresponding author.
E-mail addresses: max.koehne@uni-rostock.de (J.M. Köhne), sigrid.koehne@uni-rostock.de (S. Köhne),jiri.simunek@ucr.edu (J. Šimlnek).
0169-7722/$ - see front matter D 2006 Elsevier B.V. All rights reserved.
doi:10.1016/j.jconhyd.2006.01.001

2

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

modeling of multi-process pesticide transport through structured soils during variably saturated waterflow.
D 2006 Elsevier B.V. All rights reserved.
Keywords: Pesticides; Preferential flow; Reactive transport; Physical nonequilibrium; Chemical nonequilibrium;
Modeling; Dual-permeability; Triple-porosity; HYDRUS-1D

Abbreviations and notations
AcronymsBr dissolved bromide
BTC breakthrough curve
CNE chemical nonequilibrium
DPM dual-permeability model
DP-MIM dual-permeability mobile–immobile model
DP-KM dual-permeability kinetic modelE
statistical measure: modified coefficient of agreement (Eq. (23))
HNE hydraulic nonequilibrium
IPU isoproturon
MAE statistical measure: mean absolute deviation (Eq. (24))
MIM mobile–immobile model
PFP preferential flow path
PNE physical nonequilibrium
TER terbuthylazine
Symbolsc (M L 3) solute concentration
D (L2 T 1) dispersion coefficientD 0 (L2 T 1) diffusion coefficientf f () fraction of equilibrium sorption sites in the fracture pore region of the DP-KM,
Eqs. (19a) and (20a), or factor to increase or decrease sorption in the fracture poreregion of the DPM, Eq. (15a)f m () fraction of equilibrium sorption sites in the matrix of the DP-KM
f mo () fraction of sorption sites in contact with the mobile region
1  f mo () fraction of sorption sites in contact with the immobile regionh (L) pressure head
K (L T 1) hydraulic conductivity
K̄ m (L T 1) effective hydraulic conductivity of the matrix, Eq. (8)K s (L T 1) saturated hydraulic conductivity
k fr (Mb L3b ) Freundlich distribution coefficient
M () number of data-setsn () van Genuchten’s (1980) soil hydraulic parameter
N () number of observations
O () experimental observationq (L T 1) volumetric water flux density

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

3

s (M M 1) adsorbed solute concentrations e (M M 1) adsorbed solute concentration associated with equilibrium sorption site
fractions k (M M 1) adsorbed concentration of solute associated with kinetic sorption site fractiont (T) timev j () weighting factor for a particular data-set in the objective function, Eq. (21)w f () volumetric fraction of the fracture pore system
w i,j () weighting factor for individual data-points in the objective functionz (L) vertical spatial coordinatea (L 1) van Genuchten’s (1980) soil hydraulic parameterb () exponent of Freundlich sorption isothermk (L) dispersivityl (T 1) first-order herbicide degradation constant in the liquid phaseh (L3 L 3) water content
h r (L3 L 3) residual water content parameter in the van Genuchten’s (1980) modelh s (L3 L 3) saturated water content parameter in the van Genuchten’s (1980) modelq (M L 3) bulk densityr () variances () exponent in the van Genuchten hydraulic conductivity functionx h (L 1 T 1) first-order water transfer rate coefficientx p (T 1) first-order solute transfer rate coefficientx c (T 1) first-order kinetic sorption rate coefficientC h (T 1) first-order water transfer rate
C p (M L 3T 1) first-order solute transfer rate
U c () objective function
Subscriptscchemical nonequilibrium (in the kinetic sorption rate coefficient, x c)eequilibrium (in s e, Eqs. (18a,b))ffracture pore regionhhydraulic (in x h, Eqs. (5) and (7))
inumbering index, 1 V i V n, Eq. (21).
imimmobile regionjnumbering index, 1 V j V m, Eq. (21).
kkinetic (in s k, Eqs. (18a,b))m
matrixmo mobile regionpphysical nonequilibrium (in the solute transfer rate coefficient, x p)

1. Introduction
Pesticides are applied on arable land to prevent crop disease and to secure the crop
production. Pesticide findings in surface and ground waters have raised public concern. Diffuse
herbicide losses from agricultural fields generally dominate detections in ground and surface
waters as compared to point sources such as accidental spills and leaks (USEPA, 1993).

4

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

Fine-textured soils containing clay minerals and organic matter may act as sorption filter
against pesticide leaching. However, this filter is often perforated by soil structural elements
such as biopores (earthworm burrows, decayed plant root channels) or mechanical shrinkage
patterns (cracks or fractures). Through these preferential flow paths (PFP), up to a few percentof the surface-applied pesticide may be channeled below the root zone, particularly during
rainstorm events soon after pesticide application (e.g., Flury, 1996; Kladivko et al., 2001).
Pesticide fate and transport in structured soils was studied in the past using experimental andmodeling approaches. In batch experiments, pesticide sorption was investigated in terms ofsorption isotherms, ad- and desorption hysteresis, and kinetic sorption (e.g., Berglof et al.,
2003; Smith et al., 2003; Cooke et al., 2004). Kinetic sorption is a chemical nonequilibrium(CNE) process, where either the sorption reaction at the soil–solution interface or the diffusion
to the sorption sites is slow. Transport experiments with structured soil columns demonstrated
effects of preferential water flow on pesticide transport (e.g., Czapar et al., 1992; Meyer-Windel
et al., 1999). Preferential flow and solute transport are characterized by hydraulicnonequilibrium (HNE) and physical nonequilibrium (PNE), i.e., by distinct local differences
between flow velocities, water contents, pressure heads for HNE (Köhne and Mohanty, 2005)
and by different advective solute transport velocities and concentrations for PNE (e.g.,
Vanderborght et al., 2002) across soil regions associated with different mobility. Furthermore,
pesticide movement was also strongly affected by CNE (e.g., James and Rubin, 1979; Valocchi,
1985). To complicate the matter, PNE and CNE during transport may produce somewhat
similar effects of limited pesticide sorption, and may act concurrently (e.g., Gamerdinger et al.,
1990; Selim et al., 1999).
Models that do not account for PNE cannot accurately describe observations of pesticideleaching through structured soils in the presence of preferential flow (Pennell et al., 1990;
Parrish et al., 1992; Veeh et al., 1994). Accordingly, one line of solute transport model
development focused on the description of PNE based on the two-region concept. Common
two-region approaches include, among others, i) dual-porosity or mobile–immobile models(MIM) that assume either steady-state (e.g., van Genuchten and Wierenga, 1976) or transient(Šimlnek et al., 2003) water flow, and ii) two-flow domain or dual-permeability models (DPM)
that also assume either steady-state (Skopp et al., 1981; Ma and Selim, 1995) or transient (e.g.,
Gerke and van Genuchten, 1993; Jarvis, 1994) flow. These two-region approaches generally
assume instantaneous reversible chemical sorption and identical sorption isotherms in bothregions.
Another line of contaminant transport model development concentrated on CNE presumingthat sorption to the soil solid phase is kinetic only (Valocchi, 1985), or supposing that sorption
is characterized by simultaneous equilibrium and kinetic processes, where the sorbent in thesoil can be divided into equilibrium (type-1) and kinetic (type-2) sites (e.g., van Genuchten and
Wierenga, 1976; Cameron and Klute, 1977). Such two-site sorption models were shown toapply for many pesticides (e.g., Brusseau and Rao, 1989; Ball and Roberts, 1991).
Finally, coupled PNE – CNE transport models were developed. For instance, Brusseau et al.
(1989) proposed a model where PNE (two-region MIM) and CNE (two-site sorption model)transport descriptions were combined. Wagenet and Chen (1998) combined models of singlerate or gamma-distributed multi-rate sorption with MIM. The resulting PNE – CNE models gave
improved descriptions of atrazine breakthrough curves (BTCs) from column experiments,
particularly for the multi-rate sorption model (Wagenet and Chen, 1998). Ma and Selim (1998)
combined kinetic one-site or multiple-site adsorption models with MIM and obtained better
agreement with measured metolachlor and chromium breakthrough curves. Selim et al. (1999)

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

5linked a second-order two-site model with a MIM and yielded improved description of
metolachlor transport in soil columns. Most recently, DPM were coupled with the two-sitekinetic sorption model (Ray et al., 2004; Pot et al., 2005), which resulted into excellent
simulations of IPU and metribuzin transport in structured soil columns (Pot et al., 2005).
However, coupled PNE – CNE approaches are generally limited to steady-state water flow andcannot account for natural transient flow conditions. On the other hand, transient-flow DPM
with their simple description of equilibrium sorption may not realistically describe pesticidetransport in natural soils (Dubus and Brown, 2002). Contrary to model assumptions, differing
sorption properties were measured for PFPs and soil matrix. Organic linings of earthworm
burrows were found to increase the total organic carbon content and Freundlich distribution
coefficient, and degradation constant, as compared to the adjacent bulk soil matrix (Stehouwer
et al., 1993, 1994; Mallawatantri et al., 1996). Fracture coatings may not only result in different
hydraulic conductivity (Gerke and Köhne, 2002) and diffusion properties (Chen et al., 1997;
Köhne et al., 2002a), but may also alter herbicide sorption (Celis et al., 1997).
The overall objective of this work was to identify and simulate the processes governing theleaching of two common herbicides with contrasting sorption properties, isoproturon andterbuthylazine, in aggregated and macroporous soil columns during transient flow. For this end,
the MIM, DPM, and triple-porosity (DP-MIM) approaches within modified HYDRUS-1D(Šimlnek et al., 2003; Pot et al., 2005), were here further extended for variably saturated,transient flow conditions. Specific objectives were to analyze (i) if model parameters obtained
from inverse analysis of transient water flow and Br transport and from equilibrium or kinetic
batch herbicide sorption experiments can be applied in the forward simulation of herbicidetransport during variably saturated flow, (ii) if the experimental data were sufficient to identify
the governing HNE–PNE–CNE processes and corresponding model parameters and (iii) if there
were differences in herbicide sorption between the PFP and the matrix domains.
2. Material and methods
Part of the experimental data, i.e. Br and isoproturon transport in the Ap and Ah columns,
was published by the second author (Meyer-Windel et al., 1999). Since the terbuthylazine
transport data and the equilibrium and kinetic sorption batch experiments for terbuthylazine andisoproturon have not yet been published, the experimental procedure is described below.
2.1. Soil description
Undisturbed soil columns (14.7 cm diameter, and 15 cm long) were excavated manually froma loamy (Ap) and a sandy (Ah) topsoil (Meyer-Windel et al., 1999). Additionally, bulk soil
samples were taken from the same horizons. The soil from the Ap horizon was ploughed every
fall and had a subangular aggregate structure. The sandy Ah soil was not tilled for four years
prior to sampling and the Ah soil column contained a continuous earthworm burrow. Table 1
gives an overview of selected physical and chemical characteristics of the two soil horizons.
2.2. Batch studies: equilibrium and kinetic herbicide sorption
Equilibrium and time dependent (kinetic) sorption of the phenylurea herbicide isoproturon(IPU, 3-(4-isopropylphenyl)-1,1-dimethylurea) and the triazine herbicide terbuthylazine (TER,
N 2-tert-butyl-6-chloro-N 4-ethyl-1,3,5-triazine-2,4-diamine) as commercial products RArelon

6

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

Table 1
Selected physical and chemical properties of the Ap and Ah horizons
Horizon

Depth
(m)

Clay
(%)

Silt
(%)

Sand
(%)

q
( g cm 3 )Structure
–

C org
(%)

pH
–

Ap
Ah

0 – 0.3
0 – 0.25

2.3
9.2

43.8
19

54
71.7

1.37
1.46

Subangular to crumb
Crumb, earthworm burrows

1.83
1.55

6.7
6.1

Equivalent particle diameter (in 10 6 m): Clay, b2; Silt, 2–63; Sand, N63 – 2000. The pH was measured in 0.01 M CaCl2
solution. q = bulk density, C org = organic carbon.

and RGardoprim were tested on soil materials from the Ah and Ap horizons. The soil material
was air dried and passed through a 2-mm sieve. Adsorption rates were quantified using 0.01 M
CaCl2 solutions that carried 0.2 mg L 1 of either TER or IPU. Duplicate samples of 5 g sterilesoil and 25 mL herbicide solution were shaken overhead in inert OAK-Ridge centrifuge tubes
for 0.016, 0.33, 0.66, 1.33, 3, 9, 12, 24, 48, and 72 h. After shaking, the tubes were rotated in acentrifuge at 4600 g for 30 min whereupon 20 mL of the supernatant solution were removed
and stored at 4 8C until chemical analysis. The difference between concentrations before andafter shaking was attributed to sorption.
Additionally, we followed Streck et al. (1995) who showed that for organic pesticides
undergoing a two-site (fast and slow) sorption process, both equilibrium sorption constants andsorption rates can be derived by fitting the two-site sorption model to data of hysteretic ad- anddesorption isotherms. Ad- and desorption equilibrium isotherms were measured using a 0.01 M
CaCl2 solution containing 0.05 (0.05), 0.1 (0.1), 0.5 (0.5), 2 (2) and 10 (5) mg L 1 IPU (TER).
To obtain an adsorption isotherm, duplicate samples of 5 g soil and 25 mL herbicide solution
were equilibrated by overhead shaking for 24 h. Again tubes were rotated and solution waswithdrawn and stored as described above. Then desorption isotherms were established by adding
20 mL of herbicide-free CaCl2 solution to the soil-solution remainder of the adsorption
procedure. Again tubes were shaken (24 h), rotated, and solution was removed and stored as
described previously. The desorption procedure was repeated three more times. Thus, theadsorption–desorption procedure yielded 1 adsorption and 4 desorption isotherms per herbicidesoil combination. The adsorption isotherms were calculated from sorbed concentrations (s) as
the difference in the herbicide concentrations before (c 0) and after equilibration (c e), multipliedwith the solution / soil ratio (s = 5(s 0  c e)). For the ith-desorption isotherm of a herbicide-soil
combination (i = 1, 2, 3, 4; i = 0 is for adsorption), the sorbed concentrations (s i ) were calculatedusing the corresponding measured equilibrium concentration (c e,i ), and the previously sorbed(s i1) and dissolved (c e,i1) concentrations, as s i = s i1  [5 (c e,i  c e, i1)].
The two-site sorption concept (e.g., van Genuchten and Wagenet, 1989) was assumed torepresent the solute concentration in the solid phase, s, to be associated with two different kindsof sorption sitess ¼ se þ sk

ð1Þ

1

where s e (MM ) denotes the adsorbed concentration of a solute in instantaneous equilibriumwith the fraction of equilibrium sorption sites, and s k (MM 1) correspondingly symbolizes theremaining concentration of a solute attached to the solid phase by a slower, time-dependent
process. The first-order rate equation for kinetic solute sorption at the type-2 sites is written asBsk¼ xc ð1  f Þkfr cb  sk
Bt

ð2Þ

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

7

where x c is a first-order sorption rate constant [T 1], f is the fraction of equilibrium sorption
sites [], k fr is the Freundlich distribution coefficient [M b L3b ], and b is the exponent of theFreundlich isotherm model [] which reads s = k frc b . Four parameters of the two-site sorption
model (k fr, b, x c, f) were estimated using a Levenberg–Marquardt based optimization routine(Streck et al., 1995) to minimize the sum of squared deviations between measured and calculated
hysteretic ad- and desorption isotherms.
2.3. Transport studies
Transport experiments were conducted with natural (undisturbed, structured) soil columns
under natural (transient, variably saturated) water flow conditions. We used the soil column
set-up as described by Rambow and Lennartz (1993). The setup was placed in an air
conditioned lab (ambient temperature of 20 8C F2). At the top column boundary, a rainfall
simulator with adjustable sprinkling rate was installed. At the bottom column boundary,constant suction of 3 kPa was applied using a vacuum pump. Three magnetic valves percolumn outlet controlled the effluent collection while maintaining the constant applied suction
at the lower boundary. Tensiometers (0.5-cm diam.) were inserted into the soil columns at
depths of 2.8 and 12.8 cm for pressure head measurements. Each column was placed on abalance to register changes of column mass, from which changes in the total volumetric soil
water content were calculated based on the known column volume. A data-logging and control
unit automatically controlled the column set-up and collected data.
Initial water contents of each column were determined as the difference between initialand oven-dry soil column mass (after the experiment), divided by the soil core volume, to be
0.320 (0.316) cm3 cm 3 for the Ap (Ah) soil. The transport experiment was started byapplying 0.005 L (0.029 cm) of a solution of Br, IPU, and TER during 20 min to the surface ofboth the Ap and Ah column. For Ap (Ah), the application contained 1010 (1034) mg L 1 Br,and a whitish suspension of 254 (336) mg L 1 IPU and 290 (314) mg L 1 TER.
The solubility of IPU in RArelon (TER in RGardoprim 500) is 170 mg L 1 (8.5 mg L 1) at
20 8C (IVA, 1990). The applied doses of active ingredients (0.8–1 kg AI/ha) andconcentrations above solubility were chosen to represent typical agricultural application forcrop protection, as done previously in column experiments (e.g., Roulier and Jarvis, 2003a).
Starting 25 min after the application, a 3-cm (0.5 L) pulse of rain was applied at an intensityof 24 cm d 1. Based on solubility and applied irrigation rate, IPU and TER were completely
dissolved after 3 min and 1 h of the first irrigation, respectively. In weekly intervals, two (Ah:three) more rainfall events of the same depth and intensity were applied. An average
evaporation of 0.3 cm d 1 over the experiment was estimated to be equal to inflow minus
outflow minus change in column mass. Column effluent was continuously collected in 0.016
to 0.018 L samples by a fraction collector. Samples were stored at 4 8C until Br andherbicide analysis.
2.4. Chemical analysis
Chemical analyses for Br were performed using an ion chromatography system with
conductivity detector and for herbicides were done with on-line solid phase extraction (SPE) andreversed phase HPLC/UV–VIS spectroscopy. The detection limits were 0.2 mg L 1 for Br and
0.005 mg L 1 for herbicides in the extracted samples. Details about the chemical analysis
procedures can be found in Meyer-Windel (1998).

8

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

2.5. Hydraulic nonequilibrium (HNE) water flow models
2.5.1. Mobile–immobile model, MIM
The MIM approach assumes that the soil water content, h, is partitioned into a mobile region
water content, h mo [L3L 3], and an immobile region water content, h im [L3L 3] (van Genuchten
and Wierenga, 1976):h ¼ hmo þ him

ð3Þ

One-dimensional transient water flow in the mobile region is described with the Richards’
equation, while a source/sink term accounts for water exchange with the immobile region(Šimlnek et al., 2001, 2003):



Bhmo
B
Bhmo
K ðhmo Þ¼
 1  CMIM
h
Bz
Bt
Bz
Bhim
¼ CMIMhBt

ð4Þ

where t [T] is time, z [L] is the vertical distance positive downward, K [L T 1] is the hydraulicconductivity as a function of the pressure head in the mobile region, h mo [L], and C hMIM [T 1]
is the water transfer rate between the mobile and immobile regions (subscript h refers to
HNE):
CMIM
¼ xh ½hmo  him hð5Þwhere x h [L 1T 1] is a first-order water transfer rate coefficient.
2.5.2. Dual-permeability model, DPM
In the DPM, water flow in the preferential flow paths or fractures (subscript f) and in the soil
matrix (subscript m) are described with Richards’ equations (Gerke and van Genuchten, 1993):


Bhf
Bhf
B
CDPM
Kf
¼
 Kf  h
ð6aÞ
Bz
Bt
Bz
wf


Bhm
B
Bhm
CDPM
Km
¼
 Km þ hBz
BtBz
1  wfð6bÞ

where K m, h m, and K f, h f are the local hydraulic conductivity and water content of thematrix and fracture pore system (e.g., Köhne et al., 2002b), respectively (as opposed to
global state variables h mo, h im and K in Eq. (4) related to the total soil volume), w f is thedimensionless volume fraction of the fracture pore system for converting local to global
variables (e.g., w fh f), and C hDPM [T 1] is the rate of water exchange between the fracture andmatrix regions:P

CDPM¼ x h K m ð hf  hm Þhð7Þ

where x h [L 1T 1] is a first-order water transfer rate coefficient, and the effective hydraulicP
conductivity of the matrix, K m [L T 1], was averaged in terms of h f and h m as follows
P

K m ¼ ðKm ðhf Þ þ Km ðhm ÞÞ:

ð8ÞJ.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

9

2.6. Physical nonequilibrium (PNE) solute transport models
2.6.1. Mobile–immobile model, MIM
The MIM describes solute transport at variably saturated water flow subject to instantaneous
reversible sorption and first-order degradation in the liquid phase of each region as follows

Bhmo cmo
Bfmo qsmo
B
Bcmo
Bqmo cmohmo Dmoþ¼
 lhmo cmo  CMIM
ð9aÞ

p
Bz
Bt
Bt
Bz
Bz
Bhim cimBð1  fmo Þqsim
þ¼  lhim cim þ CMIM
p
Bt
Bt

ð9bÞ

where for the mobile (subscript mo) region, c mo is the solute concentration [M L 3], D mo is thedispersion coefficient [L2 T 1], q mo is the volumetric flux density [L T 1], f mo is the fraction ofsorption sites in the mobile region [], C pMIM is the solute transfer rate between the two regions(subscript p refers to PNE) [M L 3 T 1], and l [T 1] represents the first-order rate coefficientfor degradation in the liquid phase. Furthermore, c im is the solute concentration in the immobile
region [M L 3], and s mo and s im are the adsorbed solute concentrations [MM 1] in the mobile
and immobile regions, respectively, which were calculated using Freundlich equations:smo ¼ kfr cbmo

ð10aÞ

sim ¼ kfr cbim

ð10bÞ
3b

 b

where k fr is the Freundlich distribution coefficient [L M ], and b is the Freundlich model
exponent []. The dispersion coefficient in the mobile domain, D mo, was described with thefollowing equation:
hmo Dmo ¼ hmo D0

h7=3
moþ kmo jqmo jh2s;moð11Þ

where D 0 is the molecular diffusion [L2 T 1], k mo signifies the dispersivity in the mobile region[L], and h s,mo denotes the saturated water content in the mobile region. The rate of solute masstransfer between the mobile and immobile regions, C pMIM, is given by (Šimlnek et al., 2003):
CMIM¼ xMIMðcmo  cim Þ þ CMIM
c4
p
p
h

ð12Þ

where c* is equal to c mo for C w N 0 and c im for C w b 0, and x pMIM is the constant first-order(PNE) diffusive solute mass transfer coefficient [T 1].
2.6.2. Dual-permeability model, DPM
Solute transport in the fracture (subscript f) and matrix (subscript m) pore regions of the DPM
were described with convection–dispersion equations (Gerke and van Genuchten, 1993):

CDPM
Bhf cf
Bqsf
B
Bcf
Bqf cfp
þ
¼
 lhf cf 
hf Df

Bz
Bt
Bt
Bz
Bz
wf

ð13aÞ



CDPM
Bhm cmBqsm
B
Bcm
Bqm cmp
þ
¼
 lhm cm þ
hm Dm

Bz
Bt
Bt
Bz
Bz
1  wf

ð13bÞ

10

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

where for the fracture and matrix regions, c is the solute concentration [M L 3], D is thedispersion coefficient [L2 T 1] which was described with an equation analogous to (11), q is thevolumetric flux density [L T 1], s is the adsorbed concentration [MMv1], and solute transfer
C pDPM [M L 3T 1] is calculated as follows
CDPM¼ hm xDPM
p
p

h7=3
m
ð1  wf Þðcf  cm Þ þ CDPM
c4
h
h2s;m

ð14Þ

where x pDPM is a first-order diffusive solute mass transfer coefficient [T 1], c* is equal to c f for
C hDPM N 0 and c m for C hDPM b 0, and the adsorbed concentrations in the fracture and matrix
regions, s f and s m [MM 1], were calculated assf ¼ ff kfr cbf

ð15aÞ

sm ¼ kfr cbm

ð15bÞ

where f f is an empirical factor [] used to increase ( f f N 1) or decrease ( f f b 1) sorption in thefracture pore region.
2.6.3. Dual-permeability mobile–immobile model: DP-MIM
The DPM solute transport model described in Section 2.6.2 was, following Pot et al. (2005),
extended into a triple-porosity model by assuming that the matrix water content, h m, can befurther partitioned into mobile and immobile fractions similarly as in Eq. (3):hm ¼ hmo þ him :ð16Þ

In the DP-MIM, solute transport in the fractures was described using Eq. (13a), while transport
in the matrix was described using Eqs. (17a,b) with first-order diffusive solute transfer between
the mobile and immobile matrix regions, and assuming time-variable h mo and constant h im:
CDPM
Bhmo cmo
Bqfmo smo
B
Bcmo
Bqmo cmopþ
¼
 hmo lcmo þ
hmo Dmo

Bz
Bt
Bt
Bz
Bz
1  wf
 xMIMðcmo  cim Þ
p
him

Bcim
Bqð1  fmo Þsimþ
¼ xMIMðcmo  cim Þ  him lcim
p
Bt
Bt

ð17aÞ
ð17bÞ

where s mo and s im were described with Freundlich Eqs. (10a) and (10b). In DP-MIM, f mo is thefraction of sorption sites in the mobile region of the matrix domain [].
2.7. Physical and chemical nonequilibrium solute transport model
2.7.1. Dual-permeability kinetic sorption transport model: DP-KM
In analogy to the two-site kinetic sorption model for bulk soil (Eqs. (1) and (2)), the adsorbed
solute concentrations in the fracture pore system, s f, and the matrix pore system, s m, of the DPKM were partitioned into equilibrium and rate-limited sorption sites
sf ¼ se;f þ sk;f

ð18aÞ

sm ¼ se;m þ sk;mð18bÞ

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

11

where for fracture (subscript f) and matrix (subscript m) domains, s e (MM 1) denotes theadsorbed solute concentration at equilibrium (type-1) sorption sites, and s k (MM 1) symbolizes
the remaining adsorbed concentration of a solute attached to kinetic (type-2) sites of the solid
phase assuming the following first-order rate equations (Pot et al., 2005):
h
i
Bsk;f
¼ xc;f ð1  ff Þkfr;f cbf  sk;f
ð19aÞ
Bt


Bsk;m
¼ xc;m ð1  fm Þkfr;m cbm  sk;m
Bt

ð19bÞ

where for fracture (subscript f) and matrix (subscript m) domains, x c,f and x c,m are first-ordersorption rate constants [T 1] and f f and f m are the fractions of equilibrium sorption sites [].
The DP-KM approach describes solute transport at concurrent PNE and CNE conditions
accounting for two-site kinetic reversible nonlinear sorption as follows:
Bhf cf Bqsf
B
Bcf
hf D f
þ
¼
Bz
Bt
Bt
Bz



DPM



Cp
Bqf cf
 lhf cf 
Bz
wf

h
i
 qxc;f ð1  ff Þkfr;f cbf  sk;f
ð20aÞ



CDPM
Bhm cm Bqsm
B
Bcm
Bqm cmphm D mþ¼
 lhm cm þ

Bz
Bt
Bt
Bz
Bz
ð1  wf Þ qxc;m ð1  fm Þkfr;m cbm  sk;m :ð20bÞ

A model somewhat similar to the DP-KM, although assuming linear adsorption isotherms,
was recently presented by Ray et al. (2004).
2.8. Model parameter estimation
Model approaches in this study were parameter demanding, although to different extents(Table 2). To reduce the number of unknowns, some parameters were always fixed as follows.
The parameter s was set equal to 0.5 (Mualem, 1976). The parameters h mo,r and h im,r (MIM)
and h m,r and h f,r (DPM) were set to zero, which is consistent with zero h r values fitted to water
retention and hydraulic conductivity data for the Ap and Ah soil columns (Meyer-Windel,
1998). The diffusion coefficient, D 0, was calculated for Br (1.797 cm2 d 1) and for herbicides(0.7 cm2 d 1) according to Atkins (1990). Using dye tracer experiments with columns takenfrom the same soil horizons (Köhne, 1999), w f was set equal to average values of dye coverage(Ap: w f = 0.2; Ah: w f = 0.05). Similarly, k m (0.5 cm) was derived from Br transport
experiments performed with the Ap and Ah soil columns at unsaturated ( 3 kPa) steadystate flow conditions (Meyer-Windel, 1998). In forward herbicide transport simulations,additional parameters were fixed. The parameter f mo in the MIM and DP-MIM was assumed tobe equal to h mo,s / h s, whereas for DP-KM, f f and f m were assumed to be equal to f as
determined in the batch test for total soil (Table 3). The sorption coefficients k fr, b, and x c werealso fixed based on the batch experiments results and were as a first approximation assumed tobe identical for different regions. In DPM and DP-MIM, the empirical factor f f increasing or
decreasing sorption in PFPs was set to its default value of one. The degradation coefficient l at

12

Simulated process
Variably saturated transientwater flow andwater transfer
between regions(HNE)+Convective – dispersive
solute transportand solute transfer
between regions(PNE)+Nonlinear
instantaneous sorption
and decay+Kinetic sorption
(CNE)

Data

MIM

N
a

DPM

N
a

DP-MIM

N
a

DP-KM

N
a

q, h, h(s, h mo,r, h im,r)K s, h mo,s, a mo,
n mo, h im,s, a im,
n im, x h8(h m,r, s m, h f,r, s f)
h m,s, a m, n m, K m,s,
a f, n f, h f,s, K f,s, x h

9(h mo,r, s mo, h f,r, s f)
h mo,s, a mo, n mo, K mo,s,
h f,s, a f, n f, K f,s, x h

9(h m,r, s m, h f,r, s f)
h m,s, a m, n m, K m,s,
a f, n f, h f,s, K f,s, x h

9BrBTC

D 0a, k, x p

2(D 0, k m)a, k f, x p

2(D 0, k mo)a, h im, k f,
x DPM
, x MIMpp

4(D 0, k m)a, k f, x p

2Batch tests
IPU BTC(k fr, b)a,
f mo, l2

(k fr, b)a, f f, l

2

(k fr, b)a, f f, f mo, l

3

TER BTC

n.a.

n.a.

(k fr, b, f m, x c,m)a,
f f, x c,f, l

3

N = maximum number of fitted parameters, n.a. = not available.
aFixed at known or estimated value.

n.a.

n.a.

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

Table 2Model approaches with parameter requirements to describe transport of Br, IPU, and TER

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

13

Table 3Equilibrium (Ah-soil) and two-site kinetic (Ap-soil) sorption parameters of TER and IPU as obtained from batch tests
Soil-herbicide

k fr

b

f

xc

R2

Ap-TER
Ap-IPU
Ah-TER
Ah-IPU

11.99
6.20
3.73
1.33

0.92
0.86
0.92
0.89

0.499
0.375
–
–

0.106
0.105
–
–

0.978
0.981
0.999
0.993

20 8C was set to 0.046 d 1 for IPU (half-life of 15 d) and 0.011 d 1 for TER (half-life of 60 d),
within the typical range of measured values (e.g., PerrinGanier et al., 1996; James et al., 1998;
Walker et al., 2002).
The inverse parameter estimation was accomplished as follows. The water flow andconservative solute transport parameters (Table 2) were estimated using the experimental
observations of cumulative water flow, pressure heads, water contents, and Br effluent
concentrations. For herbicide transport simulations, the resulting hydraulic and transportparameter values were fixed, except that x pMIM and x pDPM were both multiplied by the ratio
between D 0 of herbicides and Br (0.7 / 1.797). Most reactive transport parameters were fixed attheir batch test values. Sorption parameters of the PFPs, and the degradation constant l assumed
to be similar in PFPs and matrix, were inversely estimated based on the measured herbicide
effluent concentrations. For the various approaches, f mo and f f were estimated either alone or
together with l (DP-KM: also x c,f and f f) (Table 2). Limiting the optimization of sorption
properties to only the PFP region was done for three reasons: (i) PFP sorption properties
are much more sensitive to herbicide transport during preferential flow than matrix sorption
properties (Pot et al., 2005), (ii) the batch values mainly represented the matrix properties, since
the matrix made up the bulk of soil, and (iii) the number of fitted parameters should be
kept small.
Recently, an increasing number of theoretical and experimental studies reported that thedispersivity may not be a constant for a soil, but increases from conservative to (particularly
nonlinearly) adsorbed solutes (e.g., Pot and Genty, 2005; Vanderborght et al., in press). To
obtain a first impression of the potential significance of this effect for dispersion in preferential
flow paths, k f was additionally fitted in few selected herbicide transport simulations.
The inverse parameter estimation was performed by Levenberg–Marquardt minimization of
the objective function U(Šimlnek et al., 1998):UðbÞ ¼

m
X
j¼1

vj

n
X


2
wi;j Oj ð z; ti Þ  Ej ð z; ti ; bÞ

ð21Þ

i¼1

where m represents the number of different measurement sets, n is the number of observations
in a particular measurement set, O j (z,t i ) are observations at time t i at location z, E j (z,t i ,b) are
the corresponding estimated space–time variables for the vector b of optimized parameters(Table 2), and v j and w i,j are weighting factors associated with a particular measurement set orpoint, respectively. When simulating the Br transport experiment, O j (z,t i ) represented pressure
heads in the mobile region at 2.8 and 12.8 cm depths, average soil column water contents(averaged over the column depth and over regions), cumulative water fluxes across the lower
boundary, and effluent Br concentrations. For simulating the herbicide transport experiment,O j (z,t i ) included only herbicide concentrations in the effluent. The w i,j in Eq. (21) were set

14

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

equal to one assuming similar error variances within a particular measurement set. The v j were
calculated as (Šimlnek et al., 1998):
1
vj ¼
ð22Þ
Nj r2jwhere r2j is the variance within the j-th measurement set, and N j is the number of measurements
within the set. The model performance was statistically evaluated using the dimensionlessmodified coefficient of efficiency, E (Legates and McCabe, 1999):N P
jEj ð z; ti ; bÞ  Oj ð z; ti Þj
ð23Þ
E ¼ 1:0  i¼1 P
NP
jOj ð z; ti Þ  O j
i¼1

where E varies between  l to 1.0, with higher values indicating better agreement between themodel and observations. A zero (negative, positive) value of E indicates that the model isstatistically as good as (worse, better than) the mean (Ō) to predict the data, with a unit value
indicating perfect agreement between model and data. Compared to the coefficient ofdetermination, R 2, E is a more conservative and reliable statistical measure and is less sensitive
to extreme values (Legates and McCabe, 1999). Throughout this study we will assume thatE N 0 represents an bacceptable Q simulation, while E N 0.5 represents a bgood Q simulation.
Additionally, the mean absolute error (MAE) was calculated asN PjEj ð z; ti ; bÞ  Oj ð z; ti Þj
MAE ¼ i¼1
ð24ÞN
where MAE carries the same units as the observations. The observed and simulated herbicide
and Br BTCs were compared using the relative and absolute error measures (23) and (24).

Fig. 1. Left: Measured time dependent sorption of IPU and TER with a) Ah soil and b) Ap soil, Right: data of a batch
adsorption-(S0) experiment with 4 desorption (S1 to S4) steps, and corresponding fitted log-linear (Freundlich)isotherms, c) TER with the Ah soil material, and d) IPU with the Ap soil material.

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

15

For the water flow simulations, additionally an average soil hydraulic index of agreement,
E hyd, was calculated as the average over individual E values for outflow, water contents, andpressure heads (E outflow, E water contents, E pressure heads):
Ehyd ¼

Eoutflow þ Ewater content þ Epressure head3

ð25ÞThe initial condition for water flow was represented by a hydrostatic pressure head profile
linearly interpolated between measurements of tensiometers installed at 2.8 and 12.8 cm depth.

Fig. 2. Transient water flow and bromide (Br) transport data from the Ap-soil column, and corresponding model
simulation results, a) bulk soil water contents, b) pressure heads at depths of 3 and 13 cm, c) cumulative outflow, andd) Br effluent concentrations.

16

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

The upper flow boundary condition was a variable flux-type condition equal to time-variable
rates of irrigation (24 cm d 1) and evaporation (0.3 cm d 1) set according to the experimental
conditions. For the lower boundary, an unsaturated seepage face at user-specified pressure headsof  6 kPa (Ah1) or  3 kPa (Ap1) was implemented in HYDRUS-1D. The initial Br andherbicide concentrations were set to zero. The upper and lower boundaries for solute transport
were represented by specified flux and (zero) gradient conditions, respectively. The Br andherbicide concentrations (above solubility) were directly used in the small injection pulse(0.029 cm). In the following 25 min between injection and irrigation, the models internallysimulated herbicide partitioning between soil solid and liquid phases according to the herbicide
sorption and diffusion characteristics, after which the first water irrigation was applied.
Neglecting herbicide dissolution kinetics in the boundary conditions should be negligible for
IPU, which was applied in concentrations just slightly above solubility, and may have some(unknown) effect on TER transport simulations.
3. Results
3.1. Batch experiments
Fig. 1 shows selected results for the relative concentrations for the 4 herbicide-soil
combinations during 72 h of adsorption (Fig. 1a,b) and ad- and desorption data (Fig. 1c,d). For
the Ah soil material and both IPU and TER, the early drop from unity at time zero to more or lessconstant concentrations suggests instantaneous equilibrium adsorption (Fig. 1a). By contrast, for
the Ap soil the almost log-linear decrease of dissolved concentrations of both herbicides over theentire 72 h (Fig. 1b) revealed rate-limited sorption. Accordingly, for the Ah material and themeasured ad- and desorption data together with the isotherms obtained by linear regression ofcorresponding 5 log c eq–log s-pairs were similar for ad- and desorption (Fig. 1c exemplarily
shows Ah-TER), whereas the shift of sequential isotherms revealed hysteretic ad- and desorption
relation for the Ap material and IPU or TER (Fig. 1d shows Ap–IPU as an example). Hence,
instantaneous equilibrium sorption was assumed for both IPU and TER in the Ah soil, andvalues for k fr and b were obtained by linear regression of the Freundlich isotherm to thecorresponding double-logarithmic data of adsorbed versus dissolved concentrations (Table 3).
For the Ap soil, the hysteretic sorption data were analyzed with the model of Streck et al. (1995)
to obtain parameters of the two-site sorption model (Table 3) as described in Section 2.2.
3.2. Bromide transport
3.2.1. Ap soil column
The MIM, DPM, and DP-MIM inverse simulation results for water flow and Br transport in
the Ap soil column are displayed in Fig. 2. We note that simulations were also performed with asingle-domain model based on Richards’ and convection-dispersion equations in HYDRUS-1D(Šimlnek et al., 1998) which failed to describe the observed Br concentrations (not shown).
Observed bulk soil water contents varying between 0.3 and 0.36 (Fig. 2a), pressure heads
between  60 and 0 cm (Fig. 2b), and cumulative outflow (Fig. 2c) for three intermittent
irrigation events were more or less matched by all three models. The good agreement between
model simulations and data was confirmed by the E hyd values around 0.8 (Table 4).
Simulated Br BTCs were within fluctuations of measured concentrations up to 7.5 days
(Fig. 2d). The early Br arrival with the first column effluent, the low Br peak of 5.5 mg L 1,

Table 4
Results for hydraulic and conservative solute transport parameters and corresponding (F) 95% confidence limits for the MIM, DPM, and DP-MIM, obtained by inverse estimation
using water flow and Br concentration data for the Ah and Ap soil columns
Soil

Fracture and mobile region

Matrix and immobile region

Inter-region transfer

Ap
MIM
DPM
DP-MIM

Ah
MIM
DPM
DP-MIM

h§s

a

K§s

k

(cm3 cm 3)(cm 1)(cm d 1)

(cm)

0.148
F 0.0140.159
F 0.006
0.162

0.0341
F0.00200.0312
F0.0006
0.03461.250
F 0.0421.301
F 0.028
1.250

22.55
F0.90
61.62
F2.98
54.88

2.30
F 0.892.75
F 0.402.82

F 0.007

F0.0025

F 0.039

F1.43

F 0.460.055
F 0.006
0.050
F 0.003
0.039

0.0461
F0.0066
0.0136
F0.00100.01401.348
F 0.0711.827
F 0.084
2.188

34.57
F4.60
6.71
F0.40
6.88

2.34
F 2.61
10#
F 4.810.42

F 0.001

F0.001

F 0.111

F0.52

F 0.08nBromide

Kbs

xh

x MIM
p

x DPM
p

(cm d 1)(d 1)(d 1)(d 1)

n.a.

1.336
F 0.274n.a.

n.a.

hbs (DPM–
MIM:h mo,h im)(cm3 cm 3)

a

0.214
F0.0140.212
F0.0080.201;
0.0107
F0.008;0.136

0.0162
F0.00200.0158F0.0004
0.01332.000
F0.3101.937
F0.0971.951

1.645
F 0.008
1.475

0.192
F 0.0360.092
F 0.007
0.156

F0.0005

F0.156

F 0.1380.0045
F0.00280.0011
F0.0000
0.0016

1.441
F0.323
1.362
F0.0391.250

n.a.

F0.0001

F0.0390.313
F0.1800.323
F0.0270.110;
0.222
F0.010;0.031

n(cm 1)

Goodness-of-fit

Water

E (MAE)
E hyd(mg L 1)

0.021

5.35
F 2.412.55

0.43 (0.75)
0.82
0.48 (0.72)
0.773
0.55 (0.62)

F 0.021

F 0.020

F 4.520.82

0.774
F 0.348
n.a.

n.a.

2.741
F 0.218
2.154

0.010
F 0.0080.009
F 0.001
0.010

0.074

0.606
F 0.252
3.62

0.45 (0.92)
0.82
0.64 (0.72)
0.81
0.82 (0.36)

F 0.090

F 0.001

F 0.030

F 0.602

0.81

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

Model

Fixed: h r = 0, k m = 0.5 cm, w f = 0.05 (0.2) for Ah (Ap); n.a. — not available. §For DPM and DPM–MIM, divide value by w f = 0.2 (Ap) and w f = 0.05 (Ah) to obtain thecorresponding local value for the fracture pore region, bfor DPM and DPM–MIM, divide by (1  w f) = 0.8 (Ap) and (1  w f) = 0.95 (Ah) to obtain the corresponding local value forthe matrix, #at upper constraint.

17

18

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

corresponding to only 0.5% of the applied concentration of 1010 mg L 1, and the shape of theBr BTC were delineated by all model approaches. The late-time low Br concentrations
observed after two weeks could be only accurately captured using the DP-MIM (Fig. 2d). For
the entire Br BTC with three irrigation events, the model performance was acceptable for MIM(E = 0.427) and DPM (0.476), and was good for DPM–MIM (0.549). Mean absolute deviations,
MAE, were below 0.75 mg L 1 and decreased in the order of MIM, DPM, and DPM–MIM(Table 4). For the Br BTC only up to 7.5 days, i.e., until outflow after the second irrigation
event ceased, there was a good agreement between all models and data as reflected by E N 0.5.
We conclude that for the Ap soil, the assumption of two pore domains was a suitable concept todescribe Br leaching at HNE and PNE during the first two irrigation events. However, the latetime Br transport appeared to be dominated by diffusive transfer into the microporosity
subdomain of the matrix acting as a sink for Br.
Table 4 furthermore shows estimated model parameters. First-order water and solute transfer
parameters can only approximately be compared between two-flow region and MIM approaches
as used in this study. For example, the fitted x h values controlling water and advective solutetransfer are smaller in DPM and DP-MIM than in MIM, but estimated x h values in DPM and
DP-MIM compensate for the effective matrix hydraulic conductivity, K̄m, in Eqs. (7) and (8).
Similarly, in the solute transfer term (Eq. (14)) of DPM and DP-MIM, the fitted x pDPM values are
multiplied by the variable term (h m10 / 3 / h 2s,m) (1  w f) b 1. The average reduction is compensatedby the larger optimized values for x pDPM as compared to x pDPM in Eq. (12). In DP-MIM,
where x pDPM characterizes transfer between the fracture and matrix domains and x pMIM between
the mobile and immobile regions within the matrix, estimated values of x pMIM and h im were sosmall (Table 4) that diffusive transfer between the mobile and immobile matrix regions only had
a notable effect on simulated Br transport after prolonged time.
To facilitate inter-comparison between MIM and DPM approaches, the h s and K s parameters in
Table 4 are presented as global variables related to the total soil volume. In other words, h s and K svalues listed in Table 4 for DPM and DP-MIM were multiplied by w f and (1  w f) for the fracture
and matrix regions, respectively. Region-specific parameters did not differ much among different
models, except that w f K s,f for DPM and DP-MIM were almost 3 times as large as correspondingK s values in the mobile region of MIM (Table 4). However, the fracture pore region in DPM and
DP-MIM simulations did not wet up entirely (h f =  2 cm) during rainfall applications, while themobile region in the MIM simulation fully saturated (h mo = 0 cm) as was experimentallyobserved. For DPM and DP-MIM, the inverse estimate of soil hydraulic conductivity(K = w fK f + (1  w f)K m) at h =  2 cm came to within 10% of the K mo,s value (22.55 cm d 1)
estimated with MIM. Similarly, the w fh f,s values for DPM and DP-MIM were slightly larger thanh mo,s for MIM (Table 4), whereas similar values were found for h mo,s and w fh f at h f =  2 cm.
Estimated K f,s values for DPM and DP-MIM were thus values extrapolated beyond the simulatedpressure head range and are expected to be associated with large estimation errors.
These findings can be partly explained with the over-parameterization of DPM and DPM–
MIM and may be used as an argument for always using a simpler approach (MIM)
simultaneously with more complex approaches to cross-check the physical meaning ofparameters obtained by inverse estimation. However, recognizing the overestimated K s,f values,
DPM and DP-MIM can still be used in subsequent herbicide transport simulations.
3.2.2. Ah soil column
Experimental results for water flow and Br transport in the Ah column are shown in Fig. 3
together with corresponding model results, while estimated model parameters and statistical

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

19

measures are presented in Table 4. Pressure heads between  200 and 0 cm (Fig. 3a) andcumulative outflow (Fig. 3b) show the hydraulic response of the soil system to four irrigations
and subsequent redistribution events. No water contents were measured for the Ah soil column.
Simulations were in good agreement with hydraulic data (Fig. 3a,b) with E hyd N 0.8 for all threemodels (Table 4).
The Br BTC had a highly left-skewed shape with the peak already encountered in the first
percolate sample, and a higher peak concentration of 18 mg L 1 (2% of the applied
concentration) than for the Ap soil column (Fig. 3c). This is again indicative of preferential flow
with pronounced HNE and PNE. Similarly as for the Ap soil column, both MIM and
DPM approximated the Br breakthrough well during the first week, while only the DP-MIM

Fig. 3. Transient water flow and bromide (Br) transport data from the Ah-soil column, and corresponding model
simulation results, a) pressure heads at depths of 3 and 13 cm, b) cumulative outflow, and c) Br effluent concentrations.

20

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

approach simulated the low measured Br concentrations during the second week (Fig. 3c). Thisvisual evaluation was again supported by the E-statistics that increased from 0.449 for MIM to
0.644 for DPM and 0.818 for DPM–MIM, with MAE decreasing in the same order (Table 4).
The narrow Br BTC in the Ah soil was reflected by a x h value 10 times smaller than for the
Ap soil (Table 4). In the MIM, the K s,mo was larger for Ah than for the Ap soil, as expected for the
more pronounced preferential flow. However, for the DPM and DP-MIM, values for w f K s,f were
smaller for Ah than for Ap, and this discrepancy remained even for local K s,f values obtained by
dividing w f K s,f (Table 4) by w f . However, the K s,f values for the Ap soil were overestimated as
explained above. Moreover, smaller a f and larger n f values for Ah (Table 4) mean relatively larger
unsaturated K f values than for Ap. Hence, K s,f values fitted under variably saturated conditions

Fig. 4. Transient water flow and isoproturon (IPU) transport data from the Ap-soil column, and corresponding model
simulation results, a) forward simulations, b) inverse simulation 1, c) inverse simulation 2, and d) inverse simulation
3 (DPM fl indicates the dual permeability model with fitted f and l parameters).

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

21

are less suitable indicators of the extent of preferential flow than, e.g., x h and x p. As expected,
the x p values in MIM and DPM were also smaller for Ah than for Ap. Only for the DP-MIM
approach, x pDPM (3.62 d 1), x pMIM (0.07 d 1) and h im (0.222) were larger, and k f (0.42 cm)
smaller than for Ap (Table 4). We note that for all approaches, x p and k f were to some extent
negatively correlated ( 0.7 b r 2 b  0.3, not further shown). The best fit of the Br BTC suggeststhat the DP-MIM triple-porosity approach best represented the dispersivity and solute transfer
between macropores, inter-aggregate regions, and aggregates of the Ah soil. The small k f(0.42 cm) value estimated using DP-MIM may consequently best represent the actual dispersivityof the macropore domain of the Ah soil column. The striking differences between x f and x p
values for the relatively similar model approaches (Table 4) point to the need of further improving
experimental designs for parameter identification of two- and multi-region models.
Overall, the two-region concept appeared to be a reasonable approximation of the actual HNE
and PNE solute transport in the Ah soil, with the two mobile regions of the DPM approachproviding a somewhat better description of the data. The decrease of Br concentrations at latertimes, similarly as for the Ap soil, could be reproduced only with the DP-MIM approach thatconsidered diffusive solute transfer into an immobile sub-region of the matrix.
3.3. Herbicide transport
3.3.1. Ap soil column: isoproturon transport
For the weakly adsorbed herbicide IPU, Fig. 4 shows the measured BTC together withsimulation results using various forward (Fig. 4a) and inverse (Fig. 4b,c,d) modeling approaches.
Isoproturon appeared simultaneously with Br in the effluent, but at much lower peak
Table 5Isoproturon transport parameters and corresponding (F) 95% confidence limits for the various model approaches used
for simulating the observed IPU breakthrough curve in the Ap soil column
Approachf mo, f f§

x c,f(d 1)

ForwardMIM
DPM
DP-MIM
DP-KM

0.408
1
1
0.375

–
–
–
0.105

InverseDPM f
DPM flDP-MIM f
DP-MIM fl
DP-KM f
DP-KM fx
DP-KM fl0.240 F 10.50.240 F 0.001
0.236 F 1.5
0.229 F 0.005
0.258 F 2.8
0.242 F 0.003
0.254 F 0.001

–
–
–
–
0.105
0.045 F 0.024
0.105yf m, f mo

k mo, kbf(cm)

l(d 1)E (MAE)(mg L 1)

–
0.043
0.375

2.30
2.75
2.82
2.75

0.046
0.046
0.046
0.046

0.54 (0.088)
0.54 (0.088)
0.25 (0.068)
0.03 (0.053)

–
–
0.043
0.043
0.375
0.375
0.375

2.75
2.75
2.82
2.82
2.75
2.75
2.75

0.046
0.410 F 0.126
0.046
0.761 F 0.146
0.046
0.046
0.421 F 0.221

0.30 (0.038)
0.51 (0.026)
0.01 (0.055)
0.49 (0.028)
0.40 (0.033)
0.48 (0.028)
0.45 (0.030)

Fitted parameters are in bold letters. Default parameters see in Tables 3 and 4; x MIM
and x DPM
(Table 4) were multipliedppby diffusion correction factor 0.39 (=0.7 / 1.797). §Fraction of sorption sites in the mobile region ( f mo; MIM), or factor toincrease or decrease sorption in the fracture pore system ( f f; DPM, DP-MIM), or fraction of equilibrium sorption sites in
the fracture pore system ( f f; DP-KM). Multiply f f with w f (0.2) to obtain global value ( f f w f) that can be compared withf mo. yFraction of equilibrium sorption sites in the matrix ( f m; DP-KM) or fraction of sorption sites in the mobile region of
the matrix ( f mo; DP-MIM). b Dispersivity in the mobile region (k mo; MIM) or in the fracture pore system (k f; DPM, DPMIM, DP-KM).

22

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

concentrations below 0.2 mg L 1. The IPU breakthrough was characterized by concentration
drops during redistribution times between irrigations, which were not observed in Br
breakthrough.
The MIM forward approach predicted zero concentrations for all times (not shown), and the
MIM was not further evaluated for inverse parameter estimation. While none of the forward
modeling approaches approximated the observed IPU breakthrough, only DP-KM predicted anearly concentration increase (Fig. 4). Accordingly, among the forward approaches only DP-KM
had a non-negative E statistics (Table 5). Hence, DP-KM forward simulation results confirm therelevance of kinetic sorption (CNE) for IPU transport.
However, the DPM and DP-MIM equilibrium sorption approaches apparently much
improved when fitting only f f (Fig. 4b,c), even if huge confidence intervals revealed large
uncertainty for fitting f f alone (Table 5). Note that we refer to DPM with optimized parameter f f
as DPM f in Fig. 4b and that we use similar notation for other inverse optimizations. Inverse
estimates were similar for the DPM f ( f f = 0.240) and DP-MIM f ( f f = 0.236), and were thus fourtimes smaller than assumed in the forward simulations (Table 5). The IPU concentrations afterFig. 5. Transient water flow and terbuthylazine (TER) transport data from the Ap-soil column, and corresponding model
simulation results, a) forward simulations, b) inverse simulation 1, c) inverse simulation 2, and d) inverse simulation 3.

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

23

one week were overestimated with DPM (Fig. 4b) and even more so with DP-MIM (Fig. 4c),which resulted in acceptable (E = 0.301) and poor (E =  0.011) model performance for DPM and
DP-MIM, respectively. When f f and l were both fitted in DPM and DP-MIM, the IPU
concentrations after one week were matched well (Fig. 4b,c), which increased E to values
around 0.5 (Table 5). Estimates for f f hardly changed (with associated confidence intervals
drastically decreased), while l increased from 0.046 to 0.41 d 1 (DPM fl) or 0.76 d 1
(DP-MIM fl) corresponding to very short half-lives of only 1.7 or 0.9 days for DPM fl and
DP-MIM fl, respectively.
The best BTC match (E = 0.40; Table 5) when fitting only f f (= 0.258) was achieved using theinverse DP-KM f (Fig. 4d). Gradual further improvement for late times was obtained by fitting f f
and x c,f (E = 0.48), or f f and l (E = 0.45), giving very similar curves (Fig. 4d, Table 5). The
observed IPU loss of 0.090 mg or 7.1% of 1.27 mg IPU applied was approximated by the model
approaches, e.g., the DP-KM fx approach simulated a mass loss of 0.074 mg or 5.8%.
Additionally, the hypothesis was tested that k f may be different for conservative Br andnonlinearly adsorbed IPU. Jointly optimizing k f with f f resulted in the best approximation of
the first IPU peak (Fig. 4d). However, the large positive correlation (r N 0.9) between k f (15 cm;
at upper constraint) and f f showed that these two parameters could not be identified
simultaneously.
3.3.2. Ap soil column: terbuthylazine transport
The measured and simulated TER BTCs are shown in Fig. 5 for forward (Fig. 5a) and inverse(Fig. 5b,c,d) approaches. The breakthrough of the more strongly adsorbed TER started
simultaneously with Br and IPU, but at still lower concentrations below 0.06 mg L 1. Forward

Table 6
Terbuthylazine (TER) transport parameters and corresponding (F) 95% confidence limits for various model approaches
used for simulating the observed TER breakthrough curve in the Ap soil column
Approachf mo, f f§

x c,f(d 1)

yf m, f mo

k mo, kbf(cm)

l(d 1)E (MAE)(mg L 1)

ForwardMIM
DPM
DP-MIM
DP-KM

0.408
1
1
0.499

–
–
0.106

–
–
0.043
0.499

2.30
2.75
2.82
2.75

0.011
0.011
0.011
0.011

1.36 (0.026)
1.36 (0.026)
1.36 (0.026)
1.18 (0.024)

InverseDPM f
DPM flDP-MIM f
DP-MIM fl
DP-KM f
DP-KM fx
DP-KM fl0.401 F1.6
0.191 F 0.002
0.500 F 1.6
0.185 F 0.003
0.188 F 3.4
0.189 F 0.003
0.186 F 0.002

–
–
–
–
0.106
0.095 F 0.037
0.106

–
–
0.043
0.043
0.499
0.499
0.499

2.75
2.75
2.82
2.82
2.75
2.75
2.75

0.011
1.0b F 0.070
0.011
1.0b F 0.052
0.011
0.011
0.340 F 0.113

0.93 (0.021)
0.01 (0.012)
1.04 (0.022)
0.19 (0.012)
0.21 (0.008)
0.20 (0.008)
0.37 (0.006)

Fitted parameters are in bold letters. See Tables 3 and 4 for default parameters. §Fraction of sorption sites in the mobile
region ( f mo; MIM), or factor to increase or decrease sorption in the fracture pore system ( f f; DPM, DP-MIM), or fractionof equilibrium sorption sites in the fracture pore system ( f f; DP-KM). Multiply f f with w f (0.2) to obtain global value ( f f
w f) that can be compared with f mo. yFraction of equilibrium sorption sites in the matrix ( f m; DP-KM) or fraction ofsorption sites in the mobile region of the matrix ( f mo; DP-MIM). bDispersivity in the mobile region (k mo; MIM) or in thefracture pore system (k f; DPM, DP-MIM, DP-KM).

24

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

simulations did not match observed concentrations, with only DP-KM predicting concentrations
gradually increasing after one week (Fig. 5a).
The inverse DPM (DP-MIM) approaches failed to simulate TER breakthrough when fitting
f mo ( f f) alone, but simulations improved upon concurrent inverse estimation of f mo and l( f f and l) (Fig. 5b,c). However, based on zero or negative E, neither approach reproduced TER
breakthrough any better than the mean TER concentration (Table 6).
Significant further improvement could only be attained using the DP-KM approach (Fig. 5d).
According to the E statistics, acceptable simulation results were achieved by inverse estimation off f (0.188) yielding E = 0.21, while simultaneously fitting f f and l further increased E to 0.37
(Table 6). Simultaneous estimation with f f made l reach its upper constraint of 1.0 d 1 in DPM and
DP-MIM, while l reached 0.34 d 1 in the DP-KM fl approach (Table 6). Simultaneously fitting k f
and f f in DP-KM did not significantly improve the simulation results (Fig. 5d).
Mean absolute deviations, MAE, were smaller for TER than for IPU, due to smaller TER
concentrations. The estimation of the total TER mass loss after 14.2 days using the DP-KM flapproach was quite accurate. While the observed mass loss was 0.039 mg or 2.3% of the applied
1.45 mg TER, the simulated loss was 0.031 mg or 2.1%.
3.3.3. Ah soil column: isoproturon transport
Results for the IPU breakthrough are illustrated in Fig. 6 along with corresponding forward(Fig. 6a) and inverse (Fig. 6b) model simulations. The IPU BTC displayed the same extremely
left-skewed shape as Br, peaking already in the first column effluent, which is characteristic ofpreferential flow. The peak concentration of 3 mg L 1 was over 10 times higher than the IPU
peak in the Ap soil column. Furthermore, the IPU BTC did not show concentration drops during
redistribution phases (Fig. 6), in contrast to herbicide BTCs in the Ap soil column (Figs. 4, 5).

Fig. 6. Isoproturon (IPU) transport data from the Ah-soil column, and corresponding model simulation results, a) forward
simulations, b) inverse simulation.

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

25

Table 7Isoproturon (IPU) and Terbuthylazine (TER) transport parameters and corresponding (F) confidence limits for various
model approaches used for simulating the observed IPU breakthrough curve in the Ah soil column
Approach

Isoproturon
f f, f §moTerbuthylazineyf m, f mo

E (MAE)(mg L 1)

f f, f §moy
f im, f mo

(cm)

k mo, kbf(cm)

E (MAE)(mg L 1)

k mo, kbf

ForwardMIM
DPM
DP-MIM

0.150
1
1

0.850
–
0.363

2.34
10.0
0.42

0.12 (0.41)
0.15 (0.39)
0.34 (0.61)

0.150
1
1

0.850
–
0.363

2.34
10.0
0.42

0.39 (0.19)
0.09 (0.15)
0.12 (0.16)

InverseMIM f
DPM f
DP-MIM f
DP-MIM fk

0.174 F 1.2
0.728 F 17.2
0.705 F 13.1
1.548 F 0.04

0.826
–
0.363
0.363

2.34
10.0
0.42
15#

0.17 (0.38)
0.36 (0.29)
0.21 (0.55)
0.54 (0.21)

0.084 F 12.6
1.260 F 3.3
0.855 F 100
1.825 F 0.09

0.916
–
0.363
0.363

2.34
10.0
0.42
15#

0.01 (0.14)
0.34 (0.09)
0.30 (0.18)
0.44 (0.08)

Fitted parameters are in bold letters. Fixed: b (IPU) = 0.89, l (IPU) = 0.046 d 1, b (TER) = 0.92, l (TER) = 0.011 d 1,
and parameters in Tables 3 and 4. §Fraction of sorption sites in contact with the mobile region ( f mo  MIM) or empiricalfactor to increase or decrease sorption in the fracture pore system ( f f; DPM, DP-MIM). Multiply f f with w f (0.05) toobtain global value ( f f w f) that can be compared with f mo. yFraction of equilibrium sorption sites in the matrix ( f m; DPKM) or fraction of sorption sites in contact with mobile region within the matrix ( f mo; DP-MIM). bDispersivity themobile region (k mo; MIM) or in the fracture pore system (k f; DPM, DP-MIM). #At upper constraint.

The MIM, DPM and DP-MIM forward approaches assuming instantaneous sorption provided
much better predictions of the observed IPU BTC for the Ah (Fig. 6a) than corresponding
forward models did for the Ap soil column. Furthermore, the DP-KM approach did not yieldfurther improvement, even when optimizing f f and x c,f (not shown). These results suggest that
there was sorption equilibrium within the PFP of the Ah soil column even during fast
preferential transport.
Both DPM and DP-MIM predicted secondary peaks around 0.25 days. However, around0.25 days the column outflow almost ceased (Fig. 3b), such that only one sample could be
retrieved between 0.15 and 0.3 days (Fig. 6). Due to minimal outflow, the relevance of thesecondary peaks for simulated IPU mass loss was negligible.
The best (acceptable) goodness-of-fit statistics among the forward approaches (Table 7) was
obtained for the DPM (E = 0.15). Optimizing f f improved the DPM description of the IPU BTC(Fig. 6b) and provided the best goodness-of-fit (E = 0.36). The total IPU mass loss after
20.2 days (0.54 mg or 32% of applied) was somewhat overestimated using the DPM f approach(0.81 mg or 48%). The inverse DP-MIM f and MIM f simulations could not approximate the first
IPU peak (Fig. 6b, Table 7). For DPM f and DP-MIM f approaches, inverse estimates for f f wereslightly above 0.7 (Table 7). Estimating both f f and l gave no overall improvement for any
approach (not shown), while fitted large values of f f and k (Table 7) at moderate correlation(r = 0.4) enabled the DP-MIM to simulate the early IPU BTC rise (Fig. 6b).
3.3.4. Ah soil column: terbuthylazine transport
Finally, breakthrough data and simulation results for the more strongly adsorbed herbicide
TER are depicted in Fig. 7 for forward (Fig. 7a) and inverse (Fig. 7b) model approaches. The
TER BTC appeared concurrently as IPU and Br, but with a less sharp and lower peak of
0.5 mg L 1. The agreement between models and observed TER concentrations was less goodthan for IPU (Fig. 7). The forward MIM approach severely underestimated TER concentrations.

26

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

Fig. 7. Terbuthylazine (TER) transport data from the Ah-soil column, and corresponding model simulation results, a)forward simulations, b) inverse simulation.

The forward DPM approach predicted the early TER appearance in the effluent, but
overestimated the peak concentration. The DP-MIM approach predicted a later TER
concentration increase than observed (Fig. 7a).
Inverse estimation of f mo improved the agreement between MIM results and data, but without
a prediction of the early increase of TER concentrations. Among all approaches, fitting f f in the
DPM approach resulted in the best (acceptable) agreement with the data (E = 0.338,
MAE = 0.09 mg L 1, Table 7). The total TER mass loss after 20.2 days was reproduced using
the DPM f approach. While the observed mass loss was 0.31 mg or 20% of the applied 1.57 mg
IPU, the simulated mass loss was 0.33 mg or 21%. In DP-MIM, the early peak could not be
reproduced by fitting f f, even when assigning higher weights, w i,j , to initial TER observations(Table 7, approach DPM–MIM f).
For DPM f and DP-MIM f approaches, inverse estimates for f f were above and below one,respectively (Table 7). Again, joint optimization of f f and l did not improve results (not shown),
while fitted large values of f f and k (Table 7) in the DP-MIM improved simulation of the early
TER BTC rise (Fig. 7b).
4. Discussion
4.1. Sorption
For the Ap soil column, the inverse DPM, DP-MIM, and DP-KM simulations of IPU
transport produced comparable estimates for f f (0.19 to 0.27). Slightly smaller f f values (0.16 and0.19) were found for the more strongly adsorbed TER. Overall, DPM and DP-MIM inverse
simulation results suggest small equilibrium sorption capacity in the PFPs, while the DP-KM
approach suggests rate-limited sorption in the PFPs. There is however much evidence for theJ.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

27

latter explanation: (i) IPU and TER batch tests revealed kinetic sorption in the bulk soil,(ii) among forward simulations only the DP-KM predicted fast herbicide leaching, (iii) amongmodel approaches with only f f fitted, the DP-KM best described the herbicide concentration
drops during the one-week redistribution period.
For the Ah soil column, where herbicide sorption was instantaneous according to batch test
results, the inverse herbicide transport simulations yielded larger f f values of around one (orabove one, if k f was additionally fitted). The f f values above one were in accord with reported
findings of batch tests that organic linings of earthworm burrows may increase the sorption ascompared to the soil matrix (Stehouwer et al., 1993, 1994; Mallawatantri et al., 1996). On theother hand, the predictive capacity of batch sorption tests performed with material scratchedfrom earthworm burrows may be limited for transport conditions with variable saturation, where
effects such as partial wetting of burrow walls might counterbalance the increased sorption.
Sorption often appeared to be reduced in mobile soil water regions of aggregated or structured
soils, as we found for our Ap soil column. For example, Vanderborght et al. (2002) utilized MIM
simulations of Cl and fluorescent dye tracer transport in undisturbed soil columns at steadystate flow and found a smaller sorption capacity of the mobile region, f mo, for the more strongly
adsorbed dye sulforhodamine B ( f mo b h mo / h) than for brilliant sulfaflavine ( f mo = h mo / h). Thisfinding was experimentally confirmed by measurements of lower sulforhodamine B concentrations in the mobile region (Vanderborght et al., 2002). Applications of the dual-permeability
model MACRO for simulating preferential pesticide leaching yielded values for the fraction ofsorption sites in the macropore, f, which appear rather small (although direct comparison withf mo or w f f f values is limited by the conceptual model differences). For example, f in MACRO
was assumed to be 0.02 to describe bentazone leaching (Larsson and Jarvis, 1999) and f was
fitted to values ranging from 0.008 to 0.048 to describe transport of the herbicide MCPA in
different soils (Roulier and Jarvis, 2003a). The latter authors found f to be associated with large
estimation uncertainty (Roulier and Jarvis, 2003a,b), which is in agreement with our results
showing large confidence intervals for f f when fitted alone.
In our inverse DPM and DP-MIM simulations of herbicide transport in the Ap soil column,
the strong reduction of f f (from unity down to 0.2) as compared to less reduction of f f (from batch
test values) in the DP-KM approach suggest that some previously reported small f estimates may
simply result from ignoring chemical nonequilibrium. The acceleration of pesticide leaching byslow sorption reaction rates within the preferential flow domain was also shown by comparing
dual-permeability model scenarios with kinetic and equilibrium sorption (Ray et al., 2004).
In a recent study by Pot et al. (2005), a set of experiments to study the transport of Br, IPU,
and metribuzin through structured undisturbed soil columns was performed at five different
steady-state flow rates. The effect of flow velocity on herbicide transport was evaluated using
comparable models as in our study. Herbicide sorption in the preferential flow region was foundto be kinetic, and only combined PNE–CNE models were successful in simulating the BTCs.
Similar as for our Ap soil, no model could capture the early rise of the herbicide concentrations,while the DP-KM approach yielded the best agreement with the observed BTCs. The DP-KM
approach for low and intermediate flow velocities gave similar values ( f f = 0.22 to 0.29) as for
our Ap column, while f f was lower for the highest flow velocity ( f f = 0.12).
4.2. Degradation
For IPU and TER leaching in the Ap soil column, inverse estimation gave values for the firstorder degradation rate in the liquid phase, l, in the range of 0.3–1.0 d 1, while typical values

28

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

are at least an order of magnitude smaller (e.g., James et al., 1998; Walker et al., 2002). Muchlarger magnitudes of degradation rates in laboratory column studies (particularly underunsaturated flow conditions) as compared to batch sorption experiments were found before,
and the results of several such studies were reviewed in Pot et al. (2005). For example for inverse
simulation of MCPA transport using the MACRO model, l values of 0.3 – 0.4 d 1 wereapparently well identified with very low uncertainty (Roulier and Jarvis, 2003a). In transport
experiments under unsaturated flow conditions, Pot et al. (2005) found l values up to 31 (IPU)
and 18 times (metribuzin) larger than the literature values (Agritox, 2001). However, thedegradation constant increased with the extent of nonequilibrium sorption, such that theestimated high degradation rates may have only reflected chemical irreversible or hysteretic
sorption, defined as fast attachment and much slower detachment. Irreversible or hysteretic
sorption may thus have been lumped together with degradation and mineralization into l(Pot et al., 2005). This argument can be corroborated to some extent by our own results. Whilefitting of l did not improve simulation results for the Ah soil where kinetic (irreversible)
sorption was not involved, it did apparently improve simulations of herbicide transport in the Ap
soil where kinetic sorption was important. Moreover, jointly optimizing f f and x c,f or f f and l in
DP-KM (and even f f and l in DPM and DP-MIM) produced comparable simulated BTCs.
Hence, effects of sorption rate, x c,f, and degradation rate, l, on herbicide transport cannot be
jointly identified based on BTC fitting.
4.3. Dispersion
For the Ah soil, DP-MIM was the best approach for simulating Br transport but not for
simulating herbicide transport. This discrepancy could be explained when assuming that thedispersivity was not the same for Br and herbicides. For example, Vanderborght et al. (2002)
found that MIM simulation results of dye tracer BTCs could be significantly improved whenassuming different dispersivity values for two different dyes. Recently, a number of studies
reported increased dispersivity of adsorbed as compared to conservative chemicals (e.g., Pot andGenty, 2005), this increase being more pronounced for nonlinearly than linearly adsorbed
chemicals (e.g., Vanderborght et al., in press). The same relations should be evaluated for thedispersivity within pore regions of two-flow region models. As a first indication, the descriptionof IPU and TER breakthrough through the Ap and Ah soil columns improved when k f was
optimized to larger values than were obtained for Br breakthrough.
4.4. Limitations of the approach
Despite the advanced models and the comprehensive data used in this study, there are some
limitations to our approach. To overcome these limitations will require further experimental andmodel refinement. For example even for conservative solute transport, two-region model
parameters may depend on time and length scales, as well as boundary conditions used in theexperiment (Young and Ball, 2000). For reactive nonequilibrium solute transport, the parameters
f f, x c, and l seemed to depend on flow rate, water saturation, and degree of nonequilibrium(Pot et al., 2005). Although the transient flow conditions with wetting and drying cycles used inthis study approximate better the natural flow conditions, fitted parameters may only reflect
effective values averaged over variable saturations and flow velocities. Little is known about
how natural transient flow conditions affect average parameters. Improved experimental
techniques are required to separate real functional dependences of model parameters on flow

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

29

conditions from model artefacts. For examples, incorrect sorption parameter values can beobtained when equilibrium conditions are not fully understood (Sabbah et al., 2005), andanalyzing solute breakthrough curves alone may not be sufficient to separate effects ofdegradation and irreversible sorption (Pot et al., 2005). For better parameter identification,column experiments will need to be augmented with data about local processes within thecolumn. Examples would be flow path visualization and local concentration measurements using
computer tomography and dye tracers (e.g., Vanderborght et al., 2002), or measurements ofadsorbed herbicide concentrations in preferential flow paths and matrix at the end of theexperiment. Furthermore, possibly inadequate first-order water and solute transfer terms may
cause the dependence of optimized transport parameters on time and length scales. A secondorder water transfer term was shown to give higher accuracy for the two-region model
simulation of preferential water flow (e.g., Köhne et al., 2004). Similarly, the improved accuracy
gained by using second-order (advective and) diffusive solute transfer terms on solute transport
under (transient and) steady-state flow should be evaluated.
5. Conclusions
For the aggregated Ap soil column, the BTCs of Br, IPU, and TER started simultaneously
suggesting preferential flow conditions. Peak concentrations decreased in the order of Br, IPU,
and TER, due to increasing sorption in accordance with batch test results. The inverse DPM
simulation of water flow and Br transport revealed HNE and PNE. Inverse DPM and DP-MIM
simulations could apparently approximate the herbicide breakthrough curves when fitting
reduced equilibrium sorption in the preferential flow paths. However, the DP-KM approachrevealed that herbicide transport in the fracture pore region of the Ap soil column was governedby rate-limited sorption. Furthermore, large optimized degradation rates were estimated only for
the Ap soil. However, comparable simulated BTCs were obtained when either the degradation or
sorption rates were optimized. Our results therefore suggest that it may not be possible toinversely identify the degradation rate in nonequilibrium transport models from a BTC if
transport is affected by time-dependent sorption.
The Ah soil column contained crumbs and an earthworm burrow. The BTCs of IPU and Br
featured narrow early peaks characteristic of preferential flow, while the TER breakthroughdisplayed an early but less sharp initial peak. Inverse DPM simulations of water flow and Br
displacement revealed HNE and PNE as for the Ap soil, while batch experiments revealed
instantaneous sorption (no CNE) for IPU and TER in the Ah soil. This was confirmed by themodel analysis of herbicide transport where assuming kinetic sorption did not improve themodel-data agreement. The forward simulation using the DPM (equilibrium sorption) approachgave an acceptable approximation of the IPU BTC. This suggests that in the absence of CNE,
more reliable predictions of preferential herbicide transport can be made based on standard batch
sorption tests and hydraulic and tracer transport data, than in the presence of CNE. As expectedfrom the observed triple-porosity structure (macropore, inter-aggregate region, aggregates) of the
Ah soil column, the DP-MIM was the best approach for simulating the Br BTC, with transport
after two weeks controlled by slow diffusion into the immobile part of the matrix. However, it
was the DPM and not the DP-MIM that provided the best description of the IPU and TER
breakthrough. No obvious reason was found for this discrepancy, other than the possibility thatthe dispersivity may be different for conservative and (nonlinearly) adsorbed solutes.
Overall, modeling our data permitted identification of different processes involved in
herbicide transport through variably saturated structured soil. For the aggregated Ap soil column,

30

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

herbicide transport was under complex HNE–PNE–CNE conditions with rate-limited sorption,
while herbicide transport through the macroporous Ah soil column during pronounced
preferential flow was only under HNE–PNE conditions with equilibrium sorption, which showsthat sorption kinetics in structured soils is not necessarily related to the extent of preferential
flow. However, to describe the herbicide transport in structured soils requires many parameters
for models of increasing complexity. To ensure realistic process description, as many model
parameters as possible should be identified independently, and column experiments need to befurther developed to facilitate parameter identification in complex multi-process flow-transportmodels.
Acknowledgements
We wish to thank Valérie Pot and an anonymous reviewer for helpful comments. The
development of nonequilibrium flow and transport models was partly supported by the
Terrestrial Sciences Program of the Army Research Office (Terrestrial Processes and Landscape
Dynamics and Terrestrial System Modeling and Model Integration).
References
Agritox, 2001. Toxicological Data Base on Pesticides. INRA, France.
Atkins, P.W., 1990. Physical chemistry. 995 p., 4th ed. New York.
Ball, W.P., Roberts, P.V., 1991. Long-term sorption of halogenated organic chemicals by aquifer materials: Part 2.
Intraparticle diffusion. Environ. Sci. Technol. 25 (7), 1237 – 1249.
Berglof, T., Koskinen, W.C., Duffy, M.J., Norberg, K.A., Kylin, H., 2003. Metsulfuron methyl sorption–desorption infield-moist soils. J. Agric. Food Chem. 51 (12), 3598 – 3603.
Brusseau, M.L., Rao, P.S.C., 1989. The influence of sorbate-organic matter interactions on sorption nonequilibrium.
Chemosphere 18, 1691 – 1706.
Brusseau, M.L., Jessup, R.E., Rao, P.S.C., 1989. Modeling the transport of solutes influenced by multiprocess
nonequilibrium. Water Resour. Res. 25, 1971 – 1988.
Cameron, D.R., Klute, A., 1977. Convective-dispersive solute transport with a combined equilibrium and kinetic
adsorption model. Water Resour. Res. 13, 183 – 187.
Celis, R., Cox, L., Hermosin, M.C., Cornejo, J., 1997. Sorption of thiazafluron by iron- and humic acid-coatedmontmorillonite. J. Environ. Qual. 26, 472 – 479.
Chen, S., Franklin, R.E., Johnson, A.D., 1997. Clay film effects on ion transport in soil. Soil Sci. 162 (2), 91 – 96.
Cooke, C.M., Shaw, G., Collins, C.D., 2004. Determination of solid–liquid partition coefficients (K–d) for the herbicides
isoproturon and trifluralin in five UK agricultural soils. Environ. Pollut. 132 (3), 541 – 552.
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.
Dubus, I.G., Brown, C.D., 2002. Sensitivity and first-step uncertainty analysis for preferential flow model MACRO.
J. Environ. 31, 227 – 240.
Flury, M., 1996. Experimental evidence of transport of pesticides through field soils — a review. J. Environ. Qual. 25,
25 – 45.
Gamerdinger, A.P., Wagenet, R.J., van Genuchten, M.T., 1990. Application of two-site/two-region models for studying
simultaneous nonequilibrium transport and degradation of pesticides. Soil Sci. Soc. Am. J. 54, 957 – 963.
Gerke, H.H., Köhne, J.M., 2002. Hydraulic properties of soil aggregate skins. Soil Sci. Soc. Am. J. 66, 26 – 36.
Gerke, H.H., van Genuchten, M.T., 1993. A dual-porosity model for simulating the preferential movement of water andsolutes in structured porous media. Water Resour. Res. 29, 305 – 319.
IVA-Industrieverband Agrar e.V, 1990. Wirkstoffe in Pflanzenschutz und Schädlingbekämpfungsmitteln- Physikalische,chemische und toxikologische Daten: 2. Aufl. BLV Verlags-Gesellschaft, München.
James, R.V., Rubin, J., 1979. Applicability of the local equilibrium assumption to transport through soil of solutes
affected by ion exchange. In: Jenne, E.A. (Ed.), Chemical Modeling of Aqueous Systems. American Chemical
Society, Washington, D.C., pp. 225 – 235.

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

31

James, T.K., Rahman, A., Holland, P.T., McNaughton, D.E., Heiermann, M., 1998. Degradation and movement of
terbuthylazine in soil. Proc. 51st N.Z. Plant Protection Conf., pp. 157 – 161.
Jarvis, N., 1994. The MACRO model (version 3.1)-Technical description and sample simulations. Rep. and Diss., vol.
19. Dept. Soil Sci., Swed. Univ. of Agric. Sci., Uppsala, Sweden 51 pp.
Kladivko, E.J., Brown, L.C., Baker, J.L., 2001. Pesticide transport to subsurface tile drains in humid regions of North
America. Crit. Rev. Environ. Sci. Technol. 31 (1), 1 – 62.
Köhne, J.M., 1999. Analyse präferentiellen Wasserflusses und Stofftransports in strukturierten Böden mit Hilfe eines
Dual-Porositätsmodells. Doctoral thesis (in German, with English summary). In: Roweck, H., Lennartz, B. (Ed.),
Schriftenreihe Inst. Water Management and Landscape Ecology 30, 173 pp. Christian-Albrechts-University, Kiel,
Germany.
Köhne, J.M., Mohanty, B.P., 2005. Water flow processes in a soil column with a cylindrical macropore: experiment andhierarchical modeling. Water Resour. Res. 41, W03010.
Köhne, J.M., Gerke, H.H., Köhne, S., 200a. Effective diffusion coefficients of soil aggregates with surface skins. Soil Sci.
Soc. Am. J. 66 (5), 1430 – 1438.
Köhne, J.M., Köhne, S., Gerke, H.H., 2002b. Estimating the hydraulic functions of dual-permeability models from bulk
soil data. Water Resour. Res. 38 (7).
Köhne, J.M., Mohanty, B.P., Simlnek, J., Gerke, H.H., 2004. Evaluation of a second-order water transfer term forvariably saturated dual-permeability flow models. Water Resour. Res. 40, W07409.
Larsson, M., Jarvis, N., 1999. Evaluation of a dual-porosity model to predict field-scale solute transport in a macroporous
soil. J. Hydrol. 215, 153 – 171.
Legates, D.R., McCabe Jr., G.J., 1999. Evaluating the use of bgoodness-of-fitQ measures in hydrologic and hydroclimatic
model validation. Water Resour. Res. 35 (1), 233 – 241.
Ma, L., Selim, H.M., 1995. Transport of a nonreactive solute in soils: a two-flow domain approach. Soil Sci. 159 (4),224 – 234.
Ma, L., Selim, H.M., 1998. Physical nonequilibrium in soils: modeling and application. In: Selim, H.M., Ma, L. (Eds.),
Physical Nonequilibrium in Soils. Ann Arbor Press, Chelsea, Michigan, pp. 83 – 115.
Mallawatantri, A.P., McConkey, G.B., Mulla, M.J., 1996. Characterization of pesticide sorption and degradation in
macropore linings and soil horizons of Thatuna silt loam. J. Environ. Qual. 25, 227 – 235.
Meyer-Windel, S., 1998. Transport processes of bromide and herbicides under steady-state and transient flow conditions.
PhD-dissertation, Schriftenreihe Institut für Wasserwirtschaft und Landschaftsökologie der Christian-AlbrechtsUniversität zu Kiel, 28. 113 pp.
Meyer-Windel, S., Lennartz, B., Widmoser, P., 1999. Bromide and herbicide transport under steady-state and transient
flow conditions. J. Eur. Soil Sci. 50, 23 – 33.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated soils. Research Project Report,
vol. 442. Technion, Israel Institute of technology, Haifa.
Parrish, R.S., Smith, Ch.N., Fong, F.K., 1992. Tests of the pesticide root zone model and the aggregate model fortransport and transformation of aldicarb, metolachlor, and bromide. J. Environ. Qual. 21, 685 – 697.
Pennell, K.D., Hornsby, A.G., Jessup, R.E., Rao, P.S.C., 1990. Evaluation of five simulation models for predicting
aldicarb and bromide behavior under field conditions. Water Resour. Res. 26 (11), 2679 – 2693.
PerrinGanier, C., Breuzin, C., Portal, J.M., Schiavon, M., 1996. Availability and persistence of isoproturon under field
and laboratory conditions. Ecotoxicol. Environ. Saf. 35 (3), 226 – 230.
Pot, V., Simlnek, J., Benoit, P., Coquet, Y., Yra, A., Matrinez Cordün, M.-J., 2005. Impact of rainfall intensity on thetransport of two herbicides in undisturbed grassed filter strip soil cores. J. Contam. Hydrol. 81 (1–4), 63 – 88.
Pot, V., Genty, A., 2005. Sorbing and non-sorbing solute migration in rough fractures with a multi-species LGA model:
dispersion dependence on retardation and roughness. Transp. Porous Media 59, 175 – 196.
Rambow, J., Lennartz, B., 1993. Laboratory method for studying pesticide dissipation in the vadoze zone. Soil Sci. Soc.
Am. J. 57, 1476 – 1479.
Ray, C., Vogel, T., Dusek, J., 2004. Modeling depth-variant and domain-specific sorption and biodegradation in dualpermeability media. J. Contam. Hydrol. 70, 63 – 87.
Roulier, S., Jarvis, N., 2003. Modeling macropore flow effects on pesticide leaching: inverse parameter estimation using
microlysimeters. J. Environ. Qual. 32, 2341 – 2353.
Roulier, S., Jarvis, N., 2003. Analysis of inverse procedures for estimating parameters controlling macropore flow andsolute transport in the dual-permeability model MACRO. Vadose Zone J. 2, 349 – 357.
Sabbah, I., Ball, W.P., Young, D.E., Bouwer, E.J., 2005. Misinterpretations in the modeling of contaminant
desorption from environmental solids when equilibrium conditions are not fully understood. Environ. Eng. Sci.
22 (3), 350 – 366.

32

J.M. Köhne et al. / Journal of Contaminant Hydrology 85 (2006) 1–32

Selim, H.M., Ma, L., Zhu, H., 1999. Predicting solute transport in soils: second-order two-site models. Soil Sci. Soc. Am.
J. 63, 768 – 777.
Šimlnek, J., Šejna, M., van Genuchten, M.T., 1998. The HYDRUS-1D Software Package for Simulating the OneDimensional Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 2.0, IGWMCTPS-70. International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado. 202 pp.
Šimlnek, J., Wendroth, O., Wypler, O., van Genuchten, M.T., 2001. Nonequilibrium water flow characterized from an
upward infiltration experiment. Eur. J. Soil Sci. 52 (1), 13 – 24.
Šimlnek, J., Jarvis, N.J., van Genuchten, M.T., Gärdenäs, A., 2003. Review and comparison of models for describing
nonequilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272, 14 – 35.
Skopp, J., Gardner, W.R., Tyler, E.J., 1981. Solute movement in structured soils: two-region model with smallinteraction. Soil Sci. Soc. Am. J. 45, 837 – 842.
Smith, M.C, Shaw, D.R., Massey, J.H., Boyette, M., Kingerey, W., 2003. Using nonequilibrium thin-disc and batch
equilibrium techniques to evaluate herbicide sorption. J. Environ. Qual. 32 (4), 1393 – 1404.
Stehouwer, R.C., Dick, W.A., Traina, S.J., 1993. Characteristics of earthworm burrow lining affecting atrazine sorption.
J. Environ. Qual. 22, 181 – 185.
Stehouwer, R.C., Dick, W.A., Traina, S.J., 1994. Sorption and retention in vertically oriented earthworm and artificial
burrows. J. Environ. Qual. 23, 286 – 292.
Streck, T., Poletika, N.N., Jury, W.A., Farmer, W.J., 1995. Description of simazine transport with rate-limited, two-stage,
linear and non-linear sorption. Water Resour. Res. 31, 811 – 822.
U.S. Environmental Protection Agency (USEPA), 1993. Pesticides in Ground Water Database—A Compilation of
Monitoring Studies: 1971–1991 National Summary. EPA 734-12-92-001. US EPA Office of Pesticide Programs,
Washington, DC.
Valocchi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing solute transport throughhomogeneous soils. Water Resour. Res. 21, 808 – 820.
Vanderborght, J., Gähwiller, P., Flühler, H., 2002. Identification of transport processes in soil cores using fluorescent
tracers. Soil Sci. Soc. Am. J. 66, 774 – 787.
Vanderborght, J., Kasteel, R., Vereecken, H., in press. Stochastic continuum transport equations for field-scale solutetransport: overview of theoretical and experimental results. Vadose Zone J. 5 (2). doi:10.2136/vzj2005.0024.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J. 44, 892 – 898.
van Genuchten, M.T., Wagenet, R.J., 1989. Two-site/two-region models for pesticide transport and degradation:
theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303 – 1310.
van Genuchten, M.T., Wierenga, P.J., 1976. Mass transfer studies in sorbing porous media: I. Analytical solutions.
Soil Sci. Soc. Am. J. 40, 473 – 481.
Veeh, R.H., Inskeep, W.P., Roe, F.L., Ferguson, A.H., 1994. Transport of chlorsulfuron through soil columns. J. Environ.
Qual. 23, 542 – 549.
Wagenet, R.J., Chen, W., 1998. Coupling sorption rate heterogeneity and physical nonequilibrium in soils. In: Selim,
H.M., Ma, L. (Eds.), Physical Nonequilibrium in Soils. Ann Arbor Press, Chelsea, Michigan, pp. 1 – 36.
Walker, A., Bromilow, R.H., Nicholls, P.H., Evans, A.A., Smith, V.J.R., 2002. Spatial variability in the degradation rates
of isoproturon and chlorotoluron in a clay soil. Weed Res. 42, 39 – 44.
Young, D.F., Ball, W.P., 2000. Column experimental design requirements for estimating model parameters from temporal
moments under nonequilibrium conditions. Adv. Water Resour. 23 (5), 449 – 460.