EVALUATION OF OVERSHOOTING ERRORS IN PARTICLE METHODS FOR DIFFUSION BY BIASED GLOBAL RANDOM WALK

. The adjustment of grid steps which guarantees that particles methods yield no numerical diﬀusion inevitably induces overshooting errors in the solution of the parabolic partial diﬀerential equations with space variable coeﬃcients. In this paper we give an evaluation of the overshooting errors of the “global random walk” algorithm (GRW), a computational eﬃcient method used in simulations for transport in environmental problems. The evaluation is performed by comparisons between the GRW solutions and those of the “biased global random walk” algorithm (BGRW), a cellular automaton, which is computational more expensive but is also free of overshooting errors. The reference problem was the diﬀusive transport in a random velocity ﬁeld, a model for the transport of the contaminant solutes in groundwater. The evaluation reveals that, for an optimum choice of the parameters, GRW results for time intervals of practical interest lie in ranges of acceptable precision, for both the ensemble averaged observables and for their ﬂuctuations.


INTRODUCTION
Particles methods supply stable solutions, free of numerical diffusion, for parabolic equations with highly variable parameters.These methods are successfully used in environmental problems, characterized by high Péclet numbers and irregular space variable velocity and diffusion coefficients [5].An often used method is the "particle tracking" (PT), which is a forward-time Euler scheme for the Itô stochastic equation [4].Itô equation describes the trajectory of a diffusing particle and provides a description of the advection-diffusion process which is equivalent to the Fokker-Plank equation for the corresponding probability densities [1].In particular, for the onetime density, which can be identified with the concentration c(x, t), the Fokker-Plank equation coincides with the advection dispersion parabolic equation ( 1) where V is the velocity field and D is the diffusion tensor.Note that (1) is the Itô form of the Fokker-Planck equation, the only one which is equivalent to Itô stochastic equation and for which PT provides correct numerical solutions.More familiar and used in practice is the Stratonovich form (1') which is the local mass balance equation derived from the Fick's constitutive law (i.e., concentration flux ∼ −D∇c) [3].The equation (1') takes the form (1) and can be associated to Itô equation, if the velocity V * is supplemented with a drift vector with components given by the divergence of the rows of the diffusion tensor: Due to the equivalence between the Itô and Fokker-Planck descriptions of the transport process, the PT solution of (1) at a given time t = kδt and at a grid point where N is the total number of particles and ∆ l = 2s l δx l , l = 1, 2, 3, are the lengths of symmetrical intervals centered at x l .The severe limitation of PT concerns the total number of particles, which in practice cannot be large enough to ensure the desired accuracy of c [5,6,9].
The GRW algorithm is equivalent to a superposition of many PT procedures [9].Starting with a given distribution of N particles in a computational grid, all the particles lying at a grid site are simultaneously spread, first by an advection displacement, then by unbiased diffusion jumps.For illustration we present the GRW algorithm in the case of constant D, when the two-and three-dimensional algorithms can be designed by repeating a one-dimensional procedure for all space directions.The onedimensional GRW algorithm describes the scattering of n(i, k) particles from (x i , t k ) by (3) n(j, k) = δn(j, j + v j , k) + δn(j + v j − d, j, k) + δn(j where v j are discrete displacements in a given velocity field and d is the amplitude of unbiased diffusive jumps.The quantities δn are Bernoulli random variables and describe respectively, the number of particles which remain at the same grid site j +v j reached after an advective displacement, the number of particles jumping to the left and those jumping to the right (with respect to the advected position j + v j ).The mostly used version is the so called "reduced fluctuations" version of GRW, defined by is the integer part of n/2 and θ is a random variable taking the values 0 and 1 with probability 1/2.When applied in large scale problems, the reduced fluctuations GRW algorithm has two advantages.First, a much smaller number of calls of the random number generator is required, which significantly reduces the computational costs.Then, the diffusion front does not extend beyond the limit concentration defined by one particle at a grid point, keeping a physical significant shape (compared to the finite differences schemes, where a pure diffusion front has a cubic shape of size ∼ (2Dt) 1/2 ).
The distribution of the particles at the next time (k + 1)δt is given by The averages of the terms in (1) over an ensemble of simulations are related by ( 4) where 0 ≤ r ≤ 1.The diffusion coefficient D is expressed in terms of the grid steps by the relation For r ≤ 1, there is no numerical diffusion.Moreover, since the total number of particles N contained in the grid is conserved, the GRW algorithm is stable.GRW can be regarded as a particular cellular automaton (CA), i.e. it is a stochastic process in the space of configurations, defined at a given time by the occupation numbers at each grid site.In the GRW algorithm the number of particles per grid site is not limited by an "exclusion principle".Therefore, GRW is "self-averaging", in the sense that the solution given by a single simulation is in fact the same as that obtained after averaging over large ensembles of simulations.For instance, the GRW solution of the heat equation converges to the Gaussian distribution as O(δx 2 ) +O(1/ √ N ); for large numbers of particles the convergence order is O(δx 2 ), the same as for the finite differences scheme [9].Though if GRW overcomes the limitation of PT concerning the number of particles, the correctness of its solution is still affected by overshooting errors, which characterize particles methods.Overshooting errors occur when the particles jump over more than one lattice site and the variability of the velocity field between the ends of the jumps is not accounted for.Overshooting can be somehow kept under control (with increasing computational costs) by a fine discretization but cannot be completely removed [6].Therefore, to assess the reliability of GRW solutions, an evaluation of the overshooting errors by comparisons with other methods is necessary.

BGRW ALGORITHM
To prevent overshooting, one imposes that particles jump only on the first-neighbor grid sites.This can be done by fixing to unity the amplitude of jumps in GRW, d = 1, and by modeling the advection by a bias in the random walk jumps.The new cellular automaton algorithm is called "biased global random walk" (BGRW) [7].We describe here the one-dimensional BGRW algorithm, for variable velocity and diffusion coefficients.In BGRW, the relations (3-5) are replaced by Thus unlike in the unbiased GRW algorithm, where the advective displacements v j and the unbiased diffusive jumps d are independent, in BGRW the advective flux v j n(j, k) is given, as follows from (4'), by the bias of the diffusive jumps The configuration of the grid, at the time (k + 1)δt is obtained by summing the contributions from the first-neighbor sites, Taking the averages of this relation over an ensemble of simulations and using (4') and (5') one obtains where n(x, t) = n(iδx, kδt).Thus, BGRW is equivalent to the forward-time centredspace finite difference scheme for the 1-dimensional version of the advection-diffusion equation (1).For large N , n(x, t) n(x, t), i.e.BGRW is also self-averaging.As it follows from (4'), BGRW is subjected to the restrictions The first condition (6) corresponds to the von Neumann criterion for stability, implying that there is no numerical diffusion.The second condition in (6) ensures that the Courant numbers are sub-unitary, then the algorithm also avoids the overshooting errors.Implementations of BGRW for higher dimensional cases are easily done by using vector parameters, for instance r = (r 1 , r 2 , r 3 ), v = (v 1 , v 2 , v 3 ), and accordingly modifying (3'-5') [7].

EVALUATION OF GRW BY COMPARISON WITH BGRW
We consider an isotropic two-dimensional diffusion in groundwater (D 1 = D 2 = D = 0.01 m 2 /day) in a velocity field V with mean U = 1 m/day, oriented along the x 1 axis and with a standard deviation σ = 0.2 m/day.The velocity field is generated as a realization of a periodic random field, consisting of a superposition of 64 sin modes which approximates a Gaussian field; this Gaussian field is in turn an approximation for small fluctuations of the hydraulic conductivity of the Darcy velocity field in saturated groundwater formations [8].For all simulations presented in this paper, we fixed N = 10 10 , the total number of particles which ensures the self-averaging of the GRW simulations for this transport problem [6].
Since D is constant we use in BGRW algorithm δx 1 = δx 2 = δx and r 1 = r 2 = r = 0.5.Supposing that the maximum velocity can be as large as V max = U + 5σ = 2 m/day, from the second relation (5') and the second condition (6) it follows that δx 2D/V max = 0.01 m.Correspondingly, from the first relation in (5') one obtains δt = 0.0025 days.The BGRW simulation of the transport over 100 days, for a point instantaneous injection at the origin of the lattice, requires about 15 cpu hours.For the same problem and consuming the same cpu time, the (unbiased) GRW algorithm can perform the simulation of the transport over 4 000 days, when an optimum choice of parameters, (δx = 0.1 m, δt = 0.5 day, r = 0.25) is used [8].Therefore, GRW is more appropriate for investigations on the large time behavior of transport processes which show complex scaling properties, like the fate of contaminants in groundwater [8].To assess the reliability of the computational efficient GRW simulations we shall use the overshooting-free (and more expensive in terms of cpu time) BGRW algorithm.Because of the statistical nature of the predictions for groundwater contamination, we go beyond the single realization comparisons done in [7], and proceed to an evaluation of GRW solutions for the means and fluctuations, with respect to an ensemble of 256 velocity realizations, of the mostly used observables.The latter are the moments of the solute concentration and the space average of the concentration over the crosssection of the solute plume.
The two-dimensional version of the equation ( 1) was solved by GRW and BGRW, for identical realizations of the velocity, point instantaneous injection at the origin, and no-boundary conditions.The latter were achieved by fixing the grid dimensions to L 1 = 150 m and L 2 = 20 m, so that during the total simulation time, T = 100 days, no particle reaches the boundary.The concentration c was computed by (2) for ∆ 1 = ∆ 2 = ∆ = 1 m.The 1st and the 2-nd moments of the concentration c are defined by (7) where l = 1, 2 and the integrals are computed over the support of c.For comparisons we used the derivative of the 1-st moments V cm l (t) = dµ l (t)/dt, which represent the velocity components of the center of mass of the solute body, and the rate of increase with time of the 2-nd moments, which is used to define effective diffusion coefficients D ef f ll (t) = µ ll (t)/(2t) [8].The cross-section concentration was computed by ( 8) where L 2 is the transverse dimension of the grid.For comparisons we used the crosssection concentration at the center of mass of the plume, i.e. (8) estimated at x 1 = µ 1 (t).
The evaluation of the center of mass velocity was done using the absolute errors ( 9) where Φ stands, respectively, for the expectation E(V cm l ), computed by arithmetic means over the 256 velocity realizations, and for the corresponding standard deviations . The results for l = 1 and l = 2 are presented in figure 1.In these figures we also plotted, by horizontal lines, the mean errors (9) calculated by formula , where T 1 = 1 day.
The evaluation for the effective coefficients D ef f ll and cross-section concentration C were achieved by using the percentage relative errors  where Φ stands, again, for the corresponding expectations, E(•), and standard deviations, SD(•).The results for the longitudinal and transverse effective coefficients are given in figure 2, respectively, and those for cross-section concentration in figure 3. The horizontal lines in these figures correspond to the mean errors (11), calculated similarly to (10), as ε(Φ) .

CONCLUSIONS
A comparison between figure 1 and the theoretical expectation values E(V cm 1 ) = 1 m/day and E(V cm ) = 0 m/day, shows that GRW reproduces the mean and the fluctuations of the velocity of the plume center of mass with a very good precision of the order of a few cm/day.Figures 2a and 3 also show that the longitudinal effective coefficient and cross-section concentration are simulated with GRW with a satisfactory precision of about 5%.The errors for the transverse effective coefficient are larger, mainly those for the standard deviation (figure 2b).But in modelling groundwater contamination the most important effects are described by the longitudinal effective coefficients.Moreover, the existing limit theorems for large time behavior of the transport also predict up-scaled diffusion coefficients given by asymptotic expansions truncated at the order of local diffusion coefficient D. After tens of days the transverse effective coefficient already reaches the theoretical asymptotic value D ef f 22 ∼ D and its standard deviation is one order of magnitude smaller [8].Because the relative errors presented in figure 2b are smaller than 100%, from (11) it follows that the absolute errors in GRW simulations are smaller than D, thus still acceptable for inclusion into the prediction process.Therefore, we can say that GRW is accurate enough for the purpose of stochastic analyses of large time behavior of the transport.

Fig. 1 .
Fig. 1.The evolution of the center of mass velocity.
Transverse effective coefficient.

Fig. 2 .
Fig. 2. The evolution of the effective coefficient.