Analysis of Tracer Test Breakthrough Curves in Heterogeneous

Porous Media Using Continuous Time Random Walks



1. Introduction

Groundwater movement in naturally fractured and heterogeneous porous aquifers is highly complex, due to strongly varying velocity fields. A key problem is how to describe tracer and contaminant movement in such systems. Realistic quantification of this movement is complicated by the uncertainty in characterization of aquifer properties.

Transport is almost invariably treated by using some form of the advection-dispersion equation (ADE) (e.g., Bear, 1972), or by using related approaches, such as standard particle tracking random walk techniques, that are based on the same assumptions as the ADE (e.g., Kinzelbach, 1988). These treatments include deterministic and stochastic approaches (e.g., Dagan and Neuman, 1997). However, as demonstrated frequently in the literature, these approaches often fail to capture contaminant migration even in many "homogeneous" systems. This failure is most clearly evidenced by the finding of "scale-dependent dispersion": contrary to the fundamental assumptions underlying use of the ADE, dispersivity is not constant, and the very nature of the dispersive transport seems to change as a function of time or distance traveled by the contaminant (e.g., Gelhar et al., 1992). This scale-dependent behavior (also sometimes referred to as "pre-asymptotic", "anomalous" or "non-Gaussian") is what we shall refer to as "non-Fickian" transport. Other evidence of non-Fickian behavior lies in the often observed "unusually" early breakthrough times, and "unusually" long late time tails, in measured breakthrough curves.

The generally accepted explanation for non-Fickian transport is that heterogeneities which cannot be ignored are present at all scales. One approach to address this difficulty is to attempt to resolve the hydraulic conductivity (or velocity) field at a sufficiently high level and apply a numerical code that incorporates the ADE (using either partial differential equation or particle tracking formulations) at the scale of these blocks. However, even highly discretized systems (e.g., with block sizes of the order of 10 cubic meters in large aquifers) have not adequately captured the migration patterns (e.g., Eggleston and Rojstaczer, 1998), suggesting that unresolved heterogeneities also exist at these relatively small scales. Maybe somewhat surprisingly, non-Fickian transport has been observed even in small-scale, relatively homogeneous, laboratory-scale models (Berkowitz et al., 2000). We note, parenthetically, that the all-too-frequent "patch" solution, involving use of a "functional form" for dispersion which allows the dispersivity to change with travel distance or time, is mathematically incorrect, and contradicts the fundamental assumptions used to develop the ADE (e.g., Berkowitz and Scher, 1995).

The ADE formulation assumes Fickian transport. In simple terms, this formulation assumes that the center of mass of the contaminant plume advances with the average ("macroscopic") fluid velocity, while the contaminant spread about this center of mass (due to both mechanical and chemical effects) is a purely Fickian process. While these assumptions may seem to be a reasonable working hypothesis, and the ADE has in some cases provided reasonable first-order estimates of how a contaminant plume will migrate, the underlying conceptual picture and mathematical framework is in fact valid only under highly restrictive conditions. Details of these conditions, which in simple terms require a very high degree of homogeneity in the hydraulic conductivity, can be found in Berkowitz and Scher (2000).

Thus, in summary, in spite of at least four decades of intensive research efforts, we lack proper tools to quantitatively capture contaminant plume migration behavior. Specifically, we still require an easy applicable framework to deal with non-Fickian transport in realistic, heterogeneous geological formations.

A general, physically-based approach to quantify transport is based on continuous time random walk (CTRW) theory. It was first applied to calculate impurity conduction in semiconductors (Scher and Lax, 1973a,b) and to analyze properties of amorphous semiconductors (Montroll and Scher, 1973; Scher and Montroll, 1975). An excellent pedagogical treatment of these analyses is given in Scher et al. (1991). In the context of geological materials, CTRW theory was developed and applied to numerical studies of transport in fracture networks (Berkowitz and Scher, 1995, 1997, 1998), to consideration of the oft-cited MADE tracer experiment (Berkowitz and Scher, 1998), to the analysis of a tracer test in fractured till (Kosakowski et al., 2000), and to modeling tracer transport in laboratory flow cells containing porous media (Hatano and Hatano, 1998; Berkowitz et al., 2000). Berkowitz and Scher (2000) recently demonstrated how the CTRW framework accounts for a very wide range of non-Fickian and Fickian transport behaviors, and how the ADE (as well as fractional derivative transport equations) can be derived from it under specific, well-defined conditions. Margolin and Berkowitz (2000) have further extended the theory.

The references mentioned above focus largely on theoretical and mathematical aspects of the CTRW, considering the underlying conceptual picture of transport in heterogeneous systems, and lay the groundwork for widespread application of the theory. In this contribution we demonstrate how the relatively complex mathematical framework of the CTRW approach can be practically implemented to examine tracer test measurements and analyze the transport characteristics of fractured and heterogeneous porous media. In the next section, we briefly introduce the CTRW formalism and discuss its physical meaning. We then (Section 3) present analytical expressions for tracer breakthrough curves (also referred to as "first passage time distributions"), for the full range of possible transport behaviors. Section 4 describes how to implement these solutions, and discusses an accompanying software library we have developed which can be downloaded freely from a website. We note that if desired, the reader can skip the mathematical details in Section 3 and move directly to Section 4. Example applications of these solutions to the analysis of measurements from laboratory and field tracer tests are provided.



2. The Continuous Time Random Walk (CTRW) Approach

We review here briefly the conceptual picture and mathematical language associated with the CTRW approach. Extensive discussion of the entire theory can be found in the references cited in the previous section. A major advantage of the CTRW theory is that the resulting solutions are robust over a wide range of cases, and require a minimum of fitting parameters.

We consider the movement of water and (conservative) chemical species (i.e., tracers) moving through a geological formation. Clearly, heterogeneities occur on a broad range of scales in most geological formations. These heterogeneities can consist of fractures (joints and/or faults), variations in the rock matrix (e.g., grain sizes, mineralogy, layering, and lithology), and/or large-scale geological structures. Under an external pressure gradient, the velocity and flux distributions are determined by the liquid properties and by the structure of the aquifer heterogeneities. Tracer particles (representing the contaminant mass) transported within the water move through the formation via different paths with spatially changing velocities. Different paths are traversed by different numbers of particles. Typically, heterogeneous systems show a broader distribution of velocities than homogeneous systems.

This kind of transport can in general be represented by a joint probability density function,, 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. By coupling particle migration in space and time, such a function naturally accounts for particle transitions that extend over short and long distances, and over short and long times. Similar to any probabilistic/stochastic approach, we definefor an ensemble average over many possible realizations of the medium. As such, we are assuming 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.

Identification oflies at the heart of the CTRW theory. It can be shown that the principal characteristics of tracer plume migration patterns are dominated by the behavior of at "large" (or "asymptotic") times. In the context of the CTRW theory, it turns out that "large" time is in practice rather small, and is for all intents and purposes reached almost immediately. Further discussion of this point is given in Berkowitz and Scher (1998). Simple asymptotic forms ofwhich can exist include exponential decay and power law (algebraic) decay. It can be shown that exponential decay leads to Fickian transport (Margolin and Berkowitz, 2000). More interesting is the power law decay, whose long time behavior we can approximate as , with the (constant) exponent> 0.

The exponentcontrols the particle migration behavior, and thus functionally quantifies the dispersion behavior. We stress, however, that the parameteris fundamentally more general than the dispersivity parameter in the ADE. 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:

We stress that the behavior of the ratio of the plume standard deviation to the mean is critical. For all cases where < 2, fitting the central region of the breakthrough curves (as is usually done, where the relative concentration is 0.5) with the ADE model yields poor results, particularly in the early and late time regions of the curves. Moreover, such erroneous fitting leads to the often observed, apparent increase in the effective dispersivity with time (or distance), as discussed in the Introduction.



3. First Passage Time Distribution (FPTD) Solutions

3.1 Background

Tracer test measurements often consist of one-dimensional (averaged) tracer concentration breakthrough curves, as a function of time t, at selected distances from the tracer source. The breakthrough usually refers to the plane of exiting particles, and the (non-cumulative) curve corresponds therefore to what is known historically as a "first passage time distribution" (FPTD). In terms of the CTRW formulation, the FPTD is defined as the probability per time for a tracer particle to reach the measurement plane at time t for the first time. Averaging over the coordinates parallel to the exit plane yields the average distribution of tracer per unit time. Thus, the solutions presented herein are applicable to breakthrough curves in time, i.e., to contaminant distributions measured at fixed distances downstream of the point of tracer injection. We note that this probabilistic formulation of tracer distribution is equivalent to that of flux-averaged concentration values.

As already stated in the previous sections, the solutions presented below appear, in various forms, in the literature cited in Section 1, and we refer the interested reader to these sources for the detailed derivations. Here we discuss aspects relating to the practical application of these solutions to the analysis of tracer test measurements.

3.2 FPTD solutions for 0 << 1

Anomalous transport for the case 0 << 1 was first investigated by Montroll and Scher (1973) and by Scher and Montroll (1975). They derived a solution for the average distribution of tracer per unit time, the FPTD function, for a pulse injection. Following Montroll and Scher (1973),is given by

(1)

where L is the distance between the inlet (origin) and the outlet (measurement) planes, is a dimensionless time,denotes the inverse Laplace transform, and u is the (dimensionless) Laplace variable. The constant b is defined here aswhereis the average length of a single particle transition andis a dimensionless constant defined at a fixed.

Implementation of this solution uses the formulas from Appendix C in Scher and Montroll (1975). The exact inverse Laplace transform solution of equation 1 is

(2)

where. Numerically we can evaluate equation 2 only for small to moderate values of x, because of the j exponent. For large x, we can approximate equation 1 by (Scher and Montroll, 1975)

.

(3)

Often, it is more convenient to examine breakthrough data in terms of cumulative FPTD (CFPTD) curves. As discussed in Berkowitz et al. (2000), the CFPTD curves for a pulse input are equivalent to the (non-cumulative) curves for a step input. We derive in Appendix A the following expressions for the CFPTD for 0 << 1.

Similar to the FPTD solutions, the same two cases introduced by x arise. The exact solution, which can be evaluated up to moderate values of x, is

(4)

while the approximation for large x (>>1) is

.

(5)

The crossover point between the two different solutions for the FPTD and CFPTD curves depends upon the value of. By estimating visually the optimal crossover point,, for different(or, alternatively, by defining two criteria at the crossover point: the two solutions, and their slopes, must coincide), we interpolate between these points and find that for 0.9 << 0.97,

(6)

while for 0.1 << 0.9,

.

(7)

If x is smaller than, we use the exact solution (equation 2 or 4); otherwise we use the solution for large x (equation 3 or 5).

Shown in Figure 1 are a series of typical FPTD and CFPTD curves for a range of values (0 << 1). From Figure 1, it is evident that asincreases, the cumulative breakthrough curves (CFPTD solutions) become sharper and less disperse. Note that the curves are not symmetrical, and that the long tails persist. Because the mean and standard deviation of the plume distribution both scale as, the shapes of the FPTD and CFPTD curves are functions only of, and are similar on different spatial scales.

The temporal x-axis scale in Figure 1 is given in dimensionless units of. (With respect to the FPTD, division ofbyon the x-axis scale requires multiplication of the y-axis scale by, in order to preserve normalization for.) In order to allow comparison of the FPTD and CFPTD curves to measurements, a multiplicative "shift factor",, is used to translate between dimensionless temporal units (from theory) and dimensional temporal units (from measurements). This factor, which scales the time, contains parameters that correspond to a characteristic tracer particle velocity(distinct from the average fluid velocity) for a given length scale of the medium. In the case 0 << 1,is always larger than the average fluid velocity, and we can write (e.g., Berkowitz and Scher, 1998). The derivation ofis shown in Appendix B, and its actual application to data analysis is discussed below in Section 4.3.





Figure 1: A series of typical FPTD and CFPTD curves for a range ofvalues (0 << 1). Shown here are semilog plots for a range of FPTD curves () and CFPTD curves () versus, for a range of.  (a) versus ,= 0.33, 0.50, 0.65, (b)versus,= 0.75, 0.80, 0.90 (c) versus ,= 0.33, 0.50, 0.65, and (d) M(L,t) versus,= 0.75, 0.80, 0.90 (after Berkowitz et al., 2000).



3.3 FPTD solutions for 1 << 2

The FPTD solution for 1 << 2 is similar to that for the case 0 << 1. Margolin and Berkowitz (2000) start from a simplification of the transition time density function. Considering the leading terms of the Laplace transform of this density function, denoted , they obtain for small u

(8)

wheredenotes a positive constant andrepresents a mean transition time for . Details of the solution derivations can be found in Margolin and Berkowitz (2000); here we provide only the solutions themselves. We note here that these particular solutions assume no tracer transport upstream of the point of injection; this approximation is applicable for many practical situations, and certainly valid at sufficient measurement distances from the input source.

We must first introduce a dimensionless timeand the dimensionless parametersandwith. The exact solution for the FPTD is

.

(9)

For large positive h we approximate the FPTD by

(10)

and for large negative h we have

.

(11)

Note that here we must consider three functional forms for the FPTD solution over the complete range of h, whereas only two functional forms arise in the FPTD solutions for 0 << 1.

As in the previous section, we consider also the cumulative FPTD (CFPTD) solutions. The derivation of the CFPTD solutions is shown in Appendix C. The exact solution is

.

(12)

For large positive h we use

(13)

and for large negative h

.

(14)

Similar to the 0 << 1 solutions, here the crossover points between the (three) different solutions for the FPTD and CFPTD curves depend upon the value of. We estimate visually the optimal crossover points for different, and interpolate between these points. The crossover points depend on h. We define

,

(15)

for the FPTD and CFPTD solutions, respectively, and

(16)

for both cases. Ifwe use equations 10 and 13 for the FPTD and CFPTD solutions, respectively. For equations 11 and 14 are used, and in all other cases the exact solutions 9 and 12 are applied.

The resulting FPTD and CFPTD curves for 1 << 2 are similar to those shown in Figure 1: the cumulative breakthrough curves (CFPTD solutions) become sharper and less disperse asincreases, and again, the curves are asymmetric with long tails. However, in contrast to the Figure 1 curves, and because the mean of the transition time (i.e., the first moment of) is finite, the shapes of the FPTD and CFPTD curves are functions of the actual spatial scale, as well as functions of .

In order to compare the FPTD and CFPTD curves to measurements, we must translate between the dimensionless variable h and the dimensional temporal unit, t, for breakthrough measurements at a distance L. Because of the dependence of the solutions on the spatial scale L, and as shown in Appendix D, two parameters must be determined to permit this translation. One of these parameters, denoted below by, is analogous to the multiplicative shift factor,, defined for the case 0 << 1, while the second parameter depends on. The derivation of these translation parameters is shown in Appendix D, and their actual application to data analysis is discussed below in Section 4.4.



3.4 Convolution solutions for transport through multiple layers and variable

input boundary conditions

As discussed in Section 2, the transport solutions presented in the previous sections allow determination and prediction of breakthrough curves at specified distances from the inlet boundary. These solutions are valid as long as the overall behavior of the flow field does not change, so thatremains constant. Moreover, the solutions have been developed for the case of either a pulse or step function tracer input condition. Clearly, in many experimental systems, either the input boundary condition is different (e.g., a decaying step input), or transport occurs through regions with distinctly different properties (e.g., through two or more sedimentary layers). Conveniently, convolution techniques can be used in a straightforward manner to deal with these situations.

We define first an input function, F(t), and a response function for the medium G(t), which corresponds to the FPTD solution for a pulse input. The convolution of F(t) with G(t), H(t), is defined formally as

(17)

where H(t) is the resulting breakthrough curve. This flexible and convenient formulation allows us to consider both of the situations outlined above. In the case of transport through two characteristically different layers, we can analyze tracer transport through each of the layers individually, defining two differentvalues (and thus two FPTD solutions, F(t) with G(t). The overall breakthrough curve that develops after transport through the two layers is given by the convolution of these two solutions. For the second case of a variable inlet boundary condition, G(t) represents the FPTD solution with a specificvalue (assuming a pulse inlet boundary condition), while F(t) represents the functional form of the tracer input. Applications of convolution methods in this context have been presented by Kosakowski et al. (2000), and are demonstrated in Section 4.6.


4. Application of the FPTD Solutions

4.1 Basic approach

We now describe how the solutions given in Section 3 can be applied to analyze measurements from laboratory and field tracer tests. These solutions, strictly valid for an infinite flow domain, can be used for both characterization and prediction of tracer migration. Input data for the analyses consist of breakthrough curve measurements, in the form of pairs of concentration versus time values at (one or more) given distances from the tracer inflow to the domain. The concentration measurements may be either cumulative or non-cumulative. We emphasize that we consider here averaged, one-dimensional concentration measurements, corresponding to one-dimensional or non-diverging/non-converging, averaged, two- or three-dimensional flow fields, and line (or areal, in three dimensions) tracer input sources. In addition, at sufficiently large measurement distances, the theoretical solutions can be applied to cases of point sources in two- and three-dimensional systems.

The FPTD solutions can be evaluated numerically with relative ease. We have developed a series of subroutines written in the C programming language, which are evaluated within the GRACE graphical analysis package. The solutions, along with experimental data which are to be analyzed, can then be plotted, and iterative fitting or straightforward comparison of the solutions to the data points can be achieved.

Regardless of the analysis to be performed (i.e., fitting with FPTD or ADE solutions), two important issues should be kept in mind. First, it is critical that reliable, high-resolution (flux-averaged) concentration data be available, especially at early and late times. Otherwise, it is difficult to properly distinguish the nature of the transport, whether it be Fickian or non-Fickian, and the appropriate value of(for further discussion, see Berkowitz and Scher, 2000). Second, to further aid proper analysis, the availability of two sets of breakthrough measurements, at two different distances from the tracer source, are useful (e.g., to check variations in parameter fits). Although these two requirements represent an "ideal case", their importance should not be under-estimated.

A critical aspect of any model is the number of fitting parameters used in its application. The ADE model in fact has two explicit fitting parameters - the average velocity, and the dispersivity (or, alternatively, the coefficient of dispersion). While the average velocity may not seem to be a fitting parameter, but rather a "given" of the system, we emphasize that imposition of the ADE model on a transport process assumes Fickian transport, which leads to the average tracer velocity being identical to the average fluid velocity.

The FPTD solutions for 0 << 1 involve two explicit fitting parameters:, which characterizes the dispersive process, and, which characterizes the effective tracer velocity at a given distance (i.e., the effective, average arrival time of tracer particles). In contrast, the FPTD solutions for 1 << 2 require three fitting parameters:,and. The latter two parameters are contained in the dimensionless parameter h, which is the variable of the FPTD solution (compare, e.g., equations 2 and 9). Note that in the case 0 << 1, the parameteris "absorbed" in. Although it is moot, we could argue that the ADE solution (which corresponds to the FPTD solutions for> 2) also requires three fitting parameters (,and) if we account for the implicit choice= 2.



4.2 Use of the GRACE graphical analysis package

The codes we have developed to evaluate the FPTD solutions are currently implemented as a shared library of functions for use with GRACE (http://plasma-gate.weizmann.ac.il/Grace). However, if desired, the codes can be run independently or in conjunction with other graphics software. The GRACE graphical analysis package produces a powerful array of two-dimensional data plots, working within a WYSIWYG user environment. As noted in the GRACE user manual (http://plasma-gate.weizmann.ac.il/Grace), GRACE "...runs under various (if not all) flavors of Unix with X11 and Motif (LessTif). It also runs under VMS, OS/2, and Windows (95/98/NT). Its capabilities are roughly similar to GUI-based programs like Sigmaplot or Microcal Origin plus script-based tools like Gnuplot or Genplot. Its strength lies in the fact that it combines the convenience of a graphical user interface with the power of a scripting language which enables it to do sophisticated calculations or perform automated tasks. ... GRACE can access external functions present in either system or third-party shared libraries or modules specially compiled for use with GRACE".

The ability of GRACE to support extended automated fitting of data and advanced visualization, as well as the portability across different OS platforms, are good reasons for using GRACE instead of programming our own environment. In order to work with the CTRW library as it is presented in this paper, it is necessary to download either the source file or the binary library for the relevant platform. Once downloaded, the file must be uncompressed and untarred; installation instructions appear in the README file. After successful installation of GRACE, four additional functions specific to the FPTD solutions must be downloaded (as described below and on the web homepage).

The software library can be downloaded freely from this webpage. A web link and instructions to access and download the GRACE graphical analysis package are also provided.It is of course not necessary to access this website in order to apply the FPTD solutions presented above (and the interested reader can program the solutions independently), but we note that the program names and implementation described in the following sections correspond to those in the website software library.


4.3 Application of FPTD solutions for 0 << 1

The software library contains two programs. The first program, ctrw01(xshift, beta, time), provides the FPTD solution for a pulse source (solving equations 2 and 3, and using equations 6 and 7), while the second program, ctrw01i(x_shift, beta, time), provides the FPTD solution for a step injection (equations 4 and 5, which is identical to the cumulative (integrated) FPTD, or CFPTD, for the pulse solution). The variable names "x_shift" and "beta" refer, respectively, tomeasurement in the data set being considered. It should be observed that because of computational sensitivities in the current numerical implementation, these solutions are limited to the range 0.1 << 0.97. However, our analyses to date of realistic formations have not found dispersive transport characterized byvalues less than about 0.3.

The nonlinear curve fitting option in GRACE can be used to fit the two parametersandto the experimental data. While the fitting routine is in this case relatively robust, we caution that nonlinear curve fitting routines often have difficulty differentiating among multiple local minima, and can be sensitive to the initial estimates parameter values. As such, it is important to test carefully the optimal parameter values returned by the routine, by using a variety of initial parameter estimates. For the CFPTD calculations, a convenient initial estimate ofis given by the time at which the relative concentration is about 0.5. Alternatively, trial-and-error variation of these parameters can be used to obtain a suitable fit. The programs return a value of the FPTD and CFPTD for each time specified in the input file. Once the parameter estimates (,) are available, for the tracer breakthrough at some distance L from the plane of tracer injection, predictions for the FPTD and the CFPTD can be made for the breakthrough behavior at any distancefrom this inlet by calling the program with the parameters (,) and a specified time.

To illustrate the application of these solutions, we consider briefly two examples of breakthrough curve analyses. One breakthrough curve consists of measurements in a laboratory flow cell (2.13 m in length) containing a uniformly heterogeneous packing of sand (Silliman and Simpson, 1987). Tracer (chloride) was introduced as a step concentration at the flow cell inlet, and concentration measurements were made downstream. As shown in Figure 2, we find (using the nonlinear curve fitting option) estimates ofandwhich yield a CFPTD curve that fits well the measured concentration breakthrough. Full details of the analysis are presented in Berkowitz et al. (2000).



Figure 2: Comparison of CFPTD solution (solid line), with= 0.87 and= 390, to measured concentrations (points) at a distance L = 1.37 m from the inlet (data from Silliman and Simpson, 1987).



The second breakthrough curve was measured in a field tracer experiment in a fractured till (Sidle et al., 1998). A 4 m ´ 4.8 m areal section was isolated vertically, and samplers were installed at depths of 2.5 m and 4 m. The breakthrough of chloride introduced as a step input concentration at ground surface was measured at these two depths. Figure 3 shows the fitted CFPTD curve (using the nonlinear curve fitting option) to a measured concentration breakthrough at the 2.5 m depth. Full details of the analysis are presented in Kosakowski et al. (2000).



Figure 3: Comparison of CFPTD solution (solid line), with= 0.61 and= 0.203, to measured chloride breakthrough data in a sampler at a depth of 2.5 m (sampler C5, data from Sidle et al., 1998). Note the logarithmic time scale.



We stress that, as discussed by Berkowitz et al. (2000) and Kosakowski et al. (2000), neither of the measured breakthrough curves shown in Figures 2 and 3 could be fitted adequately by using the standard ADE solution.



4.4 Application of FPTD solutions for 1 << 2

As for the case 0 << 1, the software library contains two programs. These programs, ctrw12(t_mean, beta, b_beta, time) and ctrw12i(t_mean, beta, b_beta, time), provide the FPTD solution for a pulse source, and the FPTD solution for a step injection (identical to the cumulative (integrated) FPTD, or CFPTD, for the pulse solution), solving equations 9-11 and 12-14, respectively, using terms 15 and 16 to determine the crossover points. The variable names "t_mean", "beta" and "b_beta" refer, respectively, to,, and, while "time" is a (dimensional) time corresponding to an actual experimental measurement in the data set being considered. Because of computational sensitivities in the current numerical implementation, these solutions are limited to the range 1.03 << 2. Note that for this range of, as well as for> 2, the mean velocity of the tracer plume is essentially that of the average fluid velocity, so thatis given by L divided the mean fluid velocity; this approximation is increasingly good asincreases from unity.

The nonlinear curve fitting option in GRACE can be used to fit the three parameters, andto the experimental data. Again, we caution that the existence of multiple local minima necessitates careful testing of optimal parameter values returned by the routine. In many cases, a good strategy is to obtain initial estimates ofandusing a fixed value of and then to fit all three free parameters. For the CFPTD calculations, a convenient initial estimate ofis given by the time at which the relative concentration is about 0.5. Alternatively, trial-and-error variation of these parameters can be used to obtain a suitable fit. The programs return a value of the FPTD and CFPTD for each time specified in the input file. Once the parameter estimates (,and) are available, for the tracer breakthrough at some distance L from the plane of tracer injection, predictions for the FPTD and the CFPTD can be made for the breakthrough behavior at any distancefrom this inlet by calling the program with the parameters (,,) and a specified time. Finally, as noted in Appendix C, the solutions given here are valid only for. This condition will be met in virtually all practical cases of interest, as long as sufficient breakthrough data are available at early and long times.

Example applications of these codes are provided in the next section, where we discuss solutions for both> 2 and for the ADE. We have, in the limited number of data sets analyzed to date, not found any real data sets indicating dispersive behavior characterized by 1 << 2.

4.5 Application of FPTD solutions for> 2

As shown from theoretical considerations (e.g., Margolin and Berkowitz, 2000), Fickian transport arises when³ 2. In these situations, the classical ADE can be applied. For completeness here, and to enable full data analysis within the software library that is provided in this paper, the solutions presented in Section 4.4 for 1 << 2 can also be applied using the value= 2. All solutions for³ 2 are captured by using the solution for= 2, because of an equivalence in the mathematical formulation (Margolin and Berkowitz, 2000). The program parameters are given exactly as described in Section 4.4, with= 2.

If the match between the FPTD curve and the breakthrough measurements is satisfactory, then it is recommended to continue the analysis using the standard ADE solution. For convenience, this solution is also provided in the website software library. However, we strongly caution the reader that erroneous application of the ADE is frequently made. Careful examination of the early and late time portions of the breakthrough curves is often necessary in order to distinguish between Fickian (ADE) and non-Fickian (CTRW) transport behavior.

To illustrate the application of these solutions, we consider an artificially generated set of breakthrough curve data, based on perturbed (random noise) values calculated using the ADE solution with an average velocity of 3.34 m/d and a dispersivityof 0.068 m, at a distance L = 3.39 m from the inlet. Shown in Figure 4 are the concentration points, the fitted solution using the FPTD with= 2, and the ADE solution. Clearly, the solutions are very similar, and capture the dispersive behavior demonstrated by the data points. The optimal estimate ofis 1.01, which is identical to L divided by the average fluid velocity (recall Appendix D). Similarly, the optimal value ofis 0.02; it can be shown thatis related directly to the dispersivity defined in the ADE,.



Figure 4: Comparison of CFPTD solution (solid line), with= 2, and the ADE solution (dashed line), to artificially generated breakthrough concentrations. The CFPTD and ADE solutions are essentially identical, and capture the dispersive behavior demonstrated by the data points. As discussed in the text, there is excellent agreement between the CFPTD and ADE parameter values.



4.6 Convolution solutions for transport through distinct heterogeneity layers/regions

As discussed in Section 3.4, convolution of the solutions given above can be used to deal with experimental systems in which either the input boundary condition is something other than pulse or step input, or transport occurs through regions with distinctly different properties.

Clearly, completely general, "automated", software analyses such as developed and explained in Sections 4.3-4.5 cannot be provided in such cases; more direct user involvement is required. The software library contains the basic programs to perform the convolution calculations, once the relevant FPTD solutions, F(t) and G(t), and correspondingvalues have been specified.

We illustrate the application of these solutions by considering again measured breakthrough curves in the fractured till field tracer experiment of Sidle et al. (1998); as noted previously, the dispersive behavior displayed in these breakthrough curves could not be captured by fitting the standard ADE solution. Sidle et al. (1998) suggested that a change in heterogeneity patterns at the field site might be present, leading to tracer migration through two layers with different dispersive properties. Quantification of transport through such a system can be treated easily by use of the convolution solutions.

Recalling the analysis for Figure 3, we use an estimate of the breakthrough curve at a depth of 2.5 m, based on a fittedvalue of= 0.61 (for sampler C5), as the input function to the lower layer. Clearly, this input function, which we store in the software library input data file, conv_in.dat, is distinct from the usual pulse and step input functions.

The software library contains two programs. These programs, conv01(x_shift, beta, time) and conv12(t_mean, beta, b_beta, time), can then be used to fit the breakthrough curve in the lower layer. The solutions used in the programs are, of course, based on those given in the previous sections; the variable names are defined previously, and the same limitations on implementation are relevant here. The resulting convolution solution, shown in Figure 5, is obtained using the nonlinear curve fitting option. It is clear that, in this example, both the single region and the two-region, convolution solutions provide excellent fits to the data. Thus, geological and other considerations are necessary to select the most appropriate conceptualization. Full details of the analysis are presented in Kosakowski et al. (2000).

Figure 5: Illustration of convolution solution application. Points indicate measured chloride concentrations, for a sampler at the 2.5 depth (sampler C5; triangles), and a sampler at the 4 m depth (sampler F5; circles) (data from Sidle et al., 1998). Parameter estimates for the individual breakthrough curves are= 0.61,= 0.203 (for sampler C5; long dashes) and= 0.5,= 0.724 (for sampler F5; short dashes). The convolution solution using the FPTD, shown as a solid line, is obtained by using the calculated breakthrough curve for the upper layer (= 0.61) as the input function to the lower layer. Then the fitted parameters for the lower layer are= 0.41,= 0.158, where the subscripts 1 and 2 denote the upper and lower layers, respectively. Note the logarithmic time scale.



5. Discussion and Concluding Remarks

The CTRW theory captures a broad range of dispersive transport behaviors. It can be applied to quantify the ubiquitous, non-Fickian tracer migration patterns not accounted for by the ADE, as well as Fickian (ADE) transport. Mathematically, the ADE can be derived as a special, limiting case of the CTRW. Importantly, the underlying conceptual and physical picture of tracer migration is realistic and intuitively (and mathematically) satisfying. Moreover, the resulting CTRW solutions are parsimonious, requiring a minimum number of fitting parameters.

Theparameter characterizes the functional nature of the dispersion occurring in the system. As a consequence, in fractured and heterogeneous porous media, this parameter may be expected to vary as a function of the relative degree of heterogeneity. For example, as the size of the heterogeneities becomes extremely small relative to the length scale of interest over which transport occurs, transport will tend to Fickian behavior (> 2). However, as is well-known from field studies, Fickian transport is seldom reached, even at extremely large length scales. In our analyses here, we assume thatvaries only very slowly with the length scale of interest, so that for all practical purposescan be considered constant. Our analyses to date support this assumption, and its precise limits are currently under detailed investigation.

The FPTD solutions presented here can be readily used to predict breakthrough curves for a different values of L, once the relevant parameters have been obtained by fitting the solutions to one (or preferably two) breakthrough curves. However, it is important to keep in mind that for very large L values, different scales of heterogeneity may be encountered in the field situation, which would imply a change in the value of. This caution is of course relevant also to application of the ADE (and the dispersion coefficient), when transport is Fickian.

To date, we have considered the movement of water and (conservative) chemical species migrating through (fully water-saturated) geological formations. We note that these solutions do not explicitly account for contaminant retardation (adsorption, desorption), decay and production processes. Similar non-Fickian transport can in principle be expected in some partially saturated systems, and in reactive transport systems. Moreover, solutions which apply to (1) the spatial distribution of contaminant concentrations, (2) temporal and spatial solutions valid for higher dimensional flow fields and for radially convergent/divergent flow fields, and (3) point sources, remain to be fully derived. These and other features are currently under development, and will be incorporated into future reports and versions of this software library.



References

Bear, J., 1972. Dynamic of Fluids in Porous Media, Dover Publications, New York.

Berkowitz, B., and H. Scher, 1995. On characterization of anomalous dispersion in porous and fractured media. Water Resour. Res. 31. no. 6: 1461-1466.

Berkowitz, B., and H. Scher, 1997. Anomalous transport in random fracture networks. Phys. Rev. Lett. 79, no. 20: 4038-4041.

Berkowitz, B., and H. Scher, 1998. Theory of anomalous chemical transport in fracture networks. Phys. Rev. E 57, no. 5: 5858-5869.

Berkowitz, B., H. Scher, and S.E. Silliman, 2000. Anomalous transport in laboratory- scale, heterogeneous porous media. Water Resour. Res. 36, no. 1: 149-158 (with a minor correction, published in Water Resour. Res. 36, no. 5: 1371).

Berkowitz, B., and H. Scher, 2000. The role of probabilistic approaches to transport theory in heterogeneous media., Transp. Porous Media, in press.

Dagan, G., and S. P. Neuman (Eds.), 1997. Subsurface Flow and Transport. A Stochastic Approach, Cambridge University Press, New York.

Eggleston, J., and S. Rojstaczer, 1998. Identification of large-scale hydraulic conductivity trends and the influence of trends on contaminant transport. Water Resour. Res. 34, no. 9, 2155-2168.

Gelhar, L.W., C. Welty, and K.R. Rehfeldt, 1992. A critical review of data on field-scale dispersion in aquifers. Water Resour. Res. 28, no. 7: 1955-1974.

Hatano, Y., and N. Hatano, 1998. Dispersive transport of ions in column experiments: An explanation of long-tailed profiles. Water Resour. Res. 34, no. 5: 1027-1033.

Kinzelbach, W., 1988. The random walk method in pollutant transport simulation, in Groundwater Flow and Quality Modelling, NATO ASI Ser., Ser. C Math and Phys. Sci., vol. 224, edited by E. Custodio, A. Gurgui, and J. P. Lobo Ferreria, pp. 227-246, D. Reidel, Norwell, Mass.

Kosakowski, G., B. Berkowitz, and H. Scher, 2000. Analysis of field observations of tracer transport in a fractured till. J. Contam. Hydrol., in press.

Margolin, G., and B. Berkowitz, 2000. Application of continuous time random walks to transport in porous media. J. Phys. Chem. B 104, no. 16: 3942-3947 (with a minor correction, published in J. Phys. Chem. B 104, no. 36: 8762).

Montroll, E.W., and H. Scher, 1973. Random walks on lattices. IV. Continuous-time walks and influence of absorbing boundaries. J. Stat. Phys. 9, no. 2: 101-135.

Scher, H., and M. Lax, 1973a. Stochastic transport in a disordered solid, I. Theory. Phys. Rev. B 7, no. 10: 4491-4502.

Scher, H., and M. Lax, 1973b. Stochastic transport in a disordered solid, II. Impurity conduction. Phys. Rev. B 7, no. 10: 4502-4519.

Scher, H., and E.W. Montroll, 1975. Anomalous transit-time dispersion in amorphous solids. Phys. Rev. B 12, no. 6: 2455-2477.

Scher, H., M.F. Shlesinger, and J.T. Bendler, 1991. Time-scale invariance in transport and relaxation. Physics Today, January, 26-34.

Sidle, C.R., B. Nilsson, M. Hansen, and J. Fredericia, 1998. Spatially varying hydraulic and solute transport characteristics of a fractured till determined by field tracer tests, Funen, Denmark. Water Resour. Res. 34, no. 10: 2515-2527.

Silliman, S.E. and E.S. Simpson, 1987. Laboratory evidence of the scale effect in dispersion of solutes in porous media. Water Resour. Res. 23, no. 8: 1667-1673.



PRINCIPAL SYMBOLS

Greek:

: exponent which functionally quantifies the dispersion behavior

: factor for tracer breakthrough measurement length

: standard deviation of plume location

: dimensionless time

: joint probability density function which describes each particle

"transition" over a distance and direction, s, in time, t

English:

b : dimensionless constant, defined as

: dimensionless constant for a fixed transition length

: first passage time distribution (FPTD) function for a pulse injection

h : dimensionless parameter, defined as

j : integer

: inverse Laplace transform

L : distance between the inlet (origin) and the outlet (measurement) planes

l : dimensionless parameter, defined as

: average length of a single particle transition

: mean plume location

s : distance and direction of a single particle transition

t : time

: characteristic tracer velocity

x : dimensionless constant, defined as





Appendix A: Derivation of the CFPTD for 0 << 1

We derive here expressions for the CFPTD for 0 << 1. Because these CFPTD solutions have not been published previously, we present here the full derivations. There are in fact two possible ways (at least) to derive the CFPTD. One method involves integrating the already known solutions for the (non-cumulative) FPTD. A second method is now shown. It is known that

(A1)

where c is a real, positive constant. Then

.

(A.2)

With the change of variables(A.2) becomes

(A.3)

where. In the last equality in (A.2), we used the fact that

(A.4)

which also means that the CFPTD can be calculated as

.

(A.5)

In other words, no particles arrive at zero time or before it (as is physically correct). This result does not hold for the case 1 << 2; see Appendix C).

Using the identity

(A.6)

we have that

(A.7)

where here c represents a contour integral and it is known that

.

(A.8)

Since for j ¹ 0,

,

(A.9)

equation 4 follows immediately. This is an exact solution. For relatively small times, when x >> 1, the method of steepest descent gives from (A.3) the approximate solution 5.



Appendix B: Derivation offor 0 << 1

A multiplicative "shift factor",, is used to translate between dimensionless units in the FPTD and CFPTD solutions and the dimensional temporal units of actual measurements. As seen in equation 2, the basic independent variable is. We define (Berkowitz et al., 2000; Margolin and Berkowitz, 2000)

(B.1)

where t is the dimensional time,and. Using the parameter(Margolin and Berkowitz, 2000), we have from (B.1) that

(B.2)

We can write this expression somewhat more conveniently. The asymptotic form, and the accompanying CTRW formulation, are such that the mathematics (e.g., the result) lead to expressions that can be shown to be independent of the choice of(Margolin and Berkowitz, 2000). Then choosingand defining vL as the effective, average tracer velocity over the length scale L, equation (B.2) can be simplified andcan be defined simply as a "mean effective time", , for tracer particles to reach the distance L:

(B.3)

Our fitting procedure (Section 4.3) finds the optimal pair of values forandat a given distance L from the tracer source. Since we have that(Margolin and Berkowitz, 2000), then to make a prediction of the FPTD at a distance lL, we write

(B.4)

where the subscripts L anddenote the distance from the inlet at which the parameters are defined. Thus, given the parameter values (from fitting breakthrough measurements) forand at distance L, the predicted FPTD and CFPTD curves at a distance(as functions of time) are determined by evaluating the FPTD and CFPTD with the parameter valuesand.


Appendix C: Derivation of the CFPTD for 1 << 2

We derive here expressions for the CFPTD for 1 << 2. The procedure is similar to that shown in Appendix A for 0 << 1. The additional complications here are that

(C.1)

and that CFPTD= 1 only if we evaluate the integral of FPTD fromto t. The particles arriving during the time before the pulse are an artifact of the approximations; the asymptotic approximations made in the mathematical development (Margolin and Berkowitz, 2000) require the conditionat distance L (see below: CFPTD(t=0) must be negligibly small; this condition is required by equation 13). The reason is that for short times, the "tail" approximation breaks down, and the realbehavior must be considered. Due to normalization considerations,cannot scale asfor small t, because this will lead to a divergence. (Note that this divergence does not cause any difficulties for the case 0 << 1 because the rate of divergence is slower.)

We have then that

(C.2)

where c is a real, positive constant, and

.

(C.3)

With the change of variables (C.3) becomes

(C.4)



where c in the latter integrals represents a contour integral. Using the identities in (A.8) and (A.9), equation 12 follows immediately. This is an exact solution. Note also that equation 12 is equal to 1/Ò when h = 0 (i.e., t =; see Appendix D), where for Fickian transport (= 2), equation 12 equals 0.5. The approximation for large positive h is obtained by using the method of steepest descent. From (C.3), this yields equation 13. Conversely, for large negative h (i.e.,), from (C.3), and using the change of variablesand then we have

(C.5)

Using again the identity (A.8) leads to equation 14.


Appendix D: Derivation of axis translation parameters for 1 << 2

Two parameters must be defined in order to translate between dimensionless units in the FPTD and CFPTD solutions and the dimensional temporal units of actual measurements. As seen in equations 9-14, the basic independent variable is h. Because of the form of h, the "shift factor" is not a simple multiplicative factor as it is for the case 0 << 1. Rather, as seen from the last equality in D.1 below, two parameters control the shift (i.e., functional relationship between h and t is of the formfor parametersand).

With the definitions(whereis the average tracer velocity, which for all practical purposes is equal to the average fluid velocity; in contrast to the case 0 << 1,is almost independent of the length scale L),and again simplifying the analysis by choosing(see Appendix B), we have that

(D.1)

with the "mean time",for tracer particles to reach the distance L given by

.

(D.2)

Our fitting procedure (Section 4.4) finds the optimal set of values for,andat a given distance L from the tracer source. Then at a distance, we have thatwhere the subscripts L anddenote the distance from the inlet at which the parameters are defined. Moreover, because(Margolin and Berkowitz, 2000), i.e., it follows thatwhere in the last equality, we again choose.

Thus, the predicted FPTD and CFPTD curves at a distance(as functions of time) are given by solving the FPTD and CFPTD with the parameter values,and.