This article appeared in a journal published by Elsevier. The attachedcopy is furnished to the author for internal non-commercial researchand education use, including for instruction at the authors institution
and sharing with colleagues.
Other uses, including reproduction and distribution, or selling or
licensing copies, or posting to personal, institutional or third party
websites are prohibited.
In most cases authors are permitted to post their version of thearticle (e.g. in Word or Tex form) to their personal website or
institutional repository. Authors requiring further information
regarding Elsevier’s archiving and manuscript policies areencouraged to visit:
http://www.elsevier.com/copyrightAuthor's personal copy
Journal of Contaminant Hydrology 104 (2009) 107–125

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology
j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j c o n h yd

Impact of varying soil structure on transport processes in different
diagnostic horizons of three soil types
Radka Kodešová a,⁎, Nadia Vignozzi b, Marcela Rohošková a, Tereza Hájková a, Martin Kočárek a,Marcello Pagliai b, Josef Kozák a, Jirka Šimůnek caDept. of Soil Science and Soil Protection, Czech University of Life Sciences Prague, Kamýcká 129, 16521 Prague 6, Czech Republic
CRA-ABP, Centro di ricerca per l'agrobiologia e la pedologia, Piazza M Azeglio 30, I-50121 Florence, Italyc
Dept. of Environmental Sciences, University of California Riverside, 900 University Avenue, Riverside, CA 92521, USAba r t i c l e

i n f oArticle history:
Received 12 April 2008
Received in revised form 13 October 2008
Accepted 22 October 2008
Available online 31 October 2008
Keywords:
Soil pore structure
Soil structure stability
Column experiment
Preferential ﬂow
Pesticide transport
Single-porosity model
Dual-permeability modelTwo-site sorption modela b s t r a c t
When soil structure varies in different soil types and the horizons of these soil types, it has asigniﬁcant impact on water ﬂow and contaminant transport in soils. This paper focuses on theeffect of soil structure variations on the transport of pesticides in the soil above the water table.
Transport of a pesticide (chlorotoluron) initially applied on soil columns taken from various
horizons of three different soil types (Haplic Luvisol, Greyic Phaeozem and Haplic Cambisol)
was studied using two scenarios of ponding inﬁltration. The highest inﬁltration rate andpesticide mobility were observed for the Bt1 horizon of Haplic Luvisol that exhibited a welldeveloped prismatic structure. The lowest inﬁltration rate was measured for the Bw horizon of
Haplic Cambisol, which had a poorly developed soil structure and a low fraction of large
capillary pores and gravitational pores. Water inﬁltration rates were reduced during theexperiments by a soil structure breakdown, swelling of clay and/or air entrapped in soil
samples. The largest soil structure breakdown and inﬁltration decrease was observed for the Ap
horizon of Haplic Luvisol due to the low aggregate stability of the initially well-aggregated soil.
Single-porosity and dual-permeability (with matrix and macropore domains) ﬂow models in
HYDRUS-1D were used to estimate soil hydraulic parameters via numerical inversion using
data from the ﬁrst inﬁltration experiment. A fraction of the macropore domain in the dualpermeability model was estimated using the micro-morphological images. Final soil hydraulicparameters determined using the single-porosity and dual-permeability models were
subsequently used to optimize solute transport parameters. To improve numerical inversion
results, the two-site sorption model was also applied. Although structural changes observed
during the experiment affected water ﬂow and solute transport, the dual-permeability model
together with the two-site sorption model proved to be able to approximate experimental data.
© 2008 Elsevier B.V. All rights reserved.

1. Introduction
Soil and groundwater contamination caused by pesticides
used in agriculture is an environmental problem worldwide. The
controlling processes, which include water ﬂow and contaminanttransport, strongly depend on the composition of soil structurecomponents resulting from pedogenesis and soil management.

⁎ Corresponding author.
E-mail address: kodesova@af.czu.cz (R. Kodešová).
0169-7722/$ – see front matter © 2008 Elsevier B.V. All rights reserved.
doi:10.1016/j.jconhyd.2008.10.008

Different compositions of soil structure components are reﬂectedin soil properties that affect water and solute transport, such as
the soil structure and the conﬁguration of the soil porous system,
the aggregate stability, and soil hydraulic properties. Soilaggregation is under the control of different mechanisms in
different soil types. Flocculated clay particles, or their complexes
with humus (organo-mineral complexes) and soil organic matter,
act as the main cementing agents in soil aggregates. The
cementing effect of free Fe and Al oxides is important in soils
with low organic matter content (Six et al., 2002). Generally, thelevel of aggregation and the stability of aggregates increase with

Author's personal copy
108

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

increasing organic matter content, surface area of clay minerals,
and cation exchange capacity (Bronick and Lal, 2005). In addition
to soil properties, soil management also has a very important
inﬂuence on the stability of soil structure. By tillage the topsoil ismixed and new aggregates are exposed to different breakdown
mechanisms (Six et al., 1998). Soil processing at improper soil
moisture, crossing of heavy machinery, irrigation, and use offertilizers can also lead to soil structure degradation (Pagliai et al.,2003, 2004 and Servadio et al., 2005). Soil structure, soil porous
systems and subsequently soil hydraulic properties are also
inﬂuenced by many additional factors, including the mineralogical composition of the parent material, its alteration degree,organic matter content, soil water regime, transport processes
within the soil proﬁle, climate, plant roots, soil organisms, andmanagement practices.
The shapes and sizes of soil pores can be studied by image
analysis of thin soil sections at various magniﬁcations. Porous
systems with macropores and their impact on saturated
hydraulic conductivities, Ks, were explored by Bouma et al.
(1977, 1979) using micromorphometric data. Differently
shaped pores in soils under different management practices
and their Ks values were studied by Pagliai et al. (1983, 2003,
2004) using both micromorphometric and micromorphological investigations. Finally, Kodešová et al. (2006, 2007,2008) described the effects of gravitational pores and large
capillary pores detectable in micromorphometric images on
the shape of soil hydraulics functions.
Multi-modality of the soil porous system and hierarchicalpore composition may cause non-equilibrium water ﬂow andsolute transport in soils. Numerical models that assume multimodal soil porous systems have been developed to describe nonequilibrium water ﬂow and solute transport in such soils.
Overviews of various approaches were given by Šimůnek et al.
(2003), Gerke (2006), Köhne et al. (2009a,b-this issue), orŠimůnek and van Genuchten (2008b). In these models the soil
porous system is divided into two domains, and each domain ischaracterized by its own set of transport properties and equations
describing ﬂow and transport processes. The dual-porosity
approach deﬁnes water ﬂow and solute transport in systems
consisting of domains of both mobile and immobile water. The
dual-porosity formulation is based on a set of equations
describing water ﬂow and solute transport in the mobile domain,and mass balance equations describing soil water and solutecontent in the immobile domain. On the other hand, the dualpermeability approach assumes that water ﬂow and solutetransport occur in both domains. The dual-permeability formulation is based on a set of equations that describe water ﬂow andsolute transport separately in each domain (matrix and macropore domains). Numerical models based on the multi-modality ofthe soil porous system are very useful for the assessment of soil
and groundwater contamination from point and diffused sources.
An alternative approach to dual-porosity and dual-permeability models can be to use composite bimodal soil hydraulicproperty models in the single-porosity ﬂow and transportmodels (e.g., Othmer et al., 1991; Durner, 1994; Mohanty et al.,
1997; Coppola et al., 2009-this issue) that can account for asigniﬁcant increase in the hydraulic conductivity near saturation. However, these models do not by themselves lead topreferential ﬂow, since when used in a single-domain model
based on the Richards equation, simulations will still produce
uniform wetting fronts, although possibly with some acceler-ated advance of the front for surface ponding conditions
because of the higher saturated conductivity. Such models,
however, cannot predict preferential ﬂow since there is no
mechanism to account for lateral nonequilibrium in terms ofwater moving primarily through larger pores or fractures andbypassing the matrix pore-space (Šimůnek et al., 2003).
Water ﬂow and pesticide transport in soils has been
recently studied, both experimentally and numerically, by
Köhne et al. (2005, 2006a,b), Pot et al. (2005), Gärdenäs et al.
(2006), Kodešová et al. (2005, 2008), and Dousset et al.
(2007). These studies showed that pesticide transport isfrequently affected by preferential ﬂow due to the presence ofmacropores. All these studies also considered pesticide
application on the top of humic horizons, where pesticide
adsorption and degradation are higher than in deeper
horizons due to the presence of organic matter and microorganisms. In these studies, it was often observed that
pesticides pass through the top humic horizon and contaminate horizons below relatively quickly.
The objective of the current study was to evaluate ﬂow andtransport in both the top and underlying soil horizons. For
this purpose, undisturbed soil samples were taken fromvarious diagnostic horizons of three different soil types andsubjected to various analyses. Conﬁgurations of the soil
porous system, soil aggregate stability, water and pesticidetransport were investigated on all samples to assess either
protective or non-protective functions of various soil horizons
with respect to groundwater contamination.
2. Materials and methods
2.1. Soil porous system
The study was performed on Haplic Luvisol (parent material
loess), Greyic Phaeozem (parent material loess) and Haplic
Cambisol (parent material orthogneiss). The elevations of soil
locations above sea level are 272 m for Haplic Luvisol, 276 m for
Greyic Phaeozem and 523 m for Haplic Cambisol. Averageannual temperatures and average annual precipitations (over30 years) are 8.5 °C and 618 mm for Haplic Luvisol, 8.6 °C and
578 mm for Greyic Phaeozem, and 7.4 °C and 660 mm for Haplic
Cambisol, respectively. A ﬁve year rotation system withconventional tillage was used at all locations. Winter barley
was planted at all areas when ﬁeld work was performed on
August 9 and 10, 2006. Four soil horizons (Ap 0–35 cm, Bt1 35–
62 cm, Bt2 62–102 cm, Ck 102–133 cm) were identiﬁed in Haplic
Luvisol, three horizons (Ap 0–21 cm, Bth 21–27 cm, BCk 27–
37 cm, Ck 37–80 cm) in Greyic Phaeozem, and three horizons(Ap 0–32 cm, Bw 32–62 cm, Ck 62–97 cm) in Haplic Cambisol.
Large undisturbed soil samples of 10 ×8 ×4 cm were taken fromeach horizon of each soil proﬁle except for the Bth horizon of
Greyic Phaeozem (the thickness of this horizon was only 5 cm).
Micromorphological properties characterizing the soil pore
structure were studied on thin soil sections prepared from thesoil samples. The ﬁnal thin section size was approximately
7 × 4 cm. Two thin sections were made from each soil sample.
The soil porous system was analyzed using the proceduredescribed by Pagliai et al. (1983, 2003, 2004). Two images ofeach thin section were taken at one magniﬁcation at a resolutionof 300 dpi. The size of the images was 2048 × 1536 pixels, andthe pixel size was 15.7 µm. The applied method allowedAuthor's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

studying pores larger than approximately 40 µm. The shapefactor, SF [−], was calculated according to Pagliai et al. (1983)
using the following equation:
SF =

P2
4πA

ð1Þ

where P [L] is the pore perimeter and A [L2] is the pore area.
Pores were divided into different shape groups: regular with SFof 1–2, irregular with SF of 2–5, and elongated with SF larger
than 5. These groupings correspond approximately to thoseused by Bouma et al. (1977). Pores of each shape group were
further subdivided into size classes either according to theirequivalent pore diameter (regular and irregular pores) or theirwidth (elongated pores) (Pagliai et al., 1983).
2.2. Soil structure stability

Wds
Wds + Wdw

day after the second inﬁltration, all soil samples were cut into 8
layers. The sub-samples were dried and the mass of soil in eachlayer was measured. The soil was then ground and sieved
through the 2-mm sieve. The total amount of chlorotoluron
remaining in each soil layer was ascertained as follows. 5 g of thedry soil was placed into the centrifuge curette. 5 ml of methanol
was added and the centrifuge cuvette was placed for 15 h into
the shaking apparatus. After that, analyzed soil samples were
centrifuged for 30 min at 13,800 rpm. The chlorotoluron
concentration in the methanol extract was again measured
using HPLC. The chlorotoluron concentration in the soil sample
was expressed as the total amount of solute per a unit of mass ofdry soil (µg g− 1). Finally, chlorotoluron recovery was found bycomparing the amount of applied solute with the sum of soluteleached at the bottom of the column and remaining in the soil.
2.4. Numerical evaluation of experimental data

The stability of the soil structure was evaluated using theindex of water-stable aggregates (Nimmo and Perkins, 2002).
Four grams of air dry soil aggregates (segregates) of the size of
2–5 mm was sieved for 3 min in distilled water (sieve 0.25 mm).
Aggregates remaining on the sieve were next sieved in sodium
hexametaphoshate until only sand particles remained on thesieve. The index of water-stable aggregates, WSA [−], was thendetermined as:
WSA =

109

ð2Þ

where Wds [M] is the weight of aggregates dispersed indispersing solution and Wdw [M] is the weight of aggregatesdispersed in distilled water. In addition, organic carbon
contents, clay contents and cation exchange capacity, which
all affect soil structure stability, were measured.
2.3. Water ﬂow and solute transport experiments
Plastic cylinders with a diameter of 10.5 cm and a height of
13 cm were used to collect undisturbed soil samples. Undisturbed soil samples were taken from the Ap (0–13 cm), Bt1 (42–
55 cm), and Bt2 (65–78 cm) horizons of Haplic Luvisol, the Ap (0–13 cm) and Ck (46–59 cm) horizons of Greyic Phaeozem, and the
Ap (0–13 cm) and Bw (40–53 cm) horizons of Haplic Cambisol.
Soil samples were transported into the laboratory where microtensiometers Tensior 5 (UMS GmbH, Munich, 2005) for pressure
head measurements were installed 4 and 8 cm below the soil
surface. Cementation was used to prevent water ﬂow and solutetransport between the soil core and plastic cylinder. 10 cm3 ofchlorotoluron solution (198 µg cm− 3 of chlorotoluron) wasapplied with an inﬁltration rate of 0.116 cm min− 1 on the top of
the soil samples. Two ponded inﬁltration experiments with
distilled water were then performed during the following 2 days.
The ponding depth was maintained at 1.5 cm. The duration of theinﬁltration experiment varied between 30 min and 7 h depending on the inﬁltration rate. The duration of individual experiments was set to study possible chlorotoluron leakage from soil
horizons under extreme inﬁltration conditions, and to evaluate
chlorotoluron retention in the soil. Water outﬂow and soluteconcentration from the bottom of the soil sample weremonitored in time. Solute concentrations were determinedusing High Performance Liquid Chromatography (HPLC). One

2.4.1. Single-porosity and dual-permeability water ﬂow models
Water ﬂow in the soil proﬁle may be simulated using thesingle-porosity and dual-permeability models implemented in
HYDRUS-1D (Šimůnek et al., 2008a,b). The Richards equation,
describing the one-dimensional isothermal Darcian ﬂow in avariably saturated rigid porous medium, is used in both models.
The Richards equation is used to describe water ﬂow in the
single-porosity system:


Aθ A
Ah
=
K ðhÞ
+ K ðhÞ
At Az
Az

ð3Þ

where θ is the volumetric soil water content [L3L− 3], h is thepressure head [L], K is the hydraulic conductivity [LT− 1], t is
time [T], and z is the vertical axes [L]. Eq. (3) is solved for theentire ﬂow domain using one set of soil water retention andhydraulic conductivity functions.
In the case of the dual-permeability model, the Richards
equation is applied separately to each of the two pore regions,
i.e., the macropore (fractures, a domain of larger pores) andmatrix domains (Gerke and van Genuchten, 1993):
  Ahf
  Cw
Aθf
A
Kf hf
=
+ Kf hf −
Az 
At
Az
fw

Aθm
A
Ahm
Cw
Km ðhm Þ=+ Km ðhm Þ +
Az
At
Az
1−fw

ð4Þ

where subscripts f and m refer to the macropore and matrix
domains, respectively. Γw is the mass transfer term for waterexchange between the macropore and matrix domains [T− 1] andfw is the macropore domain fraction [−], that is, the ratio of thevolume of the macropore domain and the total ﬂow domain.
The matrix domain fraction is speciﬁed as one minus fw. Eq. (4) issolved using two sets of soil water retention and hydraulicconductivity functions, which are deﬁned for each domain. The
total soil water content or hydraulic conductivity are deﬁned as asum of the soil water contents or hydraulic conductivities ofeach domain multiplied by corresponding domain fractions. The
mass exchange between the matrix and macropore regions, Γw,
is calculated using the following equation:
Cw =



b
Ka γ w hf −hm
2
a

ð5Þ

Author's personal copy
110

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

where Ka is the effective saturated hydraulic conductivity of theinterface between the two pore domains [LT− 1]. Parametersdescribing aggregate shapes are the shape factor b [−] (=15 forspherical aggregates, =3 for cubic aggregates), the characteristic
length of an aggregate a [L] (a sphere radius or half size of the
cube edge), and the dimensionless scaling factor γw [−] (=0.4).
Analytical expressions proposed by van Genuchten (1980)
for the soil water retention curve, θ(h), and the hydraulicconductivity function, K(θ), are used in both models:θe =θðhÞ−θr1=m ; hb0θs −θr
1 + jαhjn

ð6Þθe = 1; hz0h  m i2K ðθÞ = Ks θle 1− 1−θ1=m
; hb0
e

ð7Þ

where θe is the effective soil water content [−], Ks is thesaturated hydraulic conductivity [LT− 1], θr and θs are the residual
and saturated soil water contents [L3L− 3], respectively, l is thepore-connectivity parameter [−], α is reciprocal of the air-entry
pressure head, [L− 1], and n [−] is related to the slope of theretention curve at the inﬂection point, and m = 1 − 1/n [−].
2.4.2. Single-porosity and dual-permeability solute transportmodelsModels for solute transport correspond conceptually to
the water ﬂow models described above. The single advection–
dispersion equation for solute transport is used for the singleporosity system:

ð8Þ

where c [ML− 3] and s [MM− 1] are solute concentrations in theliquid and solid phases, respectively, q is the volumetric ﬂux
density [LT− 1], ρd is the soil bulk density [ML− 3], D is thedispersion coefﬁcient [L2T− 1], and Φ describes zero- and ﬁrstorder rate reactions [ML− 3T− 1].
The dual-permeability formulation for solute transport isbased on two advection–dispersion equations:Aθf cf Aρd;f sf
Acf∂qf cf
A
Γs
θf Df
+
=
−
−Φf −
Az  Az
At
At
f
∂z
w

Aθm cm Aρd;m smA
Acm ∂qm cm
Γs
θm Dm
+
=
−
−Φm +
Az
At
At
Az
1−fw
∂z

ð9Þ

where subscripts f and m refer to the macropore and matrix
domains, respectively. The mass transfer term for soluteexchange between the macropore and matrix domains, Γs[ML− 3T− 1], is described using the following equation:
Cs =



b
De ð1−fw Þθm cf −cm + Cw cT
a2

s = kF cβð10Þ

where De is the effective diffusion coefﬁcient [L2T− 1] and c⁎ isequal either to cf for Γw N 0 or to cm for Γw b 0.

ð11Þ

where kF and β are empirical coefﬁcients. This equation isvalid for both single-porosity and dual-permeability models.
Nonequilibrium solute adsorption may also be simulatedusing the two-site sorption concept (van Genuchten andWagenet, 1989). The adsorbed concentration of solute on soil
particles may be deﬁned as:s = se + sk

K ðθÞ = Ks ; hz0Aθc Aρd s AAc ∂qc
+
=
θD
−
−Φ
At
At
Az
Az
∂z

Assuming equilibrium solute adsorption, the equilibrium
isotherm relating the sorbed concentration, s, to the aqueousconcentration, c, can be described using the Freundlich
equation:

ð12Þ

where se [MM− 1] denotes the adsorbed solute concentration onequilibrium sorption sites in the solid phase and sk [MM− 1]
represents the remaining solute concentration adsorbed to thesolid phase by a slow time-dependent process. The ﬁrst-order
rate equation for kinetic solute sorption is written as:Ask= ωc ð1−fs ÞkF cβ −sk
At

ð13Þ

where ωc [T− 1] is the ﬁrst-order sorption rate constant, fs [−] isthe fraction of equilibrium sorption sites, and kF and β areparameters of the Freundlich adsorption isotherm.
2.4.3. Estimation of soil hydraulic properties
The single-porosity model in HYDRUS-1D was ﬁrst used toanalyze transient inﬁltration data (cumulative inﬁltration, cumulative outﬂow and pressure heads at two positions) to obtain theparameters of both soil hydraulic functions (soil water retention
and unsaturated hydraulic conductivity functions) that were
described using the van Genuchten model (Eqs. (6) and (7)). Only
the initial part (up to 150 min) of the ﬁrst ponding inﬁltration
experiment was analyzed to reduce the impact of changes in thesoil structure and consequently in the porous system. Initialconditions within the soil sample were set based on pressure
heads measured at the beginning of the experiments. Pressure
heads at monitored locations were set to measured values, whilepressure heads between, above, and below were linearlyinterpolated and/or extrapolated. Since the soil samples werestored vertically in plastic bags at a cold place for 4 weeks before
experiments, the initial conditions were close to hydrostatic
conditions. The upper boundary condition was deﬁned using an
initial inﬁltration pulse of 0.116 cm min− 1 and a duration of 1 min,followed by a constant pressure head of 1.5 cm. The bottom
boundary condition was deﬁned using a seepage face.
While the saturated water content, θs, was set equal to themeasured porosity, the remaining soil hydraulic parameters(i.e., θr, α, n, and Ks) were optimized. In all cases the poreconnectivity parameter was assumed to be equal to an
average value for many soils (l = 0.5) (Mualem, 1976) since it
was found (not shown) that optimization of the l parameterincreased uncertainty for all optimized parameters.
The dual-permeability model of Gerke and van Genuchten(1993), as implemented in HYDRUS-1D, was applied next. Two
scenarios (DPM1 and DPM2) that were considered are described
below. To obtain unique optimization results for this complex
model, many of its parameters must be set equal toAuthor's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

111

Fig. 1. Images of soil-pore structure (left; frame length of 32.1 mm) and pore size and shape analysis (right) for various horizons of Haplic Luvisol.

independently estimated values. In both cases, the fraction of themacropore domain was estimated using the fraction of pores
identiﬁed from micromorphological images. The following
parameters were used to deﬁne the structure of both domains:

b=5 (the shape factor characterizing angular blocky aggregates),
a=0.2 cm, and γw =0.4. Saturated water contents of the matrix
and macropore domains, θs,m and θs,f, respectively, were assumed
to be the same and equal to the measured porosity. The residualAuthor's personal copy
112

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

water content of the macropore domain, θr,f, was set to zero. In the
DPM1 scenario, it was assumed that large pores control saturated
water ﬂow and that the Ks value obtained via numerical inversion
using the single porosity model can be used to derive thesaturated hydraulic conductivity of the macropore domain, Ks,f, asa ratio of Ks and the macropore domain fraction, fw. The remaining
parameters θr,m, αm, nm, Ks,m, αf, nf, and Ka were optimized. In the
DPM2 scenario, the nf parameter was set equal to a value of 2.5 todeﬁne the step-like shape retention curve that is characteristic for
macropores. The remaining parameters θr,m, αm, nm, Ks,m, αf, Ks,fand Ka were optimized. In all cases, parameters n were restricted
to vary between 1.1 and 3.
2.4.4. Estimation of solute transport parameters
Solute transport parameters were estimated using measured
outﬂow concentrations versus time for the single-porosity anddual-permeability models. In the case of the single-porosity

model, two scenarios (SPM-S1 and SPM-S2) were simulated.
Initial and boundary water ﬂow conditions were deﬁned as
before. Soil hydraulic parameters obtained from previous
inversions of water ﬂow data were ﬁxed. In addition, a soluteconcentration of 198 µg cm− 3 was speciﬁed for the initial
inﬁltration pulse (1 min). A zero concentration gradient wasdeﬁned at the bottom of the sample. Chlorotoluron adsorption
in different soil horizons was studied on undisturbed soil
samples using the standard laboratory procedure and adsorption isotherms were described using the Freundlich Eq. (11).
Bulk density was set to the measured value. Molecular diffusion
was assumed to be zero. In the SPM-S1 scenario, only thelongitudinal dispersivity was optimized. The two-site sorption
concept (van Genuchten and Wagenet, 1989) was applied in the
SPM-S2 scenario to characterize nonequilibrium solute adsorption. Due to the fact that experiments often do not provide
enough information to optimize all remaining parameters, theFig. 2. Images of soil-pore structure (left; frame length of 32.1 mm) and pore size and shape analysis (right) for various horizons of Greyic Phaeozem.

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

113

Fig. 3. Images of soil-pore structure (left; frame length of 32.1 mm) and pore size and shape analysis for (right) various horizons of Haplic Cambisol.

ωc parameter was set to 0.0001 min− 1, as found by Köhne et al.
(2006a,b). The longitudinal dispersivity and parameter fs were
optimized using measured transient solute transport data in the
SPM-S2 scenario.

Two scenarios (DPM-S1 and DPM-S2) were also consideredusing the dual-permeability model. Initial and boundary
condition were speciﬁed in the same way as for the singleporosity model. The fraction of surface ﬂow into the macropore

Table 1
Total image porosities and portion of different pore sizes
Soil type

HorizonImage
porosity

Portion [%] of different size classes [µm]b 100100–200

200–300

300–400

400–500

500–1000

N1000

Haplic Luvisol

Ap
Bt1
Bt2
Ck
Ap
BCk
Ck
Ap
Bw
C

0.04940.0826
0.05820.04750.03860.04240.05190.09200.03540.0396

6.3
4.9
5.4
6.7
7.1
7.2
6.0
6.7
7.1
15.8

18.5
11.6
14.4
16.3
22.4
20.4
14.1
16.4
18.1
27.6

17.4
10.9
15.8
19.4
12.6
22.4
15.6
13.8
16.8
21.0

12.2
10.7
13.5
14.2
16.6
14.8
13.9
18.1
11.7
14.1

9.4
7.5
11.4
7.3
6.7
10.2
11.7
6.6
11.9
7.0

24.4
18.8
23.7
20.9
22.0
13.9
25.0
23.1
17.1
14.5

11.8
35.6
15.8
15.2
12.6
11.1
13.5
15.3
17.3
0.0

Greyic Phaeozem

Haplic Cambisol

Author's personal copy
114

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

Table 2
Index of water stable aggregates (WSA) and soil properties affecting aggregate stability
Soil type

Horizon

WSA index [−]Organic carbon content [%]Clay content [%]Cation exchange capacity [mmol+(100 g)− 1]

Haplic Luvisol

Ap
Bt1
Bt2
Ck
Ap
Bth
BCk
Ck
Ap
Bw
C

0.48
0.61
0.30
0.05
0.56
0.60
0.45
0.10
0.75
0.46
0.39

1.11
0.50
0.32
0.17
1.11
0.50
0.46
0.29
1.77
0.41
0.37

24.4
34.7
33.7
24.0
32.6
37.3
31.0
27.6
15.0
24.6
23.5

14.0
21.2
25.5
16.5
20.5
29.2
22.5
15.5
14.2
13.5
15.0

Greyic Phaeozem

Haplic Cambisol

domain was deﬁned as a ratio of Ks,f multiplied by the domain
fraction and the total saturated hydraulic conductivity (the sumof saturated hydraulic conductivities of both domains multiplied by corresponding domain fractions). Soil hydraulicparameters obtained from the previous inversions (the DPM2
scenario) were again ﬁxed. In both scenarios, the two-site
sorption concept was applied to describe the adsorption
processes in the matrix domain. Parameters kF,m and βm were
set to measured values. Parameter ωc was again set equal to
0.0001 min− 1. The longitudinal dispersivity and the fs,mparameter were optimized using the measured transient solutetransport data and the numerical inversion. In the DPM-S1
scenario, similar to MACRO (Larsbo and Jarvis, 2003), the piston
type of solute transport was considered in the macropore
domain and the longitudinal dispersivity was set to a low valueof 0.01 cm. Chlorotoluron adsorption in the macropore domain
was neglected because of a much smaller surface area in contact
with solute than in the matrix domain and because of the shortcontact time in this fast ﬂowing region. In the DPM-S2 scenario,
the longitudinal dispersivity in the macropore domain was also
optimized. Molecular diffusion was neglected in both domains
and both scenarios. Chlorotoluron degradation, which may bedeﬁned using the ﬁrst-order rate reaction, was neglected in allcases due to the short experimental duration.
3. Results and discussion
3.1. Soil porous systems
Depending on the type of pedogenesis, soils exhibit
different porous systems. Micromorphological images (onefor each soil sample) and pore-size histograms (average forfour images) of all horizons of Haplic Luvisol, Greyic

Phaeozem, and Haplic Cambisol are shown in Figs. 1–3,
respectively. Total image porosities and fractions of different
pore sizes are given in Table 1. In all cases, the results
demonstrate the multimodality of soil porous systems.
Bimodality of soil porous systems detected on these images
is more or less evident for all soils and all soil horizons. A limit
value between two size class modes is 500 µm. Anotherdomain of considerably smaller pores is expected either inhigher order aggregates (in the top humic Ap, and deeper Bt1
and Bt2 horizons), or in the homogeneous matrix structure(the Ck, C, BCk and Bw horizons). Pores of different sizes may
be classiﬁed according to their impact on water ﬂow. A
deﬁnition of macropores as pores causing the physical nonequilibrium of water ﬂow close to saturation (frequently
under ponding conditions) is given in the review paper ofJarvis (2007). He suggested that macropores are pores with anequivalent diameter larger than about 300–500 µm, corresponding to pressure heads of −10 and −6 cm. On the other
hand, Kodešová et al. (2006, 2008) showed that even smaller
pores may play a signiﬁcant role. They deﬁned gravitational
pores as pores with an equivalent diameter larger than
1470 µm, corresponding to a pressure head of −2 cm (Watsonand Luxmoore, 1986), and large capillary pores as pores with adiameter larger than 40 µm (pores detectable on micromorphological images), corresponding to a pressure head of−70 cm, and smaller than 1470 µm. This deﬁnition was also
based on the bimodal shape of the soil water retention curve(Kodešová at al., 2007) obtained from multiple outﬂow
experiments. The macropore domain in the dual-permeability
model of HYDRUS-1D was deﬁned as a domain of large
capillary pores. Interestingly, the limiting value between two
size class modes documented in Figs. 1–3 is similar to thelimiting value for macropores as summarized by Jarvis (2007).

Table 3Inﬁltration durations, total volumes of inﬁltrating water and total volumes of water captured below the soil samples
Soil type

HorizonInﬁltration duration
1st day

Total inﬁltration
1st day (cm3)

Total outﬂow
1st day (cm3)Inﬁltration duration
2nd day

Total inﬁltration
2nd day (cm3)

Total outﬂow
2nd day (cm3)

HaplicLuvisol

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

6 h 54 min
0 h 30 min
4 h 42 min
5 h 5 min3 h 2 min6 h 4 min4 h 51 min

725 a (875) b
2400 (2550)
1275 (1450)
350 (550)
1225 (1400)
1705 (1855)
300 (425)

507 a (600) b
2092 (2228)
1098 (1199)
283 (405)
1129 (1264)
1307 (1470)
92 (207)

4 h 27 min
0 h 32 min
0 h 42 min
5 h 52 min
2 h 51 min
5 h 37 min
2 h 14 min

25 a (60) b
950 (1500)
425 (550)
100 (200)675 (850)
200 (325)
50 (175)

10 a (29) b
904 (1035)
430 (550)
57 (156)409 (540)
178 (301)
29 (136)

GreyicPhaeozem
HaplicCambisolab

Observed during the monitored inﬁltration period.
Total (during and after inﬁltration).

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

115

Fig. 4. Cumulative outﬂow from the bottom of soil samples versus chlorotoluron concentrations (left) and ﬁnal chlorotoluron concentration distributions (right) for
Haplic Luvisol (top), Greyic Phaeozem (middle), and Haplic Cambisol (bottom).

caused a signiﬁcant increase in fractions of pores with sizes
larger than 500 µm. In addition, many elongated pores wereobserved in this horizon. Their presence in the soil sample may
increase water ﬂux through the soil under saturated conditions.

However, an approach proposed by Kodešová et al. (2006) isapplied in this study.
The highest porosity in Haplic Luvisol (Fig. 1) was detected in
the Bt1 horizon due to a well-developed prismatic structure that

Table 4Chlorotoluron mass balances in solutes captured below the soil samples and remaining in soils
Soil type

Horizon

Chlorotoluron content
in outﬂow (µg)

Chlorotoluron content
in soil (µg)

Recovered chlorotoluron
content (µg)

Recoveredchlorotoluron (%)

Haplic Luvisol

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

857
1387
1484
671
1443
1251
672

818
360
283
788
265
522
1025

1675
1747
1767
1459
1708
1773
1697

84.6
88.3
89.2
73.7
86.3
89.5
85.7

Greyic Phaeozem
Haplic Cambisol

Author's personal copy
116

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

Table 5
Soil sample porosities, image porosities, and fraction of large pore and macropore domains
Soil type

HorizonBulk density[g cm− 3]Speciﬁc density[g cm− 3]

Porosity of soil
sample [−]Imageporosity [−]Fraction of large pore andmacropore domains, fw [−] aHaplic Luvisol

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

1.39
1.31
1.23
1.48
1.29
1.46
1.41

2.43
2.53
2.53
2.51
2.53
2.39
2.39

0.42800.48050.51190.41170.49470.38820.41270.04940.0826
0.05820.0386
0.05190.09200.0354

0.115
0.172
0.114
0.094
0.105
0.237
0.086

Greyic Phaeozem
Haplic CambisolaA ratio of the overall soil sample porosity and image porosity.

All pores are affected by clay coatings and inﬁllings (Kodešová
et al., 2006) that control the water ﬂow interaction between thelarger pores and pores of the matrix structure. The largest
porosity in Greyic Phaeozem (Fig. 2) was detected in the Ck
horizon, again due to an increase in fractions of pores with sizes
larger than 500 µm. While a majority of the soil pores wereregular or irregular, only a limited number of elongated pores
were observed. Lower water ﬂuxes under saturated condition
are expected in these soils. Images in Fig. 2 show the relatively
homogeneous matrix structure and a system of larger pores
created by roots and soil organisms. The highest porosity in
Haplic Cambisol (Fig. 3) was detected in the Ap horizon. Despite
the fact that aggregates in this horizon are poorly developed, theporous system does not show intrapedal or interpedal pores,
and pores are developed mainly along gravel particles, images
show many elongated pores of all detectable pore sizes.
Comparing all soil horizons (Table 1), the largest image porosity
was detected in the Ap horizon of Haplic Cambisol. A slightly
lower porosity was observed in the Bt1 horizon of Haplic Luvisol,
but a greater fraction of pores larger than 500 µm was also
detected. The lowest porosity was observed in the Bw and C
horizons of Haplic Cambisol and no pores of a size larger than
1000 µm were detected in the C horizon.
3.2. Soil structure stability
The indexes of water-stable aggregates, organic carbon
contents, clay contents and cation exchange capacities of all
horizons in all soil proﬁles are shown in Table 2. The most stable
aggregates in Haplic Luvisol were found in the Bt1 horizon,
presumably because of the presence of clay coatings in thishorizon. The stability of the soil structure of the Ap horizon(where the high organic carbon content should insure a high soil
aggregate stability) was probably reduced by the agricultural

practices. The impact of tillage and clay coatings on soil structurestability within the soil proﬁle is also evident in Greyic Phaeozem.
The highest aggregate stability was found in the Bth horizon. In
the case of Haplic Cambisol, the soil aggregate stability in the Ap
horizon was higher than in deeper horizons (Bw, C), due to higherorganic carbon content. A relatively high aggregate stability was
found despite micromorphological images showing weakly
developed soil aggregates. This was probably also caused by thepresence of free iron oxides in this soil.
3.3. Water ﬂow and solute transport experiments
The inﬁltration durations, total volumes of inﬁltratingwater, and total volumes of water captured below the soil
samples are shown in Table 3. Observed inﬁltration andoutﬂow rates partly reﬂected the soil structure of various
horizons (Figs. 1–3). The highest inﬁltration rate was observed
for the Bt1 horizon of Haplic Luvisol exhibiting a welldeveloped prismatic structure. A dense network of large
capillary pores with clay coatings and gravitational porescaused signiﬁcant preferential ﬂow. The irregular, nonuniform ﬂow path geometry for inﬁltration in structured soils
has also been discussed by German and Hensel (2006). The
lowest inﬁltration rate was measured for the Bw horizon of
Haplic Cambisol, which displayed a poorly developed soil
structure and had no gravitational pores. Lower inﬁltration
rates were observed for the humic Ap horizons of Haplic
Luvisol and Greyic Phaeozem than for their subhorizons. On
the other hand, a higher inﬁltration rate was monitored forthe humic Ap horizon than for the subhorizon of Haplic
Cambisol. Inﬁltration rates for humic horizons of all three soil
types increased due to the presence of macropores as follows:
Greyic Phaeozem, Haplic Luvisol, and Haplic Cambisol. In allcases, except for the Bt2 horizon, inﬁltration rates during theTable 6
Soil hydraulic parameters (van Genuchten, 1980) for the single-porosity model
Soil type
Haplic Luvisol

Greyic Phaeozem
Haplic CambisolaHorizon
Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

Fixed at 0 during optimization.
Not optimized.
c
95% conﬁdence interval.
bθr [cm3 cm− 3]a00.0102 ± 0.0413
0.0018 ± 0.0146
0.0056 ± 0.0011
0.0680 ± 0.0079
0a0.2415 ± 0.0264

θs [cm3 cm− 3]b0.42800.4805 b
0.5119 b
0.4117 b
0.4947 b
0.3882 b
0.4127 b

α [cm− 1]c0.0124 ± 0.0010
0.0185 ± 0.0206
0.0114 ± 0.0039
0.0599 ± 0.0191
0.0169 ± 0.00154
0.1035 ± 0.1550
0.0072 ± 0.0003n [−]Ks [cm min− 1]R22.001 ± 0.173
1.101 ± 1.260
1.221 ± 0.095
1.126 ± 0.009
1.093 ± 0.009
2.276 ± 1.401
3.000 ± 0.676

0.0377 ± 0.0031
0.8141 ± 0.0951
0.0694 ± 0078
0.0178 ± 0.0011
0.0872 ± 0.0023
0.0968 ± 0.0082
0.0069 ± 0.0006

0.94540.5036
0.46990.83020.92290.66150.9007

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

second inﬁltration decreased signiﬁcantly. The breakdown ofsoil structure, swelling of clay particles, and/or air entrapmentin the soil may cause this decrease in inﬁltration rates. The
restriction of water ﬂow due to the presence of entrapped air
has already been described by Císlerová et al. (1988), as wellas others. The largest decrease in inﬁltration rate was117

observed for the Ap horizon of Haplic Luvisol, in which
breakdown of aggregates led to the formation of the lesspermeable layer at the surface of the soil sample. Thisoccurred despite the fact that the WSA of the Ap horizon(Table 2) was higher than for the Bt2 and Ck horizons of thesame soil proﬁle. A homogeneous matrix structure with a

Fig. 5. Cumulative inﬁltration (Inﬂow), cumulative outﬂow (outﬂow), pressure heads at upper (UT) and lower (LT) tensiometers measured and simulated using thesingle-porosity (SPM) and dual-permeability (DPM) models (scenario DPM2) for Ap (left), Bt1 (middle) and Bt2 (right) horizons of Haplic Luvisol.

Author's personal copy
118

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

system of larger pores created by roots and soil organisms indeeper horizons was apparently more stable than the soil
aggregates of the Ap horizon.
The relationship between the monitored outﬂow andoutﬂow chlorotoluron concentrations for all soil samples isshown in Fig. 4. Peak solute concentrations were observed in allcases immediately at or close after the initiation of outﬂow.
Faster inﬁltration rates usually resulted in lower peak concentrations and a more gradual decrease of measured values. A
dense network of large capillary and gravitational pores is likely
to have partly bypassed initial chlorotoluron deposits, whichwere then gradually released into the macropores by mass

Fig. 6. Cumulative inﬁltration (Inﬂow), cumulative outﬂow (outﬂow), pressure heads at upper (UT) and lower (LT) tensiometers measured and simulated using thesingle-porosity (SPM) and dual-permeability (DPM) models (scenario DPM2) for Ap (left) and Ck (right) horizons of Greyic Phaeozem.

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125transfer. Outﬂow concentrations observed during the secondday inﬁltration experiments decreased signiﬁcantly comparedto those measured at the end of the ﬁrst-day inﬁltration tests.
This was likely due to kinetic adsorption processes and masstransfer of chlorotoluron to less mobile zones. Interestingly,
chlorotoluron was detected even in the outﬂow from the Bw

119

horizon of Haplic Cambisol, where a very low outﬂow wasmeasured.
Final distributions of chlorotoluron concentrations in all
soil samples are shown in Fig. 4. The highest chlorotoluron
concentrations were found in the surface layers of the Ap
horizons of all three soils. A similar phenomenon was also

Fig. 7. Cumulative inﬁltration (Inﬂow), cumulative outﬂow (outﬂow), pressure heads at upper (UT) and lower (LT) tensiometers measured and simulated using thesingle-porosity (SPM) and dual-permeability (DPM) models (scenario DPM2) for Ap (left) and Bw (right) horizons of Haplic Cambisol.

Author's personal copy120

Table 7
Soil hydraulic parameters (van Genuchten, 1980) for the dual-permeability model (scenario DPM1)
Horizonθr,m [cm3 cm− 3]

θs,m = θs,f [cm3 cm− 3]αm [cm− 1]

nm [−]

Ks,m [cm min− 1]αm [cm− 1]

nf [−]

Ks,f [cm min− 1]

Ka [cm min− 1]

R2

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

0.1209 ± 0.1310 a
0.1359 ± 0.0172
0.1081 ± 0.4220
0.1593 ± 0.0269
0.2003 ± 0.3610
0.2163 ± 0.0602
0.2500 ± 0.0149

0.4280 b
0.4805 b
0.5119 b
0.4117 b
0.4947 b
0.3882 b
0.4127 b

0.0086 ± 0.0010
0.0202 ± 0.0025
0.0102 ± 0.0142
0.0296 ± 0.0092
0.0051 ± 0.0055
0.0267 ± 0.0073
0.0096 ± 0.0009

2.500 ± 1.670
1.319 ± 0.176
1.293 ± 0.295
1.100 ± 0.086
1.100 ± 0.013
2.501 ± 0.802
1.851 ± 0.165

0.0063 ± 0.0026
0.0015 ± 0.0002
0.0522 ± 0.0085
0.00252 ± 0.0002
0.0022 ± 0.0027
0.0051 ± 0.0029
0.0004 ± 0.0001

0.0100 ± 0.0026
0.052 ± 0.0053
0.0137 ± 0.1700
0.0319 ± 0.0078
0.0462 ± 0.0042
0.0101 ± 0.0148
0.0103 ± 0.0009

3.000 ± 5.080
1.190 ± 0.137
3.000 ± 2.58
2.851 ± 0.389
2.191 ± 1.800
3.000 ± 0.890
3.000 ± 0.230

0.3285 b
9.4704 b
0.6080 b
0.1900 b
0.8305 b
0.4080 b
0.0802 b

8.68E−06 ± 5.36E−06
9.49E−06 ± 1.02E−06
3.89E−06 ± 1.25E−04
1.08E−05 ± 3.00E−06
6.05E−06 ± 3.38E−06
2.03E−07 ± 1.46E−07
3.58E−05 ± 2.43E−06

0.91550.6457
0.46270.20140.96290.84210.9212

Greyic Phaeozem
Haplic Cambisolab

95% conﬁdence interval.
Not optimized.

Table 8
Soil hydraulic parameters (van Genuchten, 1980) for the dual-permeability model (scenario DPM2)
Soil type

Horizonθr,m [cm3 cm− 3]

θs,m = θs,f [cm3 cm− 3]αm [cm− 1]

nm [−]

Ks,m [cm min− 1]αm [cm− 1]

n f [−]

Ks,f [cm min− 1]

Ka [cm min− 1]

R2

Haplic Luvisol

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

0.1749 ± 0.0185 a
0.1588 ± 0.0177
0.1001 ± 0.0391
0.1875 ± 0.0083
0.1735 ± 0.0688
0.1039 ± 0.1130
0.2500 ± 0.0075

0.4280 b
0.4805 b
0.5119 b
0.4117 b
0.4947 b
0.3882 b
0.4127 b

0.0225 ± 0.0027
0.0047 ± 0.0001
0.0102 ± 0.0005
0.0467 ± 0.0024
0.0044 ± 0.0002
0.0267 ± 0.0047
0.0265 ± 0.0014

2.418 ± 0.267
1.102 ± 0.018
1.294 ± 0.297
1.100 ± 0.033
1.104 ± 0.047
2.505 ± 0.920
1.635 ± 0.088

0.0233 ± 0.0029
0.0150 ± 0.0021
0.0528 ± 0.0399
0.0031 ± 0.0001
0.0030 ± 0.0023
0.0085 ± 0.0043
0.0002 ± 0.0000

0.0107 ± 0.0015
0.051 ± 015
0.0107 ± 0.0102
0.0589 ± 0.0087
0.0226 ± 0.0004
0.0106 ± 0.0139
0.0108 ± 0.0003

2.5 b
2.5 b
2.5 b
2.5 b
2.5 b
2.5 b
2.5 b

0.1610 ± 0.0205
4.5775 ± 0.1491
0.1050 ± 0.2650
0.1370 ± 0.0026
0.7770 ± 0.025
0.3594 ± 0.0237
0.0929 ± 0.0036

5.32E−05 ± 6.45E−06
1.15E−07 ± 1.24E−08
1.93E−06 ± 7.73E−05
9.15E−06 ± 1.16E−06
3.57E−06 ± 1.03E−06
3.21E−07 ± 2.65E−07
1.86E−05 ± 9.22E−07

0.9496
0.89400.4711
0.6632
0.97100.84140.8738Greyic Phaeozem
Haplic Cambisolab

95% conﬁdence interval.
Not optimized.

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

Soil type
Haplic Luvisol

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125
Table 9Parameters of Freundlich adsorption isotherms determined in the laboratory
Soil type

HorizonFreundlich adsorption
isotherm coefﬁcient, kF[cm3β µg1− β g− 1]

Freundlich
adsorption isotherm
coefﬁcient, β [–]HaplicLuvisol

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

2.89
0.86
0.84
2.73
0.88
4.77
0.52

0.80
0.84
0.83
0.75
0.81
0.77
0.85

GreyicPhaeozem
HaplicCambisolobserved in the ﬁeld by Kočárek et al. (2005). This may be
explained by solute diffusion into the soil aggregates
immediately after the solute application, and a subsequentslow solute release from aggregates into the preferential ﬂow
paths, or even solute isolation in the aggregates. Chlorotoluron concentrations in top layers of subsurface horizons
were lower than those in the Ap horizons and chlorotoluron
was more evenly distributed with depth. Chlorotoluron
concentrations within the Bw horizon of Haplic Cambisol
were very high due to very low outﬂow from the soil sample.
Despite this low outﬂow, chlorotoluron was regularly distributed within the soil proﬁle.
Chlorotoluron amounts captured in the outﬂow from soil
samples and remaining in soils are presented in Table 4. As
expected, higher chlorotoluron contents in outﬂow were foundfor samples with higher cumulative outﬂow. The highest
chlorotoluron content remaining in the soil was found for the
Bw horizon of Haplic Cambisol, due to very low outﬂow.
Chlorotoluron recovery varied for different soil samples. The
highest chlorotoluron recovery was found for soil samples withhigh inﬁltration rates. Chlorotoluron recovery was partly
affected by pesticide degradation, which was assumed to be
between 5 and 8% of the applied mass. Pesticide losses during
sample segmentation, soil sample grinding and sieving,incomplete pesticide dissolution in methanol, heterogeneityof applied solute, and measurement errors are other likely
reasons for incomplete chlorotoluron recovery.
3.4. Numerical evaluation of experimental data
3.4.1. Estimation of soil hydraulic properties
Soil sample porosities evaluated from bulk densities(calculated for the entire 1125.1-cm3 undisturbed soil sample)
and speciﬁc densities are presented in Table 5. Soil hydraulicparameters obtained using numerical inversions of the transient ﬂow data and the single-porosity model and R2 of

121

measured and simulated data are shown in Table 6. Measured
transient outﬂow data (Figs. 5–7) were, in most cases, notsatisfactorily described by the model. In all cases, water outﬂow
was observed almost immediately after the beginning of the
ponding inﬁltration. The later wetting front breakthrough wassimulated especially in all Ap horizons. While simulated
inﬁltration and outﬂow rates remained constant after thewetting front breakthrough, observed inﬁltration and outﬂow
rates decreased with time due to the breakdown of the soil
structure and changes in particle volumes. As discussed above,the soil structure breakdown (occurred as the soil surface
collapsed) was most evident in the soil sample taken from theAp horizon of Haplic Luvisol. In several cases (especially thosewith the lowest R2), while simulated pressure heads indicated afast advance of the wetting front, the observed data demonstrated a much slower wetting front movement. Measured
pressure heads likely characterized pressure heads in thematrix domain that may represent either immobile domains
or domains with very low permeability. Conﬁdence intervals(95%) for the optimized Ks parameter show that this parametercan be estimated from the inﬁltration experiment highly
reliably. On the other hand, 95% conﬁdence intervals for someother optimized parameters were very wide, indicating thatobserved data were affected by the factors discussed above, andthat the single porosity model was not suitable to describe theﬂow process in studied soils.
The soil hydraulic parameters for all soils obtained using thenumerical inversion of the dual-permeability model and R2 ofmeasured and simulated data are shown in Tables 7 and 8 forscenarios DPM1 and DPM2, respectively. As discussed above,
fractions of large capillary pore and gravitational pore domains(Table 5) were deﬁned as a ratio of the image porosity and thesoil sample porosity. R2 values in Table 8 show that correlations
between measured and simulated data are better for Haplic
Luvisol and Greyic Phaeozem when using the DPM2 scenario
and slightly better for Haplic Cambisol when using the DPM1
scenario. However, 95% conﬁdence intervals show that optimized parameters are more reliable for the DPM2 scenario in allcases. Similarly as for the single-porosity model, dynamic
parameters (i.e., Ks,m and Ks,f) are determined better than otheroptimized parameters. Despite the large number of optimized
parameters this experiment provides enough information about
the saturation water ﬂuxes. Another factor that increases theuncertainty of optimized parameters, apart from those alreadydiscussed earlier, is the fact that the experiments were
performed under conditions close to soil saturation, andconsequently the wetting front moved very fast. Pressure
heads ranged from −250 to 0 cm during all inﬁltration
experiments, except for the Bt2 horizon of Haplic Cambisol

Table 10
Solute transport parameters obtained using the numerical optimization with the single-porosity model
Scenario SPM-S1
Soil type
Haplic Luvisol

Greyic Phaeozem
Haplic Cambisol

Horizon
Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

Longitudinal dispersivity [cm]40.000 ± 0.004
40.000 ± 0.010
40.000 ± 0.016
40.000 ± 540.0
40.000 ± 0.040
7.454 ± 0.094
26.408 ± 0.022

Scenario SPM-S2
R2
0.3916
0.84180.00590.17600.00830.22480.7728Longitudinal dispersivity [cm]40.000 ± 14.00
38.242 ± 9.942
40.000 ± 0.022
36.523 ± 321.10
17.751 ± 15.500
40.000 ± 5.850
4.084 ± 4.510

The fraction of equilibrium sorption sites, fs [−]0.140 ± 0.172
0.826 ± 0.121
0.340 ± 0.509
0.039 ± 0.419
0.109 ± 0.123
0.900 ± 0323
0.117 ± 3.102

R2
0.44200.86070.03330.17680.15270.58970.9026

Author's personal copy
122

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

Table 11
Solute transport parameters obtained using the numerical optimization with the dual-permeability model
Scenario DPM-S1

Scenario DPM-S2
2

Soil type

HorizonLongitudinal
dispersivity inmatrix [cm]The fraction of equilibrium
sorption sites in matrix,
fs,m [−]

R

HaplicLuvisol

Ap
Bt1
Bt2
Ap
Ck
Ap
Bw

30.70 ± 19.30
5.620 ± 4.740
1.828 ± 3.311
13.373 ± 66.122
40.000 ± 2.990
0.365 ± 4.520
2.592 ± 5.361

0.280 ± 0.379
0.499 ± 0.243
0.498 ± 0432
0.473 ± 0.593
0.212 ± 0.114
0.016 ± 0.009
0.334 ± 0.548

0.44920.52580.09820.01690.77320.10220.9357

GreyicPhaeozem
HaplicCambisolwhere they ranged from −700 to 0 cm. Such data may hardly
provide enough information to optimize the residual watercontent, and in some cases (having data for only relatively
narrow experimental pressure head interval without information about the dry range), even the optimization of α and nparameters (characterizing the shape of soil hydraulic functions)
may be difﬁcult.
A comparison of correlations between measured and simulated data using either the single-porosity or dual-permeability
models shows that the R2 values did not, in some cases, increase
for the dual-permeability model. The R2 value for the Ap horizonof Greyic Phaeozem due to the applied restrictions evensigniﬁcantly decreased. However, Figs. 5–7 show that outﬂows/
inﬂows simulated using the DPM2 scenario describe observed
trends better, such as earlier wetting front breakthrough. Better
correlations between measured and simulated pressure heads
were not always achieved due to the reasons previously
discussed. In addition, the bimodality of the large capillary poreand the gravitational pore fraction detected in micromorphological images, which were not considered in numerical simulations, may have had a signiﬁcant impact on water ﬂow within soil
samples. Dividing a detected pore-ﬂow domain into two domains(a domain of pores with an equivalent diameter larger than
500 µm and a domain of pores with an equivalent diametersmaller than 500 µm), each with its own soil hydraulic properties,
would most likely increase the simulated inﬁltration rate in thedomain of large pores and cause earlier wetting front break-Longitudinal
dispersivity inmatrix [cm]The fraction of equilibrium
sorption sites in matrix,
fs,m [−]

Longitudinal
dispersivity inmacropores [cm]R2

30.00 ± 24.70
4.010 ± 6.412
8.019 ± 2.000
5.670 ± 116.604
40.000 ± 35.100
0.538 ± 7.390
3.589 ± 1.990

0.165 ± 0.395
0.504 ± 0.894
0.097 ± 1.110
0.502 ± 13.135
0.215 ± 0.117
0.021 ± 0.164
0.505 ± 1.583

5.294 ± 7.143
0.018 ± 0.056
1.399 ± 8.808
0.014 ± 19.457
0.023 ± 0.007
0.005 ± 0.199
0.0003 ± 0.038

0.46650.59030.63300.01970.79560.11050.9396through. However, such a scenario cannot be simulated using theapplied simulation program. It must be also pointed out that thesize of the thin soil sections presumably allows for studying
capillary pores. Thin soil sections would have to be larger in orderto study larger gravitational pores. In addition, a thin soil section
provides only a two-dimensional image of the soil porous system.
The accidental occurrence of larger gravitational pores and theirshape makes it difﬁcult to study them using this technique. The
impact of gravitational pores on water ﬂow, which were notdetected on micromorphological images, is the most evident in
the Ap horizon of Greyic Phaeozem. The inﬂuence of gravitational
pores signiﬁcantly increases in soils with lower fractions of poresdetectable on micromorphological images. On the other hand,
the ﬂow conditions similar to ponding inﬁltration, which cause
water ﬂow through gravitational pores, are not frequently
observed in nature. Therefore, soil micromorphological images
should provide sufﬁcient information for characterizing ﬂow
domains, model selection and deﬁnition of soil hydraulicproperties under the ﬂow conditions prevailing in nature(Kodešová et al., 2006, 2008).
3.4.2. Estimation of solute transport parameters
The parameters of the Freundlich adsorption isotherm usedin all simulations were independently determined in thelaboratory using batch experiments (Table 9). Additional
optimized solute transport parameters and R2 coefﬁcientsobtained using the single-porosity and dual-permeability

Fig. 8. Outﬂow concentrations measured and simulated using the single-porosity (SPM) (scenarios SPM-S1 and SPM-S2) and dual-porosity (DP) (scenarios DPM-S1
and DPM-S2) models for Ap (left), Bt1 (middle) and Bt2 (right) horizons of Haplic Luvisol.

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

123

Fig. 9. Outﬂow concentrations measured and simulated using the single-porosity (SPM) (scenarios SPM-S1 and SPM-S2) and dual-porosity (DPM) (scenarios DPMS1 and DPM-S2) models for Ap (left) and Ck (right) horizons of Greyic Phaeozem.

models are given in Tables 10 and 11, respectively. Observed andsimulated chlorotoluron concentrations for all numericalscenarios and all soil horizons are presented in Figs. 8–10. It isevident that the SPM-S1 scenario (the single porosity modelwith optimized value of the longitudinal dispersivity) could not
simulate observed outﬂow concentrations for the Ap horizonsof all soil types due to the high prescribed sorption (Table 9),
and was able to predict only concentration tailing for deeper
horizons. Using the concept of two-site sorption improved
correspondence between measured and simulated data fordeeper horizons. However, it was possible to predict the early
chlorotoluron breakthrough only using the dual-permeability
model. As with water ﬂow, although the R2 values did notalways indicate improved correlations between measured andsimulated data when the dual-permeability model was used,
Figs. 8–10 show that the dual-permeability model better
described the general trends exhibited by measured concentrations. However, simulated peak concentrations were almost
always lower than measured values. Simulated tail concentrations were usually either similar or slightly higher than
measured values. A comparison of modelling results obtainedby the two dual-permeability scenarios reveals that the DPM-S2
scenario provided better correlations between measured andsimulated data, mainly due to the increase in the number of
optimized parameters. However, in many cases the value of the
optimized longitudinal dispersivity was very low, and similar tothe ﬁxed value in the DPM-S1 scenario. Note that several dual-

permeability models, such as MACRO (Larsbo and Jarvis, 2003),
assume piston-type solute transport in the macropore domain
and neglect dispersion in this domain altogether. Again,
conﬁdence intervals of many optimized parameters were
often very wide, indicating their nonuniqueness and mutual
correlations. Changes in soil structure and bimodality of themacropore domain hampered attempts to ﬁnd a unique set ofsolute transport parameters that would describe collectedtransient data. Since it was shown by Köhne et al. (2006b) that
adsorption in the macropore domain should be assumed in thehorizons containing organic matter, chlorotoluron adsorption
in macropores was also simulated (not shown). However, such
simulations led to either no correlation between measured andsimulated data, or the two-site sorption concept had to be
applied again with a very low value of the fraction ofequilibrium sorption sites in macropores. In addition, observed
data do not provide enough information to optimize so many
parameters. Finally, macropores in the subsurface horizons of
Haplic Luvisol and the Ap horizon of Greyic Phaeozem wereaffected by clay coatings that may decrease the pesticide
sorption (Celis et al., 1997).
This study showed that even in soils with poorly or
undeveloped soil aggregates and without larger gravitational
pores (the Bw horizon of Haplic Cambisol) the small fracture oflarge capillary pores detected in micromorphological images
signiﬁcantly inﬂuenced solute transport. The impact of thesepores was previously documented by Kodešová et al. (2008), by

Fig. 10. Outﬂow concentrations measured and simulated using the single-porosity (SPM) (scenarios SPM-S1 and SPM-S2) and dual-porosity (DPM) (scenarios
DPM-S1 and DPM-S2) models for Ap (left) and Bw (right) horizons of Haplic Cambisol.

Author's personal copy
124

R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125

comparing simulated pesticide distribution within the soil proﬁlewith pesticide distribution studied experimentally in the ﬁeld.
This assumption was proved here via numerical inversion oflaboratory experimental data.
4. Conclusions
Soil pore structure, soil aggregate stability, water ﬂow, andchlorotoluron transport were studied in various diagnostic
horizons of three soil types. Observed inﬁltration and outﬂow
rates and pesticide transport closely reﬂected the soil structureof the studied horizons. Transport processes were also affected
to a certain degree by the breakdown of soil structure. The
highest water ﬂow rates were observed in soils with welldeveloped soil structures, which were affected by clay coatingsand high aggregate (structure) stability. The lowest water ﬂow
rates were observed in soils with poorly developed soil
structures, a low fraction of large capillary pores, and theabsence of gravitational pores. Water ﬂux was affected signiﬁcantly during the experiment by the aggregate breakdown in
the well-aggregated soil with moderate aggregate stability. On
the other hand, structural changes were less signiﬁcant incompact homogeneous matrix structures with systems of larger
pores created by roots and soil organisms and low structurestability. The dual-permeability model improved the descriptionof collected data in comparison to the single-porosity model.
However, better correlations were not always obtained due tothe breakdown of soil structure and particle volume changes
during experiments. In addition, the bimodality of the large poreand gravitational pore fraction detected in micromorphologicalimages, which was not considered in numerical simulations,
may have had a signiﬁcant impact on water ﬂow in the soil
samples. Utilization of the dual-permeability model and theconcept of a two-site sorption applied to the matrix domain
made it possible to ﬁt the measured data more closely than when
the single-porosity model with equilibrium adsorption wasused. This study conﬁrmed that chlorotoluron is highly mobile
under saturated conditions in the deeper horizons of Haplic
Luvisol and Greyic Phaeozem. Interestingly, the parent material
for both soils was loess, which is usually assumed to be a lowpermeability material. On the other hand, Cambisols are usuallyassumed to be soils with high permeability due to the presenceof preferential pathways. This assumption was not conﬁrmed fordeeper horizons of the studied Cambisol. However, the prevailing fast chlorotoluron transport through large capillary pores,
and frequently through gravitational pores, was documented in
all cases. Micromorphological images and soil porous systemanalysis signiﬁcantly improved the understanding of transport
processes in studied soils.
Acknowledgement
Authors acknowledge the ﬁnancial support of the Grant
Agency of the Czech Republic (grant No. GA CR 526/08/0434)
and the Ministry of Education, Youth and Sports (grants No.
2B06095 and No. MSM 6046070901). The thin section
preparation and the image analysis were funded by the Italian
Agricultural Ministry (grant No. D.M. 2/7331/06). We also thankto Anna Žigová and Karel Němeček for their help with the ﬁeld
work. The comments of anonymous reviewers and editors
A. Coppola and E.O. Frind are gratefully acknowledged.

References
Bouma, J., Jongerius, A., Boersma, O., Jager, A., Schoonderbeek, D., 1977. The
function of different types of macropores during saturated ﬂow throughfour swelling soil horizons. Soil Sci. Soc. Am. J. 41, 945–950.
Bouma, J., Jongerius, A., Schoonderbeek, D., 1979. Calculation of saturated
hydraulic conductivity of some pedal clay soils using micromorphometric data. Soil Sci. Soc. Am. J. 43, 261–264.
Bronick, C.J., Lal, R., 2005. Soil structure and management: a review. Geoderma
124, 3–22.
Celis, R., Cox, L., Hermosin, M.C., Cornejo, J., 1997. Sorption of thiazaﬂuron by ironand humic acid-coated montmorillonite. J. Environ. Qual. 26 (2), 472–479.
Císlerová, M., Šimůnek, J., Vogel, T.,1988. Changes of steady-state inﬁltration rates
in recurrent ponding inﬁltration experiment. J. Hydrol. 104 (1–4), 1–16.
Coppola, A., et al., 2009. Darcian preferential water ﬂow and solute transport
through bimodal porous systems: Experiments and modelling, J. Contam.
Hydrol. 104, 74–83 (this issue).
Dousset, S., Thevenot, D., Pot, V., Šimůnek, J., Andreux, F., 2007. Evaluating
equilibrium and non-equilibrium transport of bromide and isoproturon in
disturbed and undisturbed soil columns. J. Contam. Hydrol. 94, 261–276.
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30, 211–233.
Gärdenäs, A., Šimunek, J., Jarvis, N., van Genuchten, M. Th., 2006. Twodimensional modelling of preferential water ﬂow and pesticide transport
from a tile-drained ﬁeld. J. Hydrol. 329, 647–660.
Gerke, H.H., 2006. Review article. Preferential ﬂow descriptions forstructured soils. J. Plant Nutr. Soil Sci. 169, 382–400.
Gerke, H.H., van Genuchten, M.Th., 1993. A dual-porosity model for
simulating the preferential movement of water and solutes in structured
porous media. Water Resour. Res. 29, 305–319.
German, P.F., Hensel, D., 2006. Poiseuille ﬂow geometry inferred fromvelocities of wetting front in soils. Vadose Zone J. 5, 867–876.
Jarvis, N., 2007. Review of non-equilibrium water and solute transport in soil
macropores: principles, controlling factors and consequences for water
quality. Eur. J. Soil Sci. 58, 523–546.
Kočárek, M., Kodešová, R., Kozák, J., Drábek, O., Vacek, O., 2005. Chlorotoluron
behaviour in ﬁve different soil types. Plant, Soil Environ. 51 (7), 304–309.
Kodešová, R., Kozák, J., Šimůnek, J., Vacek, O., 2005. Single and dualpermeability model of chlorotoluron transport in the soil proﬁle. Plant
Soil Environ. 51 (7), 310–315.
Kodešová, R., Kodeš, V., Žigová, A., Šimůnek, J., 2006. Impact of plant roots
and soil organisms on soil micromorphology and hydraulic properties.
Biologia (Bratisl.) 61 (Supl. 19), S339–S343.
Kodešová, R., Pavlů, L., Kodeš, V., Žigová, A., Nikodem, A., 2007. Impact of
spruce forest and grass vegetation cover on soil micromorphology andhydraulic properties of organic matter horizon. Biologia (Bratisl.) 62 (5),565–568.
Kodešová, R., Kočárek, M., Kodeš, V., Šimůnek, J., Kozák, J., 2008. Impact of soil
micromorphological features on water ﬂow and herbicide transport in soils.
Vadose Zone J. 7 (2), 798–809 Special Issue “Vadose Zone Modeling”.
Köhne, J.M., Köhne, S., Mohanty, B.P., Šimůnek, J., 2005. Inverse mobile–immobile modeling of transport during transient ﬂow: effect of betweendomain transfer and initial content. Vadose Zone J. 3, 1309–1321.
Köhne, J.M., Mohanty, B.P., Šimůnek, J., 2006a. Inverse dual-permeability
modeling of preferential water ﬂow in a soil column and implications forﬁeld-scale solute transport. Vadose Zone J. 5, 59–76.
Köhne, J.M., Köhne, S., Šimůnek, J., 2006b. Multi-process herbicide transport
in structured soil columns: experiments and model analysis. J. Contam.
Hydrol. 85 (1–2), 1–32.
Köhne, J.M., Köhne, S., Šimůnek, J., 2009a. A Review of Model Applications for
Structured Soils: a) Water Flow and Tracer Transport, J. Contam. Hydrol.
104, 4–35 (this issue).
Köhne, J.M., Köhne, S., Šimůnek, J., 2009b. A Review of Model Applications for
Structured Soils: b) Pesticide Transport, J. Contam. Hydrol. 104, 36–60(this issue).
Larsbo, M., Jarvis, N., 2003. MACRO 5.0. A model of water and solute transport
in macroporous soil. Technical description. Swedish Univ. of Agricult. Sci.,
Dept. of Soil Sci., Uppsala, Swedem, Emergo, 6. 48 pp.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten, M.Th., 1997.
New piecewise-continuous hydraulic functions for modeling preferential
ﬂow in an intermittent ﬂood-irrigated ﬁeld. Water Resour. Res. 33,
2049–2063.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12 (3), 513–522.
Nimmo, J.R., Perkins, K.S., 2002. Aggregate stability and size distribution. In:
Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis, Part 4 — Physical
Methods. Soil Science Society of America, Inc., Madison, USA, pp. 317–328.
Othmer, H., Diekkrüger, B., Kutílek, M., 1991. Bimodal porosity andunsaturated hydraulic conductivity. Soil Sci. 152 (3), 139–150.

Author's personal copy
R. Kodešová et al. / Journal of Contaminant Hydrology 104 (2009) 107–125
Pagliai, M., La Marca, M., Lucamante, G., 1983. Micromorphometric andmicromorphological investigations of a clay loam soil in viticulture underzero and conventional tillage. J. Soil Sci. 34, 391–403.
Pagliai, M., Marsili, A., Servadio, P., Vignozzi, N., Pellegrini, S., 2003. Changes
in some physical properties of clay soil in central Italy following thepassage of rubber tracked and wheeled tractors of medium power. Soil
Tillage Res. 73, 119–129.
Pagliai, M., Vignozzi, N., Pellegrini, S., 2004. Soil structure and the effect ofmanagement practices. Soil Tillage Res. 79, 131–143.
Pot, V., Šimůnek, J., Benoit, P., Coquet, Y., Yra, A., Martinez Cordun, M.-J., 2005.
Impact of rainfall intensity on the transport of two herbicides in undisturbed
grassed ﬁlter strip soil cores. J. Contam. Hydrol. 81 (1–4), 63–88.
Servadio, P., Marsili, A., Vignozzi, N., Pellegrini, S., Pagliai, M., 2005. Effect onsome qualities in central Italy following the passage of four wheel drive
tractor ﬁtted with single and dual tires. Soil Tillage Res. 84 (1), 87–100.
Šimůnek, J., Jarvis, N.J., van Genuchten, M. Th., Gärdenäs, A., 2003. Review andcomparison of models for describing non-equilibrium and preferential
ﬂow and transport in the vadose zone. J. Hydrol. 272, 14–35.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th., 2008a. The
HYDRUS-1D Software Package for Simulating the Movement of Water,
Heat, and Multiple Solutes in Variably Saturated Media, Version 4.0,

125

HYDRUS Software Series 3, Department of Environmental Sciences,
University of California Riverside, Riverside, California, USA, pp. 315.
Šimůnek, J., van Genuchten, M. Th., 2008b. Modeling nonequilibrium ﬂow
and transport with HYDRUS. Vadose Zone J. 7 (2), 782–797 Special Issue“Vadose Zone Modeling”.
Six, J., Elliott, E.T., Paustian, K., Doran, J.W., 1998. Aggregation and soil organic
matter accumulation in cultivated and native grassland soils. Soil Sci. Soc.
Am. J. 62, 1367–1377.
Six, J., Feller, C., Denef, K., Ogle, S.M., 2002. Soil organic matter, biota andaggregation in temperate and tropical soils — effect of no-tillage. Agronomie
22, 755–775.
UMS GmbH, Munich, 2005. T5 Miniature Pressure Transducer Tensiometer,
User Manual version 1.8, Munich.
van Genuchten, M. Th., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Genuchten, M.T., Wagenet, R.J., 1989. Two-site/two-region models forpesticide transport and degradation: theoretical development andanalytical solutions. Soil Sci. Soc. Am. J. 53, 1303–1310.
Watson, K.W., Luxmoore, R.J.,1986. Estimating macroporosity in a forest watershedby use of a tension inﬁltrometer. Soil Sci. Soc. Am. J. 50 (3), 578–582.