Published online 13 May 2005
Published in Vadose Zone J 4:360-379 (2005)
DOI: 10.2136/vzj2004.0125
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
REVIEWS AND ANALYSES
Simulating Solute Transport in Porous or Fractured Formations Using Random Walk Particle Tracking
A Review
Frédérick Delaya,*,
Philippe Ackererb and
Charles Danquignyb
a Univ. of Poitiers, Earth Sciences Building, 40 Av. du Recteur Pineau, 86022 Poitiers, Cedex, France
b Fluid and Solid Mechanics Institute of Strasbourg, 2 rue Boussingault, 67000 Strasbourg, France
* Corresponding author (fred.delay{at}hydrasa.univ-poitiers.fr)
Received 9 September 2004.
 |
ABSTRACT
|
|---|
Since the first attempts some 20 yr ago in the field of hydrology, random walk (RW) particle tracking as applied to solute transport has experienced profound changes. Concepts and mathematical techniques have improved to the point that numerically difficult problems (e.g., advection-dominated transport in highly heterogeneous media, or reactive transport) are now much easier to address. Random walk has never been widely used for multiphase flow, probably because numerical dispersion is not a major problem for modeling exercises at large scales. However, vadose zone hydrologic studies often point out very strong variations in fluid velocity over relatively short distances. Random walk methods may be well suited for such studies, a possibility which motivated us to write this review. We first give a comprehensive discussion of the theoretical context of the method. The Fokker–Planck–Kolmogorov equation (FPKE) is established for solute transport, as well as the ordinary Langevin equation and its simplifications for transport of small particles (e.g., colloids). Next, numerical methods are developed for the motion of particles in space. An important section is subsequently dedicated to recent RW concepts in the time domain, and to their application to anomalous (non-Fickian) transport and inverse problems. Adaptations of RW to transport with solute–solid reactions are also provided, as well as several numerical recipes for resolving a few computational difficulties with the RW method. We purposely did not include any comparisons with Eulerian and Lagrangian approaches. These approaches are discussed at length in several references cited in this review. We note, however, that today's computing capabilities provide new incentives to using RW methods for problems where Eulerian methods are potentially unstable or hampered by numerical diffusion.
Abbreviations: ADE, advection–dispersion equation • CTRW, continuous time random walk • ELLAM, Eulerian Lagrangian localized adjoint method • EPT, enhanced particle tracking • FPKE, Fokker–Planck–Kolmogorov equation • LEA, local equilibrium assumption • RW, random walk • TDRW, time domain random walk
|
INTRODUCTION
|
|---|
THERE IS AN INCREASING need for accurate numerical solutions to solute transport problems in porous and/or fractured formations. This need is predicated by increased concerns about soil and water quality, the longterm impact of industrial and agricultural contaminants on the environment, and the feasibility of underground repository sites for nuclear or other wastes. While considerable efforts during the past several decades have resulted in vastly improved numerical techniques for solving subsurface transport problems, accurate simulations of seemingly straightforward advection–dispersion problems often remain a challenge when transport is advection dominated. This problem becomes even more challenging when attempts are made to accurately account for the enormous heterogeneity of soil and groundwater systems (e.g., de Marsily et al., 1998, 2004).
Traditional Eulerian approaches to solute transport generally require very fine discretization of the transport domain to overcome recurrent problems of unstable numerical solutions and/or artificial diffusion. Solving such problems generally imposes a heavy computational burden. If the discretization is too coarse, oscillations and numerical diffusion of standard Eulerian methods may yield poor or even incorrect solutions, particularly when nonlinear processes such as multiphase flow or adsorption–desorption are considered in conjunction with standard advective–dispersive transport. This problem proved to be crucial in safety calculations of underground repository sites for nuclear wastes. To address these problems, several alternative Eulerian-Lagrangian schemes based on particle tracking have been developed, such as the method of characteristics (Konikow et al., 1996). Particle tracking has also been used for numerical integration in more complicated schemes, such as the Eulerian Lagrangian localized adjoint method (ELLAM) (Russell and Celia, 2002; Younès, 2004).
We provide a review of the principles of fully Lagrangian schemes based on RW particle tracking for solving advection–dispersion problems. Lagrangian approaches have long been used widely in physics for a variety of numerical applications. In the hydrologic sciences they were first applied in the late 1970s and early 1980s to transport in saturated alluvial formations (e.g., Ahlstrom et al., 1977; Prickett et al., 1981). Despite the encouraging results during the last two decades for simulating transport in saturated media, RW methods have been used only very rarely for unsaturated media. This is surprising since no inherent incompatibility exists between the RW method and unsaturated flow. On the contrary, soils are often characterized by strong contrasts in local flow velocities when a moisture front propagates during variably saturated flow. While Eulerian methods often experience numerical difficulties when simulating such velocity contrasts over short distances, Lagrangian models are much easier and more accurately to implement, as we shall demonstrate later.
In attempts to focus our review primarily on RW methods, we decided not to provide detailed comparisons between discrete Eulerian approaches (such as finite volume and finite element methods) and Lagrangian methods. Abundant literature already exists on such comparisons (we also provide several useful references). An extensive comparison between Eulerian and Lagrangian approaches would have required a much longer review. For these same reasons, we will not discuss all available RW methods, such as Boltzmann lattice gas approaches, the use of particles for numerical thermodynamics, stochastic properties of perfect gases, and molecular simulations in theoretical chemistry. Our study is aimed at solving macroscopic solute transport problems in porous or fractured media involving advection–dispersion as the principal mechanism. However, some discussions will be devoted to the ordinary Langevin equation which, as shown below, is well suited to mimic transport of solid particles subject to interaction forces in a moving fluid. This additional discussion seems justified here in that it may help the formulation of models for colloidal transport in porous media, and because the Langevin equation under some assumptions can be simplified into a slightly modified form of the advection–dispersion equation (ADE).
The main topics to be reviewed are as follows. The first section provides a comprehensive overview of the FPKE and the ordinary Langevin equations, and their relationship with the advection–dispersion equation for solute transport. The principles of stochastic physics on which these equations rely are simplified as best as possible to make this section readable without much background. A second section explains how to construct "explicit" solutions by randomly moving particles in space (classical RW particle tracking). Special emphasis is placed on algorithms that preserve solute mass balance at interfaces showing sharp contrasts in dispersion and on algorithms that accurately account for flow field heterogeneity. A third section presents recent concepts of random walk in the time domain. Two approaches are discussed: (i) a continuous time RW approach that calculates the probability density of particle residence times as a function of the fluid velocity distribution and (ii) a time domain RW method which moves particles explicitly between fixed sites of a network, but calculates random transition times between two successive locations. Applications of these methods to anomalous transport and inverse problems are suggested. A fourth section reviews adaptations of classical RW to transport with solute–solid interactions. Different methods, including several recent algorithms, are reviewed for simulating both instantaneous and first-order kinetic reactions. We conclude with a short section that provides a few numerical recipes and references to improve the RW efficiency and to implement boundary conditions.
|
FUNDAMENTALS
|
|---|
The theory and concepts in this section have been simplified as much as possible to obtain a treatise that is readable by most vadose zone hydrologists without being well-versed in theoretical physics. As a result, some statements may appear as assertions or conjectures compared with a more rigorous development. Several more extensive studies exist in current literature on stochastic physics; in particular, we refer more experienced readers to a recent review by Zaslavsky (2002).
The Fokker–Planck–Kolmogorov Equation
Let us assume that a particle moves at random and x(t) designates its location in Euclidean space 
at time t. Let e be the dimension of space. Furthermore, let 
(x, y, t) be the conditional probability density of a transition y 
x of duration t. ( ) is of dimension [L–e] and stationary in time, which means that the probability of a jump toward x, of duration t of a particle that was at y at time 
, is independent of . If the motion of the particle obeys a Markovian process, the particle "loses memories" between successive transitions. In other words, for successive times ti–1, ti, ti+1 and associated locations xi–1, xi, xi+1, (xi+1, xi, ti+1 – ti) does not depend on previous location xi–1 and lag-time ti – ti–1. Introduce now P(x, t), the probability density for a particle to be at location x at time t. P(x, t) is equivalent to (x, x0, t), that is, the probability density of a transition x0 x within duration t. If t is large enough for the transition x0 x to be performed in several steps x0 x1, ... , x
–1 x, x x because of the Markovian nature of this process, (x, x0, t) only depends on all possible locations x experienced before the last step x x. Thus, any reference to the initial location x0 may be dropped without loss of generality, and (x, x0, t) becomes P(x, t). A classical assumption of first-order derivatives allows us to write
 | [1] |
or with the equivalence between P( ) and ( ):
 | [2] |
An expression for (x, x0, t + 
t) for small values of t can be derived using the property of the Markovian process that
 | [3] |
which is the probability density to be at y at time t and to make a jump y x within duration t. This probability is summed over all possible locations y of the domain . Note that Eq. [3] is often referred to as the Chapmann–Kolmogorov equation for random particles. Introducing Eq. [3] in [2] yields
 | [4] |
When t tends to zero, no transition y x may occur except for x = y. Thus, we can write
 | [5] |
in which 
( ) is the Dirac-delta function [(0) = 1, (
0) = 0]. This property allows us to expand (x, y, t) as a series with respect to the derivatives of the ( ) function (Zaslavsky, 2002):
 | [6] |
where i( ) is the ith order derivative of () with respect to , and Mi(y, t) is the ith order moment of the transitions starting at y and being of duration t:
 | [7] |
Note that the exponent i in Eq. [7] refers to elevation to a power. To the first-order, M1(y, t) is the vector (e components) of the mean size of transitions starting at location y at any time and lasting t. To the second-order, M2(y, t) is the covariance matrix (e2 components) between the components along the e spatial directions of the above transitions. Replacing (x, y, t) in Eq. [4] by its expansion Eq. [6] up to the second order yields
 | [8] |
Knowing that the first and fourth terms of the right-hand side cancel, and using integration by parts for the second and third terms, leads to
 | [9] |
The so-called Kolmogorov assumptions state that the following limits exist:
 | [10] |
The third condition (i.e., moments of orders above 2 for transitions starting from x and of duration t cancel out) justifies the limited development of Eq. [9] up to the second order. The vector A(x) and tensor B(x) have dimensions of (L T–1) and (L2 T–1), respectively. Their common interpretation for the motion of particles by successive jumps is that A(x) corresponds to the mean of the jump velocity, while B(x) is the statistical dispersion of this velocity around its mean. Introducing Kolmogorov assumption (Eq. [10]) in Eq. [9] yields the so-called FPKE for transport:
 | [11] |
We show below that this equation has an explicit solution that for the simplest case consists of moving particles that undergo successive jumps of mean size A(x)t and covariance B(x)t, with t being the time step of the jumps.
The classical ADE for solute transport in a porous medium is given by
 | [12] |
in which C(x, t) (M L–3) is the volumetric concentration of the solute at location x and time t, u(x, t) (L T–1) is the mean fluid velocity vector, and D(x, t) (L2 T–1) is the dispersion tensor. The similarity of the ADE and FPKE should be evident by noting that
 | [13] |
Note in Eq. [13] that A(x) and B(x) are independent of time but compared with possible time-dependent values u(x, t) and D(x, t). This apparent discrepancy stems from the stationary assumption for ( ) and could mean that FPKE is not suited for transport under, for example, transient flow conditions. In fact, for practical and numerical applications, the stationary assumption can be relaxed provided transition times between successive particle locations are small compared with characteristic times of the flow field variations. Strict equivalence between the ADE and FPKE using Eq. [13] holds only if the dispersion tensor D(x, t) is constant in space; that is, 
(DC/x)/x = 2(DC)/x2 = D2C/x2. While this may the case for a homogeneous medium subject to flow with a uniform velocity u, such a condition is not frequently encountered in natural environments. However, the ADE given by Eq. [12] may be rewritten as
 | [14] |
In this form, and assuming the equivalence A(x) 
u(x, t) + D(x, t)/x instead of A(x) u(x, t), the ADE is similar to the FPKE. The importance of this modified velocity u(x, t) + D(x, t)/x will be discussed in more detail in a section below. The modified velocity is key in random walk applications to avoid mass-balance discrepancies at interfaces with sharp dispersion contrasts. Ackerer and Kinzelbach (1985) and Uffink (1985) were first in the hydrologic community to raise this problem. Ever since, much literature has been dedicated to algorithms able to solve the nonphysical excess of mass in low-dispersion areas when FPKE-based random particles are used without the modified velocity. We will revisit this topic below. While the form of A(x) u(x, t) + D(x, t)/x may appear surprising at first, we note that A(x) has the physical meaning of an estimate of the mean displacement per unit time, whereas u(x, t) is a unit flux of particles through a certain section (Uffink, 1990). This difference between A(x) and u(x, t) was also shown by Kitanidis (1988)(1994) using integration by parts of the ADE given by Eq. [12]. Note that when the porosity 
of the medium is variable in space and/or time, or when dealing with transport in unsaturated soils having a water content 
, the concentration C(x, t) in the ADE must be replaced by (x, t)C(x, t) or (x, t)C(x, t). The similarity with the FPKE is then ensured by using (x, t) C(x, t) P(x, t). The physical meaning is that CV is the solute mass in a bulk volume V of porous medium, while the mass conservation principle has a meaning equivalent to the property of probability density functions which states that their integral should be equal to one.
The Ordinary Langevin Equation
The ordinary Langevin equation of motion has been widely used to model the transport of immiscible gas bubbles or drops in fluids (Ramarao and Tien, 1992) and of atmospheric aerosols (Gupta and Peters, 1985). Its application to transport problems in natural porous media is less common. One reason is that its solution is cumbersome when dealing with complex velocity fields since the problem is addressed at the pore scale and requires, as a preliminary step, a solution of the Navier–Stokes equation. Explicit solutions for the flow field are still available for simplified systems such as the "sphere in cell" model (Happel, 1958) or periodic stacking of spherical beads. These models are useful, for instance, to simulate the interaction between suspended matter and the solid phase of the porous medium (Elimelech and O'Melia, 1990), and may find many applications involving colloid transport in soils. Nevertheless, the Langevin approach to transport requires, in theory, very small time steps for the best accuracy. While our aim here is not to discuss in detail potential applications for which the Langevin equation is relevant, it is included in our discussion for three reasons. First, the Langevin equation is not well known in current literature on groundwater. Second, its complete solution is based on a Lagrangian approach for moving particles in space and time. Third, the equation can be simplified with certain assumptions into the classical RW approach, which is the main topic of this review.
The ordinary Langevin equation is based on the conservation of momentum (mv) of a particle of mass m moving at velocity v in a fluid. In the following we assume that the Euclidean space is three-dimensional. Thus, except when specified, the variables involved are vector quantities depending on their space and time coordinates. The notation as a simple variable without reference to space coordinates is used here for the sake of simplicity. A classical form of the Langevin equation for transport of particles in porous media can be written as follows:
 | [15] |
where v(t) is the particle velocity (L T–1), u(t) is the fluid velocity (L T–1), while F(t) (M L T–2) combines interaction forces between both the particle and the solid phase and between the particles themselves. Additionally, A(t) (L T–2) is the acceleration due to Brownian motion stemming from collisions between the particle and fluid molecules, m is the mass of the particle, and ß = 6
µrp/m is the scalar Stokes coefficient of the drag force (T–1) in which µ is the dynamic viscosity of the fluid (M L–1 T–1), and rp the radius of the particle. Given the form of Eq. [15], the particle velocity is of the following general form:
 | [16] |
where K(t) is a vector of dimension (L T–1). Reintroducing Eq. [16] in [15] gives an equation for K(t):
 | [17] |
For a time step defined by the relative times 0 t, one obtains
 | [18] |
in which K0 = v0 (i.e., the value given by Eq. [16] at t = 0). Substituting Eq. [18] into [16] gives the following expression for the particle velocity:
 | [19] |
If the time step 0 t is small enough such that both the fluid velocities u() and the external forces F() can be assumed constant, Eq. [19] simplifies to yield a recursive algorithm for calculating the particle velocity at relative time t, provided its value is known at t = 0:
 | [20] |
The third term of the right-hand side of Eq. [20]; that is,
is random and expresses the Brownian component of the particle velocity. Its main characteristics will be discussed below.
If the particle velocity is available at successive time steps, then the particle trajectory can be calculated by integration of v over 0 t:
 | [21] |
where r(t) is the location of the particle at relative time t. Introducing Eq. [20] in [21] yields
 | [22] |
Integration by parts of the last term of the right-hand side gives
 | [23] |
Finally, the particle location at relative time t is
 | [24] |
The last term, that is,
is also a random vector corresponding to the Brownian component of the particle displacement during 0 t. Einstein's work between 1905 and 1908 on Brownian motion (e.g., as documented in a book edited by Furth, 1956) has shown that the acceleration A(t) is a vector of independent components A(t) ( = x, y, z directions) obeying Gaussian distributions and having the following properties:
 | [25] |
in which 
is the mathematical expectation, the Dirac-delta function as before, and q = kTß/m, where k is the Boltzmann constant (M L2 T–2 K–1) and T is absolute temperature (K). Using these properties of Brownian acceleration, Chandrasekhar (1943) showed that the vectors 
1(t) and 2(t) are correlated and obey multi-Gaussian distributions. It can be shown that more generally the vector quantities in three dimensions, that is,
have Gaussian probability density functions of the form:
 | [26] |
in which T indicates the transposition operator. Given the algebraic forms of 1(t) and 2(t) (see Eq. [20] and [24]) and Eq. [26], the first- and second-order moments of these distributions can be calculated as
 | [27] |
where i (i = 1, 2; = x, y, z) are components of the vector i(t) along direction .
Equations [20], [24], and [27] define a Lagrangian method for solving the ordinary Langevin equation. As stated above, the calculations are cumbersome since they require an accurate description of the flow field at the pore scale. Also, the assumption of having a constant fluid velocity and constant external forces during a time step forces the time steps to be small (e.g., on the order of 10–4 to 10–6 s for transport of stable clay suspensions in sand columns, as shown by Compere et al., 2001). In many cases, however, the above Lagrangian approach can be simplified to produce a slightly modified ADE. Assume that the time step 0 t is such that t >> 1/ß, and hence ßt >> 1 (common values of ß for clay particles of size 0.1 to 1 µm in water are in the range of 108 to 109 s–1). Equations [20] and [24] then simplify to
 | [28] |
Knowing that u0 + F0/ßm has the magnitude of v0, the term v0/ß is negligible compared with (u0 + F0/ßm)t. Moreover, 1(t) = 0 (see Eq. [27]), which allows Eq. [28] to be rewritten in the form
 | [29] |
The term 2(t) in Eq. [29] is known to be a vector with independent Gaussian components, while simplifications of Eq. [27] stemming for ßt >> 1 lead to 
= 2qt/ß2. With D = kT/6µrp (L2 T–1) designating the Stokes–Einstein diffusion coefficient, we have: 
2Dt and 1/ßm = D/kT. Equation [29] can now be rewritten as follows:
 | [30] |
where z is a vector having independent random components drawn from normal deviates of zero mean and unit variance. As further discussed below, Eq. [30] expresses a classical random walk procedure for advective–diffusive transport with a unit mass flux, J (M L–2 T–1): J = [u + (FD/kT)]C – D ·
C, which is the usual form when external forces have a nonnegligible influence on solute transport.
The same simplifications of the ordinary Langevin equation can be implemented in another manner, and probably with a more physical basis. Assume that the particles are very small and very light in weight, such that inertia in the Langevin equation can be neglected (the extreme would be particles with properties of fluid molecules, i.e., a "perfect" solute). Setting dv/dt in Eq. [15] to zero and assuming constant fluid velocities and external forces during a time step 0 t yields
 | [31] |
Since A(t) = 0 (see Eq. [25]), the mean velocity of the particles becomes v = u0 + F0/ßm. Calculation of the particle displacement from Eq. [31] is now straightforward to give
 | [32] |
Applying the general Eq. [26] to the term 
= 1/ß
t0 A
d results in
 | [33] |
which is a tri-Gaussian distribution of independent components (t) (L) with zero mean and a variance of 2qt/ß2. Using the earlier definition of the Stokes–Einstein diffusion coefficient D in the expressions for v and r(t) leads again to the classical random walk method given by Eq. [30].
|
NUMERICAL IMPLEMENTATION OF RANDOM WALK PARTICLE TRACKING IN SPACE
|
|---|
We have already shown that both the FPKE given by Eq. [11] and the simplified form of the Langevin equation given by Eq. [30] relate the transport mechanism to the motion of particles with well-established probability distributions for transitions between two successive locations. Basically, these equations are valid from the Darcy scale (say, a few ten thousand pores) up to large scales (from meters to kilometers). On the other hand, a full Langevin approach (Eq. [20], [24], and [27]) to transport in porous media is only applicable at the pore scale but proceeds in the same way, that is, moving particles with known transition probabilities. For the sake of simplicity we focus the following on mimicking the FKPE with particles. Such an approach assumes that the location of a particle evolves in time as follows:
 | [34] |
where x(t) is the vector of e coordinates defining the location of the particle at time t, e is the dimension of Euclidean space, and A(x) and B(x) are the vector and the tensor given by Eq. [10]. As discussed earlier, A(x)t is the mean displacement of the particle during time step t and B(x)t is its covariance, while Z is a vector of e independent random numbers drawn from a normal deviate (with zero mean and unit variance). As for the Langevin equation as discussed above, the multi-Gaussian distribution of stochastic jumps [B(x)t]0.5Z is inspired by the work of Einstein on Brownian motion and molecular diffusion. Since the vector Z is drawn independently for each particle at each time step, the particle moves randomly while following a general direction as dictated by vector A(x)t. This feature renders the name random walk to the method.
The first numerical implementations of RW in subsurface hydrology for solving the ADE were performed by Ahlstrom et al. (1977) and Prickett et al. (1981). They used a simplified RW based on the equivalence of FPKE and ADE as given by Eq. [13]. This yields the algorithm
 | [35] |
in which u[x(t)] and D[x(t)] are the mean pore-water velocity vector and the dispersion tensor, respectively, at location x(t). This standard equation is often referred to as the discrete form of the Itô (1951) stochastic differential equation. An alternative way of computing the dispersive term is given by Stratonovich (1966), where dispersion is taken at an intermediate location, that is:
 | [36] |
The Advection Step
Because of the random displacement of the particles, the velocity components must be known everywhere. The velocity field is usually obtained with a numerical flow model based on a finite volume or finite element method, and calculated on edges or faces of elements or at the element centroids. Various interpolation schemes to obtain the velocity at any location were compared by Goode (1990), Semra (1994), and LaBolle et al. (1996). Schafer-Perini and Wilson (1991) suggested that the interpolation scheme should preserve a null divergence of the velocity field. When a two-dimensional finite volume scheme with rectangular elements was used to compute the heads, they showed that fluid mass balance would be preserved only when the following interpolator within one element is used
 | [37] |
where x and y now are the two spatial coordinates, and a, b, c, and d are determined from the velocities at the edges of the element and the size of the element. For an element with its center located at the ith row and jth column of a regular grid, x0 = (i – 1/2, j) and y0 = (i, j – 1/2) are the coordinates of the lower left corner. The constants a, b, c, and d are given by
 | [38] |
The constants b and d correspond to the velocity gradient into the element along the x and y directions. Note that the method can be extended to three-dimensional flow and parallelepiped elements by simply adding the velocity in the z direction at z0 = (i, j, k – 1/2) and the velocity gradient along z. For a particle initially located at (x1, y1), integration of Eq. [37] gives its new location:
 | [39] |
after time step t (Pollock, 1988; Schafer-Perini and Wilson, 1991). Of course, the algorithm needs to be slightly modified when the time step t is fixed and the particle presumably leaves the element within the time step. The time step is then split such that when the particle reaches an element boundary, the exit point coordinates are stored, the constants for the velocity interpolation are updated with values calculated from the characteristics of the neighboring element, and the rest of the time step is allotted to motion of the particle within the neighboring element (e.g., see Pollock, 1988, for details of the algorithm). This interpolation is consistent with the scheme used in the mixed finite element method for rectangular elements (e.g., Hoteit et al., 2002), except that the scheme with the mixed finite element is given a priori. The method is therefore self-consistent, and no additional assumptions are required for interpolating the velocities, unlike with finite volume and finite difference methods.
The Dispersion Step
As shown by many authors (e.g., Ackerer and Kinzelbach, 1985; Uffink, 1987; Kinzelbach, 1988; Tompson and Gelhar, 1990; Kinzelbach and Uffink, 1991) and as discussed above when we established the equivalence between FPKE and ADE, the RW algorithm used to mimic the ADE requires a correction term that depends on the derivatives of the dispersion coefficients when the dispersion tensor varies significantly in space. Equivalence between the mean displacement of the FPKE and the modified velocity of the ADE (see Eq. [14]) allows algorithm [35] to be rewritten as
 | [40] |
This modification is straightforward if variations of the dispersion tensor in space are continuous enough to allow a first-order approximation. However, the algorithm needs to be modified when the derivatives of the dispersion coefficients are not defined, for example because of abrupt changes in the flow velocities. Abrupt changes may stem either from heterogeneities in the porous material, or from the discrete nature of the velocity field computed with a numerical model. For a nonuniform flow field, velocities may not always be unique at the interface of two adjacent elements. The discontinuity in dispersion due to a jump in the velocity at an element interface must be accounted for to avoid physically unrealistic results. This was clearly shown by Hoteit et al. (2002), who performed simulations with and without the correction term. The effect of discontinuous dispersion tensors has been widely studied for stratified media. The correction terms may be evaluated using two approaches: a reflection principle (Uffink, 1985; Cordes et al., 1991; Semra et al., 1993; LaBolle et al., 1998) and an interpolation technique (LaBolle et al., 1996; Bunzl, 2002).
The reflection principle is easy to implement in an efficient manner, even for three-dimensional problems. This may be done by generating an additional random number when a particle is expected to cross an interface with dispersion contrasts. If the random number is greater than the reflection coefficient, the particle crosses the interface. If not, the particle is reflected. The different reflection techniques have been compared by LaBolle et al. (1998) and discussed by Ackerer and Mosé (2000) and LaBolle and Fogg (2000). At the interface between two layers with different dispersion coefficients, a fraction R of the mass flux must be reflected to preserve mass balance in the dispersive flux on both sides of the interface. Following Uffink (1990), we consider here a one-dimensional diffusion process and an interface located at x = 0 between two layers denoted by and 
. Let the x coordinate be positive in layer and negative in layer , and assume that the particle is at x0 in layer . Using image theory and the superposition principle for linear processes, the probability density function of the particle displacement by diffusion in layer is given by
 | [41] |
where R is the reflected fraction of the mass flux. The first term of the right-hand side corresponds to diffusion up to a location x from a source at x0, and the second term to diffusion up to x from an image source that is symmetric of x0 with respect to the interface, and therefore located at –x0. The probability density function for layer is
 | [42] |
where the coefficient ß serves to modify the location of the source to account for the fact that diffusion is calculated for a homogeneous medium , whereas in reality diffusion occurs first in medium . Considerations on continuity at the interface (see Appendix) enable R and ß to be calculated as
 | [43] |
Semra et al. (1993) suggested the following numerical procedure to reflect particles. When a particle may cross the interface during a certain time step t, its displacement is split into two parts. During the first time step t1 the particle moves to the interface, and during the second time step t2 (=t – t1) the particle moves from the interface to either layer or layer . The displacement is then computed using the transport properties of the medium into which the particle entered. To preserve mass balance when the particle is at the interface, the probabilities that the particle enters layers and are P = 
/
and P = 1 – p = 
/, respectively. This algorithm must be used whether the particle originates from one side of the interface or the other one. It can be shown that this splitting is equivalent to the reflection principle suggested by Uffink (see Appendix).
An alternative to reflection methods is the interpolation technique. LaBolle et al. (1996) suggested to interpolate the velocities in the dispersion tensor (D = 
u, with the dispersivity) to smooth the dispersion tensor across the interface, and hence to dampen or even completely eliminate the discontinuities. However, as stated by the authors, an unbiased solution to transport requires not only very small time steps when the particle reaches the discontinuity, but also a fine spatial discretization for the interpolation scheme. Particle displacement during each time step must be significantly smaller than the width of the transition area over which the interpolation is performed. The size of this transition area is unknown a priori, but should be large enough to depict the change in the dispersion coefficient and small enough to avoid perturbations of the transport properties on either side of the discontinuity. The interpolation technique of LaBolle et al. (1996) hence leads to small time steps and high computational costs. A smoother interpolation was used by Bunzl (2002), who suggested the following equation:
 | [44] |
where D and D are the dispersion coefficients on both sides of the interface, xf is the location of the interface and w is the width of the transition zone.
We note here that LaBolle et al. (2000) suggested still another method for handling abrupt changes in dispersion. That method is based neither on reflection nor on strict interpolation. They changed the RW Eq. [35] for this purpose into
 | [45] |
where
 | [46] |
LaBolle et al. (2000) suggested to evaluate dispersion at location x + x, where the increment x is the dispersive step that would have been calculated in x(t) without correction. Their scheme was found to compare successfully with analytical solutions (LaBolle et al., 2000; Hassan and Mohamed, 2003).
An Alternative Using Cellular Automata
Delay et al. (1996) proposed a numerical alternative for solving RW in space with particles managed as variables over a regular grid. For simplicity we assume here that a two-dimensional domain is discretized into square cells of size x on a side. Our reasoning is the same for three-dimensional transport. Each cell contains a numerical variable representing the number of particles that can move by advection and dispersion during a certain time step t. If the RW algorithm with a uniform distribution of the dispersion jumps (Uffink, 1985) is applied to a particle located at (x, y) at time t, its position at t + t is
 | [47] |
where (x, y) represents a random displacement within a dispersion rectangle of size 2
, 2
centered at [x(t) + uxt,y(t) + uyt] and with its principal axis along the direction of the fluid velocity u. We used DL and DT for the longitudinal and transverse dispersion coefficients, respectively. Let us assume that all particles in cell j are assembled at its center (x, y) and represented by a single variable nj: their number in cell j. The jump t (t + t) will spread the particles uniformly over cells k that overlap the dispersion rectangle defined above. Each cell k receives a number of particles corresponding to the number nj of particles in cell j multiplied by the ratio of the area of cell k within the dispersion rectangle to the total area of the dispersion rectangle. With this method, all advection–dispersion jumps during time step t are performed sequentially, with j varying from 1 to the number of cells in the domain.
A classical scenario is as follows: the jump of particles from cell is performed before the jumps from cell while cell spreads its particles over . To avoid counting these particles in the jumps of cell , which would result in the anomaly of having two displacements during one time step, the number of dispersed particles that arrive in a cell are first stored and added to form an intermediate variable. When all jumps have been completed, the particles of the intermediate variable are assigned cell by cell to the variable representing the mean concentration that will move during the next time step. This enhanced particle tracking (EPT) scheme updates at each time step N variable (N being the number of cells) instead of moving particles individually and is therefore faster than classical RW methods. Moreover, as long as the velocity field remains constant in time, the jumps from each cell keep the same settings of the particle distribution. These settings in each cell (i.e., the neighboring cells under the dispersion rectangle and the fraction of mass assigned to those cells) can be stored in a table. The table is constructed only once at the beginning of the simulation; subsequent calculations at each time step then only need to follow the distribution rules of the table. A drawback of the method is that particles supposedly spread evenly over a cell are "reconcentrated" at the center of the cell before each time step. This may yield artificial dispersion. In most cases, it was observed that reconcentration at the cell center was able to balance eventual excess of spreading occurring when just a small portion of a cell is beneath the dispersion rectangle. In the end, artificial dispersion remains small.
Basically, the size of the dispersion rectangle is calculated with DL and DT values from the starting cell, which may become a problem when particles are supposed to cross boundaries with dispersion contrast. This was solved by improving the algorithm in using a constant displacement scheme. The concept was first proposed by Wen and Gomez-Hernandez (1996) for classical RW. These authors correctly reasoned that RW calculations should improve when keeping the size of the advection jumps constant instead of the time step. In classical calculations, the advection jump is of size ut. To correctly sample the velocity field, the time step t must be small enough so that, irrespective of the velocity u, the jump ut is smaller than the space step x of the grid over which the velocity field is depicted. This results in numerous jumps for low velocity areas, whereas only a few jumps are necessary for high velocity situations. The computation is improved with particles holding both their time and space position. For a jump at velocity ui, the local time step ti is adjusted so that the jump uiti is constant for any ui (e.g., uiti = x/2). After the jump (to which is also added a classical random motion due to dispersion), the location of the particles and their resting time are updated.
With the EPT method using dispersion rectangles, a constant displacement scheme may be generated as follows. An elementary time step t is defined such that for the maximum velocity of the flow field, umaxt is for instance equal to x/2. Suppose that a cell j with velocity uj requires nj elementary time steps for a jump ujnjt x/2. A series of nj intermediate variables is assigned to the cell and managed by means of a "first in–first out" sequence or queue. For each elementary time step t, the standard variable corresponding to the moving particles is spread over the other cells within a dispersion rectangle which size is given by Eq. [47] for a time step njt. The outlet of the queue is then flushed into the standard variable, while the ranking of the particles in the queue are decreased such that particles that arrive from other cells by advection–dispersion are stored at the inlet of the queue. Proceeding in this manner enables the particles to be stored during njt in the cell, and then to be moved by advection over an almost constant distance irrespective of the velocity in the cell. As compared with the algorithm of Wen and Gomez-Hernandez (1996), more jumps may be required since the EPT moves particles at each elementary time step t. However, the aforementioned principle of moving sets of particles with a table that stores the particle distribution remains unchanged and saves considerable computational time. The EPT method is rapid, even with a constant displacement scheme, while the use of constant jumps avoids the atavistic problems of mass conservation at the interfaces between layers with contrasting dispersion properties.
A similar approach for moving groups of particles was described by Vamos et al. (2003). The particles belonging to one cell are for this purpose first gathered at its center. They are subsequently moved to a neighboring cell following a Bernoulli distribution, depending on the advective fluxes between the starting cell and its neighbors. As stated by the authors, random fluctuations are dampened, but the scheme is equivalent to standard finite differences for rectangular cells, and therefore often generates artificial numerical diffusion.
|
APPROACHES TO RANDOM WALK IN THE TIME DOMAIN
|
|---|
In the preceding section we showed that classical RW methods are based on the motion of particles in space. They are displaced during each time step over the domain according to local properties of the flow field (see e.g., Eq. [40]). In very heterogeneous media such as fractured rocks and macroporous soils where flow occurs in both the matrix and the fractures, RW calculations may become very time-consuming. Several theoretical studies have shown the absence of a homogenization scale for the connectivity and the flow properties of such media (e.g., de Dreuzy et al., 2001a, 2001b, for fractured networks). Very heterogeneous media are often modeled as discrete networks involving tens of thousands to even a few millions bonds (Bour and Davy, 1997; Adler and Thovert, 1999; Rivard and Delay, 2004). In addition, flow velocities may span several orders of magnitude, which imposes very small time steps for accurate sampling of the flow field using a classical RW. A large number of jumps are required in low velocity areas to move the particles significantly, which may cause the calculations to become inefficient in terms of computational costs. These constraints have motivated the development of new approaches based on random particles managed in the time domain. There are two reasons for this. One is to provide a framework for describing nonGaussian transport often observed in field and laboratory studies (Scher et al., 2002; Bromly and Hinz, 2004; Cortis and Berkowitz, 2004). A second reason is to produce efficient and rapid methods for calculating transport problems, for example for wide bond networks, while accounting for processes other than pure advection.
One assumption of the use of time domain approaches applied to subsurface hydrology is that the flow field can be represented by a set of spatially distributed sites connected by bonds that ensure the transition of particles between sites. This is an efficient abstraction for fracture networks (e.g., Dershowitz and Fidelibus, 1999) and for heterogeneous porous media, which may be viewed as regular percolation networks with varying hydraulic conductivities in the bonds (e.g., Sahimi and Mukhopadhyay, 1996). The particles are also assumed to follow a Markovian process in that there is no correlation between successive transitions from one site to the other.
The Continuous Time Random Walk
The continuous time random walk (CTRW) theory for lattices was first developed by Montroll and Weiss (1965), Montroll and Scher (1973), and Scher and Lax (1973) for problems involving electron hop conductivities in semiconductors. The theory was recently enhanced and extended to anomalous transport in fracture networks (Berkowitz and Scher, 1998) and more generally to groundwater hydrology in heterogeneous media (Berkowitz et al., 2001; Berkowitz and Scher, 2001). Contrary to classical RW moving particles through a given domain, CTRW is an up-scaled approach that does not use particles explicitly. The theory is based on calculations of P(s, t), the probability density for a particle to be at site s at time t. The aim is to depict the macroscopic behavior of the system, with the first abstraction considering that the real domain can be averaged by sites located on a regular lattice with N sites on a side and involving periodic boundary conditions. To avoid any loss of generality, the jumps are not limited to adjacent sites, which means that particles can move from site s' to site s with a distance s – s' spanning all sizes available in the lattice. This is an important feature when CTRW is applied to random fracture networks where the distance between adjacent sites at fracture intersections may vary from almost zero to the size of the network. If R(s, t) designates the probability density for a particle to arrive at site s and time t, and 
(d, t) is the probability density of making a jump of size d of duration t, the Markovian process followed by the particles allows us to write
 | [48] |
which is the probability of just arriving at s' at time multiplied by the probability of moving from s' to s within time t – . This local probability is extended to each possible site s' and each time between 0 t. The (s – 0) and (t – 0) Dirac-delta functions represent the initial conditions to ensure that the probability for a particle to be in 0 at time 0 is one. We emphasize that the form of Eq. [48] does not presume any particular transport process in the bonds. Potentially several processes can occur simultaneously, but they need to be merged in the probability (d, t). In other words, small-scale transport details are supposed to be encapsulated in ( ) before averaging them when calculating P( ). The probability P(s, t) can be written as:
 | [49] |
which is the probability of arriving at site s at time between 0 t, multiplied by the probability (t – ) of staying in place during t – . The variable (t) in Eq. [49] may be expressed as the probability that no jump occurs during t:
 | [50] |
Note that () is the probability of doing a jump of any size and of duration . Equations [48] and [49] are convolution products in time and space. Their calculation can be simplified if Laplace–Fourier transforms are used. For two functions f(t) and g(s), the Laplace and discrete Fourier transforms are defined as
 | [51] |
The Laplace and Fourier transforms of convolution products are the simple algebraic products of the Laplace and Fourier transforms. Therefore, in using the properties L[(t)] = 1 and F[(x)] = 1, Eq. [48] yields
 | [52] |
Integration Eq. [50] by parts and taking the Laplace transform of (t) leads to
 | [53] |
where L(p) is the Laplace transform of () = 
d(d, ). The probability density P(s, t) given by Eq. [49] can be calculated as the inverse Fourier transform of PF(w, t) as follows:
 | [54] |
in which L–1 is the inverse Laplace transform, N the number of sites, and e the dimension of Euclidean space. P(s, t) expresses the evolution in time of the concentration plume at location s. Another probability density relevant to solute transport is that of the first arrival time at a location s. This probability is equivalent to a breakthrough curve that would be sampled at the outlet s of a laboratory column in response to a Dirac injection at the inlet. Under these conditions and for transitions 0 s, we have
 | [55] |
where F(s, t) is the probability of arriving for the first time at s at time t. Note here that Dirac-delta function ( ) are not required for initial conditions since they are already in R(0, t) (see Eq. [48]). Equation [55] expresses the fact that a particle just arriving at s at t may have visited the site for a first time at , and then may have moved back and forth during t – . In the Laplace domain, one obtains
 | [56] |
Assuming that the Laplace and Fourier inversions can be calculated either numerically or analytically, the central notion of the CTRW theory is the probability density (d, t) of doing a jump of size d during time t. As stated above, ( ) is assumed to enclose small-scale details of transport mechanisms in the real domain. For transport by advection in two-dimensional fracture networks, Berkowitz and Scher (1998) decomposed this probability as follows:
 | [57] |
in which K is a normalization constant, 
(u) is the probability density of the fluid velocity u in the bonds, and p(d|u) is the conditional probability of a transition of size d knowing that it occurs at velocity u. In other words, since t = d/u, p is the probability of a transition of size d knowing that it occurs within time t. Also, f(u) in Eq. [57] is the probability for a particle to experience the velocity u in the network. Assuming perfect mixing of mass fluxes at bond intersections, this probability at the scale of the network is the ratio of the total flow rate through bonds with velocity u to the total flow rate in the network.
Berkowitz and Scher (1998) performed calculations of flow through several realizations of two-dimensional random fracture networks having exponential law length distributions [i.e., the number of fractures of length l are given by n(l) 
exp(–l/l0)]. They found the best fits for p(d|u) and (u) to be
 | [58] |
where d0, u0, and ß are fitting parameters. Whereas d0 and u0 are the characteristic distance and velocity that normalize the probability densities, ß is a key parameter that expresses the typical feature of fracture networks to exhibit strong asymmetry in the velocity distribution. With u = d/t, Eq. [57] becomes
 | [59] |
Introducing Eq. [59] into Eq. [54] and [56] with a value of ß in the range 0.5 to 0.9, Berkowitz and Scher (1998) demonstrated the persistent nature of anomalous (non-Gaussian) behavior of solute plumes and breakthrough curves, even for long residence times. For ß values of approximately 0.5, the plumes were highly asymmetric with a peak remaining close to the injection point and a flat forward front due to only a few particles experiencing very rapid transitions. On the other hand, the breakthrough curves showed flat tails for large residence times. For ß values of about 0.8 to 0.9, the peak moved away from the injection point, with a forward tail in space becoming thinner and shorter than the backward tail. The centers of mass of the plumes along the main flow 
TOP
ABSTRACT
INTRODUCTION
