Published in Soil Sci. Soc. Am. J. 69:291-300 (2005).
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
Division S-1—Soil Physics
Equation for Describing Solute Transport in Field Soils with Preferential Flow Paths
Young-Jin Kima,
Christophe J. G. Darnaultb,
Nathan O. Baileyd,
J.-Yves Parlangec and
Tammo S. Steenhuisc,*
a Dep. of Biological and Environmental Engineering, Riley-Robb Hall, Cornell Univ., Ithaca, NY 14853, currently at Environmental Sciences Division, Oak Ridge National Lab., Oak Ridge, TN 37831
b Dep. of Civil and Materials Engineering, Univ. of Illinois at Chicago, Chicago, IL 60607
c Dep. of Biological and Environmental Engineering, Riley-Robb Hall, Cornell Univ., Ithaca, NY 14853
d Biological and Agricultural Systems Engineering, Florida A&M Univ., Tallahassee, FL 32307
* Corresponding author (TSS1{at}cornell.edu)
 |
ABSTRACT
|
|---|
Modeling the transport of chemicals is challenging due to the presence of preferential flow paths and the input data needed to describe these paths. We propose a simple equation that can predict the breakthrough of solutes without excessive data requirements under steady state flow conditions. In our generalized preferential flow model (GPFM), the soil is conceptually divided into a distribution zone and a conveyance zone. The distribution zone acts as a linear reservoir resulting in an exponential loss of solutes through preferential flow paths from this zone. In the conveyance zone, the transport of solutes is described with the convective-dispersive (CD) equation. Input data required are the apparent water content of the distribution zone, and solute velocities and dispersivities in the conveyance zone. The model is tested with chloride breakthrough curves (BTCs) resulting from three sets of experiments performed with steady state artificial acid rainfall at different intensities applied sequentially using duplicate columns filled with undisturbed soil or sand both exhibiting preferential flow. The model is able to describe the breakthrough of solutes with physically meaningful parameters through coarse sand with fingered flow and through undisturbed soil cores. In coarse sand, water and solutes flow only through the fingered flow paths in the conveyance zone. In undisturbed columns there is both flow through the matrix and preferential flow paths, and two velocities in the conveyance zone are needed to simulate the breakthrough curves.
Abbreviations: AIC, Akaike information criteria • BTC, breakthrough curve • CD, convective-dispersive • GPFM, generalized preferential flow model • MCE, mean cumulative error
|
INTRODUCTION
|
|---|
SINCE THE DISCOVERY of pesticide contamination of the Long Island aquifers in the early 1980s, both the effect on water quality of preferential flow and their occurrence have been researched widely. The term preferential flow refers to the rapid, nonuniform transport of solutes and water through preferred pathways in the subsoil (Stagnitti et al., 1994). Three categories of preferential flow have been distinguished: macropore flow in well-structured soils (Quisenberry and Phillips, 1976; Beven and Germann, 1982), fingered flow in granular soils and water repellent soils as a consequence of unstable wetting fronts (Hill and Parlange, 1972; Bauters et al., 1998), and funnel flow (Kung, 1990). Dekker and Ritsema (1994) note that preferential flow is more the rule than the exception.
Solute movement prediction in soils was initiated by van der Molen (1956) to predict the rate of desalinization of the Dutch polder soils after inundation by the sea. He derived the CD equation from chromatography theory, based on the assumption that all water percolating through the soil moves approximately with the same velocity. The solute disperses around the solute front that moves down with the average velocity and is described by a dispersion coefficient. It is generally assumed that the dispersion coefficient varies linearly with the average solute velocity (Jury et al., 1991). Subsequently, the CD equation was tested many times successfully with repacked soil columns (e.g., Wierenga and van Genuchten, 1989; Huang et al., 1995; Brush et al., 1999). In the early 1980s, it became clear that under field conditions groundwater contamination of toxic chemicals could not be explained by the CD equation, leading to the development of models that include preferential flow (Ahuja et al., 1993; Steenhuis et al., 1994; Ritsema et al., 1998; Hendriks et al., 1999; Faybishenko et al., 2000). One of earliest type of preferential flow models, formulated by Jury and co-workers, was the transfer function (Jury and Roth, 1990). In this model, the solute flow input response at a certain depth can be calculated from solute flow input response in the layer above when the correlation is known between the points (Nissen et al., 2000). In all cases, the model development was usually ahead of the experimental work by requiring input parameters that were not readily available.
A novel set of experiments was performed by Kung et al. (2000) at four different sites across the USA. They noted that under steady state conditions, the solute transport behavior was surprisingly similar between sites (Kung et al., 2000). This could mean that there is one equation for describing the solute transport in field soils in which water and solutes can move through preferential transport paths. However, no generally accepted analytical expression for describing this type of solute transport has been agreed on. In this paper, we suggest and test such an expression.
|
Generalized Preferential Flow Model Development
|
|---|
One of the ways to model preferential solute movement is by dividing up the soil profile (Fig. 1) into a distribution zone (or an induction zone) near the soil surface and a conveyance zone below (Ahuja and Lehman, 1983; Jarvis et al., 1991; Steenhuis et al., 1994; Ritsema and Dekker, 1995). The distribution zone funnels both water and solutes into flow paths of the conveyance zone (Steenhuis et al., 1994). The thickness of the distribution zone depends on land use and tillage practices. For simplicity, we will assume that initially the conveyance zone has only one set of flow paths through which the water flows with approximately the same velocity. Without loss of generality, the theory can be extended equally well for two (matrix and preferential) or more flow regimes.
For the case when the initial concentration in the distribution zone is C0(kg/m), and the rainfall is solute free, the concentration, C (kg/m), in the percolating water out of the distribution zone can be described as (Jury and Roth, 1990; Steenhuis et al., 1994)
 | [1] |
where t (s) is time, 
(s–1) is the coefficient equal to q/w, q (m/s) is the steady state flow rate, and w (m) is the apparent water content in the distribution zone. For nonadsorbed chemicals, w is simply defined as the volume water in the distribution zone per unit area. For adsorbed chemicals, w = d(
kd + 
s); d (m) is the depth of the distribution zone, s (m3/m3) is the saturated moisture content, (kg/m3) is the bulk density of the soil, and kd (m3/m3) is the desorption partition coefficient. The depth of the distribution zone is dependent on soil type. For plowed soil, it is approximately equal to the depth of plowing. For repacked soil columns, it is the saturated layer near the surface, which is usually around 5 cm deep.
When water is added with a solute concentration, C0, and no solutes initially present, the concentration in the water leaving the distribution zone is
 | [2] |
Equations [1] and [2] are equivalent to those used for sludge by the USEPA (1992) in predicting the loss of metals from the incorporation zone. We further assume that the transport in the preferential flow paths of the conveyance zone can be described with the CD equation (Parker and van Genuchten, 1984)
 | [3] |
where D (m/s) is the dispersion coefficient, and v is the velocity of the solute and equals q/[ß(kd + )], where ß is the mobile region fraction. Using Laplace transforms, Eq. [3] can be integrated for a semi-infinite column subject to the boundary condition described in Eq. [1] and initially no solutes in the column, for 4D/v2 < 1 as (see also Toride et al., 1995)
 | [4] |
where 
= 
. The last term can usually be neglected when x or t are sufficiently large, that is, (x + vt)/(4Dt)1/2 > 3. For the boundary condition in Eq. [2], we find by superposition for sufficient large x or t that the concentration in the conveyance zone equals
 | [5] |
Equation [5] is valid for continuous application of the solute. For a pulse application for which at time t = 
the solute application is stopped and only rainwater is applied, the solute concentration for t > can be found by subtracting the concentration calculated with Eq. [5] at time t – from the concentration calculated at time t. In case there is more than one flow path with different velocities and assuming no mixing between the flow paths, the solute concentration can simply be expressed by summing the contributions of the different flow paths,
 | [6] |
where ai is the fraction of water moving at a velocity of vi through the different flow paths i. Cvi can be found by substituting the velocity vi in either Eq. [4] or Eq. [5], depending on the boundary condition.
|
MATERIALS AND METHODS
|
|---|
Model performance was evaluated with chloride BTCs for three sets of experiments under steady state artificial acid rainfall. The first experiments were performed using duplicate columns of undisturbed soil which exhibited preferential flow, whereas the second and third experiments were performed on duplicate homogenous sand columns (the sand exhibited unstable fingered flow). Steady state rainfall rate was the independent variable (treatment) within each experiment. In Exp. 1 and 2, solute was applied as solution, whereas in Exp. 3, solute was applied as a solid. Acid rain was used with ionic concentrations close to natural rain in the northeastern USA to minimize ionic imbalances in the undisturbed soil cores that had been exposed to this type of rainfall before being brought into the laboratory. The experimental setup is shown in Fig. 2. The artificial rainfall was applied with a rainfall simulator containing six moving needles, rotating in two directions to randomize raindrop distribution (Andreini and Steenhuis, 1990). Chloride was applied either in the artificial rainfall or to the surface after steady state conditions were reached. Drainage water was collected at 10-min time intervals and analyzed for chloride with a digital chloridometer (Buchler Instruments). The artificial acid rain had a pH of 4.1. Its composition is given in Table 1.
For the first set of experiments, two undisturbed Hudson silt loam/silty clay loam (fine, illitic, mesic Glossaquic Hapludalfs) columns (30 cm in diam. and 40 cm high) were hand-excavated in the Cornell University Orchard, in an area with a sod cover with roots in the upper 20 cm and with cracks below the 20-cm depth. The soil structure is described as firm, moderate coarse prismatic parting to medium, subangular blocky (Akhtar et al., 2003). Earthworms and root holes followed the crack network between the hexagonal 25-cm-wide peds. The columns were brought to the laboratory and each subjected to two sequential rainfall intensity treatments, hereafter referred to as cycles, starting with a low rainfall rate cycle of approximately 0.18 cm/h, followed by a second cycle at a rate of 1.5 cm/h (Table 2). Each of the two cycles consisted of three separate parts: first, regular acid rain was applied. After steady state flux conditions were established, acid rain with a chloride solution (1 g/L of CaCl2) was applied at the same intensity. After the chloride concentration in the drainage water equaled the input concentration, regular acid rain was added again in the last part of the cycle. The experiment was stopped or switched to the next intensity in the next cycle after the chloride concentration in the drainage water was undetectable.
View this table:
[in this window]
[in a new window]
|
Table 2. Rainfall intensities and chloride pulse duration applied in units of time and cumulative drainage for the three experiments.
|
|
In the second set of experiments, two homogeneous sand columns (14 cm in diam. and 40 cm high) were used. The columns were filled with dry and sieved commercial 12/20 sand (1.1 mm average diam., Union Corporation) and packed using a vibrator. The water and chloride application was similar to the first experiments, with the difference that three cycles with different intensities were used. In the first cycle, the intended flow rate was 0.4 cm/h, but one of the pumps was misadjusted and the application rate was 0.05 cm/h, which was too low for the pump to function well. In the second cycle the high application rate was around 1.7 cm/h, and in the third cycle the low intensity was approximately the same as the intended rate for the first cycle (Table 2).
For the third set of experiments, six columns, two for each of the three rainfall intensity cycles, were used as part of an experiment on the movement of Cryptosporidium parvum oocysts and are of interest here because chloride was applied in the feces. The type of columns (14 cm in diam. and 40 cm high) and the sand (1.1-mm average diam.) were the same as in Exp. 2. In contrast to Exp. 1 and 2, chloride was applied in the feces. For each of the three rainfall intensity cycles, a different set of two columns was used. Regular acid rain was applied until steady state flux conditions were reached, after which cow manure mixed with 2 g of NaCl was added to the columns and the rainfall was continued. At the end of the experiment, a 0.5% FD&C No. 1 blue dye was added to the rain and the columns were segmented and the flow pattern observed. The rainfall intensities employed were approximately 0.3, 1, and 2 cm/h (Table 2).
Model Comparisons
To compare the experimental results with the model data, graphical and statistical criteria can be used. Graphical comparisons between model and observed data are subjective; however, no best statistical-quality criterion has been identified for hydrologic models (Weglarczyk, 1998). Therefore, three assessment criteria were used to compare the calibrated output of both the GPFM and the CD model with the observed values; this was done in accordance with the standard R2 correlation method, the mean cumulative error (MCE) described by Perrin et al. (2001), and the Akaike information criteria (AIC) that takes into account the number of parameters fitted in the regression model (Hippel, 1981; Russo and Jury, 1987). The R2 ranges from 0 to 1, while the MCE ranges from –
to 1, where 1 represents a perfect fit. The AIC is a relative measure between two models with the lowest AIC value indicating the best fit. The first criterion, the R2 correlation method, can be found in most text books and is not given here.
The second criterion, transformed to a scale of – to 1, was derived from the mean absolute model error (Perrin et al., 2001). While the R2 yield is an assessment of how well each of the individual points fit the observed data, MCE measures the ability of a model to correctly reproduce the concentrations during the entire experimental period. The MCE is given as (Perrin et al., 2001):
 | [7] |
The third criterion is the AIC and can be expressed as (Webster and McBratney, 1989)
 | [8] |
where L is the maximum likelihood estimate and p is the number of parameters in the regression equation. According to Larget (2003), Eq. [8] can be given in a regression context as:
 | [9] |
where n is the number of observed and predicted concentration pairs (Cobs,i and Csim,i, respectively). There are other forms of the AIC that have different constants, but since the AIC value of the two models is compared, this will not affect the outcome.
The numerical evaluation criteria were applied to the fitted data for each of the columns. In addition, the parameter values were evaluated both for consistency and agreement with measured values.
|
RESULTS AND DISCUSSION
|
|---|
Chloride BTCs for all experiments are given in Fig. 3. The chloride concentration is plotted both as a function of time and as a function of cumulative outflow. In the first experiments for the two undisturbed cores, the same pattern emerged when plotted as a function of the cumulative outflow (Fig. 3a). For both flow rates, the chloride concentration in the drainage water rose rapidly during the first 2 to 3 cm after application, followed by a plateau in the chloride concentration (Fig. 3a). After approximately 6 cm of drainage, this was then followed by a more gradual increase in chloride concentration until the outflow chloride concentration was close to the input concentration. The difference between the flow rates was the chloride concentration of the plateau (Fig. 3a). The plateau for the high rainfall application rate was at 0.38 Cl/Cl0, and for the low rate at about 0.15 Cl/Cl0, where Cl is the concentration of chloride. On the basis of experiments done by Camobreco et al. (1996), we hypothesize that the early fast breakthrough and plateau were directly caused by the chloride breakthrough in the macropore network, while the concentration of chloride in the matrix flow was negligible. The second, more gradual increase in chloride concentration in the outflow water was the result of the breakthrough of chloride in the matrix. The BTCs plotted as a function of time (Fig. 3b) show that, despite the similarity when plotted as a function of the amount of water drained, the velocities for high and low rates were distinctly different.
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 3. Chloride breakthrough curves (BTCs) for undisturbed soil columns: (a) as a function of cumulative drainage, and (b) as a function of time (Exp. 1). The BTCs rising limbs for sand columns: (c) as a function of cumulative drainage, and (d) as a function of time. The BTCs falling limbs for sand columns: (e) as a function of cumulative drainage, and (f) as a function of time (Exp. 2). The BTCs for sand columns with chloride applied on the surface of the column: (g) as a function of cumulative drainage, and (h) as a function of time (Exp. 3). For (a) and (b), S represents the point where the rain with Cl was switched to regular acid rain. For (e) and (f), the y axis indicates the time after switching to the regular rain. C1 and C2 refer to Columns 1 and 2, respectively.
|
|
The BTCs for the two sand columns (the second experiments) are shown in Fig. 3c to 3f. The BTCs' rising limbs (Fig. 3c and 3d) for the first two rainfall intensity cycles (low and high application rates) coincided reasonably well when plotted as a function of time (Fig. 3d). The BTCs for the third cycle with a low application rate was delayed significantly. Conversely, the amount of water needed for initial breakthrough was directly dependent on the rate of application (Fig. 3c): the smallest amount for the 0.05 cm/h rate and the greatest for the 1.8 cm/h. The BTCs falling limbs (Fig. 3e and 3f) mirror the rising limbs in Fig. 3c and 3d. As before, the BTCs for the first two rainfall intensity cycles nearly coincided when plotted as a function of time (Fig. 3f) and varied with application rate when cumulative drainage was the independent variable (Fig. 3e).
These results differ from the standard CD equation where the BTCs are expected to nearly coincide when plotted as a function of the application rate, and greatly differ when charted against time. As we will see later, the different behavior of the BTCs in Fig. 3c and 3d is a direct consequence of the preferential flow paths in the conveyance zone. In sandy soils with fingered flow paths, the solute velocity through these fingers is dependent on the quotient of conductivity and equilibrium moisture content, and thus almost independent of the flow rate (Selker et al., 1992). The number of flow paths is dependent on the highest application rate experienced (Selker et al., 1996), that is, as many fingers form as needed to carry the flow. Thus, for the second low application rate, the wetted area was larger compared with the first application rate, resulting in lower chloride velocities (i.e., same flow rate over a larger area).
The results of the third set of experiments (Fig. 3g and 3h) where chloride was added to the surface of the column are, in many respects, the same as for the second set of experiments. The BTCs' rising limbs for the three intensity cycles coincided when plotted vs. time, and are independent of the flow rate. Note that in this set of experiments, different columns were used for each application rate, and thus each column was only subjected to one rainfall intensity rate and not several, as was the case for Exp. 1 and 2. The number of fingers increased, as expected, for increasing flow rates (Fig. 4; Table 3). As before, plotting the Cl concentration vs. the cumulative drainage separated the BTCs (Fig. 3g). Since the amount of chloride was fixed, the falling limbs were dependent on the initial breakthrough of the chloride.
View larger version (87K):
[in this window]
[in a new window]
|
Fig. 4. The number of fingers created at a depth of 14 cm with different flow rates. The left picture was taken after the rainfall rate of 0.3 cm/h during the third set of experiments (three fingers were delineated with a marker for convenience). The picture on the right is for the 1.8 cm/h of rainfall intensity in the second set of experiments. Under the higher rate of rain, fingers merged in the center of the column (dark area), but the outer edge of the column was still dry (bright area).
|
|
Model Results
In general, to calibrate the GPFM, three parameters are required for the fingered flow columns (Exp. 2 and 3). These are the amount of water in the distribution zone, the solute velocity, v, and the dispersion coefficient in the conveyance zone, D. On the basis of these parameters, the mobile region fraction, ß, can be calculated as:
 | [10] |
