| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Dep. of Environmental Sciences and Energy Research, Weizmann Institute of Science, 76100 Rehovot, Israel
* Corresponding author (brian.berkowitz{at}weizmann.ac.il)
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONCONTINUOUS TIME RANDOM WALK...MATERIALS AND METHODSCONCLUSIONSREFERENCES |
|---|
Abbreviations: ADE, advection-dispersion equation • BTC, breakthrough curve • CTRW, continuous time random walk
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONCONTINUOUS TIME RANDOM WALK...MATERIALS AND METHODSCONCLUSIONSREFERENCES |
|---|
In the usual one-dimensional form, the ADE is simply
 | [1] |
= D/v is usually referred to as the dispersivity (Bear, 1972), and is considered a property only of the porous medium microgeometry.
In the majority of cases, the actual BTCs do not appear in the published literature. Rather, only the final results of fitting the experimental data with the ADE—that is, the estimated values of D or —have been presented. As such, direct and independent checks of the quality of the fits are not possible. In fact, as early as 1959, it was already noticed that there were systematic errors in fitting BTCs using the classical ADE. In a series of careful and well-documented column experiments, Scheidegger (1959) observed that deviations in fits of the ADE to the BTCs could not be explained simply by the usual variability (error) in experimental measurements:
Thus, the experiments confirm the diffusivity equation up to this standard error. The deviations are systematic which appears to point toward an additional, hitherto unknown, effect.
Aronofsky and Heller (1957) also reported, in an analysis of published tracer experiments, that systematic deviations arise between measurements and predictions using the ADE.
Indeed, over the last two decades, several studies have pointed out specific and serious inadequacies in the applicability of the ADE, even at the laboratory (small) scale. Silliman and Simpson (1987), for example, demonstrate convincingly in laboratory experiments the scale dependency of the dispersivity coefficient; this is in stark contrast to the fundamental assumption that the dispersivity is a constant. It is also critical to note that the ADE assumes complete homogeneity of the porous medium, at least at the relevant scale of measurement. However, it has been shown that even carefully packed laboratory-scale flow cells and columns containing porous media are in fact "heterogeneous". For example, studies using magnetic resonance imaging to visualize flow conditions within "homogeneous" geologic materials in laboratory-scale column experiments report the existence of preferential flow paths, which strongly influence both water flow and tracer transport (e.g., Hoffman et al., 1996; Oswald et al., 1997). These paths occur due to macrostructures, caused for example from bridging effects and microstructures, reflecting grain-size heterogeneities. As a consequence, in contrast to BTCs predicted by the ADE, tracer test measurements tend to display anomalous early arrival times and long late time tails (Levy and Berkowitz, 2003). Our use of the term "anomalous" refers mainly to the fact that the tracer arrivals are consistently found in such experiments—at both early and late times—to be slower than those described by the (Fickian-based) ADE. This behavior is often referred to as non-Fickian; it is sometimes called pre-asymptotic or pre-ergodic. Thus even at the laboratory scale, serious inadequacies in the applicability of the ADE are evident.
And yet, in spite of the early (and entirely correct) observation of Scheidegger (1959) and these other studies of anomalous dispersion behavior, the ADE has remained the principal means by which dispersion in porous media is quantified. In this context, tracer transport in both fully and partially saturated porous media has been analyzed.
Only more recently have efforts to quantify anomalous (non-Fickian) behavior been undertaken. Stochastic perturbative approaches, employing an ensemble-averaged ADE, have been applied (as reviewed in Dagan and Neuman, 1997; Dagan, 1989; Gelhar, 1993) to quantify transport in mildly heterogeneous porous media. Development of these approaches, however, has been motivated primarily by non-Fickian transport behaviors observed in field experiments (e.g., Boggs et al., 1992; Adams and Gelhar, 1992) and in numerical transport simulations (e.g., Burr and Sudicky, 1994; Naff et al., 1998; Salandin and Fiorotto, 1998; Pannone and Kitanidis, 2001; McLaughlin and Ruan, 2001; Dentz et al., 2002). However, to the best of our knowledge, application of these methods in the context of (re)analyzing non-Fickian BTC measurements in laboratory-scale flow systems has not been considered. With regard to transport in heterogeneous media, an alternative (non-perturbative) probabilistic modeling framework is offered by CTRW theory. The CTRW theory has proven highly successful in quantifying non-Fickian transport in hydrogeological systems at a variety of scales (for a review see Berkowitz and Scher, 2001). In particular, transport behavior predicted by the CTRW theory is well suited to capturing the anomalous early and late time tailing observed so frequently in BTC measurements (Levy and Berkowitz, 2003).
Specific analysis of anomalous transport in homogeneous porous media has received little attention (Levy and Berkowitz, 2003). In light of recent developments in applying the CTRW theory, and the accompanying new conceptual understanding of transport behavior in porous media, we re-examine here some classical BTCs on Berea sandstones by Scheidegger (1959), on various kinds of Aiken clay loamy soils and Oakley sands by Nielsen and Biggar (1962), and on undisturbed saturated and unsaturated soils by Jardine et al. (1993). Originally, these BTCs were analyzed in terms of the classical ADE model to derive dispersion coefficients.
In the Continuous Time Random Walk Theory section below, we present a brief overview of the essential conceptual and quantitative aspects of CTRW theory, and some developments of specific relevance to analysis of the BTC experiments under consideration here. We then provide, in the Laboratory Experiments section below, a description of the experiments by Scheidegger (1959), Nielsen and Biggar (1962), and Jardine et al. (1993). In the Quantitative Analysis of Breakthrough Curve Data Using CTRW and ADE Models section below, we fit these historical BTC measurements with solutions from the CTRW theory, and contrast the fits with those provided by the ADE solution.
|
CONTINUOUS TIME RANDOM WALK THEORY |
|---|
|
TOPABSTRACTINTRODUCTIONCONTINUOUS TIME RANDOM WALK...MATERIALS AND METHODSCONCLUSIONSREFERENCES |
|---|
Berkowitz and Scher (1995)(2001) note that a broad distribution of characteristic time scales—notably power-law transport time distributions—can result from a broad distribution of heterogeneity length scales. This phenomenological picture is the basis of CTRW models. Significantly, forms of the ADE, as well as a variety of mobile-immobile and multirate mass transfer models (e.g., Carrera et al., 1998; Roth and Jury, 1993; Cunningham et al., 1997; Haggerty and Gorelick, 1995; Haggerty and Gorelick, 1998) can be derived as special cases within the CTRW framework (e.g., Berkowitz et al., 2002; Dentz et al., 2004). Fractional-in-time derivative formulations (e.g., Metzler and Klafter, 2000; Schumer et al., 2003) of the transport equation are special forms of the CTRW, and their application is limited to a narrow subset of transport behaviors and to large times (Berkowitz and Scher, 2001; Berkowitz et al., 2002; Scher et al., 2002). Fractional-in-space derivatives have also been considered in the literature (e.g., Pachepsky et al., 2000), though some questions regarding their applicability to soil science and hydrological problems have arisen recently (Lu et al., 2002; Zhou and Selim, 2003). Moreover, neither fractional derivative formulation (space or time) can allow transport to evolve from non-Fickian to Fickian over long spatial/temporal scales.
We present here a short summary of the conceptual picture associated with CTRW theory and the derivation of its formulation in terms of partial differential equations. Comprehensive explanations and a full accounting of the mathematical development associated with the CTRW theory have been presented elsewhere (Berkowitz and Scher, 1998; Margolin and Berkowitz, 2000, 2002; Berkowitz et al., 2002; Dentz and Berkowitz, 2003; Cortis et al., 2004). A practical users' guide to application of CTRW theory solutions employed here is given in Berkowitz et al. (2001) (freely accessible software is available at www.weizmann.ac.il/ESER/People/Brian/CTRW, verified 12 May 2004).
In the CTRW framework, contaminant migration in a variable velocity field is envisioned as particles executing a series of steps, or transitions, through the formation via different paths with spatially changing velocities. The sporadic interaction of particles in high and moderate velocity paths with low velocity regions often leads to non-Fickian transport behaviors, that is, non-Fickian transport can arise if the encounter-range relationship between particles and the velocities produces a wide spread of different sequences in the flow paths of migrating particles. This kind of transport can in general be represented by a joint probability density function, 
(s, t), which describes each particle transition over a distance and direction, s, in time, t. Of course, particle movement occurs along continuous paths; our definition of discrete transitions here refers to a conceptual discretization of these paths, which can be made at as high or as low a resolution as desired. Similar to any probabilistic/stochastic approach, we define (s, t) for an ensemble average over many possible realizations of the medium. As such, we assume here that the formation properties are stationary (i.e., statistical properties are the same at any location in the system), although the system itself is not homogeneous. The non-stationary case has been studied by Cortis et al. (2004).
We stress that, as we are introducing new information in the model on the microscopic geometry, new parameters are needed to precisely quantify this disorder. Also, the classical parameters of the ADE model (namely, velocity and dispersion) assume a new significance in light of the CTRW description.
Identification of (s, t) lies at the basis of the CTRW theory. In the following we will focus on a decoupled transition length and time (s, t):
 | [2] |
