The Global Random Walk algorithm (GRW) is equivalent to a superposition of many particle tracking procedures. Starting with a given distribution of $N$$N$NN$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, as shown in Fig. 1. GRW is 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 lattice site. In the GRW algorithm the number of particles per grid site is not limited by an "exclusion principle" and there are no limitations as to the total number of particles. Therefore, GRW is "self-averaging"in the

sense that the solution given by a single simulation is practically 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 $\mathcal{O}\left(\delta x^2\right)+\mathcal{O}(1/\sqrt{N})$$\mathcal{O}\left(\delta x^2\right)+\mathcal{O}(1/\sqrt{N})$O(deltax^(2))+O(1//sqrtN)\mathcal{O}\left(\delta x^{2}\right)+\mathcal{O}(1 / \sqrt{N})$\mathcal{O}\left(\delta x^2\right)+\mathcal{O}(1/\sqrt{N})$, where $\delta x$$\delta x$delta x\delta x$\delta x$ is the lattice parameter and $N$$N$NN$N$ is the total number of particles. Thus, for large $N$$N$NN$N$ the convergence order is $\mathcal{O}\left(\delta x^2\right)$$\mathcal{O}\left(\delta x^2\right)$O(deltax^(2))\mathcal{O}\left(\delta x^{2}\right)$\mathcal{O}\left(\delta x^2\right)$, the same as that for the finite differences scheme. By working with integers, GRW is free of round-of errors, avoids numerical diffusion and it is inherently stable [1]. However, for variable drift and diffusion coefficients overshooting errors occur when the particles jump over more than one lattice site (see Fig. 1).
To get rid of overshooting errors, we impose that particles jump only to the nearest sites (Fig. 2). In this procedure the advection will be simulated by a bias in the random walk jumps. Therefore, we call it "biased global random walk" (BGRW) algorithm. Since BGRW moves all the particles of one lattice site in a single numerical procedure, $N$$N$NN$N$ can be as large as necessary to ensure the self-averaging, which is the main difference with respect to other CA for diffusion without exclusion principle [2].

The 2-dimensional BGRW is defined by the CA rule

where $n(i,j,k)$$n(i,j,k)$n(i,j,k)n(i, j, k)$n(i,j,k)$ is the number of particles at the site $(x,y)=(i\delta x,j\delta y)$$(x,y)=(i\delta x,j\delta y)$(x,y)=(i delta x,j delta y)(x, y)=(i \delta x, j \delta y)$(x,y)=(i\delta x,j\delta y)$ at the time $t=k\delta t$$t=k\delta t$t=k delta tt=k \delta t$t=k\delta t$. All the terms of (1) are Bernoulli random variables, computed as in the "reduced fluctuations GRW" presented in ref. [1]. Corresponding to the components of the drift (velocity) and diffusion coefficients of the transport problem, $V_x(x,y,t),V_y(x,y,t),D_x(x,y,t)$$V_x(x,y,t),V_y(x,y,t),D_x(x,y,t)$V_(x)(x,y,t),V_(y)(x,y,t),D_(x)(x,y,t)V_{x}(x, y, t), V_{y}(x, y, t), D_{x}(x, y, t)$V_x(x,y,t),V_y(x,y,t),D_x(x,y,t)$ and
$D_y(x,y,t)$$D_y(x,y,t)$D_(y)(x,y,t)D_{y}(x, y, t)$D_y(x,y,t)$, we define the dimensionless parameters

The averages of the terms in (1) over an ensemble of simulations are related by

Defining the particle density $\rho (x,y,t)=\overline{n}(i,j,k)$$\rho (x,y,t)=\overline{n}(i,j,k)$rho(x,y,t)= bar(n)(i,j,k)\rho(x, y, t)=\bar{n}(i, j, k)$\rho (x,y,t)=\overline{n}(i,j,k)$ and summing the contributions from the first neighbors to a lattice site, from (1-3) one obtains

which is just the forward-time centred-space finite difference scheme for the advectiondiffusion (Fokker-Plank) equation

The equation which corresponds to Fick's law

is equivalent to (5) if the new drift coefficients are given by the relations $V_x^{\ast }=V_x-{\partial }_x{D}_x$$V_x^{\ast }=V_x-{\partial }_x{D}_x$V_(x)^(**)=V_(x)-del_(x)D_(x)V_{x}^{*}=V_{x}-\partial_{x} D_{x}$V_x^{\ast }=V_x-{\partial }_x{D}_x$ and $V_y^{\ast }=V_y-{\partial }_y{D}_y$$V_y^{\ast }=V_y-{\partial }_y{D}_y$V_(y)^(**)=V_(y)-del_(y)D_(y)V_{y}^{*}=V_{y}-\partial_{y} D_{y}$V_y^{\ast }=V_y-{\partial }_y{D}_y$, of which the corresponding BGRW can be easily derived.
As it follows from (3), BGRW is subject to the restrictions

Adding the conditions $r_x\le 0.5$$r_x\le 0.5$r_(x) <= 0.5r_{x} \leq 0.5$r_x\le 0.5$ and $r_y\le 0.5$$r_y\le 0.5$r_(y) <= 0.5r_{y} \leq 0.5$r_y\le 0.5$, the von Neumann criterion for stability is satisfied, implying that there is no numerical diffusion. The last two inequalities in (6) ensures that the Courant numbers are sub-unitary, thus the algorithm also avoids the overshooting errors.

As a direct consequence of (6), we can see that removing overshooting errors requires high computational costs. Let us consider an isotropic two-dimensional diffusion in groundwater ( $D_x=D_y=D=0.01{\mathrm{m}}^2/$$D_x=D_y=D=0.01{\mathrm{m}}^2/$D_(x)=D_(y)=D=0.01m^(2)//D_{x}=D_{y}=D=0.01 \mathrm{~m}^{2} /$D_x=D_y=D=0.01{\mathrm{m}}^2/$ day ) in a mean flow of $U=1\mathrm{m}/$$U=1\mathrm{m}/$U=1m//U=1 \mathrm{~m} /$U=1\mathrm{m}/$ day , oriented along the $x$$x$xx$x$ axis and with a standard deviation $\sigma =0.2\mathrm{m}/\mathrm{day}$$\sigma =0.2\mathrm{m}/\mathrm{day}$sigma=0.2m//day\sigma=0.2 \mathrm{~m} / \mathrm{day}$\sigma =0.2\mathrm{m}/\mathrm{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 [4]. Admitting that the maximum velocity can be as large as $V^{\text{max}}=U+5\sigma =2\mathrm{m}/\mathrm{day}$$V^{\text{max }}=U+5\sigma =2\mathrm{m}/\mathrm{day}$V^("max ")=U+5sigma=2m//dayV^{\text {max }}=U+5 \sigma=2 \mathrm{~m} / \mathrm{day}$V^{\text{max }}=U+5\sigma =2\mathrm{m}/\mathrm{day}$, from (2) and the second condition (6) it follows that $\delta x\le 2D/V_x^{\text{max}}=0.01\mathrm{m}$$\delta x\le 2D/V_x^{\text{max }}=0.01\mathrm{m}$delta x <= 2D//V_(x)^("max ")=0.01m\delta x \leq 2 D / V_{x}^{\text {max }}=0.01 \mathrm{~m}$\delta x\le 2D/V_x^{\text{max }}=0.01\mathrm{m}$. Since this space step also fulfils the third condition (6), in the following we take $\delta y=\delta x$$\delta y=\delta x$delta y=delta x\delta y=\delta x$\delta y=\delta x$. Correspondingly, from (2), $\delta t=0.0025day$$\delta t=0.0025day$delta t=0.0025 day\delta t=0.0025 d a y$\delta t=0.0025day$ (the case represented in Fig. 2). The simulation of the transport over 50 days, for a point instantaneous injection at the origin of the lattice, requires about 3 cpu hours. For the same problem and consuming the same cpu time, the unbiased GRW algorithm (UGWR) with $\delta x=0.1m$$\delta x=0.1m$delta x=0.1 m\delta x=0.1 m$\delta x=0.1m$ and $\delta t=0.5day$$\delta t=0.5day$delta t=0.5 day\delta t=0.5 d a y$\delta t=0.5day$ (Fig. 1) can perform the simulation of the transport over 1000 days. Nevertheless, the BGRW simulations are very helpful for the evaluation of other numerical methods, mainly, as in the case presented here, when no analytical solutions are available. We computed the 1 st and 2 -nd centered moments of the density $\rho$$\rho$rho\rho$\rho$, defined by

where $\alpha$$\alpha$alpha\alpha$\alpha$ stands for $x$$x$xx$x$ or $y$$y$yy$y$ and the integrals are computed over the support of $\rho$$\rho$rho\rho$\rho$. Further, using (7), we computed the derivatives of the 1 -st moments ${\mathcal{V}}_{\alpha }=d{\mu }_{\alpha }/dt$${\mathcal{V}}_{\alpha }=d{\mu }_{\alpha }/dt$V_(alpha)=dmu_(alpha)//dt\mathcal{V}_{\alpha}=d \mu_{\alpha} / d t${\mathcal{V}}_{\alpha }=d{\mu }_{\alpha }/dt$, which represent the velocity components of the center of mass of the solute body, and the rates of increase with time of the 2-nd moments ${\mathcal{D}}_{\alpha \alpha }={\mu }_{\alpha \alpha }/(2t)$${\mathcal{D}}_{\alpha \alpha }={\mu }_{\alpha \alpha }/(2t)$D_(alpha alpha)=mu_(alpha alpha)//(2t)\mathcal{D}_{\alpha \alpha}=\mu_{\alpha \alpha} /(2 t)${\mathcal{D}}_{\alpha \alpha }={\mu }_{\alpha \alpha }/(2t)$, which in the large time limit define the effective diffusion coefficients for this transport problem.
The self-averaging of the GRW simulations for the transport problem considered in this paper is ensured if the total number of particles is of the order $N={10}^{10}$$N={10}^{10}$N=10^(10)N=10^{10}$N={10}^{10}$ [4]. Using this value of $N$$N$NN$N$ in all cases, the numerical solution $\rho =\overline{n}$$\rho =\overline{n}$rho= bar(n)\rho=\bar{n}$\rho =\overline{n}$ was estimated by the actual number of particles $n$$n$nn$n$ at the lattice sites.
The moments (7) were computed with BGRW for the parameters $\delta x=0.01m$$\delta x=0.01m$delta x=0.01 m\delta x=0.01 m$\delta x=0.01m$ and $\delta t=0.0025day$$\delta t=0.0025day$delta t=0.0025 day\delta t= 0.0025 d a y$\delta t=0.0025day$ (case b1) and for a finer discretization, $\delta x=0.005m$$\delta x=0.005m$delta x=0.005 m\delta x=0.005 m$\delta x=0.005m$ and $\delta t=0.000625day$$\delta t=0.000625day$delta t=0.000625 day\delta t=0.000625 d a y$\delta t=0.000625day$ (case b2), with $r_x=r_y=r=0.5$$r_x=r_y=r=0.5$r_(x)=r_(y)=r=0.5r_{x}=r_{y}=r=0.5$r_x=r_y=r=0.5$ in both cases. The errors of BGRW simulation for the case (b1) are estimated by

where $\mathrm{\Delta }{\mathcal{V}}_{\alpha }$$\mathrm{\Delta }{\mathcal{V}}_{\alpha }$DeltaV_(alpha)\Delta \mathcal{V}_{\alpha}$\mathrm{\Delta }{\mathcal{V}}_{\alpha }$ and $\mathrm{\Delta }{\mathcal{D}}_{\alpha \alpha }$$\mathrm{\Delta }{\mathcal{D}}_{\alpha \alpha }$DeltaD_(alpha alpha)\Delta \mathcal{D}_{\alpha \alpha}$\mathrm{\Delta }{\mathcal{D}}_{\alpha \alpha }$ are the deviations of the corresponding quantities computed in case (b1) with respect to those obtained in case (b2) and $T$$T$TT$T$ is the simulation duration.

The estimations presented in Table 1 are orders of magnitude smaller than the fluctuations of the first two moments of the density $\rho$$\rho$rho\rho$\rho$ (governed by the physical parameters $D=0.01m^2/day$$D=0.01m^2/day$D=0.01m^(2)//dayD=0.01 m^{2} / d a y$D=0.01m^2/day$ and $\sigma =0.2m/day$$\sigma =0.2m/day$sigma=0.2 m//day\sigma=0.2 m / d a y$\sigma =0.2m/day$ ). A numerical investigation on the convergence of BGRW by comparisons with analytical solutions has not yet been done. However, since there are no overshooting errors, it is expected that the convergence order for BGRW is the same as in the case of the genuine diffusion (which was shown in [1] to be $\mathcal{O}\left(\delta x^2\right)$$\mathcal{O}\left(\delta x^2\right)$O(deltax^(2))\mathcal{O}\left(\delta x^{2}\right)$\mathcal{O}\left(\delta x^2\right)$ ). Since, due to conditions (6), this order is much smaller than for the particles methods with overshooting, BGRW solutions can serve as reference to evaluate the faster (but coarser) unbiased algorithms.
As an illustration, we compare in Figs. 3 and 4 the deviations $\mathrm{\Delta }{\mathcal{V}}_{\alpha }$$\mathrm{\Delta }{\mathcal{V}}_{\alpha }$DeltaV_(alpha)\Delta \mathcal{V}_{\alpha}$\mathrm{\Delta }{\mathcal{V}}_{\alpha }$ and $\mathrm{\Delta }{\mathcal{D}}_{\alpha \alpha }$$\mathrm{\Delta }{\mathcal{D}}_{\alpha \alpha }$DeltaD_(alpha alpha)\Delta \mathcal{D}_{\alpha \alpha}$\mathrm{\Delta }{\mathcal{D}}_{\alpha \alpha }$ with respect to BGRW (case b1) of the results given by UGRW for the sets of parameters $\delta x=0.1m$$\delta x=0.1m$delta x=0.1 m\delta x=0.1 m$\delta x=0.1m$, $\delta t=0.5day,r=0.25$$\delta t=0.5day,r=0.25$delta t=0.5 day,r=0.25\delta t=0.5 d a y, r=0.25$\delta t=0.5day,r=0.25$ (case u1) and $\delta x=0.01m,\delta t=0.1day,r=0.408$$\delta x=0.01m,\delta t=0.1day,r=0.408$delta x=0.01 m,delta t=0.1 day,r=0.408\delta x=0.01 m, \delta t=0.1 d a y, r=0.408$\delta x=0.01m,\delta t=0.1day,r=0.408$ (case u2). The corresponding error estimations via (8) are given in Table 2.

Even if the coarser discretization (u1) yields errors $\epsilon \left({\mathcal{D}}_{xx}\right)$$\epsilon \left({\mathcal{D}}_{xx}\right)$epsi(D_(xx))\varepsilon\left(\mathcal{D}_{x x}\right)$\epsilon \left({\mathcal{D}}_{xx}\right)$ of the order of $D$$D$DD$D$ it is still accurate enough to reproduce the behavior of the expectations (averages over ensambles of velocity fields). In this case UGRW can be successfully used in investigations on the large time behavior and self-averaging properties of the transport process [4]. Since for (u2) the errors are one order of magnitude smaller, in this case UGRW can be used to simulate the behavior in single realizations. But when higher accuracy is necessary (smaller times, reaction-diffusion processes) BGRW should be used.

The local CA character of the rule (1) makes the BGRW algorithm naturally parallel. Therefore, large simulations are possible on massively parallel computers. Moreover, since only a minimal amount of communication between physically neighboring processors is necessary, the parallel computing implementation of BGRW could result in a considerable reduction of the computing time.
The discrete stochastic process governing the movement of the particles on the BGRW lattice is a ratchet-like mechanism which induces advection from asymmetric fluctuations [3]. This provides a plausible physical description of the transport in natural porous media: at the pores scale the solute molecules move along the erratic stream lines of the water flow through the void space, under the forcing effect of an almost constant hydraulic pressure. The equivalence of BGRW with the finite difference scheme (4) of the partial derivative equation (5) also agrees with the experimental finding that, at a macroscopic scale, the transport can be described by an advection-diffusion equation.
Owing to the simplicity of the reaction-diffusion CA [2] and to the fact that the number of particles can be as large as the real number of molecules of various species involved in chemical reactions [1], the extension of the BGRW algorithm to reaction-diffusion processes appears to be a promising research direction in contaminant hydrology.

Acknowledgments: The research reported in this paper has been supported in part by the Deutsche Forschungsgemeinschaft grant SU 415/1-1 awarded to the first author.

[1] Vamoş, C., N. Suciu, and H. Vereecken: Generalized random walk algorithm for the numerical modeling of complex diffusion processes. Comp. Phys. 186(2) (2003), p. 527-544.
[2] Karapiperis, T. and B. Blankleider: Cellular automaton model of reaction-transport processes. Physica D 78 (1994), p. 30-64.
[3] Doering, C. R., W. Horsthemke, and J. Riordan: Nonequilibrium fluctuation-induced transport. Phys. Rev. Lett. 72(19) (1994), p. 2984-2987.
[4] Suciu N., C. Vamoş, J. Vanderborght, H. Hardelauf, and H. Vereecken: Numerical modeling of large scale transport of contaminant solutes using the global random walk algorithm. Monte Carlo Methods and Appl. 10(2) (2004), p. 155-179.