A value for = 0.15 cm3/cm3 was used as the equilibrium volumetric moisture content in the finger (Selker et al., 1992). In addition to these three parameters, the proportion of flux through the preferential flow path, ai, is required for the undisturbed soil columns in Exp. 1. The fitted and calculated model parameter values are given in Table 3, and the goodness of fit between observed and modeled data in Table 4. For both the R2 correlation and the MCE methods, unity indicates the best agreement between observed and predicted values, while the AIC is comparative indicator with the smaller AIC value indicating the better model.
View this table:
[in this window]
[in a new window]
|
Table 4. Statistical measures of generalized preferential flow model (GPFM) and convective-dispersive model (CD) performance for the three experiments.
|
|
For the undisturbed soil (Exp. 1 in Table 3), the proportion of flux through the preferential flow path was derived directly from the plateau of the preferential flow in Fig. 3a and 3b. For example, the relative concentration of 0.38 Cl/Cl0 was assumed to imply that 38% of the flux through the preferential flow paths at the input concentration and the remaining flux through the matrix with zero concentration. As long as the flux is below the saturated conductivity—as it was here—higher flow rates resulted in more pores taking part in the transport of the water and solutes. The velocity of water through the preferential flow paths was affected minimally by the number of pores transporting water and dissolved solutes, and was assumed independent of the application rate. A velocity 1.2 m/h for the preferential flow region fitted the data well. For the matrix flow, an eightfold increase in the flux increased the solute velocity by a factor of five to six times (Table 3). It is likely, therefore, that at a greater flow rate larger pores participated in the flow. The dispersivity, computed by dividing the dispersion coefficient by the pore water velocity, was between 1.5 and 2.5 cm and nearly independent of flux. The amount of water in the distribution layer increased from 1.8 cm at the low rate to 3 cm at the high rate. The fitted and observed BTCs for the high application rate for Column 1 are shown, as an example, in Fig. 5. Overall, the model produced a reasonable fit (Table 4) with realistic fitting parameters, although the increase in the apparent water contents of the distribution zone was initially not expected.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 5. Fitting results of the Generalized Preferential Flow Model (GPFM) and the Convective-Dispersive (CD) model for undisturbed soil Column 1 with the high rate (Exp. 1).
|
|
The BTCs from Fig. 3c to 3f for the second experiments (Exp. 2 in Table 4) were fitted to Eq. [5]. The fitted solute velocities in the conveyance zone for the first two rainfall intensity cycles were independent of application rates except for the very low application rate of 0.05 cm/h. This might have been caused by uneven water distribution associated with this low rate so that the whole distribution layer did not wet. In the third rainfall intensity cycle, the velocities were significantly lower because the high flux in the second cycle increased the number of flow paths which were maintained when the flow rate decreased. Hence, there was a greater proportion of the column wetted (Table 3). The apparent water contents in the distribution zone ranged from 0.6 to 1.5 cm (excluding the lowest application rate) and corresponded well to the observed distribution zone depth in the columns of 3 to 5 cm. A dispersivity value of 2 cm fitted all data equally well. Experimental observations and predicted curves are shown for Column 1 in Fig. 6. Generally, a good fit was obtained with R2 values of 0.97 or higher, with the exception of the 0.05 cm/h rate (Table 4).
View larger version (13K):
[in this window]
[in a new window]
|
Fig. 6. Fitting results of the Generalized Preferential Flow Model (GPFM) and the Convective-Dispersive (CD) model for sand Column 1 with the high rate (Exp. 2).
|
|
The third set of experimental BTCs were visually fitted to Eq. [4] (Exp. 3 in Table 3). The velocities were slightly less, and water contents of the distribution zone were greater than in Exp. 2. This might have been caused by the manure that plugged some of the pores. A good example was that in one of the columns under the high flow rate, the drainage rate decreased after application of the manure and then increased slowly but did not reach its initial value. As a result, some water ponded at the surface for this column. There was also a greater variability in the fitted parameters, especially for the dispersivity that varied from a value of 1.1 to 5. Example fitted and observed results are shown graphically for Column 1 with 1 cm/h rain in Fig. 7. Statistical analysis is given in Table 4.
View larger version (16K):
[in this window]
[in a new window]
|
Fig. 7. Fitting results of Generalized Preferential Flow Model (GPFM) and the Convective-Dispersive (CD) model for sand Column 1 with 1 cm/h rain (Exp. 3).
|
|
The amount of water in the distribution zone per unit area is plotted as a function of the flow rate in Fig. 8 for all three experiments. The water content values were consistent between the sand and the undisturbed soil columns and increased with increasing flow rates. The dependence of water contents on flow rate merits further investigation.
Comparison with the Convective-Dispersive Equation
Although the GPFM fitted the data reasonably well with parameter values that were physically meaningful, the question remains whether we could have obtained the same fit with the standard CD equation without the distribution zone. We fitted, therefore, all the BTCs to the standard CD equation (Eq. [3]) using the CXTFIT computer program (Parker and van Genuchten, 1984). The flux type boundary condition was applied with a retardation coefficient set to 1 and a decay rate set to 0. Assuming that the water was uniformly distributed through the columns, the moisture content was calculated by dividing the flow rate by the velocity. As shown in Table 5, for undisturbed columns the moisture content of the column decreased for the increased flow rate which is counterintuitive. Also, the dispersivity values for the high application rate were out of the expected range of 5 to 20 cm for field soils (Jury et al., 1991). Both models having approximately the same values for goodness of fit when both the R2 and the MCE tests are used (Table 4). The AIC, which penalize for the greater number of fitting parameters, indicate that the standard CD equation is the best model for the low application rate but not for the high application rate. Thus, for the undisturbed columns, despite that the standard CD model gave a statistically slightly better fit for the low application rate, the model parameters are not consistent for the two application rates and will give unrealistic results when used outside the calibration range.
For the sand columns in the second set of experiments, it was more difficult to judge from the parameter values in Table 5 whether the CD equation was valid. The R2, AIC, and MCE values were generally high (Table 4). The AIC clearly favors the standard CD model. The moisture contents were generally lower than would be expected. It is interesting to note when comparing the parameter values in Table 5 with Table 3 that the estimated velocities were lower and the dispersion coefficient higher than for the GPFM. Dispersivity values of 0.5 to 2 cm are expected for packed columns (Jury et al., 1991). Thus, again although the standard CD equation gave a statistically better fit, especially when the test penalized for the number of parameters used, it assumed that the water flowed uniformly through the subsoil, which was different than observed (Fig. 4).
The standard CD equation did not fit the third experimental data overall (Table 4). Especially the R2 and MAE values were low, indicating that the individual points fitted poorly. In addition, the AIC values were lower for the GPFM. Most of the parameters used were also physically unrealistic. The moisture contents were unreasonably low, similar to the second set of experiments. Dispersivities were also too large compared with the suggested experimental values (Jury et al., 1991). Thus, unlike the GPFM (Table 3), the standard CD equation was incapable of predicting solute transport applied on the soil surface with the feces. Since in Exp. 2 the CD equation fitted the observed data of a fingered flow experiment, the main culprit in Exp. 3 was likely the surface boundary condition. The surface boundary condition for the standard one-dimensional CD equation is a constant concentration which was imposed for Exp. 2 but was not met in Exp. 3.
The final question is if the two-region CD equation (Parker and van Genuchten, 1984) can predict the observed BTCs. For the sandy soils, the soil remained dry between the fingers, and the assumed exchange in the two-region CD equation between the matrix and the preferential flow paths could not take place. Thus, the sand columns will not be considered further. We also fitted the parameter values for the undisturbed soil column experiment. Although the best-fitting results were achieved with reasonable values for the velocity and dispersivity of 0.6 cm/h and 5.7 cm, respectively, with the high R2 value of 0.98, the dimensionless variables estimated, ß and 
, were unrealistic. The product ßR (where R is the retardation factor and set to 1 for chloride) of 0.2 is significantly higher than the observed data, and of 3732 is awkward compared with the values shown as examples (0.424–14.2) by Parker and van Genuchten (1984). Therefore, we concluded that although the predictions fitted the observed data well, the two-region CD equation is not satisfactory to predict the solute transport through preferential flow paths. Earlier we showed that the two-region CD equation and the GPFM can fit the BTCs for undisturbed columns equally well with parameter values that are physically acceptable (Steenhuis et al., 2001).
In summary, the GPFM is an alternative to predict BTCs in field soils where preferential flow is a common occurrence. Overall, it gives the best fit to the data and might, under some circumstances, provide a better physical representation of solute transport processes within the soil than either the standard two- or one-region CD equation. We do not reject the validity of the CD equation itself. We just interpret the parameters such as velocity and cross-section of flow differently. Moreover, the surface boundary condition has an exponentially changing solute concentration and not a constant concentration as is the case for the standard CD equation. The exact representation of flow patterns might prove critical when transport of adsorbed chemicals and colloids are considered.
Received for publication March 17, 2003.
|
REFERENCES
|
|---|
- Ahuja, L.R., D.G. Decoursey, B.B. Barnes, and K.W. Rojas. 1993. Characteristics of macropore transport studied with the ARS root zone water quality model. Trans. ASAE 36:369–380.
- Ahuja, L.R., and O.R. Lehman. 1983. The extent and nature of rainfall–soil interaction in the release of soluble chemicals to runoff. J. Environ. Qual. 12:34–40.[Abstract/Free Full Text]
- Akhtar, M.S., T.S. Steenhuis, P.A. Medrano, M. deGroot, and B.K. Richards. 2003. Dissolved phosphorus losses from undisturbed soil cores: Related to adsorption strength, flow rate, or soil structure? Soil Sci. Soc. Am. J. 67:458–470.[Abstract/Free Full Text]
- Andreini, M.S., and T.S. Steenhuis. 1990. Preferential paths of flow under conventional and conservation tillage. Geoderma 46:85–102.
- Bauters, T.W.J., D.A. DiCarlo, T.S. Steenhuis, and J.-Y. Parlange. 1998. Preferential flow in water-repellent sands. Soil Sci. Soc. Am. J. 62:1185–1190.[Abstract/Free Full Text]
- Beven, K., and P. Germann. 1982. Macropores and water flow in soils. Water Resour. Res. 18:1311–1325.
- Brush, C.F., W.C. Ghiorse, L.J. Anguish, J.-Y. Parlange, and H.G. Grimes. 1999. Transport of Cryptosporidium parvum oocysts through saturated columns. J. Environ. Qual. 28:809–815.[Abstract/Free Full Text]
- Camobreco, V.J., B.K. Richards, T.S. Steenhuis, J.H. Peverly, and M.B. McBride. 1996. Movement of heavy metals through undisturbed and homogenized soil columns. Soil Sci. 161:740–750.
- Dekker, L.W., and C.J. Ritsema. 1994. Fingered flow: The creator of sand columns in dune beach sands. Earth Surf. Processes Landforms 19:153–164.
- Faybishenko, B., C. Doughty, M. Steiger, J.C.S. Long, T.R. Wood, J.S. Jacobsen, J. Lore, and P.T. Zawislanski. 2000. Conceptual model of the geometry and physics of water flow in a fractured basalt vadose zone. Water Resour. Res. 36:3499–3520.
- Hendriks, R.F.A., K. Oostindie, and P. Hamminga. 1999. Simulation of bromide tracer and nitrogen transport in a cracked clay soil with the FLOCR/ANIMO model combination. J. Hydrol. (Amsterdam) 215:94–115.
- Hill, D.E., and J.-Y. Parlange. 1972. Wetting front instability in layered soils. Soil Sci. Soc. Am. Proc. 36:697–702.
- Hipel, K.W. 1981. Geophysical model discrimination using the akaike information criterion. IEEE Trans. Autom. Control 26(2):358–378.
- Huang, K., N. Toride, and M.Th. van Genuchten. 1995. Experimental investigation of solute transport in large, homogeneous and heterogeneous, saturated soil columns. Trans. Porous Media 18:283–302.
- Jarvis, N.J., P.E. Jansson, P.E. Dik, and I. Messing. 1991. Modeling water and solute transport in macroporous soil: 1. Model description and sensitivity analysis. J. Soil Sci. 42:59–70.
- Jury, W.A., W.R. Gardner, and W.H. Gardner. 1991. Soil physics. 5th ed. John Wiley & Sons, New York.
- Jury, W.A., and K. Roth. 1990. Transfer functions and solute movement through soil: Theory and applications. Birkhäuser, Verlag, Basel.
- Kung, K.-J.S. 1990. Preferential flow in a sandy vadose zone: 2. Mechanisms and implications. Geoderma 46:59–71.[ISI]
- Kung, K.-J.S., T.S. Steenhuis, E. Kladivko, T.J. Gish, G. Bubenzer, and C.S. Helling. 2000. Impact of preferential flow on the transport of adsorptive and non-adsorptive tracers. Soil Sci. Soc. Am. J. 64:1290–1296.[Abstract/Free Full Text]
- Larget. 2003. Applied regression analysis statistics lecture handouts: AIC and BIC. Class notes Statistics 333. Univ. of Wisconsin, Madison.
- Nissen, H.H., P. Moldrup, and R.G. Kachanoski. 2000. Time domain reflectometry measurements of solute transport across a soil layer boundary. Soil Sci. Soc. Am. J. 64:62–74.[Abstract/Free Full Text]
- Parker, J.C., and M.Th. van Genuchten. 1984. Determining transport parameters from laboratory and field tracer experiments. Bull. 84–3. Virginia Ag. Exp. Station, Blacksburg, VA.
- Perrin, C., C. Michel, and V. Andreassian. 2001. Does a large number of parameters enhance model performance? Comparative assessment of common catchment model structures on 429 catchments. J. Hydrol. (Amsterdam) 242:275–301.
- Quisenberry, V.L., and R.E. Phillips. 1976. Percolation of surface-applied water in the field. Soil Sci. Soc. Am. J. 40:484–489.[Abstract/Free Full Text]
- Ritsema, C.J., and L.W. Dekker. 1995. Distribution flow—A general process in the top layer of water repellent soils. Water Resour. Res. 31:1187–1200.
- Ritsema, C.J., L.W. Dekker, J.L. Nieber, and T.S. Steenhuis. 1998. Modeling and field evidence of finger formation and finger recurrence in a water repellent sandy soil. Water Resour. Res. 34:555–567.
- Russo, D., and W.A. Jury. 1987. A theoretical-study of the estimation of the correlation scale in spatially-variable fields: 1. Stationary fields. Water Resour. Res. 23:1257–1268.
- Selker, J.S., T.S. Steenhuis, and J.-Y. Parlange. 1992. Wetting front instability in homogeneous sandy soils under continuous infiltration. Soil Sci. Soc. Am. J. 56:1346–1350.[Abstract/Free Full Text]
- Selker, J.S., T.S. Steenhuis, and J.-Y. Parlange. 1996. An engineering approach to fingered vadose pollutant transport. Geoderma 70:197–206.
- Stagnitti, F., J.-Y. Parlange, T.S. Steenhuis, B. Nijssen, and D. Lockington. 1994. Modeling the migration of water soluble contaminants through preferred paths in the soils. p. 367–379. In K. Kovar and J. Soveri (ed.) Groundwater quality management. IAHS Publ. No. 220. Int. Assoc. of Hydrological Sci. Press, Oxfordshire, UK.
- Steenhuis, T.S., J. Boll, G. Shalit, and I.A. Merwin. 1994. A simple equation for predicting preferential flow solute concentration. J. Environ. Qual. 23:1058–1064.[Abstract/Free Full Text]
- Steenhuis, T.S., Y.-J. Kim, J.-Y. Parlange, M.S. Akhtar, B.K. Richards, K.-J.S. Kung, T.J. Gish, L.W. Dekker, C.J. Ritsema, and S.A. Aburime. 2001. An equation for describing solute transport in field soils with preferential flow paths. p. 137–140. In D.D. Bosch and K.W. King (ed.) Preferential flow, water movement and chemical transport in the environment. Proc. ASAE 2nd Int. Symp., Honolulu, HI. 3–5 Jan. 2001. Am. Soc. of Agricultural Engineers, St. Joseph, MI.
- Toride, N., F.J. Leij, and M. Th. van Genuchten. 1995. The CXTFIT code for estimating transport parameters from laboratory or field tracer experiments. v. 2.0. Res. Rep. 137. U.S. Salinity Lab., Riverside, CA.
- USEPA. 1992. Technical support document for land application of sewage sludge. Vol. I. Publ. 822/R-93–00a. USEPA, Washington, DC.
- van der Molen, W.H. 1956. Desalinisation of saline soils as a column process. Soil Sci. 81:19–27.
- Webster, R., and A.B. McBratney. 1989. On the akaike information criterion for choosing models for variograms of soil properties. J. Soil Sci. 40(3):493–496.
- Weglarczyk, S. 1998. The interdependence and applicability of some statistical quality measures for hydrological models. J. Hydrol. (Amsterdam) 206:98–103.
- Wierenga, P.J., and M.Th. van Genuchten. 1989. Solute transport through small and large unsaturated soil columns. Ground Water 27:35–42.
This article has been cited by other articles:
|
|
|
|
 
L. Luo, H. Lin, and P. Halleck
Quantifying Soil Structure and Preferential Flow in Intact Soil Using X-ray Computed Tomography
Soil Sci. Soc. Am. J.,
June 18, 2008;
72(4):
1058 - 1069.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
K.-J.S. Kung, E. J. Kladivko, C. S. Helling, T. J. Gish, T. S. Steenhuis, and D. B. Jaynes
Quantifying the Pore Size Spectrum of Macropore-Type Preferential Pathways under Transient Flow
Vadose Zone J.,
August 24, 2006;
5(3):
978 - 989.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