(t) denotes the marginal density of (s, t) (e.g., Berkowitz et al., 2002). The p(s) is assumed to be sharply peaked around its mean (i.e., it is assumed to have finite first and second moments). According to the Central Limit Theorem, after a sufficiently large number of steps the sum of the displacements will be Gaussian distributed. It is therefore reasonable to assume p(s) as Gaussian-distributed from the outset. Under this uncoupling assumption, anomalous behavior can arise with a power-law, transit-time (or event time) distribution:
 | [3] |
Many choices of the transition-time distribution (t) are possible (Cortis et al., 2004), depending on the nature of the porous medium. Here we use a (t) with a truncated power-law distribution (Dentz et al., 2004):
 | [4] |
2 
t2/t1 and 
(a, x) is the incomplete Gamma function (Abramowitz and Stegun, 1970). For transition times t1 << t << t2, (t) behaves according to a power law 
(t/t1)–1–ß, while for t >> t2, (t) decreases exponentially. This simple, yet general form of (t) allows for a power law (algebraic) decay (whose intermediate time behavior we can approximate as (t) 
t–1–ß, with the constant exponent 0 < ß < 2) to evolve to an exponential decay (characteristic of Fickian transport). The time for crossover from power law to exponential decay is governed by t2. The characteristic time t1 denotes a typical median transition time, so that the time range of interest is t > t1; below this cut-off, the system has not yet reached the ergodic limit and the upscaling is not meaningful. Further discussion of application of this (t) is given in Dentz et al. (2004). Using functional forms such as Eq. [3] or [4] (with t1 << t << t2), the very nature of the dispersion (e.g., Fickian or non-Fickian) can be characterized by the value of ß and falls into three possible ranges: ß > 2, 1 < ß < 2, and 0 < ß < 1 (Shlesinger, 1974; Dentz et al., 2004). For ß > 2, the CTRW model yields the classical Fickian behavior described by the ADE model. In this case, the tracer plume center of mass (mean location of the plume) travels at the average fluid velocity (and therefore scales as time t), while the standard deviation scales as t1/2. This case is equivalent to the ADE and thus the BTCs are Fickian in shape. For 1 < ß < 2, the mean of the tracer plume moves with a constant velocity. The standard deviation of the plume, however, scales as t(3–ß)/2. The curves are asymmetric with long late time tails and as ß increases the resulting BTCs become sharper and less disperse. For 0 < ß < 1, the BTCs display the most anomalous behavior. The curves are not symmetrical and delayed early and long late time tails exist. The mean and standard deviation each scale as tß, the shapes of the BTCs are functions of ß, and are similar on different spatial scales.
The Transport Equation
It is well known that the diffusion equation can be derived from the so-called Brownian motion (random walk) model. In this model the residence time of a particle in a particular point in space is independent of the position itself. This is a reasonable assumption for a completely homogeneous system such as a still fluid. Porous media are, however, heterogeneous systems even at the scale of the diffusion length, with regions of different local velocity. In this case the assumption of a uniform residence time breaks down. The pertinent model is therefore one in which the residence time of the walker at a given point in space is itself a random variable.
We note here that the CTRW formulation is particularly well suited to describe anomalous dispersion in unsaturated porous media. The complex maze of preferential channels and stagnant zones typical of such situations is described by the conceptual picture of different resident times of the tracer at different locations. As such a mass balance formulation in terms of a Master Equation is in order. The air phase simply creates a more tortuous pathway for the tracer as it migrates through the liquid phase. This can be incorporated directly in the form of (s, t). In fact, in Quantitative Analysis of Breakthrough Curve Data Using CTRW and ADE Models below, we will see that while the classical ADE model fails to describe the BTCs of unsaturated porous media, the CTRW formulation provides an excellent description of the early short times, and late large time arrivals typical of these situations.
In the following, it will be more expedient to work in the Laplace transformed space of the time variable t: we denote the Laplace variable by u and the Laplace transformed quantities by a ~.
In CTRW, the underlying mass balance equation for the (dimensionless) bulk concentration, cb, is given by a Generalized Master Equation
 | [5] |
(s, u) is a memory function that can be expressed in terms of 
(u) by
 | [6] |
The Laplace transform of the truncated power-law Eq. [4] is given by:
 | [7] |
As discussed in the Basic Concepts section above, we adopt the uncoupling assumption, (s, u) = p(s)(u)—Eq. [2]—and we require p(s) to have finite first and second moments. Further assuming smooth behavior of the concentration field, we can then obtain a Fokker-Planck equation with memory for the Laplace transformed concentration 
(s, u) = b(s, u)/n of the form (in one dimension)
 | [8] |
 | [9] |
is some characteristic time, and
 | [10] |
We renormalize the space variable x by the column length L = 1, such that v v/L, and D D/L2. By definition, the transport velocity v is generally different from the average pore velocity v, which enters the original form of the ADE model (Eq. [1]). It is convenient to define a dimensionless dispersivity, ,
 | [11] |
The definition of dispersivity in Eq. [11] is completely different from the classical definition of dispersivity as a characteristic pore length (e.g., Bear, 1972). In particular, need not necessarily be constant, but can vary, for instance, with the velocity field. It is important to recognize that the transport velocity, v, is distinct from the average pore velocity, v, whereas in the classical ADE picture, these velocities are identical. These two quantities do not scale linearly with each other, because the particle jumps are influenced differently than the fluid by the mobile–immobile zones. This implies that the first two moments of p(s) do not scale necessarily in a linear fashion as a function of the fluid velocity v. Also, the fact that the two distributions are uncoupled does not mean that the inherent non-linearities associated with the heterogeneities of the microscopic fluid flow do not affect the first two moments of the spatial distribution. We note that when 
(u) = 1, Eq. [8] reduces to a form analogous to Eq. [1].
We now define the flux-averaged concentration as 
= v ( – 
/x) (Kreft and Zuber, 1978). Given a unit steady-state flux boundary condition = u–1 at the inlet (x = 0) and a natural boundary condition /x = 0 at the outlet (x = 1), the exact analytical solution of Eq. [8] for (u) at x = 1 is
 | [12] |
 | [13] |
