Numerical Modeling of Large Scale Transport of Contaminant Solutes Using the Global Random Walk Algorithm

Nicolae Suciu
(original), global random walk, Monte Carlo simulations, paper

Abstract

The methods which track particles to simulate concentrations are successfully used in problems for large scale transport in groundwater, mainly when the aquifer properties are spatially heterogeneous and the process is advection dominated. These methods, sometimes called “analogical Monte Carlo methods”, are not concerned with the numerical diffusion occurring in finite difference/element schemes. The limitations of classical random walk methods are due to large computation time and memory necessary to achieve statistically reliable results and accurate concentration fields. To overcome these computational limitations a “global random walk” (GRW) algorithm was developed. Unlike in the usual approach where the trajectory of each particle is simulated and stored, in GRW all the particles from a given grid node are scattered, at a given time, using a single numerical procedure, based on Bernoulli distribution, partial-deterministic, or deterministic rules. Because no restrictions are imposed for the maximum number of particles to be used in a simulation, the Monte Carlo repetitions are no longer necessary to achieve the convergence. It was proved that for simple diffusion problems GRW converges to the finite difference scheme and that for large scale transport problems in groundwater, GRW produces stable and statistically reliable simulations. A 2-dimensional transport problem was modeled by simulating local diffusion processes in realizations of a random velocity field. The behavior over 5000 days of the effective diffusion coefficients, concentration field and concentration fluctuations were investigated using 2560 realizations of the velocity field and 1010 particles in every realization. The results confirm the order of magnitude of the effective diffusion coefficients predicted by stochastic theory but the time needed to reach the asymptotic regime was found to be thousands times larger. It is also underlined that the concentration fluctuations and the dilution of contaminant solute depend essentially on local diffusion and boundary conditions.

 

Authors

N. Suciu
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

C. Vamoş
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

H. Hardelauf
Research Center Julich (Germany), ICG-IV: Institut of Agrosphere

J. Vanderborght
Research Center Julich (Germany), ICG-IV: Institut of Agrosphere

H. Vereecken
Research Center Julich (Germany), ICG-IV: Institut of Agrosphere

Keywords

?

Cite this paper as

N. Suciu, C. Vamoş, J. Vanderborght, H. Hardelauf, H. Vereecken (2004), Numerical modeling of large scale transport of contaminant solutes using the global random walk algorithm, Monte Carlo Methods and Applications, 10(2), 153-177, doi: 10.1163/156939604777303235

References

see the expansion block below.

PDF

2004 Vamos Numerical Modeling of Large Scale Transport of Contminant Solutes Using the Global Random Walk Algorithm

About this paper

Journal

Monte Carlo Methods and Applications

Publisher Name
DOI

10.1163/156939604777303235

Print ISSN

0929-9629

Online ISSN

1569-3961

Google Scholar Profile

google scholar profile

[1] Ames, W. F., Numerical Methods for Partial Differential Equations, 2nd ed., Academic Press, New York, 1977.

[2] Bellin, A., P. Salandin and A. Rinaldo, Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations, Water Resour. Res. 28(9), 2211-2227, 1992.
CrossRef (DOI)

[3] Chorin, A. J., Vortex sheet approximation of boundary layers, J. Comput. Phys. 27, 428-442 (1978).
CrossRef (DOI)

[4] Dagan, G., Flow and Transport in Porous Formations, Springer, New York, 1989.
CrossRef (DOI)

[5] Dagan, G. and A. Fiori, The influence of pore-scale dispersion on concentration statistical moments in transport through heterogeneous aquifers, Water Resour. Res., 33, 1595-1650, 1997.
CrossRef (DOI)

[6] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach, Temporal behavior of a solute cloud in a heterogeneous porous medium 3. Numerical simulations, Water Resour. Res., 38, 23:1-12, 2002.
CrossRef (DOI)

[7] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach, Numerical studies of the transport behavior of a passive solute in a two-dimensional incompressible random flow field, Phys. Rev. E 67, 046306:1-10, 2003.
CrossRef (DOI)

[8] Fiori, A., Finite Peclet extension of Dagan’s solutions to transport in anisotropic formations, Water Resour. Res., 32, 193-198, 1996.
CrossRef (DOI)

[9] Kapoor, C. and L. W. Gelhar, Transport in three-dimensionallity heterogeneous aquifers 1. Dynamics of concentration fluctuations, Water Resour. Res., 30, 1775-1789, 1994.
CrossRef (DOI)

[10] Kapoor, C. and P. K. Kitanidis, Concentration fluctuation and dilution in aquifers, Water Resour. Res., 34, 1181-1193, 1998.
CrossRef (DOI)

[11] Lundgren, T. S., Turbulent pair dispersion and scalar diffusion, J. Fluid Mech., 111, 27-57, 1981.
CrossRef (DOI)

[12] Pannone, M. and P. K. Kitanidis, Large-time behavior of concentration variance and dilution in heterogeneous media, Water Resour. Res., 35, 623-634, 1999.
CrossRef (DOI)

[13] Roth, K. and K. Hammel, Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with steady state flow, Water Resour. Res., 32, 1653-1663, 1996.
CrossRef (DOI)

[14] Sabelfeld, K. K. and G. A. Mikhailov, On the numerical simulation of the particle’s diffusion in random velocity fields, Izv. AN SSSR, Physics of Atmosphere and Ocean, 16 (3), 229-235 (in Russian), 1980.

[15] Sabelfeld, K. K., Monte Carlo methods in boundary value problems, Springer, New York – Heidelberg – Berlin, 1991.

[16] Schwarze, H., U. Jaekel and H. Vereecken, Estimation of macrodispersivity by different approximation methods for flow and transpot in randomly heterogeneous media, Transport in Porous Media 43, 265-287, 2001.
CrossRef (DOI)

[17] Suciu, N., C. Vamos, H. Vereecken, and J.Vanderborght, pp. 111-140 in G. Marinovschi and I. Stelian (Eds.) Proceedings of the first Workshop on Mathematical Modeling of Environmental Problems, Rom. Acad. Publishing House, Bucharest 2002.
CrossRef (DOI)

[18] Suciu N., C. Vamos, H. Vereecken, J. Vanderborght, and H. Hardelauf, pp. 215-221 in T. Maghiar, A. Georgescu, M. Balaj, I. Dzitac, I. Mang (Eds.), Proc. of the 11th Conf. of Appl. and Ind. Math., Vol. 1, Oradea Univ. Publishing House, Oradea, 2003.

[19] Suciu, N. and C. Vamos, Effective diffusion in heterogeneous media, Interner Bericht ICG-IV 00303, Forschungszentrum Julich, 2003. Sun, Ne-Z., Mathematical Modeling in Groundwater Pollution, Springer, New York, 1996.

[20] Vamos, C., N. Suciu, H. Vereecken, J. Vanderborht and O. Nitzsche, Path decomposition of discrete effective diffusion coefficient, Internal Report ICG-IV.00501, Forschungszentrum Julich, 2001.

[21] Vamo¸s, C., N. Suciu, and H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186, 527-544, 2003.
CrossRef (DOI)

[22] van Kampen, N. G., Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1981.

Monte Carlo Methods and Appl., Vol. 10, No. 2, pp. 153 – 177 (2004) c VSP 2004 Numerical Modeling of Large Scale Transport of Contminant Solutes Using the Global Random Walk Algorithm ∗ N. Suciu and C. Vamo¸s Romanian Academy, ”T. Popoviciu” Institute of Numerical Analysis H. Hardelauf, J. Vanderborght, and H.Vereecken Research Center J¨ ulich (Germany), ICG-IV: Institut of Agrosphere Abstract — The methods which track particles to simulate concentrations are successfully used in problems for large scale transport in groundwater, mainly when the aquifer properties are spatially heterogeneous and the process is advection dominated. These methods, sometimes called “analogical Monte Carlo methods”, are not concerned with the numerical diffusion occur- ring in finite difference/element schemes. The limitations of classical random walk methods are due to large computation time and memory necessary to achieve statistically reliable results and accurate concentration fields. To overcome these computational limitations a “global random walk” (GRW) algorithm was developed. Unlike in the usual approach where the trajectory of each particle is simulated and stored, in GRW all the particles from a given grid node are scattered, at a given time, using a single numerical procedure, based on Bernoulli distribution, partial-deterministic, or deterministic rules. Because no restrictions are imposed for the maxi- mum number of particles to be used in a simulation, the Monte Carlo repetitions are no longer necessary to achieve the convergence. It was proved that for simple diffusion problems GRW converges to the finite difference scheme and that for large scale transport problems in ground- water, GRW produces stable and statistically reliable simulations. A 2-dimensional transport problem was modeled by simulating local diffusion processes in realizations of a random veloc- ity field. The behavior over 5000 days of the effective diffusion coefficients, concentration field and concentration fluctuations were investigated using 2560 realizations of the velocity field and 10 10 particles in every realization. The results confirm the order of magnitude of the effective diffusion coefficients predicted by stochastic theory but the time needed to reach the asymptotic regime was found to be thousands times larger. It is also underlined that the concentration fluctuations and the dilution of contaminant solute depend essentially on local diffusion and boundary conditions. 1 Introduction The solutions of the parabolic equations can be numerically obtained using random walkers to describe the contribution of the second order space partial derivatives [Chorin, 1978]. For simple diffusion processes, the concentration is obtained by simulating the random walk trajectories for a collection of particles and by counting the number of * Paper presented at the IVth IMACS Seminar on Monte Carlo Methods MCM-2003, 15-19 September 2003, Berlin.
154 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken particles in a reference volume. This algorithm, equivalent to the finite difference scheme (FD), is a analogical Monte Carlo method, where the accuracy of solutions is improved by increasing the number of particles [Ames, 1977]. The “particle tracking method” (PT) which generalizes the algorithm for advection-diffusion processes is suitable to simulate the transport in heterogeneous velocity fields, as in problems related to monitoring the fate of contaminant solutes in groundwater [Roth and Hammel, 1996]. With respect to FD and finite element methods, PT has the advantage to be structurally stable and without numerical diffusion but in many practical problems its applicability is limited by the computing resources necessary to simulate the trajectories of collections of particles large enough to ensure the required accuracy [Sun, 1996]. In [Vamo¸s et al., 2003] it was shown that the “global random walk” (GRW) is an alternative tool to describe the transport in heterogeneous velocity fields. This new par- ticles method scatters all particles lying at a given grid point simultaneously. Unlike in usual random walk algorithms where the trajectories of the particles are simulated indi- vidually and stored, the GRW provides the whole concentration field at each time step and it is not concerned with limitations related to the number of tracked particles. In the present work the GRW is used to obtain statistically reliable results in simulation of the large time behavior of a typical transport process in groundwater. Following the usual stochastic approach of the problem, the large scale heterogeneities of the natural porous medium are described by a random velocity space field and the effects of erratic structure at local scales are modeled as a diffusive process, called “pore scale dispersion”. Collec- tion of particles are used to obtain the concentration and other quantities of interest in given realizations of the field and averages over realizations are taken to obtain statistical estimations. The accuracy of results depends on the choice of space and time steps of the grid used in the GRW algorithm, the number of particles used to simulate the transport in a given realization and the number of realizations used to obtain the statistical estimations. The grid steps give the lower bound of the numerical errors in simulation of the diffusion process and, mainly, of the advective transport. Even stable, the GRW, as well as PT, yields “overshooting” errors when particles jump over grid points with different assigned velocities. These errors distort the advection and cannot be completely eliminated. To reduce the overshooting, refined grids are necessary implying high computational costs. A trade off between accuracy and the aim to obtain statistically reliable large scale simula- tions, describing the transport during 5000 days, was achieved by choosing an acceptable overshooting which keeps the numerical results close to those given by first order approx- imations from stochastic theory. The stochastic convergence of the simulations of the transport in given realization of the velocity field, towards the threshold corresponding to the fixed diffusion and overshooting errors, was obtained using N = 10 10 particles. Accu- rate statistical estimations were obtained by simulating the transport in 2560 realizations of the velocity field. 2 The global random walk algorithm The GRW algorithm was proposed as a generalization of the random walk and PT methods, aiming to increase the speed of computations and to increase the accuracy [Vamo¸s et al., 2003]. The one-dimensional GRW algorithm describes the scattering of the
Numerical Modeling of Large Scale Transport of Contminant Solutes 155 n(i, k) particles from (x i ,t k ) by n(j, k)= δn(j, j + v j ,k)+ δn(j + v j − d, j, k)+ δn(j + v j + d, j, k), where v j are discrete displacements in a given velocity field, d describes the diffusive jumps, and δn(j + v j ± d, j, k) are Bernoulli random variables. The distribution of the particles at the next time (k + 1)δt is given by n(i, k + 1) = j δn(i, j, k). For two and three-dimensional cases, the same procedure is repeated for all space di- rections. The average number of particles undergoing diffusive jumps and the average number of particles particles remaining at the same node after the displacement v j are given by the use of the parameter r, r 1, trough the relations δn(j + v j ± d, j, k)= 1 2 r n(j, k) and δn(j, j + v j ,k) = (1 − r) n(j, k). The diffusion coefficient D is related to the grid steps through D = r (dδx) 2 2δt . (1) Because the total number of particles N contained in the grid is conserved the GRW algorithm is stable and (1), with condition r 1, ensures that there is no numerical diffusion. The GRW can be implemented in several variants. The fast algorithm, which we call simply GRW, approximates the Bernoulli distribution for N> 10 20 with that correspond- ing to 10 20 particles without producing significant errors. GRW0 is the algorithm where no approximations are made. GRWG uses the Gauss distribution (the limit of Bernoulli distribution for N −→ ∞). GRWD, the determinist variant obtained when n are real numbers and δn(j + v j ± d, j, k)= 1 2 rn(j, k), is shown to be equivalent to the FD scheme for constant velocity. The “reduced fluctuations algorithm” GRWR is defined by δn(j + v j − d, j, k)= n/2 if n is even [n/2] + θ if n is odd , where n = n(j, k) − δn(j, j + v j ,k), [n/2] is the integer part of n/2 and θ is a variable taking the values 0 and 1 with probability 1/2. GRWR is suitable for large scale problems because the diffusion front do not extend beyond the limit concentration defined by one particle at a grid point, keeping a physical significant shape (unlike in FD where a pure diffusion front has a cubic shape of side ∼ √ Dt), and it also requires a minimum number of calls of the random number generator. The main difference with respect to the PT algorithm is that while the output of PT is a single trajectory of a particle undergoing advection and random walk, the GRW algorithm yields paths containing the trajectories of all the N particles used in simulation. Because at every time step the whole probability density (i.e. the normalized concentration inferred from the number of particles at grid sites) is available, the GRW, as a mathematical object is an approximation of a discrete diffusion process. An illustration for the 2-dimensional case is presented in Fig. 1. For a given velocity field u, from all points (i, j ) containing particles develop paths P (s) α (i, j ; u) of length α = t/δt so that at the intermediate time step s the n(i, j ) particles are spread in 4 s different grid points [Vamo¸s et al. 2001].
156 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken Fig. 1. The path P (2) α (i, j ; u), after 2 time steps, of particles starting at (i, j ) for a realization u of the velocity field. In the following some results regarding the properties and performances of GRW are quoted from [Vamo¸s et al., 2003]. A comparison with the PT code ”ParTrace”, on Cray T3E, for a small scale problem (constant diffusion coefficient, no velocity, diffusion of N particle starting at the center of a cubic grid the side of which was 20 space steps) show that while in GRW algorithm there is practical no limitation, N> 10 9 particles becomes prohibitive for PT (Fig.2). Fig. 2. Comparison between the computing times required by PT and GRW for increasing N . (c 2003 Elsevier Science B. V. All rights reserved) The results obtained with GRW algorithm for a one-dimensional Wiener process, de- scribed by ∂ t c = D∂ 2 x c, lim t−→0+ c = Nδ (x), where δ (x) is the Dirac function, were
Numerical Modeling of Large Scale Transport of Contminant Solutes 157 compared with the Gaussian solution trough the norm ‖c − c Gauss ‖. The results from Fig. 3 show that the numerical solution converges to the Gaussian as O(δx 2 )+O(1/ √ N ), i.e. for large numbers of particles, N , the convergence order is O(δx 2 ), the same as for FD. Fig. 3. The convergence as function of N for δx =0.1 (a) and δx =0.01 (b).(c 2003 Elsevier Science B. V. All rights reserved) The dependence of the norm on r was studied with the GRWD algorithm, equivalent to FD scheme (Fig. 4). It was found that the best precision is obtained for r ≃ 0.3. Fig. 4. Dependence on r. (c 2003 Elsevier Science B. V. All rights reserved) When N = 1 (PT procedure) the value ‖c − c Gauss ‖ ∼ =0.01, corresponding to FD scheme, was reached for S ∼ = 10 6 repetitions of the simulation and the computation time was of about 17 minutes. The same result is obtained in a single simulation, S = 1, for N = 10 6 , in 0.5 seconds (Fig. 5).
158 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken Fig. 5. Convergence as function of S .( c 2003 Elsevier Science B. V. All rights reserved) The conclusion is that for large enough N the Monte Carlo repetitions are not necessary to obtain the required precision in GRW algorithm. 3 Transport in heterogeneous aquifers The transport of non-reactive solutes in heterogeneous aquifers depends on the interplay of the following factors: 1) the local dispersion with constant coefficient D, described as a Wiener process w, 2) the initial position X 0 of solute molecules, and 3) the large scale space variability of the advection velocity field V. The width of diffusion fronts is described by the mean square displacements σ 2 ij (t)= 〈(X i −〈X i 〉)(X j −〈X j 〉)〉 , where 〈···〉 = 〈···〉 w,X 0 denotes the average over the Wiener process and initial positions. When the effective coefficients D eff i,j = lim t→∞ σ 2 ij (t) / (2t) can be defined, then an ”up- scaled” description similar to the Wiener process is possible: σ 2 ij ∼ 2D eff ij t, for t →∞. The structure of the effective coefficients is given by D eff ij = D + D adv ij − D cm ij + M ij , (2) where D adv ij = lim t→∞ t 0 〈V i V j 〉 (t, t ′ )dt ′ and D cm ij = 1 2 lim t→∞ 〈V i 〉 (t) t 0 〈V j 〉 (t ′ )dt ′ + 〈V j 〉 (t) t 0 〈V i 〉 (t ′ )dt ′ are the ”advective” and ”center of mass” components of the effective coefficients and M ij = 1 2 lim t→∞ [〈X 0 i V j + X 0 j V i 〉− (〈X 0 i 〉〈V j 〉 + 〈X 0 j 〉〈V i 〉)] (t) is a ”momentum” term which vanishes for a singular initial concentration [Vamo¸s et al., 2001; Suciu and Vamo¸s, 2003; Suciu et al., 2003].
Numerical Modeling of Large Scale Transport of Contminant Solutes 159 The stochastic approach considers realizations of a random field V = U + u, with mean U and fluctuation u, and defines the effective coefficients as averages over the realizations of the velocity field of the effective coefficients in single realizations, D eff ij = D eff ij V = D + D adv ij V − D cm ij V . The ergodic hypothesis [Dagan, 1989] assumes that the effective coefficients in a given realization equal their estimations, and D cm ii = 0, i.e. D eff ij = D + D adv ij . The 1-st order approximation [Dagan, 1989; Fiori, 1996] yields u = O(σ y ) and, for mean flow U =(U, 0, 0), D eff 11 = D + Uλ y σ 2 y ,D eff 22 = D eff 33 = D, (3) where σ 2 y is the variance of the logarithm of the hydraulic conductivity y and λ y is the corresponding correlation length. 4 A typical large scale transport problem Two-dimensional velocity fields with constant mean U =(U, 0), U =1 m/day , ex- ponential correlated normal y , with correlation length λ y =1 m and variance σ 2 y =0.1 and σ 2 y = 1, were generated using 640 Fourier modes by the Kraichnan procedure (as in [Schwarze et al., 2001; Dentz et al. 2003]) on a Cray T3E computer. Because the results of the stochastic theory in first order approximation are obtained for unbounded domains, the grid was chosen to be larger than the maximum extension of the plume. All the N particles were initially located at the origin of the grid and the computations were conducted for dimensionless times t/Uλ y corresponding to 5000 correlation lengths, using the reduced fluctuations algorithm GRWR. The transport depends on both heterogeneity of the field, described by σ 2 y , and the local dispersion coefficient D. For N = 10 10 , as suggested by Fig. 3, and choosing the space and time steps δx 1 = δx 2 =0.25 m and δt =1 day , several combinations of σ 2 y and D were investigated.
160 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken 0 10 20 30 40 50 60 70 80 90 100 300 400 500 600 700 800 900 1000 1100 1200 1300 x 2 (m) x 1 (m) σ 2 y =1.0 and D=0.001 σ 2 y =0.1 and D=0.001 σ 2 y =1.0 and D=0.01 σ 2 y =0.1 and D=0.01 Fig. 6. Diffusion fronts at t = 1000days. Comparing the results presented in Fig. 6, the transport parameters σ 2 y =0.1 and D =0.01 m 2 /day were chosen, because they lead to symmetric plumes, as in first order approximation. In this case, the theoretical values of the effective diffusion coefficients given by (3) are D eff 11 =0.11 m 2 /day and D eff 22 =0.01 m 2 /day . 5 Computational parameters In the following the number of particles is fixed at N = 10 10 and because D is constant and isotropic, one chooses δx 1 = δx 2 = δx. To avoid the overshooting due to the diffusive jumps one takes d = 1. The parameters δx, δt, and d are related to r and D through the relation (1). The structure of the velocity field is correctly reproduced when the space step and all the displacements in a time step are smaller than the correlation length λ y . The following constraint are imposed: δx < λ y (4a) dδx < λ y (4b) Uδt λ y (4c) The space step δx =0.25 m fulfils the condition (4a). Previous published works also recommend to choose at least four velocity values per correlation length to obtain an accurate description of the velocity field [Bellin et al., 1992]. The condition (4b) being met for d = 1, it remains to choose the time step δt satisfying (4c). The test for different values Uδt/δx in Fig. 7 shows that the minimum time step which gives effective coefficients close to the theoretical values is δt =1 day , corresponding to Uδt/δx = 4. From (1) one obtains a value r =0.32, which is closed to the optimum in Fig 4.
Numerical Modeling of Large Scale Transport of Contminant Solutes 161 0 1 2 3 4 5 6 7 1 2 3 4 5 6 7 8 9 10 m 2 /day D adv 11 D eff 11 D cm 11 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 1 2 3 4 5 6 7 8 9 10 m 2 /day Uδt/δx D adv 22 D eff 22 D cm 22 Fig. 7. Asymptotic coefficients, defined as time averages between t = 4900 and t = 5000days, for δx =0.25 m. Larger values, Uδt/δx > 4, infringes the condition (4c). A choice Uδt/δx < 4 yields unrealistic effective coefficients and do not describe the diffusion fronts (Fig. 8). Thus for fixed δx =0.25 m, the only possible choice is δt =1 day .
162 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken Fig. 8. Effective coefficients and diffusion fronts for δx =0.25m andUδt/δx = 1. For fixed Uδt/δx = 4, the choice δx > 0.25 (i.e. less than 4 velocity values per correlation length), do not reproduce the structure of the velocity field (Uδt > λ) and the effective coefficients become larger than the theoretical values (Fig. 9).
Numerical Modeling of Large Scale Transport of Contminant Solutes 163 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 m 2 /day D adv 11 D eff 11 D cm 11 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 m 2 /day δx D adv 22 D eff 22 D cm 22 Fig. 9. Effective coefficients, defined as time averages between t = 4900 and t = 5000 days, for Uδt/δx = 4. Because smaller space and time steps become prohibitive as computing resources, the choice δx =0.25 m, δt =1 day , and d = 1 is a trade off on accuracy which permits the simulation of the asymptotic effective diffusion coefficients and the comparison with theoretical results. To estimate the numerical errors produced with these parameters, the effect of over- shooting should be analyzed. Even if there is no recognized non-dimensional parameter for overshooting, a meaningful description can be obtained using the ratio between the mean square differences of the velocity over a mean displacement in a time step and the square mean velocity. For statistically homogeneous velocity fields, this parameter, hereinafter called Ov, can be computed as it follows: Ov = [V(x) − V(x + Uδt)] 2 U 2 = 2σ 2 V U 2 [1 − ρ(Uδt)] , (5)
164 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken where σ 2 V is the velocity variance and ρ V (Uδt)=[〈V(x)V(x + Uδt)〉− U 2 ]/σ 2 V the cor- responding correlation function evaluated at Uδt. For the exponential correlation of the logarithm of the hydraulic conductivity considered here, and in the 1-st order approxima- tion of the flow, one obtains σ 2 V = U 2 σ 2 y and ρ V (|x|) = exp(−|x| /λ y ). Considering (4c) and the first order of the Taylor expansion in |x| /λ y , one obtains Ov =2 σ 2 y Uδt λ y . (6) Thus the overshooting parameter introduced here depends on mean velocity, U , veloc- ity variance σ 2 V , and correlation ρ V . The shape of the correlation function is impor- tant for the value of the overshooting coefficient. For instance, when the logarithm of the hydraulic conductivity has a Gaussian correlation, then in 1-st approximation ρ V (|x|) = exp(−|x| 2 /λ y ), and instead of (6) one obtains Ov =4σ 2 y Uδt/λ y . It is also to be noticed that for large d the diffusive displacements becomes important and the maximum overshooting is described when ρ(Uδt + dδx) is considered in (5). For the previously chosen parameters, δx =0.25 m, δt =1 day , and d = 1, the overshooting coefficient is Ov =0.2. One remarks that even if the overshooting errors are reduced when reducing δt, for fixed δx, the quality of numerical results do not improve, moreover it is worse (see Fig. 7 and 8). This can happen, first, because the parameter r also modifies, according to (1), and departure from the optimum r =0.3 increases the numerical errors related to diffusion, and, second, due to an inaccurate representation of the variability of the velocity for small values Uδt/δx. The numerical errors can be reduced by reducing both δt and δx, choosing d so that (1) gives a value r close to 0.3, and keeping a small Ov and a large enough Uδt/δx. Results of preliminary tests are presented in Fig. 10, where computations of the effective coefficients, for a fixed realization of the velocity field, are compared for the set of chosen parameters (δx =0.25 m, δt =1 day , Uδt/δx = 4, d = 1) and for three sets of parameters corresponding to δx =0.1 m. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 D eff 11 (m 2 /day) Ut/λ y δx=0.25m δt=1.0day Uδt/δx=4 r=0.32 Ov=0.20 δx=0.10m δt=0.4day Uδt/δx=4 r=0.80 Ov=0.08 δx=0.10m δt=0.5day Uδt/δx=5 r=1.00 Ov=0.10 δx=0.10m δt=0.6day Uδt/δx=6 r=0.30 Ov=0.12 Fig. 10a. Comparison of longitudinal effective coefficients for different δx and δt.
Numerical Modeling of Large Scale Transport of Contminant Solutes 165 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 1 10 100 1000 D eff 22 (m 2 /day) Ut/λ y δx=0.25m δt=1.0day Uδt/δx=4 r=0.32 Ov=0.20 δx=0.10m δt=0.4day Uδt/δx=4 r=0.80 Ov=0.08 δx=0.10m δt=0.5day Uδt/δx=5 r=1.00 Ov=0.10 δx=0.10m δt=0.6day Uδt/δx=6 r=0.30 Ov=0.12 Fig. 10b. Comparison of transverse effective coefficients for different δx and δt. One remarks that, for δx =0.1 m and almost the same overshooting Ov, the numerical results are closer to the theoretical predictions (i.e. D cm 11 is smaller and D eff 22 closer to D =0.01 m/day ) in the last two cases in Fig. 10a and 10b, corresponding to the values Uδt/δx > 4. Then the accuracy is less sensitive to the parameter r and mainly depends on the discretization of the velocity field. It was also checked that correcting the overshooting by replacing the actual velocity with its mean over a time step (as proposed in [Roth and Hammel, 1996]), do not significantly change the results for δx =0.1 m and Uδt/δx =5 [Suciu et al., 2002]. Thus the difference between the first and the last two cases in Fig. 10 could be considered as an estimation of the order of magnitude for the numerical errors due to the discretization used in the following computations. 6 Stochastic convergence for given velocity realiza- tion In the case of the Wiener process, presented in Fig. 3, the stochastic convergence of GRW simulations towards the Gaussian solution was quantified by a norm (defined by means of the global variance with respect to the exact solution) and compared to the convergence of the FD scheme, quantified by the same norm. It was found that, for the same space-time discretization, the norm of the GRW solution cannot be smaller than the norm of the FD scheme, regardless how large N is. Thus, the accuracy of the FD solution defines a threshold for the GRW, corresponding to a given discretization. For the transport in realizations of a random velocity field, considered here, there is neither analytical solution nor other numerical solution available. Thus the stochastic convergence cannot be investigated in the same way. For the fixed parameters δx =0.25 m, δt =1 day , Uδt/δx = 4 and d = 1, the number of particles N was progressively increased until the effective coefficients reach stable values. Once the effective coefficients remain unchanged when N is increased, the threshold corresponding to the fixed diffusion and overshooting errors is reached. It was
166 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken found that for N 10 10 the values of all the coefficients (2) are identical in the limit of double precision (Fig. 11). Thus, for the large scale transport problem considered, the convergence of simulations for a given realization of the velocity field requires numbers of particles of the order N = 10 10 . 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 1 100000 1e+10 1e+15 1e+20 1e+25 m 2 /day D adv 11 D eff 11 D cm 11 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 1 100000 1e+10 1e+15 1e+20 1e+25 m 2 /day Number of particles D adv 22 D eff 22 D cm 22 Fig.11. Convergence as function of N of the effective coefficients, defined as time averages between t = 4900 and t = 5000 days. The Fig. 12 shows that small numbers of particles induce large numerical errors in simulation of time behavior of the coefficients.
Numerical Modeling of Large Scale Transport of Contminant Solutes 167 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 D eff 11 (N) (m 2 /day) Ut/λ y N=4 N=10e02 N=10e06 N=10e10 N=10e24 Fig. 12. Time behavior for different N . Then, one can assess that the large scale computations with PT procedure, reported in the literature, did not ensure the convergence because a too small number of particles was used (N = 4 in [Schwarze et al., 2001] and N ∼ 100 in [Dentz et al., 2002 and 2003]). Moreover, as shown by Fig. 2, the PT algorithm is not able to handle N = 10 10 particle even if efficient parallel computes are used. The convergence of large scale transport simulation in given realizations of the velocity field seems to be a new result, for the first time presented here. 7 Ensemble averages of the effective coefficients The convergence of the ensemble averages (more precisely, of the advective coefficients <D adv ii > V , neglecting the local dispersion, using the PT procedure with a single particle in every realization of the velocity field and 1500 realizations) was obtained in the past, for tens of correlation lengths [Bellin et al., 1992]. Here, the convergence of the statistical estimation of effective coefficients (2) is investigated for large scales, corresponding to 5000 correlation lengths. For every realization of the transport (i.e. diffusion process in a fixed realization of the velocity field) the statistical convergence was ensured using N = 10 10 particles. With the chosen grid steps, δx =0.25 m and δt =1 day , it was possible to obtain a convergent transport computation on a single processor of the parallel machine Cray T3E and 256 realizations of the transport in a single run, using a GRW code. Repeating the job 10 times, an ensemble of 2560 realizations was obtained. The mean values and the standard deviations of the longitudinal and transverse coef- ficients D eff 11 and D eff 22 given by (2), as well as for the corresponding advective and center of mass coefficients, were computed at every time step. The relative errors presented in Fig. 13, computed as averages over the time interval between 4900 days and 5000 days, show that the estimation of the stochastic averages with a precision between 1 and 5% requires more than S = 2500 realizations of the velocity field.
168 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 500 1000 1500 2000 2500 σ D adv 11 / <D adv 11 > v σ D adv 22 /< D adv 22 > v 0 0.005 0.01 0.015 0.02 0.025 0.03 500 1000 1500 2000 2500 σ D eff 11 / <D eff 11 > v σ D eff 22 / <D eff 22 > v 0 0.1 0.2 0.3 0.4 0.5 0.6 500 1000 1500 2000 2500 Number of realizations σ D cm 11 / <D cm 11 > v σ D cm 22 / <D cm 22 > v Fig. 13. Convergence of the stochastic averages of the effective coefficients as function of number of realizations S . The averaged longitudinal coefficients are of the same order as predicted by the first order approximation (3) but the ergodic behavior does not occur (the center of mass coefficient <D cm 11 > V does not vanish) (Fig. 14a). The ergodicity of the transverse coefficients occurs at travel times of thousands days and their asymptotic values are larger than D (Fig. 14b).
Numerical Modeling of Large Scale Transport of Contminant Solutes 169 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 m 2 /day Ut/λ y D eff 11 D cm 11 D+D adv 11 Fig. 14a. Ensemble averaged longitudinal effective diffusion coefficients. 0 0.005 0.01 0.015 0.02 0.025 1 10 100 1000 m 2 /day Ut/λ y D eff 22 D cm 22 D+D adv 22 Fig. 14b. Ensemble averaged transverse effective diffusion coefficients. The non-ergodicity of the longitudinal coefficients is due to the nonvanishing value of the center of mass coefficient in the average of (2). The average value of this coefficient monotonically decreases with S , for S> 512, as shown in Fig. 15. After S = 2048 it reaches a quasi constant value <D cm 11 > V ≃ 0.02 m 2 /day estimated, according to Fig. 13, with a precision of about 5%.
170 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken 0.019 0.02 0.021 0.022 0.023 0.024 0.025 0.026 3000 3500 4000 4500 5000 1/S Σ s=1,S ( D cm 11 ) (m 2 /day) Ut/λ y S=512 S=767 S=1024 S=1280 S=1536 S=1792 S=2048 S=2304 S=2560 Fig. 15a. Time behavior of the averaged center of mass coefficient D cm 11 for different numbers of realizations. 0 0.005 0.01 0.015 0.02 0.025 0 500 1000 1500 2000 2500 3000 (m 2 /day) S 1/S Σ s=1,S ( D cm 11 ) 1/S Σ s=1,S ( D cm 22 ) x 10 Fig. 15b. Dependence of the averaged center of mass coefficients D cm 11 and D cm 22 at t = 5000 days, on number of realizations. The average <D cm 22 > V reaches a constant value two orders of magnitude smaller than the local coefficient D, after S = 1024 realizations (Fig. 15b), with the same precision of 5%. The effective coefficient <D eff 22 > V ≃ D+ <D adv 22 > V ≃ 0.014 m 2 /day is estimated, according to Fig. 13, with a precision ≃ 1%. This result is at variance with the 1-st order approximation (3) which predicts a vanishing advective coefficient, <D adv 22 > V = 0, and an effective coefficient equal to the local coefficient, <D eff 22 > V = D =0.01 m 2 /day .
Numerical Modeling of Large Scale Transport of Contminant Solutes 171 These results are statistically reliable but the discretization errors, estimated from Fig. 10, are of the same order of magnitude as the differences with respect to the 1-st approximation. Even if the prediction of the 1-st order stochastic theory concerning the ergodicity and the values of the asymptotic coefficients cannot yet be verified, the results on ensemble averages of the effective coefficients presented in Fig. 14 clarify some other aspects. First, because the GRW numerical approach permits the computation of the terms composing the effective coefficients (2), the explanation of the significant difference in the behavior at intermediate times of effective and advective transverse coefficients, already noted in [Dentz et al., 2002], becomes evident in Fig. 14b. The difference < D adv 22 > − <D eff 22 > is just the <D cm 22 > coefficient, due to the randomness of the plume center of mass. The maximum of the advective coefficient <D adv 22 > at about ten dimensionless times, predicted by stochastic theory [Dagan, 1989], is an artifact due to the ensemble averaging and do not describe the behavior of the transport in given aquifers. The estimated effective coefficient <D eff 22 > behaves almost monotonically. A second conclusion, statistically reliable and not affected by the discretization errors is that the upscaled Fickian behavior can be expected to occur only after thousands of dimensionless times. This reinforces previous numerical results, obtained by averaging over transport realizations that were not convergent [Schwarze et al., 2001; Dentz et al. 2002, 2003], and contradicts the theoretical predictions saying that the asymptotic coefficients are reached after the plume have traveled tens correlations lengths [Dagan, 1989; Fiori, 1996]. 8 Concentration fluctuations and dilution The normalized concentration is the probability density of molecules moving on the trajectories X (t, x 0 ) of the transport process described in Section 3 and it is given by the general definition c(x,t)= 〈δ (x − X (t, x 0 ))〉 w,x 0 [van Kampen, 1981]. The average of the concentration over the realizations of the velocity field can be written 〈c(x,t)〉 V = 〈δ (x − X (t, x 0 ))〉 w,x 0 V = V 0 〈δ (x − X (t, x 0 ))〉 w,V c(x 0 )dx 0 , (7) where V 0 is the domain occupied by the initial plume and c(x 0 ) is the initial concentration field. The average of the square concentration is c 2 (x,t) V = V 0 V 0 〈δ (x − X (t, x 0 1 ))〉 w 〈δ (x − X (t, x 0 2 ))〉 w V c(x 0 1 )c(x 0 2 )dx 0 1 dx 0 2 (8a) = V 0 V 0 〈δ (x − X (t, x 0 1 ))δ (x − X (t, x 0 2 ))〉 w,V c(x 0 1 )c(x 0 2 )dx 0 1 dx 0 2 (8b) = V 0 V 0 p(x,t; x 0 1 , 0; x,t; x 0 2 , 0)c(x 0 1 )c(x 0 2 )dx 0 1 dx 0 2 , (8c) where (8a) and (8b) are identical because the trajectories of the Wiener process starting in two different points x 0 1 and x 0 2 are independent. The function p in (8c) is the joint density of four events (x,t) laying on two trajectories, sometimes called “two particles probability density” [see e.g. Lundgren, 1981]. This joint density is used to build up two- trajectory estimators for correlation function and variance, introduced in [Sabelfeld and
172 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken Mikhailov, 1980], which are useful tools in Monte Carlo methods for turbulent transport in random velocity fields [Sabelfeld, 1991]. From (7) and (8), the variance of the concentration is σ 2 c = 〈c 2 (x,t)〉 V − (〈c(x,t)〉 V ) 2 . In the Lagrangian approach [Dagan and Fiori, 1997] first order results are provided, for which probability density factorizes, p = p V p w , and both p V and p w are Gaussian. In the Eulerian approach [Kapoor and Gelhar, 1994, Kapoor and Kitanidis, 1998], σ 2 c is derived by taking the moments of an advection-diffusion equation supposed to be valid at the Darcy scale. The concentration statistics is usually described by the variance σ 2 c , the concentration coefficient of variation CV = σ c / 〈c(x,t)〉 V , computed at the plume center of mass, and the global variance ||σ 2 c || = σ 2 c dx 1 dx 2 . Besides the averages over the ensemble of realizations given by (7-8), the risk of contamination in a given aquifer should also be quantified by the space mean concentration and variance [Kapoor and Gelhar, 1994]. The dilution of the contaminant solute in given realization of the transport is described by the “dilution index” E [Kapoor and Kitanidis, 1998], E = exp − c(x,t) ln c(x,t)dx , (9) where the expression under the exponential is the entropy of the process and the integral extends over the entire problem domain. The simulations of the two-dimensional transport in unbounded domain and point source from previous sections were supplemented with simulations for extended initial plumes, with and without local dispersion, and with a simulation of the transport in a confined aquifer (for initial plume uniformly distributed over the vertical direction and no flux at upper and lower boundaries, similar to [Kapoor and Kitanidis, 1998]). In all computations the initial plume contained 10 10 particles and the concentration was computed as the relative number of particles in a square meter. Space statistics were inferred from averages over the transverse dimension of the computation domain and, for the point source case, averages over 500 realizations of the velocity were used to estimate the ensemble averages. The results for unbounded domains from Fig. 16 are similar to ones published in the past and underline the crucial role played by the local dispersion (unbounded increase of global variance an CV in absence of local dispersion). In Fig. 16b the ensemble averaged numerical values σ 2 c / 〈c(x,t)〉 V are represented by the points “+”. They correspond to the theoretical CV inferred from the statistics over realizations given by (7-8) and are smaller than the space average CV . To avoid the underestimation of the incertitude of the spatial distribution of the contaminant in a given realization of the aquifer the space averages should be considered. Further, the statistical CV accounts for the incertitude of the aquifer’s parameters (hydraulic conductivity).
Numerical Modeling of Large Scale Transport of Contminant Solutes 173 0 1e-05 2e-05 3e-05 4e-05 5e-05 6e-05 7e-05 8e-05 9e-05 0.0001 100 1000 10000 ||σ c 2 || Ut/λ y Point source Band, D=0, H=100 m, L=10 m Band, D=0.01, H=10 m, L=10 m Fig. 16a. Time behavior of the global variance. 0 2 4 6 8 10 12 14 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 CV Ut/λ y Point source Band, D=0, H=100 m, L=10 m Band, D=0.01, H=10 m, L=10 m Point source, σ C / <C> v Fig. 16b. Time behavior of the concentration coefficient of variation. The dilution index (9), normalized by E max , corresponding to Gaussian distribution, for unbounded domains, and Gaussian distribution on longitudinal direction and uniform on transverse direction, for confined aquifer, is presented in Fig. 17. For both, the behavior is similar, towards E/E max = 1 (occurrence which is called “complete dilution” in [Kapoor and Kitanidis, 1998]), while in the absence of a local dispersion mechanism there is no dilution.
174 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 1 10 100 1000 10000 E / E max Ut/λ y ~[Kapoor and Kitanidis, 1998], no flux at H=0 and 100 m Point source Extended, D=0, σ y =12 m Linear, D=0, H=100 m, L=0 m Band, D=0, H=100 m, L=10 m Band, D=0.01, H=10 m, L=10 m Fig. 17. Normalized dilution index. Fig. 18 shows that for both boundary problems the dilution index E is proportional to the apparent dimension of the plume σ x 1 σ x 2 . 0 1000 2000 3000 4000 5000 6000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 E Ut/λ y mumeric fit with E=at 0.5 10σ x 1 σ x 2 Fig. 18a. Time behavior of the dilution index E compared to the apparent dimension of the plume, for unbounded domain.
Numerical Modeling of Large Scale Transport of Contminant Solutes 175 0 2000 4000 6000 8000 10000 12000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 E Ut/λ y mumeric fit with E=at 0.5 10σ x 1 σ x 2 Fig. 18b. Time behavior of the dilution index E compared to the apparent dimension of the plume, for confined aquifer. Since the apparent plume behaves at large times as ∼ t for unbounded domains (Fig. 18a) and as ∼ t 0.5 for confined aquifer (Fig. 18b), the observation of Pannone and Kitanidis [1999] that the results in [Dagan and Fiori, 1997] are in variance with those in [Kapoor and Kitanidis, 1998] is not really surprising. The two results are not comparable because they are derived from different boundary problems. For both transport processes there is no thermodynamic equilibrium, because the entropy increases unboundedly. 9 Conclusions The GRW algorithm is able take over the limitations of classic particles methods and to track very large numbers of particles necessary to achieve the statistical convergence and to compute accurate concentrations in large scale transport simulations. The computation parameters used made possible the simulation of the asymptotic diffusive behavior of the transport in given realizations of the velocity field and the com- putation of statistically reliable averages of the effective coefficients over large ensembles of realizations. The simulated effective coefficients are close to the prediction of the first order approx- imation in the stochastic theory of transport. Reinforcing previously published results, the time necessary to reach the asymptotic Fickian regime of transport was found to be thousands times larger than the theoretical prediction. The numerical errors are still to large to test the validity of the ergodic hypothesis and the prediction of vanishing advec- tive contribution to the transverse effective coefficient in the first order theory. It was shown that the results can be improved by decreasing the space and time steps so that the overshooting is kept at small values and the numerical resolution of the velocity field is improved. The mean concentration field and fluctuations were accurately described by space and ensemble averages. The comparison between simulations in different initial and boundary
176 N. Suciu, C. Vamo¸s, H. Hardelauf, J. Vanderborght, and H.Vereecken conditions shows the determinant role of the local dispersion mechanism in fluctuations extinction. As shown by the long time behavior of the entropy, the transport process has no equilibrium state and the contamination approaches a state of complete dilution. References Ames, W. F., Numerical Methods for Partial Differential Equations, 2nd ed., Academic Press, New York, 1977. Bellin, A., P. Salandin and A. Rinaldo, Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations, Water Re- sour. Res. 28(9), 2211-2227, 1992. Chorin, A. J., Vortex sheet approximation of boundary layers, J. Comput. Phys. 27, 428-442 (1978). Dagan, G., Flow and Transport in Porous Formations, Springer, New York, 1989. Dagan, G. and A. Fiori, The influence of pore-scale dispersion on concentration statistical moments in transport through heterogeneous aquifers, Water Resour. Res., 33, 1595-1650, 1997. Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach, Temporal behavior of a solute cloud in a heterogeneous porous medium 3. Numerical simulations, Water Resour. Res., 38, 23:1-12, 2002. Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach, Numerical studies of the transport behavior of a passive solute in a two-dimensional incompressible random flow field, Phys. Rev. E 67, 046306:1-10, 2003. Fiori, A., Finite Peclet extension of Dagan’s solutions to transport in anisotropic forma- tions, Water Resour. Res., 32, 193-198, 1996. Kapoor, C. and L. W. Gelhar, Transport in three-dimensionallity heterogeneous aquifers 1. Dynamics of concentration fluctuations, Water Resour. Res., 30, 1775-1789, 1994. Kapoor, C. and P. K. Kitanidis, Concentration fluctuation and dilution in aquifers, Water Resour. Res., 34, 1181-1193, 1998. Lundgren, T. S., Turbulent pair dispersion and scalar diffusion, J. Fluid Mech., 111, 27-57, 1981. Pannone, M. and P. K. Kitanidis, Large-time behavior of concentration variance and dilution in heterogeneous media, Water Resour. Res., 35, 623-634, 1999. Roth, K. and K. Hammel, Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with steady state flow, Water Resour. Res., 32, 1653-1663, 1996.
Numerical Modeling of Large Scale Transport of Contminant Solutes 177 Sabelfeld, K. K. and G. A. Mikhailov, On the numerical simulation of the particle’s diffusion in random velocity fields, Izv. AN SSSR, Physics of Atmosphere and Ocean, 16 (3), 229-235 (in Russian), 1980. Sabelfeld, K. K., Monte Carlo methods in boundary value problems, Springer, New York - Heidelberg - Berlin, 1991. Schwarze, H., U. Jaekel and H. Vereecken, Estimation of macrodispersivity by different approximation methods for flow and transpot in randomly heterogeneous media, Transport in Porous Media 43, 265-287, 2001. Suciu, N., C. Vamo¸s, H. Vereecken, and J.Vanderborght, pp. 111-140 in G. Marinovschi and I. Stelian (Eds.) Proceedings of the first Workshop on Mathematical Modeling of Environmental Problems, Rom. Acad. Publishing House, Bucharest 2002. Suciu N., C. Vamo¸s, H. Vereecken, J. Vanderborght, and H. Hardelauf, pp. 215-221 in T. Maghiar, A. Georgescu, M. Balaj, I. Dzit ¸ac, I. Mang (Eds.), Proc. of the 11th Conf. of Appl. and Ind. Math., Vol. 1, Oradea Univ. Publishing House, Oradea, 2003. Suciu, N. and C. Vamo¸s, Effective diffusion in heterogeneous media, Interner Bericht ICG-IV 00303, Forschungszentrum J¨ ulich, 2003. Sun, Ne-Z., Mathematical Modeling in Groundwater Pollution, Springer, New York, 1996. Vamo¸s, C., N. Suciu, H. Vereecken, J. Vanderborht and O. Nitzsche, Path decom- position of discrete effective diffusion coefficient, Internal Report ICG-IV.00501, Forschungszentrum J¨ ulich, 2001. Vamo¸s, C., N. Suciu, and H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186, 527-544, 2003. van Kampen, N. G., Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1981.
2004

Related Posts