SSSAJ Journal of Natural Resources and Life Sciences Education
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 Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (22)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Pachepsky, Y.
Right arrow Articles by Rawls, W.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Pachepsky, Y.
Right arrow Articles by Rawls, W.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Pachepsky, Y.
Right arrow Articles by Rawls, W.
Soil Science Society of America Journal 64:1234-1243 (2000)
© 2000 Soil Science Society of America

DIVISION S-1-SOIL PHYSICS

Simulating Scale-Dependent Solute Transport in Soils with the Fractional Advective–Dispersive Equation

Yakov Pachepskya, David Bensonb and Walter Rawlsa

a USDA-ARS Hydrology Lab., Bldg. 007, Rm.104 BARC-WEST, Beltsville, MD 20705 USA
b Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512-1095 USA

ypachepsky{at}hydrolab.arsusda.gov


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 Theory
 Materials and methods
 NOTES
 Results
 Discussion
 REFERENCES
 
Solute dispersivity defined from the classical advective–dispersive equation (ADE) was found to increase as the length of a soil column or the soil depth increased. The heterogeneity of soil is a physical reason for this scale dependence. Such transport can be described assuming that the random movement of solute particles belongs to the family of so-called Lévy motions. Recently a differential solute transport equation was derived for Lévy motions using fractional derivatives to describe advective dispersion. Our objective was to test applicability of the fractional ADE, or FADE, to solute transport in soils and to compare results of FADE and ADE applications. The one-dimensional FADE with symmetrical dispersion included two parameters: the fractional dispersion coefficient and the order of fractional differentiation {alpha} , 0 < {alpha} <= 2. The FADE reduces to the ADE when the parameter . Analytical solutions of the FADE and the ADE were fitted to the data from experiments on Cl- transport in sand, in structured clay soil, and in columns made of soil aggregates. The FADE simulated scale effects and tails on the breakthrough curves (BTCs) better than, or as well as, the ADE. The fractional dispersion coefficient did not depend on the distance. In the clay soil column, the parameter {alpha} did not change significantly when the flow rate changed provided the degree of saturation changed only slightly. With the FADE, the scale effects are reflected by the order of the fractional derivative, and the fractional dispersion coefficient needs to be found at only one scale.

Abbreviations: ADE, advective–dispersive equation • BTC, breakthrough curve • CTRW, continuous time random walks • FADE, fractional advective–dispersive equation • FPE, Fokker-Planck equation • MIM, mobile–immobile model


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
Theory
 Materials and methods
 NOTES
 Results
 Discussion
 REFERENCES
 
SOLUTE TRANSPORT in soil is a phenomenon of the utmost importance and of challenging complexity. The first mathematical models of solute transport in soils were proposed {approx}40 yr ago (Day and Forsyth, 1957; Gardner and Brooks, 1957) and the search for appropriate models continues.

The advective–dispersive equation was the first solute transport model widely used for soils. For one-dimensional steady state flow of a nonreactive solute in a homogeneous soil layer, the ADE is:

(1)

Here and below c is the solute concentration (M L-3), Ds is the solute dispersion coefficient (L2 T-1), v is the average pore-water velocity (L T-1), x is distance, and t is time.

Both continuum mechanics and statistical physics have been used to derive the ADE. Continuum mechanics uses the solute concentration as the variable under study. The ADE is derived using the mass conservation law of solute. It is assumed that (i) convective solute transport occurs in all parts of the pore space and (ii) the solute flux deviates from purely convective or piston flow according the solute dispersion hypothesis defining the dispersion term, which is proportional to the concentration gradient (Hillel, 1980).

With the statistical physics approach, the variable under study is the probability distribution function of individual solute particle motions in soil solution. The ADE is derived from the assumption that solute particles undergo Brownian motion B(t), which is an addition of successive increments that are independent, identically distributed random variables having finite variance. This derivation is based on the fact that the distribution of the sum of such increments is a normal distribution. A solute transport equation is obtained by considering an ensemble of particles. The Fokker-Planck equation (FPE) describes the changes in the probability of a particle to be found in a particular spatial location at a particular time. For Brownian motion in one dimension, this equation is (Bhatacharya and Gupta, 1990):

(2)
where P is the probability, v is the drift of the process or mean velocity vector, and the coefficient D is the variance growth rate of B(t). If the particles are released simultaneously and do not affect each other, then the probability of a particle being found in a particular spatial location at a particular time is interchangeable with concentration, and mass conservation is equivalent to requiring that the sum of probabilities equals unity. Then Eq. [2] transforms into ADE for a nonreactive solute.

As results of ADE applications accumulated, it appeared that the ADE could not satisfactorily describe several important features of solute transport in soils. First, the dispersivity (defined as the ratio Ds/v) tended to increase as the length of the soil column or the soil depth increased. Second, BTCs of nonreactive solutes had shapes different from those suggested by the ADE, so that the solute remained in soil longer or left soil faster than the ADE predicted.

Some published observations of the dependence of the solute dispersivity on scale are summarized in Fig. 1 . Data in this figure show results of column and field experiments in which parameters of Eq. [1] have been estimated from measurements made at several depths in a soil profile or soil column. Porro et al. (1993) used a 5-m-long soil column, and tritium as a tracer. Ellsworth et al. (1996) reported results of a field plot study in which solutes were extracted at 25- and 65-cm depths during an unsaturated steady-state flow experiment, and the plot was excavated to the 2-m depth after the experiment. Yasuda et al. (1994) described a field plot experiment in which solute samples were taken at five 0.3- to 1.5-m depths during a ponded infiltration experiment. Zhang et al. (1994) worked with a 12.5-m-long horizontal sandy soil column and collected electrical conductivity data at 50-cm increments along the column during a steady state saturated flow experiment with NaCl as a tracer. Toride et al. (1995) reported results of Shiozawa who used four-electrode sensors to measure NaCl concentrations at depths of 11, 17, and 23 cm in sand columns with saturated and unsaturated flow. Snow et al. (1994) sampled a Br- solution at four depths from a 1-m-long lysimeter with a layered soil that was sprinkler-irrigated for 60 d. Khan and Jury (1990) applied pulses of CaCl2 solution to loamy sand soil columns of 20-, 40-, and 80-cm length with three different flow rates. Jury (1988) summarized a field study in which a plot was irrigated every second day after application of KCl. Tracer and solute concentrations were monitored at depths from 0.3 to 3 m. Data in Fig. 1 show that the solute dispersivity as found from application of the ADE can be considered as a power function of the distance traveled by the solute. The companion Table 1 shows the slopes of regression lines in Fig. 1. The slopes vary widely, from small values of 0.2 to a relatively large value of 1.7. Data indicating an increase in the solute dispersivity as the sampling depth increases have also been published by Jaynes et al. (1988), Butters and Jury (1989), Hamlen and Kachanoski (1992), and Jaynes and Rice (1993).



View larger version (17K):
[in this window]
[in a new window]
  		Fig. 1 Observed dependencies of the solute dispersivity on distance in soils. Data sources: Porro et al., 1993 (open circle); Ellsworth et al., 1996 (open square); Yasuda et al., 1994 (open triangle); Zhang et al, 1994 (open inverted triangle); Toride et al., 1995 (open diamond); Snow et al., 1994 (filled circle); Khan and Jury, 1990 (filled square); Khan and Jury, 1990 (filled triangle); Jury, 1988 (filled inverted triangle)

 

View this table:
[in this window]
[in a new window]
  		Table 1 Slopes of the log(dispersivity)-log(distance) dependencies found from field plot and soil column experiments

 
The ADE predicts an ideal sigmoidal shape for BTCs from soil columns that are sufficiently long. van Genuchten and Wierenga (1976) pointed out the deviations from this shape in experimental BTCs. They used the term tail to describe the last part of a nonsigmoidal BTC that approaches its asymptotic value. Heavy tails of BTC, that is approaching the asymptotic value more slowly than Eq. [1] predicts, were observed by a number of authors (Nielsen et al., 1986; Vachaud et al., 1990).

Continuum mechanics addressed the problems with the ADE applications by replacing the hypothesis of the convective solute mobility in the whole pore space with the hypothesis of the existence of pore space zones with different mobilities of solutes. The physical reason for the scale dependence of solute transport parameters is the heterogeneity of soil. It is established in subsurface hydrology that the scale dependence of the dispersivity arises from the way in which an individual solute particle will gradually sample more and more of the velocity fluctuations associated with the heterogeneity of an aquifer (Mishra and Parker, 1990; Beven et al., 1993; Liu and Dane, 1996). As Ellsworth et al. (1996) remarked, theoretically for unsaturated soil water flow systems, scale-dependent dispersion would apply when the dominant mechanism of effective dispersion is the variation in local pore water velocity. van Genuchten and Wierenga (1976) developed the first model that explicitly described the variations in solute mobility and divided soil pore space into two zones, mobile and immobile. This mobile–immobile (MIM) model was remarkably successful in describing column experiments and explaining tails on BTCs (Nielsen et al., 1986). Moreover, whereas the ADE required the depth-dependent dispersivity to simulate the solute transport in a soil profile, the MIM model worked with a constant (or, at least, not distance-dependent) dispersivity. Determining parameters of this model is a nontrivial task because the parameters appear to be strongly correlated; the MIM model seems to overparameterize the solute transport, although it represents the actual distribution of the flow velocities in soil in a piecewise fashion with only two intervals, zero and non-zero (Beven et al., 1993; Koch and Fluhler, 1993; Bronswijk et al., 1995; Montas et al., 1997). Introducing more intervals is possible and has been explored by several research groups (Gwo et al., 1995; Steenhuis et al., 1990; Hutson and Wagenet, 1995). However, as more regions are used, correlation among parameters increases, and it is more difficult to estimate parameters reliably, as the problem becomes yet more overparameterized.

The occurrence of heavy-tailed BTCs and scaledependent dispersivity can be addressed in a statistical physics framework. Montroll and Weiss (1965) and Scher and Lax (1973) defined a broad framework for particle transport that can incorporate very general motions called continuous time random walks (CTRW). Each motion that a particle undertakes is described by the probability of moving a random distance in a random amount of time. The CTRW reduce to standard Brownian motion unless random increments in the particle travel time or distance are heavy-tailed, meaning that the probability decays like a power law for large values (see Mandelbrot, 1983). When only the individual jump sizes have a heavy-tailed probability density, excursions of particles converge to a Lévy motion (sometimes called Lévy flights; Shlesinger et al., 1982; Klafter et al., 1987; Compte, 1996; Benson, 1998). Particles undergoing Lévy motion can be simply characterized as behaving mostly like in Brownian motion except for occasional large jumps. These large jumps are clustered (Hughes et al., 1981). In chaotic or turbulent flows, this type of transport has been described as being trapped in relatively stagnant vortices for periods of time with occasional travel within "jets" of high velocity fluid (see Weeks et al., 1995). The Lévy motion predicts heavier tails in the BTC than those produced by Brownian motion. It also predicts the growth of the apparent variance faster than in Brownian motion. These features make Lévy motion an attractive generalization of Brownian motion when describing scale-dependent solute transport in soil (Pachepsky et al., 1997).

A differential solute transport equation derived for Lévy motions would facilitate applications in the same way as the differential ADE facilitated applications of the Brownian motion model. Recently, Zaslavsky (1994) suggested a procedure to derive such an equation using fractional derivatives that in effect account for the spatial memory of solute particles. Saichev and Zaslavsky (1997), Benson (1998), and Chaves (1998) modified Zaslavsky's procedure to explicitly describe all one-dimensional Lévy motions.

The objective of this study was to test the applicability of the FADE to solute transport in soils and to compare results of FADE and ADE applications.


   Theory
TOP
 ABSTRACT
 INTRODUCTION
 Theory
Materials and methods
 NOTES
 Results
 Discussion
 REFERENCES
 
The movement of solute particles in soils may not follow Brownian motion because high velocity regions in soils tend to be spatially continuous at all scales. A particle traveling faster (slower) than the mean is much more likely to do so over a large distance. One can say that the particles have a spatial memory, and this feature is absent in Brownian motion. Motions with the persistence in movements can be simulated with Lévy statistics (Bouchard, 1995; Weeks et al., 1995). Lévy motions produce particle excursions in which significant deviations from the mean are much more likely to occur than in Brownian movement.

Lévy's random variables arise from the generalization of the central limit theorem (Gnedenko and Kolmogorov, 1954). If X1, X2, ..., Xn are independent identically distributed random variables, then

(3)
where 0 < {alpha} <= 2 is the parameter of the Lévy distribution, µ is the shift (the mean for {alpha} > 1), {sigma} is the scale parameter (the standard deviation for ), the sign -> denotes convergence in probability, and Y is a standard Lévy random variable. If the increments Xn have finite variance, for example they obey the Gaussian or uniform distributions, then the exponent and the random variable Y in Eq. [3] has the Gaussian distribution. Then Eq. [3] expresses the classic central limit theorem (Feller, 1966). If the increments Xn have an infinite variance (e.g., they obey power law distributions with the exponent {alpha} between zero and two), then the generalization of the central limit theorem (Eq. [3]) is valid with {alpha} < 2. A complete set of assumptions about the properties of random variables in Eq. [3] can be found in (Feller, 1966, Ch. XVII).

The family of Lévy distributions is parameterized with the value {alpha} showing, in particular, how heavy are the tails of the distribution. The heavier the tail is, the higher the probability is of large deviations from average. Values of {alpha} are within the range 0 < {alpha} <= 2. The Gaussian distribution is a specific case of a Lévy distribution with . The standard symmetric Lévy probability distribution function F{alpha}(y) for parameter {alpha} can be computed as (McCulloch, 1996):

(4)
where {phi} is the integration variable, and C and U{alpha} are auxilary functions of {alpha}.

(5)

The integral in Eq. [4] can be evaluated numerically with standard integration procedures. The asymmetrical Lévy distributions have been also studied (Samorodnitsky and Taqqu, 1994).

Figure 2 shows several symmetrical probability distribution functions from the Lévy family computed using Eq. [4]. As parameter {alpha} decreases from two, Lévy distributions differ more and more from the Gaussian distribution as the tails of the distributions become heavier (Fig. 2). The probability density functions of the Lévy distributions with {alpha} < 2 decrease slower than Gaussian probability density as the deviation from average becomes large.



View larger version (15K):
[in this window]
[in a new window]
  		Fig. 2 Probability distribution functions for symmetric Lévy distributions with (a) linear scaling and (b) probability scaling. The curves differ in values of the distribution parameter {alpha}

 
The governing (Fokker-Planck) equation of a particle undergoing Lévy motion uses a fractional-order dispersion derivative that accounts for the heavy-tailed and scaling nature of the motion (Saichev and Zaslavsky, 1997; Chaves, 1998; Benson, 1998; Meerschaert et al., 1999). For an ensemble of particles, the probability of particle location is replaced by the solute concentration similar to the derivation of the ADE in Eq. [2]. An uncommon feature of the resultant solute transport equation is a fractional derivative in the dispersion term. The simplest form of the one-dimensional equation assumes symmetrical dispersion (Benson, 1998; Meerschaert et al., 1999):

(6)

Here Dsf is the fractional dispersion coefficient (L{alpha} T-1) and the superscript {alpha} is the order of fractional differentiation, 0 < {alpha} <= 2, corresponding with the parameter of the distributions shown in Fig. 2. Fractional derivatives are integro-differential operators defined as (Samko et al., 1993):

(7)


(8)

Here {alpha} is the order of the fractional derivative, {alpha} > -1, {Gamma} is the gamma function, m is the smallest integer number larger than {alpha}. If the value of {alpha} is actually an integer, then fractional derivatives reduce to well-known ordinary derivatives. For example, it is easy to see that, for , the right-hand side expressions in Eq. [7] and [8] will reduce to usual second derivatives. Then Eq. [6] will reduce to the common ADE given in Eq. [1]. This is expected since the Lévy distribution with is the Gaussian distribution, and the particle motion is Brownian. We will use the term FADE for the Eq. [6] hereafter, noting that only the dispersion term is fractional.

Equation [6] has analytical, or closed form, solutions obtained using the Fourier transform (Benson, 1998). When the solute concentration at the top of soil layer changes from 0 to , the dependence of the solute concentration on x and t relative to the average position of the displacing front vt is reasonably well described with the probability distribution function of the Lévy distribution:

(9)
where F{alpha}(y) is given by Eq. [4] and [5], and . The solution (Eq. [9]) appears similar to the solution of the ADE given by the probability distribution function of the normal distribution (Ogata and Banks, 1961):

(10)

Equation [6] can also be approximated using finite differences (Benson, 1998). The resulting system of the difference equations can be solved with various boundary and initial conditions.


   Materials and methods
TOP
 ABSTRACT
 INTRODUCTION
 Theory
 Materials and methods
NOTES
 Results
 Discussion
 REFERENCES
 
Experimental Data
Published Cl- BTCs were used to test the applicability of the fractional differential equation to the solute transport in soils. Shiozawa (cited in Toride et al., 1995) measured Cl- BTCs with four electrode EC sensors at depths of 11, 17, and 23 cm. The experiments included (i) leaching with solute-free water during unsaturated condition at volumetric water content , and (ii) continuous application of 0.01 M NaCl solution to an initially solute-free saturated sand at a . The data in digital form are included in the CXTFIT software and data package (Toride et al., 1995).

Dyson and White (1987) studied Cl- transport in structured clay soil (Evesham series; UK classification) irrigated under flow rates of 0.28 and 2.75 cm h-1. Soil cores 16.4 ± 1.5 cm long were irrigated from 16 evenly spaced hypodermic needles set above the soil surface. A steady-state near-saturated flow was created. Initial volumetric water content was 0.52 ± 0.07 cm3 cm-3, saturated water content was estimated as 0.67 ± 0.02 cm3 cm-3, and the steady-state water content in soil columns was 0.59 to 0.62 cm3 cm-3. Soil was preirrigated with 10 mM CaSO4 solution to reach the steady-state water flow, and the step input of CaCl2 was applied at the same intensity afterwards. The data points were obtained by digitizing graphs found in the aforementioned publication. The digitizing was made in triplicate. Coefficient of variation within the replications did not exceed 0.1%.

Biggar and Nielsen (1962) studied Cl- transport in 30-cm-long columns filled with aggregates of Aiken (fine, parasesquic, mesic Xeric Haplohumult) clay loam. Aggregate sizes were in the ranges of 0.25 to 0.5, 0.5 to 1., and 1.0 to 2.0 mm. The water content and volume flux were controlled using hanging water columns and negative air pressure imposed on porous plates. All BTCs were determined on columns at slight negative pressures (values not reported). The columns were leached with the tracer-free 0.005 M CaSO4 solution, and the step input of 0.05 M CaCl2 solution was applied at the same flow rate afterwards. Experiments were carried out with flow rates of 2.0 and 0.04 cm h-1. We digitized graphs to obtain data points.

Fitting and Comparison of ADE and FADE
The lack-of-fit mean square (Whitmore, 1991)

(11)
was used to estimate parameters Ds and v in the ADE, and parameters {alpha}, Dsf, and v in the FADE. For the originally solute-free columns, is defined as the ratio of the effluent concentration c to the influent concentration c0. For the columns originally containing the solute at concentration ci and leached with the solute-free water, is defined as the ratio (ci - c)/ci. The value of N is the number of observation points, K is the number of parameters to estimate, two for the ADE and three for the FADE. A version of the Marquardt-Levenberg algorithm published by van Genuchten (1981) was used to minimize the lack-of-fit square errors for the FADE and the ADE. The procedure allowed us also to estimate the standard errors of parameters along with their average values.1 The BTCs were simulated using Eq. [10] for ADE and Eq. [9] for FADE. Trapezoidal integration with 10000 nodes was used to calculate the integral in Eq. [4]. Value of the root mean square error, , was computed along with the mean square lack-of-fit.

Two statistical criteria were used to compare performance of the ADE and the FADE. The statistical significance of difference between lack-of-fit square values was assessed using Fisher's statistics, . The hypothesis of equality between s2r, FADE and s2r, ADE was rejected when the Fisher's statistics exceeded the critical value FN-2,N-3; the latter was computed using the FDF routine from the IMSL library (Microsoft, 1994).

To compare the models' performance using paired residuals, we applied the Williams and Kloot (1953) criterion. Let two models (the FADE and the ADE) be fitted to measured data and produce residuals Ri, FADE and Ri, ADE. To apply the criterion, two auxiliary values were formed and the slope {lambda} of the regression was estimated. The hypothesis of the similar models' performance could not be rejected if the value of {lambda} did not differ significantly from zero. The FADE (ADE) model performed better when {lambda} was significantly less (greater) than zero.

The Student's t test was used to study the significance of differences in values of estimated parameters. The 0.05 significance level was used in all statistical comparisons.


   Results
TOP
 ABSTRACT
 INTRODUCTION
 Theory
 Materials and methods
 NOTES
 Results
Discussion
 REFERENCES
 
Experiments with Sand Columns
Measured (Toride et al., 1995) and fitted BTCs are shown in Fig. 3 . Both ADE and FADE capture shapes of BTCs. However, inspection of the graphs shows that FADE better simulates the tails on BTCs in unsaturated sand. Statistics of the model performance shown in Table 2 support this observation. The root mean square errors are markedly lower for FADE than for ADE in unsaturated samples but almost the same for FADE and ADE in saturated samples. The ratio of lack-of-fit squares is greater than the critical value in unsaturated samples but lower than the critical value in saturated samples. The Williams-Kloot regression slope {lambda} is significantly less than zero, indicating significantly better performance of FADE than that of ADE in unsaturated samples.



View larger version (18K):
[in this window]
[in a new window]
  		Fig. 3 Comparison of measured (filled circles) and calculated with the fractional advective–dispersive equation (FADE) (solid line) and with the advective–dispersive equation (ADE) (dotted line) Cl- breakthrough curves in (a) unsaturated and (b) saturated sand. Measured values are from Toride et al. (1995)

 

View this table:
[in this window]
[in a new window]
  		Table 2 Estimated parameters and performance statistics for the advective–dispersive equation (ADE) and fractional advective dispersive equation (FADE) applied to data on Cl- breakthrough curves from soil columns

 
Estimated parameter values are also shown in Table 2. Estimates of the pore water velocity do not differ between FADE and ADE in the same data sets. The ADE dispersion coefficient Ds increases significantly with the length of the column in experiments with unsaturated sand. In the same experiments, the FADE dispersivity parameter Dsf does not change statistically significantly with the length, although there is a trend for the estimated average Dsf to decrease as the length grows. The differences between estimates of the FADE parameter {alpha} are not statistically significant for consecutive lengths.

In experiments with saturated sand, values of the parameter {alpha} are close to two, although they significantly differ from two. There is no significant difference between these values among columns of different lengths. Estimates of the FADE parameter Dsf are very close to estimates of the ADE parameter Ds. Both Dsf and Ds tend to decrease as column length increases.

We also fitted ADE and FADE to combined data for all column lengths (data not shown). For the data from the unsaturated columns, FADE appeared to be superior to ADE, the F ratio was significantly larger than the critical value, and the regression slope {lambda} was significantly less than zero. No significant difference was found between FADE and ADE performances using the data from the experiments with saturated columns.

Experiments with Structured Clay Soil
Measured (Dyson and White, 1987) and calculated BTCs are shown in Fig. 4 . Inspection of graphs suggests better FADE performance than that of ADE. The model performance statistics F and {lambda} shown in Table 2 support this observation for the low flow rate experiment with . In the high flow rate experiment, the value of F appears to be less then the critical value, and {lambda} does not significantly differ from zero. Therefore, the hypothesis of similarity in performance between FADE and ADE cannot be rejected in the high flow rate experiment.



View larger version (13K):
[in this window]
[in a new window]
  		Fig. 4 Comparison of measured (filled circles) and calculated with the fractional advective–dispersive equation (FADE) (line) and with the advective–dispersive equation (ADE) (dotted line) Cl- breakthrough curves in near-saturated clay soil; water flux was (a) 0.28 cm h-1 and (b) 2.75 cm h-1. Measured values are from Dyson and White (1987)

 
Estimated values of the parameter {alpha} do not differ significantly for the two flow rates and are similar to those in unsaturated sand. Both of the dispersion coefficients Ds in the ADE and Dsf in the FADE increase approximately 37 times as the flow rate increases 10 times.

Experiments with Soil Aggregates
Estimated parameter values are shown in Table 3 . The FADE dispersivity parameter Dsf increases as the aggregate size increases at both flux rates. The differences between estimates of the FADE parameter {alpha} are not statistically significant for three aggregate diameters at the same flow rate. Values of the parameter {alpha} are significantly different from 2.0 at the high flow rate of 2.0 cm h-1 (Table 3) and the FADE RMSE values are {approx}20% less than the corresponding ADE values. Values of {alpha} do not differ significantly from two at the low flow rate, and no significant difference was found between the FADE and the ADE performances at this rate.


View this table:
[in this window]
[in a new window]
  		Table 3 Parameters of the fractional convective–dispersive equation found from Cl- breakthrough curves obtained by Biggar and Nielsen (1062) in columns filled with soil aggregates

 
Using FADE Parameters Estimated at One Depth to Simulate Breakthrough Curves at Other Depths
The FADE dispersivity parameter Dsf and the FADE parameter {alpha} found from data at the 17-cm depth in experiments reported by Toride et al. (1995) were used to generate the BTCs at other depths where the observations were available. Results shown in Fig. 5 illustrate the ability of FADE to capture scale effects in solute dispersion. The RMSE values of the FADE were 0.094 and 0.104 for the BTCs at 11 and 23 cm, respectively. The ADE individually fitted to each of BTCs provided significantly lower accuracy with RMSE values of 0.0128 and 0.0187, respectively (Table 2).



View larger version (17K):
[in this window]
[in a new window]
  		Fig. 5 Using fractional advective–dispersive equation (FADE) parameters found from the breakthrough curve at 17 cm to simulate breakthrough curves at 11 and 23 cm in unsaturated sand. Lines show the simulations, symbols show measured values from Toride et al. (1995)

 
Consequences of Using a Wrong Model
We carried out a numerical experiment. Averaged FADE parameters found for unsaturated sand and FADE parameters found for the clay soil in the low velocity experiment were used to generate BTCs at 10-, 50-, and 100-cm depths using the solution of the FADE given by Eq. [9]. Then the solution (Eq. [10]) of the ADE was fitted to these BTCs. Results are shown in Fig. 6 . The ADE fits simulated BTCs with an accuracy that would be considered satisfactory if the BTCs had a realistic experimental variation. The fitted ADE dispersivity Ds increases with depth. The relative increase in Ds is larger for the column with a larger Dsf value: the fitted Ds doubles when and almost triples when . The dependencies of Ds on depth follow a power law in the form Ds{propto}(depth)b, where for the sand parameters and 0.37 for the soil parameters. This type of scaling was also shown in Fig. 1.



View larger version (25K):
[in this window]
[in a new window]
  		Fig. 6 Fitting advective–dispersive equation (ADE) to breakthrough curves generated with the fractional advective–dispersive equation (FADE). Optimized values of the ADE dispersivity Ds, (cm h-1) are shown for three depths x (cm); (a) FADE parameters for unsaturated sand, (b) FADE parameters for the clay structured soil

 

   Discussion
TOP
 ABSTRACT
 INTRODUCTION
 Theory
 Materials and methods
 NOTES
 Results
 Discussion
REFERENCES
 
The FADE model performed as expected, able to simulate tails on BTCs and column length effects better than the ADE model. However, these advantages were not obvious in saturated sand as compared with unsaturated sand and at high flow rates in clay soil as compared with slow flow. The relative differences between velocities in different parts of the saturated pore space in sand columns might not be high enough to obey a Lévy distribution with parameter {alpha} values significantly different from two. In such a case, the motion of the particles would be close to Brownian motion and the ADE would be a satisfactory model. Data on columns with soil aggregates (Table 3) illustrate the role of the convective transport in applicability of the FADE. The FADE worked better than the ADE when the flow rate was relatively high and transport was dominated by the convective dispersion. The ADE was a good model when the flow rate was low, the molecular diffusion was a dominating transport mechanism, and the relative differences between velocities in different parts of pore space were not substantial.

It seems to be important that the FADE includes the classical ADE as a specific case when the parameter {alpha} does not differ significantly from two. Some authors did not find a tangible dependence between the dispersivity and depth. Khan and Jury (1990) observed an increase of the dispersivity with depth with undisturbed soil but not in columns with repacked soil. In such case, the value of {alpha} needs to be set at two to simulate solute transport in repacked columns and to a value less than two to simulate solute transport in undisturbed soil. Then Eq. [6] remains unchanged, but its parameters change.

The advantage of using the FADE to describe solute transport in soils is the separation of scale effects from the values of the transport coefficients. The scale effects are reflected by the exponent of the fractional derivative, and the transport coefficients need to be found at only one scale as the example in Fig. 5 shows. Both the exponent and the transport coefficient reflect the structure and heterogeneity of the pore space. In the clay column, the transport coefficient exhibited the dependence of the flow rate, whereas the parameter {alpha} did not change significantly when the flow rate changed, provided the saturation degree changed only slightly. However, the range of scales within which a single value of the exponent {alpha} is applicable is currently unknown. Experimental setups limit the range of scales within which the scaling of the dispersivity shown in Fig. 1 and 3 is valid. Germann (1991) argued that the scale-dependence of solute transport should exist within and beyond the depth of a heterogeneous structured soil. A number of researchers found several distinct ranges of scaling for various soil properties. Avnir et al. (1985) found at least two ranges of radii with different surface fractal dimensions in studied soils. Potential sources for the dependence of fractal dimensions on pore radii were reported (Wu et al., 1993; Perfect et al., 1993). Pachepsky et al. (1995a) found three or four distinct scaling intervals for pore surfaces in the range of pore radii from 4 nm to 5 µm. Further testing and validating the FADE should shed light on factors affecting the range of scales where value of {alpha} can be set constant.

The BTCs calculated with the symmetrical version of the FADE are more or less symmetric (i.e., that the very flat early concentration increase is similar to the slow approach of the asymptotic value; see Fig. 3). A heavier tailing was observed where the BTCs were highly asymmetric and had a fast early concentration increase and a slow approach of the final value. This could be attributed to the asymmetric solute dispersion. It would be interesting to explore fitting this type of data with the asymmetrical FADE (Meerschaert et al., 1999):

(12)
that describes asymmetrical Lévy distributions and has an additional parameter ß explicitly accounting for the skewness of the distributions. This equation reduces to Eq. [6] when the parameter . Figure 7 shows the effect of the skewness parameter ß on the BTCs computed using the analytical solution of Eq. [12] derived in Benson (1998). Introducing non-zero ß values renders an opportunity to simulate asymmetrical dispersion.



View larger version (16K):
[in this window]
[in a new window]
  		Fig. 7 Breakthrough curves simulated using asymmetrical Lévy distributions. Values of parameters: 1, {alpha} = 2, ß = 0; 2, {alpha} = 1.3, ß = 0; 3, {alpha} = 1.3, ß = 0.9; 4, {alpha} = 1.3, ß = -0.9; {sigma} = (Dfs|cos ({pi}{alpha}/2)|)1/{alpha}

 
Observations of the scale-dependent solute dispersion in soils correspond well with theoretical predictions of non-Fickian mass transport in fractal media and accumulated evidence of the fractal nature of soil structure. The scale-dependent dispersivity has long been known for geologic materials. Neuman (1990, 1995) summarized results of several studies of solute dispersion in geologic materials and found solute dispersivity to increase with scale approximately as L3/2, where L is the apparent length scale. Wheatcraft and Tyler (1988) attributed scale-dependent dispersivities in aquifers to the fractal geometry of the pore space. Fractal properties of soil pore space are well documented (Giménez et al., 1998). However, non-fractal porous media can show scale dependence of the dispersivity as well, caused by the heterogeneous structure of the flow field. Soil heterogeneity in general presents a probable reason for the dispersivity scaling as shown in Fig. 1.

Assuming that the solute transport is governed by FADE leads to the scale dependence in the ADE dispersivity when the ADE is fitted to the FADE BTCs (Fig. 6). This dependence could be expressed with a power law Ds{propto}(depth)b similarly to those shown in Fig. 1. Values of the exponent b were around 0.3; that is smaller than many of values observed in Fig. 1 (Table 1). Such moderate dependence on scale happened because values of {alpha} used in simulations were close to 1.6. Additional simulations have revealed that the exponent b grows both with an increase in Dsf and a decrease in {alpha} (data not shown). Both increases in Dsf and decreases in {alpha} correspond with an increase in soil heterogeneity, and that causes an increase in fitted ADE dispersivity values.

The spatial arrangement of particles in soils has been successfully described using mass, volume, and surface fractal dimensions. However, these dimensions may not be relevant for estimating scale dependencies of solute transport in soils. Particles undergoing a Lévy motion are known to leave self-similar trails (Hughes et al., 1981). However, scaling properties of these trails must be related to scaling of pore connectivity, which does not have a direct relationship with mass or volume fractal dimensions. Lemaitre and Adler (1990) studied saturated viscous flows through two Menger sponges having similar fractal dimensions and found quite different water transport properties in these sponges. To estimate soil hydraulic conductivity, Crawford (1994) has suggested the use of both mass fractal dimensions and the spectral dimension. One method of estimating the spectral dimension was recently suggested by Anderson et al. (1996). It remains to be seen whether this dimension can be related to the parameter {alpha} in FADE.

The FADE has three parameters as compared with two in the ADE. This difference has been taken into account in the comparison of models in Table 2, because the F criterion employed to distinguish models uses root mean square pure errors (Eq. [11]) that correct the root mean square errors with respect to the number of parameters in the model. The ADE could have three parameters after an introduction of an empirical power law dependence of the dispersivity on the distance as shown in Fig. 1. Such an amended ADE probably could reproduce data from Fig. 1. However, a physically based solute transport model should be able to simulate more complex solute movement than just a unidirectional solute propagation from a stepwise input source. For example, it should be possible to simulate the solute transport with alternating infiltration, redistribution, and evaporation. Both the original ADE (Eq. [1]) and the FADE (Eq. [6]) are physically based and therefore are not unidirectional, contrary to the ADE amended with an empirical dependence of dispersivity on distance from a boundary. Admittedly, to simulate complex solute movement, both the FADE and the ADE will require a concept that captivates the dependency of their parameters on water content and velocity.

Data in Fig. 1 show that scale-dependent dispersivity was observed both in saturated and in unsaturated soils. The structural heterogeneity effect on the distribution of the particle velocities can be both larger and smaller as water contents decrease. The FADE was significantly better than the ADE when sand was unsaturated (Fig. 3 and Table 2). However, this result cannot be generalized. Seyfried and Rao (1987) observed solute transport with heavy tails on BTC curves in soil with capillary pressure greater than -0.1 kPa, whereas the ADE was sufficient at pressures less than -0.1 kPa. We realize that much more data are needed to evaluate the FADE as a solute transport model in unsaturated soils. The analytical solution (Eq. [9]) that we have assumes a step input, whereas the majority of published data on solute transport in unsaturated soils (Fig. 1) are obtained in experiments with the finite pulse input. The work to apply the FADE to the finite pulse input experiments is under way.

The scale-dependent solute dispersion can be an important phenomenon to consider in estimations of the fate of agricultural chemicals. As the dispersivity grows with depth, the spread of solute concentrations can be larger than that estimated from data within a soil profile. There may be a potential for chemicals to travel into and within the vadose zone faster than expected from their movement in the upper part of soil profiles. The FADE is a new promising tool to tackle this problem.


   ACKNOWLEDGMENTS
 
D.A. Benson acknowledges support from DOE-BES grant DE-FG03-98ER14885.


   NOTES
TOP
 ABSTRACT
 INTRODUCTION
 Theory
 Materials and methods
 NOTES
Results
 Discussion
 REFERENCES
 
1 The FORTRAN program to estimate parameters of the FADE model is available from the corresponding author upon request. Back

Received for publication April 22, 1999.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 Theory
 Materials and methods
 NOTES
 Results
 Discussion
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
J. Vanderborght and H. Vereecken
One-Dimensional Modeling of Transport in Soils with Depth-Dependent Dispersion, Sorption and Decay
Vadose Zone J., January 24, 2007; 6(1): 140 - 148.
[Abstract] [Full Text] [PDF]

Home page
 Soil Sci.Home page
A. Cortis and B. Berkowitz
Anomalous Transport in "Classical" Soil and Sand Columns
Soil Sci. Soc. Am. J., September 1, 2004; 68(5): 1539 - 1548.
[Abstract] [Full Text] [PDF]

Home page
 Soil Sci.Home page
L. Zhou and H. M. Selim
Application of the Fractional Advection-Dispersion Equation in Porous Media
Soil Sci. Soc. Am. J., July 1, 2003; 67(4): 1079 - 1084.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited