VZJ
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 27 May 2008
Published in Vadose Zone J 7:798-809 (2008)
DOI: 10.2136/vzj2007.0079
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Kodesová, R.
Right arrow Articles by Kozák, J.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Kodesová, R.
Right arrow Articles by Kozák, J.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Kodesová, R.
Right arrow Articles by Kozák, J.
Related Collections
Right arrow Solute Transport Models
Right arrow Water Flow Models

SPECIAL SECTION: VADOSE ZONE MODELING

Impact of Soil Micromorphological Features on Water Flow and Herbicide Transport in Soils

Radka Kodesováa,*, Martin Kocáreka, Vít Kodesb, Jirí Simunekc and Josef Kozáka

a Czech Univ. of Life Sciences in Prague, Dep. of Soil Science and Geology, Kamycká 129, 16521 Prague 6, Czech Republic
b Czech Hydrometeorological Institute, Dep. of Water Quality, Na Sabatce 17, 14306 Prague 4, Czech Republic
c Univ. of California Riverside, Dep. of Environmental Sciences, 900 University Ave., A135 Bourns Hall, Riverside, CA 92521, USA
* Corresponding author (kodesova{at}af.czu.cz).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.


Received 23 April 2007.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 Mathematical Models
 Materials and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
The impact of varying soil micromorphology on soil hydraulic properties and, consequently, on water flow and herbicide transport observed in the field is demonstrated on three soil types. The micromorphological image of a humic horizon of 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, separating higher-order peds with small intrapores that possibly formed zones with immobile water. Herbicide was regularly distributed in this soil. The majority of detectable large capillary pores in a humic horizon of Greyic Phaeozem were separated and affected by clay coatings and fillings. The herbicide transport in this soil was highly affected by preferential flow. Macropores corresponding to pressure heads higher than –2 cm were detected in a humic horizon of Haplic Cambisol. However, preferential flow only slightly influenced the herbicide transport in this soil. Single-porosity and either dual-porosity or dual-permeability flow and transport models in HYDRUS-1D were used to estimate the soil hydraulic parameters from laboratory multistep outflow and ponded infiltration experiments via numerical inversion and to simulate the herbicide transport experimentally studied in the field. Appropriate models were selected on the basis of the soil micromorphological study.

Abbreviations: WF, weighting factor


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
Mathematical Models
 Materials and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
WATER FLOW and solute transport are traditionally described by assuming a monomodal soil porous system and a single continuum approach. However, soil-porous systems are often bimodal or multimodal with a hierarchical composition of pores. A hierarchical pore composition on a microscale was shown, for instance, by Rösslerová-Kodesová and Kodes (1999). Complicated soil structure was also studied on a larger 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 flow and contaminant transport.

