JEQ Grow Your Career With ASA
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

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 Related articles in JEQ
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (13)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Roulier, S.
Right arrow Articles by Jarvis, N.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Roulier, S.
Right arrow Articles by Jarvis, N.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Roulier, S.
Right arrow Articles by Jarvis, N.
Related Collections
Right arrow Pesticides
Right arrow Inverse Procedures/Parameter Estimation
Right arrow Dual Porosity/Permeability Models
Right arrow Nonequilibrium Flow
Published in J. Environ. Qual. 32:2341-2353 (2003).
© ASA, CSSA, SSSA
677 S. Segoe Rd., Madison, WI 53711 USA

TECHNICAL REPORTS

Vadose Zone Processes and Chemical Transport

Modeling Macropore Flow Effects on Pesticide Leaching

Inverse Parameter Estimation Using Microlysimeters

Stéphanie Roulier* and Nicholas Jarvis

Department of Soil Sciences, SLU, Box 7014, 750 07 Uppsala, Sweden

* Corresponding author (stephanie.roulier{at}mv.slu.se).

Received for publication October 16, 2002.

   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Macropore flow is a key factor determining pesticide fate, but models accounting for this process need parameters that cannot be easily measured. This study was conducted to investigate the use of inverse techniques to estimate parameters controlling macropore flow and pesticide fate in the dual-permeability model MACRO. Undisturbed columns were sampled at three landscape positions (hilltop, slope, hollow) with contrasting texture and organic carbon content. Transient leaching experiments were performed for an anionic tracer and the herbicide MCPA (4-chloro-2methylphenoxy acetic acid) during a 4-mo period, first under natural rainfall, and then under controlled irrigation in the laboratory. The tracer breakthrough for the finer-textured soil from the hilltop showed strong evidence of macropore flow, resulting in a rapid leaching of MCPA, while leaching was minimal from the organic-rich hollow soil, since macropore flow was weaker and adsorption stronger. The MACRO model was linked to the inverse modeling program SUFI (Sequential Uncertainty Fitting) to enable calibration of nine key model parameters. Based on calculated model efficiencies, MACRO–SUFI gave generally good predictions of water movement and tracer and pesticide transport, although some errors were attributed to difficulties in simulating the effects of soil moisture on degradation and the timing of water outflows. Even after calibration, significant uncertainties remained for some key parameters controlling macropore flow. Nevertheless, the parameter estimates were significantly different between landscape positions and could also be related to basic soil properties. The posterior uncertainty ranges could probably be reduced with a more exhaustive sampling of the parameter space and improved experimental designs.

Abbreviations: EF, model efficiency • SUFI, Sequential Uncertainty Fitting


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
IN STRUCTURED SOILS, water flowing in macropores may quickly bypass the less permeable soil matrix. Preferential movement of water and solutes in macropores is a key process for the fate of agrochemicals such as pesticides, because it reduces the residence time of the solute in the topsoil, where the chemical reactivity and biological activity is the greatest. Consequently, a fraction of the pesticide moves faster and to greater depths in soils affected by macropore flow and the risk of ground water contamination usually increases (Jarvis, 2002).

In recent years, experimental work investigating the mechanisms that generate and sustain preferential movement of water and pesticides has resulted in significant improvements in our understanding of these complex processes. For example, the extent of preferential flow and transport in soil macropores is known to be influenced by factors such as the initial water content, rainfall intensity, and the timing and method of application of the solute (Coles and Trudgill, 1985; Flury et al., 1994; Magesan et al., 1995; Kätterer et al., 2001). This knowledge and understanding has underpinned the development of a large number of models that include treatments of macropore flow (Feyen et al., 1998; Jarvis, 1998; Simunek et al., 2003). The most widely adopted approach has been to divide the total pore space into two parts. Early approaches assumed that one part of the pore system contained "mobile" water, while the remaining fraction of the water content could be considered stagnant (van Genuchten and Wierenga, 1976). This concept is still widely used, but is best suited to steady water flow under controlled laboratory conditions. In recent years, dual-permeability models have been developed that are better suited to transient field conditions. These models divide the soil porosity into one part characterized by a large storage capacity and small flow capacity (matrix) and another part (macropores) with a small storage capacity and a large flow capacity (Gerke and van Genuchten, 1993; Jarvis, 1994). This type of model has shown promise in recent field tests (Larsson and Jarvis, 1999), and is now starting to be used for management purposes. For example, one widely tested dual-permeability model (MACRO; Jarvis, 1994) is one of four leaching models recommended for use in European Union pesticide registration procedures [Forum for the Coordination of Pesticide Fate Models and Their Use (FOCUS), 1995].

Nevertheless, there are some unresolved issues that hamper the widespread adoption of dual-permeability models in the policy and management arena even though the critical importance of macropore flow for pesticide leaching is generally acknowledged. One difficulty is the general lack of acceptable and reliable methods for estimating model parameters regulating macropore flow [Forum for the Coordination of Pesticide Fate Models and Their Use (FOCUS), 1995], especially since several of these cannot easily be directly measured. Fitting to solute breakthrough curves measured under steady state flow from soil columns is a widely used technique to estimate the parameters of the mobile–immobile model (Jardine et al., 1998). However, this method cannot provide information on the parameters of dual-permeability models that simulate macropore flow under unsaturated, transient conditions in the field. Transient experiments performed in microlysimeters characterized by a short column length (e.g., <30 cm) maintain the advantages of experiments performed on undisturbed soil under controlled conditions, and also allow interpretation of the derived parameters in relation to measured properties of individual soil horizons. Analytical solutions are available for solute transport models based on steady water flow, whereas numerical "inverse" modeling is required to derive the parameters of dual-permeability models from transient microlysimeter experiments. Recently, Schwartz et al. (2000) and Kätterer et al. (2001) have, with mixed success, investigated the use of inverse techniques for estimating the parameters of dual-permeability models. Using generated data, Roulier and Jarvis (2003) demonstrated the crucial role of data quantity and quality in determining the likelihood of achieving reliable inverse parameter estimation for dual-permeability pesticide leaching models, and proposed a methodology for microlysimeter experiments based on the SUFI program (Abbaspour et al., 1997).

