Journal of Contaminant Hydrology 81 (2005) 63 – 88www.elsevier.com/locate/jconhyd

Impact of rainfall intensity on the transportof two herbicides in undisturbed grassed
filter strip soil cores
V. Pot a,*, J. Šimlnek b, P. Benoit a, Y. Coquet a,A. Yra a,c, M.-J. Martı́nez-Cordón a,da

UMR INRA-INA P-G, Unité Environnement et Grandes Cultures, BP01, 78850 Thiverval-Grignon, Franceb
Department of Environmental Sciences, University of California, Riverside, California, USAc
UMR CNRS 8508-Laboratoire Transferts Ecoulements Fluides Energétiques, site de l’ENSAM,
33405 Talence Cedex, FrancedDepartamento de Quı́mica Agrı́cola y Edafologı́a, Campus Universitario de Rabanales, 14071 Córdoba, Spain
Received 9 July 2004; received in revised form 13 June 2005; accepted 19 June 2005
Available online 16 September 2005

Abstract
Two series of displacement experiments with isoproturon and metribuzin herbicides were
performed on two undisturbed grassed filter strip soil cores, under unsaturated steady-state flow
conditions. Several rainfall intensities (0.070, 0.147, 0.161, 0.308 and 0.326 cm h 1) were used. A
water tracer (bromide) was simultaneously injected in each displacement experiment. A descriptive
analysis of experimental breakthrough curves of bromide and herbicides combined with a modeling
analysis showed an impact of rainfall intensity on the solute transport. Two contrasting physical nonequilibrium transport processes occurred. Multiple (three) porosity domains contributed to flow at
the highest rainfall intensities, including preferential flow through macropore pathways. Macropores
were not active any longer at intermediate and lowest velocities, and the observed preferential
transport was described using dual-porosity-type models with a zero or low flow in the matrix
domain. Chemical non-equilibrium transport of herbicides was found at all rainfall intensities.
Significantly higher estimated values of degradation rate parameters as compared to batch data werecorrelated with the degree of non-equilibrium sorption. Experimental breakthrough curves wereanalyzed using different physical and chemical equilibrium and non-equilibrium transport models:* Corresponding author. Tel.: +33 1 30 81 54 02; fax: +33 1 30 81 53 96.
E-mail address: vpot@grignon.inra.fr (V. Pot).
0169-7722/$ - see front matter D 2005 Elsevier B.V. All rights reserved.
doi:10.1016/j.jconhyd.2005.06.013

64

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

convective–dispersive model (CDE), dual-porosity model (MIM), dual-permeability model (DP),triple-porosity, dual permeability model (DP-MIM); each combined with both chemical instantaneous and kinetic sorption.
D 2005 Elsevier B.V. All rights reserved.
Keywords: Pesticide transport; Nonequilibrium transport; Undisturbed soil columns; Breakthrough curves;
Numerical model; HYDRUS-1D

1. Introduction
Grassed buffer strips were introduced as one of the best management practices toreduce soil losses due to erosion of cultivated soils. Located downslope of cropped fields
and adjacent to streams or drainage ditches, grassed buffer strips are effective sedimenttraps and, moreover, have been shown to effectively remove soluble and sediment-boundchemicals from cropland runoff (Dillaha et al., 1989; Arora et al., 1996). Reviewing
published data, Baker et al. (1995) found that herbicide runoff was reduced by an averageof 48% by different types of buffer systems with reductions ranging from 9% to 91%.
Many parameters such as soil type, soil water content, runoff volume, buffer width, andbuffer vegetation can explain this variability (Patty et al., 1997). Vertical infiltration
through the grass is the major reduction mechanism (Mersie et al., 1999). Physicalfiltration by the grass and sedimentation of the particle-bound pesticides at the soil surface
further reduce the pesticide concentration in runoff water. The high potential for sorption
and mineralization of pesticides in surface layers of the grassed soil has been pointed out
in different studies (Hurto et al., 1979; Benoit et al., 1999). It is related to the high organic
matter content of the topsoil, originating from the high biological activity in grassland soils(Merckx et al., 1985; Haynes and Francis, 1993). On the other hand, root growth andearthworm activity can lead to development of preferential pathways for soluble andsediment-bound pesticides that will then by-pass the surface horizon where pesticides canbe retained and degraded. Therefore, there is a potential risk of pesticide contamination ofshallow groundwater under grassed buffer strips. Field studies have reported preferential
transport of herbicides through grassed buffer strips (Kanchanasut and Scotter, 1982;
Delphin and Chapot, 2001). Runoff simulations performed at the local scale (i.e. 3 m2) andtension infiltrometer measurements at a smaller scale (50 cm2) showed the importance of
infiltration in the reduction of pesticide runoff, but could not identify the particularmechanisms involved within the soil. Lateral sub-surface flow or vertical preferential
transport through macropores were possible explanations (Souiller et al., 2002).
Although the high sorption capacities and mineralization rates in the grassed strips arewell documented, the fate of pesticides inside the grassed buffer strip is still poorly
understood. Especially, the role of infiltration needs to be better quantified. Some studies
about vertical infiltration of pesticides in grassed buffer strip soils have been performed in
the laboratory, using column displacement experiments. Mersie et al. (1999) showed thatswitchgrass roots caused preferential movement of bromide and atrazine herbicide inrepacked soil columns. Reungsang et al. (2001) observed preferential transport of bromide
and atrazine in undisturbed soil columns, although they found little difference in transport

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

65

characteristics between the switchgrass and adjacent cultivated soils. Other authors foundcontrasting results. Vogeler et al. (1997) did not observed preferential flow of anions inundisturbed soil columns sampled under pasture. Benoit et al. (2000) showed using
undisturbed soil columns that leaching of isoproturon herbicide was significantly higher in
the adjacent cropped soil compared to the grassed soil. They measured an enhanced
isoproturon degradation due to the formation of nonextractable (bound) residues in thegrassed soil. Leaching also depended on isoproturon aging in the soil before percolation in
the columns.
This study focuses on the identification of the main flow and transport mechanisms in a
grassed buffer strip soil and on the fate of two herbicides, isoproturon and metribuzin. A
series of displacement experiments was performed at increasing rainfall intensities on two
undisturbed soil columns, sampled from the topsoil of a grassed buffer strip. The modified
one-dimensional transport model HYDRUS-1D (Šimlnek et al., 1998, 2003) was used to
identify and quantify different physical and chemical non-equilibrium transport processes
acting in the grassed soil.

2. Material and methods
2.1. Experimental
2.1.1. Soil description and core sampling
The studied grassed buffer strip is located at the experimental site of ARVALIS at La
Jaillière in the western part of France (Loire-Atlantique). The 20-m-wide grassed strip,
planted with perennial rye-grass (Lolium perenne) in 1992, is cut twice a year. The soil is ahydromorphic silt loam (24% clay, 40% silt, 36% sand) developed on altered schists of theMassif Armoricain. A perched aquifer often forms during the winter period at a depth of
25–60 cm and can occasionally reach the soil surface (Patty, 1997).
Two undisturbed soil cores (14-cm i.d. and 30-cm height) were sampled in June 2002
from the soil surface layer. To minimize soil compaction, a hydraulic pressure device wasused to extract the cores. Bulk densities were sampled in the immediate vicinity. The
arithmetic mean of the bulk density over the 0–30 cm depth was 1.25 g cm 3.
2.1.2. Displacement experiments
Displacement experiments were conducted at steady-state flow under unsaturated
conditions.
A CaCl2 solution (5.10 3 M CaCl2) was supplied to the soil surface using a rainfall
simulator consisting of 47 hypodermic needles (0.3-mm i.d.) at a 25-mm spacing. A
membrane pump delivered a constant flux to the rainfall simulator. A stainless steel
mesh (25 Am porosity) was placed at the soil base and covered with glass beads (95–150 Am diameter) to establish a full contact between the soil base and the mesh. To
maintain unsaturated conditions, a negative pressure of  15 cm was applied at themesh level by fixing the height between the base of the core and the constant water level
outlet of the column apparatus. Two time domain reflectometry (TDR) probes (8-cm
long uncoated mini-probes) were horizontally inserted 6 and 20 cm below the soil

66

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

surface to measure the volumetric water content inside the soil columns. Threetensiometers (2-cm long and 0.6-cm diameter) were located at 6, 15 and 25 cm below
the soil surface to measure the soil water pressure head. The effluent was collected atfixed time intervals of 0.5 h at the highest velocity and 1.5 h at the intermediate andlowest velocities to obtain the breakthrough curves. Note that the length of the tubing
below the base of the core was about 80 cm long and 4.8-mm i.d. The delay in thetubing was taken into account when the experimental breakthrough curves wereconstructed by correcting the arrival time of the effluents.
The displacement experiments consisted of injecting a pulse of the CaCl2 solution
containing actual concentrations of 66.1 F 2.1 mg l 1 Br, 7.6 F 1.5 mg l 1 isoproturon(IPU) and 8.5 F 0.2 mg l 1 metribuzin (Met). Three and two successive displacement
experiments were performed at different rainfall intensities for columns I and II,
respectively. The solute pulse duration was fixed so that an approximately equal massof Br, IPU and Met was injected at each velocity. Differences in flow rates between thepulse solution application and the rainfall were less than 10% for all experiments. Between
the successive solute pulses, the columns were thoroughly rinsed with the CaCl2 solution
for several days until no solutes were detected in the effluent.
2.1.3. Herbicides properties
Isoproturon and metribuzin were purchased from CIL Cluzeau (France). Water
solubilities for IPU and Met are 70 and 1200 mg l 1; K oc values range from 36 to 241
l kg 1 and from 9 to 95 l kg 1, and DT50s (half-lives) range from 12 to 33 days and from23 to 180 days, respectively (Agritox, 2001). Met therefore shows higher mobility andpersistence in soils than IPU. Madrigal et al. (2002) measured the sorption properties of
IPU in the first 60-cm of the grassed strip at La Jaillière. The distribution coefficient, K d,
showed a strong vertical gradient (decreasing with the soil depth) that was correlated withthe gradient of organic matter (also decreasing with the soil depth). Sorption anddegradation properties were also measured by Benoit et al. (1999) in the same grassed
strip. For IPU, a mean distribution coefficient was found to be equal to 2.4 l kg 1, a mean
K oc 143 l kg 1, a mean degradation rate 0.001203 h 1, and a mean DT50 24 days, when
considering the 0–30 cm horizon as one single layer.
2.1.4. Chemical analysis
Effluent was filtered on polyvinylidene difluoride syringe filters (0.45 Am) and stored at
4 8C until analysis. Bromide and herbicide concentrations were measured using a highperformance liquid chromatography (Waters HPLC appliance) with an Allsep A-2 anion
column for Br and a Nova Pack silica-C18 column (Waters) for IPU and Met. UV
detection (Waters 996 Photodiode Array Detector) was performed at 194, 240 and 298 nm
for Br, IPU and Met, respectively. Detection limits were 1 mg l 1 for Br and 0.001 to 0.01
mg l 1 for the herbicides.
2.2. Transport models
We used the HYDRUS-1D model that simulates water flow and transport ofconservative and reactive solutes in soils (Šimlnek et al., 1998, 2003). The latest version

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

67

of the code implements different equilibrium and non-equilibrium flow and transportmodels in both direct and inverse mode. The following transport processes can be used
to analyze experimental breakthrough curves (BTCs): a) physical and chemical
equilibrium transport described using the classical convection–dispersion transportmodel (Lapidus and Amundson, 1952); b) chemical non-equilibrium transport using thetwo-site chemical model (Cameron and Klute, 1977); c) physical non-equilibriumtransport using the two region model (mobile and immobile water) (van Genuchten and
Wierenga, 1976) and the dual-permeability model (Gerke and van Genuchten, 1993);
and finally, e) simultaneous physical and chemical non-equilibrium transport using thedual-permeability model (Gerke and van Genuchten, 1993) combined with the two-sitechemical model.
Although we describe below only models that use the linear sorption isotherm,
HYDRUS-1D can also consider nonlinear sorption using Freundlich or Langmuir sorption
isotherms. Similarly, HYDRUS-1D considers transient flow and transport processes(Šimlnek et al., 2001), while below we describe only models for steady-state water flow
condition.
2.2.1. Equilibrium transport model — CDE
Assuming that decay is active only in the liquid phase (assumed in all models
described below) and linear relationship between liquid and solid phase concentrations:s = K d c, where c (ML 3) and s (MM 1) are concentrations in the liquid and solid
phases, respectively, and K d is an empirical distribution constant (L3M 1), thetransport of a reactive solute for steady-state water flow conditions can be written as
follows:
hR

Bc
B2 c
Bc
¼ hD 2  q
 lw hc
Bt
BZ
Bz

ð1Þ

where h is the volumetric water content (L3L 3), D is the dispersion coefficient(L2T 1), q is the Darcian water flux (LT 1), l w is first-order degradation constant forthe solute in the liquid phase (T 1), t (T) and z (L) are the temporal and spatial
coordinates, respectively, and R is the retardation factor ():
R¼1þ

qKd
h

ð2Þ

where q is the bulk density (ML 3).
2.2.2. Chemical non-equilibrium transport model
The concept of two-site sorption (Selim et al., 1977; van Genuchten and Wagenet,1989) permits consideration of non-equilibrium adsorption–desorption reactions. The twosite sorption concept assumes that the sorption sites can be divided into two fractions:adsorbed concentration, s e (MM 1), is associated with one fraction of the sites (the type-1
sites) where sorption is assumed to be instantaneous, while adsorbed concentration, s k(MM 1), is associated with the remaining (type-2) sites where sorption is considered to be

68

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

a time-dependent, first-order kinetic process. The governing convection–dispersion
equation is then written as follows:
hRe

Bc
B2 c
Bc
¼ hD 2  q
 lw hc  qach ½ð1  fe ÞKd c  sk 
Bt
Bz
Bz

ð3Þ

where f e is the fraction of exchange sites assumed to be in equilibrium with the solution
phase (), a ch is the first-order mass transfer rate constant (T 1) for chemical nonequilibrium and R e is the retardation factor for the equilibrium phase ():Re ¼ 1 þ

qfe Kd:h

ð4Þ

2.2.3. Physical non-equilibrium transport models
2.2.3.1. Dual-porosity (mobile–immobile water) model — MIM. The concept of tworegion, dual-porosity type solute transport (van Genuchten and Wierenga, 1976) permits
consideration of physical non-equilibrium transport (MIM model). The two-region
concept assumes that the liquid phase can be partitioned into mobile (flowing), h m(L3L 3), and immobile (stagnant), h im (L3L 3), flow regions and that solute exchange
between the two liquid regions can be modeled as a first-order process. The governing
convection–dispersion equation is then written as follows:h m Rm

Bcm
B2 cm
Bcm
¼ hm D m
 lw hm cm  aph ðcm  cim Þ
q
2
Bt
Bz
Bz

ð5Þ

where the subscript m and im refer to mobile and immobile regions, respectively; a ph is themass transfer rate coefficient (T 1) for physical non-equilibrium; and R m is the retardation
factor in the mobile zone ():Rm ¼ 1 þ

qfm Kdhm

ð6Þwhere f m is the fraction of sorption sites in contact with mobile region.
2.2.3.2. Dual-permeability model — DP. The dual-permeability (DP) model (Gerke andvan Genuchten, 1993; Šimlnek et al., 2003) assumes that the porous medium is composedof two mobile flow regions, each with its own flow and transport properties. The tworegions are often alternatively referred to as matrix and fracture regions, micropores andmacropores, or intra- and inter-porosity. Contrary to the dual-porosity model, dualpermeability model assumes that flow can occur in both domains. The total water content
and water flux density are then defined as follows:
h ¼ wf hf þ ð1  wf Þhmaq ¼ wf qf þ ð1  wf Þqma

ð7Þ
3 3

where h ma and h f are volumetric water contents in the two pore regions (L L ), q ma andq f are corresponding water fluxes (LT 1), subscripts ma and f refer to the matrix andV. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

69

fracture regions, and w f is the volumetric weighting factor () expressing the ratio of thefracture domain to total volume.
The governing convection–dispersion equations for the matrix and fracture regions are
written as follows:
hma Rma

hf R fBcma
B2 cmaBcma
Cs¼ hma Dma lw hma cma þ
 qma
Bt
Bz2
Bz
1  wf

B2 cf
Bcf
Bcf
Cs
¼ hf Df
 lw h f c f 
 qf
2
Bt
Bz
Bz
wf

ð8Þ

ð9Þ

where the retardation coefficients are written as:
Rma ¼ 1 þ

Rf ¼ 1 þ

qKd
hma

qKd
:
hf

ð10Þ

Although in general reaction coefficients may be different in the two regions, in theformulation above we assume that both reaction coefficients (K d and l w) and bulk
densities are the same in both regions.
The transfer rate, C s, for solutes between the fracture and matrix regions (ML 3 T 1)
is usually given as the sum of diffusive and convective fluxes. For steady-state water flow
condition convective fluxes are equal to zero and C s can thus be written as (Gerke and van
Genuchten, 1996):
Cs ¼ as ð1  wf Þðcf  cma Þ
where a s is the first-order solute mass transfer rate coefficient (T
as ¼

b
Da
a2

ð11Þ
1

) of the form:
ð12Þ

in which D a is an effective diffusion coefficient (L2T 1), which represents the diffusion
properties of the fracture–matrix interface as well as other parameters, b is a shape factor
that depends on the geometry (the value of b ranges from 3 for rectangular slabs to 15 forspherical aggregates), and a is the characteristic length of the aggregate (L) (e.g. the radiusof a spherical or solid cylindrical aggregate, or half the width of a rectangular aggregate).
In the following, we assumed spherical aggregates (b = 15) and a = 0.1 cm as thecharacteristic length of the aggregate (Kay, 1990).
2.2.3.3. Dual-permeability mobile–immobile water model — DP-MIM. Dual-permeability
mobile–immobile water model (DP-MIM) assumes that the liquid phase in the matrix region
can be further partitioned into mobile (flowing), h m,ma (L3L 3), and immobile (stagnant),
h im,ma (L3L 3), regions:hma ¼ hm;ma þ him;ma

ð13Þ

70

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

where h ma is the volumetric water content of the matrix pore system (L3L 3). The governing
convection–dispersion equation for transport in the matrix region (8) is then replaced withmodified Eq. (5):hm;ma Rm;ma

Bcm;maB2 cm;ma
Bcm;ma¼ hm;ma Dm;ma lw hm;ma cm;ma
 qma
Bt
Bz2
Bz

Cs
þ
 ama cm;ma  cim;ma
1  wfð14Þ

where c im,ma and c m,ma are solute concentrations in the immobile and mobile zones of thematrix region (ML 3), respectively; and a ma is the mass transfer rate coefficient between
mobile and immobile zones of the matrix region (T 1). The solute exchange between the two
liquid regions in the matrix is modeled as a first-order process, i.e.,him;ma Rim;ma
Bcim;ma¼ ama cm;ma  cim;ma  lw him;ma cim;ma
Bt

ð15Þ

with retardation coefficients R m,ma and R im,ma for the mobile and immobile domains of thematrix zone, respectively, given as:
Rm;ma ¼ 1 þqfm;ma Kdhm;ma
q 1  fm;ma Kd
Rim;ma ¼ 1 þhim;mað16Þ

where f m,ma is the fraction of sorption sites in contact with mobile region in the matrix.
2.2.4. Physical and chemical non-equilibrium model
Simultaneous physical and chemical non-equilibrium processes are implemented in
HYDRUS-1D using the dual-permeability model (Gerke and van Genuchten, 1993;
Šimlnek et al., 2003). The governing convection–dispersion equations for the two regions
are written as follows:
hma Rma

Bcma
B2 cmaBcma
Cs¼ hma Dma lw hma cma þ
 qma
Bt
Bz2
Bz
1  wf qach ½ð1  fma ÞKd cma  skma 

hf Rf

Bcf
B 2 cfBcf
Cs¼ hf D f
 lw h f c f  qf
 qach ½ð1  ff ÞKd cf  skf 
Bt
Bz2
Bz
wf

ð17Þ
ð18Þ

where the adsorbed concentrations, s kma and s kf are associated with the type-2 sites in thematrix and fracture domains, respectively; f ma and f f are fractions of exchange sites
assumed to be in equilibrium with the solution phase () in the matrix and fracture

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

71

domains, respectively; and a ch is again the first-order mass transfer rate constant [T 1].
The retardation coefficients are written again as:
Rma ¼ 1 þ

Rf ¼ 1 þ

qfma Kdhma

qff Kd:hf

ð19ÞNote that in all formulations above we assume that the bulk density, reaction (K d andl w) and transfer (a) coefficients are the same in both regions. The dual permeability
model, however, allows for different dispersivities in the two pore regions.

3. Results
3.1. Descriptive analysis
3.1.1. Flow
A satisfactory steady-state flux was obtained for all cases (Table 1). For column I, themean water content, derived from TDR measurements, slightly decreased by nearly 2% asrainfall intensity decreased, while for column II water contents remained about the samefor both fluxes. For both columns, the water content in the upper part of the column wasalways slightly higher compared to the lower part. Pressure head measurements showed

Table 1
Measured flow properties and leaching chemical parametersq(cm h 1)h(cm3 cm 3)

T0
()

Br
M
(%)

IPU
T min
()

Column I
q max 0.308 F 0.016 0.3980 F 0.0212 0.1548 82.6a 0.03
q mid 0.147 F 0.059 0.3865 F 0.0276 0.1521 95.6
q min 0.070 F 0.018 0.3766 F 0.0086 0.168 98.9

0.13
0.09

Column II
q max 0.326 F 0.007 0.3745 F 0.0078 0.1741 81.6

0.05

q mid 0.161 F 0.003 0.3750 F 0.0092 0.1712 83.3f 0.13

Met

T peak
()

M
(%)

T min
()

T peak
()

M
(%)

T min
()

T peak
()

0.10
0.21
0.49
0.60

50.4b 0.03

0.10
0.21
0.75
0.92

75.9b 0.03

0.10
0.21
n.d.
0.80

0.20
0.41
0.51

35.4c 0.32
13.0d 0.14

67.5e 0.06
58.2g 0.02

0.23
0.41
0.26

n.d. n.d.
48.5d 0.09

75.8e 0.05
49.7g 0.09

0.23
0.41
0.39

q is the Darcian water flux measured from the effluent weights, q max, q mid and q min correspond to the highest,
intermediate and lowest rainfall intensities, respectively, h is the volumetric water content, T 0 is the pulse
duration, M is mass recovery, T min is the time of the first appearance of the chemical in the effluent, and T peak is
the time when the maximum chemical concentration was reached in the effluent. T 0, T min and T peak are expressedin pore volumes. Note that at q max, the arrival times of two peaks are presented. Measured after 2.8a, 10.5b, 7.9c,
5.3d, 10.3e, 4f, and 6.6g pore volumes. n.d. is not determined. Standard errors are given after the Fsign.

72

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

also slightly higher pressure heads (between 0 and  5 cm) close to the soil surface,compared to the other two locations (pressure heads ranged from  5 to  15 cm).
Observed moisture gradients can be partially explained by the pressure head distribution
and by the vertical variability of the hydraulic conductivity in this soil. Souiller et al.
(2002) showed using in situ infiltrometer measurements that the surface layer (0–5 cm) has
a higher hydraulic conductivity compared to the underlying layer (a factor from 3 up to 25
was observed) and they hypothesized that the vertical infiltration capacity of the grassed
strip for conditions close to saturation was controlled by the hydrodynamic properties of
the underlying layer. The hypothesis of a solely gravitational flow in a homogeneous soil
was therefore not fully achieved in our column displacement experiments. However,considering that the water content did not change much with depth (less than 2%) when
the flux rate was increased, the hypothesis of a gravity-driven flow was used as a first
approximation.
3.1.2. Bromide and herbicides BTCs
The BTCs are plotted as the measured relative concentration (ratio of the effluent
concentration, c, to the influent concentration, c 0) versus the number of pore volumes
eluted (flux-averaged volume of solution eluted, V, per the column pore volume, V 0) (Figs.
1–3). Due to chemical analysis problems, the metribuzin BTC was not available in the caseof the intermediate velocity of column I. From now on, for the sake of simplicity, we will
use the following notation for the effluent velocities: the highest, intermediate and lowest
rainfall intensities are, respectively, q max, q mid and q min as defined in Table 1. Since for all
experiments we used approximately the same pulse duration, a qualitative comparison
between particular BTCs can be easily made.
For both columns, the Br BTCs for the q max case were markedly shifted to the left,
exhibited double fronts or peaks and a significant asymmetry with a long tailing elution
curve compared to the other velocities (Fig. 1). Arrival times of the two peaks were
notably different for the two columns: at 0.1 and 0.21 pore volumes for column I with arelative concentration of 0.04 and 0.22, respectively, and at 0.23 and 0.41 pore volumes forcolumn II with a relative concentration of 0.17 and 0.14, respectively. At q mid and q min,comparable peak positions of the BTCs were observed (Table 1). Early arrivals of thepeaks at q max indicated that a rapid transport of the tracer caused by a possible by-pass ofthe main part of the soil porosity and the use of macropore pathways probably occurred in
the grassed strip soil at the largest rainfall intensity. Visual observation of several large
macropores made by earthworm activity (of about 3 to 5 mm in diameter) and numerous
root channels at the bottom of the soil columns (Fig. 4) suggested potential candidates for
macroporous pathways.
At q max, the two herbicides elution curves displayed identical peak positions as thosefor Br (Figs. 2ad and 3ac and Table 1). No delay originating from soil-herbicide sorption
was thus observed. These results supported the observation that preferential flow throughmacropores dominated at this high rainfall intensity.
At q mid and q min, the comparison of arrival times of the concentration peaks of IPU andMet showed a delay as compared to Br for column I (Figs. 2bc and 3b), while earlier
concentration peak arrivals than Br were found for column II (Figs. 2e and 3d). Note that a
non-zero IPU concentration was found in the first effluent for column II at q mid.

C/C0 (-)

0.10

0.25

Column II - Br - q=0.326 cm/h

0.20

Experimental
CDE
MIM
DP
DP-MIM

0.05
0.00

0.00
0.01

0.05

0.50

5.00

0.01

0.05

0.50

5.00

0.10

C/C0 (-)

0.15

0.25
0.20

(b)

Column II - Br - q=0.161 cm/h
Experimental
CDE
MIM
DP

(e)

0.15

Experimental
CDE
MIM
DP

0.00

0.00

0.05

0.10

0.25

Column I - Br - q=0.147 cm/h

0.20

V/V0 (-)

0.05

C/C0 (-)

V/V0 (-)

0.01

0.05

0.50

5.00

V/V0 (-)

0.01

0.05

0.50

5.00

0.25

Column I - Br - q=0.070 cm/h

0.20

V/V0 (-)

Experimental
CDE
MIM
DP

0.10

0.15

(c)

0.00

0.05

C/C0 (-)