An excellent approximation of Eq. [12] for large Peclet numbers (say Pe –1 > 10) is
 | [14] |
Different boundary conditions can be implemented just as easily. Also note that our initial condition was c0(x) = 0: the solution in Eq. [12] can be modified easily to account for a different initial condition.
The expression for the flux-averaged concentration in Eq. [12] is explicit in Laplace space, but to obtain BTCs as a function of time t we need to numerically invert Eq. [12]. This is accomplished by means of the so-called de Hoog algorithm (de Hoog et al., 1982), which makes use of complex-valued u Laplace parameters. As a rule of thumb, the inversion of a complete BTC requires 10 to 20 solutions of Eq. [8] for each decade of time considered. The expression for in Eq. [12] must remain bounded for all values of the u variable, such that Re(u) > 
, some positive constant defining the axis of convergence in the complex u-plane. It is not possible to find analytically the zeros of the denominator. However, numerical simulations indicate that no singularities are present for Eq. [12] for the (u) in Eq. [7] (for all possible combinations of parameter values). Some analytical consideration can, in addition, be made for the approximate solution in Eq. [14], as that expression is in fact unbounded only when (u) = 0. It is seen easily that the memory function constructed from the (u) in Eq. [7] can never equal zero.
In Fig. 1, we plot a series of BTCs obtained from Eq. [12] for different values of the parameters ß and 2. We note that for > 2 the BTC decays quickly (in a Gaussian fashion, as for the ADE solution). On the other hand, power-law behavior is evident for late arrival times of the BTCs smaller than 2.
|
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONCONTINUOUS TIME RANDOM WALK...MATERIALS AND METHODSCONCLUSIONSREFERENCES |
|---|
Nielsen and Biggar (1961)(1962) subsequently performed a series of miscible displacement experiments in both saturated and unsaturated soils. The material used was an air-dry soil screened through a 2-mm sieve. The samples were packed in columns 30 cm long (the column area and the total porosity of the samples are not available from the experiment description). Nielsen and Biggar (1961) designed a special apparatus to keep stationary flow and constant water content in which the tracer-free water (0.005 M in CaSO4) was replaced by water containing the dissolved tracer (0.1 M CaCl2). The analysis of the Cl– concentration was made by titration with AgNO3. We revisit here three typical BTCs presented in Nielsen and Biggar (1962). Two BTCs correspond to experiments using columns filled with Aiken clay loam 0.23- to 0.50-mm aggregates, in saturated conditions, with imposed microscopic velocities of 3.40 and 0.058 cm h–1, respectively. The third BTC corresponds to a partially saturated (volumetric water content 
= 0.27) Oakley sand with an imposed microscopic velocity equal to 1.03 cm h–1.
More recently, (Jardine et al., 1993) presented results of miscible displacement experiments on undisturbed cylindrical soil columns (8.5-cm diam., 24-cm length). The soil columns were saturated with 0.05M CaCl2 from the bottom and were then allowed to drain. Bromide was used as the non-reactive, passive tracer. We re-examine here three of the typical BTCs obtained for different degrees of saturation. The BTC for a fully saturated column corresponds to a pressure head of h = 0 cm with = 0.549 (cm3 cm–3) and a pore-water velocity v = 14.66 cm h–1. Breakthrough curves for two partially saturated columns correspond to volumetric water contents of = 0.533 (h = –10 cm, v = 2.82 cm h–1) and = 0.513 (h = –15 cm, v = 0.354 cm h–1), respectively.
Quantitative Analysis of Breakthrough Curve Data Using CTRW and ADE Models
We now present and contrast the CTRW and ADE solution fits to the various measured BTCs described in the previous section. In all cases, best-fit curves (in the least squares sense, with all parameters fit simultaneously) were obtained from solutions of the CTRW formulation, using Eq. [8], and from the ADE [1] as a reference. A summary of the fitted parameter values for all of the BTCs is given in Table 1. Error estimates for parameters represent 95% confidence intervals. Also shown in Table 1 are the standard deviations of the fits.
|
|
|
|
|
|
|
|
We stress that the information contained in the early and late time regions of a BTC is critical, and should not be underestimated. Accurate characterization of early arrival time behavior is particularly significant for issues related to, for example, escape of contaminants from subsurface waste repositories, while late time tailing behavior is of key importance with respect to ground water remediation problems.
In their study of tracer transport in saturated and partially saturated columns, Nielsen and Biggar (1962) also reported systematic deviations of the calculated values using the ADE from the experimental data. Here, we apply the CTRW description to fit three of the measured BTCs, all of which display non-Fickian transport behavior. As before, the CTRW solutions capture the BTC behavior more completely than those provided by the ADE, as seen in Fig. 3, 4, and 5.
Analysis of the experiments on saturated and unsaturated soil columns in Jardine et al. (1993) further confirm the occurrence of anomalous transport. Once again, as seen in Fig. 6, 7, and 8, the CTRW solutions capture the full evolution of the BTCs. With respect to Fig. 2 through 8, it is noted that in three cases, the last two points on each graph are not matched by the CTRW solutions; although no error bars are available, it is suggested (due to the abrupt changes in slope) that this is due to measurement uncertainty.
From the analyses of the three independent sets of experiments, with parameter values from the fits shown in Fig. 2 through 8 summarized in Table 1, a number of comments are in order. First, as discussed above in the introduction and Continuous Time Random Walk Theory section above, the definition of dispersivity takes on different meanings in the CTRW and ADE formulations. Comparison of the dispersivity values for these two models for the BTCs in Fig. 3 and 4 provides a clear explanation for this difference in the definitions. In the ADE model the value of the dispersivity should remain constant as the pore water velocity changes, especially over such small-scale columns. However, decreasing the velocity from v = 3.4 cm h–1 (Fig. 3) to v = 0.058 cm h–1 (Fig. 4) yields a change in the ADE value of by a factor of more than four. The CTRW model, on the other hand, does not require that = D/v remain constant.
Next, with respect to application of the CTRW theory, we observe from Table 1 that the ß values range from 0.94 to 1.67. As discussed in the Continuous Time Random Walk Theory section above, these values are well below the range in which Fickian transport occurs (ß > 2). Moreover, for a given water flux, the disorder might be expected to increase (i.e., ß should decrease) with decreasing saturation. The ß values shown in Table 1 do not, in fact, follow this pattern; compare, for example, the results for the BTCs shown in Fig. 3 and 4 and for Fig. 6 to 8. We suggest that this is due to the complex interplay between the change in saturation and velocity, with the subsequent changes in paths of tracer transport. Recall that the values of the velocities in the experiments referred to in Fig. 3, 4, and 6 to 8 are different. Finally, we note also that the key parameter in Eq. [4] is the exponent ß. The t2 parameter was extremely large for all the BTCs and therefore plays no significant role in these experiments. The parameter t1 represents the characteristic time for the scaling of the equations and sets the lower limit from which a power-law behavior begins to occur. The underlying flow process is what ultimately governs the precise shape of the transition probability (t), whose validity must be assessed for each particular physical situation.
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONCONTINUOUS TIME RANDOM WALK...MATERIALS AND METHODSCONCLUSIONSREFERENCES |
|---|
The CTRW framework offers a viable and general means to model tracer transport, which is more accurate than descriptions obtained with the conventional ADE and other approaches. This general framework can guide the design of new experiments on tracer displacement. Systematic experimental studies and analyses on the effects of pore structure, fluid velocity, column length, and water saturation on BTCs will, in particular, help to shed additional light on the relationship between characteristics and specific parameters of the CTRW theory.
|
ACKNOWLEDGMENTS |
|---|
Received for publication September 25, 2003.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONCONTINUOUS TIME RANDOM WALK...MATERIALS AND METHODSCONCLUSIONSREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  Y.-q. Tao, X. Jiang, Y.-r. Bian, X.-l. Yang, and F. Wang Transport of Malathion in Homogeneous Soil Liquid Chromatographic Columns: Influence of Nonequilibrium Sorption Vadose Zone J., February 1, 2009; 8(1): 42 - 51. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
L. Bolshov, P. Kondratenko, K. Pruess, and V. Semenov Nonclassical Transport Processes in Geologic Media: Review of Field and Laboratory Observations and Basic Physical Concepts Vadose Zone J., November 1, 2008; 7(4): 1181 - 1190. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
F. Delay, P. Ackerer, and C. Danquigny Simulating Solute Transport in Porous or Fractured Formations Using Random Walk Particle Tracking: A Review Vadose Zone J., May 12, 2005; 4(2): 360 - 379. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| The SCI Journals | Agronomy Journal | Crop Science | |||
| Journal of Natural Resources and Life Sciences Education |
Vadose Zone Journal | ||||
| Journal of Plant Registrations | Journal of Environmental Quality |
The Plant Genome | |||
N–1. Error estimates represent 95% confidence intervals.
