Published online July 18, 2005

Water and Solute Transport in a Cultivated Silt Loam Soil: 2. Numerical Analysis
Y. Coquet,* J. Šimůnek, C. Coutadeur, M. Th. van Genuchten, V. Pot, and J. Roger-Estrade

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT

soil profile based on pedological investigations (Snowet al., 1994; Hammel et al., 1999; Vanderborght et al.,2001) and assuming that each soil horizon has its own
soil hydraulic parameters. It is then generally assumed
that the soil horizons are horizontally homogeneous andthat any heterogeneity within the horizon can be neglected. Previous studies, however, have shown that individual soil horizons can also be highly heterogeneous(e.g., Mallants et al., 1996, 1997). In the case of the upper
tilled layer of an agricultural soil, such heterogeneities
can be affected significantly, and explained deterministically, by the invoked tillage practices (Manichon, 1982;
Roger-Estrade et al., 2000).
Several investigations have shown that the soil hydraulic properties are strongly affected by tillage practices (Logsdon et al., 1993; Azevedo et al., 1998; Coutadeur
et al., 2002). Such practices can generate heterogeneities
within the soil in both the horizontal and vertical directions, and may produce various compacted and noncompacted zones and clods (Manichon and Roger-Estrade,1990). The compacted zones within a tilled layer generally
have a much lower hydraulic conductivity than noncompacted zones (Meek et al., 1989; Ankeny et al., 1990).
This paper reports results of a numerical analysis of adetailed field experiment designed to evaluate the effects of tillage on water and solute transport dynamics
in an agricultural field. Experimental details of the fieldstudy are given in Part 1 of this two-part series (Coquet
et al., 2005). In this paper, hydraulic properties of thevarious structural components of a tilled silt loam will
be analyzed. The properties were measured in situ witha tension disc infiltrometer (Coutadeur et al., 2002), as
well as using a variety of laboratory methods, including
Wind’s evaporation method, a constant head method,
and the pressure-plate method. These independently
determined soil hydraulic properties were subsequently
used in HYDRUS-2D (Šimůnek et al., 1999) to deterministically simulate the infiltration of water and bromide (Br⫺) into the tilled soil. The approach assumes explicit knowledge of the geometry (i.e., size and shape)and location of the compacted zones within the soil
profile. By comparing field-measured and numerically
predicted water contents, pressure heads, and tracer
concentrations, we will determine whether one can use
a standard numerical model in combination with fieldcharacterization of the heterogeneities involved, including their hydraulic properties, to successfully describe
flow and transport at the plot scale.

A field experiment was performed to study the effects of soil structure heterogeneity generated by farming practices (i.e., compaction
by wheel traffic, plowing, surface tillage) on plot-scale water flow andsolute transport. The experiment involved a 4 by 2 m2 field plot thatwas uniformly sprinkle irrigated with water and bromide for about
6 h. Independently measured soil hydraulic functions were used tosimulate the experiment with a numerical flow and transport model(HYDRUS-2D) using a fully deterministic approach for describing soil
heterogeneity. The numerical model reproduced observed flow andtransport processes only after adjustments were made to the soil hydraulic functions. Adjustments were needed to account for increased
flow and transport into and through the soil between the compacted zones
below the wheel tracks, and to predict double concentration peaks
caused by the umbrella (or shadow) effects of compacted soil clods.
Global optimization of the soil hydraulic parameters produced a satisfactory description of the very heterogeneous flow patterns, with theresulting hydraulic parameters showing only limited correlation among
each other. We demonstrate that double-peak concentration profiles
can result from the presence of tillage-induced soil heterogeneities.

H

eterogeneity in vadose zone flow and transport
properties is usually modeled using either stochastic or deterministic approaches. In the stochastic approach, either a set of stochastic equations is formulated
and solved for a particular problem (e.g., Yeh et al.,
1985; Mantoglou and Gelhar, 1987) or the governing
deterministic equations are solved for the transport domain using randomly distributed soil properties (e.g.,Roth et al., 1991; Feyen et al., 1998). The latter approach
uses random generators that produce stochastic fieldsof either the hydraulic properties, or linear scaling factors with certain statistical properties such as a mean,variance, and correlation length (Destouni, 1992; Rothand Hammel, 1996; Vanderborght et al., 1997; Hammel
et al., 1999).
By contrast, a process-based deterministic approachassumes that soil heterogeneity in space is precisely
known, and that the soil properties within the transportdomain can be specified exactly at each location. Thismay be done by describing successive layers within a

Y. Coquet, C. Coutadeur, and V. Pot, UMR INRA/INAPG Environment and Arable Crops, Institut National de la Recherche Agronomique/Institut National Agronomique Paris-Grignon, B.P. 01, 78850
Thiverval-Grignon, France; J. Šimůnek, Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521; M.Th. van Genuchten, USDA-ARS, George E. Brown, Jr. Salinity Lab., 450 West Big
Springs Road, Riverside, CA 92507; J. Roger-Estrade, UMR INRA/
INAPG Agronomy, Institut National de la Recherche Agronomique/
Institut National Agronomique Paris-Grignon, B.P. 01, 78850 Thiverval-Grignon, France. Received 7 Oct. 2004. *Corresponding author(coquet@grignon.inra.fr).

MATERIALS AND METHODS
Description of the Tilled Soil

Published in Vadose Zone Journal 4:587–601 (2005).
Original Researchdoi:10.2136/vzj2004.0153© Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

The soil was a Calcic Cambisol (FAO classification) inan agricultural field at the Institut National de la Recherche
Abbreviations: TDR, time domain reflectometry.

587

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

588

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 1. Soil profile perpendicular to the tillage direction showing the different structural components of the plow layer (the untilled soil does
not appear on the photo). Compacted zones are delineated by solid white lines. Also shown in (a) are locations of the TDR probes (circles)
and tensiometers (squares), while (b) shows the spatial distribution of the different soil structure types used in the HYDRUS-2D model.

Agronomique (INRA) Experimental Station at Grignon (Yvelines, France). The field had been plowed in November 2000,
and harrowed in May 2001 to a depth of 15 cm using a 3-m
wide rotary harrow to prepare a seed bed for maize (Zea maysL.). The structure of the soil was described according to amethod proposed by Manichon (1982) based on observationsof the face of a vertical trench, 70 cm deep and 3.1 m wide,
perpendicular to the tillage direction (Fig. 1a). First, several
vertical soil compartments were identified from top to bottom:
the seed bed, the plow layer undisturbed by surface tillage, andthe unplowed subsurface. Next, lateral compartments wererecognized according to the location of the wheel tracks below
which the soil was compacted. The structure of the soil compartment corresponding to the plow layer between the wheel
tracks was further divided into compacted clods (delineatedby white lines in Fig. 1a), called ⌬ clods (Manichon, 1982),
and the remaining soil having a macroporous noncompacted
structure, called ⌫ structure. According to Manichon (1982),
these two types of structure may be characterized as follows:• ⌬ structure, corresponding to soil with no structural porosity, and having smoothly breaking faces. This type ofsoil structure is created by compaction under the wheel
tracks of tractors or other heavy farm machinery. This
⌬ structure is also characteristic of clods located in theplow layer between the wheel tracks. Such ⌬ clods werecreated by plowing of compacted soil formed under wheel
tracks the previous years.
• ⌫ structure, resulting from the coalescence of macroporous

aggregates or clods, with structural porosity clearly visible
to the eye.
In the following, we will use the term “structure type” torefer to either the structure of entire compartments that havea homogeneous structure (e.g., the seed bed, or deeper soil
layer unaffected by tillage), or to the structure of single clods(⌬ or ⌫ ) within heterogeneous compartments (e.g., the compacted clods within the plow layer between the wheel tracks).
Additional details about the soil and its structure can be foundin Coquet et al. (2005).

Measurements of the Soil Hydraulic Properties
Water contents (␪) and hydraulic conductivities (K ) as afunction of the pressure head (h ) were determined for eachsoil structure type: the seed bed, the compacted ⌬ soil in theplow layer below the wheel tracks, the macroporous ⌫ soil inthe plow layer between the wheel tracks, compacted ⌬ clods
in the plow layer between the wheel tracks, and the untilled
soil. A range of field and laboratory methods were used forthis purpose as explained below.
The hydraulic conductivity of the soil at ⫺10, ⫺5, ⫺3, ⫺2,
and 0 cm had previously been measured in the field using atension disc infiltrometer (Coutadeur et al., 2002). One to 18
replicate measurements per pressure head were performed
for each structure type to account for spatial variability. Thisdifference in the number of replicates was due to the fact thatthe volume of soil of each structure type available for the measurements was quite variable. For example, few well-defined

589

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org⌬ clods could be identified in the plow layer between thewheel tracks in 1999 when the infiltrometer measurements
were performed (Coutadeur et al., 2002).
In May 2000, the year before the field experiment took
place, the seed bed, the untilled soil, the macroporous plowed
soil (⌫ structure), and the plow layer compacted by the wheel
tracks (⌬ structure) were sampled immediately after surface
tillage by collecting large cylinders (15 cm diam., 7 cm height),
as well as large blocks (≈1 dm3). The cylinders were collectedby manually pushing and/or hammering stainless steel sleeves
into the soil and then extracting them by removing soil fromaround the cylinders. The blocks were removed manuallywhile using a knife to separate the blocks from the surrounding
soil. One of the large cylinders was used to measure the saturated hydraulic conductivity (Ks) of each soil structure typein the laboratory using the upward constant-head method(Klute and Dirksen, 1986). The other cores, one per soil structure type, were used to simultaneously measure ␪(h ) and K(h )using Wind’s evaporation method (Wind, 1968; Tamari et al.,
1993; Bertuzzi et al., 1999). Finally, water retention values
were measured using Richards pressure plates (Klute, 1986)at ⫺1000, ⫺3300, ⫺10 000, and ⫺15 000 cm on relatively small
clods (≈5 cm3) extracted from the larger blocks.

Field Experiment
The field experiment is described in detail in Part 1 of thisstudy (Coquet et al., 2005). After having been covered by plastic sheets for more than 2 wk, a 4 by 2 m2 area was uniformly
sprinkle irrigated with tap water for 4.33 h at an intensity of
2.1 cm h⫺1. The longer dimension of the plot was perpendicular
to the path of the tractor pulling the harrow for seed bed
preparation. The plot contained one entire spatial unit of theperiodic soil structural pattern created by tillage, and hencecomprised one passage of the tractor pulling the harrow. Thirty
time domain reflectometry (TDR) probes (20 cm long) and30 mini-tensiometers (2 cm long, 0.6 cm outside diam.) were
installed to monitor the soil water dynamics within the tilled
soil (Fig. 1a). However, only 29 TDR probes and 22 tensiometers were found to give reliable results (Coquet et al., 2005).
After quasi-steady state infiltration was reached in the seed
bed and the plow layer (as reflected by the TDR probes andtensiometers not showing further changes in the water contents and pressure heads), a potassium bromide (KBr) solutionof 850 mg Br ⫺ L⫺1 was applied for 2 h at a rainfall intensityof 2.6 cm h⫺1. The structure of the seed bed was stable andno surface sealing or local water ponding or run-off occurred
throughout the experiment. After infiltration, the surface wascovered with plastic sheet to allow water to redistribute within
the soil. Twelve hours after the end of the tracer experiment,
small core samples (4 cm diam., 2 cm height) were taken fromthe face of a trench, 3 m wide and 0.66 m deep, and analyzed
for Br ⫺ (Coquet et al., 2005).

Numerical Modelingfor saturated-unsaturated water flow and the Fickian-based
advection-dispersion equation

冢

冣

⳵␪c
⳵c
⳵q c
⳵
␪Dij
⫺ i
⫽
⳵t
⳵xi
⳵xj
⳵xi

[2]

for solute tracer transport. In Eq. [1] and [2] ␪ is the volumetricwater content (L3 L⫺3), h is the pressure head (L), xi (i ⫽ 1,2)
are the spatial coordinates (L), t is time (T), K Aij are components of a dimensionless anisotropy tensor (KA) in the twomain spatial directions, K is the unsaturated hydraulic conductivity function (L T ⫺1), c is the solute concentration in theliquid phase (M L⫺3), q i is the i th component of the volumetricflux density (L T ⫺1), and Dij are the components of the dispersion coefficient tensor (L2 T ⫺1) for the liquid phase. The
above governing equations were solved with HYDRUS-2D
using Galerkin type linear finite elements applied to an unstructured mesh involving a network of triangular elementswith irregular geometries (Šimůnek et al., 1999).
The unsaturated soil hydraulic properties in this study were
described using van Genuchten (1980) type analytical functions that involve the statistical pore-size distribution modelof Mualem (1976) to obtain a predictive equation for theunsaturated hydraulic conductivity function in terms of soil
water retention parameters:␪ s ⫺ ␪rh⬍0
⎧ ␪r ⫹[1 ⫹ |␣h|n]m
␪(h) ⫽ ⎨
hⱖ0
⎩ ␪s
2

K(h) ⫽ ksS le[1 ⫺ (1 ⫺ S e1/m)m]

[3]
[4]

in which ␪r and ␪s denote the residual (r) and saturated (s)water content (L3 L⫺3), respectively, Ks is the saturated hydraulic conductivity (L T ⫺1), ␣ (L⫺1) and n are shape parameters,m ⫽ 1 ⫺ 1/n , S e is effective saturation, and l is a pore-connectivity parameter estimated to be about 0.5 as an averagefor many soils (Mualem, 1976). A modified van Genuchten(1980) model with an air-entry value of ⫺2 cm (Vogel et al.,2001) was used in all simulations. This model implements avery minor change in the shape of the water retention curvenear saturation, but significantly affects and improves predictions of the hydraulic conductivity function, especially for finetextured soils with small n values (Vogel et al., 2001).
The longitudinal dispersivity for all simulations was assumed
to be equal to 5 cm, the transverse dispersivity was taken to
be one-tenth of the longitudinal dispersivity, and the molecular
diffusion coefficient 0.0675 cm2 h⫺1. The value for the longitudinal dispersivity is consistent with dispersivities typicallyfound for field transport studies (e.g., Jury et al., 1991; Warrick,2003), as well as with observations (e.g., Anderson, 1984)
indicating that the longitudinal dispersivities generally areabout one-tenth of the spatial scale (or distance) of a transport
experiment. Few vadose zone field studies exist where detailed
measurements of both the longitudinal and transverse dispersivities have been made; their ratio is generally assumed
to be about 10, with values typically ranging between about
6 and 20 (e.g., Fetter, 1979).

Governing Flow/Transport Equations
The field experiment was simulated using the HYDRUS-2D
software package (Šimůnek et al., 1999) that numerically predicts two-dimensional water flow and solute transport in variably saturated porous media. The program numerically solves
the Richards equation

冤 冢

⳵␪
⳵
⳵h
K KijA⫹ KizA
⫽
⳵t
⳵xi
⳵xj

冣冥

[1]

Transport Domain
Water flow and solute transport were simulated for a rectangular transport domain, 300 cm wide and 100 cm deep, perpendicular to the path of the tractor pulling the harrow. Thistransport domain may be regarded as the basic spatial unit ofa periodic pattern created by tillage and being repeated overthe field. The plow layer between the wheel tracks was centered in the middle of the transport domain (Fig. 1b). We used

590

VADOSE ZONE J., VOL. 4, AUGUST 2005

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

an unstructured triangular mesh with 16 145 elements to spatially discretize the transport domain. Triangular elements ofsmaller sizes were generated closer to the soil surface and in
the plow layer than near the bottom of the soil profile. The
average size (i.e., the length of the longest side of a triangle)of elements close to the soil surface was about 1.3 cm, andabout 4 cm at or near the bottom of the profile.
Material Distribution
The soil was characterized in terms of four different materials that corresponded to the four soil structure types: theseed bed, the ⌫-structured soil, the compacted ⌬-unstructured
material, and the nontilled subsoil (Manichon, 1982; Coquet
et al., 2005). These materials were distributed in accordancewith field observation (Fig. 1): a 14-cm-thick seed bed overlaying a 16-cm-thick plow layer and a nontilled subsoil below
30 cm. While both the seed bed and the subsoil were assumed
to be homogeneous horizontally, the plow layer was extremely
heterogeneous. Two 80-cm wide compartments of compacted⌬ soil were located below the wheel tracks between depthsof 14 and 28 cm on both sides of the transport domain. Soil
between the wheel tracks consisted of a mixture of compacted⌬ clods and macroporous ⌫-structured soil, whose exact locations were based on visual inspection of the soil trenches. The
continuous ⌬ zones were created by compaction under thewheels of the tractor pulling the harrow. The ⌬ clods between
those wheel tracks resulted from compacted soil under previous wheel tracks that had been cut and displaced by plowing(Roger-Estrade et al., 2000). The ⌬ zones under the wheel
tracks and the older ⌬ clods between the wheel tracks wereassumed to have identical hydraulic properties (Coutadeur
et al., 2002) and hence were modeled using HYDRUS-2D asa single material.
Initial Conditions
HYDRUS-2D allows initial conditions to be specified either
in terms of pressure heads or water contents. While water contents were expected to vary significantly in the soil profile due
to spatial heterogeneity, pressure heads should be relatively
uniform since the soil profile was covered for a substantial
length of time with plastic sheet before initiating the infiltration experiment. Although measured pressure heads decreased somewhat toward the soil surface, presumably becauseof limited evaporation through the plastic cover (Coquet et al.,2005), pressure heads were nearly constant in the horizontal
direction in most or all parts of the soil profile. Minor differences in the initial pressure head between tensiometers locatedat approximately the same depth within the seed bed were
probably due to some wetting by local condensation below
the plastic sheet used to cover the soil before infiltration.
Initial conditions were obtained by averaging measured
initial pressure head values from tensiometers located at approximately the same depths (⫾5 cm). Pressure heads at thebottom of the soil profile were set to a value of ⫺225 cm,
which was the average value recorded by the tensiometers in
the untilled soil before infiltration (Fig. 5 in Coquet et al.,2005). This value was kept constant throughout the untilled
subsoil. Pressure heads were assumed to decrease linearly
from the top of the untilled subsoil toward the soil surface
where they reached a value of ⫺500 cm. The initial Br ⫺ concentration in the solute transport model was set to zero.
Boundary Conditions
The upper boundary condition for water flow corresponded
to a variable flux. The water flux varied from 2.1 cm h⫺1 during

the first 4.33 h when water without any solute was sprayedover the soil surface, to 2.6 cm h⫺1 for the next 2 h during
Br ⫺ application at a concentration of 850 mg L⫺1, and thento a zero flux for the remainder of the experiment until about1000 h. A free drainage (zero pressure head gradient) boundary condition was implemented at the bottom of the soil profileand no flow boundary conditions along both sides. The freedrainage condition at the bottom of the transport domain wasconsistent with the initial soil moisture status of the untilled
bottom layer, while the zero flux conditions along the sides
reflected symmetry in the transport domain (which consistedof a basic unit that was repeated laterally). A third-type boundary condition for solute transport was applied to the soil surface, using a zero input concentration during the initial 4.33 h,
and a concentration of 850 mg L⫺1 during the following 2 h.

RESULTS AND DISCUSSION
Hydraulic Properties of the Tilled Soil
Slightly different hydraulic properties were found forthe four soil structure types (i.e., the seed bed, the compacted ⌬ soil, the noncompacted ⌫ soil, and the untilledsoil). Measured soil water retention [␪(h)] and hydraulicconductivity [K(h)] relationships are given in Fig. 2,
together with the Mualem–van Genuchten functions fitted simultaneously to the retention and conductivity
data using the RETC code (van Genuchten et al., 1991).
Soil hydraulic parameters are given in Table 1. R2 values
were 0.995, 0.958, 0.988, and 0.985, for the seed bed, the⌫ soil, the ⌬ soil and the untilled soil, respectively. The
main differences between the four structure types area higher near-saturated hydraulic conductivity for theseed bed than for the ⌫-structured soil and the untilled
soil, and a lower near-saturated hydraulic conductivity
for the ⌬-structured soil.
The agreement between the tension infiltrometerdata and the Wind data was reasonable for K, except forthe compacted ⌬ soil (Fig. 2c, bottom). Infiltrometer measurements for the ⌬ soil were performed in the upper
part of the compacted soil below the wheel tracks. The
presence of several small shrinkage cracks in this zone
(Coutadeur et al., 2002) may explain the somewhat higher
hydraulic conductivity values measured with the tension
infiltrometer as compared with Wind’s method. The Ks
values measured with the constant head method in thelaboratory and using tension infiltrometry in the field
were very similar. Water retention data in the dry range
measured using Wind’s method generally did not closely
match the measurements on small clods using pressure
plates, except for the seed bed (Fig. 2).

Numerical Analysis Using Independently
Measured Hydraulic Properties
Water Flow
The infiltration experiment was first simulated using
the soil hydraulic parameters determined independently
from the field experiment (Table 1). Figures 3 and 4compare measured and simulated water contents andpressure heads, respectively. Because of the manner inwhich the TDR probes and tensiometers were installed
into the face of the transect (Coquet et al., 2005), some un-Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

591

Fig. 2. Soil water retention, ␪ (h ) (top) and unsaturated hydraulic conductivity, K(h ) (bottom) curves for the four soil structure types: (a) seed
bed, (b) macroporous ⌫ structure, (c) compacted ⌬ structure, and (d) subsoil as determined using tension infiltrometry (open triangle), Wind’s
evaporation method (open diamond), a constant-head method (open square), and using pressure plates (open circle). Solid curves represent
the fitted Mualem-van Genuchten functions. Measured values for Ks were plotted on the left axes at log(⫺h ) ⫽ ⫺2, where h is in cm (or hPa).

certainty existed about their exact location in the plowlayer relative to the local structures (⌬ or ⌫). This was
especially true for the 20-cm-long TDR probes. For example, TDR Probe 19 (Fig. 3b) was located at the edgeof a ⌬ block next to a zone with high organic matter content. The soil structure around each probe or tensiome-

ter was determined more precisely when the instruments
were removed.
Although arrival times of the moisture front were reasonably well predicted for all four soil structure types,
absolute values of the water content were significantly
overestimated with the model (Fig. 3). Overestimation

592

VADOSE ZONE J., VOL. 4, AUGUST 2005

Table 1. Soil hydraulic parameters fitted to data collected using
tension-disc infiltrometry, a laboratory evaporation method, aconstant-head method, and pressure plates (Fig. 2).
Material

␪r

␪s

␣

n

Ks

0.42
0.44
0.40
0.42

cm⫺1
0.104
0.07520.0287
0.107

⫺
1.18
1.12
1.11
1.14

cm h⫺1
42.2
24.4
4.09
29.6

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

m 3 m ⫺3
Seed bed⌫ structure⌬ structure
Subsoil

0.06
0.00
0.00
0.06

was by far the worst for the seed bed where measured
water contents did not exceed 0.3 m3 m⫺3, while predicted values were close to 0.4 m3 m⫺3. Even the initialwater content, calculated from measured pressure heads
and the corresponding water retention curve for theseed bed, was significantly larger than the measured initial water content. During the redistribution phase of theexperiment, the calculated water content of the topsoil
never decreased below 0.2 m3 m⫺3, while the measured
initial water content was only about 0.05 m3 m⫺3. Several
explanations are possible for these significant differences. One possibility is poor contact between the TDR

probes and the relatively loose soil of the seed bed.
However, since all five replicated TDR measurements
produced similar data (Fig. 3a), we considered this possibility to be either unlikely or not significant. Instead, we
believe that a combination of various factors, including
hysteresis, entrapped air, and soil aggregation, caused
the differences between the data and the model predictions. Soil water hysteresis may have played an importantrole. While the evaporation experiment represented adrying branch, our field experiment followed a wetting
branch. The extent of hysteresis in the seed bed may
have been further enhanced by the strongly aggregatednature of the seed bed.
Entrapped air may also have played an important role
during the infiltration phase. It is widely accepted (Klute,
1986; van Genuchten et al., 1991; Šimůnek et al., 1998)
that field-measured ␪s values are generally about 20 to
30% lower than the total porosity, while laboratory values are usually somewhere between field measurements
and the total porosity. The fact that higher ␪s values
were measured in the laboratory as compared with theFig. 3. Measured (symbols) and simulated (lines) water contents of the four different soil structure types: (a) seed bed, (b) compacted ⌬ structure,(c) untilled subsoil below the wheel tracks, (d) macroporous ⌫ structure, and (e) untilled subsoil between the wheel tracks. Vertical grid lines
indicate time period during which bromide was applied with the infiltrating water. Numerical simulations were performed using the independently
estimated soil hydraulic parameter values in Table 1.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

593

Fig. 4. Measured (symbols) and simulated (lines) pressure heads of the four different soil structure types: (a) seed bed, (b) compacted ⌬ structure,(c) macroporous ⌫ structure, and (d) untilled subsoil. Numerical simulations were performed using the independently estimated soil hydraulicparameter values in Table 1.

field may have been due to allowing for more time for
equilibration in the laboratory (thus yielding steady statevalues with less entrapped and dissolved air), the useof de-aired water, and/or having had water flow constrained one-dimensionally within the soil cylinders. For
example, the soil cylinders used for the measurementsof K(h) and ␪(h) with Wind’s method were left for atleast 3 wk to resaturate in the laboratory. Only for longterm saturated conditions (e.g., groundwater) will theeffective field saturation value approach porosity.
Retention curves determined independently from thefield experiment were compared with those obtainedby relating field-measured water contents and pressure
heads at similar soil depths (Fig. 5). Tensiometers andTDR probes were heuristically paired based on theirproximity, provided that they belonged to the same soil
structure type. For instance, TDR probe 5 in the seed
bed was paired with Tensiometer 25 rather than Tensiometer 26, which was located closer to TDR 5 (Fig. 1a),since tensiometer 25 was located at approximately thesame depth as TDR 5. The field-determined retention
curve for the seed bed was shifted in its entirety by
about 0.10 to 0.15 m3 m⫺3 toward lower water contents
as compared with the laboratory-based curve used in
the simulations. Water contents deviated far less (onlyabout 0.03–0.05 m3 m⫺3) for the nontilled, ⌬-compacted,
and macroporous ⌫ zones (Fig. 3). Differences between
measured initial and final water contents of the ⌬ zones
and clods during infiltration were much smaller thanthose predicted using only laboratory parameters. Thisagain can be explained with Fig. 5, which shows thatmeasured field water contents of the ⌬ clods (Fig. 5c)decreased little with increasing tension. It seems thatthe field ⌬ clods were much more compacted than thoseused for the laboratory measurements, and hence hadnarrower pore-size distributions. It is conceivable thatthe compaction was slightly disturbed during sample
collection and transport to the laboratory. Differences between measured and simulated water contents of thenontilled and ⌫-structure soils corresponded closely withthe differences in water contents close to saturation between the field-constructed retention curves and thosemeasured independently (Fig. 5). It is, however, important to note that data presented in Fig. 5 need to beinterpreted with some caution because of the different response times of TDRs and tensiometers (Coquet
et al., 2005).
Interpretation of the tensiometer measurements wasmore difficult than of the TDR data. In almost all cases
the model predicted much earlier arrivals of the moisture front than recorded with the tensiometers (Fig. 4).
This suggests that the quantity of water that can be
stored in the seed bed and/or the plow layer during
infiltration (i.e., the difference between the initial watercontent before infiltration and the maximum water content during infiltration) was underestimated with themodel as compared to what actually occurred in thefield. This interpretation is supported by the TDR measurements (Fig. 3), where the simulated change in watercontent was much lower than the measured change.
With the exception of the seed bed and the plow layer
between the wheel tracks (Fig. 4a, 4c), most tensiometers recorded slow and much more gradual increases in
the pressure heads over time periods of up to 6 h. This
is not consistent with the water content measurements,
which indicated passage of the moisture front in ⬍1 h.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

594

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 5. Comparison of retention curves measured in the laboratory using the evaporation method, and field-measured water contents and pressure
heads for four types of soil structures: (a) seed bed, (b) macroporous ⌫ structure, (c) compacted ⌬ structure, and (d) untilled subsoil. Whitesquares in (d) are data from an independent internal drainage experiment (Nicole, 2000).

Highly conductive materials such as those described inthis study (with the exception of the ⌬ clods), usuallydisplay relatively sharp moisture fronts. The later arrivals of the wetting front as recorded by the tensiometers
relative to those recorded by the TDRs and the numerical model may reflect the different response times of
the tensiometers and TDRs. It seems that the majorityof the tensiometers either did not have sufficient contact
with the surrounding soil or that their response times
were relatively large. A tensiometer connected to a mercury manometer requires water to physically move intoand out of the ceramic cup, thus resulting in a slower
response time than when using a TDR. The soil slurry
added around the tensiometer cups during their installation may also have slowed flow between the surrounding
undisturbed soil and the cup, leading to what Klute andGardner (1962) referred to as soil-limiting tensiometer
response times.
Solute Transport
Numerical calculations produced a relatively uniform
concentration front (Fig. 6 and 7). This was contrary to
the field data, which displayed significantly deeper Br ⫺penetration in the area between the two wheel tracks.
Field measurements showed that water and solute could
not easily enter the compacted zones under the wheel
tracks, and that water largely bypassed these two zones

in favor of the more permeable area between the wheel
tracks. The TDR probes in the compacted zones did
not indicate any significant changes in the water content(Fig. 3b), while also little or no Br ⫺ was measured inthese areas (Fig. 6). Mass balance calculations of theamount of Br ⫺ recovered from the profile showed significant lateral redistribution from the ⌬ zones below
the wheel tracks (having an average mass balance recovery of 70%), toward profiles next to or between thewheel tracks where mass balance recoveries were between 107 and 118%. These findings are contrary to thenumerical simulations which showed water flow andsolute transport to be mostly vertical in the entire plot.
The numerical results were expected since the Ks of all
four structure types were larger than the applied irrigation fluxes (Table 1). While the K(h) curve of the⌬-structured soil was consistently lower than that ofthe ⌫-structured soil for pressure heads above ⫺100 cm(Fig. 2), this difference was not sufficient to induce
strong contrasts in the water flow and Br ⫺ transport
rates within the soil.
Simulated concentration profiles also showed significantly more vertical movement below the wheel tracks
as compared to the measured profiles (Fig. 7a), and did
not reproduce the multiple peaks clearly visible inFig. 7b for the area between the wheel tracks. Measured
Br ⫺ penetrated to a depth of only about 15 cm below

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

595

Fig. 6. Contours of measured (a) and calculated (b) bromide concentrations 12 h after irrigation. Simulated values were obtained using theindependently estimated hydraulic parameter values in Table 1.

the wheel tracks, compared to about 30 cm for the model
simulation. Between the wheel tracks, especially in thearea immediately next to the compacted ⌬ bands, up to
100 mg L⫺1 of Br ⫺ was measured at a depth of 64 cm.
This was contrary to the model simulation, which did
not show deep Br ⫺ penetration in areas between thewheel tracks and immediately to the side of the wheel
tracks.
The hydraulic conductivity of the compacted soil needs
to be smaller than the applied irrigation flux in order
for water to seek alternative pathways. A lower conductivity would lead to higher water contents or even ponding (i.e., perched water) above the compacted ⌬ zones,
thus creating horizontal gradients and forcing water andsolutes to move laterally around the compacted zones
into neighboring areas with higher permeabilities. Having a lower Ks for the ⌬ soil than was estimated in 1999
with tension infiltrometry is supported by the fact thatsome drying cracks in the ⌬ soil were visible before thetension infiltrometer measurements. It should be notedthat the field tension-infiltrometer measurements, thelaboratory soil hydraulic measurements, and the field

tracer experiment were performed in different (successive) years from 1999 to 2001. Although performed at
the same time each season (i.e., immediately after seed
bed preparation), and although also the same cropping
and tillage practices had been implemented each year,
it is possible that the soil hydraulic properties were
not completely identical in the 3 yr as exemplified bychanges in the near-saturated hydraulic conductivity of
the compacted ⌬ structure.

Adjustment of Soil Hydraulic Properties
The discussion above suggests that the independently
measured hydraulic properties may not have been fully
representative of the in situ flow processes at the site,and hence that some modifications to the hydraulic parameters are needed to better describe the measured
water content and Br ⫺ concentration data. One approachto obtain more representative values of the soil hydraulic parameters would be to carry out a complete numerical inversion (Hopmans et al., 2002; Šimůnek et al.,2002) of the field data. This approach is followed here.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

596

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 7. Measured (symbols) and simulated (lines) bromide concentration profiles at two locations with the shallowest (a) and deepest (b) simulated
solute fronts 12 h after irrigation (x ⫽ 78 cm is below the wheel tracks, and x ⫽ 190 cm is between the wheel tracks). Simulated values wereobtained using the independently estimated hydraulic parameter values in Table 1.

However, since the complexity of the two-dimensional
problem with four different materials, and the considerable heterogeneity within the soil profile, make it difficult to formulate an inverse problem with a unique solution, we decided to first modify only a few selectedparameters to obtain more suitable initial estimates ofthe optimized parameters.
The most important modifications made were thosefor the soil hydraulic parameters of the seed bed andthe compacted ⌬ structure (Table 2). Soil hydraulic parameters of the first layer were modified to reflect hysteresis of the seed bed and to allow for the much lower
initial water contents that were measured (i.e., 0.05 m3m⫺3 at a pressure head of about ⫺500 cm). The saturated
and residual water contents for this purpose were decreased and the parameter n was increased from 1.18 to
1.60 to produce more curvature in the ␪(h) relationship.
For reasons discussed previously, the Ks of the compacted ⌬ structure was decreased by one order of magnitude to a value smaller than the simulated rainfall rate.
The corrected Ks value of the ⌬ zones was now muchcloser to the value estimated using Wind’s method (0.07
Table 2. Soil hydraulic parameters heuristically adjusted (upperhalf) and optimized (lower part) to better describe observed watercontent and concentration data.
Material␪r␪s3

Seed bed⌫ structure⌬ structure
Subsoil
Seed bed⌫ structure⌬ structure
Subsoil†† Not optimized.

m ⫺3

␣

n

Ks

cm⫺1

⫺

cm h⫺1

1.60
1.12
1.11
1.14

42.2
24.4
0.1
29.6

1.61
1.10
1.16
1.14

19.44
83.9
0.112
29.6

m
Adjusted hydraulic parameters
0.02
0.30
0.104
0.00
0.38
0.07520.00
0.39
0.00287
0.06
0.38
0.107
Optimized parameters
0.02
0.272
0.01230.00
0.368
0.04910.00
0.361
0.00173
0.06
0.38
0.107

cm h⫺1). The amount of water between a pressure headof ⫺250 cm and saturation was also adjusted by decreasing the value of ␣ by one order of magnitude. We additionally decreased the saturated water content of all
four structure types slightly to account for air entrapment and air dissolution in the irrigating water. The factthat the field ␪s was slightly lower than the ␪s measured
in the laboratory is supported by the plateau values
reached by the water contents during the infiltration experiment (Fig. 3). The decrease in ␪s was kept relatively
small for all soil structure types, except for the seed bed
as discussed earlier. The adjusted parameters that were
found to describe the observed flow and transport datamuch better are listed in Table 2.
Inverse Optimizations
In attempts to further improve the predictions, additional inverse simulations were performed to calibrate
the model against the measured TDR data. Two inverse
simulations with different initial estimates for the optimized parameters were performed. For the first inverse
optimization we used as initial estimates the soil hydraulic parameters determined by fitting the independent
laboratory data (Table 1), while for the second optimization we used the heuristically adjusted parameters (Table 2). For both optimizations we defined the objectivefunction, ⌽, as the sum of the weighted squared deviations between the measured and calculated water contents as follows (Clausnitzer and Hopmans, 1995; Hopmans et al., 2002):m nj⌽ ⫽ 兺 兺 w 关␪*(xj ,zj ,ti ) ⫺ ␪(␤,xj ,zj ,ti )兴

2j⫽1i⫽1w⫽

1m

␴ 兺 nj2

j⫽1[5]

www.vadosezonejournal.org

Table 3. Values of the objective function, ⌽, and regression coefficients, R2, for direct and inverse simulations. Values of ⌽were defined in terms of either the water content, ␪, or thepressure head, h.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Description
Direct simulation using the soil hydraulic parametersfrom Table 1 (⌽ defined in terms of water content)
Direct simulation using the soil hydraulic parametersfrom Table 1 (⌽ defined in terms of pressure head)
Direct simulation using the adjusted soil hydraulicparameters given in Table 2 (⌽ defined in termsof water content)
Direct simulation using the adjusted soil hydraulicparameters given in Table 2 (⌽ defined in terms
of pressure head)
Inverse simulation using the soil hydraulic parametersfrom Table 1 as initial estimates for optimized
parameters (⌽ defined in terms of water content)
Inverse simulation using the adjusted soil hydraulicparameters as initial estimates for optimized
parameters (⌽ defined in terms of water content)

⌽

R2

1.11

0.30

1.77

0.11

0.476

0.74

1.72

0.12

0.344

0.70

0.280

0.74

where m is the number of locations with water content(or pressure head) measurements, nj is the number ofmeasurements at location j, ␪* and ␪ are measured andmodeled water contents at location j and time i, ␤ is theset of optimized parameters, and w is a weight associatedwith the measured data that causes the objective function to become the average weighted squared deviation
normalized by the measurement variances ␴ 2 (Clausnitzer and Hopmans, 1995).
The parameters ␪s , n, ␣, and Ks of the seed bed, themacroporous ⌫ structure, and the compacted ⌬ structure
were optimized simultaneously. We found that the second optimization using the heuristically adjusted initial
hydraulic parameter values produced a much lower final
value of the objective function, ⌽ (Table 3). Results forthis optimization are shown in Table 2. Selected statistical information such as values of the objective function
and the regression coefficients for selected direct and inverse simulations are given in Table 3. Values are given
for both the direct and inverse simulations to provide
a measure of the goodness of fit for particular simulations. They show how the various adjustments and calibrations improved the description of the experimental
data. Values of the objective function for selected simulations were evaluated separately for the measured
water contents and pressure heads (Table 3).
The data in Table 3 show that while heuristic adjustment of the soil hydraulic parameters produced much
lower values of ⌽ if defined in terms of the water content,
the objective function changed only minimally when defined in terms of the pressure head. The inverse analyses, starting with initial estimates using both the original
and adjusted soil hydraulic parameters, in all cases produced substantially lower values of the objective function. Interestingly, little correlation between the optimized parameters was found, with only two values of thecorrelation coefficients being larger than 0.60 and all
smaller than 0.8 (data not shown here). Although it ishighly unlikely that, because of the complexity of the
field problem, the calibrated parameters represent uniquesolutions of the inverse problem, they did provide a very
good description of the measured water contents.

597

Water Flow
After calibration of the hydraulic properties, the modelgave much more realistic descriptions of the initial watercontents, the infiltration front arrival times, the maximum water contents reached in the seed bed during infiltration (Fig. 8a), and the water contents throughout
the soil profile during the field experiment. While not always giving the exact absolute values, the adjusted
model now also closely matched the relatively smallchanges in water contents of the ⌬ clods between the
beginning and end of the infiltration experiment (Fig.
8b). We assumed that all ⌬ clods could be characterizedwith the same soil hydraulic properties, and hence did
not attempt to describe the variable compaction status
of individual clods. Simulated water contents for the remaining structural types (Fig. 8c through 8e) were mostly
within the range of measured values and correctly predicted the arrival of the moisture front.
Solute Transport
The adjusted soil hydraulic properties produced dramatically different concentration profiles than the independently measured values (Fig. 6 and 9). The lower
hydraulic conductivity of the compacted ⌬ soil now
caused water and bromide to largely bypass the ⌬ zones
below the wheel tracks and to funnel water into areas
between the wheel tracks and around individual ⌬ clods.
Figure 9 shows that the simulated concentration map is
now much closer to the measured concentration profile.
Bromide no longer enters the compacted bands under
the wheel tracks, while individual ⌬ clods between thewheel tracks serve as umbrellas (Coquet et al., 2005)
by forcing water and solutes to flow around them, thus
bypassing much of the soil immediately below the ⌬ clods.
The umbrella or shadow effects of the ⌬ clods now produces double concentration peaks in many vertical concentration profiles, as was observed also in the field(Fig. 10). Most of the flow and the deepest Br ⫺ penetration depths were observed and simulated in areas
immediately to the sides of the compacted zones below
the wheel tracks (Fig. 9) at horizontal coordinates between x ⫽ 80 and 100 cm and between x ⫽ 170 and
200 cm (where x is the abscissa corresponding to thelateral position within the soil). We found the simulated
Br ⫺ concentration pattern to be very sensitive to theKs value of the ⌬ structure type, but insensitive to thelateral dispersivity value.
Figure 10 compares measured and simulated vertical
Br ⫺ concentration distributions at several locations in
the transport domain. Since the concentration profiles
below the two (left and right) wheel tracks were very
similar, they were grouped into the same graphs (Fig.
10a and 10k). Numerical simulations for these two areas
did not show much horizontal variation, and closely
matched the measurements. Several vertical profiles between the wheel tracks showed multiple Br ⫺ concentration peaks that were relatively well described with themodel, especially for the verticals at x ⫽ 82, 122, 194,and 210 cm (Fig. 10c, 10f, 10i, and 10j). Bromide concentration profiles for other locations (i.e., x ⫽ 90, 146,

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

598

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 8. Measured (symbols) and simulated (lines) water contents in the four types of soil structure: (a) seed bed, (b) compacted ⌬ structure, (c)untilled subsoil below the wheel tracks, (d) macroporous ⌫ structure, and (e) untilled subsoil between the wheel tracks. Vertical grid lines
indicate time period during which bromide was applied with the infiltrating water. Simulated values were obtained using the optimized soil
hydraulic parameter values in Table 2.

and 186, Fig. 10d, 10g, and 10h) were less well reproduced with the model, although the simulated profiles
qualitatively showed similar features and shapes as themeasured distributions.

Preferential Flow
This study shows that preferential flow of water andBr ⫺ can be initiated without ponding at the soil surface.
However, our numerical analysis also indicates that locally saturated conditions within the soil profile are necessary for lateral flow to develop. These local saturated
or near-saturated conditions occurred above the lowconductivity ⌬ soil structures created by compaction. The
low-conductivity regions were not only located below themost recent wheel tracks, but also within the plow layer
between the wheel tracks where ⌬ clods of various sizes
were present as a result of plowing and related fragmentation of compacted zones that had formed below former
wheel tracks in previous years.
The abundance and size of the compacted soil zones
can be explicitly related to adopted tillage and cropping
practices (Roger-Estrade et al., 2000), and may vary ac-

cording to the invoked cropping system. The numericalanalyses presented here show that the abundance andsizes of the compacted zones within tilled soils will determine the extent of preferential flow. In our example,
the presence of large compacted ⌬ zones was sufficient
to move Br ⫺ down to a depth of 64 cm after only 5.2 cm
of a Br ⫺ solution was added to the soil surface at a
relatively moderate flow rate of 2.6 cm h⫺1. The moreabundant and the larger the ⌬ clods and zones, the more
preferential flow should be expected.

SUMMARY AND CONCLUSIONS
A detailed field experiment was performed to study
the effects of subsurface heterogeneities generated byfarming practices, such as compaction by wheel traffic,
plowing, and superficial tillage, on plot-scale water flow
and solute transport in a silt loam soil. The experiment
involved uniform sprinkler application of water and Br ⫺for a period of 6 h over a surface area of 4 by 2 m2. The
soil profile was instrumented with 30 TDR probes and
30 mini-tensiometers to monitor water dynamics withinReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

599

Fig. 9. Contours of measured (a) and calculated (b) bromide concentrations 12 h after irrigation. Simulated values were obtained using theadjusted hydraulic parameter values in Table 2.

the soil profile. Soil hydraulic properties of the four different structure classes observed in the soil profile were
measured using a combination of field and laboratory
methods. These hydraulic functions were subsequently
used to simulate the experiment.
The HYDRUS-2D numerical code could not reproduce
the observed flow and transport processes when independently measured soil hydraulic properties were used
as input parameters. When such independently measured properties were used, the model predicted relatively uniform Br ⫺ concentration profiles that were inconsistent with the highly heterogeneous distributions
observed in the field, while predicted water contents
during infiltration were also larger than the observed
values.
Much better predictions for both water flow and Br ⫺transport were obtained after making relatively minor
modifications to the soil hydraulic properties to accountfor possible air entrapment during infiltration, hysteresis
and reduced permeability of the compacted soil, andsubsequent model calibration. The adjusted model described the absolute values of the observed water contents, as well as the wetting and solute front arrival times,
reasonably well. The model was then also successful indescribing spatial patterns of water fluxes and soluteconcentrations, and in predicting double-peak concen-

tration distributions vs. depth for several vertical profiles containing compacted ⌬ clods.
Double-peak concentration profiles are invariably explained in the literature by the presence of macropore
flow (Jury et al., 1986; Ghodrati & Jury, 1990; Rothet al., 1991). In this paper we show that double-peak behavior can be caused also by the presence of compacted
clods that have a lower hydraulic conductivity than otherparts of the tilled layer. Such compacted clods are created by compaction under the wheels of farming machinery and subsequent fragmentation and displacement by tillage, especially during moldboard plowing.
The compacted clods can significantly affect the waterand solute pathways. Water is diverted by the less permeable compacted parts to more permeable soil areas,
with soil immediately below the clods being protectedfrom infiltrating water and solutes (referred to as umbrella and shadow effects in this study). Lateral diffusionof water and solute below the compacted clods can thenproduce the secondary concentration peak.
Our results highlight the importance of quantifying
the scale of structural heterogeneities that may exist ina soil. The main idea of Manichon’s method was to study
and quantify the soil profile along a trench while also
considering previous tillage and cropping operations at
the site. The method is based on identification of the mainReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

600

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 10. Measured and simulated concentration profiles versus depth at selected locations. Simulated values were obtained using the adjusted
hydraulic parameter values listed in Table 2.

structural units created by tillage at macroscopic scales
between approximately 1 and 100 cm. The challenge thenis to characterize the K(h) and ␪(h) properties of theseheterogeneities at the right scale. Using field instrumentation or sampling the soil only according to visible surface wheel tracks or the rows of a crop, if present, may
not be sufficient because of the presence of other structural features such as compacted clods within the tilled
soil. More efforts are needed to fully characterize theseheterogeneities, possibly using noninvasive imaging techniques (e.g., Vereecken et al., 2004).
ACKNOWLEDGMENTS
This work was supported in part by the Research and Education Department of the French Ministry of Agriculture,

Alimentation, Fisheries and Rural Affairs. Support from
OECD (a fellowship under the OECD Co-operative Research
Programme: Biological Resource Management for Sustainable Agriculture Systems) and INRA for J. Šimůnek’s summer
stay at the UMR INRA/INAPG EGC-Grignon is greatly acknowledged. This study was supported in part also by SAHRA(Sustainability of Semi-Arid Hydrology and Riperian Areas)
under the STC Program of the National Science Foundation,
Agreement No. EAR-9876800.

REFERENCES
Anderson, M.P. 1984. Movement of contaminants in groundwater:
Groundwater transport–advection and dispersion. p. 37–45. In
Groundwater contamination. Studies in geophysics. National Academy Press, Washington, DC.
Ankeny, M.D., T.C. Kaspar, and R. Horton. 1990. Characterization

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

of tillage and traffic effects on unconfined infiltration measurements.
Soil Sci. Soc. Am. J. 54:837–840.
Azevedo, A.S., R.S. Kanwar, and R. Horton. 1998. Effect of cultivation
on hydraulic properties of an Iowa soil using tension infiltrometers.
Soil Sci. 163:22–28.
Bertuzzi, P., D. Morath, L. Bruckler, J.C. Gaudu, and M. Bourlet.
1999. Wind’s evaporation method: Experimental equipment anderror analysis. p. 323–328. In M.Th. van Genuchten et al. (ed.)
Characterization and measurement of the hydraulic properties of
unsaturated porous media. Part 1. Univ. of California, Riverside.
Clausnitzer, V., and J.W. Hopmans. 1995. Non-linear parameter estimation: LM_OPT. General purpose optimization code based on theLevenberg-Marquardt algorithm. Land, Air and Water Resources
Paper No. 100032, Univ. of California, Davis.
Coquet, Y., C. Coutadeur, C. Labat, P. Vachier, M.Th. van Genuchten,
J. Roger-Estrade, and J. Šimůnek. 2005. Water and solute transport
in a cultivated silt loam soil: I. Field observations. Available atwww.vadosezonejournal.org. Vadose Zone J. 4:573–586 (this issue).
Coutadeur, C., Y. Coquet, and J. Roger-Estrade. 2002. Variation of
hydraulic conductivity in a tilled soil. Eur. J. Soil Sci. 53:619–628.
Destouni, G. 1992. The effect of vertical soil heterogeneity on field
scale solute flux. Water Resour. Res. 28:1303–1309.
Fetter, C.W. 1979. Contaminant hydrogeology. MacMillan Publ. Co.,New York.
Feyen, J., D. Jacques, A. Timmerman, and J. Vanderborght. 1998.
Modelling water flow and solute transport in heterogeneous soils:
A review of recent approaches. J. Agric. Eng. Res. 70:231–256.
Ghodrati, M., and W.A. Jury. 1990. A field study using dyes to characterize preferential flow of water. Soil Sci. Soc. Am. J. 54:1558–1563.
Hammel, K., J. Gross, G. Wessolek, and K. Roth. 1999. Two-dimensional simulation of bromide transport in a heterogeneous fieldsoil with transient unsaturated flow. Eur. J. Soil Sci. 50:633–647.
Hopmans, J.W., J. Šimůnek, N. Romano, and W. Durner. 2002. Simultaneous determination of water transmission and retention properties. Inverse methods. p. 963–1008. In J.H. Dane and G.C. Topp(ed.) Methods of soil analysis. Part 4. Physical methods. SSSA
Book Ser. 5. SSSA, Madison, WI.
Jury, W.A., H. Elabd, and M. Resketo. 1986. Field study of napropamide movement through unsaturated soil. Water Resour. Res.
22:749–755.
Jury, W.A., W.R. Gardner, and W.H. Gardner. 1991. Soil Physics.
John Wiley & Sons, New York.
Klute, A. 1986. Water retention: Laboratory methods. p. 635–662.
In A. Klute (ed.) Methods of soil analysis. Part 1. Physical andmineralogical methods. 2nd ed. Agronomy Monogr. 9. ASA and
SSSA, Madison, WI.
Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity:
Laboratory methods. p. 687–734. In A. Klute (ed.) Methods ofsoil analysis. Part 1. Physical and mineralogical methods. 2nd ed.
Agronomy Monogr. 9. ASA and SSSA, Madison, WI.
Klute, A., and W.R. Gardner. 1962. Tensiometer response time. Soil
Sci. 93:204–207.
Logsdon, S.D., J.L. Jordahl, and D.L. Karlen. 1993. Tillage and crop
effects on ponded and tension infiltration rates. Soil Tillage Res.
28:179–189.
Mallants, D., B.P. Mohanty, D. Jacques, and J. Feyen. 1996. Spatial
variability of hydraulic properties in a multi-layered soil profile.
Soil Sci. 161:167–181.
Mallants, D., P.H. Tseng, N. Toride, A. Timmerman, and J. Feyen.
1997. Evaluation of multimodal hydraulic functions in characterizing a heterogeneous field soil. J. Hydrol. (Amsterdam) 195:172–199.
Manichon, H. 1982. Influence of cropping systems on the soil ‘tillage
profile’: Development of a diagnosis method based on morphological observation. (In French.) Thèse de doctorat de l‘Institut National Agronomique Paris-Grignon, INA-PG, Paris.
Manichon, H., and J. Roger-Estrade. 1990. Characterization of soil
structure and its evolution at short and medium term due to cropping system effects. (In French.) p. 27–55. In D. Picard and L.
Combe (ed.) Les systèmes de culture. INRA Editions, Paris.
Mantoglou, A., and L.W. Gelhar. 1987. Stochastic modeling of large-

601

scale transient unsaturated flow systems. Water Resour. Res. 23:
37–46.
Meek, B.D., E.A. Rechel, L.M. Carter, and W.R. DeTar. 1989.
Changes in infiltration under alfalfa as influenced by time andwheel traffic. Soil Sci. Soc. Am. J. 53:238–241.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:
513–522.
Nicole, J. 2000. In situ determination of the hydraulic properties of a
Grignon silt loam soil by the internal drainage method. (In French.)Unpublished M.Sc. thesis. Hydrologie, Hydrogéologie, Géostatistique et Géochimie. Univ. of Paris XI.
Roger-Estrade, J., G. Richard, H. Boizard, J. Boiffin, J. Caneill, and
H. Manichon. 2000. Modelling structural changes in tilled topsoil
over time as a function of cropping systems. Eur. J. Soil Sci. 51:
455–474.
Roth, K., and K. Hammel. 1996. Transport of conservative chemical
through an unsaturated two-dimensional Miller-similar medium
with steady-state flow. Water Resour. Res. 32:1653–1663.
Roth, K., W.A. Jury, H. Flühler, and W. Attinger. 1991. Transportof chloride through an unsaturated field soil. Water Resour. Res.
27:2533–2541.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1999. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat, and multiple solutes in variably saturated
media. Version 2.0. IGWMC- TPS- 53. International Ground Water
Modeling Center, Colorado School of Mines, Golden, CO.
Šimůnek, J., M.Th. van Genuchten, D. Jacques, J.W. Hopmans, M.
Inoue, and M. Flury. 2002. Solute transport during variably-saturated flow-Inverse methods. p. 1435–1449. In J.H. Dane and G.C.
Topp (ed.) Methods of soil analysis. Part 4. Physical methods. SSSA
Book Ser. 5. SSSA, Madison, WI.
Šimůnek, J., D. Wang, P.J. Shouse, and M.Th. van Genuchten. 1998.
Analysis of a field tension disc infiltrometer experiment by parameter estimation. Int. Agrophysics 12:167–180.
Snow, V., B.E. Clothier, D.R. Scotter, and R.E. White. 1994. Solutetransport in a layered field soil: Experiments and modelling using
the convection-dispersion approach. J. Contam. Hydrol. 16:339–358.
Tamari, S., L. Bruckler, J. Halbertsma, and J. Chadoeuf. 1993. A
simple method for determining soil hydraulic properties in thelaboratory. Soil Sci. Soc. Am. J. 57:642–651.
Vanderborght, J., D. Jacques, D. Mallants, P.H. Tseng, and J. Feyen.
1997. Comparison between field measurements and numerical simulation of steady-state solute transport in a heterogeneous soil
profile. Hydrol. Earth Syst. Sci. 4:853–871.
Vanderborght, J., M. Vanclooster, A. Timmerman, P. Seuntjens, D.
Mallants, D.J. Kim, D. Jacques, L. Hubrechts, C. Gonzalez, J.
Feyen, J. Diels, and J. Deckers. 2001. Overview of inert tracer
experiments in key Belgian soil types: Relation between transportand soil morphological and hydraulic properties. Water Resour.
Res. 37:2873–2888.
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.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC
code for quantifying the hydraulic functions of unsaturated soils.
EPA/600/2–91–065, USEPA, Office of Research and Development,
Washington, DC.
Vereecken, H., S. Hubbard, A. Binley, and T. Ferré (ed.) 2004. Hydrogeophysics. Spec. Sect., Vadose Zone J. 3(4):1060–1239.
Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2001. Effect of
the shape of the soil hydraulic functions near saturation on variablysaturated flow predictions. Adv. Water Resour. 24:133–144.
Warrick, A.W. 2003. Soil water dynamics. Oxford University Press.
Wind, G.P. 1968. Capillary conductivity data estimated by a simple
method. p. 181–191. In P.E. Rijtema and H. Wassink (ed.) Waterin the unsaturated zone. Proc. of the Wageningen Symp. June 1986.
Vol. 1. Gentbrugge/UNESCO, Paris.
Yeh, T.C., L.W. Gelhar, and A.L. Gutjahr. 1985. Stochastic analysis
of unsaturated flow in heterogeneous soils, 1. Statistically isotropic
media. Water Resour. Res. 21:447–456.