(d)

0.15

0.20
0.15

(a)

73

0.10

0.25

Column I - Br - q=0.308 cm/h
Experimental
CDE
MIM
DP
DP-MIM

0.05

C/C0 (-)

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

0.01

0.05

0.50

5.00

V/V0 (-)

Fig. 1. Bromide BTCs measured and calculated for columns I and II using four different physical transportmodels. Symbols are the experimental data, dotted lines, dashed lines, dashed-and-dotted lines and long-dashedlines are the fits given by the CDE, MIM, DP and DP-MIM transport models, respectively. Relative
concentrations of the effluents are presented against the number of pore volumes.

74

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

Column II - IPU - q=0.326 cm/h

Column I - IPU - q=0.308 cm/h
0.12

(d)

0.08

Experimental
DP-MIM - K d
DP-(K d ,αch)

0.04

C/C0 (-)

0.08

Experimental
DP-MIM - K d
DP-(K d ,αch)

0.00

0.00

0.04

C/C0 (-)

0.12

(a)

0.01

0.05

0.50

0.01

5.00

0.05

0.05

Column II - IPU - q=0.161 cm/h

0.04

(e)
Experimental
DP-Kd
DP-(K d ,αch)

0.01

0.02

C/C0 (-)

0.03

Experimental
DP-Kd
DP-(K d ,αch)

0.03

0.04

(b)

0.00

0.00

0.01

C/C0 (-)

5.00

0.02

0.05

Column I - IPU - q=0.147 cm/h

0.01

0.05

0.50

5.00

0.05

0.01

0.05

0.50

5.00

V/V0 (-)

V/V0 (-)

Column I - IPU - q=0.070 cm/h

0.04

(c)

0.02

0.03

Experimental
DP-Kd
DP-(K d ,αch)

0.00

0.01

C/C0 (-)

0.50
V/V0 (-)

V/V0 (-)

0.01

0.05

0.50

5.00

V/V0 (-)

Fig. 2. Comparison of measured and calculated isoproturon concentrations for columns I and II. The DP-MIM
transport model was combined with the equilibrium chemical model (instantaneous sorption) for the highest
rainfall intensities (DP-MIM-K d transport model) and the DP transport model was combined with either theequilibrium (DP-K d transport model) or non-equilibrium chemical models (instantaneous and first-order kinetic
sorption) (DP-(K d, a ch) transport model) for the other rainfall intensities. Symbols represent the experimental
data, dotted and dashed lines the DP-MIM-K d or DP-K d and DP-(K d, a ch) transport models, respectively. Relative
concentrations of the effluents are presented against the number of pore volumes.

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

Column II - Met - q=0.326 cm/h

Column I - Met - q=0.308 cm/h
0.12

(c)

0.08

Experimental
DP-MIM - Kd
DP-(K d ,αch)

0.04

C/C0 (-)

0.08

Experimental
DP-MIM - Kd
DP-(K d ,αch)

0.00

0.00

0.04

C/C0 (-)

0.12

(a)

0.01

0.05

0.50

0.01

5.00

0.05

5.00

0.05

Column II - Met - q=0.161 cm/h

(d)

0.04

(b)

0.03

Experimental
DP-Kd
DP-(K d ,αch)

0.00

0.00

0.01

0.01

0.02

C/C0 (-)

0.03

Experimental
DP-Kd
DP-(K d ,αch)

0.02

0.05

Column I - Met - q=0.070 cm/h

0.04

0.50
V/V0 (-)

V/V0 (-)

C/C0 (-)

75

0.01

0.05

0.50
V/V0 (-)

5.00

0.01

0.05

0.50

5.00

V/V0 (-)

Fig. 3. Comparison of measured and calculated metribuzin concentrations for columns I and II. The DP-MIM
transport model was combined with the equilibrium chemical model (instantaneous sorption) for the highest
rainfall intensities (DP-MIM-K d transport model) and the DP transport model was combined with either theequilibrium (DP-K d transport model) or non-equilibrium chemical models (instantaneous and first-order kinetic
sorption) (DP-(K d, a ch) transport model) for the other rainfall intensities. Symbols represent the experimental
data, dotted and dashed lines the DP-MIM-K d or DP-K d and DP-(K d, a ch) transport models, respectively. Relative
concentrations of the effluents are presented against the number of pore volumes.

For both herbicides, a decreasing mass recovery rate was found with decreasing
velocity. This result was expected since the residence time of herbicides increased as
velocity decreased and thus more time was left for sorption and degradation to occur inside
the soil column. No degradation products were detected in the effluent. Since we had noexperimental data about the fate of herbicides inside the soil columns, it was difficult tomake conclusions about the potential decay processes in the soil column. Degradation withthe formation of metabolites of lower mobility than the parent compound, mineralization,or irreversible sorption of either the parent compound or metabolites could have occurredin our soil columns. Benoit et al. (2000) showed in previous IPU percolation studies on thesame grassed soil (i) the existence of metabolites in the soil, (ii) an increase of nonextractable C-14 residues suggesting that metabolites could be irreversibly sorbed, and (iii)
an increase of K d with time suggesting that IPU could be irreversibly sorbed. The

76

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

Fig. 4. Photo of the soil structure at the bottom of column I.

comparison of mass recovery rates between IPU and Met showed a higher leaching of Met
in column I, while it showed contrasting results for column II with a higher leaching of
Met at q max but a lower leaching of Met at q mid. Except for the latter case, these results arein good agreement with the literature data (Agritox, 2001) on these two herbicides thatsuggests a higher mobility for Met than for IPU.
3.2. Numerical model analyses
Physical transport parameters were first estimated by inverting the Br BTCs for bothcolumns using the CDE, MIM, DP and/or DP-MIM transport models. Then, chemicalparameters for IPU and Met were estimated by inverting the herbicide BTCs using thephysical transport model that best described Br transport combined with equilibrium or
non-equilibrium chemical transport models. Note that non-equilibrium chemical transport
is currently not yet available in the HYDRUS-1D model for the physical non-equilibrium
DP-MIM transport model. The soil profile was assumed to be initially solute free. A
prescribed solute flux (third-type boundary condition) and a zero concentration gradientwere used as the upper and lower boundary conditions, respectively. Various models were
calibrated against measured BTCs using flux concentrations at the lower boundary. For all
simulations, the average pore water velocity, v = q / h, was fixed.
3.2.1. Physical transport parameters (analyses of the Br BTCs)
Results obtained by parameter optimization are presented in Fig. 1 and Table 2. The
molecular diffusion coefficient of Br was taken as 0.066 cm2 h 1 (Weast, 1985). The
following physical transport parameters were fitted for different models: D and R for CDE;
the dispersivity k m of the mobile phase (L), h im, and a ph for MIM; q ma, k f, w f, and D a forDP, and q ma, k ma, k f, w f, q im,ma, D a, and a ma for DPMIM. The following physicaltransport parameters were calculated from the fitted ones: the dispersivity k = D / v (L) for
CDE; h m = h  h im and the pore water velocity in the mobile phase v m = vh / h m (LT 1) for
MIM; the contribution to the overall water content from the fracture and matrix domains

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

77

h f * = w f h f and h ma * = (1  w f )h ma , respectively, a s calculated from Eq. (12),q f = ( q  (1  w f)q ma) / w f, pore water velocities in the fracture and matrix domainv f = q f / h f and v ma = q ma / h ma, respectively, and the relative contribution of the fracture
and matrix domain to the total flux Q f = w fq f / q and Q ma = (1  w f)q ma / q, respectively, for
DP and DP-MIM; and additionally, h m,ma and the pore water velocity of the mobile
domain of the matrix v m,ma = q ma / h m,ma (LT 1) for DP-MIM. Note that for DP, thedispersivity in the matrix domain, k ma, was fixed. Fitting this parameter resulted indramatically high confidence intervals for k ma and q ma, suggesting that these parameters
were not very sensitive due to the fact that the contribution of the matrix to the transport
was rather small. Thus, we selected to fix k ma at a value of 1 cm, which is the commonly
observed value for homogeneous soil columns (e.g., van Genuchten and Wierenga, 1977).
3.2.1.1. Physical equilibrium transport model — CDE. The CDE model did not describe
satisfactorily the Br data for either column. Dispersivity values were significantly larger
than the column length ( q max cases) or the retardation factors were larger than 1. In the
CDE concept, solutes are transported through the entire pore space and CDE can notdescribe two arrival fronts as observed for both columns at q max. Thus, large dispersivities
accounted for physical non-equilibrium processes not considered in the CDE model, whileretardation factors larger than 1 suggested that the CDE model used the sorption process todescribe a potential storage of Br in immobile or low flow regions.
Only at q mid, a good fit was obtained for column I but with a retardation factor
significantly smaller than 1 (R = 0.76). Such low R value could not be explained by anion
exclusion alone. Considering a mean clay surface area of 130 m2 g 1 (Bower and Hatcher,1966) and an average double layer thickness of 5 Å (El-Swaify et al., 1967; Veeh et al.,1994), an exclusion volume of 0.04 cm3 g 1 of clay was calculated. The total exclusion
volume for the grassed soil column amounted thus only to about 3% of the total porevolume of the column. The low R value suggested an acceleration of the solute transport
due to physical non-equilibrium processes.
From these results, it was concluded that the CDE transport model was not successful
in modeling the bromide data for both columns and at all rainfall intensities.
3.2.1.2. Physical non-equilibrium transport model — MIM. The MIM model could notdescribe two arrival fronts in both columns at q max (Fig. 1a and d). This result was
expected since the MIM model can predict only one front since it has only one mobilezone.
For column I, at q mid an excellent correspondence between the MIM model and data
was obtained (Fig. 1b) with an immobile water content of 0.122 cm3 cm3. The masstransfer rate coefficient between mobile and immobile zones, a ph, corresponded to acharacteristic transfer time t ph = 1 / a ph ~909 h, which is much longer than the bromide
residence time, t r = LR / v ~79 h. This high value of t ph suggested that not enough time wasleft for Br to diffuse into the microporosity during the experiment and that the mass transfer
rate could be neglected at this time scale (a ph ~ 0), resulting thus in the symmetrical shapeof the experimental BTC (van Genuchten and Wierenga, 1976). At q min a good fit was
found with an immobile water content of 0.245 cm3 cm 3 and t ph = 81 h compared tot r = 160 h (Fig. 1c).

78
Table 2Estimated physical transport parameters using four different physical transport models: the CDE, MIM, DP, and DP-MIM (for the highest velocity experiment only)transport models

Column I

q max
q mid
q min

Column II

q max
q mid

k
cm

R
–

r2
–

0.308
(0.0162)0.147(0.0589)0.0704(0.01781)
0.326(0.0072)0.1605(0.0027)119.01(27.50)8.17(0.47)
14.69(0.83)
69.94(5.95)
28.00(2.10)

1.00(349.80)0.76(0.01)
1.10(0.26)
1.00(380.54)1.20(0.22)

0.83710.98690.87040.87040.9697

MIM

Column I

q max
q mid
q min

Column II

q max
q mid

q
cm h 1

km
cm

h ma
cm3/cm 3

h im
cm3/cm 3

a ph
h 1

v ma
cm h 1

r2
–

0.308
(0.0162)0.147(0.0589)0.0704(0.01781)
0.326(0.0072)0.1605(0.0027)

4.13(1.33)
5.48(0.39)
4.94(0.61)
43.78(25.23)10.28(2.58)

0.088(0.030)0.264(0.033)0.132(0.016)0.238(0.087)0.174(0.035)0.310(0.009)0.122(0.005)0.245(0.008)0.136(0.079)0.201(0.026)0.0200(0.0006)
0.0011(0.0002)0.0123(0.0013)0.0445(0.0138)
0.0266(0.0098)

3.517(1.780)0.557(0.372)0.533(0.225)1.368(0.586)0.923(0.246)

0.9422
0.99370.98800.98140.9876

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

CDEqcm h 1

DP

Column I

q max
q mid

Column II

q max
q mid

k mab
cm

kf
cm

h +maacm3/cm 3

h +fa
cm3/cm 3

a sv mah 1v maa
cm h 1v fa
cm h 1

Q maa
–

Q fa
–

r2
–

0.308
(0.0162)0.147(0.0589)0.0704(0.01781)
0.326(0.0072)0.1605(0.0027)

1.00

5.74
(1.58)
5.40(0.36)
12.03(4.22)
20.85(6.98)
11.18(14.10)0.316(0.033)0.140(0.065)0.205(0.037)0.219(0.048)0.197(0.067)0.082(0.012)0.247(0.037)0.171(0.028)0.156(0.040)0.178(0.058)0.0443(0.0150)
0.0059(0.0015)
0.0212(0.0015)0.1829(0.0600)
0.1263(0.1500)

0.156(0.076)0.051(0.059)0.106(0.034)0.001(0.128)0.027(0.473)3.150(1.189)0.567(0.474)2.284(0.302)2.094(1.023)0.872(0.793)0.1597(0.1339)0.0482(0.0813)0.3078(0.7274)0.0004(0.0005)
0.0330(0.1538)0.8403(0.7045)0.9518(1.6049)
0.6922(1.6360)
0.9996(1.4104)
0.9670(4.5081)0.9408

1.00
1.00
1.00
1.00

0.99350.9881
0.97850.9878

DP-MIM
q
cm h

k ma
cm

kf
cm

h +maah +fa
a s
cm3 cm 3 cm3 cm 3 h 1

Column I q max 0.308
0.16
0.03 0.39(0.0162) (0.13) (2.41) (0.023)
Column II q max 0.326
8.30
7.84 0.347(0.0072) (12.64) (7.87) (0.084)

0.002(0.002)0.027(0.077)h m,maah im,maa ma
cm3 cm 3 cm3 cm 3 h 1

0.0030 0.063(0.0150) (0.025)
0.0180 0.223(0.1200) (0.199)

0.334(0.002)
0.124(0.115)v maa
v fa
cm h 1 cm h 1

v m,maacm h 1

Q maa
–

Q fa
–

r2
–

0.0196 0.715 16.043
4.519
0.92050.07950.9551(0.0006) (0.054) (397.499) (117.961) (47.3264) (4.0865)
0.0062 0.590 4.475
0.918
0.62830.37170.9878(0.0077) (0.719) (67.535) (17.226) (22.4665) (13.2911)

Standard errors are given in parenthesis, r 2 is the coefficient of determination for the regression of observed versus fitted concentrations, v m, v ma, v f, and v m,ma are averagepore velocities in the mobile region of the MIM model, in matrix and fracture domain of the DP model, and in the mobile region of the matrix domain of the DP-MIM
model, respectively, and Q ma and Q f are the relative contributions of the matrix and fracture region to the total flux, respectively.
aParameters calculated from the fitted parameters.
b
Fixed parameters.

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

q min

q
cm h 179

80

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

For column II, at q mid, the MIM model overestimated the experimental data, probably due
to the low mass recovery of Br (83.3%) (Fig. 1e). The immobile water content was 0.201
cm3 cm 3 and t ph = 38 h compared to t r = 70 h.
Different transport mechanisms were thus found for the two columns, as alreadydifferent shapes of the experimental BTCs suggested. For column I, a larger immobile
domain for q min compared to q mid suggested that decreasing the velocity and the total
water content created isolated zones in the soil that acted as immobile zones for the water
flux (van Genuchten and Wierenga, 1977; Gaudet et al., 1977). At the q mid case, a higher
value of a ph was found in column II compared to column I. This higher value is inagreement with the longer elution tail of the experimental BTC of column II compared tocolumn I. Complex transport processes that depended on water flux (differences within acolumn) and on spatial variability (differences between columns) thus occurred in our twocolumns. It was concluded that the MIM model correctly simulates bromide transport at
the intermediate and low rainfall intensities, although no explanation was found for the
much lower value of a ph at the intermediate rainfall intensity for column I. The MIM
model was not well suited to model rapid preferential flow with two arrival fronts
occurring at the largest rainfall intensities.
3.2.1.3. Physical non-equilibrium transport model — DP. For both columns, nosignificant differences were found between the optimized BTCs produced by the MIM
and DP models. For all cases, the pore water velocities in the fracture domain weresignificantly higher than in the matrix domain. The relative contribution of the matrix
domain to the total flux, Q ma, differed between columns and flow rates. It was 15.97%,
4.82% and 30.78% for q max, q mid and q min cases of column I and 0.04% and 3.33% for
q max and q mid cases of column II. Thus, non-negligible contributions of the matrix
domain were found for the q max and q min cases of column I. Note that large standard
errors of estimated flux v ma were obtained, compared to their mean values, when thecontribution of the matrix domain to flow was negligible. Except for the q max case ofcolumn II, the water content in the matrix domain was similar to the water content in theimmobile domain found with the MIM model. Higher values of the mass transfer rate
coefficient, a s, between the fracture and matrix domains were found compared to a ph of
the MIM model (Table 2).
In general the dual-permeability model should be able to predict a double front solute
migration when the transfer rate between the two domains is not too fast. However, in ourcase, the DP model, similarly to the MIM model, was not able to describe two arrival
fronts that occurred during the Br transport at the highest rainfall intensities. Results foundwith MIM and DP models were quite similar since instead of an immobile zone consideredin the MIM model, the DP model predicted very slow water flow in the matrix zone. This
slow flow, together with higher fitted transfer rates, resulted into a single solute front
predicted by the DP model. The effect of the slow flow zone would be more pronounced atsmaller water fluxes than those used in our study.
3.2.1.4. Physical non-equilibrium transport model — DP-MIM. The dual-permeability
mobile–immobile model was used to describe the transport of Br at the highest velocities.
This model can well describe a double front solute migration due to faster flow in the

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

81

fracture domain and slower flow in the mobile zone of the matrix domain. While both
MIM and DP models could relatively well describe the arrival of the main part of Br as
well as the long tailing, they could not predict its very early arrival (time between 2 (0.05)
and 6 (0.15) h (pore volume)). This could be achieved with the DP-MIM model, which
added additional flexibility by providing additional pore region, i.e., fracture domain,through which Br could travel quickly, by-passing the matrix region and arriving early atthe end of the column.
For column I, the DP-MIM model succeeded in describing the early solute front arrival,
although the tail of the BTC was slightly overestimated, similarly as with the othertransport models. For column II, the DP-MIM model was able to simulate well two solutepeaks. Differences in pore velocities and water content in the fracture domain (Table 2)
explained rather well observed differences between two experimental BTCs. For column I,
the fastest zone (fracture domain) consisted of a very small pore domain (h*f = 0.002 cm3
cm 3 only) combined with a large pore velocity (v f = 16.043 cm h 1). This resulted in theearly front arrival of low concentrations after about 2 h of infiltration (0.05 pore volume),which corresponds with the time needed to travel through the 30-cm long column. The
second large front was due to the matrix domain with a mobile water content, h m,ma, equal
to 0.063 cm3 cm 3 and a pore velocity, v m,ma, equal to 4.519 cm h 1. For column II, thefastest zone (fracture domain) consisted of a larger, although still small, pore domain(h f* = 0.027 cm3 cm 3) combined with a smaller pore velocity (v f = 4.475 cm h 1). This
combination of the water content and pore velocity led to the arrival of the first peak ofhigher concentrations after about 8 h (0.23 pore volume), which is close to the time needed
to travel through the 30-cm long soil column.
Different results were found for two mass transfer rate coefficients (a s and a ma) in thetwo columns. The mass transfer rate in column I was negligible between the fracture andmatrix domains, and dominant between mobile and immobile regions of the matrix
domain. Opposite results were found for column II (Table 2). These mass transfer rates
reflect the shape of the BTCs in the two columns that also appeared very different (Fig. 1aand d). The BTC for column I shows two solute fronts of very different concentrations.
The pore velocity in the fracture domain is very high and the mass transfer rate between
fracture and matrix domains is small. The long concentration tail is then explained by large
exchange between mobile and immobile regions of the matrix domain. The BTC forcolumn II, on the other side, shows two solute peaks of similar magnitude and relatively
smooth increases and decreases in concentrations, which requires faster exchange between
the fracture and matrix domains. The latter mass transfer rate can explain the slightretardation observed between the measured arrival time of the first peak (about 8 h (0.23
pore volume)) and the arrival time calculated from v f (6.7 h (0.19 pore volume)). Note thatpore velocities in the fracture domain of column II are of similar magnitude as porevelocities in the mobile region of the matrix domain of column I, and so are the masstransfer rate coefficients between these and slower, or stagnant, regions.
Relative contributions of the matrix domain to the total flow was 92.05% and 62.83%
for columns I and II, respectively. Although larger pores (fracture domain) are responsible
for an early arrival of the solute, they do not contribute a large part of the total flux. Largemacropores generally account for the main part of the flow rate at full saturation, but areinactive for unsaturated flow. Under experimental conditions in our study, macropores

82

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

were inactive for medium and low rainfall intensities, but became partly active for thehighest fluxes. Their effect would start dominating flow if full saturation was reached.
3.2.1.5. Intermediate conclusions about bromide transport modeling. Physical nonequilibrium models of different types (MIM, DP, DP-MIM) better described bromidetransport than physical equilibrium model (CDE). A preferential flow of a macropore typewas suggested based on the experimental data for the highest rainfall intensities (twoarrival fronts and marked shift to the left of peak positions). The DP-MIM model bettermodeled the experimental data than the other non-equilibrium transport models (MIM andDP) and supported the hypothesis of multiple (three) porosity domains with multiple (two)permeabilities at this rainfall intensity. However, an important note has to be madeconcerning the confidence one can have in these simulations. Non-negligible crosscorrelations between a large number of parameters and non-unique solutions were foundfor both columns. These problems were already discussed by many authors (e.g., Brusseau
et al., 1989; Hopmans et al., 2002). It must be kept in mind that the DP-MIM simulations
made here must be considered as a working hypothesis for testing the concept of multiple
porosities and permeabilities on our data rather than a definite answer. A preferential flowof a different type, close to a MIM type, was suggested for experiments at the othervelocities. No significant differences were found in the fits given by the MIM and DP
models, confirming therefore the hypothesis of a two-domain transport. However, in thecase of the lowest rainfall intensity, the DP model showed that the contribution of thematrix to the total flux can be significant. From these results, it was concluded that atintermediate and low rainfall intensities, the preferential transport observed is of a dual
porosity type with a zero or low flow in the matrix domain.
3.2.2. Chemical transport parameters (analyses of the herbicide BTCs)
Results of the model analyses of the herbicide BTCs are presented in Figs. 2 and 3 andTable 3. In the following analyses with chemical equilibrium and non-equilibrium models,the first-order degradation rate, l w, was estimated although no degradation products weredetected. As mentioned, the fact that no degradation products were detected in theeffluents did not mean that no decay processes occurred. The decay processes
hypothesized in Section 3.1 were (i) degradation with the formation of metabolites oflower mobility than the parent compound, (ii) mineralization or (iii) irreversible sorptionof either the parent compound or the metabolites. The models used in this study assumed a
total mass recovery of the chemicals so that all processes responsible for the loss of soluteand not included in the governing equations (degradation, mineralization, irreversible
sorption) are lumped into the first-order degradation rate l w.
3.2.2.1. Chemical equilibrium transport model — (DP-MIM and DP)-K d. The dualpermeability mobile–immobile model was used as the transport model for the highest
velocity experiments, while the dual-permeability model was considered for the othervelocities. The physical transport parameters (k ma, k f, h*ma , h*f , v ma, and v f for DP, andadditionally h m,ma, h im,ma, and v m,ma for DP-MIM) were fixed to values estimated for theBr transport. The mass transfer rate coefficients for the chemicals (a s for DP plus a ma for
DP-MIM) were calculated from the ones estimated for Br, corrected for the molecularV. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

83

Table 3
Estimated chemical transport parameters for isoproturon and metribuzin using DP transport model combined with
non-equilibrium chemical model (instantaneous and first-order kinetic retention) (DP-(K d, a ch) transport model)
IPU
q

Kd
1

Column I

q max
q mid
q min

Column II

q max
q mid

f ma
1

ff

a ch

lw

1

1

DT50a

r2

cm h

l kg

–

–

h

h

d

–

0.308
(0.0162)0.147(0.0589)0.0704(0.01781)
0.326(0.0072)0.1605(0.0027)1.31(0.16)
1.70(0.12)
2.24(0.39)
1.28(0.07)
1.50(0.16)

0.131(7.136)0.991(0.001)0.183(1.319)0.994(0.000)0.117(1.872)0.115(0.024)0.219(0.025)0.286(0.049)0.036(0.008)0.003(0.007)0.0547(0.0060)
0.0124(0.0015)
0.0067(0.0376)
0.0746(0.0040)0.0397(0.0019)0.0755(0.0036)0.0190(0.006)0.0178(0.0166)
0.0107(0.0019)0.0040(0.0011)0.38(0.02)
1.52(0.05)
1.62(1.51)
2.69(0.47)
7.26(1.99)

0.9739

q

Kd

f ma

ff

a ch

lw

DT50a

r2

cm h 1l kg 1

–

–

h 1

h 1

d

–

0.308
(0.0162)0.0704(0.01781)
0.326(0.0072)0.1605(0.0027)

1.17(0.19)
0.73(0.09)
0.64(0.10)
0.68(0.03)

0.137(4.480)0.161(1.427)0.321(2.863)0.000(0.848)0.082(0.022)0.289(0.051)0.004(0.014)0.015(0.013)0.0575(0.0038)
0.0199(0.0028)0.1296(0.0189)
0.0698(0.0026)

0.0255(0.0027)0.0040(0.005)0.0065(0.0020)0.0123(0.0004)

1.29(0.15)
7.16(0.80)
4.45(1.35)
2.36(0.07)

0.97430.93410.74540.98240.9600

Met

Column I

q max
q min

Column II

q max
q mid

0.97050.97990.9902Standard errors are given in parenthesis and r 2 is the coefficient of determination for the regression of observed
versus fitted concentrations.
aParameters calculated from the fitted parameters.

diffusion coefficient of pesticides, taken as 0.0179 cm2 h 1 (Jury et al., 1983). The
number of chemical parameters to be estimated was therefore limited to K d and l w for the
DP model, and additionally f m,ma for the DP-MIM model.
Results showed that for all rainfall intensities and regardless of the used physicaltransport model the chemical equilibrium transport model was not well suited to describe
the transport of herbicides in the grassed soil columns. Systematically, the model strongly
underestimated the observed elution tails and could not correctly predict the arrival time of
the first herbicide non-zero concentrations and peak positions (Figs. 2 and 3).
3.2.2.2. Chemical non-equilibrium transport model — DP-(K d, a ch). Chemical nonequilibrium is not currently available in the dual-permeability mobile–immobile transportmodel of HYDRUS-1D. Therefore, we used the dual-permeability model as the physicaltransport model also for the highest velocity cases in order to estimate the importance of
non-equilibrium chemical transport for the two herbicides. The dual-permeability model
was thus viewed as an approximation of the physical transport occurring at the highest

84

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

rainfall intensities. Similarly as for the chemical equilibrium model, the physical transportparameters (k ma, k f, h ma*, h*f , v ma, v f ) were fixed to the values estimated for Br transport
and the chemicals mass transfer rate coefficients (a s) were calculated from the ones
obtained for Br, corrected for the molecular diffusion coefficient of pesticides. The
chemical parameters to estimate were therefore K d, a ch, f ma, f f and l w.
For both herbicides, the non-equilibrium chemical model described well theexperimental BTCs. Particularly, the elution tails of BTCs were well captured. Only theearly arrival of isoproturon for the q min case of column I was not well described (Fig. 2c).
The early arrival front at q max for column I was not described (Figs. 2a and 3a).
A comparison of the sorption parameters (K d and a ch) between the two herbicides
showed higher K d values for IPU and comparable to higher a ch values for Met (Table 3). A
higher sorption capacity of the soil for IPU compared to Met is in agreement with theliterature data that show that Met is more mobile than IPU (Agritox, 2001). The absolute
K oc values, calculated as K oc = K d / CO with CO being the organic carbon content (equal to17.29 g kg 1 from Madrigal et al., 2002) were in the range of 74 to 130 l kg 1 for IPU and
37 to 68 l kg 1 for Met. These values were thus also in the range of the literature data
(Agritox, 2001). For IPU, the K d values were lower than those derived from batch sorption
experiments with the same soil data (Madrigal et al., 2002). Note that good fits could not
be found when the K d parameter was fixed to the latter batch value of 2.4 l kg 1.
For IPU, K d values increased when rainfall intensities decreased, with comparable
values between the two columns. For Met, no significant trend was observed for K d values.
For both herbicides, a significant decrease in the sorption kinetic rate coefficient, a ch, wasobserved as velocity decreased. Maraqa (2001), in a review of column displacement
experiments, reported similar trends. A possible explanation is that for lower velocities
more sites seem to be at relative equilibrium with the solution ( f f usually increased whenvelocity decreased), and given the lower sorption kinetic rate the overall rate of reaction isthen more or less the same. Note that the values of f ma were all accompanied by very large
errors or were found to be close to 1. A series of additional simulations was performed
with this parameter fixed to different values (0.1, 0.5 and 0.9). No significant differenceswere found in simulated BTCs, showing that f ma is not a sensitive parameter for the
herbicide transport. Since the flow in the matrix region is slower, there is more time toreach equilibrium in the matrix domain and thus the fraction of exchange sites assumed to
be in equilibrium with the solution phase of the matrix domain is not a relevant parameter.
Equilibrium sorption in the matrix domain is thus a suitable model. By contrast, for thefracture domain where there is not enough time to reach equilibrium, due to faster flow,
non-equilibrium sorption is an appropriate model to explain the observed long elution tails.
Sensitivity analyses performed with f f fixed to different values (0.1, 0.5 and 0.9) showedthat the elution tails can not be described with large f f values (0.5 and 0.9).
Using soil samples originating from the same grassed strip, Benoit et al. (1999) haveshown an increase in IPU sorption during soil incubation at 18 8C and the water content
close to the soil water holding capacity. This was done by considering extractions withcalcium chloride and methanol at different dates and soil depths (0–2; 2–13 and 13–30 cm)
where the herbicide concentration was measured by HPLC. From these data, we estimateda mean sorption kinetic rate parameter for IPU using a simple first order equation to fit theincrease in sorbed IPU versus time and considering the 0–30 cm horizon as one single

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

85

layer. We calculated a value of 0.0122 h 1 which is of the same order of magnitude as the
IPU a ch values estimated by the DP-(K d, a ch) model (Table 3).
A longer half-life for IPU compared to Met for q mid case of column II and lower DT50s
for all other cases are in agreement with the IPU mass recovery rates in the effluent that are
higher for q mid case of column II and smaller for all other cases compared to Met (Table 1).
However, the DT50s values were dramatically lower than the lower bounds of theliterature ranges for both pesticides. Fitted degradation rates were up to 31 and 18 times
larger for IPU and Met, respectively, than the literature values (Agritox, 2001). Althoughattempts have been made to fit IPU BTCs with l w fixed to the measured batch value of0.001203 h 1 in the grassed strip soil (Benoit et al., 1999), no good fits could be found.
Several other authors have also found markedly faster degradation rates in column
displacement studies compared to batch data. Up to 3 times larger degradation rates for
saturated experiments and up to 14 times for unsaturated experiments were found by Rocio
Estrella et al. (1993); and up to 7 times larger rates in saturated experiments were found by
Guo and Wagenet (1999). Other authors found higher or comparable degradation rates(Gaston et al., 2003; Das et al., 2004) and/or smaller degradation rates (Langner et al.,
1998; Simoni et al., 2001). Higher degradation rates in column displacement experiments
compared to batch data were usually attributed to differences in aeration and mixing in thebatch and column systems, to decrease of substrate concentrations in batch experiments
compared with constant influent substrate concentrations in column experiments, or to ahigh degree of non-equilibrium sorption (Guo and Wagenet, 1999). Smaller degradation
rates were attributed to rate-limited diffusion of the solutes to the microbial cells (Simoni
et al., 2001).
The first-order degradation rate, l w, was found to increase when velocity increased,
except for Met of column II. Angley et al. (1992) and Guo and Wagenet (1999) foundsimilar results for alkylbenzenes and alachlor, respectively. Our modeling analysis
suggested that the degree of non-equilibrium sorption increased as velocity increased ( f fusually decreased). Thus, our results suggest that the estimated degradation rates are
positively related to the degree of non-equilibrium sorption. It supports the hypothesis thatchemical irreversible sorption, defined as fast attachment and much slower detachment,
occurred most probably on the sorption sites where kinetic sorption occurs. The
irreversible sorption acts as a disappearance process at the time scale of our experiments
and is therefore accounted for by apparently higher degradation rates by the models. This
could explain the apparent contradiction of higher degradation rates when higher mass
recovery rates in the effluent are measured.

4. Conclusions
Two series of Br and herbicides (IPU and Met) displacement experiments were
performed at different rainfall intensities (0.070, 0.147, 0.161, 0.308 and 0.326 cm h 1) in
two soil columns sampled in a 10-year-old grassed strip at La Jaillière. The experimental
Br BTCs, and their descriptive and numerical analyses, suggested two contrasted physicalnon-equilibrium processes depending on the rainfall intensity. For the highest rainfall
intensities (0.308 and 0.326 cm h 1), multiple porosity domains and multiple

86

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

permeabilities preferential flow were probably active: a) rapid flow through macropore
pathways, b) slower flow through a mesoporosity, and c) no-flow in the remaining
micropores. The total mobile porosity was suggested to be less than 0.07 cm3 cm 3 interms of the water content. For the lower rainfall intensities, macropore flow was not
active anymore.
The experimental BTCs of both herbicides confirmed the hypothesis of a rapid
preferential flow at the highest rainfall intensities. The main conclusion drawn from themodeling analyses is that non-equilibrium chemical transport processes occurred for all
velocities and both herbicides and that the degree of non-equilibrium sorption increased
with velocity. IPU was found less mobile than Met in accordance with batch data.
Significantly higher values of degradation rate parameters as compared to batch data werepositively related with the degree of non-equilibrium sorption. This work suggests that thesorption and/or degradation rates obtained from batch experiments can not be used in nonequilibrium transport models to describe the transport experiments unless the actual
underlying mechanisms measured by the two techniques are fully identified.
From an environmental point of view, these series of displacement experiments allowed
to better understand the functioning of the grassed strip at La Jaillière. Particularly, theseresults emphasized the potential role of macropores in establishing rapid preferential
transport of herbicides and in by-passing the surface layer of the grassed soil at a relatively
low rainfall intensity (0.3 cm h 1). This work should be completed with displacement
experiments under ponding conditions, as the grassed strip soils are mainly used asprotective management practices for the reduction of pesticides from agricultural runoff
waters.
Acknowledgments
We would like to thank Christophe Labat for his help in sampling at the experimental
site. This study was supported financially by the INRA-Cemagref bAQUAEQ program and
the Pesticide Research program of the French Ministry of the Ecology and Sustainable
Development. Support from OECD (a fellowship under the OECD Co-operative Research
Program: Biological Resource Management for Sustainable Agriculture Systems) and
INRA for J. Šimlnek’s summer stay at the UMR INRA/INAPG EGC-Grignon is greatly
acknowledged.
References
Agritox: Toxixological Data Base on Pesticides, 2001. INRA, France.
Angley, J.T., Brusseau, M.L., Lamar Miller, W., Delfino, J.J., 1992. Nonequilibrium sorption and aerobic
biodegradation of dissolved alkylbenzenes during transport in aquifer material: column experiments andevaluation of a coupled-process model. Environ. Sci. Technol. 26, 1404 – 1410.
Arora, K., Mickelson, S.K., Baker, J.L., Tierney, J.L., Peters, C.J., 1996. Herbicide retention by vegetative bufferfrom runoff under natural rainfall. Trans. ASAE 39, 2155 – 2162.
Baker, J.L., Mickelson, S.K., Hatfield, J.L., Fawcett, R.S., Hoffman, D.W., Franti, T.G., Tierney, D.P., 1995.
Reducing herbicide runoff: role of best management practices. In British Crop Port. Conf. — Weeds. BCPC
Publications, Farnham, pp. 479 – 487.

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

87

Benoit, P., Barriuso, E., Vidon, P., Réal, B., 1999. Isoproturon sorption and degradation in a soil from grassed
buffer strip. J. Environ. Qual. 28, 121 – 129.
Benoit, P., Barriuso, E., Vidon, P., Réal, B., 2000. Isoproturon movement and dissipation in undisturbed soil coresfrom a grassed buffer strip. Agronomie 20, 297 – 307.
Bower, C.A., Hatcher, J.T., 1966. Simultaneous determination of surface area and cation-exchange capacity. Soil
Sci. Soc. Am. Proc. 30, 525 – 527.
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 – 188.
Das, B.S., Lee, L.S., Rao, P.S.C., Hultgren, R.P., 2004. Sorption and degradation of steroid hormones in soils
during transport: column studies and model evaluation. Environ. Sci. Technol. 38, 1460 – 1470.
Delphin, J.-E., Chapot, J.-Y., 2001. Leaching of atrazine and deethylatrazine under a vegetative filter strip.
Agronomie 21, 461 – 470.
Dillaha, T.A., Reneau, R.B., Mostaghimi, S., Lee, D., 1989. Vegetative filter strips for agricultural nonpoint
source pollution control. Trans. ASAE 32, 513 – 519.
El-Swaify, S., Coleman, N.T., Bredell, G., Arca, M., 1967. Negative adsorption by vermiculite: salt exclusionfrom interlayer volumes. Soil Sci. Soc. Am. Proc. 31, 464 – 466.
Gaston, L.A., Boquet, D.J., Bosch, M.A., 2003. Fluometuron sorption and degradation in cores of silt loam fromdifferent tillage and cover crop systems. Soil Sci. Soc. Am. J. 67, 747 – 755.
Gaudet, J.-P., Jégat, H., Vachaud, G., Wierenga, P.J., 1977. Solute transfer, with exchange between mobile andstagnant water, through unsaturated sand. Soil Sci. Soc. Am. J. 41, 665 – 671.
Gerke, H.H., van Genuchten, M.Th., 1993. A dual-porosity model for simulating the preferential movement ofwater and solutes in structured porous media. Water Resour. Res. 29, 305 – 319.
Gerke, H.H, van Genuchten, M.Th., 1996. Macroscopic representation of structural geometry for simulating waterand solute movement in dual-porosity porosity media. Adv. Water Resour. 19, 343 – 357.
Guo, L., Wagenet, R.J., 1999. Evaluation of alachlor degradation under transport conditions. Soil Sci. Soc. Am. J.
63, 443 – 449.
Haynes, R.J., Francis, G.S., 1993. Changes in microbial biomass C, soil carbohydrate composition andaggregate stability induced by growth of selected crop and forage species under field conditions. J. Soil Sci.
44, 665 – 675.
Hopmans, J.W., Šimlnek, J., Romano, N., Durner, W., 2002. Inverse modeling of transient water flow. In: Dane,
J.H., Topp, G.C. (Eds.), Methods of Soil Analysis: Part 1. Physical Methods, Chapter 3.6.2, Third edition.
SSSA, Madison, WI, pp. 963 – 1008.
Hurto, K.A., Turgeon, A.J., Cole, M.A., 1979. Degradation of benefin and DCPA in thatch and soil from
Kentucky Bluegrass (Poa pratensis) turf. Weed Sci. 27, 154 – 157.
Jury, W.A., Spencer, W.F., Farmer, W.J., 1983. Behavior assessment model for trace organics in soil: I. Modeldescription. J. Environ. Qual. 12, 558 – 564.
Kanchanasut, P., Scotter, D.R., 1982. Leaching patterns in soil under pasture and crop. Aust. J. Soil Res. 20,
193 – 202.
Kay, B.D., 1990. Rates of changes of soil structure under differents cropping systems. In: Steward, B.A. (Ed.),
Advances in Soil Science. Springer-Verlag, pp. 1 – 51.
Langner, H.W., Inskeep, W.P., Gaber, H.M., Jones, W.L., Das, B.S., Wraith, J.M., 1998. Pore watervelocity and residence time effects on the degradation of 2,4-D during transport. Environ. Sci. Technol.
32, 1308 – 1315.
Lapidus, L., Amundson, N.R., 1952. Mathematics of adsorption in beds: VI. The effect of longitudinal diffusion
in ion exchange and chromatographic columns. J. Phys. Chem. 56, 84 – 988.
Madrigal, I., Benoit, P., Barriuso, E., Etiévant, V., Souiller, C., Réal, B., Dutertre, A., 2002. Capacités de stockage
et d’épuration des sols de dispositifs enherbés via-à-vis des produits phytosanitaires. Deuxième partie:
propriétés de rétention de deux herbicides, l’isoproturon et le diflufénicanil dans différents sols de bandes
enherbées. Etude et Gestion des Sols 9, 287 – 302.
Maraqa, M.A., 2001. Prediction of mass-transfer coefficient for solute transport in porous media. J. Contam.
Hydrol. 50, 1 – 19.

88

V. Pot et al. / Journal of Contaminant Hydrology 81 (2005) 63–88

Merckx, R., Den Hartog, A., Van Veen, J.A., 1985. Turnover of root-derived material and related biomass
formation in soils of different texture. Soil Biol. Biochem. 17, 565 – 569.
Mersie, W., Seybold, C., Tsegaye, T., 1999. Movement, adsorption and mineralization of atrazine in two soils
with and without switchgrass (Panicum virgatum) roots. Eur. J. Soil Sci. 50, 343 – 349.
Patty L., 1997. Limitation du transfert par ruissellement vers les eaux superficielles de deux herbicides(isoproturon et diflufénicanil). Méthodologie analytique et étude de l’efficacité de bandes enherbées. PhD
Thesis. Université Joseph Fourier, Grenoble, France.
Patty, L., Réal, B., Gril, J.-J., 1997. The use of grassed buffer strips to remove pesticides, nitrate and soluble
phosphorus compounds from runoff water. Pest. Sci. 49, 243 – 251.
Reungsang, A., Moorman, T.B., Kanwar, R.S., 2001. Transport and fate of atrazine in midwestern riparian buffer
strips. J. Am. Water Resour. Assoc. 37, 1681 – 1692.
Rocio Estrella, M., Brusseau, M.L., Maier, R.S., Pepper, I.L., Wierenga, P.J., Miller, R.M., 1993. Biodegradation,
sorption, and transport of 2,4-Dichlorophenoxyacetic acid in saturated and unsaturated soils. Appl. Environ.
Microbiol. 59, 4266 – 4273.
Selim, H.M., Davidson, J.M., Rao, P.S.C., 1977. Transport of reactive solutes through multilayered soils. Soil.
Sci. Soc. Am. J. 41, 3 – 10.
Simoni, S.F., Schäfer, A., Harms, H., Zehnder, A.J.B., 2001. Factors affecting mass transfer limited
biodegradation in saturated porous media. J. Contam. Hydrol. 50, 99 – 120.
Šimlnek, J., Šejna, M., van Genuchten, M.Th., 1998. The HYDRUS-1D software package for simulating the
one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0,
IGWMC-TPS-70. International Ground Water Modeling Center, Colorado School of Mines, Golden,
Colorado. 202 pp.
Šimlnek, J., Wendroth, O., Wypler, N., van Genuchten, M.Th., 2001. Nonequilibrium water flow characterizedfrom an upward infiltration experiment. Eur. J. Soil Sci. 52 (1), 13 – 24.
Šimlnek, J., Jarvis, N.J., van Genuchten, M.Th., Gärdenäs, A., 2003. Review and comparison of models fordescribing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272, 14 – 35.
Souiller, C., Coquet, Y., Pot, V., Benoit, P., Réal, B., Margoum, C., Laillet, B., Labat, C., Vachier, P., Dutertre, A.,
2002. Capacités de stockage et d’épuration des sols de dispositifs enherbés via-à-vis des produits
phytosanitaires. Première partie: dissipation des produits phytosanitaires à travers un dispositif enherbé;
mise en évidence des processus mis en jeu par simulation de ruissellement et infiltrométrie. Etude et Gestion
des Sols 9, 269 – 285.
van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer studies in sorbing porous media I. Analyticalsolutions. Soil Sci. Soc. Am. J. 40, 473 – 480.
van Genuchten, M.Th., Wierenga, P.J., 1977. Mass transfer studies in sorbing porous media II. Experimental
evaluation with tritium (3H2O). Soil Sci. Soc. Am. J. 41, 272 – 278.
van Genuchten, M.Th., 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.
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.
Vogeler, I., Scotter, D.R., Green, S.R., Clothier, B.E., 1997. Solute movement through undisturbed soil columns
under pasture during unsaturated flow. Aust. J. Soil Res. 35, 1153 – 1163.
Weast, R.C., 1985. Handbook of Chemistry and Physics, 66th edition. CRC Press, Boca Raton, Florida.