Numerical models that assume bimodal soil porous systems have been developed in the past to describe nonequilibrium water flow and solute transport in such soils (Durner, 1994; Gerke and 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 flow and transport processes. The dual-porosity approach defining water flow and solute transport in systems consisting of domains of mobile and immobile water was presented by Phillip (1968) and Simunek et al. (2003). The dual-porosity formulation is based on a set of equations describing water flow 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 flow and solute transport occur in both domains. The dual-permeability formulation is based on a set of equations that describe water flow and solute transport separately in each domain (matrix and macropore domains). Different equations may be used to simulate water flow in the macropore domains (Simunek et al., 2003). For example, the kinematic wave approach was used by Germann (1985), Germann and Beven (1985), and Jarvis (1994) to describe flow in preferential flow paths. Alternatively, Gerke and van Genuchten (1993, 1996) used the Richards equation to describe flow in both domains. Other approaches can be based on the Poiseuille equation (Ahuja and Hebson, 1992) and the Green-Ampt or Philip infiltration equations (Ahuja and Hebson, 1992; Chen and Wagenet, 1992). The overview of various approaches was previously given by Gerke (2006). Single-porosity, dual-porosity, and dual-permeability models based on the numerical solution of the Richards equation in all domains were implemented into the HYDRUS-1D software package by Simunek et al. (2003, 2005) and successfully applied to simulate water flow and solute transport at both laboratory and field scales by Köhne et al. (2004, 2006a,b), Pot et al. (2005), and Kodesová et al. (2005). The same approach was also used by Vogel et al. (2000).

Soil porous systems, and subsequently, their soil hydraulic properties, are influenced by many factors, including the mineralogical composition, stage of disintegration, organic matter, soil water content, transport processes within the soil profile, weather, plant roots, soil organisms, and management practices. Shapes and sizes of soil pores may be studied on images of thin soil sections taken at various magnifications. Pore systems with macropores and their impact on saturated hydraulic conductivities, Ks, were previously explored by Bouma et al. (1977, 1979). Differently shaped pores in soils under different management practices and their Ks values were studied by Pagliai et al. (1983, 2003, 2004). Kodesová et al. (2006) described the effects of detectable macropores and larger capillary pores on the shape of soil hydraulics functions. Additionally, they used the dual-permeability model to obtain soil hydraulic properties for separate flow domains using numerical inversion. Kodesová et al. (2006) also discussed the restricting impact of clay coatings on water exchange between both domains.

Soil micromorphological images were used in this study to analyze soil porous systems, to distinguish the possible characters of water flow, and to select appropriate numerical models for the three different soil types. Soil hydraulic properties were estimated using a multistep outflow and ponded infiltration experiments, and numerical inversion. Resulting soil hydraulic properties were then applied to simulate the water flow and reactive solute transport that were studied in the field.


   Mathematical Models
TOP
 ABSTRACT
 INTRODUCTION
 Mathematical Models
Materials and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
Single-Porosity, Dual-Porosity, and Dual-Permeability Water Flow Numerical Models
Water flow in the soil profile may be simulated using the single-porosity, dual-porosity, and dual-permeability models implemented in HYDRUS-1D (Simunek et al., 2005; Simunek and van Genuchten, 2008). The Richards equation, describing the one-dimensional isothermal Darcian flow in a variably saturated rigid porous medium, is used in all models.

The single Richards equation is used for the single-porosity system:

Formula 1[1]
where {theta} is the volumetric soil water content [L3 L–3], h is the pressure 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 flow domain using one set of soil water retention and hydraulic conductivity curves.

The dual-porosity formulation for water flow can be based on a mixed formulation of the Richards equation to describe water flow in mobile zones (flowing water, interaggregate pores) and a mass balance equation to describe soil water content dynamics in immobile zones (stagnant water, intra-aggregates pores):

Formula 2[2]
where subscripts mo and im refer to the mobile and immobile domains, respectively, {theta}mo and {theta}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 hydraulic conductivity [L T–1], Smo and Sim are the sink terms [T–1],and {Gamma}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 defined 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:

Formula 3[3]
where {omega}w [T–1] and {omega}w* [L–1 T–1] are first-order rate coefficients and {theta}e,mo and {theta}e,im are effective soil water contents in the mobile and immobile domains, respectively. The left part of Eq. [3] is valid 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:

Formula 4[4]
where {theta}f and {theta}m are volumetric soil water contents in the macropore and matrix pore domains [L3 L–3], respectively, hf and hm are pressure heads in both domains [L], Kf and Km are the hydraulic conductivities [L T–1], Sf and Sm are the sink terms [T–1], {Gamma}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 flow domain [–]. Equation [4] is solved using two sets of soil water retention and hydraulic conductivity curves that are defined for each domain. The mass exchange between the matrix and macropore regions is calculated using

Formula 5[5]
where Ka is the effective saturated hydraulic conductivity of the interface between the two pore domains [L T–1]. Parameters describing aggregate shapes are the shape factor β [–] (= 15 for spherical 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 {gamma}w [–] (= 0.4).

Analytical expressions proposed by van Genuchten (1980) for the soil water content retention curve, {theta}(h), and the hydraulic conductivity function, K({theta}), are used in all models:

Formula 6[6]
and

Formula 7[7]
where {theta}e is the effective soil water content [–], Ks is the saturated hydraulic conductivity [L T–1], {theta}r and {theta}s are the residual and saturated soil water contents [L3 L–3], respectively, l is the pore-connectivity parameter [–], {alpha} is reciprocal of the air entry pressure, [L–1], n [–] is related to the slope of the retention curve at the inflection point, and m = 1 – 1/n [–].

Single-Porosity, Dual-Porosity, and Dual-Permeability Solute Transport Numerical Models
Concepts of models for solute transport correspond to water flow models described above. The single convection-dispersion equation for solute transport is used for the single-porosity system:

Formula 8[8]
where c [M L–3] and s [M M–1] are solute concentrations in the liquid and solid phases, respectively, q is the volumetric flux density [L T–1], {rho}d is the soil bulk density [M L–3], D is the dispersion coefficient [L2 T–1], and {Phi} describes zero- and first-order rate reaction [M L–3 T–1].

The dual-porosity formulation for solute transport is based on the convection–dispersion and mass balance equations:

Formula 9[9]
where fs is the dimensionless fraction of sorption sites in contact with the mobile water and {Gamma}s is the solute transfer rate between the two regions [M L–3 T–1]. The solute transfer is described using the following equation:

Formula 10[10]
where {omega}s [T–1] is the first-order rate coefficient and c* is equal either to cmo for {Gamma}w > 0 or to cim for {Gamma}w < 0.

The dual-porosity formulation for solute transport is based on two convection-dispersion equations:

Formula 11[11]
where subscripts f and m refer to the macropore and matrix domains, respectively. The solute transfer is described as follows:

Formula 12[12]
where De is the effective diffusion coefficient [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:

Formula 13[13]
where kA, β, and {eta} are empirical coefficients.

The zero- and first-order rate reaction term, {Phi}, is defined in all cases as follows:

Formula 14[14]
where µw and µs are the first-order solute degradation rate constants in water and on solid [T–1], respectively, and {gamma}w [M L–3 T–1] and {gamma}s [T–1] are the zero-order solute degradation rate constants in water and on solid, respectively.


   Materials and Methods
TOP
 ABSTRACT
 INTRODUCTION
 Mathematical Models
 Materials and Methods
Results and Discussion
 Conclusions
 REFERENCES
 
Definition 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 five-year rotation system with conventional tillage is used at all locations. The winter barley (Hordeum vulgare L.) was planted at all areas when soil samplings and field 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 identified 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 taken from each horizon of all soil profiles.

Micromorphological properties characterizing the soil pore structure were studied in thin soil sections prepared from large soil aggregates. These thin sections were prepared according to methods presented by Catt (1990). The final thin section size was 1.5 by 2 cm. The soil pore structure was analyzed using the procedure presented by Kodesová et al. (2006). Images were taken at one magnification at a resolution of 300 dpi using a Nikon CoolPix 4500 camera. The size of the image was 1280 x 860 pixels; the size of the pixel side was 9.7 µm. Image-processing filters were used to detect soil pores. True color images were first 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-and-white images, compared with original images, and manually revised by either deleting regions that were not pores (usually mineral grains) or adding pore areas that had different color using the image-processing filters. Potential disturbances were eliminated by smoothing the image (removing regions smaller then 4 x 4 pixels). Black-and-white images were finally converted into grids using ArcGIS software (ESRI, Redlands, CA). Individual pores were identified 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 parameters of centroids that define pore areas, positions, and orientations. The total image porosities were determined as ratios between total pore areas and entire image areas. Since small size regions were eliminated, detected pores represent pores with radii larger then approximately 20 µm.

Determination of Soil Hydraulic Properties
Soil hydraulic properties were studied in the laboratory on undisturbed 100-cm3 soil samples (soil core height of 5.10 cm and a cross-section area of 19.60 cm2) placed in Tempe cells. First, the saturated hydraulic conductivity was measured using the constant head method, followed by the multistep outflow 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 step of the experiment. The single-porosity model in HYDRUS-1D was first used to analyze multistep outflow 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]). Points of 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 outflow data. While the {theta}s value was set to the measured value, the remaining soil hydraulic parameters, {theta}r, {alpha}, n, and Ks, were optimized. The pore connectivity parameter was in all cases assumed to be equal to an average value for many soils (l = 0.5) (Mualem, 1976). In the second set of inversions, {theta}s and Ks were set equal to the measured values, and remaining parameters, {theta}r, {alpha}, n, and l, were optimized. Depending on the soil porous structure and the character of observed outflow 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 simulations of solute transport monitored in the field. Additionally, a double-ring ponded infiltration 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 was used to obtain hydraulic properties. The applied procedure is discussed in the "Results and Discussion" section.

Field and Numerical Study of the Herbicide Transport
The transport of chlorotoluron in all three soil profiles was studied under field 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–1 of 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 a sampling probe 5, 13, and 35 d after the chlorotoluron application to study the residual chlorotoluron distribution in the soil profile. Chlorotoluron concentrations within the soil profile were expressed as the total amount of solute per unit mass of the dry soil (µg g–1). Laboratory and mathematical procedures were described in detail by Kocárek et al. (2005), who published experimental 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. Chlorotoluron was more regularly distributed in Haplic Luvisol than Haplic Cambisol and was significantly affected by preferential flow in Greyic Phaeozem. The maximum depths of significant 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 field conditions was first 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 profile was considered in numerical simulations. Each soil profile 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 defined using daily precipitation rates (Fig. 1 ). The root water uptake was simulated assuming a time-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 coefficient, fr(t):

Formula 15[15]
which is defined using the Verhulst–Pearl logistic growth function (Simunek and Suarez 1993):

Formula 16[16]
where L0 is the initial value of the rooting depth at the beginning of the growing season [L] and, r is the growth rate [T–1]. The growth rate r was calculated from the assumption that 50% of the rooting 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 profile. Numerical simulations were always started on April first, to reach a realistic pressure head distribution in the soil profile when herbicide was applied. Initial pressure heads in the soil profile were set to a constant value of –200 cm.


Figure 1
View larger version (26K):
[in this window]
[in a new window]

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

 
Parameters used in simulations of the reactive solute transport are presented in Table 1 . Bulk densities were measured for each horizon. Chlorotoluron adsorption isotherms were determined for the humic and subsurface horizons of each soil profile using standard laboratory procedures. Degradation rates were estimated from the herbicide dissipations observed in the field assuming the first-order rate reaction:

Formula 17[17]
where ct is the actual concentration in time t, c0 is the initial concentration, and µ is the first-order solute degradation rate constant [T–1]. Degradation was assumed to occur in both liquid and solid phases (Kodesová et al., 2004). The degradation rate constants (the same values in both phases) were evaluated from the amount of pesticide that was applied on the soil surface (25 µg cm–2) and the herbicide content 35 d after its application (13,0, 24.2 and 17.4 µg cm–2 in Haplic Luvisol, Greyic Phaeozem, and Haplic Cambisol, respectively). No root solute uptake was assumed. Longitudinal dispersivities (Table 1) were set to values suggested for various soil textures, experimental scales, and transport distances by Vanderborght and Vereecken (2007). The molecular diffusion was neglected. Depending on the soil porous structure and the results of numerical inversions, either dual-porosity or dual-permeability models from HYDRUS-1D (Simunek et al., 2003; Simunek and van Genuchten, 2008) were used to improve the correspondence between simulated and measured chlorotoluron concentrations under field conditions. Regression coefficients were evaluated to assess agreements between measured and simulated data. Since the measured data represented average chlorotoluron concentrations within each 2-cm-thick layer, simulated herbicide concentrations were integrated over each layer and divided by its thickness.


View this table:
[in this window]
[in a new window]

  		TABLE 1. Parameters used for reactive solute transport simulations.

 

   Results and Discussion
TOP
 ABSTRACT
 INTRODUCTION
 Mathematical Models
 Materials and Methods
 Results and Discussion
Conclusions
 REFERENCES
 
Definition of Soil Porous Systems
Depending on the type of pedogenesis, soils exhibit different porous systems. Micromorphological images of humic Ap horizons of all soil types are shown in Fig. 2 . The micromorphological image of the humic horizon of Haplic Luvisol (Fig. 2, left) shows higher-order aggregates. The majority of detectable pores with radii larger than 20 µm and corresponding to a pressure head higher 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, and Ck horizons, respectively. All pore radii were smaller than the pore radius (= 740 µm) corresponding to the pressure head of –2 cm that is usually considered to be a limit between capillary pores and macropores (Watson and Luxmoore, 1986). The image of the 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 of larger pores affected to varying degrees by clay coatings and infillings (Bt1, Bt2, BC). Porosities detected on micromorphological images 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).


Figure 2
View larger version (57K):
[in this window]
[in a new window]

  		FIG. 2. Micromorphological images of soil samples characterizing humic horizons of Haplic Luvisol (left), Greyic Phaeozem (middle), and Haplic Cambisol (right): A, pores; B, isolated aggregates; C, clay coatings and fillings; D, grains.

 
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 and a system of larger pores. Radii of the largest circle that could be drawn within each pore were 92, 108, and 107 µm for the Ap, Bth, and Ck horizons, respectively. All pore radii were again smaller than 740 µm. Detectable larger pores, corresponding to a pressure head range between –2 and –70 cm, are separated and affected by clay coatings and infillings that control water and solute interactions between larger pores and pores of the matrix structure. Preferential flow through larger capillary pores may appear in such pore systems. Images of other soil horizons show similar areas of structural pores, either affected (Bth) or not affected (Ck) by clay coatings and infillings. Image porosities were 6.2% for Ap, 5.9% for Bth, and 3.4% for the Ck horizons.

Aggregates in the humic Ap horizon of Haplic Cambisol (Fig. 2, right) are poorly developed. The pore system does not show intrapedal or interpedal pores, pores developed mainly along gravel particles. Radii of the largest circle that could be drawn within the largest image pore were 753, 273, and 318 µm for the Ap, Bw, and C horizons, respectively. In this case, several pores of the Ap horizon had radii larger than 740 µm. Macropores, corresponding to pressure heads higher than –2 cm, may create preferential pathways. Images of other soil horizons show similar features. Image porosities were 32.4% for Ap, 11.8% for Bw, and 12.3% for the C horizon. Detected porosity of the Ap image is significantly influenced by the presence of macropores.

Determination of Soil Hydraulic Properties
Soil hydraulic parameters for the Ap1, Ap2, and Bt1 soil horizons of Haplic Luvisol, shown in Table 2 (which will be used for numerical simulation of the field experiment), were obtained using the numerical inversion of the multistep outflow experiment and the single-porosity model. Measured and calibrated outflow data, as well as optimized and experimental soil water retention curves (not shown here), were satisfactorily similar. Optimized Ks values were significantly higher than those measured using the constant head method on the same soil samples. The 95% confidence intervals for the optimized Ks and l values were very wide, indicating that outflow experiments provide only limited information about the conductivity function. An agreement between measured and simulated outflow data was better when Ks rather than l was optimized.


View this table:
[in this window]
[in a new window]

  		TABLE 2. Single-porosity model parameters of soil hydraulic functions (van Genuchten, 1980) for horizons of Haplic Luvisol.{dagger}

 
Since the micromorphological study of the Ap1 horizon indicated a possible occurrence of immobile zones, the dual-porosity model was also applied in the inverse analysis. Soil hydraulic parameters for the dual-porosity model obtained using the numerical inversion of cumulative outflow and soil-water 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, {theta}r,im, was set equal to {theta}r of the single-porosity model, while the residual water content of the mobile zone, {theta}r,mo, was assumed to be equal to zero, the sum of the saturated water contents of the two pore regions ({theta}s,im and {theta}s,mo) was set equal to {theta}s of the single-porosity model, and the shape parameters {alpha}mo (= {alpha}im) and nmo (= nim) were also considered to be equal to {alpha} and n of the single-porosity model. Additionally, {theta}s,im was defined as the saturated soil water content, {theta}s, minus the measured image porosity. Only the saturated hydraulic conductivity, Ks, and the water transfer coefficient between the immobile and mobile domains, {omega}w, were thus optimized. No significant improvement in the agreement between the measured and calibrated outflow data, or the optimized and experimental soil water retention curves, was observed. Only the difference between optimized and measured Ks values significantly decreased. The 95% confidence 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 significantly decreased for the Bt1 horizon, and its 95% confidence interval remained large. This is probably due to differences between micromorphological features of the Bt1 horizon and the other two horizons.


View this table:
[in this window]
[in a new window]

  		TABLE 3. Dual-porosity model parameters of soil hydraulic functions for horizons of Haplic Luvisol.{dagger}

 
Soil hydraulic parameters for the Ap, Bth, and Ck soil horizons of Greyic Phaeozem shown in Table 4 were obtained using the numerical inversion of the multistep outflow experiment while assuming the single-porosity model. Measured and calibrated outflow data shown for the Ap horizon in Fig. 3 (left) corresponds well. However, better correlation would be desired near saturation and below the pressure head of –400 cm. Worse agreement was again obtained between measured and simulated outflow data when the l parameter was optimized (not shown). The soil water retention curve fit the experimental data well (Fig. 3, right). Optimized Ks values were again significantly higher than those measured. The 95% confidence intervals for the optimized Ks and l values were similarly, as before, very wide.


View this table:
[in this window]
[in a new window]

  		TABLE 4. Single-porosity model parameters of the soil hydraulic functions for horizons of Greyic Phaeozem.{dagger}

 

Figure 3
View larger version (13K):
[in this window]
[in a new window]

  		FIG. 3. Multistep outflow data (left) and soil water retention curves (right) obtained on the soil sample characterizing the humic horizon of Greyic Phaeozem.

 
The dual-permeability model was then applied to improve numerical inversion results and to obtain parameters for two flow domains referred to as matrix (m) and large pore (f) domains (the bicontinuum approach). To obtain unique optimization results for this complex model, many of its parameters must be set equal to independently estimated values. The fraction of the large pore domain (wf = 0.15) was estimated from the fraction of pores identified from micromorphological images. The following parameters were used to define the structure of both domains: β = 8 (the shape factor characterizing angular blocky aggregates), a = 0.5 cm (the characteristic length of an aggregate defining elongated pathways), and {gamma}w = 0.4. Saturated water contents of the matrix and large pore domains, {theta}s,m and {theta}s,f, respectively, were set assuming that the sum of their values multiplied by domains fractions (wm = 1 – wf = 0.85, wf = 0.15) was equal to the measured {theta}s value. Residual water contents of the matrix and large pore domains, {theta}r,m and {theta}r,f, were set equal to {theta}r of the single-porosity model and zero, respectively. Parameters {alpha}f and nf (Table 5 ) for the large pore domain were obtained by fitting the effective soil water retention data characterizing the large pore domain that were defined as soil water contents for particular pressure heads minus the saturated water content of the matrix multiplied by the domain fraction. The lower weighting factor was used for data where the impact of the matrix soil water retention data was expected. Zero soil water contents were set for pressure heads below –120 cm. Considering that large pores control mainly saturated flow, the saturated hydraulic conductivity of the large pore domain, Ks,f, was defined as the ratio of Ks measured using the constant head experiment on the same soil sample and domain fraction (wf). The effective saturated hydraulic conductivity, Ke, for the mass transfer between the matrix and macropore domains was set equal to a small value of 2.4 x 10–3 cm d–1 (10–4 cm h–1) due to the presence of clay coatings. Remaining parameters {alpha}m, nm, and Ks,m were optimized (Table 5) by minimizing the objective function defined using measured cumulative outflow and retention data. Optimized parameter values are shown in Table 5. The resulting total soil water retention curve was obtained as the sum of soil water retention curves for the matrix and macropore domains multiplied by their corresponding fractions.


View this table:
[in this window]
[in a new window]

  		TABLE 5. Dual-permeability model parameters of the soil hydraulic functions for horizons of Greyic Phaeozem.{dagger}

 
The optimized total soil water retention curve and simulated outflow data are presented in Fig. 3. Although the correlation coefficient slightly decreased, Fig. 3 (left) shows that agreement between measured and optimized outflow data for pressure heads close to saturation and below –400 cm exhibited improvement when the dual-permeability model was used. The total saturated hydraulic conductivity (10.6 cm d–1), evaluated as the sum of Ks values for each pore domain multiplied by their corresponding domain fractions, is not substantially higher than the measured Ks value and is considerably lower than Ks obtained using the single-porosity model with l = 0.5. The 95% confidence interval for the optimized Ks,f is also narrower. Numerical inversions using the dual-permeability model were similarly performed for the other two soil samples. The same values of {alpha}f, nf, and Ke were used for both horizons. Final optimized and fixed parameters are shown in Table 5. The total saturated hydraulic conductivity (3.32 cm d–1) for the Bth horizon is again similar to the measured Ks value and is considerably lower than the Ks obtained using the single-porosity model. On the other hand, the total saturated hydraulic conductivity (23.1 cm d–1) for the C horizon is similar to the Ks value obtained using the single-porosity model (for l = 0.5). These differences are likely due to different micromorphological features of the matrix domain of the C horizon compared with the other two horizons.

Soil hydraulic parameters of the single-porosity model for the Ap, Bw, and C soil horizons of Haplic Cambisol obtained by numerical inversion of the outflow data are shown in Table 6 . Calibrated outflow data and optimized soil water retention curves satisfactorily fitted the experimental data. Optimized Ks values were similar to or higher than those measured. The 95% confidence intervals for the Ks value are narrower than for previous soils. Optimization of the l parameter again did not improve agreement between measured and simulated data. The 95% confidence intervals for the l value are again relatively wide.


View this table:
[in this window]
[in a new window]

  		TABLE 6. Single-porosity model parameters of the soil hydraulic functions for horizons of Haplic Cambisol optimized using the multistep outflow data.{dagger}

 
The capillary pore systems of Haplic Cambisol are not bimodal as the capillary pore systems of previous soils. The micromorphological study and field observations for the Ap horizon indicated the presence of gravitational macropores (with pore radii larger than 740 µm). Since the multistep outflow experiment is not particularly suitable for characterization of macropore flow since almost the entire experiment is performed under unsaturated conditions, and since an additional ponded infiltration experiment was performed at this field site, the soil hydraulic parameters for the dual-permeability model were determined in the following way. The soil hydraulic parameters obtained by the optimization process of the multistep outflow experiment using the single-porosity model were assumed to characterize only the matrix pore domain. The additional macropore domain with pores corresponding to pressure heads larger than –2 cm was represented using a steplike shape retention curve with parameters {alpha}f = 0.07 cm–1 and nf = 3, ensuring that macropores were filled with water only for pressure heads close to zero. The macropore domain fraction (wf = 0.05) for all horizons was set based on field observations. The saturated water contents of the macropore domain, {theta}s,f, was set assuming that the sum of {theta}s,m and {theta}s,f multiplied by their corresponding domain fractions (wm = 0.95, wf = 0.05) was equal to the average measured porosities of 0.489, 0.456, and 0.386 cm3 cm–3 for the Ap, Bw, and C soil horizons. The following parameters defining the domain's structure were used: β = 8, a = 0.2 cm, and {gamma}w = 0.4. The Ks values for the macropore domain of all horizons were obtained by numerical inversion of the ponded infiltration experiment that was performed at the experimental plot. The same space discretization was used as for the pesticide transport simulation. Initial pressure head conditions were determined using the matrix soil water retention curve and soil water contents measured in the field. The constant pressure head of 3.5 cm (ponded depth) was used as the upper boundary condition. Free drainage boundary condition was used at the bottom of the soil profile. Similar to Greyic Phaeozem, the effective saturated hydraulic conductivity, Ke, was set equal to 2.4 x 10–3 cm d–1. The resulting parameters are shown in Table 7 .


View this table:
[in this window]
[in a new window]

  		TABLE 7. Single-porosity model parameters of the soil hydraulic functions for horizons of Haplic Cambisol evaluated using the infiltration experiment data.{dagger}

 
Field and Numerical Study of the Herbicide Transport
The chlorotoluron transport in Haplic Luvisol was simulated using the single-porosity and dual-porosity models. In addition to input data described above, the solute transfer coefficient, {omega}s, between the mobile and immobile domains was set equal to 10–3 d–1 (a similar {omega}s value was obtained for the Ap horizon of sandy loam by Köhne et al., 2006b). Cumulative surface water fluxes, root water uptakes and bottom water fluxes simulated between the herbicide application (5 May) and the end of the experiment (10 June) using both models are shown in Fig. 4 . Cumulative surface fluxes are the same since no evaporation was considered and all precipitated water infiltrated into the soil profile. The simulated final root depth was 75 cm in both cases. The cumulative root water uptake simulated using the dual-porosity model is lower than that simulated using the single-porosity model because of the presence of immobile domains and lower hydraulic conductivity, resulting in the higher retention ability of the dual-porosity system. This phenomenon is also reflected in simulated cumulative bottom outflows. Observed concentrations in the soil profiles were expressed as total amounts of solute (present in soil water and adsorbed on soil particles) per mass unit of the dry soil and compared to simulated chlorotoluron concentrations. Measured and calculated chlorotoluron concentrations 5, 13, and 35 d after the application are shown in Fig. 5 . Regression coefficients describing the correlation between the measured and simulated concentrations are presented in Table 8 . Chlorotoluron was regularly distributed in the soil. Although the simulated data using the single-porosity model (Fig. 5, middle) approximated experimental results (Fig. 5, left) reasonably well, the dual-porosity model (Fig. 5, right) provided a slightly better description of the chlorotoluron distribution in the soil profile. The slightly lower herbicide leaching simulated using the dual-porosity model than that simulated using the single-porosity model was caused by the presence of immobile domains and lower hydraulic conductivity.


Figure 4
View larger version (10K):
[in this window]
[in a new window]

  		FIG. 4. Cumulative surface fluxes (CSF), cumulative root water uptakes (CRWU), and cumulative bottom fluxes (CBF) simulated between the herbicide application (5 May) and the end of the experiment (10 June) using the single-porosity (SPM) and dual-porosity (DPM) models in Haplic Luvisol.

 

Figure 5
View larger version (8K):
[in this window]
[in a new window]

  		FIG. 5. Distribution of chlorotoluron in Haplic Luvisol measured (left) and simulated using single-porosity (middle) and dual-porosity (right) models.

 

View this table:
[in this window]
[in a new window]

  		TABLE 8. Regression coefficients describing correlation between measured and simulated concentrations for Haplic Luvisol.

 
The chlorotoluron transport in Greyic Phaeozem was simulated using the single-porosity and dual-permeability models. Additional solute transport parameters required by the dual-permeability model included longitudinal dispersivities in macropores, which were set to 3, 4, and 5 cm for the Ap, Bth, and Ck horizons, respectively, and sorption properties and degradation rates of the large pore domain, which were assumed to be the same as in the matrix domain. Although the soil hydraulic characteristics of the single-porosity model were very similar to the total soil hydraulic characteristics of the dual-permeability model (they were obtained using numerical inversion of the same multistep outflow and soil water retention data), the simulated field water flow and solute transport were significantly different. Cumulative surface water fluxes, root water uptakes, and bottom water fluxes between the herbicide application and the end of the experiment simulated using single-porosity and dual-permeability models are shown in Fig. 6 . Cumulative surface fluxes are again the same. The cumulative root water uptake simulated using the dual-permeability model is higher than that simulated using the single-porosity model, because of the faster saturation of the root zone through preferential pathways, and correspondingly no or low water stress reduction simulated using the dual-permeability model. The simulated final root depth was in both cases 75 cm. The dual-permeability model predicted higher cumulative bottom outflow than the single-porosity model due to again faster saturation of the entire soil profile. The dual-permeability model simulated very low preferential flow at the bottom of the soil profile. The simulated cumulative outflow from the large pore domain was 2 orders of magnitude lower than the simulated outflow from the matrix domain. Measured and calculated chlorotoluron concentrations 5, 13, and 35 d after the application in Greyic Phaeozem are shown in Fig. 7 . Even though regression coefficients (Table 9 ) for the 5th and 35th days are slightly lower, a comparison of plotted simulated and measured chlorotoluron concentrations reveals that the dual-permeability model (Fig. 7, right) described the chlorotoluron transport (Fig. 7, left) better than the single-porosity model (Fig. 7, middle), mainly because the observed herbicide transport was strongly affected by preferential flow. While solute moved only to a depth of 8 cm in the single-porosity system, in the dual-permeability system solute moved to a depth of 26 cm, compared with an observed depth of 22 cm. Observed oscillations in chlorotoluron concentrations between 6 and 22 cm depths were probably caused by preferential flow, that is, a fast solute penetration to greater depths, and a solute accumulation at the end of preferential pathways (due to their disconnection). This phenomenon was not considered in the numerical model.


Figure 6
View larger version (11K):
[in this window]
[in a new window]

  		FIG. 6. Cumulative surface fluxes (CSF), cumulative root water uptakes (CRWU), and cumulative bottom fluxes (CBF) simulated between the herbicide application (5 May) and the end of the experiment (10 June) using the single-porosity (SPM) and dual-permeability (DPM) models in Greyic Phaeozem.

 

Figure 7
View larger version (9K):
[in this window]
[in a new window]

  		FIG. 7. Distribution of chlorotoluron in Greyic Phaeozem measured (left) and simulated using single-porosity (middle), and dual-permeability (right) models.

 

View this table:
[in this window]
[in a new window]

  		TABLE 9. Regression coefficients describing correlation between measured and simulated concentrations for Greyic Phaeozem.

 
The chlorotoluron transport in Haplic Cambisol was again simulated using the single-porosity and dual-permeability models. Similar to Greyic Phaeozem, sorption properties and degradation rates for the large pore domain were assumed to be the same as those in the matrix domain, and the longitudinal dispersivity in macropores was set to 3, 4, and 5 cm for the Ap, Bw, and C horizons, respectively. As expected, the additional macropore domain considerably affected simulated water flow and solute transport. Cumulative water fluxes simulated between the herbicide application and the end of the experiment using the single-porosity and dual-permeability models are shown in Fig. 8 . Surface cumulative fluxes 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 than that simulated using the single-porosity model due to the significantly higher cumulative bottom outflow simulated using the dual-permeability model than that simulated using the single-porosity model. The dual-permeability model simulated almost no preferential flow at the bottom of the soil profile. The simulated cumulative outflow from the large pore domain was four orders of magnitude lower than the simulated outflow from the matrix domain. Measured and calculated chlorotoluron concentrations 5, 13, and 35 d after the application in Haplic Cambisol are shown in Fig. 9 . Regression coefficients describing the correlation between measured and simulated concentrations are presented Table 10 . The numerical simulation that used the single-porosity model (Fig. 9, middle) again underestimated observed herbicide mobility (Fig. 9, left). On the other hand, the dual-permeability model (Fig. 9, right) predicted herbicide movement down to the depth of 40 cm in the macropore domain (not obvious in presented figures due to the very low soil water content and low wf) and consequently large herbicide leaching in the entire flow region.


Figure 8
View larger version (10K):
[in this window]
[in a new window]

  		FIG. 8. Cumulative surface fluxes (CSF), cumulative root water uptakes (CRWU), and cumulative bottom fluxes (CBF) simulated between the herbicide application (5 May) and the end of the experiment (10 June) using the single-porosity (SPM) and dual-permeability (DPM) models in Haplic Cambisol.

 

Figure 9
View larger version (8K):
[in this window]
[in a new window]

  		FIG. 9. Distribution of chlorotoluron in Haplic Cambisol measured (left) and simulated using single-porosity (middle) and dual-permeability (right) models.

 

View this table:
[in this window]
[in a new window]

  		TABLE 10. Regression coefficients describing correlation between measured and simulated concentrations for Haplic Cambisol.

 
The dual-permeability models in both cases (for Greyic Phaeozem and Haplic Cambisol) simulated insignificant bottom drainage through the preferential pathways. However, their presence caused fast matrix saturation in the entire soil profile and consequently, the larger cumulative matrix outflow at the bottom than that simulated using the single-porosity model. The larger chlorotoluron leaching was caused mainly by herbicide transport through the preferential pathways, from which herbicide penetrated into the matrix.


   Conclusions
TOP
 ABSTRACT
 INTRODUCTION
 Mathematical Models
 Materials and Methods
 Results and Discussion
 Conclusions
REFERENCES
 
Collected field data and numerical simulations demonstrated that the multiporous nature of soils has a varying impact on water flow and transport processes at the field scale. The chlorotoluron mobility in the monitored soils increased from Haplic Luvisol to Haplic Cambisol and to Greyic Phaeozem. The pesticide mobility reflects 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 this site compared to the other two locations. The highest mobility of chlorotoluron in Greyic Phaeozem was caused by larger capillary pore pathways and sufficient infiltration fluxes that occasionally filled up these pores. The presence of clay coatings in Greyic Phaeozem that restrict water flow and contaminant transport between the macropore and matrix domains is an additional cause for this preferential transport that produces chlorotoluron penetration into deeper depths. Chlorotoluron was less regularly distributed in Haplic Cambisol. Despite the highest infiltration rate, preferential flow only slightly affected the herbicide transport. Large gravitational pores that may dominate water flow and solute 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 single- and 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.


   ACKNOWLEDGMENTS
 
This work has been supported by the grants No. GA CR 103/05/2143 and MSM 6046070901. The authors acknowledge A. Zigová, K. Nemecek, and O. Drábek for helping with the field and laboratory work, and anonymous reviewers for their valuable remarks and suggestions.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 Mathematical Models
 Materials and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
J. Simunek and S. A. Bradford
Vadose Zone Modeling: Introduction and Importance
Vadose Zone J., May 27, 2008; 7(2): 581 - 586.
[Full Text] [PDF]

Home page
 Vadose Zone JHome page
J. Simunek, M. Th. van Genuchten, and M. Sejna
Development and Applications of the HYDRUS and STANMOD Software Packages and Related Codes
Vadose Zone J., May 27, 2008; 7(2): 587 - 600.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Kodesová, R.
Right arrow Articles by Kozák, J.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Kodesová, R.
Right arrow Articles by Kozák, J.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Kodesová, R.
Right arrow Articles by Kozák, J.
Related Collections
Right arrow Solute Transport Models
Right arrow Water Flow Models

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome