VZJ sign up for etocs
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 18 July 2005
Published in Vadose Zone J 4:587-601 (2005)
DOI: 10.2136/vzj2004.0153
© 2005 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 Similar articles in ISI Web of Science
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 ISI Web of Science (4)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Coquet, Y.
Right arrow Articles by Roger-Estrade, J.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Coquet, Y.
Right arrow Articles by Roger-Estrade, J.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Coquet, Y.
Right arrow Articles by Roger-Estrade, J.
Related Collections
Right arrow Hydraulic Conductivity/Relative Permeability
Right arrow Spatial Variability
Right arrow Soil Physics
Right arrow Tillage

ORIGINAL RESEARCH

Water and Solute Transport in a Cultivated Silt Loam Soil

2. Numerical Analysis

Y. Coqueta,*, J. Simunekb, C. Coutadeura, M. Th. van Genuchtenc, V. Pota and J. Roger-Estraded

a 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
b Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521
c USDA-ARS, George E. Brown, Jr. Salinity Lab., 450 West Big Springs Road, Riverside, CA 92507
d UMR INRA/INAPG Agronomy, Institut National de la Recherche Agronomique/Institut National Agronomique Paris-Grignon, B.P. 01, 78850 Thiverval-Grignon, France
* Corresponding author (coquet{at}grignon.inra.fr)

Received 7 October 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
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 and solute transport. The experiment involved a 4 by 2 m2 field plot that was uniformly sprinkle irrigated with water and bromide for about 6 h. Independently measured soil hydraulic functions were used to simulate 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 and transport 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 the resulting 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.

Abbreviations: TDR, time domain reflectometry


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
HETEROGENEITY 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 fields of either the hydraulic properties, or linear scaling factors with certain statistical properties such as a mean, variance, and correlation length (Destouni, 1992; Roth and Hammel, 1996; Vanderborght et al., 1997; Hammel et al., 1999).

By contrast, a process-based deterministic approach assumes that soil heterogeneity in space is precisely known, and that the soil properties within the transport domain can be specified exactly at each location. This may be done by describing successive layers within a soil profile based on pedological investigations (Snow et 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 and that 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 a detailed field experiment designed to evaluate the effects of tillage on water and solute transport dynamics in an agricultural field. Experimental details of the field study are given in Part 1 of this two-part series (Coquet et al., 2005). In this paper, hydraulic properties of the various structural components of a tilled silt loam will be analyzed. The properties were measured in situ with a 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 (Simunek 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 field characterization of the heterogeneities involved, including their hydraulic properties, to successfully describe flow and transport at the plot scale.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Description of the Tilled Soil
The soil was a Calcic Cambisol (FAO classification) in an agricultural field at the Institut National de la Recherche 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 mays L.). The structure of the soil was described according to a method proposed by Manichon (1982) based on observations of 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, and the unplowed subsurface. Next, lateral compartments were recognized 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 (delineated by white lines in Fig. 1a), called {Delta} clods (Manichon, 1982), and the remaining soil having a macroporous noncompacted structure, called {Gamma} structure. According to Manichon (1982), these two types of structure may be characterized as follows:



View larger version (88K):
[in this window]
[in a new window]
  		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.

 
In the following, we will use the term "structure type" to refer to either the structure of entire compartments that have a homogeneous structure (e.g., the seed bed, or deeper soil layer unaffected by tillage), or to the structure of single clods ({Delta} or {Gamma}) 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 found in Coquet et al. (2005).

Measurements of the Soil Hydraulic Properties
Water contents ({theta}) and hydraulic conductivities (K) as a function of the pressure head (h) were determined for each soil structure type: the seed bed, the compacted {Delta} soil in the plow layer below the wheel tracks, the macroporous {Gamma} soil in the plow layer between the wheel tracks, compacted {Delta} clods in the plow layer between the wheel tracks, and the untilled soil. A range of field and laboratory methods were used for this 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 a tension 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. This difference in the number of replicates was due to the fact that the volume of soil of each structure type available for the measurements was quite variable. For example, few well-defined {Delta} clods could be identified in the plow layer between the wheel 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 ({Gamma} structure), and the plow layer compacted by the wheel tracks ({Delta} structure) were sampled immediately after surface tillage by collecting large cylinders (15 cm diam., 7 cm height), as well as large blocks ({approx}1 dm3). The cylinders were collected by manually pushing and/or hammering stainless steel sleeves into the soil and then extracting them by removing soil from around the cylinders. The blocks were removed manually while 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 type in 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 {theta}(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 ({approx}5 cm3) extracted from the larger blocks.

Field Experiment
The field experiment is described in detail in Part 1 of this study (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 the periodic soil structural pattern created by tillage, and hence comprised one passage of the tractor pulling the harrow. Thirty time domain reflectometry (TDR) probes (20 cm long) and 30 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 and tensiometers not showing further changes in the water contents and pressure heads), a potassium bromide (KBr) solution of 850 mg Br– L–1 was applied for 2 h at a rainfall intensity of 2.6 cm h–1. The structure of the seed bed was stable and no surface sealing or local water ponding or run-off occurred throughout the experiment. After infiltration, the surface was covered 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 from the face of a trench, 3 m wide and 0.66 m deep, and analyzed for Br– (Coquet et al., 2005).

Numerical Modeling
Governing Flow/Transport Equations
The field experiment was simulated using the HYDRUS-2D software package (Simunek 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

[1]
for saturated-unsaturated water flow and the Fickian-based advection-dispersion equation

[2]
for solute tracer transport. In Eq. [1] and [2] {theta} is the volumetric water content (L3 L–3), h is the pressure head (L), xi (i = 1,2) are the spatial coordinates (L), t is time (T), KijA are components of a dimensionless anisotropy tensor (KA) in the two main spatial directions, K is the unsaturated hydraulic conductivity function (L T–1), c is the solute concentration in the liquid phase (M L–3), qi is the ith component of the volumetric flux 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 elements with irregular geometries (Simunek 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 model of Mualem (1976) to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of soil water retention parameters:

[3]

[4]
in which {theta}r and {theta}s denote the residual (r) and saturated (s) water content (L3 L–3), respectively, Ks is the saturated hydraulic conductivity (L T–1), {alpha} (L–1) and n are shape parameters, m = 1 – 1/n, Se is effective saturation, and l is a pore-connectivity parameter estimated to be about 0.5 as an average for 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 a very minor change in the shape of the water retention curve near saturation, but significantly affects and improves predictions of the hydraulic conductivity function, especially for fine-textured 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 typically found 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 are about 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).

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. This transport domain may be regarded as the basic spatial unit of a periodic pattern created by tillage and being repeated over the field. The plow layer between the wheel tracks was centered in the middle of the transport domain (Fig. 1b). We used an unstructured triangular mesh with 16145 elements to spatially discretize the transport domain. Triangular elements of smaller 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, and about 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: the seed bed, the {Gamma}-structured soil, the compacted {Delta}-unstructured material, and the nontilled subsoil (Manichon, 1982; Coquet et al., 2005). These materials were distributed in accordance with 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 {Delta} soil were located below the wheel tracks between depths of 14 and 28 cm on both sides of the transport domain. Soil between the wheel tracks consisted of a mixture of compacted {Delta} clods and macroporous {Gamma}-structured soil, whose exact locations were based on visual inspection of the soil trenches. The continuous {Delta} zones were created by compaction under the wheels of the tractor pulling the harrow. The {Delta} 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 {Delta} zones under the wheel tracks and the older {Delta} clods between the wheel tracks were assumed to have identical hydraulic properties (Coutadeur et al., 2002) and hence were modeled using HYDRUS-2D as a 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 because of 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 located at 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 the bottom 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 sprayed over 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 then to a zero flux for the remainder of the experiment until about 1000 h. A free drainage (zero pressure head gradient) boundary condition was implemented at the bottom of the soil profile and no flow boundary conditions along both sides. The free drainage condition at the bottom of the transport domain was consistent 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 consisted of 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
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
SUMMARY AND CONCLUSIONS
 REFERENCES
 
Hydraulic Properties of the Tilled Soil
Slightly different hydraulic properties were found for the four soil structure types (i.e., the seed bed, the compacted {Delta} soil, the noncompacted {Gamma} soil, and the untilled soil). Measured soil water retention [{theta}(h)] and hydraulic conductivity [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 {Gamma} soil, the {Delta} soil and the untilled soil, respectively. The main differences between the four structure types are a higher near-saturated hydraulic conductivity for the seed bed than for the {Gamma}-structured soil and the untilled soil, and a lower near-saturated hydraulic conductivity for the {Delta}-structured soil.



View larger version (26K):
[in this window]
[in a new window]
  		Fig. 2. Soil water retention, {theta} (h) (top) and unsaturated hydraulic conductivity, K(h) (bottom) curves for the four soil structure types: (a) seed bed, (b) macroporous {Gamma} structure, (c) compacted {Delta} 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).

 

View this table:
[in this window]
[in a new window]
  		Table 1. Soil hydraulic parameters fitted to data collected using tension-disc infiltrometry, a laboratory evaporation method, a constant-head method, and pressure plates (Fig. 2).

 
The agreement between the tension infiltrometer data and the Wind data was reasonable for K, except for the compacted {Delta} soil (Fig. 2c, bottom). Infiltrometer measurements for the {Delta} 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 the laboratory 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 4 compare measured and simulated water contents and pressure heads, respectively. Because of the manner in which the TDR probes and tensiometers were installed into the face of the transect (Coquet et al., 2005), some uncertainty existed about their exact location in the plow layer relative to the local structures ({Delta} or {Gamma}). This was especially true for the 20-cm-long TDR probes. For example, TDR Probe 19 (Fig. 3b) was located at the edge of a {Delta} block next to a zone with high organic matter content. The soil structure around each probe or tensiometer was determined more precisely when the instruments were removed.



View larger version (40K):
[in this window]
[in a new window]
  		Fig. 3. Measured (symbols) and simulated (lines) water contents of the four different soil structure types: (a) seed bed, (b) compacted {Delta} structure, (c) untilled subsoil below the wheel tracks, (d) macroporous {Gamma} 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.

 


View larger version (30K):
[in this window]
[in a new window]
  		Fig. 4. Measured (symbols) and simulated (lines) pressure heads of the four different soil structure types: (a) seed bed, (b) compacted {Delta} structure, (c) macroporous {Gamma} structure, and (d) untilled subsoil. Numerical simulations were performed using the independently estimated soil hydraulic parameter values in Table 1.

 
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 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 initial water content, calculated from measured pressure heads and the corresponding water retention curve for the seed bed, was significantly larger than the measured initial water content. During the redistribution phase of the experiment, 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 important role. While the evaporation experiment represented a drying branch, our field experiment followed a wetting branch. The extent of hysteresis in the seed bed may have been further enhanced by the strongly aggregated nature 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; Simunek et al., 1998) that field-measured {theta}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 {theta}s values were measured in the laboratory as compared with the field may have been due to allowing for more time for equilibration in the laboratory (thus yielding steady state values with less entrapped and dissolved air), the use of de-aired water, and/or having had water flow constrained one-dimensionally within the soil cylinders. For example, the soil cylinders used for the measurements of K(h) and {theta}(h) with Wind's method were left for at least 3 wk to resaturate in the laboratory. Only for long-term saturated conditions (e.g., groundwater) will the effective field saturation value approach porosity.

Retention curves determined independently from the field experiment were compared with those obtained by relating field-measured water contents and pressure heads at similar soil depths (Fig. 5) . Tensiometers and TDR probes were heuristically paired based on their proximity, 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 the same 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 (only about 0.03–0.05 m3 m–3) for the nontilled, {Delta}-compacted, and macroporous {Gamma} zones (Fig. 3). Differences between measured initial and final water contents of the {Delta} zones and clods during infiltration were much smaller than those predicted using only laboratory parameters. This again can be explained with Fig. 5, which shows that measured field water contents of the {Delta} clods (Fig. 5c) decreased little with increasing tension. It seems that the field {Delta} clods were much more compacted than those used for the laboratory measurements, and hence had narrower pore-size distributions. It is conceivable that the compaction was slightly disturbed during sample collection and transport to the laboratory. Differences between measured and simulated water contents of the nontilled and {Gamma}-structure soils corresponded closely with the differences in water contents close to saturation between the field-constructed retention curves and those measured independently (Fig. 5). It is, however, important to note that data presented in Fig. 5 need to be interpreted with some caution because of the different response times of TDRs and tensiometers (Coquet et al., 2005).



View larger version (21K):
[in this window]
[in a new window]
  		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 {Gamma} structure, (c) compacted {Delta} structure, and (d) untilled subsoil. White squares in (d) are data from an independent internal drainage experiment (Nicole, 2000).

 
Interpretation of the tensiometer measurements was more 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 water content before infiltration and the maximum water content during infiltration) was underestimated with the model as compared to what actually occurred in the field. This interpretation is supported by the TDR measurements (Fig. 3), where the simulated change in water content 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. Highly conductive materials such as those described in this study (with the exception of the {Delta} clods), usually display 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 majority of 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 into and 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 and Gardner (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 in these areas (Fig. 6). Mass balance calculations of the amount of Br– recovered from the profile showed significant lateral redistribution from the {Delta} zones below the wheel tracks (having an average mass balance recovery of 70%), toward profiles next to or between the wheel tracks where mass balance recoveries were between 107 and 118%. These findings are contrary to the numerical simulations which showed water flow and solute 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 {Delta}-structured soil was consistently lower than that of the {Gamma}-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.



View larger version (45K):
[in this window]
[in a new window]
  		Fig. 6. Contours of measured (a) and calculated (b) bromide concentrations 12 h after irrigation. Simulated values were obtained using the independently estimated hydraulic parameter values in Table 1.

 


View larger version (27K):
[in this window]
[in a new window]
  		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 were obtained using the independently estimated hydraulic parameter values in Table 1.

 
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 in Fig. 7b for the area between the wheel tracks. Measured Br– penetrated to a depth of only about 15 cm below the wheel tracks, compared to about 30 cm for the model simulation. Between the wheel tracks, especially in the area immediately next to the compacted {Delta} 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 the wheel 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 {Delta} zones, thus creating horizontal gradients and forcing water and solutes to move laterally around the compacted zones into neighboring areas with higher permeabilities. Having a lower Ks for the {Delta} soil than was estimated in 1999 with tension infiltrometry is supported by the fact that some drying cracks in the {Delta} soil were visible before the tension infiltrometer measurements. It should be noted that the field tension-infiltrometer measurements, the laboratory 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 by changes in the near-saturated hydraulic conductivity of the compacted {Delta} 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 approach to obtain more representative values of the soil hydraulic parameters would be to carry out a complete numerical inversion (Hopmans et al., 2002; Simunek et al., 2002) of the field data. This approach is followed here. 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 selected parameters to obtain more suitable initial estimates of the optimized parameters.

The most important modifications made were those for the soil hydraulic parameters of the seed bed and the compacted {Delta} 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 m3 m–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 {theta}(h) relationship. For reasons discussed previously, the Ks of the compacted {Delta} structure was decreased by one order of magnitude to a value smaller than the simulated rainfall rate. The corrected Ks value of the {Delta} zones was now much closer to the value estimated using Wind's method (0.07 cm h–1). The amount of water between a pressure head of –250 cm and saturation was also adjusted by decreasing the value of {alpha} 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 fact that the field {theta}s was slightly lower than the {theta}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 {theta}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 data much better are listed in Table 2.


View this table:
[in this window]
[in a new window]
  		Table 2. Soil hydraulic parameters heuristically adjusted (upper half) and optimized (lower part) to better describe observed water content and concentration data.

 
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 objective function, {Phi}, 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):

[5]
where m is the number of locations with water content (or pressure head) measurements, nj is the number of measurements at location j, {theta}* and {theta} are measured and modeled water contents at location j and time i, ß is the set of optimized parameters, and w is a weight associated with the measured data that causes the objective function to become the average weighted squared deviation normalized by the measurement variances {sigma}2 (Clausnitzer and Hopmans, 1995).

The parameters {theta}s, n, {alpha}, and Ks of the seed bed, the macroporous {Gamma} structure, and the compacted {Delta} 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, {Phi} (Table 3). Results for this 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).


View this table:
[in this window]
[in a new window]
  		Table 3. Values of the objective function, {Phi}, and regression coefficients, R2, for direct and inverse simulations. Values of {Phi} were defined in terms of either the water content, {theta}, or the pressure head, h.

 
The data in Table 3 show that while heuristic adjustment of the soil hydraulic parameters produced much lower values of {Phi} 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 the correlation coefficients being larger than 0.60 and all smaller than 0.8 (data not shown here). Although it is highly unlikely that, because of the complexity of the field problem, the calibrated parameters represent unique solutions of the inverse problem, they did provide a very good description of the measured water contents.

Water Flow
After calibration of the hydraulic properties, the model gave much more realistic descriptions of the initial water contents, 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 small changes in water contents of the {Delta} clods between the beginning and end of the infiltration experiment (Fig. 8b). We assumed that all {Delta} clods could be characterized with 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.



View larger version (37K):
[in this window]
[in a new window]
  		Fig. 8. Measured (symbols) and simulated (lines) water contents in the four types of soil structure: (a) seed bed, (b) compacted {Delta} structure, (c) untilled subsoil below the wheel tracks, (d) macroporous {Gamma} 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.

 
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 {Delta} soil now caused water and bromide to largely bypass the {Delta} zones below the wheel tracks and to funnel water into areas between the wheel tracks and around individual {Delta} 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 {Delta} clods between the wheel 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 {Delta} clods. The umbrella or shadow effects of the {Delta} 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 the lateral position within the soil). We found the simulated Br– concentration pattern to be very sensitive to the Ks value of the {Delta} structure type, but insensitive to the lateral dispersivity value.



View larger version (46K):
[in this window]
[in a new window]
  		Fig. 9. Contours of measured (a) and calculated (b) bromide concentrations 12 h after irrigation. Simulated values were obtained using the adjusted hydraulic parameter values in Table 2.

 


View larger version (36K):
[in this window]
[in a new window]
  		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.

 
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 the model, 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, 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 the measured distributions.

Preferential Flow
This study shows that preferential flow of water and Br– 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 low-conductivity {Delta} soil structures created by compaction. The low-conductivity regions were not only located below the most recent wheel tracks, but also within the plow layer between the wheel tracks where {Delta} 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 according to the invoked cropping system. The numerical analyses presented here show that the abundance and sizes of the compacted zones within tilled soils will determine the extent of preferential flow. In our example, the presence of large compacted {Delta} 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 more abundant and the larger the {Delta} clods and zones, the more preferential flow should be expected.


   SUMMARY AND CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
REFERENCES
 
A detailed field experiment was performed to study the effects of subsurface heterogeneities generated by farming 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 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 account for possible air entrapment during infiltration, hysteresis and reduced permeability of the compacted soil, and subsequent 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 in describing spatial patterns of water fluxes and solute concentrations, and in predicting double-peak concentration distributions vs. depth for several vertical profiles containing compacted {Delta} clods.

Double-peak concentration profiles are invariably explained in the literature by the presence of macropore flow (Jury et al., 1986; Ghodrati & Jury, 1990; Roth et 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 other parts 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 water and solute pathways. Water is diverted by the less permeable compacted parts to more permeable soil areas, with soil immediately below the clods being protected from infiltrating water and solutes (referred to as umbrella and shadow effects in this study). Lateral diffusion of water and solute below the compacted clods can then produce the secondary concentration peak.

Our results highlight the importance of quantifying the scale of structural heterogeneities that may exist in a 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 structural units created by tillage at macroscopic scales between approximately 1 and 100 cm. The challenge then is to characterize the K(h) and {theta}(h) properties of these heterogeneities 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 these heterogeneities, 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. Simunek'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
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES