MSC 2000. Primary 65C35, 76R50; Secondary 82B41, 37B15.
Keywords. Overshooting, global random walk, groundwater contamination.

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(\mathbf{x},t)$$c(\mathbf{x},t)$c(x,t)c(\mathbf{x}, t)$c(\mathbf{x},t)$, the Fokker-Plank equation coincides with the advection dispersion parabolic equation

where $\mathbf{V}$$\mathbf{V}$V\mathbf{V}$\mathbf{V}$ is the velocity field and $\mathbf{D}$$\mathbf{D}$D\mathbf{D}$\mathbf{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
$$\begin{array}{}\text{(}\prime \text{)}& {\partial }_{t}c+\mathrm{\nabla }\left({\mathbf{V}}^{\ast }c\right)=\mathrm{\nabla }(\mathbf{D}\mathrm{\nabla }c)\end{array}$$$$\begin{array}{}\text{(}\prime \text{)}& {\partial }_{t}c+\mathrm{\nabla }\left({\mathbf{V}}^{\ast }c\right)=\mathrm{\nabla }(\mathbf{D}\mathrm{\nabla }c)\end{array}$${:('")"del_(t)c+grad(V^(**)c)=grad(Dgrad c):}\begin{equation*} \partial_{t} c+\nabla\left(\mathbf{V}^{*} c\right)=\nabla(\mathbf{D} \nabla c) \tag{$\prime$} \end{equation*}$$\begin{array}{}\text{(}\prime \text{)}& {\partial }_{t}c+\mathrm{\nabla }\left({\mathbf{V}}^{\ast }c\right)=\mathrm{\nabla }(\mathbf{D}\mathrm{\nabla }c)\end{array}$$
which is the local mass balance equation derived from the Fick's constitutive law (i.e., concentration flux $\sim -\mathbf{D}\mathrm{\nabla }c$$\sim -\mathbf{D}\mathrm{\nabla }c$∼-Dgrad c\sim-\mathbf{D} \nabla c$\sim -\mathbf{D}\mathrm{\nabla }c$ ) [3]. The equation (1') takes the form (1) and can be associated to Itô equation, if the velocity ${\mathbf{V}}^{\ast }$${\mathbf{V}}^{\ast }$V^(**)\mathbf{V}^{*}${\mathbf{V}}^{\ast }$ is supplemented with a drift vector with components given by the divergence of the rows of the diffusion tensor: $\mathbf{V}={\mathbf{V}}^{\ast }+\mathrm{\nabla }\mathbf{D}$$\mathbf{V}={\mathbf{V}}^{\ast }+\mathrm{\nabla }\mathbf{D}$V=V^(**)+gradD\mathbf{V}=\mathbf{V}^{*}+\nabla \mathbf{D}$\mathbf{V}={\mathbf{V}}^{\ast }+\mathrm{\nabla }\mathbf{D}$ [2].

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\delta t$$t=k\delta t$t=k delta tt=k \delta t$t=k\delta t$ and at a grid point $(x_1,x_2,x_3)=(i_1\delta x_1,i_2\delta x_2,i_3\delta x_3)$$\left(x_1,x_2,x_3\right)=\left(i_1\delta x_1,i_2\delta x_2,i_3\delta x_3\right)$(x_(1),x_(2),x_(3))=(i_(1)deltax_(1),i_(2)deltax_(2),i_(3)deltax_(3))\left(x_{1}, x_{2}, x_{3}\right)=\left(i_{1} \delta x_{1}, i_{2} \delta x_{2}, i_{3} \delta x_{3}\right)$(x_1,x_2,x_3)=(i_1\delta x_1,i_2\delta x_2,i_3\delta x_3)$ is obtained by counting the number $n(i_1,i_2,i_3,k)$$n\left(i_1,i_2,i_3,k\right)$n(i_(1),i_(2),i_(3),k)n\left(i_{1}, i_{2}, i_{3}, k\right)$n(i_1,i_2,i_3,k)$ of "tracked particles" in a reference volume,

where $N$$N$NN$N$ is the total number of particles and ${\mathrm{\Delta }}_l=2s_l\delta x_l,l=1,2,3$${\mathrm{\Delta }}_l=2s_l\delta x_l,l=1,2,3$Delta_(l)=2s_(l)deltax_(l),l=1,2,3\Delta_{l}=2 s_{l} \delta x_{l}, l=1,2,3${\mathrm{\Delta }}_l=2s_l\delta x_l,l=1,2,3$, are the lengths of symmetrical intervals centered at $x_l$$x_l$x_(l)x_{l}$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$$c$cc$c$ 5, 6, 9.

The GRW algorithm is equivalent to a superposition of many PT procedures 9. 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. For illustration we present the GRW algorithm in the case of constant $D$$D$DD$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)$$n(i,k)$n(i,k)n(i, k)$n(i,k)$ particles from ( $x_i,t_k$$x_i,t_k$x_(i),t_(k)x_{i}, t_{k}$x_i,t_k$ ) by

where $v_j$$v_j$v_(j)v_{j}$v_j$ are discrete displacements in a given velocity field and $d$$d$dd$d$ is the amplitude of unbiased diffusive jumps. The quantities $\delta n$$\delta n$delta n\delta n$\delta n$ are Bernoulli random variables and describe respectively, the number of particles which remain at the same grid site $j+v_j$$j+v_j$j+v_(j)j+v_{j}$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$$j+v_j$j+v_(j)j+v_{j}$j+v_j$ ). The mostly used version is the so called "reduced fluctuations" version of GRW, defined by

where $n=n(j,k)-\delta n(j,j+v_j,k),[n/2]$$n=n(j,k)-\delta n\left(j,j+v_j,k\right),[n/2]$n=n(j,k)-delta n(j,j+v_(j),k),[n//2]n=n(j, k)-\delta n\left(j, j+v_{j}, k\right),[n / 2]$n=n(j,k)-\delta n(j,j+v_j,k),[n/2]$ is the integer part of $n/2$$n/2$n//2n / 2$n/2$ and $\theta$$\theta$theta\theta$\theta$ is a random variable taking the values 0 and 1 with probability $1/2$$1/2$1//21 / 2$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 $\sim (2Dt{)}^{1/2})$$\left.\sim (2Dt{)}^{1/2}\right)${:∼(2Dt)^(1//2))\left.\sim(2 D t)^{1 / 2}\right)$\sim (2Dt{)}^{1/2})$.

The distribution of the particles at the next time $(k+1)\delta t$$(k+1)\delta t$(k+1)delta t(k+1) \delta t$(k+1)\delta t$ is given by

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

where $0\le r\le 1$$0\le r\le 1$0 <= r <= 10 \leq r \leq 1$0\le r\le 1$. The diffusion coefficient $D$$D$DD$D$ is expressed in terms of the grid steps by the relation

For $r\le 1$$r\le 1$r <= 1r \leq 1$r\le 1$, there is no numerical diffusion. Moreover, since the total number of particles $N$$N$NN$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 $\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})$; for large numbers of particles 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 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.

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$$d=1$d=1d=1$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$$v_j$v_(j)v_{j}$v_j$ and the unbiased diffusive jumps $d$$d$dd$d$ are independent, in BGRW the advective flux $v_j\overline{n}(j,k)$$v_j\overline{n}(j,k)$v_(j) bar(n)(j,k)v_{j} \bar{n}(j, k)$v_j\overline{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)\delta t$$(k+1)\delta t$(k+1)delta t(k+1) \delta t$(k+1)\delta 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 $\overline{n}(x,t)=\overline{n}(i\delta x,k\delta t)$$\overline{n}(x,t)=\overline{n}(i\delta x,k\delta t)$bar(n)(x,t)= bar(n)(i delta x,k delta t)\bar{n}(x, t)=\bar{n}(i \delta x, k \delta t)$\overline{n}(x,t)=\overline{n}(i\delta x,k\delta 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,\overline{n}(x,t)\simeq n(x,t)$$N,\overline{n}(x,t)\simeq n(x,t)$N, bar(n)(x,t)≃n(x,t)N, \bar{n}(x, t) \simeq n(x, t)$N,\overline{n}(x,t)\simeq 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 $\mathbf{r}=(r_1,r_2,r_3),\mathbf{v}=(v_1,v_2,v_3)$$\mathbf{r}=\left(r_1,r_2,r_3\right),\mathbf{v}=\left(v_1,v_2,v_3\right)$r=(r_(1),r_(2),r_(3)),v=(v_(1),v_(2),v_(3))\mathbf{r}=\left(r_{1}, r_{2}, r_{3}\right), \mathbf{v}=\left(v_{1}, v_{2}, v_{3}\right)$\mathbf{r}=(r_1,r_2,r_3),\mathbf{v}=(v_1,v_2,v_3)$, and accordingly modifying ( $3^{\prime }-5^{\prime }$$3^{\prime }-5^{\prime }$3^(')-5^(')3^{\prime}-5^{\prime}$3^{\prime }-5^{\prime }$ ) [7].

We consider an isotropic two-dimensional diffusion in groundwater ( $D_1=D_2=D=0.01{\mathrm{m}}^2/$$D_1=D_2=D=0.01{\mathrm{m}}^2/$D_(1)=D_(2)=D=0.01m^(2)//D_{1}=D_{2}= D=0.01 \mathrm{~m}^{2} /$D_1=D_2=D=0.01{\mathrm{m}}^2/$ day) in a velocity field $\mathbf{V}$$\mathbf{V}$V\mathbf{V}$\mathbf{V}$ with mean $U=1\mathrm{m}/$$U=1\mathrm{m}/$U=1m//U=1 \mathrm{~m} /$U=1\mathrm{m}/$ day, oriented along the $x_1$$x_1$x_(1)x_{1}$x_1$ axis and with a standard deviation $\sigma =0.2\mathrm{m}/$$\sigma =0.2\mathrm{m}/$sigma=0.2m//\sigma=0.2 \mathrm{~m} /$\sigma =0.2\mathrm{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}$$N={10}^{10}$N=10^(10)N=10^{10}$N={10}^{10}$, the total number of particles which ensures the self-averaging of the GRW simulations for this transport problem [6].

Since $D$$D$DD$D$ is constant we use in BGRW algorithm $\delta x_1=\delta x_2=\delta x$$\delta x_1=\delta x_2=\delta x$deltax_(1)=deltax_(2)=delta x\delta x_{1}=\delta x_{2}=\delta x$\delta x_1=\delta x_2=\delta x$ and $r_1=r_2=r=0.5$$r_1=r_2=r=0.5$r_(1)=r_(2)=r=0.5r_{1}=r_{2}= r=0.5$r_1=r_2=r=0.5$. Supposing that the maximum velocity can be as large as $V^{\text{max}}=U+5\sigma =2\mathrm{m}/$$V^{\text{max }}=U+5\sigma =2\mathrm{m}/$V^("max ")=U+5sigma=2m//V^{\text {max }}=U+5 \sigma= 2 \mathrm{~m} /$V^{\text{max }}=U+5\sigma =2\mathrm{m}/$ day, from the second relation (5') and the second condition (6) it follows that $\delta x⩽2D/V^{\text{max}}=0.01\mathrm{m}$$\delta x⩽2D/V^{\text{max }}=0.01\mathrm{m}$delta x <= 2D//V^("max ")=0.01m\delta x \leqslant 2 D / V^{\text {max }}=0.01 \mathrm{~m}$\delta x⩽2D/V^{\text{max }}=0.01\mathrm{m}$. Correspondingly, from the first relation in (5') one obtains $\delta t=0.0025$$\delta t=0.0025$delta t=0.0025\delta t=0.0025$\delta 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 4000 days, when an optimum choice of parameters, ( $\delta x=0.1\mathrm{m},\delta t=0.5$$\delta x=0.1\mathrm{m},\delta t=0.5$delta x=0.1m,delta t=0.5\delta x=0.1 \mathrm{~m}, \delta t=0.5$\delta x=0.1\mathrm{m},\delta t=0.5$ day, $r=0.25$$r=0.25$r=0.25r=0.25$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\mathrm{m}$$L_1=150\mathrm{m}$L_(1)=150mL_{1}=150 \mathrm{~m}$L_1=150\mathrm{m}$ and $L_2=20\mathrm{m}$$L_2=20\mathrm{m}$L_(2)=20mL_{2}=20 \mathrm{~m}$L_2=20\mathrm{m}$, so that during the total simulation time, $T=100$$T=100$T=100T=100$T=100$ days, no particle reaches the boundary. The concentration $c$$c$cc$c$ was computed by (2) for ${\mathrm{\Delta }}_1={\mathrm{\Delta }}_2=\mathrm{\Delta }=1\mathrm{m}$${\mathrm{\Delta }}_1={\mathrm{\Delta }}_2=\mathrm{\Delta }=1\mathrm{m}$Delta_(1)=Delta_(2)=Delta=1m\Delta_{1}=\Delta_{2}=\Delta=1 \mathrm{~m}${\mathrm{\Delta }}_1={\mathrm{\Delta }}_2=\mathrm{\Delta }=1\mathrm{m}$. The 1st and the 2-nd moments of the concentration $c$$c$cc$c$ are defined by

where $l=1,2$$l=1,2$l=1,2l=1,2$l=1,2$ and the integrals are computed over the support of $c$$c$cc$c$. For comparisons we used the derivative of the 1 -st moments $V_l^{cm}(t)=d{\mu }_l(t)/dt$$V_l^{cm}(t)=d{\mu }_l(t)/dt$V_(l)^(cm)(t)=dmu_(l)(t)//dtV_{l}^{c m}(t)=d \mu_{l}(t) / d t$V_l^{cm}(t)=d{\mu }_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_{ll}^{eff}(t)={\mu }_{ll}(t)/(2t)$$D_{ll}^{eff}(t)={\mu }_{ll}(t)/(2t)$D_(ll)^(eff)(t)=mu_(ll)(t)//(2t)D_{l l}^{e f f}(t)=\mu_{l l}(t) /(2 t)$D_{ll}^{eff}(t)={\mu }_{ll}(t)/(2t)$ [8]. The cross-section concentration was computed by

where $L_2$$L_2$L_(2)L_{2}$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={\mu }_1(t)$$x_1={\mu }_1(t)$x_(1)=mu_(1)(t)x_{1}= \mu_{1}(t)$x_1={\mu }_1(t)$.

The evaluation of the center of mass velocity was done using the absolute errors

where $\mathrm{\Phi }$$\mathrm{\Phi }$Phi\Phi$\mathrm{\Phi }$ stands, respectively, for the expectation $E\left(V_l^{cm}\right)$$E\left(V_l^{cm}\right)$E(V_(l)^(cm))E\left(V_{l}^{c m}\right)$E\left(V_l^{cm}\right)$, computed by arithmetic means over the 256 velocity realizations, and for the corresponding standard deviations $SD\left(V_l^{cm}\right)={\{E\left[{\left(V_l^{cm}\right)}^2\right]-{\left[E\left(V_l^{cm}\right)\right]}^2\}}^{1/2}$$SD\left(V_l^{cm}\right)={\left\{E\left[{\left(V_l^{cm}\right)}^2\right]-{\left[E\left(V_l^{cm}\right)\right]}^2\right\}}^{1/2}$SD(V_(l)^(cm))={E[(V_(l)^(cm))^(2)]-[E(V_(l)^(cm))]^(2)}^(1//2)S D\left(V_{l}^{c m}\right)=\left\{E\left[\left(V_{l}^{c m}\right)^{2}\right]-\left[E\left(V_{l}^{c m}\right)\right]^{2}\right\}^{1 / 2}$SD\left(V_l^{cm}\right)={\{E\left[{\left(V_l^{cm}\right)}^2\right]-{\left[E\left(V_l^{cm}\right)\right]}^2\}}^{1/2}$. The results for $l=1$$l=1$l=1l=1$l=1$ and $l=2$$l=2$l=2l=2$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$$T_1=1$T_(1)=1T_{1}=1$T_1=1$ day.
The evaluation for the effective coefficients $D_{ll}^{eff}$$D_{ll}^{eff}$D_(ll)^(eff)D_{l l}^{e f f}$D_{ll}^{eff}$ and cross-section concentration $C$$C$CC$C$ were achieved by using the percentage relative errors

where $\mathrm{\Phi }$$\mathrm{\Phi }$Phi\Phi$\mathrm{\Phi }$ stands, again, for the corresponding expectations, $E(\cdot )$$E(\cdot )$E(*)E(\cdot)$E(\cdot )$, and standard deviations, $SD(\cdot )$$SD(\cdot )$SD(*)S D(\cdot)$SD(\cdot )$. 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 $\Vert \epsilon (\mathrm{\Phi })\Vert$$\Vert \epsilon (\mathrm{\Phi })\Vert$||epsi(Phi)||\|\varepsilon(\Phi)\|$\Vert \epsilon (\mathrm{\Phi })\Vert$.

A comparison between figure 1 and the theoretical expectation values $E\left(V_1^{cm}\right)=1\mathrm{m}/$$E\left(V_1^{cm}\right)=1\mathrm{m}/$E(V_(1)^(cm))=1m//E\left(V_{1}^{c m}\right)= 1 \mathrm{~m} /$E\left(V_1^{cm}\right)=1\mathrm{m}/$ day and $E\left(V_2^{cm}\right)=0\mathrm{m}/$$E\left(V_2^{cm}\right)=0\mathrm{m}/$E(V_(2)^(cm))=0m//E\left(V_{2}^{c m}\right)=0 \mathrm{~m} /$E\left(V_2^{cm}\right)=0\mathrm{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\mathrm{\%}$$5\mathrm{\%}$5%5 \%$5\mathrm{\%}$. 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$$D$DD$D$. After tens of days the transverse effective coefficient already reaches the theoretical asymptotic value $D_{22}^{eff}\sim D$$D_{22}^{eff}\sim D$D_(22)^(eff)∼DD_{22}^{e f f} \sim D$D_{22}^{eff}\sim 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$$D$DD$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.

Acknowledgement. The research reported in this paper is a contribution to Deutsche Forschungsgemeinschaft grant SU 415/1-1 and to "Interdisciplinary Programme for the Prevention of the Major Risk Phenomena at National Level" of the Romanian Academy.

[1] Doob, J. L., Stochastic Processes, John Wiley & Sons, London, 1953.
[2] Gardiner, C. W. (1985), Handbook of Stochastic Methods (for Physics, Chemistry and Natural Science), Springer, New York.
[3] Kitanidis, P. K., Particle-tracking equations for solution of the advection-dispersion equation with variable coefficients, Water Resour. Res., 30(11), 3225-3227, 1994.
[4] Kloeden, P. E. and E, Platten, Numerical solutions of stochastic differential equations, Springer, Berlin, 1995.
[5] SUN, Ne-Z., Mathematical Modeling in Groundwater Pollution, Springer, New York, 1996.
[6] 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), 153-177, 2004.
[7] Suciu, N., C. Vamoş, P. Knabner and U. Rüde, Biased global random walk, a cellular automaton for diffusion, in Simulationstechnique, $18{}^{\text{th}}$$18{}^{\text{th }}$18^("th ")18{ }^{\text {th }}$18{}^{\text{th }}$ Symposium in Erlangen, September 2005, F. Hülsemann, M. Kowarschik and U. Rüde (eds.), pp. 562-567, SCS Publishing House e. V., Erlangen, 2005.
[8] Suciu N., C. Vamoş, J. Vanderborght, H. Hardelauf and H. Vereecken, Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi:10.1029/2005WR004546, 2006.
[9] Vamoş, C., N. Suciu and H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186(2), 527-544, 2003.

Received by the editors: February 26, 2006.