Another unresolved issue is the lack of information on the significance of macropore flow at larger scales (e.g., the catchment or landscape), as a function of soil type and topographic position. Surprisingly, hardly any attempts have been made to relate the observed variability in leaching in the presence of macropore flow to variations in fundamental soil properties, although the few studies that have attempted to do so have met with some success. For example, Lennartz (1999) showed that the parameters in a convective log–normal transfer function model of solute transport derived from breakthrough experiments performed on 99 core samples of loamy moraine soil were significantly correlated to clay content. Similarly, Shaw et al. (2000) showed that solute transport was related to the texture and structural development of contrasting soil horizons and developed simple estimation routines for the parameters of the mobile–immobile model based on clay content, cation exchange capacity, and structural development.

The objectives of this study were twofold: (i) to further investigate the feasibility of using inverse modeling, specifically the procedure proposed by Roulier and Jarvis (2003), to derive estimates for key parameters governing the strength of macropore flow and pesticide leaching in dual-permeability models, and (ii) to see the extent to which these derived parameters can be related to variation in fundamental soil properties. To these ends, transient leaching experiments for a tracer (Cl-) and a mobile herbicide (MCPA) were performed on microlysimeters containing undisturbed soil sampled from three landscape positions (hilltop, slope, and hollow) from a field in the loamy moraine catchment of Vemmenhög in southern Sweden. The dual-permeability model MACRO was linked to the inverse modeling program SUFI (Abbaspour et al., 1997) to enable estimation of model parameters related to macropore flow, solute transport, and pesticide leaching by calibrating against the measured data. Differences in the estimated parameter values for the three landscape positions are interpreted in relation to variations in fundamental soil characteristics such as texture and organic carbon content, and direct measurements (Holman and Hollis, 2001) are used to "validate" the estimates where possible.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Soil Properties, Sampling, and Experimental Setup
The soil columns used for this study came from Näsbygård, an intensively instrumented and studied experimental field located in the Vemmenhög catchment (Kreuger, 1999) in southern Sweden (55°26' N, 13°27' E). The field has an undulating topography, and texture and organic carbon content are strongly related to landscape position (Holman and Hollis, 2001). Clay content is larger at the hilltop locations (Eutrudept, according to USDA classification), whereas the hollows are rich in organic carbon (Hapludoll). Four replicate undisturbed soil columns encased in PVC pipes (20 cm in diameter, 20 cm deep) were excavated in early July 2000 from the Ap horizon from each of three representative landscape positions (hilltops, slopes, and hollows) at Näsbygård using a mechanical jack. The columns were transported back to the laboratory, where they were saturated with filtered natural rainwater by slowly raising the water table. The water table was then lowered to the base of the column to reach drainage equilibrium before the start of the tracer and pesticide leaching experiments. The experiments were performed in two steps. On 28 July, doses of KCl solution (31.4 mL at 5000 mg L-1 Cl- = 50 kg ha-1) and the herbicide MCPA (0.628 mL at 5000 mg L-1 a.i. = 1 kg ha-1) were applied to the soil surface of each column, using a hand-held sprayer. They were then placed in the outdoor lysimeter station at SLU Soil Science Department, and solute breakthrough under natural rainfall was monitored during the following 3.5-mo period, with samples collected on a daily basis. This part of the experiment was performed until slightly more than 1 pore volume of drainage (approximately 100 mm) had flowed through the columns. The accumulated potential evaporation was 132 mm, and the total amount of rainfall at the end of the first stage of the experiment was 173 mm (Fig. 1) . The summer and the beginning of the autumn were rather dry, with only small rainfall events in August (Fig. 1), whereas most of the rain (151 mm) came in October and in the first half of November. The columns were then moved to the laboratory on 15 November, and a second dose of KCl (314 mL at 5000 mg L-1 Cl- = 500 kg ha-1) and MCPA (3.14 mL at 5000 mg L-1 a.i. = 5 kg ha-1) was applied to the soil surface. The columns were then irrigated with filtered rain water using a hand-held sprayer. Each irrigation consisted of 8 to 9 mm of water applied every third day, given in intermittent 3-min pulses during a 5-h period, each pulse consisting of an average rate of 25 mm h-1. Outflows were again collected on a daily basis, until approximately half a pore volume of outflow had passed through the columns. For both stages of the experiment, the Cl- concentration in the eluted solution was measured every day, whereas the MCPA concentration was measured on bulked samples for 3-d periods. Finally, the soil was excavated from five soil layers within the columns, each 4 cm in thickness, for extraction of chloride and MCPA to obtain information on the depth profile of resident concentrations at the termination of the experiment. The particle size distribution, the organic carbon content, and pH were also measured in the soil samples. Particle size distribution (<2, 2–20, 20–60, 60–200, 200–600, 600–2000 µm, and >2000 µm) was measured by the sedimentation (pipette) method, and for the larger particle sizes, by wet sieving. The organic carbon content was measured by combustion at 1250°C (CN 2000; LECO Corp., St. Joseph, MI) after removal of carbonate carbon by HCl. The pH was measured in a water and 1 mol L-1 solution of KCl suspension of soil.



View larger version (38K):
[in this window]
[in a new window]
  		Fig. 1. (a) Daily precipitation and (b) potential evaporation (calculated with the Penman equation) during the field phase of the experiment.

 
Analyses
Chloride was extracted from the excavated soil samples by end-over-end shaking in distilled water for 24 h, followed by filtration. Chloride in the water samples was measured by flow injection analysis (FIAstar 5010; FOSS Tecator, Höganäs, Sweden).

The MCPA was extracted from soil samples first by acetone and dichloromethane (1:2) with addition of 4% phosphoric acid and water solution (1:1) and then by an additional 50 mL of dichloromethane. An internal standard of 0.6 µg 2,4,5-TP was added to each sample and to standard solutions to ensure the accuracy and precision of the analytical procedure. After acidification to pH of <2, the solution was extracted by solid-phase extraction (SPE), centrifuged, and filtered before final analysis by liquid chromatography–mass spectrometry (LC–MS). The MCPA analyses were made for samples from only three columns, one representative column from each landscape position. This selection was based on an inspection of the results of the chloride breakthrough for each column.

Simulation Model
The MACRO model (Jarvis, 1994; Larsson and Jarvis, 1999) was used to simulate and analyze the experiments. The MACRO model is a comprehensive physically based dual-permeability model, where the porosity is divided in two domains, characterized by different flow rates and solute concentrations. Only brief descriptions of the most relevant aspects of the model are given here.

The division between the two flow domains is defined by a given "boundary" water potential {psi}b (m) and the corresponding saturated micropore water content {theta}b (m3 m-3), and saturated micropore hydraulic conductivity Kb (m s-1). Water flow in the micropores is governed by the Richards' equation:

[1]
where the subscript mi refers to micropores, {theta} is the water content (m3 m-3), t is the time (s), z is depth (m), K is the hydraulic conductivity (m s-1), {psi} is the soil water potential (m), and Sw is a source–sink term for lateral water exchange between micro- and macroporosity (s-1). In the micropores, the water retention characteristic {psi}({theta}mi) is given by the Brooks and Corey (1964) function, whereas the hydraulic conductivity Kmi({theta}) follows Mualem's (1976) model. Water flow in the macropores is modeled with a modified kinematic wave approach (Germann and Beven, 1985), where macropores are assumed to drain by gravity, and the macropore hydraulic conductivity Kma (and hence flow rate) is expressed as a power law function of the degree of saturation Sma:

[2]
where Ks is the saturated hydraulic conductivity of the total pore system (m s-1) and n* is an empirical "kinematic" exponent accounting for pore size distribution and tortuosity in the macropores.

Water exchange from macro- to micropores is treated as a first-order approximation to the water "diffusion" equation that results from Richards' equation when the influence of gravity is neglected:

[3]
where d is an effective diffusion pathlength (m), Dw is an effective water diffusivity (m2 s-1), and {gamma}w is a scaling factor introduced to match the approximate and exact solutions to the diffusion problem (set to 0.8; for details see Jarvis, 1994). Water "flow" in the reverse direction (micropores to macropores) occurs instantaneously if the micropores become oversaturated.

Solute transport in micropores is calculated using the convection–dispersion equation with source–sink terms representing mass exchange between flow domains Ue and biodegradation Ud:

[4]
where C and S are the solute concentrations in the liquid (kg m-3) and solid phases (kg kg-1), {rho} is the bulk density (kg m-3), f is the fraction of sorption sites in contact with water in macropores, q is the water flow rate (m s-1), and D is the dispersion coefficient (m2 s-1), calculated as the sum of an effective diffusion coefficient and a dispersion term that is linearly dependent on the pore water velocity. An equivalent expression is used for the macropores, except that dispersion is not explicitly considered.

The mass transfer term Ue (kg m-3 s-1) between the macro- and micropores accounts for both diffusion and mass flow:

[5]
where De is the effective diffusion coefficient, Cma and Cmi are the solute concentrations (kg m-3) in the macropores and micropores, respectively, and C' indicates either Cma or Cmi depending on the direction of the water flow. The solute concentration in water routed into the macropores at the soil surface is calculated assuming instantaneous local equilibrium and complete mixing of infiltrating water with the water stored in a shallow layer, or "mixing depth," at the soil surface, zd (m) (Jarvis, 1994).

Pesticide degradation is assumed to follow first-order kinetics and, in this study, to proceed at the same rate in solution and sorbed phases, so that the degradation sink term (kg m-3 s-1) can be expressed as:

[6]
where µref is a user-defined reference first-order rate constant (s-1) adjusted by a temperature response function Ft given by a modified form of the Arrhenius equation (Boesten and van der Linden, 1991) and a soil moisture response function Fw given by a modified form of Walker's function (Walker, 1974):

[7]
where B is an empirical exponent. An instantaneous sorption equilibrium is assumed and a Freundlich isotherm is used to partition the pesticide between solution and adsorbed phases. However, in this study, the Freundlich exponent n was fixed at unity, because a literature survey suggested that MCPA sorption from batch experiments may often be described by a linear isotherm (e.g., Helweg, 1987; Riise et al., 1994; Socías-Viciana et al., 1999) and because Roulier and Jarvis (2003) showed that the two parameters describing the Freundlich adsorption isotherm could not be uniquely identified from flux concentrations and a single resident concentration profile. The linear isotherm is given by:

[8]
where kd is the sorption coefficient (m3 kg-1).

Initial and Boundary Conditions
For the first phase of the experiment in the lysimeter station, hourly rainfall data recorded close to the site (approximately 300 m) were used as driving data in the model, together with daily air temperatures and potential evaporation calculated with the Penman equation. In the laboratory phase of the experiment, evaporation and air temperature were assumed constant at 0.5 mm d-1 and 20°C, respectively, while the irrigation regime was reproduced exactly in the model simulations.

