M
Z

Impact of Soil Micromorphological
Features on Water Flow andHerbicide Transport in Soils

S

S

:V

Radka Kodešová,* Mar n Kočárek, Vít Kodeš, Jiří Šimůnek, and Josef Kozák
The impact of varying soil micromorphology on soil hydraulic proper es and, consequently, on water flow and herbicidetransport observed in the field is demonstrated on three soil types. The micromorphological image of a humic horizonof Haplic Luvisol showed higher-order aggregates. The majority of detectable pores corresponding to the pressure
head interval between −2 and −70 cm were highly connected, separa ng higher-order peds with small intrapores thatpossibly formed zones with immobile water. Herbicide was regularly distributed in this soil. The majority of detectablelarge capillary pores in a humic horizon of Greyic Phaeozem were separated and aﬀected by clay coa ngs and fillings.
The herbicide transport in this soil was highly aﬀected by preferen al flow. Macropores corresponding to pressure
heads higher than −2 cm were detected in a humic horizon of Haplic Cambisol. However, preferen al flow only slightly
influenced the herbicide transport in this soil. Single-porosity and either dual-porosity or dual-permeability flow andtransport models in HYDRUS-1D were used to es mate the soil hydraulic parameters from laboratory mul step ou low
and ponded infiltra on experiments via numerical inversion and to simulate the herbicide transport experimentallystudied in the field. Appropriate models were selected on the basis of the soil micromorphological study.

A: WF, weighting factor.

W

and solute transport are traditionallydescribed by assuming a monomodal soil porous
system and a single continuum approach. However, soil-porous
systems are often bimodal or multimodal with a hierarchicalcomposition of pores. A hierarchical pore composition on a
microscale was shown, for instance, by Rösslerová-Kodešová and
Kodeš (1999). Complicated soil structure was also studied on alarger scale by Císlerová et al. (2002), Císlerová and Votrubová(2002), and Votrubová et al. (2003). High variability of the soil
porous structure in such soils may cause nonequilibrium waterﬂow and contaminant transport.
Numerical models that assume bimodal soil porous systems
have been developed in the past to describe nonequilibrium waterﬂow and solute transport in such soils (Durner, 1994; Gerke andR. Kodešová, M. Kočárek, and J. Kozák, Czech Univ. of Life Sciences inPrague, Dep. of Soil Science and Geology, Kamýcká 129, 16521 Prague 6,
Czech Republic; V. Kodeš, Czech Hydrometeorological Ins tute, Dep. of
Water Quality, Na Šabatce 17, 14306 Prague 4, Czech Republic; J. Šimůnek,Univ. of California Riverside, Dep. of Environmental Sciences, 900 University Ave., A135 Bourns Hall, Riverside, CA 92521, USA. Received 23 Apr.
2007. *Corresponding author (kodesova@af.czu.cz).
Vadose Zone J. 7:798–809
doi:10.2136/vzj2007.0079© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmi ed in any form or by any means, electronic or mechanical,
including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng from the publisher.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

798

van Genuchten, 1993). The soil porous system in these models is
divided into two domains, and each domain is characterized by its
own set of transport properties and equations describing ﬂow andtransport processes. The dual-porosity approach deﬁning waterﬂow and solute transport in systems consisting of domains ofmobile and immobile water was presented by Phillip (1968) andŠimůnek et al. (2003). The dual-porosity formulation is based ona set of equations describing water ﬂow and solute transport in
the mobile domain, and mass balance equations describing moisture dynamic and solute content in the immobile domain. The
dual-permeability 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). Diﬀerent equations may be used to simulate water ﬂow
in the macropore domains (Šimůnek et al., 2003). For example,
the kinematic wave approach was used by Germann (1985),
Germann and Beven (1985), and Jarvis (1994) to describe ﬂow in
preferential ﬂow paths. Alternatively, Gerke and van Genuchten(1993, 1996) used the Richards equation to describe ﬂow inboth domains. Other approaches can be based on the Poiseuilleequation (Ahuja and Hebson, 1992) and the Green-Ampt orPhilip inﬁltration equations (Ahuja and Hebson, 1992; Chen
and Wagenet, 1992). The overview of various approaches waspreviously given by Gerke (2006). Single-porosity, dual-porosity,and dual-permeability models based on the numerical solutionof the Richards equation in all domains were implemented into
the HYDRUS-1D software package by Šimůnek et al. (2003,
2005) and successfully applied to simulate water ﬂow and solutetransport at both laboratory and ﬁeld scales by Köhne et al. (2004,
2006a,b), Pot et al. (2005), and Kodešová et al. (2005). The same
approach was also used by Vogel et al. (2000).
Soil porous systems, and subsequently, their soil hydraulicproperties, are inﬂuenced by many factors, including the mineralogical composition, stage of disintegration, organic matter,soil water content, transport processes within the soil proﬁle,weather, plant roots, soil organisms, and management practices.
Shapes and sizes of soil pores may be studied on images of thinsoil sections taken at various magniﬁcations. Pore systems withmacropores and their impact on saturated hydraulic conductivities, Ks, were previously explored by Bouma et al. (1977, 1979).
Diﬀerently shaped pores in soils under diﬀerent management
practices and their Ks values were studied by Pagliai et al. (1983,
2003, 2004). Kodešová et al. (2006) described the eﬀects ofdetectable macropores and larger capillary pores on the shape ofsoil hydraulics functions. Additionally, they used the dual-permeability model to obtain soil hydraulic properties for separate
ﬂow domains using numerical inversion. Kodešová et al. (2006)
also discussed the restricting impact of clay coatings on waterexchange between both domains.
Soil micromorphological images were used in this study toanalyze soil porous systems, to distinguish the possible charactersof water ﬂow, and to select appropriate numerical models for thethree diﬀerent soil types. Soil hydraulic properties were estimatedusing a multistep outﬂow and ponded inﬁltration experiments,
and numerical inversion. Resulting soil hydraulic properties werethen applied to simulate the water ﬂow and reactive solute transport that were studied in the ﬁeld.

where subscripts mo and im refer to the mobile and immobile
domains, respectively, θmo and θim are volumetric soil water contents in the mobile and immobile pore domains [L3 L−3], hmo is
the pressure head in the mobile domain [L], Kmo is the hydraulicconductivity [L T−1], Smo and Sim are the sink terms [T−1],and Γw
is the mass exchange between the mobile and immobile regions[T−1]. Equation [2] is solved using one set of soil water retention
and hydraulic conductivity curves deﬁned for the mobile domain
and another soil water retention curve for the immobile domain.
The mass exchange between the mobile and immobile regions is
calculated using the following equation:
Γ w = ω w (θ e,mo − θ e,im ) or Γ w = ω *w (hmo − him )where ωw [T−1] and ωw* [L−1 T−1] are ﬁrst-order rate coeﬃcientsand θe,mo and θe,im are eﬀective soil water contents in the mobile
and immobile domains, respectively. The left part of Eq. [3] isvalid when the same soil water retention curve shape parameters
are used in both domains.
In the case of the dual-permeability model, the Richards
equation is applied separately to each of the two pore regions
macropore (fractures, domain of larger pores) and matrix
domains:
⎤
∂θ f
∂ ⎡
∂h
Γ
=
⎢ K f (hf ) f + K f (hf )⎥ − S f − w⎥⎦∂t∂ z ⎢⎣
∂z
wf
⎤
∂θ m
∂ ⎡
∂h
Γ
=
⎢ K m (hm ) m + K m (hm )⎥ − S m + w⎥⎦∂t∂ z ⎢⎣∂z
1− w fSingle-Porosity, Dual-Porosity, and Dual-Permeability
Water Flow Numerical Models

Water ﬂow in the soil proﬁle may be simulated using thesingle-porosity, dual-porosity, and dual-permeability models
implemented in HYDRUS-1D (Šimůnek et al., 2005; Šimůnek
and van Genuchten, 2008). The Richards equation, describing
the one-dimensional isothermal Darcian ﬂow in a variably saturated rigid porous medium, is used in all models.
The single Richards equation is used for the single-porosity system:

β[5]
K aγ w (hf − hm )a2where Ka is the eﬀective saturated hydraulic conductivity of theinterface between the two pore domains [L T−1]. Parametersdescribing aggregate shapes are the shape factor β [–] (= 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 content retention curve, θ(h), and the hydraulicconductivity function, K(θ), are used in all models:θ ( h )− θ r
1
θe =
=
h <0
[6]
n m
θ s − θr
1+ αh
Γw =

[1]

where θ is the volumetric soil water content [L3 L−3], h is thepressure head [L], K is the hydraulic conductivity [L T−1], S is
the sink term [T−1], t is the time [T], z is the vertical axes [L].
Equation [1] is solved for the entire ﬂow domain using one set ofsoil water retention and hydraulic conductivity curves.
The dual-porosity formulation for water ﬂow can be based ona mixed formulation of the Richards equation to describe waterﬂow in mobile zones (ﬂowing water, interaggregate pores) and amass balance equation to describe soil water content dynamics inimmobile zones (stagnant water, intra-aggregates pores):
⎤
∂θ mo
∂ ⎡
∂h
=
⎢ K mo ( hmo ) mo + K mo ( hmo )⎥ − S mo −Γ w⎥⎦∂t
∂ z ⎢⎣∂z∂θ im= −S im + Γ w
θ = θ mo + θ im∂t
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

[4]

where θf and θm are volumetric soil water contents in the macropore and matrix pore domains [L3 L−3], respectively, hf and hm arepressure heads in both domains [L], Kf and Km are the hydraulicconductivities [L T−1], Sf and Sm are the sink terms [T−1], Γw
is the mass exchange between the matrix and macropore regions[T−1], and wf is the ratio between volume of the macropore
domain and the total ﬂow domain [–]. Equation [4] is solved
using two sets of soil water retention and hydraulic conductivity
curves that are deﬁned for each domain. The mass exchange between
the matrix and macropore regions is calculated using

MathemaƟcal Models

⎤
∂θ
∂ ⎡
∂h
=
⎢ K ( h ) + K ( h )⎥ − S∂ t ∂ z ⎣⎢
∂z
⎦⎥

[3]

(

θe = 1

)

h≥0

and
2
m m⎤
K ( θ ) = K s θ el ⎡⎢1−(1− θ1/
e ) ⎥
⎣
⎦
K (θ ) = K s
h ≥0

[2]

799

h<0

[7]

where θe is the eﬀective soil water content [–], Ks is the saturated hydraulic conductivity [L T−1], θr and θs are the residual
and saturated soil water contents [L3 L−3], respectively, l is thepore-connectivity parameter [–], α is reciprocal of the air entry
pressure, [L−1], n [–] is related to the slope of the retention curveat the inﬂection point, and m = 1 − 1/n [–].

s=

Φ = μ w θ c + μ sρ d s − γ w θ − γ sρ d

[8]Materials and Methods

where c [M L−3] and s [M M−1] are solute concentrations in
the liquid and solid phases, respectively, q is the volumetric ﬂux
density [L T−1], ρd is the soil bulk density [M L−3], D is the dispersion coeﬃcient [L2 T−1], and Φ describes zero- and ﬁrst-order
rate reaction [M L−3 T−1].
The dual-porosity formulation for solute transport is basedon the convection–dispersion and mass balance equations:∂θ moc mo ∂ f sρ d smo
+
∂t
∂t
∂ ⎛⎜
∂ c mo ⎞⎟ ∂ q moc mo=−Φ mo −Γ s⎟−
⎜ θ mo Dmo∂z ⎝∂ z ⎠⎟∂z

Defini on of Soil Porous Systems

The study was performed on Haplic Luvisol (parent material loess), Greyic Phaeozem (parent material loess), and Haplic
Cambisol (parent material orthogneiss). The ﬁve-year rotation
system with conventional tillage is used at all locations. The
winter barley (Hordeum vulgare L.) was planted at all areas whensoil samplings and ﬁeld experiments were performed. Six soil
horizons (Ap1 0–29 cm, Ap2 29–40 cm, Bt1 40–75 cm, Bt2
75–102 cm, BC 102–120 cm, Ck 120–145 cm) were identiﬁed in Haplic Luvisol, three horizons (Ap 0–25 cm, Bth 25–44
cm, Ck 44–125 cm) in Greyic Phaeozem, and three horizons(Ap 0–29 cm, Bw 29–62 cm, C 62–84 cm) in Haplic Cambisol.
Large undisturbed soil aggregates of an approximate size of 3 by
4 by 1.5 cm and undisturbed 100-cm3 soil samples were takenfrom each horizon of all soil proﬁles.
Micromorphological properties characterizing the soil
pore structure were studied in thin soil sections preparedfrom large soil aggregates. These thin sections were preparedaccording to methods presented by Catt (1990). Th e ﬁ nal
thin section size was 1.5 by 2 cm. The soil pore structure wasanalyzed using the procedure presented by Kodešová et al.
(2006). Images were taken at one magniﬁcation at a resolution of 300 dpi using a Nikon CoolPix 4500 camera. Th e
size of the image was 1280 × 860 pixels; the size of the pixel
side was 9.7 μm. Image-processing ﬁlters were used to detectsoil pores. True color images were ﬁrst converted to negative
images and then to 16-color images. Pixels of a particular color
range were then considered to represent possible pores using
thresholding. Resulting images were converted to black-andwhite images, compared with original images, and manuallyrevised by either deleting regions that were not pores (usually mineral grains) or adding pore areas that had diﬀerent
color using the image-processing ﬁlters. Potential disturbanceswere eliminated by smoothing the image (removing regions
smaller then 4 × 4 pixels). Black-and-white images were ﬁnallyconverted into grids using ArcGIS software (ESRI, Redlands,
CA). Individual pores were identiﬁed using a “regiongroup”
function of ArcGIS raster processing tools. Using the zonal
geometry function, the following information was obtained:pore areas, perimeters, thicknesses (the radius of the largest
circle that can be drawn within each pore), and parametersof centroids that deﬁ ne pore areas, positions, and orientations. Th e total image porosities were determined as ratios[9]∂θ imc im ∂(1− f s )ρ d s im+
= −Φ im + Γ s∂t∂twhere fs is the dimensionless fraction of sorption sites in contact
with the mobile water and Γs is the solute transfer rate between
the two regions [M L−3 T−1]. The solute transfer is describedusing the following equation:Γ s = ω s (c mo − c im )+ Γ w c *[10]where ωs [T−1] is the ﬁrst-order rate coeﬃcient and c* is equal
either to cmo for Γw > 0 or to cim for Γw < 0.
The dual-porosity formulation for solute transport is basedon two convection-dispersion equations:∂θ f c f ∂ρ d,f s f∂ ⎛⎜∂c ⎞ ∂q cΓ+=θ f Df f ⎟⎟⎟− f f −Φ f − s
⎜
wf
∂t
∂t
∂z ⎝
∂z ⎠
∂z

[11]

f

where subscripts f and m refer to the macropore and matrix
domains, respectively. The solute transfer is described as follows:
Γs =

[14]

where μw and μs are the ﬁrst-order solute degradation rate constants in water and on solid [T−1], respectively, and γw [M L−3
T−1] and γs [T−1] are the zero-order solute degradation rate constants in water and on solid, respectively.

Concepts of models for solute transport correspond to waterﬂow models described above. The single convection-dispersion equation for solute transport is used for the single-porosity system:∂θ mc m ∂ρ d,m s m
+
∂t
∂t
∂ ⎛⎜
∂ c m ⎞⎟ ∂ q mc m
Γs=−Φ m +⎟−⎜ θ m Dm1− w∂z ⎝∂ z ⎠⎟∂z[13]

where kA, β, and η are empirical coeﬃcients.
The zero- and ﬁrst-order rate reaction term, Φ, is deﬁned in
all cases as follows:

Single-Porosity, Dual-Porosity, and Dual-Permeability
Solute Transport Numerical Models∂θ c ∂ρ d s∂ ⎛⎜∂ c ⎞ ∂ qc+=
−Φ⎜ θ D ⎟⎟−∂t∂t
∂ z ⎝ ∂ z ⎠⎟ ∂ z

kAc β
1+ ηc ββDe (1− w f )θ m (c f − c m )+ Γ w c *
a2

[12]

where De is the eﬀective diﬀusion coeﬃcient [L2 T−1].
The adsorption isotherm relating in all models adsorbed concentration of solute on soil particles, s, and solute concentration,c, is described by the following general equation:www.vadosezonejournal.org · Vol. 7, No. 2, May 2008800

between total pore areas and entire image areas. Since small
size regions were eliminated, detected pores represent poreswith radii larger then approximately 20 μm.

Chlorotoluron was more regularly distributed in Haplic Luvisol
than Haplic Cambisol and was signiﬁcantly aﬀected by preferential ﬂow in Greyic Phaeozem. The maximum depths of signiﬁcant
chlorotoluron concentrations on the 35th day were observed at
6, 22, and 10 cm below the soil surface in Haplic Luvisol, Greyic
Phaeozem, and Haplic Cambisol, respectively.
Chlorotoluron transport under ﬁeld conditions was ﬁrst simulated using the single-porosity model. Since the chlorotoluron
was not detected below the depth of 80 cm, only the upper 80 cm
of the soil proﬁle was considered in numerical simulations. Eachsoil proﬁle was divided into soil layers: 0–29, 29–40, and 40–80
cm for Haplic Luvisol; 0–25, 25–44, and 44–80 cm for Greyic
Phaeozem; and 0–29, 29–62, and 62–80 cm for Haplic Cambisol.
Top boundary conditions were deﬁned using daily precipitation
rates (Fig. 1). The root water uptake was simulated assuming atime-variable root zone depth, the Feddes stress response function
model with parameters for wheat, and daily potential transpiration rates (Fig. 1). The actual root depth, LR, was simulated in
HYDRUS-1D as the product of the maximum rooting depth, Lm[L], and a root growth coeﬃcient, fr(t):Determina on of Soil Hydraulic Proper es

Soil hydraulic properties were studied in the laboratory on
undisturbed 100-cm3 soil samples (soil core height of 5.10 cm anda cross-section area of 19.60 cm2) placed in Tempe cells. First, thesaturated hydraulic conductivity was measured using the constanthead method, followed by the multistep outﬂow test. Applied pressure heads were −10, −30, −50, −80, −120, −200, −400, −650, and−1000 cm. Soil water retention data points were obtained by calculating water balance in the soil sample at each pressure head stepof the experiment. The single-porosity model in HYDRUS-1D
was ﬁrst used to analyze multistep outﬂow data and to obtain
the parameters of both soil hydraulic properties (soil water retention and unsaturated hydraulic conductivity functions) that were
described using the van Genuchten model (Eq. [6] and [7]). Pointsof the soil water retention curve were also utilized in the inversion.
While large values of the weighting factor (WF = 10) were used
for all soil water retention data points, standard values (WF = 1)
were used for the multistep outﬂow data. While the θs value was
set to the measured value, the remaining soil hydraulic parameters,
θr, α, n, and Ks, were optimized. The pore connectivity parameter
was in all cases assumed to be equal to an average value for manysoils (l = 0.5) (Mualem, 1976). In the second set of inversions,θs and Ks were set equal to the measured values, and remaining
parameters, θr, α, n, and l, were optimized. Depending on the soil
porous structure and the character of observed outﬂow and calculated water balance data, either dual-porosity or dual-permeability
models in HYDRUS-1D were applied to improve optimization
results and obtain soil hydraulic properties suitable for simulationsof solute transport monitored in the ﬁeld. Additionally, a doublering ponded inﬁltration experiment was performed on the Haplic
Cambisol site. The inner-ring diameter was 35.7 cm. The ponded
depth of 3.5 cm was applied. The dual-permeability model wasused to obtain hydraulic properties. The applied procedure is discussed in the “Results and Discussion” section.

LR ( t ) = Lm f r ( t )[15]which is deﬁned using the Verhulst–Pearl logistic growth function(Šimůnek and Suarez 1993):
f r (t ) =L0
L0 + (Lm − L0 )e −rt[16]

Field and Numerical Study of the Herbicide Transport

The transport of chlorotoluron in all three soil proﬁles wasstudied under ﬁeld conditions. The herbicide Syncuran (Synthesia,
Semtin, Czech Republic), containing 80% a.i., was applied on a
4-m2 plot on 5 May 2004 at an application rate of 2.5 kg ha−1of the active ingredient. Two liters of Syncuran solution (0.755
g L−1, e.g., 0.5 g L−1 of chlorotoluron) were applied to the soil
surface followed by one liter of fresh water. Chlorotoluron dissolves in water, adsorbs on soil particles, and degrades with time.
Soil samples from layers 2 cm thick (to the total depth of 30 cm)
were taken at three positions at each experimental plot using asampling probe 5, 13, and 35 d after the chlorotoluron application to study the residual chlorotoluron distribution in thesoil proﬁle. Chlorotoluron concentrations within the soil proﬁlewere expressed as the total amount of solute per unit mass ofthe dry soil (μg g−1). Laboratory and mathematical procedures
were described in detail by Kočárek et al. (2005), who publishedexperimental results obtained 35 d after the herbicide application. The chlorotoluron mobility in the monitored soils increased
from Haplic Luvisol to Haplic Cambisol and to Greyic Phaeozem.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 1. Precipitation (top) and potential transpiration (bottom)
rates between the first (1 April) and the last (10 June) day ofthe experiment.
801

Results and Discussion

where L0 is the initial value of the rooting depth at the beginningof the growing season [L] and, r is the growth rate [T−1]. The
growth rate r was calculated from the assumption that 50% of therooting depth was reached after 50% of the growing season. The
initial rooting depth was 20 cm in all cases. The maximum rooting depth was 120, 120, and 70 cm for Haplic Luvisol, Greyic
Phaeozem, and Haplic Cambisol, respectively, while the harvest
time was considered to be 130, 130, and 140 d after the beginning of numerical simulations for the three soils, respectively.
The evaporation was neglected since the soil surface was almost
fully covered with plants from the beginning of the experiment.
Free drainage was considered at the bottom of the soil proﬁle.
Numerical simulations were always started on April ﬁrst, to reacha realistic pressure head distribution in the soil proﬁle when herbicide was applied. Initial pressure heads in the soil proﬁle were
set to a constant value of −200 cm.
Parameters used in simulations of the reactive solute transport are presented in Table 1. Bulk densities were measured foreach horizon. Chlorotoluron adsorption isotherms were determined for the humic and subsurface horizons of each soil proﬁleusing standard laboratory procedures. Degradation rates were
estimated from the herbicide dissipations observed in the ﬁeld
assuming the ﬁrst-order rate reaction:c t = c 0e −μt

Defini on of Soil Porous Systems

Depending on the type of pedogenesis, soils exhibit diﬀerent
porous systems. Micromorphological images of humic Ap horizons of all soil types are shown in Fig. 2. The micromorphologicalimage of the humic horizon of Haplic Luvisol (Fig. 2, left) showshigher-order aggregates. The majority of detectable pores withradii larger than 20 μm and corresponding to a pressure headhigher than −70 cm are highly connected. Higher-order peds,
with intrapedal small pores separated by interpedal larger pores,
represent possible zones of immobile water. Radii of the largest
circle that could be drawn within the largest image pore were 75,
71, 280, 190, and 239 μm for the Ap1, Ap2, Bt1, Bt2, BC, andCk horizons, respectively. All pore radii were smaller than thepore radius (= 740 μm) corresponding to the pressure head of −2
cm that is usually considered to be a limit between capillary poresand macropores (Watson and Luxmoore, 1986). The image ofthe sample from the Ap2 soil horizon (not shown here) displays
a similar but more compact structure than the Ap1 horizon. Bt1,
Bt2, BC, and Ck horizons display relatively homogeneous matrix
structure with many pores smaller than 100 μm and a system oflarger pores aﬀected to varying degrees by clay coatings and inﬁllings (Bt1, Bt2, BC). Porosities detected on micromorphologicalimages were as follows: 7.75 (for the Ap1 horizon), 3.7 (Ap2),
12.2 (Bt1), 10.1 (Bt2), 13.9 (BC), and 7.3% (Ck).
The micromorphological image of the humic Ap horizon of
Greyic Phaeozem (Fig. 2, middle) shows a relatively homogeneous
matrix structure with many pores of radii smaller than 50 μm anda system of larger pores. Radii of the largest circle that could
be drawn within each pore were 92, 108, and 107 μm for theAp, Bth, and Ck horizons, respectively. All pore radii were again[17]

where ct is the actual concentration in time t, c0 is the initialconcentration, and μ is the ﬁrst-order solute degradation rate
constant [T−1]. Degradation was assumed to occur in both
liquid and solid phases (Kodešová et al., 2004). The degradation
rate constants (the same values in both phases) were evaluatedfrom the amount of pesticide that was applied on the soil
surface (25 μg cm−2) and the herbicide
1. Parameters used for reac ve solute transport simula ons.
content 35 d after its application (13,0, T−2
24.2 and 17.4 μg cm in Haplic Luvisol,
Freundlich
Freundlich
Bulkadsorp on
adsorp on
Degrada on
Greyic Phaeozem, and Haplic Cambisol,
Soil type
Horizon
Dispersivity
density
isotherm
isotherm
raterespectively). No root solute uptake wascoeﬃcient kA
coeﬃcient β
assumed. Longitudinal dispersivitiesg cm−3cmd−1
cm3β μg1−β g−1(Table 1) were set to values suggested for Haplic Luvisol
Ap1
1.66
1.5
2.89
0.80
0.019
various soil textures, experimental scales,
Ap2
1.62
2.5
2.89
0.80
0.019
and transport distances by Vanderborght
Bt1
1.51
4.0
0.86
0.84
0.019
1.54
1.5
2.73
0.75
0.001
and Vereecken (2007). The molecular dif- Greyic Phaeozem Ap
Bth
1.45
2.5
0.88
0.81
0.001
fusion was neglected. Depending on the
Ck
1.45
4.0
0.88
0.81
0.001
soil porous structure and the results of Haplic Cambisol
Ap
1.39
1.5
4.77
0.77
0.010
numerical inversions, either dual-porosBw
1.55
2.0
0.52
0.85
0.010
ity or dual-permeability models from
C
1.66
3.5
0.52
0.85
0.010
HYDRUS-1D (Šimůnek et al., 2003;
Šimůnek and van Genuchten, 2008)
were used to improve the correspondence between simulated and measured
chlorotoluron concentrations under ﬁeld
conditions. Regression coeﬃcients were
evaluated to assess agreements between
measured and simulated data. Since themeasured data represented average chlorotoluron concentrations within each
2-cm-thick layer, simulated herbicide F . 2. Micromorphological images of soil samples characterizing humic horizons of Haplic
concentrations were integrated over each Luvisol (le ), Greyic Phaeozem (middle), and Haplic Cambisol (right): A, pores; B, isolated
aggregates; C, clay coa ngs and fillings; D, grains.
layer and divided by its thickness.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008802

2. Single-porosity model parameters of soil hydraulic func ons (van Genuchten,
smaller than 740 μm. Detectable larger pores, T
1980)
for horizons of Haplic Luvisol.†corresponding to a pressure head range between−2 and −70 cm, are separated and aﬀected by
Horizon
n
l
Ks
R2
θr
θs
α
3
−3
−1
−1
clay coatings and inﬁllings that control water———— cm cm ————
cm
–
–
cm d
and solute interactions between larger pores Ap1
0.21680.3629§0.0072
1.758
0.5§
0.449
0.9985± 0.0105‡± 0.0004± 0.112± 0.956
and pores of the matrix structure. Preferential
0.2166
0.3629§0.00721.756
0.00030.180§
0.9984ﬂow through larger capillary pores may appear Ap1±0.0106± 0.0004± 0.113± 0.0209in such pore systems. Images of other soil Ap20.00010.3977§0.0221
1.161
0.5§
8.52
0.9989horizons show similar areas of structural pores,± 0.0019± 0.0014± 0.005± 16.58
0.00010.3977§0.02181.162−0.00213.34§
0.9987either aﬀected (Bth) or not aﬀected (Ck) by Ap2± 0.0102± 0.0014± 0.007± 1.626
clay coatings and inﬁllings. Image porosities
Bt1
0.27290.4125§0.0508
1.530
0.5§
34.37
0.9970were 6.2% for Ap, 5.9% for Bth, and 3.4% for± 0.0112± 0.0092± 0.132± 202.11
the Ck horizons.
Bt1
0.27280.4125§0.05081.529
0.003214.04§
0.9970Aggregates in the humic Ap horizon of± 0.0112± 0.0093± 0.132± 0.6524
Haplic Cambisol (Fig. 2, right) are poorly † θ , residual soil water content; θ , saturated soil water content; α, n, l, van Genuchten (1980)rsdeveloped. The pore system does not show
parameters; Ks, saturated hydraulic conduc vity.
intrapedal or interpedal pores, pores developed ‡ 95% confidence interval.
mainly along gravel particles. Radii of the larg- § Not op mized.
est circle that could be drawn within the largest T
3. Dual-porosity model parameters of soil hydraulic func ons for horizons ofimage pore were 753, 273, and 318 μm for the Haplic Luvisol.†
Ap, Bw, and C horizons, respectively. In thisHorizon
nmo
Ks,mo
R2
θr
θs
θ s,imα mo
ω w × 10 −4
case, several pores of the Ap horizon had radii———— cm3 cm−3 ————
cm−1–
cm d−1d−1larger than 740 μm. Macropores, correspondAp1
0.2168‡0.3629‡0.2854‡0.0072‡
1.758‡
0.213
10.69
0.9985ing to pressure heads higher than −2 cm, may± 1.090§± 90.72
create preferential pathways. Images of other Ap2
116.660.99900.0001‡ 0.3977‡ 0.3607‡ 0.0221‡ 1.161‡
1.130
soil horizons show similar features. Image± 4.601± 510.34
porosities were 32.4% for Ap, 11.8% for Bw, Bt1
0.2729‡ 0.4125‡ 0.2905‡ 0.0508‡ 1.530‡
24.26
1.36
0.9970± 162.06± 14.04
and 12.3% for the C horizon. Detected porosity of the Ap image is signiﬁcantly inﬂuenced † θ r, total residual soil water content; θ s, total saturated soil water content; θ s,im, saturated soil
water content of the immobile domain; α mo, nmo, van Genuchten (1980) parameters of theby the presence of macropores.
Determina on of Soil Hydraulic Proper es

mobile domain; Ks,mo, saturated hydraulic conduc vity of the mobile domain.
‡ Not op mized.
§ 95% confidence interval.

Soil hydraulic parameters for the Ap1, Ap2,and Bt1 soil horizons of Haplic Luvisol, shownin Table 2 (which will be used for numerical simulation of theﬁeld experiment), were obtained using the numerical inversion of
the multistep outﬂow experiment and the single-porosity model.
Measured and calibrated outﬂow data, as well as optimized andexperimental soil water retention curves (not shown here), weresatisfactorily similar. Optimized Ks values were signiﬁcantly higher
than those measured using the constant head method on the same
soil samples. The 95% conﬁdence intervals for the optimized Ks
and l values were very wide, indicating that outﬂow experiments
provide only limited information about the conductivity function.
An agreement between measured and simulated outﬂow data wasbetter when Ks rather than l was optimized.
Since the micromorphological study of the Ap1 horizon
indicated a possible occurrence of immobile zones, the dualporosity model was also applied in the inverse analysis. Soil
hydraulic parameters for the dual-porosity model obtainedusing the numerical inversion of cumulative outﬂow and soilwater retention curve data are shown in Table 3. To minimize
the number of optimized parameters in the dual-porosity model,
the residual water content of the immobile zone, θr,im, was set
equal to θr of the single-porosity model, while the residual watercontent of the mobile zone, θr,mo, was assumed to be equal tozero, the sum of the saturated water contents of the two poreregions (θs,im and θs,mo) was set equal to θs of the single-porosity model, and the shape parameters αmo (= αim) and nmo (=www.vadosezonejournal.org · Vol. 7, No. 2, May 2008nim) were also considered to be equal to α and n of the singleporosity model. Additionally, θs,im was deﬁned as the saturated
soil water content, θs, minus the measured image porosity. Only
the saturated hydraulic conductivity, Ks, and the water transfer
coeﬃcient between the immobile and mobile domains, ωw, werethus optimized. No signiﬁcant improvement in the agreement
between the measured and calibrated outﬂow data, or the optimized and experimental soil water retention curves, was observed.
Only the diﬀerence between optimized and measured Ks values
signiﬁcantly decreased. The 95% conﬁdence interval for the optimized Ks value is narrower than for the single-porosity model.
Horizons below the surface layer were analyzed similarly. While
the results for the Ap2 horizon exhibited similar improvements,
Ks does not signiﬁcantly decreased for the Bt1 horizon, and its
95% conﬁdence interval remained large. This is probably due
to diﬀerences between micromorphological features of the Bt1
horizon and the other two horizons.
Soil hydraulic parameters for the Ap, Bth, and Ck soil
horizons of Greyic Phaeozem shown in Table 4 were obtainedusing the numerical inversion of the multistep outﬂow experiment while assuming the single-porosity model. Measured andcalibrated outﬂow data shown for the Ap horizon in Fig. 3 (left)
corresponds well. However, better correlation would be desirednear saturation and below the pressure head of −400 cm. Worseagreement was again obtained between measured and simulated803

T
4. Single-porosity model parameters of the soil hydraulic func ons for horizons of Greyic Phaeozem.†

were deﬁned as soil water contents for particular pressure heads minus the saturated water content of the2matrixmultiplied by the domain fraction. The lower
Horizon
n
l
K
R
θr
θs
α
s
weightingfactor was used for data where the impact
3
−3
−1
−1
—— cm cm ——
cm
–
–
cm d
of
the
matrixsoil water retention data was expected.
Ap0.0470 0.4182§ 0.0374
1.122
0.5§
55.63
0.9987
Zero soil water contents were set for pressure heads± 0.1627‡± 0.0063 ± 0.072± 234.8Ap
0.0821 0.4182§ 0.0363
1.140
0.0097 10.26§
0.9985
below −120 cm. Considering that large pores control± 0.1218± 0.0058 ± 0.071± 1.1597
mainly saturated ﬂow, the saturated hydraulic conBth0.2760 0.4072§ 0.0281
1.626
0.5§
7.11
0.9987
ductivity of the large pore domain, Ks,f, was deﬁned± 0.0071± 0.0026 ± 0.101± 26.01as the ratio of Ks measured using the constant headBth
0.27490.4072§ 0.0283
1.617
0.00032.95†
0.9987± 0.0072± 0.0027 ± 0.101± 0.1138
experiment on the same soil sample and domain
Ck0.1520 0.4423§ 0.0272
1.348
0.5§
22.80
0.9988
fraction (wf ). The eﬀective saturated hydraulic con± 0.0242± .0020 ± 0.057± 57.58
ductivity, Ke, for the mass transfer between the matrix
Ck
0.1520 0.4423§ 0.0272
1.349
0.0158 13.20§
0.9988and macropore domains was set equal to a small
0.02400.00200.057
0.8446valueof 2.4 × 10−3 cm d−1 (10−4 cm h−1) due to the† θ r, residual soil water content; θ s, saturated soil water content; α, n, l, van Genuchtenpresence of clay coatings. Remaining parameters αm,(1980) parameters; Ks, saturated hydraulic conduc vity.
‡ 95% confidence interval.
nm, and Ks,m were optimized (Table 5) by minimiz§ Not op mized.
ing the objective function deﬁned using measured
cumulative outﬂow and retention data. Optimized
outﬂow data when the l parameter was optimized (not shown).
parameter values are shown in Table 5. The resulting total soil
The soil water retention curve ﬁt the experimental data well (Fig.
water retention curve was obtained as the sum of soil water reten3, right). Optimized Ks values were again signiﬁcantly higher thantion curves for the matrix and macropore domains multiplied bythose measured. The 95% conﬁdence intervals for the optimized
their corresponding fractions.
Ks and l values were similarly, as before, very wide.
The optimized total soil water retention curve and simulated
The dual-permeability model was then applied to improveoutﬂow data are presented in Fig. 3. Although the correlation
numerical inversion results and to obtain parameters for two ﬂow
coeﬃcient slightly decreased, Fig. 3 (left) shows that agreement
domains referred to as matrix (m) and large pore (f ) domains
between measured and optimized outﬂow data for pressure heads(the bicontinuum approach). To obtain unique optimization
close to saturation and below −400 cm exhibited improvement
results for this complex model, many of its parameters must be
when the dual-permeability model was used. The total saturated
set equal to independently estimated values. The fraction of thehydraulic conductivity (10.6 cm d−1), evaluated as the sum of Ks
large pore domain (wf = 0.15) was estimated from the fraction ofvalues for each pore domain multiplied by their corresponding
pores identiﬁed from micromorphological images. The following
domain fractions, is not substantially higher than the measured
parameters were used to deﬁne the structure of both domains: β
Ks value and is considerably lower than Ks obtained using the= 8 (the shape factor characterizing angular blocky aggregates), a
single-porosity model with l = 0.5. The 95% conﬁdence interval= 0.5 cm (the characteristic length of an aggregate deﬁning elonfor the optimized Ks,f is also narrower. Numerical inversions usinggated pathways), and γw = 0.4. Saturated water contents of thethe dual-permeability model were similarly performed for thematrix and large pore domains, θs,m and θs,f, respectively, were set
other two soil samples. The same values of αf, nf, and Ke wereassuming that the sum of their values multiplied by domains fracused for both horizons. Final optimized and ﬁxed parameters
tions (wm = 1 − wf = 0.85, wf = 0.15) was equal to the measured
are shown in Table 5. The total saturated hydraulic conductivity
θs value. Residual water contents of the matrix and large pore
(3.32 cm d−1) for the Bth horizon is again similar to the measured
domains, θr,m and θr,f, were set equal to θr of the single-porosity
Ks value and is considerably lower than the Ks obtained using themodel and zero, respectively. Parameters αf and nf (Table 5) for
single-porosity model. On the other hand, the total saturated
the large pore domain were obtained by ﬁtting the eﬀective soil
hydraulic conductivity (23.1 cm d−1) for the C horizon is simiwater retention data characterizing the large pore domain that
lar to the Ks value obtained using the single-porosity model (forl = 0.5). These diﬀerences are likely
due to diﬀerent micromorphologicalfeatures of the matrix domain of the
C horizon compared with the othertwo horizons.
Soil hydraulic parameters of thesingle-porosity model for the Ap,
Bw, and C soil horizons of Haplic
Cambisol obtained by numericalinversion of the outflow data areshown in Table 6. Calibrated outflow data and optimized soil water
retention curves satisfactorily ﬁtted
F . 3. Mul step ou low data (le ) and soil water reten on curves (right) obtained on the soil
the experimental data. Optimized Ks
sample characterizing the humic horizon of Greyic Phaeozem.
values were similar to or higher thanwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

804

T

5. Dual-permeability model parameters of the soil hydraulic func ons for horizons of Greyic Phaeozem.†

Horizon

Ap
Bth
Ck

θ r,m

θ s,m

αm

nm

Ks,m

wf

θ r,f

θ s,f

αf

nf

Ks,f

— cm3 cm−3 —

cm−1

–

cm d−1–— cm3 cm−3 —

cm−1

–

cm d−10.0038 ± 0.0005§
0.0095 ± 0.0013
0.0162 ± 0.0012

1.185 ± 0.016
1.821 ± 0.151
1.424 ± 0.021

0.3847 ± 0.2784
0.4272 ± 0.3912
10.72 ± 5.58

0.15
0.14
0.08

0‡
0‡
0‡

0.047‡
0.047‡
0.047‡

2.03‡
2.03‡
2.03‡

68.42‡
21.07‡
165, 0‡0.0553‡0.3209‡0.1854‡0.4182‡0.407‡
0.442‡

0.418‡
0.407‡
0.442‡

R2

0.9947
0.99330.9948† θ r,m, θ r,f, residual soil water content of the matrix (m) and large pore (f) domains; θ s,m, θ s,f, saturated soil water content of the matrix and large pore domains; α m, α f, nm, nf, van Genuchten (1980) parameters of the matrix and large pore domains; Ks,m, Ks,f, saturated hydraulic conduc vity of the matrix
and large pore domains; w f, ra o between volume of the large pore domain and the total flow domain.
‡ Not op mized.
§ 95% confidence interval.

T
6. Single-porosity model parameters of the soil hydraulic func ons for horizons of Haplic Cambisol op mized using the mul step outflow data.†
Horizon
Ap
Ap
Bw
Bw
C
C

θr

θs

———— cm3 cm−3 ————
0.2911 ± 0.0064‡
0.4648§0.2945 ± 0.0156
0.4648§0.2991 ± 0.0057
0.4336§
0.3060 ± 0.0043
0.4336§0.1460 ± 0.0081
0.3673§0.1503 ± 0.0080
0.3673§

α

n

l

cm−1
0.0239 ± 0.0010
0.0243 ± 0.0012
0.0309 ± 0.0016
0.0348 ± 0.0019
0.0232 ± 0.0009
0.0232 ± 0.0009–
1.636 ± 0.059
1.639 ± 0.084
1.900 ± 0.105
1.880 ± 0.093
1.456 ± 0.037
1.458 ± 0.038–
0.5†
0.0137 ± 0.1392
0.5†−2.2021 ± 0.5763
0.5
0.0003 ± 0.0239

Ks
cm d−110.50 ± 3.75
10.61§
13.48 ± 5.58
2.82§
7.49± 10.22
2.60†

R2
0.99620.99650.9936
0.99480.99910.9991† θ r, residual soil water content; θ s, saturated soil water content; α, n, l, van Genuchten (1980) parameters; Ks, saturated hydraulic conduc vity.
‡ 95% confidence interval.
§ Not op mized.

those measured. The 95% conﬁdence intervals for the Ks value areat the experimental plot. The same space discretization was usednarrower than for previous soils. Optimization of the l parameteras for the pesticide transport simulation. Initial pressure headagain did not improve agreement between measured and simuconditions were determined using the matrix soil water retention
lated data. The 95% conﬁdence intervals for the l value are again
curve and soil water contents measured in the ﬁeld. The constantrelatively wide.
pressure head of 3.5 cm (ponded depth) was used as the upper
The capillary pore systems of Haplic Cambisol are notboundary condition. Free drainage boundary condition was usedbimodal as the capillary pore systems of previous soils. The microat the bottom of the soil proﬁle. Similar to Greyic Phaeozem, themorphological study and ﬁeld observations for the Ap horizon
eﬀective saturated hydraulic conductivity, Ke, was set equal to 2.4
indicated the presence of gravitational macropores (with pore× 10−3 cm d−1. The resulting parameters are shown in Table 7.
radii larger than 740 μm). Since the multistep outﬂow experiment
Field and Numerical Study of the Herbicide Transport
is not particularly suitable for characterization of macropore ﬂow
since almost the entire experiment is performed under unsatuThe chlorotoluron transport in Haplic Luvisol was simulated
rated conditions, and since an additional ponded inﬁltration
using the single-porosity and dual-porosity models. In addition
experiment was performed at this ﬁeld site, the soil hydraulicto input data described above, the solute transfer coeﬃcient, ωs,parameters for the dual-permeability model were determined in
between the mobile and immobile domains was set equal to
the following way. The soil hydraulic parameters obtained by
10−3 d−1 (a similar ωs value was obtained for the Ap horizon of
the optimization process of the multistep outﬂow experiment
sandy loam by Köhne et al., 2006b). Cumulative surface waterusing the single-porosity model were assumed to characterize only
ﬂuxes, root water uptakes and bottom water ﬂuxes simulated
the matrix pore domain. The additional macropore domain with
between the herbicide application (5 May) and the end of thepores corresponding to pressure heads larger than −2 cm was
experiment (10 June) using both models are shown in Fig. 4.
represented using a steplike shape retention curve with paramCumulative surface ﬂuxes are the same since no evaporation waseters αf = 0.07 cm−1 and nf = 3, ensuring that macropores wereconsidered and all precipitated water inﬁltrated into the soil proﬁlled with water only for pressure heads close to zero. The
macropore domain fraction (wf = 0.05) for all horizons was T
7. Single-porosity model parameters of the soil hydraulic func ons forset based on ﬁeld observations. The saturated water con- horizons of Haplic Cambisol evaluated using the infiltra on experiment data.†
tents of the macropore domain, θs,f, was set assuming that
Horizon
n
l
Ks
R2
θr
θs
α
the sum of θs,m and θs,f multiplied by their corresponding——cm3 cm−3 ——cm−1––
cm d−1domain fractions (wm = 0.95, wf = 0.05) was equal to theAp0‡0.4740‡0.07‡
3‡0.5‡ 1911.2 ± 195§
0.994833
average measured porosities of 0.489, 0.456, and 0.386 cm
Bw
0‡0.4520‡0.07‡
3‡0.5‡ 205.5 ± 65.6
−3
cm for the Ap, Bw, and C soil horizons. The following C
0‡0.3780‡0.07‡
3‡0.5‡ 86.4 ± 7.9
parameters deﬁning the domain’s structure were used: β = 8,
† θ r, residual soil water content; θ s, saturated soil water content; α, n, l, van Genua = 0.2 cm, and γw = 0.4. The Ks values for the macropore
chten (1980) parameters; Ks, saturated hydraulic conduc vity.
domain of all horizons were obtained by numerical inversion ‡ Not op mized.
of the ponded inﬁltration experiment that was performed § 95% confidence interval.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008805

F . 4. Cumula ve surface fluxes(CSF), cumula ve root water
uptakes (CRWU), and cumula vebo om fluxes (CBF) simulated
between the herbicide applicaon (5 May) and the end of theexperiment (10 June) using thesingle-porosity (SPM) and dualporosity (DPM) models in Haplic
Luvisol.

F . 5. Distribu on ofchlorotoluron in Haplic
Luvisol measured (le )and simulated using
single-porosity (middle)
and dual-porosity (right)
models.

ﬁle. The simulated ﬁnal root depth was 75 cm in both cases. The
were assumed to be the same as in the matrix domain. Although
cumulative root water uptake simulated using the dual-porosity
the soil hydraulic characteristics of the single-porosity model weremodel is lower than that simulated using the single-porosity
very similar to the total soil hydraulic characteristics of the dualmodel because of the presence of immobile domains and lower
permeability model (they were obtained using numerical inversion
hydraulic conductivity, resulting in the higher retention abilityof the same multistep outﬂow and soil water retention data), theof the dual-porosity system. This phenomenon is also reﬂected
simulated ﬁeld water ﬂow and solute transport were signiﬁcantly
in simulated cumulative bottom outﬂows. Observed concentradiﬀerent. Cumulative surface water ﬂuxes, root water uptakes,
tions in the soil proﬁles were expressed as total amounts of soluteand bottom water ﬂuxes between the herbicide application and(present in soil water and adsorbed on soil particles) per mass
the end of the experiment simulated using single-porosity andunit of the dry soil and compared to simulated chlorotoluron
dual-permeability models are shown in Fig. 6. Cumulative surconcentrations. Measured and calculated chlorotoluron concenface ﬂuxes are again the same. The cumulative root water uptake
trations 5, 13, and 35 d after the application are shown in Fig.
simulated using the dual-permeability model is higher than that
5. Regression coeﬃcients describing the correlation between thesimulated using the single-porosity model, because of the fastermeasured and simulated concentrations are presented in Table 8.
saturation of the root zone through preferential pathways, andChlorotoluron was regularly distributed in the soil. Although theT
8. Regression coeﬃcients describing correla on between
simulated data using the single-porosity model (Fig. 5, middle)
measured and simulated concentra ons for Haplic Luvisol.
approximated experimental results (Fig. 5, left) reasonably well,
Chlorotoluron distribu on in the soil profile Day 5
Day 13 Day 35
the dual-porosity model (Fig. 5, right) provided a slightly better
R2 measured vs. simulated using the single- 0.9964 0.9992 0.9708
description of the chlorotoluron distribution in the soil proﬁle.
porosity model
The slightly lower herbicide leaching simulated using the dualR2 measured vs. simulated using the dual0.9964 0.9992 0.9743
porosity model than that simulated using the single-porosity
porosity model
model was caused by the presence of immobile domains andlower hydraulic conductivity.
The chlorotoluron transport in Greyic Phaeozem wassimulated using the singleporosity and dual-permeability
models. Additional solute transport parameters required by
the dual-permeability model
included longitudinal dispersivities in macropores, whichwere set to 3, 4, and 5 cm for
the Ap, Bth, and Ck horizons,
respectively, and sorption prop- F . 6. Cumula ve surface fluxes (CSF), cumula ve root water uptakes (CRWU), and cumula ve bo om
erties and degradation rates of fluxes (CBF) simulated between the herbicide applica on (5 May) and the end of the experiment (10the large pore domain, which June) using the single-porosity (SPM) and dual-permeability (DPM) models in Greyic Phaeozem.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008806

F . 7. Distribu on ofchlorotoluron in Greyic
Phaeozem measured(le ) and simulatedusing single-porosity(middle), and dualpermeability (right)
models.

correspondingly no or low water stress reduction simulated using
the signiﬁcantly higher cumulative bottom outﬂow simulated
the dual-permeability model. The simulated ﬁnal root depth wasusing the dual-permeability model than that simulated using
in both cases 75 cm. The dual-permeability model predicted
the single-porosity model. The dual-permeability model simuhigher cumulative bottom outﬂow than the single-porosity model
lated almost no preferential ﬂow at the bottom of the soil proﬁle.
due to again faster saturation of the entire soil proﬁle. The dualThe simulated cumulative outﬂow from the large pore domain
permeability model simulated very low preferential ﬂow at thewas four orders of magnitude lower than the simulated outﬂow
bottom of the soil proﬁle. The simulated cumulative outﬂow fromfrom the matrix domain. Measured and calculated chlorotoluron
the large pore domain was 2 orders of magnitude lower thanconcentrations 5, 13, and 35 d after the application in Haplic
the simulated outﬂow from the matrix domain. Measured andCambisol are shown in Fig. 9. Regression coeﬃcients describing
calculated chlorotoluron concentrations 5, 13, and 35 d after thethe correlation between measured and simulated concentrations
application in Greyic Phaeozem are shown in Fig. 7. Even thoughare presented Table 10. The numerical simulation that used theregression coeﬃcients (Table 9) for the 5th and 35th days are
single-porosity model (Fig. 9, middle) again underestimatedslightly lower, a comparison of plotted simulated and measured
observed herbicide mobility (Fig. 9, left). On the other hand,chlorotoluron concentrations reveals that the dual-permeability
the dual-permeability model (Fig. 9, right) predicted herbicide
model (Fig. 7, right) described the chlorotoluron transport (Fig. 7,
movement down to the depth of 40 cm in the macropore domainleft) better than the single-porosity model (Fig. 7, middle), mainly(not obvious in presented ﬁgures due to the very low soil water
because the observed herbicide transport was strongly aﬀected bycontent and low wf) and consequently large herbicide leaching
preferential ﬂow. While solute moved only to a depth of 8 cm inin the entire ﬂow region.
the single-porosity system, in the dual-permeability system solute
The dual-permeability models in both cases (for Greyic
moved to a depth of 26 cm, compared with an observed depth
Phaeozem and Haplic Cambisol) simulated insigniﬁcant bottomof 22 cm. Observed oscillations in chlorotoluron concentrations
drainage through the preferential pathways. However, their presbetween 6 and 22 cm depths were probably caused by preferenence caused fast matrix saturation in the entire soil proﬁle andtial ﬂow, that is, a fast solute penetration to greater depths, andconsequently, the larger cumulative matrix outﬂow at the bottom
a solute accumulation at the end of preferential pathways (duethan that simulated using the single-porosity model. The larger
to their disconnection). This phenomenon was not considered
chlorotoluron leaching was caused mainly by herbicide transport
in the numerical model.
The chlorotoluron transport in Haplic Cambisol was again
T
9. Regression coeﬃcients describing correla on between
simulated using the single-porosity and dual-permeability models.
measured and simulated concentra ons for Greyic Phaeozem.
Similar to Greyic Phaeozem, sorption properties and degradation
Chlorotoluron distribu on in the soil profile Day 5
Day 13 Day 35
rates for the large pore domain were assumed to be the same as
R2 measured vs. simulated using the single- 0.9973 0.9886 0.8575
those in the matrix domain, and the longitudinal dispersivity
porosity modelin macropores was set to 3, 4, and 5 cm for the Ap, Bw, and
0.9897 0.9902 0.8548
R2 measured vs. simulated using the dualC horizons, respectively. As expected, the additional macropore
permeability model
domain considerably aﬀected simulated water ﬂow and solutetransport. Cumulative water
fluxes simulated between theherbicide application and theend of the experiment using
the single-porosity and dualpermeability models are shownin Fig. 8. Surface cumulative
ﬂuxes are again the same. The
simulated final root depth
was in both cases 45 cm. The
cumulative root water uptake
simulated using the dual-permeability model is slightly lower F . 8. Cumula ve surface fluxes (CSF), cumula ve root water uptakes (CRWU), and cumula ve bo om
than that simulated using the fluxes (CBF) simulated between the herbicide applica on (5 May) and the end of the experiment (10single-porosity model due to June) using the single-porosity (SPM) and dual-permeability (DPM) models in Haplic Cambisol.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008807

F . 9. Distribu on ofchlorotoluron in Haplic
Cambisol measured (le )and simulated using
single-porosity (middle)
and dual-permeability(right) models.

References

T
10. Regression coeﬃcients describing correla on between
measured and simulated concentra ons for Haplic Cambisol.
Chlorotoluron distribu on in the soil profile
R2 measured vs. simulated using the singleporosity modelR2 measured vs. simulated using the dualpermeability modelAhuja, L.R., and C. Hebson. 1992. Root zone water quality model. GPSR Technical Rep. No. 2. USDA–ARS, Fort Collins, CO.
Bouma, J., A. Jongerius, O. Boersma, A. Jager, and D. Schoonderbeek. 1977.
The function of diﬀerent types of macropores during saturated ﬂow
through four swelling soil horizons. Soil Sci. Soc. Am. J. 41:945–950.
Bouma, J., A. Jongerius, and D. Schoonderbeek. 1979. Calculation of saturated
hydraulic conductivity of some pedal clay soils using micromorphometricdata. Soil Sci. Soc. Am. J. 43:261–264.
Catt, J.A. 1990. Paleopedology manual. Quaternary Int. 6:1–95.
Chen, C., and R.J. Wagenet. 1992. Simulation of water and chemicals in macropore soils: Part 1. Representation of the equivalent macropore inﬂuenceand its eﬀect on soil water ﬂow. J. Hydrol. 130:105–126.
Císlerová, M., T. Vogel, J. Votrubová, and A. Robovská. 2002. Searching below thresholds: Tracing the origins of preferential ﬂow within undisturbed
soil samples. p. 265–274. In Environmental mechanics: Water, mass andenergy transfer in the biosphere. Geophysical Monogr. 129. AGU, Washington, DC.
Císlerová, M., and J. Votrubová. 2002. CT derived porosity distribution and
ﬂow domains. J. Hydrol. 267:186–200.
Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30:211–233.
Gerke, H.H. 2006. Preferential ﬂow descriptions for structured soils. J. Plant
Nutr. Soil Sci. 169:382–400.
Gerke, H.H., and M.Th. van Genuchten. 1993. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous
media. Water Resour. Res. 29:305–319.
Gerke, H.H., and M.Th. van Genuchten. 1996. Macroscopic representation of
structural geometry for simulating water and solute movement in dual-porosity media. Adv. Water Resour. 19:343–357.
Germann, P.F. 1985. Kinematic wave approach to inﬁltration and drainage intoand from soil macropores. Trans. ASAE 28:745–749.
Germann, P.F., and K. Beven. 1985. Kinematic wave approximation to inﬁltration into soils with sorbing macropores. Water Resour. Res. 21:990–996.
Jarvis, J. 1994. The MACRO model: Technical description and sample simulation. Reports and Dissertations 19. Dep. of Soil Science, Swedish Univ. ofAgricultural Science, Uppsala.
Kočárek, M., R. Kodešová, J. Kozák, O. Drábek, and O. Vacek. 2005. Chlorotoluron behaviour in ﬁve diﬀerent soil types. Plant Soil Environ. 51:304–
309.
Kodešová, R., V. Kodeš. A. Žigová, and J. Šimůnek. 2006. Impact of plant roots
and soil organisms on soil micromorphology and hydraulic properties.
Biologia Bratislava 61 (Suppl. 19):S339–S343.
Kodešová, R., J. Kozák, J. Šimůnek, and O. Vacek. 2005. Single and dual-permeability model of chlorotoluron transport in the soil proﬁle. Plant Soil
Environ. 51:310–315.
Kodešová, R., J. Kozák, and O. Vacek. 2004. Field and numerical study of chlorotoluron transport in the soil proﬁle. Plant Soil Environ. 50:333–338.
Köhne, J.M., S. Köhne, B.P. Mohanty, and J. Šimůnek. 2004. Inverse mobile–immobile modeling of transport during transient ﬂow: Eﬀects of betweendomain transfer and initial water content. Vadose Zone J. 3:1309–1321.
Köhne, J.M., S. Köhne, and J. Šimůnek. 2006a. Multi-process herbicide transport in structured soil columns: Experiments and model analysis. J. Contam. Hydrol. 85:1–32.
Köhne, J.M., B.P. Mohanty, and J. Šimůnek. 2006b. 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.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of

Day 5
Day 13 Day 35
0.8142 0.8732 0.7791
0.98490.99590.8212through the preferential pathways, from which herbicide penetrated into the matrix.

Conclusions
Collected ﬁeld data and numerical simulations demonstrated
that the multiporous nature of soils has a varying impact on waterﬂow and transport processes at the ﬁeld scale. The chlorotoluron
mobility in the monitored soils increased from Haplic Luvisol to
Haplic Cambisol and to Greyic Phaeozem. The pesticide mobility
reﬂects the soil porous system and atmospheric boundary conditions. Chlorotoluron was regularly distributed in the highly
connected domain of larger pores of Haplic Luvisol, from which
it penetrated into the soil aggregates, that is, zones of immobile
water. The lowest daily precipitation rates were recorded at thissite compared to the other two locations. The highest mobility ofchlorotoluron in Greyic Phaeozem was caused by larger capillarypore pathways and suﬃcient inﬁltration ﬂuxes that occasionally
ﬁlled up these pores. The presence of clay coatings in Greyic
Phaeozem that restrict water ﬂow and contaminant transport
between the macropore and matrix domains is an additionalcause for this preferential transport that produces chlorotoluron
penetration into deeper depths. Chlorotoluron was less regularly
distributed in Haplic Cambisol. Despite the highest inﬁltration
rate, preferential ﬂow only slightly aﬀected the herbicide transport. Large gravitational pores that may dominate water ﬂow andsolute transport under saturated conditions were inactive during
the monitored period. As a result of complex interactions between
meteorological conditions and the soil pore structure, the singleand dual-porosity models describe the herbicide behavior in
Haplic Luvisol well, while the dual-permeability model performs
better in simulating the herbicide transport in Greyic Phaeozem
and Haplic Cambisol.
A
This work has been supported by the grants No. GA CR
103/05/2143 and MSM 6046070901. The authors acknowledge A.
Žigová, K. Němeček, and O. Drábek for helping with the ﬁeld andlaboratory work, and anonymous reviewers for their valuable remarks
and suggestions.

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

808

unsaturated porous media. Water Resour. Res. 12:513–522.
Pagliai, M., M. La Mace, and G. Lucamante. 1983. Micromorphometric andmicromorphological investigations of a clay loam soil in viticulture underzero and conventional tillage. J. Soil Sci. 34:391–403.
Pagliai, M., A. Marsili, P. Servadio, N. Vignozzi, and S. Pellegrini. 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., N. Vignozzi, and S. Pellegrini. 2004. Soil structure and the eﬀect ofmanagement practices. Soil Tillage Res. 79:131–143.
Phillip, J.R. 1968. The theory of adsorption in aggregated media. Aust. J. Soil
Res. 6:1–19.
Pot, V., J. Šimůnek, P. Benoit, Y. Coquet, A. Yra, and M.-J. Martinez Cordun.
2005. Impact of rainfall intensity on the transport of two herbicides inundisturbed grassed ﬁlter strip soil cores. J. Contam. Hydrol. 81:63–88.
Rösslerová-Kodešová, R., and V. Kodeš. 1999. Percolation model for interpretation of moisture retention curves for mono-modal and bi-modal soil porous systems. p. 81–91. In M.Th. van Genuchten, F.J. Leij, and L. Wu (ed.)
Characterization and measurement of the hydraulic properties of unsaturated porous media. Univ. of California, Riverside.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review
and comparison of models for describing non-equilibrium and preferential
ﬂow and transport in the vadose zone. J. Hydrol. 272:14–35.
Šimůnek, J., and D.L. Suarez. 1993. Modelling of carbon dioxide transportand production in soil: 1. Model development. Water Resour. Res.
29:487–497.
Šimůnek, J., and M.Th. van Genuchten. 2008. Modeling nonequilibrium ﬂow
and transport processes using HYDRUS. Vadose Zone J. 7:782–797(this issue).
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2005. The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media. Version 3.0. Dep.
of Environmental Sciences, Univ. of California Riverside.
Vanderborght, J., and H. Vereecken. 2007. Review of dispersivities for transportmodeling in soils. Vadose Zone J. 6:29–52.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
Vogel, T., H.H. Gerke, R. Zhang, and M.Th. van Genuchten. 2000. Modelingﬂow and transport in a two-dimensional dual-permeability system withspatial variable hydraulic properties. J. Hydrol. 238:78–89.
Votrubová, J., M. Císlerová, M.H.G. Amin, and L.D. Hall. 2003. Recurrent
ponded inﬁltration into structured soil: A magnetic resonance imaging
study. Water Resour. Res. 39(12):1371, doi:10.1029/2003WR002222.
Watson, K.W., and R.J. Luxmoore. 1986. Estimating macroporosity in aforest watershed by use of a tension inﬁltrometer. Soil Sci. Soc. Am. J.
50(3):578–582.

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008809