Abstract
A generalized form of the random walk algorithm to simulate diffusion processes is introduced.
Unlike the usual approach, at a given time all the particles from a grid node are simultaneously scattered using the Bernoulli repartition. This procedure saves memory and computing time and no restrictions are imposed for the maximum number of particles to be used in simulations.
We prove that for simple diffusion the method generalizes the finite difference scheme and gives the same precision for large enough number of particles.
As an example, simulations of diffusion in random velocity field are performed and the main features of the stochastic mathematical model are numerically tested.
Authors
Călin Vamoş
Tiberiu Popoviciu Institute of Numerical Analysis (Romanian Academy)
Nicolae Suciu
Tiberiu Popoviciu Institute of Numerical Analysis (Romanian Academy)
Harry Vereecken
Keywords
Diffusion; random walk; groundwater; contaminant transport
Cite this paper as
C. Vamoş, N. Suciu, H. Vereecken (2003), Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186(2), 527-544, doi: 10.1016/S0021-9991(03)00073-1
About this paper
Print ISSN
Not available yet.
Online ISSN
Not available yet.
Google Scholar Profile
[1] W.F. Ames, Numerical Methods for Partial Differential Equations, second ed., Academic Press, New York, 1977.
[2] S. Attinger, M. Dentz, H. Kinzelbach, W. Kinzelbach, Temporal behavior of a solute cloud in a chemical heterogeneous porous medium, J. Fluid Mech. 386 (1999) 77.
CrossRef (DOI)
[3] M. Avellaneda, A.J. Majda, Superdiffusion in nearly stratified flows, J. Stat. Phys. 69 (3/4) (1992) 689.
CrossRef (DOI)
[4] A. Bellin, P. Salandin, A. Rinaldo, Simulation of dispersion in heterogeneous porous formations: statistics, first-order theories, convergence of computations, Water Resour. Res. 28 (9) (1992) 2211.
CrossRef (DOI)
[5] C.W. Gardiner, Handbook of Stochastic Methods, Springer, New York, 1983. C. Vamos et al. / Journal of Computational Physics 186 (2003) 527–544 543
[6] S.K. Godunov, V.S. Ryabenkii, Difference Schemes: An Introduction to the Underlying Theory, North-Holland, Amsterdam, 1987.
[7] A.J. Chorin, Vortex sheet approximation of boundary layers, J. Comput. Phys. 27 (1978) 428.
CrossRef (DOI)
[8] A.F. Ghoniem, F.S. Sherman, Grid-free simulation of diffusion using random walk methods, J. Comput. Phys. 61 (1985) 1.
CrossRef (DOI)
[9] W. Horsthemke, R. Lefever, Noise-induced Transitions. Theory and Applications in Physics, Chemistry and Biology, Springer, Berlin, 1984.
[10] C. Kapoor, L.W. Gelhar, Transport in three-dimensionallity heterogeneous aquifers 1. Dynamics of concentration fluctuations, Water Resour. Res. 30 (6) (1994) 1775.
CrossRef (DOI)
[11] C. Kapoor, L.W. Gelhar, Transport in three-dimensionallity heterogeneous aquifers 2. Prediction and observations of concentration fluctuations, Water Resour. Res. 30 (6) (1994) 1789.
CrossRef (DOI)
[12] W. Kinzelbach, Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser, Oldenbourg Verlag, Munchen, 1992. €
[13] G.L. Moltyaner, M.H. Klukas, C.A. Willis, R.W.D. Killey, Numerical simulations of twin lake natural-gradient tracer tests: A comparison of methods, Water Resour. Res. 29 (10) (1992) 3443.
CrossRef (DOI)
[14] O. Neuendorf, Numerische 3D-Simulation des Stofftransport in einem heterogenen Aquifer. Ph.D. Thesis, Jul-3421 Forschungs- € zentrum Julich, 1997. €
[15] A. Papoulis, Probability, Random Variables and Stochastic Processes, McGraw-Hill, New York, 1991.
[16] K. Roth, K. Hammel, Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with steady state flow, Water Resour. Res. 32 (6) (1996) 1653.
CrossRef (DOI)
[17] P. Salandin, V. Fiorotto, Solute transport in highly heterogeneous aquifers, Water Resour. Res. 34 (5) (1998) 949.
CrossRef (DOI)
[18] H. Schwarze, U. Jaekel, H. Vereecken, Estimation of macrodispersivity by different approximation methods for flow and transport in randomly heterogeneous media, Transp. Porous Media 43 (2001) 265.
[19] Ne.-Z. Sun, Mathematical Modeling in Groundwater Pollution, Springer, New York, 1996.
CrossRef (DOI)
[20] A.F.B. Tompson, L.W. Gelhar, Numerical simulation of solute transport in three-dimensional, randomly heterogeneous porous media, Water Resour. Res. 26 (10) (1990) 2541.
CrossRef (DOI)
[21] A.F.B. Tompson, R.B. Knapp, Reactive Geochemical Transport Problems in Nuclear Waste Analyses, Preprint UCRL-00552, Lawrence Livermore National Laboratory, Livermore, 1989.
[22] A.F.B. Tompson, R.D. Falgout, S.G. Smith, W.J. Bosl, S.F. Asby, Analysis of subsurface contaminant migration and remediation using high performance computing, Adv. Water Resour. 22 (3) (1998) 203.
CrossRef (DOI)
[23] C. Vamos, N. Suciu, H. Vereecken, O. Nitzsche, H. Hardelauf, Global random walk simulations of diffusion, in: W. Kramer, J.W.V. Gudenberg (Eds.), Scientific Computing, Validated Numerics, Interval Methods, 343, Kluwer Academic Publishers, Dordrecht, 2001.
CrossRef (DOI)
[24] R. Zhang, K. Huang, M.T. van Genuchten, An efficient Eulerian–Lagrangian method for solving solute transport problems in steady and transient flow fields, Water Resour. Res. 29 (12) (1993) 4131.
CrossREf (DOI)
' d='M65 980.1h54.8' class='g2'/%3E%0A%3C/svg%3E)
Generalized random walk algorithm for the
numerical modeling of complex diffusion processes
Calin Vamos ß
a
, Nicolae Suciu
a,
*
, Harry Vereecken
b
a
Romanian Academy, ‘‘T. Popoviciu’’ Institute of Numerical Analysis, Republicii Str. 37, P.O. Box 68, 3400 Cluj-Napoca 1, Romania
b
Forschungszentrum J€ulich GmbH, Institut f€ur Agrosph€are (ICG-IV), D-52425, J€ulich, Germany
Received 11 December 2001; received in revised form 10 January 2003; accepted 25 January 2003
Abstract
A generalized form of the random walk algorithm to simulate diffusion processes is introduced. Unlike the usual
approach, at a given time all the particles from a grid node are simultaneously scattered using the Bernoulli repartition.
This procedure saves memory and computing time and no restrictions are imposed for the maximum number of
particles to be used in simulations. We prove that for simple diffusion the method generalizes the finite difference scheme
and gives the same precision for large enough number of particles. As an example, simulations of diffusion in random
velocity field are performed and the main features of the stochastic mathematical model are numerically tested.
Ó 2003 Elsevier Science B.V. All rights reserved.
Keywords: Diffusion; Random walk; Groundwater; Contaminant transport
1. Introduction
It is well known that diffusion processes can be numerically simulated with random walk (RW) algo-
rithm. For simple diffusion processes the RW algorithm is equivalent with the finite difference (FD) scheme
[1] but, as we shall discuss in the following, this equivalence is not valid for more complex diffusion pro-
cesses. The RW algorithm can be used to model the transport of arbitrary physical quantities if parts of the
transported quantity are associated with fictitious particles obeying the RW law. To reduce the compu-
tational effort and to improve the smoothness of the numerical solution, the gradient of the transported
quantity can be associated with the particles. Integrating the gradient transported by each particle over the
computation domain the simulated field is obtained with a higher accuracy [8]. The ‘‘gradient random
walk’’ algorithm was first developed by Chorin [7] for the simulation of turbulence, the transported
quantity being the vorticity. In other applications, mainly for transport in porous media, to reduce memory
and to avoid boundary effects, a grid free algorithm called ‘‘particle tracking’’ (PT) is used [22]. The
Journal of Computational Physics 186 (2003) 527–544
www.elsevier.com/locate/jcp
*
Corresponding author.
E-mail addresses: cvamos@ictp.acad.ro (C. Vamos ß), n.suciu@fz-juelich.de (N. Suciu), h.vereecken@fz-juelich.de (H. Vereecken).
0021-9991/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved.
doi:10.1016/S0021-9991(03)00073-1
algorithm consists in building trajectories in continuous space for each particle, by performing at discrete
time steps an advection displacement and a random Gaussian one.
The application of RW algorithm to solve some practical transport problems is relatively limited. A
good estimation of the concentration field requires a large number of particles at each grid point. For every
jump of each particle a random number should be generated implying a certain number of numerical
operations. Therefore, ‘‘if we require that the numerical solution should be identical to the analytic solution
to three significant figures, the number of tracer particles simulated must be enormous’’ [19, p. 95]. In this
paper we propose an improvement of the RW algorithm aiming to eliminate this limitation: all the particles
from a given grid node are moved simultaneously, not individually. This procedure is possible since the
number of the particles jumping from a given node to a single neighboring node obeys a Bernoulli dis-
tribution. In this way a great number of particles can be distributed generating only a single random
number and the necessary random number generations are significantly reduced. We call this algorithm
‘‘Global Random Walk’’ (GRW). A more general form of the GRW algorithm is obtained when a part of
the particles remain at the initial node and only the rest of them are scattered to the neighboring nodes
according to the Bernoulli distribution.
The number of particles jumping from a given node to a single neighboring node according to GRW
fluctuates about the mean value. These fluctuations can be eliminated if we allow the particles to be
divided into parts and exactly half of them to jump to the left and the other half to the right. We show
that for simple diffusion processes this ‘‘determinist’’ GRW algorithm (GRWD), without fluctuations, is
identical with FD algorithm. Thus, GRW can be thought as a generalization of FD algorithm. If we do
not intend to give up the particle indivisibility, then the fluctuations cannot be completely canceled. The
minimum magnitude of the particles number fluctuations is equal to a single particle. In this case a form
of GRW with reduced fluctuations, denoted GRWR, is obtained: if at a node there are an odd number of
particles, then one of the particles is randomly distributed to the left or right and the rest of them is
divided to half.
To show that the RW algorithm gives a numerical solution of the diffusion equation, let us consider an
infinite lattice given by a grid with nodes at x
i
¼ idx, i 2 Z, where dx > 0 is the lattice constant or the space
step. At time t
k
¼ kdt, k 2 N, where dt > 0 is the time step, the particle jumps to the left or right neighboring
node with equal probabilities. If we denote by P ðx
i
; t
k
Þ the probability distribution at t
k
and by P ðx
i
; t
k
jx
j
; t
l
Þ
the transition probability, then the consistency property of finite dimensional probabilities for discrete
Markov processes gives
P ðx
i
; t
k
Þ¼
X
j
P ðx
i
; t
k
jx
j
; t
l
ÞP ðx
j
; t
l
Þ; ð1:1Þ
where t
k
> t
l
. According to the RW law, the transition probabilities for successive time steps are given by
P ðx
i
; t
lþ1
jx
j
; t
l
Þ¼
1
2
for i ¼ j 1;
0 for i 6¼ j 1:
ð1:2Þ
The Chapmann–Kolmogorov equation and (1.2) give the evolution equation for the transition probability
P ðx
i
; t
kþ1
j0; 0Þ¼
1
2
P ðx
i1
; t
k
j0; 0Þ ½ þ P ðx
iþ1
; t
k
j0; 0Þ: ð1:3Þ
One can prove (see for instance [5]) that if dx ! 0 and dt ! 0, and if the limit
lim
dx;dt!0
dx
2
2dt
¼ D ð1:4Þ
is finite, then the limit of the solution of (1.3) is the Gaussian transition probability with the density
528 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
pðx; tj0; 0Þ¼ð4pDtÞ
1=2
exp
x
2
4Dt
: ð1:5Þ
The Gaussian process is a continuous Markov process and the probability density that a particle has the
position x at time t,
pðx; tÞ¼
Z
R
pðx; tjx
0
; 0Þpðx
0
; 0Þ dx
0
; ð1:6Þ
is the solution of the diffusion equation
o
t
p ¼ Do
2
x
p; ð1:7Þ
with the initial condition p
0
ðxÞ¼ pðx; 0Þ [5].
The connection between the RW process, described by (1.1) and (1.2), and the solution (1.6) of the
diffusion equation is used to build up the RW algorithm to solve the diffusion equation (1.7). One considers
N fictitious particles moving on the previously described infinite grid, their movement being governed by
the RW law. At the initial time the distribution of the N particles approximates the values of the initial
probability density p
0
ðxÞ. For each particle and each time step a random number taking the values )1 and
+1 with a probability equal to 1/2 is generated. If the random number is )1, then the particle moves to the
left and if it is +1, the particle moves to the right. In this way the distribution of the N particles at the time t
k
approximately describes the solution pðx; t
k
Þ of the diffusion equation (1.7). This approximation can be
improved by increasing N and reducing dx and dt.
In Section 2 we present the GRW algorithm and discuss the connections with FD algorithm. In Sections 3
and 4 we analyze GRW simulating the Gaussian solution of the one-dimensional diffusion equation far from
the boundaries and with several boundary conditions. We show that for a large enough number of particles
N , the numerical solution obtained with GRW coincides with the FD solution, for all spatial resolutions.
In Section 5, we simulate the diffusion in random velocity fields. In this case GRW is no more equivalent
with FD algorithm. Instead of (1.7), we consider advection–diffusion equations of the form
o
t
c þ V ðxÞo
x
c ¼ Do
2
x
c; ð1:8Þ
where the advection velocity V ðxÞ is a space function and the diffusion coefficient D is a constant. The
relation between the concentration field c in (1.8) and the probability distribution p is given by
cðx; tÞ¼ N
0
pðx; tÞ; ð1:9Þ
where N
0
is the total number of diffusing particles,
N
0
¼
Z
1
1
cðx; tÞ dx: ð1:10Þ
If the advection velocity has large spatial variations then V ðxÞ is modeled as a random velocity field. Such
problems occur in studies of transport processes in heterogeneous porous materials. Many authors agree
that numerical simulations of diffusion in random fields are better achieved by RW algorithms than by
usual FD or finite element algorithms [13,16,20–22]. We shall show that GRW allows the simulation of the
two-dimensional diffusion in random velocity fields using moderate computing resources.
2. The GRW numerical algorithm
The GRW algorithm is identical with the RW algorithm in its mathematical principles, only the nu-
merical implementation of the particles displacement is different. Hence, first we describe the RW algorithm
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 529
for Eq. (1.8) using the same lattice as for (1.1)–(1.3). According to (1.4), we relate the lattice constant dx to
the time step dt, for a given constant diffusion coefficient D, by
D ¼
dx
2
2dt
: ð2:1Þ
The particles in x
i
jump either in ðx
i
dxÞ or in ðx
i
þ dxÞ. Similarly, the realization of the random velocity
field is described on the grid by a set of integers v
i
, defined by
V ðx
i
Þ¼
dx
dt
v
i
; ð2:2Þ
so that the movement of particles in x
i
due to advection field is given by ðx
i
þ v
i
dxÞ. The total displacement
of particles is obtained, similarly to PT algorithm [20], as the sum of advective and diffusive displacements.
After a time step dt, the particles starting from the node i reach either the node ði þ v
i
1Þ or ði þ v
i
þ 1Þ.
The shortcoming of this approach is that V ðx
i
Þ has only discrete values and if V ðxÞ has large variations, then
dx=dt must be small imposing a small space step. Therefore, in this case very large grids are needed. Also,
the velocity V ðx
i
Þ must be replaced by its average over a time step. For smooth variations, it was shown that
ðV ðx
i
Þþ V ðx
i
þ dtV ðx
i
ÞÞÞ=2 approximates the average up to second order in dt [16].
At a given time t
k
¼ kdt, the N particles are distributed on the grid so that nði; kÞ is the number of
particles at the node x
i
¼ idx. Moving each particle according to the RW law, the distribution of the
particles at the next time step nði; k þ 1Þ is obtained. Repeating S times the simulation under the same initial
condition, we obtain for each node i and time k a sequence of values n
s
ði; kÞ,1 6 s 6 S . When S is large
enough, we can approximate the solution cðx; tÞ of the advection–diffusion equation (1.8) using the average
over the S realizations of n
s
ði; kÞ, denoted nði; kÞ. Assuming that the particles in the node x
i
are assigned to
an interval consisting of l space steps, the numerical solution of (1.8) is given by
cðx
i
; t
k
Þ¼
1
ldx
nði; kÞ: ð2:3Þ
Usually l ¼ 1, but it is possible that only half of the grid nodes can contain particles and then l ¼ 2 (see
Section 3).
The GRW algorithm moves the particles in large groups, not individually. We denote by dnði; j; kÞ the
number of particles which at time t
k
jump from x
j
to x
i
. Then the distribution of the particles at the next
time is obtained from the relation
nði; k þ 1Þ¼ dnði; i; kÞþ
X
j6¼i
dnði; j; kÞ; ð2:4Þ
where dnði; i; kÞ is the number of particles remaining at t
kþ1
at the same node x
i
. For every j and k, the
scattering of the particles is given by the relation
nðj; kÞ¼ dnðj; j; kÞþ dnðj þ v
j
d ; j; kÞþ dnðj þ v
j
þ d ; j; kÞ; ð2:5Þ
where the positive integer d describes the jumps of particles due to diffusion. From (2.5) it follows that the
contributions dnði; j; kÞ, j 6¼ i, in (2.4) come from all points j satisfying j þ v
j
d ¼ i. The GRW algorithm
is defined if a procedure to calculate the values of the quantities in the right hand part of (2.5) is given.
We want that at a given time step only a fraction r of the number of particles jump in the neighboring
nodes, the rest of them remaining at the same node. To avoid the division of particles, r must be a positive
rational number r 6 1, such that ð1 rÞN be an integer and equal to the total number of particles remaining
at the same node at a time step. For increasing index j, we determine the number of particles remaining at a
node x
j
by means of the formula
530 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
dnðj; j; kÞ¼ ð1
2
4
rÞ
X
j
0
6 j
nðj
0
; kÞ
3
5
ð1
2
4
rÞ
X
j
0
<j
nðj
0
; kÞ
3
5
; ð2:6Þ
where ½ is the integer part of the expression in the brackets. Taking the average over a great number of
Monte Carlo realizations we obtain
dnðj; j; kÞ¼ð1 rÞ nðj; kÞ: ð2:7Þ
Since dnðj; j; kÞ is known from (2.6), (2.5) relates the random variables dnðj þ v
j
d ; j; kÞ and
dnðj þ v
j
þ d ; j; kÞ and only one of them has independent values.
As a consistency requirement, for a given diffusion process, the GRW algorithm must give the same
mean square displacement as the RW algorithm. If dx
RW
is the space step for the RW algorithm, then for a
time step dt the mean square displacement of the particles in the node j is nðj; kÞdx
2
RW
, because all the
particles jump at the first neighbors. For an equal time step, in GRW algorithm only the fraction r of the
particles in the node j jump at the nodes j d and, from (2.5) and (2.7), the mean square displacement is r
nðj; kÞðd dxÞ
2
. Imposing the condition nðj; kÞdx
2
RW
¼ r nðj; kÞðd dxÞ
2
and using (2.1) for dx
RW
, the parameter r
is given by
r ¼
2Ddt
ðd dxÞ
2
: ð2:8Þ
Usually the values of dx and dt are mainly determined by the necessity to describe enough accurately the
spatial variations of the velocity field V ðxÞ and by the limitations of the available memory and computing
speed. Then for any value of the diffusion coefficient D, even if it has spatial variations, the relation (2.8)
allows us to choose the values of d and r. This is not possible using only the parameter d , because it takes
discrete values. On the other hand r is a continuous parameter, but its value is restricted by the condition
r 6 1. Combining the unboundedness of d and the continuity of r we can obtain any value of D.
The GRW algorithm performs the evaluation of the random variables dnðj þ v
j
d ; j; kÞ directly, not as
a sum of the individual jumps of the n nðj; kÞ dnðj; j; kÞ particles. Since each of the n particles can reach
the node j þ v
j
d with a probability equal to 1=2, it follows that the probability for dnðj þ v
j
d ; j; kÞ to
take the value m,0 6 m 6 n, is given by the Bernoulli distribution b
n
ðmÞ¼ 2
n
C
m
n
. To assign to dnðj þ
v
j
d ; j; kÞ a random value satisfying the Bernoulli distribution, at each time step, a random number g with
a uniform distribution in the interval ½0; 1 is generated. If we denote by F
n
ðmÞ¼
P
m
l¼0
b
n
ðlÞ,0 < F
n
ðmÞ 6 1,
the Bernoulli repartition, then dnðj þ v
j
d ; j; kÞ takes the value m satisfying the condition F
n
ðm 1Þ 6
g < F
n
ðmÞ, where we use the convention F
n
ð1Þ¼ 0.
To analyze the relation between the GRW and the FD algorithms, we consider the centered differences
and time explicit scheme for diffusion equation obtained from (1.8) when V ðxÞ 0. Considering the ap-
proximation in the order Oðdx
2
Þ of o
2
x
c in (1.8), using finite difference between d dx space steps and the
parameter r defined by (2.8), the explicit FD scheme can be written as
cði; k þ 1Þ¼
r
2
cði þ d ; kÞþð1 rÞcði; kÞþ
r
2
cði d ; kÞ: ð2:9Þ
The solution of (2.9) is stable if the von Neumann stability criterion, r 6 1, is fulfilled. Since from (2.8) we also
have dt ¼ Oðdx
2
Þ, the FD scheme (2.9) is a consistent approximation of the exact partial differential equation
within the approximation order Oðdx
2
Þ. The stability and consistency imply the convergence of the order
Oðdx
2
Þ for the initial value problem attached to (1.8) with V ðxÞ 0 [6]. For v
i
0, the relation (2.4) becomes
nði; k þ 1Þ¼ dnði; i; kÞþ dnði; i þ d ; kÞþ dnði; i d ; kÞ:
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 531
Taking into account (2.3) and (2.7), the average of this relation is identical with (2.9) if
dnðj d ; j; kÞ¼
1
2
r nðj; kÞ: ð2:10Þ
But this is the RW law stating that the average number of particles jumping in a direction is equal to half
the total number of particles. This proves that the FD solution is identical with the ensemble average of the
GRW solutions. The parameter (2.8) defining the particles fraction jumping to the neighboring nodes in
(2.6) is the same as the stability parameter of the FD scheme (2.9).
One can define a modified GRW algorithm which is identical with the FD algorithm for V ðxÞ 0, if the
particles can be divided and nðj; kÞ is a real number, not an integer. Instead of (2.6) we introduce
dnðj; j; kÞ¼ð1 rÞnðj; kÞ ð2:11Þ
and in analogy with (2.10) we consider
dnðj þ v
j
d ; j; kÞ¼
1
2
rnðj; kÞ: ð2:12Þ
Then (2.5) is identically satisfied and all the quantities in (2.4) are defined. In this case, dnðj þ v
j
d ; j; kÞ is
not anymore a random variable but its value is uniquely determined by (2.12) and coincides with the mean
value of the corresponding random variable in GRW. Therefore we call this modified algorithm as a
‘‘deterministic’’ GRW (GRWD).
Another form of the GRW algorithm can be obtained by both preforming a deterministic scattering and
preserving the particles indivisibility. We use (2.6) and instead of (2.12) we introduce
dnðj þ v
j
d ; j; kÞ¼
n=2 if n is even;
½n=2þ h if n is odd;
ð2:13Þ
where n ¼ nðj; kÞ dnðj; j; kÞ, ½n=2 is the integer part of n=2 and h is a random variable taking the values 0
and 1 with probability 1/2. The quantity dnðj þ v
j
þ d ; j; kÞ is determined by (2.5). In comparison with
GRW, this algorithm reduces the particles number fluctuations at only one particle and we call it GRWR.
Since the fluctuations do not vanish, only the average of the GRWR solution is identical with the FD
solution. The algorithms without fluctuations (GRWD) and with reduced fluctuations (GRWR) can be
used to obtain numerical solutions for the advection–diffusion equation (1.8), as well as the stochastic
algorithm GRW. The GRW algorithm is expected to be more accurate when the fluctuations significantly
influence the simulated process [9].
The GRW algorithm and its modified forms GRWD and GRWR use the relation (2.4) where dnði; j; kÞ is
nonvanishing for every j satisfying j þ v
j
d ¼ i. Therefore, if V ðxÞ varies in space, the evolution of the
concentration in a node is obtained, unlike in (2.9), by contributions from more than the first neighboring
nodes. The terms in (2.4) are not apriori known, because they depend on the value of V in x
j
. In this case,
the GRW algorithm is no more equivalent with a FD scheme.
The implementation of the GRW algorithm as a computer code, encounters some problems related to
the computation of the Bernoulli distribution b
n
ðmÞ¼ 2
n
C
m
n
and of the corresponding repartition
F
n
ðmÞ¼
P
m
i¼0
b
n
ðiÞ. When the number of particles is of order 10
6
, the computation of b
n
ðmÞ or F
n
ðmÞ takes
too much time to be performed at each computation step. Therefore the values of F
n
ðmÞ are computed only
once and stored in files for values n ¼ 2
k
with 1 6 k 6 20. Due to the symmetry of F
n
ðmÞ with respect to
m ¼ n=2, only the values F
n
ðmÞ 6 0:5 are stored. If n < 2
21
, a binary representation n ¼
P
20
l¼0
aðlÞ2
l
is used.
The 2
l
particles of a group with aðlÞ 6¼ 0 are scattered in dn
l
ðj þ v
j
d ; j; kÞ and dn
l
ðj þ v
j
þ d ; j; kÞ, as
previously described, using a random number g uniformly generated in the interval ½0; 1. The final result is
obtained from
532 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
dnðj þ v
j
d ; j; kÞ¼
X
20
l¼0
aðlÞdn
l
ðj þ v
j
d ; j; kÞ:
If n P 2
21
, then there are several groups consisting of 2
20
particles and for each group the procedure from
above can be used. This method, referred to as GRW0 (first used in [23]), becomes time expensive for very
large n (see the ‘‘GRW0’’ curve in Fig. 9). In this paper a different method is used. The reduced variable
n ¼ðm n=2Þ=
ffiffiffiffiffiffiffi
n=4
p
and the repartitions F
n
ðnÞ are introduced. For n !1, the repartition F
n
ðnÞ tends to
the normal Gaussian repartition, according to De Moivre–Laplace theorem [15]. A good approximation is
obtained when for every n P 2
21
one uses the repartition corresponding to n ¼ 2
20
as function of the re-
duced variable n. For instance, the relative error of the values dn obtained using F
2
20 instead of F
2
30 is of the
order 10
9
. In this way, GRW can handle a number of particles equal to the maximum number of particles
that can be represented in the internal memory of the computer. For two and three-dimensional problems,
the GRW algorithms are implemented by performing the one-dimensional global scattering procedures
described in this section on x
1
, x
2
and x
3
space axes, according to the values of velocity components and
diffusion coefficients.
3. Simulation of the one-dimensional Gaussian diffusion
We verify the GRW algorithm described in the previous section for the solution of Eq. (1.8) corre-
sponding to V ðxÞ 0; DðxÞ 0:5 and initial condition lim
t!0þ
cðx; tÞ¼ N
0
dðxÞ, where dðxÞ is the Dirac
function. In this case the solution has a Gaussian analytic form and, from (1.5) and (1.9), it is given by
c
Gauss
ðx; tÞ¼ N
0
ð2ptÞ
1=2
exp
x
2
2t
: ð3:1Þ
Different numerical solutions obtained by GRW are quantitatively compared with the analytical solution
(3.1) in the space interval x 2 ½1; 1. The comparison is achieved at the time t
f
, when the number of
particles which left the interval ½1; 1 is 1% of the total number of particles. From the condition
1
N
0
Z
1
1
c
Gauss
ðx; t
f
Þ dx ¼ erf
1
ffiffiffiffiffi
2t
f
p
¼ 0:99;
we have t
f
¼ 0:15.
The numerical solution is obtained using the GRW, GRWD and GRWR algorithms, for v
i
¼ 0 and
d ¼ 1. Initially N fictitious particles are introduced in the origin of the space grid. We consider a sequence
of grids with the space steps dx ¼ð10gÞ
1
, where g ¼ 1; 2; ... ; 10. Since D ¼ 0:5, from (2.8) it follows that
the corresponding time steps are dt ¼ rdx
2
¼ rð10gÞ
2
and the numerical simulation contains k
f
¼
t
f
=dt ¼ 15g
2
=r time steps. To eliminate the boundaries effect on the numerical solution we must choose a
large enough grid so that no particles may reach the boundary at time t
f
. A particle covers the maximum
distance if it makes all the k
f
jumps in the same direction. Therefore the space grid must contain at least k
f
nodes on a side and on the other of the origin where all the N particles initially are located. Hence for
different values of r and g (dx) we obtain different numerical algorithms for D ¼ 0:5 and V ðxÞ¼ 0.
We want to compare c defined by (2.3) with the analytic solution (3.1), at time t
f
and over the spatial
interval ½1; 1. For r ¼ 1 some additional problems occur. Since initially all the N particles are at a single
node, if r ¼ 1 the simulation takes place on a single numerical mode, i.e., at even time steps the particles are
distributed only on the even nodes 2pdx and the odd nodes (2p þ 1Þdx are empty, and vice versa. So, at a
given time only half of the grid nodes contains particles and in (2.3) we must consider l ¼ 2, i.e., c ¼ n=ð2dxÞ.
The nodes corresponding to x ¼1 are always even, 1=dx ¼10g, and we introduce the quantity
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 533
i
f
¼
10g if k
f
is even;
10g 1 if k
f
is odd;
such that x
f
¼i
f
dx should be the node in the interval ½1; 1 which contains particles at t
f
and is the closest
to x ¼1. Then for r ¼ 1 we characterize the accuracy of the solution c with respect to the analytical
solution c
Gauss
by the norm kc c
Gauss
k, defined as
c k c
Gauss
k
2
¼
1
i
f
þ 1
X
i
f
i¼0
1
N
cðð2i
i
f
Þdx; t
f
Þ
1
N
0
c
Gauss
ðð2i i
f
Þdx; t
f
Þ
2
; ð3:2Þ
where only the nodes containing particles at t
f
are taken into account. If r < 1, then the two numerical
modes are mixed by the particles remaining at the same node. In this case one can use a formula analogous
to (3.2) with i
f
¼ 10g, summing over all ð2i
f
þ 1Þ nodes of the grid and with c ¼ n=dx, defined by (2.3) for
l ¼ 1.
First we analyze the numerical solution obtained with the GRWD algorithm described by (2.11) and
(2.12) with nði; kÞ real number. In the previous section it was proved that, for constant D and V 0,
GRWD is identical with the explicit FD scheme (2.9), hence we denote its solution by c
FD
ðx
i
; t
k
Þ. In Fig. 1,
the norm c
FD
c
Gauss
k k is represented as function of dx
2
for several values of the parameter r. The linear
behavior corresponds to convergence in dx
2
order of the FD scheme. From Fig. 2 it follows that the
maximum precision of the finite difference scheme is obtained for r 0:3.
In the following, besides d ¼ 1 and D ¼ 0:5, we also fix the parameter r ¼ 1. We denote by cðx
i
; t
k
Þ the
numerical solution (2.3) obtained with the GRW algorithm described in the previous section. In fact this
solution also depends on the spatial resolution dx ¼ð10gÞ
1
, the total number of the fictitious particles N ,
Fig. 1. Convergence of order dx
2
.
Fig. 2. Dependence on r.
534 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
and the number of simulations S , i.e., c ¼ cðx
i
; t
k
; dx; N ; S Þ. We successively analyze the dependence of the
solution accuracy on these parameters.
First for S ¼ 1 fixed, dx ¼ 0:1 (Fig. 3(a)) and dx ¼ 0:01 (Fig. 3(b)) we represent the norm (3.2) for the
numerical solution obtained with GRW and GRWR algorithms. Because S ¼ 1, there is a single simula-
tion, n is identical with n and because r ¼ 1, from (2.3), we have c ¼ n=ð2dxÞ. For a large enough number of
particles N , both GRW and GRWR approximates the analytical solution as well as the FD scheme. Since
the fluctuations in GRWR are reduced to minimum, this algorithm becomes equivalent to the FD scheme
for a smaller value of N . For this simple diffusion process the FD scheme is always more efficient, but for
complex diffusion processes as those presented in Section 5 other algorithms are needed.
To study the dependence of the precision on N , in Fig. 4 we represent for dx ¼ 0:1 the norm kc c
FD
k,
defined similarly to (3.2), where c corresponds to GRW, as a function of 1=
ffiffiffiffi
N
p
. One notices that the
dependence is linear. Then, Figs. 1 and 4 show that GRW converges to the analytical solution as Oðdx
2
Þ
þOð1=
ffiffiffiffi
N
p
Þ. For large N , when the condition 1=
ffiffiffiffi
N
p
¼ Oðdx
2
Þ is satisfied, the convergence order of GRW is
Oðdx
2
Þ, the same as for the finite differences scheme.
In Fig. 5 we analyze the dependence of norm (3.2) on N and S , for dx ¼ 0:1. When N ¼ 1 (the points
represented by } in Fig. 5), then the value kc c
Gauss
kffi 0:01, corresponding to FD scheme, is reached only
if S ffi 10
6
repetitions of the simulation are performed. In this case the GRW simulation is similar to those
produced by PT method, where the trajectory of one particle is generated S times. The same result is
Fig. 3. Dependence on N , for dx ¼ 0:1 (a) and dx ¼ 0:01 (b).
Fig. 4. GRW converges as 1=
ffiffiffiffi
N
p
to FD.
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 535
obtained in a single realization, S ¼ 1, for N ffi 10
6
particles handled with GRW (points represented by in
Fig. 5). This took only 0.5 s while for N ¼ 1 and S ¼ 10
6
(PT procedure) the computation time was of about
17 min. PT may be thought as an ‘‘analogical Monte Carlo method’’ where the diffusion is modeled by
averages over the diffusive behavior of individual particles. Since GRW algorithm can use large number of
particles, Monte Carlo repetitions are not necessary to obtain the required precision, and the computation
time is thousands times smaller.
4. Numerical boundary conditions
The boundary conditions for GRW algorithm depend on the values of d and v
i
. In this section we discuss
only the simplest case d ¼ 1 and v
i
¼ 0 analyzed in the previous section. In more complicated cases the
boundary conditions can be similarly derived by means of the methods presented in the following.
To formulate the boundary conditions we use the numerical flux of particles J ðx; tÞ defined as the number
of particles crossing at time t the coordinate x. We evaluate the numerical flux during a time step dt, so that
the value obtained should be assigned to the middle of the time step. For d ¼ 1 and v
i
¼ 0 the particles
jump only between the neighboring nodes, so that the numerical flux should be assigned to the middle of the
space step,
J ði þ 1=2; k þ 1=2Þ¼
1
dt
dnði ½ þ 1; i; kÞ dnði; i þ 1; kÞ: ð4:1Þ
For the GRWD algorithm, using (2.12) and (2.3) with l ¼ 1, the relation (4.1) becomes
J ði þ 1=2; k þ 1=2Þ¼ r
dx
2dt
cðx
i
; t
k
Þ ½ cðx
iþ1
; t
k
Þ:
From (2.8) it follows that this is the usual FD form of the Fick law
J ðx; tÞ¼Do
x
cðx; tÞ: ð4:2Þ
Let us consider a finite grid with 2L þ 1 nodes, i ¼L; L þ 1; ... ; L. Since v
i
¼ 0 and d ¼ 1, the
boundary conditions imply only the nodes i ¼L. We discuss only the boundary i ¼ L, the case i ¼L
being similar. A Dirichlet boundary condition can be formulated fixing the number of particles at the
boundary nðL; kÞ¼ n
b
ðkÞ, with n
b
ðkÞ a given function of time. For other boundary conditions, including
those of von Neumann type, we must evaluate the boundary flux
J ðL þ 1=2; k þ 1=2Þ¼
1
dt
dnðL ½ þ 1; L; kÞ dnðL; L þ 1; kÞ: ð4:3Þ
Fig. 5. Dependence on S.
536 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
In this formula dnðL þ 1; L; kÞ is determined by means of the GRW algorithm from nðL; kÞ, but the number
of particles dnðL; L þ 1; kÞ jumping from outside in the node i ¼ L is unknown. Therefore the boundary
condition can be formulated by calculating dnðL; L þ 1; kÞ such that the flux (4.3) would have the requested
value. Consider a von Neumann boundary condition
J ðL þ 1=2; k þ 1=2Þ¼ J
b
ðk þ 1=2Þ; ð4:4Þ
where J
b
is a given function of time. From (4.3) it follows the boundary condition
dnðL; L þ 1; kÞ¼ dnðL þ 1; L; kÞ dtJ
b
ðk þ 1=2Þ: ð4:5Þ
For J
b
ðk þ 1=2Þ¼ 0, from (4.5) we obtain the boundary condition for impermeable walls. For numerical
simulations with r ¼ 1 which use a single numerical mode, we must take into account that (4.5) mixes the
numerical modes. One can avoid the mixing of the numerical modes summing up (4.5) on two time steps.
Absorbing boundary condition corresponds to the case when all the particles leaving the grid are re-
moved, no particles being introduced from the exterior of the grid. In this case,
dnðL; L þ 1; kÞ¼ 0 ð4:6Þ
and the flux from the grid towards its exterior has the maximum value. The stationary boundary condition
imposes the equality of the fluxes on a side and the other of the boundary
J ðL þ 1=2; k þ 1=2Þ¼ J ðL 1=2; k þ 1=2Þ:
Using (4.1) and (2.5) we obtain the boundary condition
dnðL; L þ 1; kÞ¼ nðL; kÞ dnðL; L; kÞ dnðL; L 1; kÞ: ð4:7Þ
Applied to GRWD (which is equivalent to the FD algorithm), this condition becomes, by means of (2.11)
and (2.12), the ‘‘transmission boundary condition’’ [12]
nðL þ 1; kÞ¼ 2nðL; kÞ nðL 1; kÞ:
If we perform a numerical simulation on a finite grid of a nonstationary diffusion process in an un-
bounded domain after the particles reached the grid boundary, then we must formulate special boundary
conditions. A method to obtain such nonstationary conditions is to express the time derivative of the
concentration at the boundary by means of the time derivatives of the inside neighboring nodes. For
V ðxÞ¼ 0, it follows from (1.8)
o
t
cðx dx; t dtÞ¼ Do
2
x
cðx dx; t dtÞ¼ Do
2
x
½cðx; tÞ o
x
cðx; tÞdx o
t
cðx; tÞdt þ Oðdx
2
Þ
¼ o
t
cðx; tÞ Do
3
x
cðx; tÞdx þ Oðdx
2
Þ;
where we used the relation dt ¼ Oðdx
2
Þ derived from (2.8). Then we have
o
t
cðx; tÞ¼ o
t
cðx dx; t dtÞþ OðdxÞ: ð4:8Þ
Repeating the same argument for o
t
cðx 2dx; t dtÞ we also obtain
o
t
cðx; tÞ¼ 2o
t
cðx dx; t dtÞ o
t
cðx 2dx; t dtÞþ Oðdx
2
Þ: ð4:9Þ
For x ¼ Ldx and t ¼ðk þ 1Þdt, these relations written in finite difference, using (2.3), give the nonstationary
boundary condition
nðL; k þ 1Þ¼ nðL 1; k þ 1Þ nðL 1; kÞþ nðL; kÞ ð4:10Þ
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 537
and
nðL; k þ 1Þ¼ 2nðL 1; k þ 1Þ 2nðL 1; kÞ nðL 2; k þ 1Þþ nðL 2; kÞþ nðL; kÞ: ð4:11Þ
These conditions are expected to give an useful approximation when the particles distribution near the
boundary is a good approximation of the solution. But when the first particles approach the boundary,
there are significant fluctuations of the particles number. Therefore at the value obtained from (4.10) and
(4.11) supplementary conditions are imposed: (a) nðL; k þ 1Þ > 0; (b) nðL; k þ 1Þ must be smaller than the
value obtained from the impermeable wall condition (4.5) with J
b
¼ 0.
As an illustration of these boundary conditions we continue the simulation of the previous section for
temporal interval three times larger than t
f
. We use the GRWD algorithm defined by (2.11) and (2.12) with
the parameter r ¼ 0:3. In Fig. 6 we represent the time evolution of the norm (3.2) for four different
Fig. 6. Accuracy comparison.
Fig. 7. Concentrations at the boundary.
Fig. 8. Concentrations at t
f
.
538 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
boundary conditions. In these simulations dx ¼ 0:1 and the computational interval contains the spatial
interval ½1; 1. The points BC1 correspond to the impermeable wall boundary condition given by (4.5) with
J
b
¼ 0, BC2 and BC3 to the nonstationary boundary conditions (4.10) and (4.11), and BC4 to the absorbing
condition (4.6). The nonstationary condition (4.10) do not improve the accuracy with respect to (4.5) and
(4.6). The condition (4.11) keeps the norm under 0.02, which proves that it is suitable to be used in sim-
ulation of nonstationary diffusion. The same conclusion can be drawn from the evolution of the boundary
concentration (Fig. 7) and from the spatial variation of the concentration at the final time (Fig. 8).
5. Numerical simulation of diffusion in random velocity field
The transport in a random velocity field is a complex process consisting in diffusive movements of
particles and their transport along the stream lines of the velocity field. A mathematical description of this
process is given by Eq. (1.8) where D is a local diffusion coefficient and V ðxÞ a random field. Such models
are used to describe the transport of pollutants in natural porous media. In this section we consider the
numerical simulation of the transport of a contaminant substance in a saturated aquifer, for a punctual
injection case. The heterogeneity of the advection velocities is described by realizations of a random field.
Under these conditions, traditional methods (finite difference/element) are restricted at simplified aquifer
models [21]. Better results in simulation of field experiments were obtained by means of the stochastic
models based on PT method [20,22]. For instance, when the simulations of the contaminant transport at
field scale are performed with PT, the numerical diffusion and dispersion problems occurring in finite el-
ement/difference methods are completely eliminated [13].
In PT algorithm, the diffusion process is described by the movement of an ensemble of fictitious particles
in continuous space. For each particle the change of the position x in the time interval ðt
kþ1
t
k
Þ, due to the
realization V ðxÞ of the random velocity field and the local diffusion with coefficient D, is described, in one-
dimensional case by the discrete form of the It^ o equation
xðt
kþ1
Þ xðt
k
Þ¼ V ðxðt
k
Þ; t
k
Þdt þ G
ffiffiffiffiffiffiffiffiffiffi
2Ddt
p
; ð5:1Þ
where G is a Gaussian random variable with mean zero and unit variance. For large number of particles
one expects that their number density gives an approximation of the concentration field cðx; t
k
Þ satisfying an
advection–diffusion equation [22]. This assertion is based on the relation between the It^ o equation and the
Fokker–Planck equation for the probability density of particle position [5]. According to (1.9) the prob-
ability density is proportional to the concentration, so that (1.8) is in fact the Fokker–Planck equation. The
accuracy of the solution strongly depends on the number of tracked particles. As it is mentioned in liter-
ature [19], and we also have seen in Fig. 5, where the first curve (N ¼ 1) corresponds to PT procedure, a
trade off should be made to reduce the computation time without affecting the accuracy. Although im-
provements of the algorithm have been made, the high computational costs of PT simulations is still an
open problem [8]. That is why sometimes a ‘‘semianalytical’’ evaluation of diffusive movement is used [17].
In [24] for instance, on the basis of some numerical tests, the last term in (5.1) was taken as 4
ffiffiffiffiffiffiffi
Ddt
p
.
The computational effort in the PT method is due to the fact that every particle is separately displaced
and all the trajectories must be stored. That is why the GRW method, where groups of particles are si-
multaneously displaced, saves time and memory. In the following we show that GRW allows a faster and
more complete simulation of the diffusion in random fields. To illustrate the advantage of our method when
large numbers of particles are necessary, the computing time for GRW was compared with the computing
time for a PT method (‘‘ParTrace’’ code described in [14]). The same problem was solved using the same
computing platform (Cray T3E). In Fig. 9 we present the simulation of an isotropic diffusion, with D ¼ 0:5,
into a cube the side of which consisted of 21 nodes, for 10 time steps and for different number of particles
injected at the initial moment into the center of the cube.
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 539
' fill-rule='evenodd' d='M433.9 983l-1.2-.2c.1-.5 .3-.9 .7-1.1c.3-.3 .8-.4 1.5-.4c.6 0 1.1 .1 1.4 .2c.3 .2 .5 .4 .6 .6c.1 .2 .2 .6 .2 1.2v1.6c0 .5 0 .8 .1 1c0 .2 .1 .5 .2 .7H436c0-.1 0-.2-.1-.4q0-.1 0-.1c-.3 .2-.5 .4-.8 .5c-.3 .1-.5 .1-.8 .1c-.6 0-1-.1-1.3-.4c-.3-.3-.5-.7-.5-1.1c0-.3 .1-.6 .2-.8c.2-.2 .4-.4 .6-.5c.3-.1 .6-.2 1.1-.3c.7-.1 1.1-.3 1.4-.4v-.1c0-.3-.1-.5-.2-.6c-.2-.1-.4-.1-.8-.1c-.2 0-.4 0-.6 .1c-.1 .1-.2 .3-.3 .5Zm1.9 1.1c-.2 .1-.5 .2-.9 .2c-.4 .1-.6 .2-.7 .3c-.2 .1-.3 .3-.3 .5c0 .2 .1 .3 .2 .5c.2 .1 .3 .2 .6 .2c.2 0 .4-.1 .7-.3c.1-.1 .2-.2 .3-.4c0-.1 .1-.4 .1-.7v-.3Z' class='g3'/%3E%0A%3Cpath fill-rule='evenodd' d='M318.9 888.8h-7.2v-3.1c0-.7 .1-1.3 .2-1.6c.1-.4 .4-.7 .7-.9c.3-.2 .7-.3 1.1-.3c.5 0 1 .2 1.3 .5c.4 .3 .6 .8 .7 1.4c.2-.3 .4-.6 .6-.8c.2-.2 .6-.5 1.2-.8l1.4-.9v1.8l-1.6 1c-.5 .4-.9 .6-1 .8c-.2 .1-.3 .2-.3 .4c-.1 .2-.1 .4-.1 .7v.3h3v1.5Zm-4.2-1.5v-1c0-.7 0-1.2-.1-1.3c0-.2-.1-.4-.3-.5c-.1-.1-.3-.1-.5-.1c-.2 0-.4 .1-.6 .2c-.1 .1-.2 .3-.3 .5c0 .2 0 .5 0 1.1v1.1h1.8Z' class='g3'/%3E%0A%3Cpath fill-rule='evenodd' d='M318.9 878.6V880h-5.2c.5 .5 .8 1.1 1 1.8h-1.2c-.1-.4-.4-.8-.7-1.2c-.3-.4-.7-.7-1.1-.9v-1.1h7.2Z' class='g3'/%3E%0A%3Cpath fill-rule='evenodd' d='M318.9 874v1.4h-5.2c.5 .5 .8 1.1 1 1.8h-1.2c-.1-.4-.4-.8-.7-1.2c-.3-.4-.7-.7-1.1-.9V874h7.2Z' class='g3'/%3E%0A%3C/svg%3E)
GRW needs less than one second and only one Cray computing node while the computing time for
‘‘ParTrace’’ linearly increases with the number of particles and more Cray nodes are required (we stopped
the computation at 10
9
particles and 256 Cray nodes). The middle curve in Fig. 9 corresponds to GRW0
algorithm, where no approximations of binary repartitions are used (see Section 2). In this case the
computing time still remains orders of magnitude smaller than that of PT but, for N > 10
9
the time in-
creases with N .
For the simulation of two-dimensional diffusion in random fields we consider an advection–diffusion
equation similar to (1.8) and use the same procedure to generate the field as in [18]. The velocity field is
divergence free and has the form Vðx
1
; x
2
Þ¼ U þ uðx
1
; x
2
Þ, where U ¼ðU ; 0Þ is the mean velocity and
uðx
1
; x
2
Þ the fluctuation. The correlation coefficients
R
ll
¼hu
l
ðx
1
; x
2
Þu
l
ðx
1
þ a; x
2
Þi=hu
l
ðx
1
; x
2
Þu
l
ðx
1
; x
2
Þi;
where l ¼ 1; 2 and hi denotes the average over the realizations of the random field, decay as shown in Figs.
10 and 11, corresponding to an exponential correlated logarithm of the hydraulic conductivity, with cor-
relation length k ¼ 1 m. For our simulation, the parameters are the same as those used in simulations
represented in Figs. 11 and 14 of [18], i.e., U ¼ 1 m/day, and local diffusion coefficient D ¼ 0:01 m
2
=day.
Because the movement due to velocity fluctuations is much larger than diffusive motion, only the
spreading of the plume over regions with different velocities is important and not the fluctuations of the
number of particles. That is why the simulations were performed with the reduced fluctuations algorithm
Fig. 9. Comparison of computing times
Fig. 10. Longitudinal correlation.
540 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
' fill-rule='evenodd' d='M428.1 767.3l-1.3-.2c.2-.5 .4-.9 .8-1.1c.3-.3 .8-.4 1.5-.4c.6 0 1.1 .1 1.4 .2c.3 .2 .5 .3 .6 .6c.1 .2 .2 .6 .2 1.2v1.6c0 .5 0 .8 0 1c.1 .2 .2 .5 .3 .7h-1.4c0-.1-.1-.2-.1-.4c0-.1-.1-.1-.1-.2c-.2 .3-.4 .4-.7 .5c-.3 .2-.6 .2-.9 .2c-.5 0-.9-.1-1.2-.4c-.3-.3-.5-.7-.5-1.1c0-.3 .1-.6 .2-.8c.2-.2 .4-.4 .6-.5c.3-.2 .6-.3 1.1-.3c.6-.2 1.1-.3 1.3-.4v-.1c0-.3 0-.5-.2-.6c-.1-.1-.3-.2-.7-.2c-.2 0-.4 .1-.6 .2c-.1 .1-.2 .2-.3 .5Zm1.8 1.1c-.1 .1-.4 .1-.8 .2c-.4 .1-.6 .2-.8 .3c-.1 .1-.2 .2-.2 .4c0 .2 .1 .4 .2 .5c.1 .2 .3 .3 .5 .3c.3 0 .5-.1 .8-.3c.1-.1 .2-.3 .3-.4c0-.2 0-.4 0-.7v-.3Z' class='g0'/%3E%0A%3Cpath fill-rule='evenodd' d='M311.5 682.3h-7.2v-3c0-.8 .1-1.4 .2-1.7c.2-.4 .4-.6 .7-.9c.3-.2 .7-.3 1.1-.3c.6 0 1 .2 1.4 .5c.3 .3 .6 .8 .6 1.4c.2-.3 .4-.5 .6-.7c.3-.3 .7-.5 1.2-.9l1.4-.8v1.7l-1.6 1c-.5 .4-.9 .7-1 .8c-.2 .1-.2 .3-.3 .4c-.1 .2-.1 .4-.1 .8v.3h3v1.4Zm-4.1-1.4v-1.1c0-.7-.1-1.1-.1-1.3c-.1-.2-.2-.3-.3-.4c-.2-.1-.3-.2-.6-.2c-.2 0-.4 .1-.5 .2c-.2 .1-.3 .3-.3 .6c0 .1 0 .4 0 1v1.2h1.8Z' class='g0'/%3E%0A%3Cpath fill-rule='evenodd' d='M310.2 671h1.3v4.8c-.5 0-.9-.2-1.4-.4c-.4-.3-1-.8-1.7-1.6c-.6-.6-1-1-1.2-1.1c-.2-.2-.5-.3-.8-.3c-.3 0-.5 .1-.7 .3c-.2 .1-.3 .3-.3 .6c0 .3 .1 .6 .3 .7c.2 .2 .5 .3 .9 .3l-.2 1.4c-.7-.1-1.3-.4-1.6-.8c-.3-.4-.5-.9-.5-1.6c0-.7 .2-1.3 .6-1.7c.4-.4 .8-.6 1.4-.6c.3 0 .6 .1 .9 .2c.3 .1 .6 .3 .9 .6c.3 .1 .6 .4 1 .9c.4 .4 .6 .7 .8 .8c.1 .1 .2 .2 .3 .3V671Z' class='g0'/%3E%0A%3Cpath fill-rule='evenodd' d='M310.2 665.6h1.3v4.9c-.5-.1-.9-.3-1.4-.5c-.4-.3-1-.8-1.7-1.6c-.6-.6-1-1-1.2-1.1c-.2-.2-.5-.3-.8-.3c-.3 0-.5 .1-.7 .3c-.2 .1-.3 .4-.3 .6c0 .3 .1 .6 .3 .7c.2 .2 .5 .3 .9 .3l-.2 1.4c-.7-.1-1.3-.3-1.6-.8c-.3-.4-.5-.9-.5-1.6c0-.7 .2-1.3 .6-1.7c.4-.4 .8-.6 1.4-.6c.3 0 .6 .1 .9 .2c.3 .1 .6 .3 .9 .6c.3 .1 .6 .4 1 .9c.4 .4 .6 .7 .8 .8c.1 .1 .2 .2 .3 .3v-2.8Z' class='g0'/%3E%0A%3C/svg%3E)
GRWR with r ¼ 1. The GRWR algorithm requires reduced computing resources with respect to GRW and
GRWD algorithms. Unlike in GRW algorithm, only one random number will be generated (at every time
step and when at a given node there is an odd number of particles). Because the indivisibility of particles is
preserved, the plume has a smaller extension than in the case of GRWD algorithm and, consequently,
smaller grids are necessary.
The computation domain is a rectangular grid with 10
7
nodes, 0 6 i 6 10 000 and 0 6 j 6 1000, the space
step is dx ¼ dy ¼ 0:1 m and the time step dt ¼ 0:5 day. In two different realizations of the velocity field,
N ¼ 10
10
particles were released at the point ði
0
; j
0
Þ¼ð50; 500Þ. We checked that for greater values of N the
simulation results remain unchanged. A period of 800 days was simulated, so that in this interval the
particles travel a mean distance of 800 correlation lengths. The computations were performed with a PC
(Pentium III, 600 MHz, 64 Mb RAM) and lasted about 3 h for each simulation. In Fig. 12 we represent, at
different times, the number of particles summed up over transversal sections of the plume, nðx
1
; tÞ, for one
realization of the velocity field.
In PT approach the diffusion in random fields is analyzed indirectly by means of the so called ‘‘effective
diffusion coefficient’’. For a given realization of the random field, the average of ðx
l
x
l
Þ
2
over the real-
izations of the local diffusion process, described by the second term from (5.1), defines the variance of
displacements r
2
ll
ðtÞ¼ x
2
l
x
l
2
. Let hr
2
ll
ðtÞi be its average over the realizations of the random field. The
process has asymptotic diffusive behavior if there is a finite limit called effective diffusion coefficient
lim
t!1
r
2
ll
ðtÞ
t
¼ 2D
ll
; ð5:2Þ
where D < D
ll
< þ1. When the limit (5.2) is infinite one says that the asymptotic behavior is superdiffusive.
Processes with diffusive behavior also have the self-averaging property
Fig. 12. Distribution of particles across the plume.
Fig. 11. Transversal correlation.
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 541
lim
t!1
r
2
ll
ðtÞ
t
¼ 2D
ll
; ð5:3Þ
i.e., the coefficient defined by a single realization is the same as that defined as ensemble average [3].
In Figs. 13 and 14 the evolution of D
11
ðtÞ and D
22
ðtÞ, for the two-dimensional GRWR simulations, are
represented for the two realizations of the random field. The coefficients D
11
and D
22
are closed to the
constant values predicted by theory presented in [18].
The basic assumption of stochastic approach of diffusion in random media is that for large times the
process behaves asymptotically diffusive and, correspondingly, the self-averaging property (5.3) holds. Then
one expects that the stochastic model may have predictive force for individual realizations of the porous
medium [2]. ‘‘A major motivation of stochastic analysis of transport in a heterogeneous porous medium has
been to derive an effective equation for the concentration that may be used to make decision in real life
contamination events’’ [10]. First, the existence of effective diffusion coefficients (5.2) and self-averaging
property (5.3) must be checked. Further, the fluctuations of concentration in individual realizations should
be studied, in order to find the superior bound which delimits the ‘‘unsafe zone’’ [11]. Both problems are
difficult to be solved and they imply high computing resources. In [18], to obtain the approach to as-
ymptotic value of longitudinal effective diffusion coefficient, 3200 particle trajectories were simulated with
PT for times corresponding at 5000 correlation lengths. In other studies [17], the reliability of the com-
putations is achieved using 20,000 trajectories (500 realizations of the field and 40 particles in each reali-
zation), up to 20 correlations lengths but the asymptotic regime is not reached. The computational effort
needed to obtain statistically significant results ‘‘has serious implication’’ when the behavior in a single
realization is studied [4]. The limitations imposed to maximum number of particles also make it difficult to
obtain simulations of concentration fields and concentration fluctuations. Although the use of GRW
Fig. 13. Longitudinal effective diffusion coefficient.
Fig. 14. Transversal effective diffusion coefficient.
542 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544
method requires large grids, unlike PT which is a grid free method, such simulations can be obtained using
large number of particles. The results from Figs. 12–14 indicates that GRW could be an alternative tool in
modeling contaminant transport in groundwater.
6. Conclusions
The GRW algorithm moves groups of particles, unlike the usual random walk algorithm where for each
particle a random number generation is necessary. Hence, important saving of memory and computing time
are achieved. The number of particles is limited only by the maximum number that can be represented on
each computing machine. For simple diffusion processes, GRW generalizes the FD scheme and for great
enough number of particles their precision is the same. In this case, the determinist algorithm GRWD is
identical with FD algorithm.
The GRWR algorithm, where the fluctuations of the number of particles are reduced to a single particle,
can be used, as well as the GRWD algorithm, when the fluctuations do not influence the process essentially.
A version of GRWR also can be used to avoid the approximation of the binary repartitions for n > 2
20
,
described in Section 2, if only groups with n 6 2
20
particles are scattered using binary repartitions while the
group exceeding this number is scattered by the determinist rule given by (2.12). Another GRW algorithm
can be obtained when the Gaussian repartition and reduced variable is used instead of binary repartitions.
The GRW method is useful for complex diffusion processes, when the usual finite difference/element
methods become inefficient. In this paper we have presented the simulation of diffusion in random velocity
field. While the usual simulations uses tens of thousand of particles, the GRWR presented in this paper
contains 10
10
particles. This enables us to get informations about the asymptotic diffusive behavior and
concentration fluctuations. When no effective diffusion coefficients are computed and we are mainly in-
terested in concentration simulations, the boundary conditions discussed in Section 4 make it possible to
save more computing time and memory, because smaller grids can be used.
GRW can describe other complex diffusion processes. As an example, the simulation presented in
Section 5 can be completed for 3-dimensions, random space variable coefficient DðxÞ, chemical reactions
between several species of particles and radioactive decay. The modeling of chemical reactions and ra-
dioactive decay can be achieved if, at each time step, the variation of the number of particles is introduced
as described, respectively, by the kinetic formula for chemical reaction and the law of radioactive decay.
Acknowledgements
The authors are grateful to Dr. Jan Vanderborght and Dr. Olaf Nitzsche for fruitful discussions. C.V.
and N.S. are also grateful to Prof. M. Peculea for the initial impulse in the numerical simulation of dif-
fusion.
References
[1] W.F. Ames, Numerical Methods for Partial Differential Equations, second ed., Academic Press, New York, 1977.
[2] S. Attinger, M. Dentz, H. Kinzelbach, W. Kinzelbach, Temporal behavior of a solute cloud in a chemical heterogeneous porous
medium, J. Fluid Mech. 386 (1999) 77.
[3] M. Avellaneda, A.J. Majda, Superdiffusion in nearly stratified flows, J. Stat. Phys. 69 (3/4) (1992) 689.
[4] A. Bellin, P. Salandin, A. Rinaldo, Simulation of dispersion in heterogeneous porous formations: statistics, first-order theories,
convergence of computations, Water Resour. Res. 28 (9) (1992) 2211.
[5] C.W. Gardiner, Handbook of Stochastic Methods, Springer, New York, 1983.
C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544 543
[6] S.K. Godunov, V.S. Ryabenkii, Difference Schemes: An Introduction to the Underlying Theory, North-Holland, Amsterdam,
1987.
[7] A.J. Chorin, Vortex sheet approximation of boundary layers, J. Comput. Phys. 27 (1978) 428.
[8] A.F. Ghoniem, F.S. Sherman, Grid-free simulation of diffusion using random walk methods, J. Comput. Phys. 61 (1985) 1.
[9] W. Horsthemke, R. Lefever, Noise-induced Transitions. Theory and Applications in Physics, Chemistry and Biology, Springer,
Berlin, 1984.
[10] C. Kapoor, L.W. Gelhar, Transport in three-dimensionallity heterogeneous aquifers 1. Dynamics of concentration fluctuations,
Water Resour. Res. 30 (6) (1994) 1775.
[11] C. Kapoor, L.W. Gelhar, Transport in three-dimensionallity heterogeneous aquifers 2. Prediction and observations of
concentration fluctuations, Water Resour. Res. 30 (6) (1994) 1789.
[12] W. Kinzelbach, Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser, Oldenbourg Verlag,
M€ unchen, 1992.
[13] G.L. Moltyaner, M.H. Klukas, C.A. Willis, R.W.D. Killey, Numerical simulations of twin lake natural-gradient tracer tests: A
comparison of methods, Water Resour. Res. 29 (10) (1992) 3443.
[14] O. Neuendorf, Numerische 3D-Simulation des Stofftransport in einem heterogenen Aquifer. Ph.D. Thesis, J€ ul-3421 Forschungs-
zentrum J€ ulich, 1997.
[15] A. Papoulis, Probability, Random Variables and Stochastic Processes, McGraw-Hill, New York, 1991.
[16] K. Roth, K. Hammel, Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with
steady state flow, Water Resour. Res. 32 (6) (1996) 1653.
[17] P. Salandin, V. Fiorotto, Solute transport in highly heterogeneous aquifers, Water Resour. Res. 34 (5) (1998) 949.
[18] H. Schwarze, U. Jaekel, H. Vereecken, Estimation of macrodispersivity by different approximation methods for flow and
transport in randomly heterogeneous media, Transp. Porous Media 43 (2001) 265.
[19] Ne.-Z. Sun, Mathematical Modeling in Groundwater Pollution, Springer, New York, 1996.
[20] A.F.B. Tompson, L.W. Gelhar, Numerical simulation of solute transport in three-dimensional, randomly heterogeneous porous
media, Water Resour. Res. 26 (10) (1990) 2541.
[21] A.F.B. Tompson, R.B. Knapp, Reactive Geochemical Transport Problems in Nuclear Waste Analyses, Preprint UCRL-00552,
Lawrence Livermore National Laboratory, Livermore, 1989.
[22] A.F.B. Tompson, R.D. Falgout, S.G. Smith, W.J. Bosl, S.F. Asby, Analysis of subsurface contaminant migration and
remediation using high performance computing, Adv. Water Resour. 22 (3) (1998) 203.
[23] C. Vamos ß, N. Suciu, H. Vereecken, O. Nitzsche, H. Hardelauf, Global random walk simulations of diffusion, in: W. Kramer,
J.W.V. Gudenberg (Eds.), Scientific Computing, Validated Numerics, Interval Methods, 343, Kluwer Academic Publishers,
Dordrecht, 2001.
[24] R. Zhang, K. Huang, M.T. van Genuchten, An efficient Eulerian–Lagrangian method for solving solute transport problems in
steady and transient flow fields, Water Resour. Res. 29 (12) (1993) 4131.
544 C. Vamos ß et al. / Journal of Computational Physics 186 (2003) 527–544