In accordance with the experimental setup, the initial condition was defined in the model as a drainage equilibrium with the water table at the basis of the soil profile, while a "lysimeter" bottom boundary condition was used (zero pressure head, with no inflow).

Model Parameterization and Calibration
The values of nine key parameters in the model were estimated by inverse modeling (Table 1). These parameters were selected on the basis of fulfilling one or more of the following criteria: (i) difficulty and/or impossibility of direct measurement, (ii) large uncertainty in deriving parameter values from highly variable measured data, and/or the uncertainty involved in extrapolating laboratory derived values to the field (e.g., degradation and sorption parameters) (Boesten, 2000), and (iii) large inherent model sensitivity to the parameter (Dubus and Brown, 2002). The parameters selected for calibration on this basis can be grouped into three classes: (i) hydrological parameters related to partitioning of the flow between macropores and micropores and macropore flow itself (saturated micropore hydraulic conductivity Kb; saturated micropore volumetric water content {theta}b; kinematic exponent n*), (ii) parameters related to solute transport (dispersivity Dv; mixing depth zd; effective diffusion pathlength d), and (iii) parameters related to pesticide fate and mobility (degradation rate coefficient µref; the fraction of sorption sites in the macropores f; the sorption coefficient kd).


View this table:
[in this window]
[in a new window]
  		Table 1. Parameters to be determined by inverse modeling, with the initial (prior) uncertainty domains assumed.

 
The inverse modeling procedure was used to estimate the values of these parameters for each column, by calibrating against the measured accumulated water percolation, and the measured tracer and pesticide resident and flux concentrations. For this purpose, the simulation model MACRO was linked to the inverse modeling program SUFI (Abbaspour et al., 1997). The SUFI program is a forward, sequential, and iterative parameter estimation procedure. It is Bayesian in nature, since the procedure starts with user-defined prior uncertainty domains for the parameters to be fitted. Each uncertainty domain is divided into equidistant strata and parameter values are defined by the first moment of each stratum. In this study, four strata were used, based on practical considerations, in particular the number of runs it was possible to achieve bearing in mind the execution speed of the model. The MACRO model was then run for all combinations of parameter values (exhaustive stratified sampling), and the results of the simulations compared with observed variables. The deviation between the observed and corresponding simulated values is quantified by a user-defined objective or goal function. A critical value of the goal function, or tolerance, is also defined. Any parameter combination that gives values of the objective function larger than the tolerance is eliminated. This results in reduced uncertainty domains for each parameter. The next iteration consists of repeating the above steps with the reduced uncertainty domains. The procedure stops either when the goal function reaches a global minimum, or when it is not possible to reduce the uncertainty domains for the next iteration.

The SUFI program is a preferred method to estimate parameters in MACRO, as it employs a global search of the parameter uncertainty domains, and is therefore less likely to fall into local minima of the goal function, an important consideration for highly complex simulation models. Poorly constrained (ill-posed) problems simply leave the prior uncertainty domain for nonsensitive or correlated parameters unchanged (Roulier and Jarvis, 2003).

The way SUFI was applied in this study followed the procedure described in detail by Roulier and Jarvis (2003), and is therefore only briefly explained here. Parameters related to macropore flow and nonreactive solute transport were first calibrated simultaneously on the accumulated water percolation and the tracer resident and flux concentrations. In a second step, the pesticide flux and resident concentrations were used to derive estimates of the parameters related to MCPA sorption and degradation. The root mean square error was used as the goal function, calculated separately for each type of measurement (e.g., water outflow, flux concentration, resident concentration) and then combined by multiplication.

Response Surface Analysis
Constructing response surfaces is a useful way to investigate parameter sensitivity and the uniqueness of the solution in inverse modeling problems. Theoretically, the diffusion pathlength d and the dispersivity Dv might show positive correlation, since both parameters will tend to increase dispersion. Conversely, d and the kinematic exponent n* might show negative correlation since they have opposing effects on macropore flow and transport. As this can affect the reliability of the parameter estimates, we decided to investigate in more detail the behavior of these three parameters using response surfaces of the objective function in the (n*, d) and (Dv, d) planes. The analysis was performed for the data from one hilltop column (the one that showed the strongest evidence of macropore flow). The prior uncertainty domains of d, n*, and Dv (Table 1) were divided into 100, 40, and 35 discrete intervals, respectively. The model was then run for each combination of d and n*, and d and Dv, keeping the remaining parameters constant at the values estimated from the inverse procedure. The goal function was calculated as before.

Pedotransfer Functions
Local pedotransfer functions derived from a database of hydraulic properties and soil characteristics for the Vemmenhög catchment (Lindahl, 2001) were used to estimate those parameters deemed to be less sensitive for the kind of column experiment performed in this study, and easier in principle to measure directly. The bulk density {rho} (g cm-3), the saturated volumetric water content {theta}s, and the Brooks–Corey pore size distribution index in the micropores {lambda} were thus estimated from the particle size distribution and organic carbon content foc (%) measured in the individual soil columns using the relationships:

[9]

[10]

[11]
where {rho}s (2.65 g cm-3) is the particle density and the particle size distribution index µ was derived from fitting the van Genuchten (1980) equation to the accumulative particle size distribution measured in each column (Haverkamp and Parlange, 1986). Finally, the saturated hydraulic conductivity Ks was estimated from the macroporosity as:

[12]
where k = 24300 mm h-1 is an empirical constant calculated from the soil database for the Vemmenhög catchment (Lindahl, 2001). The saturated conductivity Ks was updated for each new value of {theta}b for each simulation in the inverse procedure. The measured soil characteristics and the results of the model parameterization exercise using the pedotransfer functions are shown in Table 2.


View this table:
[in this window]
[in a new window]
  		Table 2. Average values and standard deviations (in parentheses) of soil characteristics and soil hydraulic parameters estimated by pedotransfer functions.

 
To simplify the parameterization problem, the boundary water potential {psi}b was fixed at -10 cm, a pragmatic choice based on experience. Both Seyfried and Rao (1987) and Jardine et al. (1993) demonstrated that steady state solute breakthrough experiments performed at pressure heads larger than -10 cm showed early breakthrough and "tailing," while experiments run at pressure heads less than or equal to -10 cm did not. Also, in applications of dual-permeability models to field data collected under natural rainfall boundary conditions, we have obtained best results when setting {psi}b in this range (i.e., -10 cm) rather than closer to saturation at, for example, -3 or -4 cm (e.g., Larsson and Jarvis, 1999). Other parameters also assumed to be constant are shown in Table 3.


View this table:
[in this window]
[in a new window]
  		Table 3. Constant parameter values for all simulations.

 
Model Evaluation
The model efficiency EF (Loague and Green, 1991) was used to evaluate the performance of the model:

[13]
where N is the number of observations, Oi and Pi are the observed and simulated values, respectively, and is the average value of the observations. If all observed and predicted values are identical, EF will be one (maximum and ideal value), while a negative value of EF indicates a poor fit, meaning that the average value of the observations is a better predictor than the model estimations.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Water Flow and Tracer Leaching
As the replicate columns from each landscape position showed rather similar behavior, only the results for the columns where the pesticide data were available are shown graphically in this paper. The accumulated water percolation from the columns is shown in Fig. 2 . As a consequence of the dry summer (Fig. 1a) almost the entire outflow was observed after the beginning of October. The total amount of percolation was well simulated by MACRO. Nevertheless, the model failed to predict the water outflows from the slope and hollow columns following two large rain events in the beginning of October, perhaps because of an overestimation of soil evaporation during the preceding dry period.



View larger version (20K):
[in this window]
[in a new window]
  		Fig. 2. Comparison between measured and simulated accumulated water percolation for one example column from each landscape element: (a) hilltop, (b) slope, and (c) hollow.

 
The high content of fine material (clay plus silt) in the hilltop columns (Table 2) suggested a soil prone to macropore flow. Indeed, the flux concentrations for the tracer in the hilltop columns (Fig. 3a) highlighted behavior typical of macropore flow. The concentration peak on 16 November, just after the second application of KCl, and the subsequent rapid decrease of concentrations are characteristic of a flow regime dominated by the macroporosity. The tracer resident concentrations also showed distinct evidence of nonequilibrium preferential flow, in that the tracer was detected throughout the profile, but with the largest amounts retained near the surface (Fig. 3b). The slope columns showed a weaker effect of macropore flow: the initial tracer breakthrough was also rapid after the second application of KCl, but the maximum concentration was much smaller, and was only attained following the third water application (Fig. 3c). The resident concentration profile for the tracer was more consistent with convective–dispersive transport, with a maximum concentration at the 8- to 12-cm depth. In the hollow columns, preferential flow behavior was even less evident. The tracer eluted gradually from the columns, with the flux concentration increasing with time throughout the experiment (Fig. 3e).



View larger version (38K):
[in this window]
[in a new window]
  		Fig. 3. Measured and simulated tracer flux concentration for (a) hilltop, (c) slope, and (e) hollow and tracer resident concentration for (b) hilltop, (d) slope, and (f) hollow. The term EF is the model efficiency calculated using Eq. [13].

 
The simulations of water flow and tracer concentrations were globally in good agreement with the measured data. The model efficiencies shown in Table 4 were all positive, except for the resident concentration in one hilltop column (Fig. 3b), where MACRO failed to predict the shape of the profile of chloride content in soil (EF = -0.63), and the tracer flux concentrations in one hollow column (EF = -1.2). The dynamics of the flux concentration were otherwise generally well simulated, especially in the case of strong macropore flow (i.e., the hilltop columns) where MACRO matched well the maximum concentration peak after the second application of KCl, and also the subsequent secondary peaks (Fig. 3a). The model also accurately matched the small flux concentrations measured in all three columns in mid-October to early November resulting from leaching of the first chloride dose (Fig. 3a,c,e). The inability of MACRO to predict the breakthrough of chloride at the beginning of October for the hollow and slope columns (Fig. 3c,e) was due to the fact that no percolation was simulated at this time (Fig. 2).


View this table:
[in this window]
[in a new window]
  		Table 4. Model efficiencies.{dagger}

 


View larger version (39K):
[in this window]
[in a new window]
  		Fig. 6. Measured and simulated pesticide flux concentration for (a) hilltop, (c) slope, and (e) hollow and resident concentrations for (b) hilltop, (d) slope, and (f) hollow. The term EF is the model efficiency.

 
Response Surfaces
Figure 4 shows the results of the response surface analysis for the parameters d, n*, and Dv for one hilltop column. The goal function showed a good sensitivity to the diffusion pathlength d, in both the (n*, d) and (Dv, d) planes, and the response surfaces showed no evidence of parameter correlation, converging toward a minimum in both cases. For the (n*, d) response surface, the minimum lay within the estimated uncertainty domain, although its location did not exactly match with the values estimated from the inverse modeling procedure (Fig. 4a). In addition, the value of the goal function corresponding to the minimum in the (n*, d) plane was 40% smaller than the minimum obtained from the inverse modeling procedure. This was related to the stratification used to build the response surfaces (40 strata for n* and 100 for d), which gave more accurate results than the one used in the inverse modeling procedure (four strata for each parameter and iteration). Using a finer stratification in the inverse procedure would allow a better estimation of the global optimum and reduce the posterior uncertainty domain of the parameters. This is especially true for the last few iterations of the procedure, where the uncertainty domain for each parameter becomes small. At this stage, the number of strata should be larger to maintain the sensitivity of the goal function to the parameters to be estimated and to increase the accuracy of these estimates. Alternatively, the posterior uncertainty domains could be used to define initial values in a subsequent gradient-type optimization procedure to obtain true local minima (Vrugt et al., 2001).



View larger version (61K):
[in this window]
[in a new window]
  		Fig. 4. Response surfaces of the objective function in the (a) (n*, d) and (b) (Dv, d) planes for one hilltop column. The cross identifies the parameter values estimated by the inverse procedure and the hatched area defines the posterior uncertainty domains for these estimates. The shaded scale indicates the values of the goal function.

 
In the (Dv, d) plane, the prior uncertainty domain apparently did not cover the entire area of the minimum of the goal function (Fig. 4b). This problem was also noted by Roulier and Jarvis (2003) who showed that a large range of values in the (Dv, d) plane could provide a reasonably good prediction of the data. It is also apparent that Dv was not such a sensitive parameter under these experimental conditions, where macropore flow dominated the transport process. A different experimental setup (i.e., changed boundary conditions) might be needed to improve the identifiability of this parameter.

Despite an excellent fit to the data (Table 4), the results obtained by inverse modeling for one of the four replicate slope columns appeared to be questionable: the estimated values of the kinematic exponent n* and the diffusion pathlength d were considerably larger for this column (n* = 5.65, d = 186 mm) than for the other three replicates (mean values of 3.13 and 26.1 mm for n* and d, respectively, Table 5). The very large estimate for d seemed particularly suspicious, as the pattern of tracer flux concentrations for this column did not indicate any evidence of macropore flow (results not shown). Thus, a (n*, d) response surface was also built for this specific column, using the same procedure as described previously. The calculated minimum was obtained for n* = 4.25 and d = 3.5 mm, a result that contrasts strongly with the outcome of the SUFI procedure (Fig. 5) , but is much more in accordance with the subjective impression of the absence of macropore flow gained from an inspection of the data. The (n*, d) surface showed several minima, which may explain why the inverse procedure apparently failed for this column. The reason for this is unclear, but the results from this column have been excluded from further analysis. This is a useful reminder that expert judgement must never be neglected when using inverse modeling procedures, even for global search procedures such as SUFI.


View this table:
[in this window]
[in a new window]
  		Table 5. Mean water flow and solute transport parameters estimated by inverse modeling, with ranges given in parentheses. Results from one slope column have been excluded.

 


View larger version (92K):
[in this window]
[in a new window]
  		Fig. 5. Response surface of the objective function in the (n*, d) plane for one slope column. The cross identifies the parameter values estimated by the inverse procedure and the hatched area defines the posterior uncertainty domains for these estimates. The shaded scale indicates the values of the goal function.

 
Parameter Values and Their Uncertainty
The estimated water flow and solute transport parameters are shown in Table 5, while an impression of the uncertainty of these estimates is given in Table 6, which presents the average posterior uncertainty domain for each parameter expressed in dimensionless form (i.e., relative to the mean value). A preliminary check confirmed that this normalized uncertainty domain was independent of the parameter estimate, which means that this measure can be used to compare the level of uncertainty of different parameters. Table 6 shows that large differences exist between parameters that are relatively well identified by the inverse procedure (e.g., {theta}b with an average coefficient of uncertainty of 14%) and those that are poorly defined (e.g., the mixing depth zd with an average coefficient of uncertainty of 197%). It was impossible to reduce the prior uncertainty domain for the mixing depth in one of the hilltop columns. Large posterior uncertainty domains result from either a poor sensitivity or parameter correlation. In some cases the insensitivity may be inherent (Dubus and Brown, 2002) while in others it may result from inadequacies in the experimental design. Relatively large uncertainties in the saturated matrix conductivity Kb and the effective diffusion pathlength d (Table 6) are of particular concern in this respect because they are generally thought to be sensitive for long-term leaching losses (e.g., Dubus and Brown, 2002). Holman and Hollis (2001) reported values of Kb directly measured by tension disc infiltrometer at Näsbygård, close to the sampling locations of the columns. These values are similar to the calibrated values (Table 5), which gives confidence in the inverse procedure used to estimate Kb in this study.


View this table:
[in this window]
[in a new window]
  		Table 6. Uncertainties of the parameter estimates for water flow and tracer transport. Mean values of the coefficient of uncertainty are given, with ranges in parentheses.

 
Relationships with Fundamental Soil Properties
Table 5 shows that the hilltop columns had smaller values for the saturated micropore hydraulic conductivity (mean Kb = 0.61 mm h-1), which will lead to a more frequent activation of the macropores. This is presumably related to the finer texture at the hilltop landscape position. Both Smettem and Bristow (1999) and Jarvis et al. (2002) have demonstrated relationships between near-saturated hydraulic conductivity measured by tension infiltrometer and the fraction of finer soil particles. Furthermore, the hilltop columns also showed large values for the diffusion pathlength (mean d = 96 mm), which implies a slow exchange of solute between the macroporosity and the soil matrix, leading to a stronger physical nonequilibrium. That the diffusion pathlength d also seems closely related to clay content is an important result, implying that soil texture exerts a significant control on soil structural development and thus macropore flow and transport. Similarly, Shaw et al. (2000) found a significant negative relationship between clay content and the first-order mass transfer coefficient in the mobile–immobile solute transport model and attributed this to stronger grades of structure in clayey soil horizons. In contrast to the hilltop columns, the hollow columns had smaller values of the diffusion pathlength (mean d = 13.2 mm) and larger values of the saturated micropore hydraulic conductivity (mean Kb = 2.4 mm h-1), both of which will reduce the significance of macropore flow. Student's t test showed that the values of these two key parameters controlling the strength of macropore flow (d and Kb) were significantly different (P = 0.05) between the hilltop and hollow columns. Holman and Hollis (2001) reported a similar result. In their study, measured Kb was significantly smaller at the hilltop location (0.4 mm h-1) than at the slope and hollow locations (1.11 and 0.82 mm h-1 respectively; Table 5).

Interestingly, the estimated average dispersivity was also largest in the hilltop columns (4 cm) and was significantly larger than in the hollow columns (P = 0.05). One speculative explanation is that the larger content of fines has not only induced stronger macropore structures and therefore macropore flow (as noted above) but also a better developed secondary mesopore system, which results in a greater solute dispersion in the matrix. The saturated matrix water content {theta}b was much larger in the hollow columns with smaller bulk density and the derived estimates were very similar to independent direct measurements made at the site (Table 5). The two remaining parameters (mixing depth zd and the kinematic exponent n*) did not vary significantly between the three landscape positions.

MCPA Leaching
The macropore flow in the hilltop columns highlighted by the tracer data also resulted in a rapid transfer of the pesticide through the column, with a large peak flux concentration observed just after the application of the second dose of pesticide followed by a marked decrease in the rate of pesticide leaching (Fig. 6a) . Moreover, the resident concentration profile also suggests macropore flow, with maximum concentrations near the surface, while some residues of MCPA were also found in the lower layers of the column (Fig. 6b). Despite dry antecedent soil conditions, heavy rainstorms on two occasions during the first phase of the experiment in late August and early October also resulted in large flux concentrations of MCPA (500–1000 µg L-1) but in relatively small amounts of outflow (Fig. 2 and 6a). Although the pattern of leaching was similar, the maximum pesticide flux concentrations in leachate from the slope columns were more than 10 times smaller than from the hilltop columns (see Fig. 6c). This suggests that macropore flow was weaker, in accordance with the results from the tracer experiment. In contrast to the hilltop and slope columns, very little pesticide leached through the hollow columns during the experiment, with flux concentrations less than 2 µg L-1 observed during the laboratory phase of the experiment following the second MCPA dose (Fig. 6e). This is probably due to the combined effects of the absence of macropore flow, as shown by the tracer experiments, and a more effective sorption resulting from the significantly larger organic carbon content in the hollow columns (Table 1). A strong adsorption was also indicated by the fact that no MCPA residues were detected below 4 cm in the hollow column (Fig. 6f). A maximum flux concentration of 4 µg L-1 was observed in the small outflow event in early October during the field phase of the experiment following the first MCPA dose, which might be explained by the presence of shrinkage cracks in the dry soil at this time.

The pesticide parameters estimated by the inverse procedure are shown in Table 7. Despite significant differences in organic carbon content and texture (Table 1), the estimated degradation rates were similar in the three columns, although, of course, no definite conclusions can be drawn because of lack of replication. Degradation was very fast in all columns with the derived rate coefficients being equivalent to first-order half-lives of between 1 and 2 d (Table 7). Widely varying degradation half-lives in the range from 2 to 60 d have been reported for MCPA (e.g., Wauchope et al., 1992; Nicholls, 1994). The relatively fast dissipation in our columns can probably be explained by microbial adaptation, a phenomenon that is often reported for MCPA and other phenoxy-acid herbicides (e.g., Helweg, 1987; Smith and Aubin, 1991). The estimated degradation rate coefficients are similar to those obtained from direct incubation experiments on soil samples collected at Näsbygård in early 2002 (average "pseudo" half-life of 3.5 d, unpublished data, 2002). This gives confidence in the results of the inverse procedure adopted in our study. As expected, the organic-rich hollow column had a sorption coefficient approximately 12 and 2 times larger than the hilltop and slope columns, respectively (Table 7). Adsorption of the weak acid MCPA is affected by soil pH (Riise et al., 1994; Shang and Arshad, 1998). However, since the variation in soil pH was small (Table 1) and pH was much higher than the compound pKa (3.07; Nicholls, 1994), this effect was ignored and effective koc values were estimated from the calibrated sorption coefficient kd and the organic carbon content foc in the columns, following the relation:

[14]


View this table:
[in this window]
[in a new window]
  		Table 7. Values of the pesticide parameters estimated by inverse modeling with the posterior uncertainty domains given in parentheses.

 
The estimated koc values varied between 35 and 147 cm3 g-1, with an average of 89 cm3 g-1, which is within the range of values reported for MCPA in the literature (Helweg, 1987; Riise et al., 1994; Shang and Arshad, 1998; Socías-Viciana et al., 1999). Encouragingly, the degradation rate coefficient seemed to be well identified by the procedure, with the smallest coefficient of uncertainty of all estimated parameters (mean of 6.6%). In contrast, the parameter related to the fraction of sorption sites in the macropores was more uncertain indicating a greater inherent degree of insensitivity for this parameter, as also highlighted by Roulier and Jarvis (2003).

The fit between measured and simulated concentrations of pesticide was generally good (Fig. 6). However, two main failures of MACRO were noted. First, the model was not able to match the observed leaching of pesticide through the hollow column at small flux concentrations (<4 µg L-1, EF = 0; Fig. 6), especially the leaching taking place during the last days of the experiment. The reason for this was that it was not possible to simulate both the presence of MCPA in the leachate and an absence of pesticide residues in the bottom layers of the hollow column (Fig. 6e,f). This apparent contradiction in the observed values is due to the detection limit in the measurement of pesticide concentrations, which is more than one order of magnitude larger for the measurement of the resident concentration (0.01 µg g-1 = approximately 2 µg L-1) than for the flux concentration (0.1 µg L-1). In the future, it is recommended that the influence of the detection limit should be accounted for in the inverse procedure.

Second, MCPA leaching due to macropore flow occurred on a few occasions during the first phase of the experiment, which could not be predicted by the model. In some cases, this was due to a failure to accurately simulate the commencement of water outflow at the beginning of autumn, as already noted. However, even when water outflows were simulated, the model could not predict large transient peak MCPA flux concentrations in late August and early October, especially in hilltop and slope columns (see Fig. 6a; note that small water outflows of <2 mm were simulated on these two occasions, but these are not visible in Fig. 2a). This is probably attributable to differences in degradation rate of MCPA between the first and second doses. With a reference half-life of 1 to 2 d, the first dose of MCPA is predicted to be completely dissipated by late August. The degradation rate coefficient estimated by the inverse procedure is dominated by the observed rate of dissipation of the second dose, both because leaching losses were much larger and because the observed resident concentration profile results primarily from the second dose. A test simulation keeping all other parameters fixed at the calibrated values showed that a value of the reference degradation rate coefficient almost 100 times smaller (approximately 0.006 d-1, or a half-life of 115 d) was needed to be able to match the leaching observed in August and early October. It is possible that microbial adaptation to the second MCPA dose occurred resulting in very fast degradation (half-life of 1–2 d) while the potential degradation rate of the first dose was much smaller. The second dose of MCPA in our study was five times larger than the first, which may also have been a contributing factor. Helweg (1987) showed that fast metabolic degradation of MCPA occurred at high doses on previously untreated soil, but smaller doses failed to trigger this response. However, a reference degradation rate of 0.006 d-1 does not seem at all reasonable for MCPA, even for previously untreated soil (Wauchope et al., 1992; Nicholls, 1994). Furthermore, as noted earlier, independent incubation experiments on soil sampled from Näsbygård in early 2002 confirmed the potential for fast degradation even for the first dose, since the field had not been treated with MCPA since spring 1999 (unpublished data, 2002). A likelier explanation for this apparent anomaly is related to the water response function for degradation (Eq. [7]) in the model (Walker, 1974; Jarvis, 1994). The soil was very dry during the first phase of the experiment due to a large excess of potential evaporation over rain, but the soil was maintained close to saturation during the second (laboratory) phase. A default value of Walker's exponent (B = 0.7) was assumed throughout. It is possible that a larger B value (i.e., a greater sensitivity to drought) would have been more appropriate. Boesten (1986) reported an average value of 0.7 for this exponent from a large number of experiments although he also reported many values larger than unity and one extreme value of 2.8. Another explanation is that the form of the moisture response function may be inappropriate. Helweg (1987) reported that degradation of MCPA was essentially zero at water contents below wilting point, an observation that cannot be reproduced by Eq. [7]. In addition to uncertainty in the parameterization, for numerical reasons the model may also have failed to predict sufficiently dry soil in the critical zone at and very close to the soil surface where the pesticide resides soon after application (Boesten, 2000) and thus overestimated degradation. Clearly, more research is needed on the effects of soil water conditions on pesticide degradation, especially in shallow surface soil layers in the critical period soon after application.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
In this paper we have reported pesticide and tracer leaching experiments performed in columns sampled from three landscape positions (hilltop, slope, hollow) with contrasting soil properties. The analysis of the tracer leaching data showed strong evidence of macropore flow in the fine-textured soil from the hilltop, with consequences for rapid leaching of the herbicide MCPA through the columns. In contrast, macropore flow was much less evident in the organic-rich soil with lower clay content from the hollow, while the large organic carbon content also led to a strong adsorption of the pesticide.

The parameters related to water flow and tracer and pesticide transport were obtained by calibration using an inverse modeling tool based on the dual-permeability MACRO model and the inverse modeling program SUFI. This tool gave generally very good predictions of the measured data based on calculated model efficiencies, and reliable estimation of parameters when these could be checked by independent measurements. Some failures of the model to match the observed data were attributed to errors in simulating the timing of recommencement of water outflows on rewetting dry soil and in the effects of soil moisture on MCPA degradation. A failure to account for analytical detection limits also introduced small errors when calculating the fit between model simulations and measurements in one column. In some cases, the uncertainties in the estimates for some of the key parameters controlling macropore flow remained large even after conditioning MACRO–SUFI against the measured data. This may have been partly due to inherent parameter insensitivity but a more exhaustive sampling of the parameter space (i.e., more sampling strata) would probably have reduced the uncertainty. It is also suggested that further research is needed to identify improved experimental designs that would facilitate estimation of parameters regulating macropore flow and pesticide sorption with less uncertainty. This might include sampling outflows at a higher temporal resolution, the use of multiple tracers of differing characteristics, and successive destructive sampling of replicate soil columns. It would also be important to extend this work in the future to test other types of preferential flow and transport models.

Parameters controlling nonequilibrium macropore flow appeared strongly related to the soil texture and landscape position. The links demonstrated here between the observed variability in leaching characteristics (and the associated key model parameters) and fundamental soil properties suggest that some potential exists to develop pedotransfer functions for the estimation of parameters in dual-permeability models that regulate macropore flow.


   ACKNOWLEDGMENTS
 
This work was funded by a post-doctoral grant to Stéphanie Roulier from the Swedish Natural Sciences Research Council (NFR, now VR), and also partially from the EU project CAMSCALE ("Upscaling, predictive models and catchment water quality") under the contract ENV4-CT97-0439. The authors are also grateful to Dr. Thomas Kätterer (Dep. of Soil Sciences, SLU, Uppsala) and Dr. Karim Abbaspour (EAWAG, Zurich) for helpful discussions concerning SUFI.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES