Global random walk modeling of transport in complex systems

Nicolae Suciu
(original), paper

Abstract

The Global random walk algorithm performs simultaneously the tracking of large collections of particles and permits massive simulations at reasonable costs. Applications were developed for transport in systems with anisotropic, non-homogeneous, and randomly distributed parameters. As a first illustration we present simulations for diffusion in human skin. Further, a case study for contaminant transport in groundwater shows that the realizations of the transport process converge in mean square limit to a Gaussian diffusion. This investigation also indicates that the use of the Kraichnan routine, based on periodic random fields, yields reliable simulations of transport in Gaussian velocity fields.

Authors

N. Suciu
Tiberiu Popoviciu Institute of Numerical Analysis, Cluj Napoca branch of the Romanian Academy

C. Vamoş
Tiberiu Popoviciu, Institutue of Numerical and analysis, Romanian Academiy

I. Turcu
National R&D Institute for Isotopic and Molecular Technologies, Cluj-Napoca, Romania

C.V.L. Pop
National R&D Institute for Isotopic and Molecular Technologies, Cluj-Napoca, Romania

L.I. Ciortea
National R&D Institute for Isotopic and Molecular Technologies, Cluj-Napoca, Romania

Keywords

Global random walk; lattice gas; human skin; groundwater contamination

References

see the expansion block below.

Cite this paper as:

N. Suciu, C. Vamoş, I. Turcu, C.V.L. Pop, L.I. Ciortea (2007), Global random walk modeling of transport in complex systems, Computing and Visualization in Science, doi: 10.1007/s00791-007-0077-6

PDF

2009 Vamos Global random walk modeling of transport in complex systems

About this paper

Journal

Computing and Visualization in Science

Publisher Name
DOI

10.1007/s00791-007-0077-6

Print ISSN

Not available yet.

Online ISSN

Not available yet.

Google Scholar Profile

google scholar link

  • [1] Eberhard, J., Approximations for transport parameters and selfaveraging properties for point-like injections in heterogeneous media. J. Phys. A Math. Gen. 37, 2549–2571 (2004)[2] Eberhard, J., Suciu N., Vamos, On the self-averaging of dispersion for transport in quasi-periodic random media. J. Phys. A: Math. Gen. 40, 597–610, doi:10.1088/1751-8113/40/4/002 (2007)[3] Jaekel, U., Vereecken, H.: Renormalization group analysis of macrodispersion in a directed random flow. Water Resour. Res. 33, 2287–2299 (1997)
    [4] Karapiperis, T., Blankleider, B., Cellular automaton model of reaction-transport processes. Physica D 78, 30–64 (1991)

    [5] Kesten, H., Papanicolaou, G.C., A limit theorem for turbulent diffusion. Commun. Math. Phys. 65, 97–128 (1979)

    [6] Kinzelbach, W., Uffink, G. , The random walk method and extensions in groundwater modelling. In: Bear, J., Corapcioglu, M.Y. (eds.) Transport Processes in Porous Media, pp. 761–787. Kluwer, Norwell (1991)

[7] Johnson, M.E., Blankschtein, D., Langer, R., Evaluation of solute permeation through the stratum corneum: lateral bilayer diffusion as the primary transport mechanism. J. Pharm. Sci. 86(10), 1162– 1172 (1997)

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

[9] Schwarze, H., Jaekel, U., Vereecken, H., Estimation of macrodispersivity by different approximation methods for flow and transport in randomly heterogeneous media. Transp. Porous Media 43, 265– 287 (2001)

[10] Suciu, N., Vamos, C., Vanderborght, J., Hardelauf, H., Vereecken, H., Numerical modeling of large scale transport of contminant solutes using the global random walk algorithm. Monte Carlo Methods Appl. 10(2), 153–177 (2004)

[11] Suciu, N., Vamos, C., Knabner, P., Rüde, U., Biased global random walk, a cellular automaton for diffusion. In: Hüsemann, F., Kowarschik, M., Rüde, U. (eds.) Simulationstechnique, 18th Symposium in Erlangen, September 2005., pp. 562–567. SCS Publishing House e. V., Erlangen (2005)

[12] Suciu, N., Vamos, C., Evaluation of overshooting errors in particle methods for diffusion by biased global random walk. Rev. Anal. Num. Th. Approx. 35, 119–126 (2006)

[13] Suciu, N., Vamos, C., Vanderborght, J., Hardelauf, H., Vereecken, H., Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res. 42, W04409. doi:10.1029/2005WR004546 (2006)

[14] Suciu, N., Vamos, C., Eberhard, J., Evaluation of the firstorder approximations for transport in heterogeneous media. Water Resour. Res. 42, W11504. doi:10.1029/2005WR004714 (2006)

[15] Vamos, C., Suciu, N., Vereecken, H., Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys. 186(2), 527–544 (2003

[16] Vamos, C., Suciu, N., Turcu, I., Pop, C.V.L., Ciortea, L. I.,  Programme BIOTECH, Project No. 01-8-CPD-042/19.10.2001. Research report (2004, in Romanian)

[17] Vereecken, H., Döring, U., Hardelauf, H., Jaekel, U., Hashagen, U., Neuendorf, O., Schwarze, H., Seidemann, R.,  Analysis of solute transport in a heterogeneous aquifer: the Krauthausen field experiment. J. Contam. Hydol. 45, 329–358 (2000)

[18] Wilke, J., Pohl, T., Kowarschik, M., Rüde, U., Cache performance optimizations for parallel lattice Boltzmann codes. In: Proceedings of the EuroPar-03 Conference. Lecture Notes in Computer Science, vol. 2790, pp. 441–450. Springer, Heidelberg (2003)

[1] Eberhard, J., Approximations for transport parameters and selfaveraging properties for point-like injections in heterogeneous media. J. Phys. A Math. Gen. 37, 2549–2571 (2004)

[2] Eberhard, J., Suciu N., Vamos, On the self-averaging of dispersion for transport in quasi-periodic random media. J. Phys. A: Math. Gen. 40, 597–610, doi:10.1088/1751-8113/40/4/002 (2007)

[3] Jaekel, U., Vereecken, H., Renormalization group analysis of macrodispersion in a directed random flow. Water Resour. Res. 33, 2287–2299 (1997)

[4] Karapiperis, T., Blankleider, B., Cellular automaton model of reaction-transport processes. Physica D 78, 30–64 (1991)

[5] Kesten, H., Papanicolaou, G.C., A limit theorem for turbulent diffusion. Commun. Math. Phys. 65, 97–128 (1979)

[6] Kinzelbach, W., Uffink, G., The random walk method and extensions in groundwater modelling. In: Bear, J., Corapcioglu, M.Y. (eds.) Transport Processes in Porous Media, pp. 761–787. Kluwer, Norwell (1991)

[7] Johnson, M.E., Blankschtein, D., Langer, R., Evaluation of solute permeation through the stratum corneum: lateral bilayer diffusion as the primary transport mechanism. J. Pharm. Sci. 86(10), 1162– 1172 (1997)

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

[9] Schwarze, H., Jaekel, U., Vereecken, H., Estimation of macrodispersivity by different approximation methods for flow and transport in randomly heterogeneous media. Transp. Porous Media 43, 265– 287 (2001)

[10] Suciu, N., Vamos, C., Vanderborght, J., Hardelauf, H., Vereecken, H., Numerical modeling of large scale transport of contminant solutes using the global random walk algorithm. Monte Carlo Methods Appl. 10(2), 153–177 (2004)

[11] Suciu, N., Vamos, C., Knabner, P., Rüde, U., Biased global random walk, a cellular automaton for diffusion. In: Hüsemann, F., Kowarschik, M., Rüde, U. (eds.) Simulationstechnique, 18th Symposium in Erlangen, September 2005., pp. 562–567. SCS Publishing House e. V., Erlangen (2005)

[12] Suciu, N., Vamos, C., Evaluation of overshooting errors in particle methods for diffusion by biased global random walk. Rev. Anal. Num. Th. Approx. 35, 119–126 (2006)

[13] Suciu, N., Vamos, C., Vanderborght, J., Hardelauf, H., Vereecken, H., Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res. 42, W04409. doi:10.1029/2005WR004546 (2006)

[14] Suciu, N., Vamos, C., Eberhard, J., Evaluation of the firstorder approximations for transport in heterogeneous media. Water Resour. Res. 42, W11504. doi:10.1029/2005WR004714 (2006)

[15] Vamos, C., Suciu, N., Vereecken, H., Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys. 186(2), 527–544 (2003)

[16] Vamos, C., Suciu, N., Turcu, I., Pop, C.V.L., Ciortea, L. I., Programme BIOTECH, Project No. 01-8-CPD-042/19.10.2001. Research report (2004, in Romanian)

[17] Vereecken, H., Döring, U., Hardelauf, H., Jaekel, U., Hashagen, U., Neuendorf, O., Schwarze, H., Seidemann, R., Analysis of solute transport in a heterogeneous aquifer: the Krauthausen field experiment. J. Contam. Hydol. 45, 329–358 (2000)

[18] Wilke, J., Pohl, T., Kowarschik, M., Rüde, U., Cache performance optimizations for parallel lattice Boltzmann codes. In: Proceedings of the EuroPar-03 Conference. Lecture Notes in Computer Science, vol. 2790, pp. 441–450. Springer, Heidelberg (2003)

Comput Visual Sci (2009) 12:77–85 DOI 10.1007/s00791-007-0077-6 REGULAR ARTICLE Global random walk modelling of transport in complex systems N. Suciu · C. Vamo¸ s · I. Turcu · C. V. L. Pop · L. I. Ciortea Received: 9 February 2005 / Accepted: 8 December 2006 / Published online: 5 October 2007 © Springer-Verlag 2007 Abstract The Global random walk algorithm performs simultaneously the tracking of large collections of particles and permits massive simulations at reasonable costs. Appli- cations were developed for transport in systems with anisotropic, non-homogeneous, and randomly distributed parameters. As a first illustration we present simulations for diffusion in human skin. Further, a case study for conta- minant transport in groundwater shows that the realizations of the transport process converge in mean square limit to a Gaussian diffusion. This investigation also indicates that the use of the Kraichnan routine, based on periodic random fields, yields reliable simulations of transport in Gaussian velocity fields. Keywords Global random walk · Lattice gas · Human skin · Groundwater contamination Communicated by G. Wittum. N. Suciu (B) Institute of Applied Mathematics, Friedrich-Alexander University of Erlangen-Nuremberg, Martensstrasse 3, 91058 Erlangen, Germany e-mail: suciu@am.uni-erlangen.de; n.suciu@fz-juelich.de N. Suciu · C. Vamo¸ s “Tiberiu Popoviciu” Institute of Numerical Analysis, Cluj Napoca Branch of the Romanian Academy, P.O. Box 68-1, 400320 Cluj-Napoca, Romania I. Turcu · C. V. L. Pop · L. I. Ciortea National R&D Institute for Isotopic and Molecular Technologies, P.O. Box 700, 400293 Cluj-Napoca, Romania L. I. Ciortea School of EEIE, South Bank University, 103 Borough Road, London SE1 0AA, UK 1 Introduction The global random walk (GRW) is a recently proposed method for the simulation of transport processes which generalizes the classical Particle Tracking algorithm (PT). The new approach considerably increases the accuracy and reduces the computation time. Non-isotropic and space variable diffusion coefficients as well as various boundary and initial conditions can be easily implemented in GRW [15]. Unlike PT, GRW acts “globally” since, starting with a given distribution of particles in a computational grid, all the particles lying at a grid site are simultaneously spread, first by an advection displacement, then by diffusion jumps, as shown in Fig. 1. The numbers of particles undergoing dif- fusion jumps are Bernoulli random variables, describing the number of jumps in a given direction of a discrete random walk process. Thus, the GRW procedure is equivalent to a superposition of many PT procedures and can manage large numbers of particles. GRW is also equivalent to a finite- difference scheme, for genuine diffusion processes, and it is appropriate to solve more complex problems for para- bolic equations. The “overshooting” cannot be avoided in GRW procedure but, for given overshooting, the discretiza- tion errors can be controlled and kept under a desirable thre- shold by an accurate discretization of the velocity [10, 13]. The “global change of state” of the whole system of com- putational particles, as a cellular automaton, is the principle on which are built the “lattice gas” and “lattice Boltzmann” methods, and it allows high performance parallel computing applications [18]. In this respect, GRW is a particular lat- tice gas model, i.e. it is a stochastic process in the space of configurations, defined at a given time by the number of par- ticles at a grid site. Since the particles do not interact, there is no collision step, as for lattice Boltzmann. The “particle distribution function” along all the possible positions 123
78 N. Suciu et al. 0 0.1 0.2 0.3 0.4 x 1 (m) -0.1 0 0.1 x 2 (m) 1e+09 1e+10 n(x 1 ,x 2 ) Fig. 1 The advection displacement and diffusive jumps of 10 10 par- ticles starting at (0,0), for r = 0.85; the blue columns and lines represent the succession of the one-dimensional procedures occupied by particles after a “propagation” step consisting of an advection displacement and diffusion jumps is given by the Bernoulli distribution or its approximations. Similarly with gas lattice, the next configuration of the system is obtai- ned by simple summations of the number of particles at each grid site. Since in GRW algorithm, unlike in lattice gas models, the number of particles per grid site is not limited by an “exclu- sion principle” it is possible to obtain statistically reliable estimations of the macroscopic quantities, in a single simu- lation, without averaging over ensembles of simulations. The required number of particles ranges between 10 6 , for small size lattices of the order of 100 nodes [15] to 10 10 for lat- tice size of ∼10 6 nodes [10, Fig. 11]. In some simplified implementations of the algorithm there are no limitations as concerning the number of particles. Therefore, GRW can simulate the behavior of the real number of molecules invol- ved in microscopic processes (diffusion, chemical reactions and radioactive decay) as well as in large-scale environmen- tal problems (contaminant transport in water or atmosphere). The paper is organized as follows. Section 2 presents the principle of the GRW algorithm. A two-dimensional algo- rithm for space variable and anisotropic diffusion tensor is described in Sect. 3, and in Sect. 4 it is applied for simu- lations of diffusion in human skin. Section 5 deals with two-dimensional simulation of large-scale transport in groundwater and Sect. 6 presents some conclusions. 2 The GRW algorithm The GRW algorithm simulates the random walk of an ensemble of fictitious particles to solve parabolic partial derivative equations. As an illustration, we first present the GRW algorithm for the case of the one-dimensional diffu- sion whose probability density (normalized concentration of the diffusing substance) is a solution of the Fokker–Planck equation with space variable drift coefficient V (x ) and constant diffusion coefficient D, ∂ t c + V (x )∂ x c = D∂ 2 x c. (1) Let us consider a lattice with constant δ x , defined by the nodes x i = i δ x , and the discrete time t k = k δt , where i and k , k ≥ 0, are integers. The configuration of N particles evolving on the lattice is given by the number n(i , k ) of particles lying at x i at the time step t k . At every time step the particles from every lattice site are grouped into three subsets: n( j , k ) = δn( j + v j | j , k ) + δn( j + v j − d | j , k ) + δn( j + v j + d | j , k ). (2) The term δn( j + v j | j , k ) is the number of particles which are moved from x j to x j + v j , where v j = V (x i )δt /δ x . The terms δn( j + v j ± d | j , k ) in Eq. (2) count the rest of particles, which undergo diffusive right–left jumps over d lattice nodes and are Bernoulli random variables. The new configuration of the particles after a time step is given by the sum of the particles which at the time (k + 1)δt are at the same lattice site, n(i , k + 1) = j δn(i | j , k ). (3) By repeating this procedure for different sets of realiza- tions of the Bernoulli variables δn one obtains a statistical ensemble of simulations. The ensemble average of the num- ber of particles undergoing diffusive jumps and of the number of particles remaining at the same site after the displacement v j are given by the relations δn( j + v j ± d | j , k ) = 1 2 r n( j , k ), δn( j + v j | j , k ) = (1 − r ) n( j , k ), (4) where r ,0 ≤ r ≤ 1 is a real parameter. The diffusion coef- ficient D is related to the amplitude of diffusive jumps d , the parameter r , the lattice constant and the time step by the relation D = r (d δ x ) 2 2δt . (5) The concentration simulated by the one-dimensional GRW algorithm is given by c(x i , t k ) = n(i , k )/δ x . For vani- shing or constant drift coefficients V , the relations (3–5) imply that the average of c(x i , t k ) over an ensemble of GRW simulations coincides with the solution of the stable finite dif- ference scheme for the diffusion equation [15]. For variable V , GRW is no longer equivalent to a finite difference scheme. This is due to overshooting errors occurring when particles jump over grid sites with different V (x i ) and is common to most of the random walk methods. Nevertheless, an evalua- tion by comparisons with a method free of overshooting [11] 123
Global random walk modelling of transport in complex systems 79 has shown that, even for highly variable V , the GRW solu- tion remains accurate in a range of 5% for suitable chosen parameters δ x , δt , d , and r [12]. In the “exact” version of the GRW the quantities δn( j + v j ± d | j , k ) are picked up from Bernoulli repartition func- tions, using computer generated random numbers. This pro- cedure requires the expression of the number of jumping particles in a series of powers of 2, the computation and the storage of the Bernoulli repartitions for increasing power of 2, and the generation of some random numbers at every time step and lattice site. Accordingly, the computing time increases for large N . Several approximations of the exact algorithm can be used to avoid this inconvenience. For ins- tance, if the number of particles at a lattice site is larger than 2 20 , the Bernoulli repartition can be successfully approxima- ted by the Gaussian one (i.e. with the erf-function). Another useful and robust approximation is that leading to the “redu- ced fluctuations GRW algorithm”. In both cases there are no limitations for the maximum number of particles, which can be as large as the maximum double precision value of the computing platform. The reduced fluctuations GRW algorithm 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 + v j | j , k ), [n/2] is the inte- ger part of n/2 and θ is a variable taking the values 0 and 1 with probability 1/2. This algorithm minimizes the num- ber of calls of the random number generator and makes the computing time independent of N . Therefore the reduced fluctuations GRW algorithm is mainly useful for large-scale problems. Moreover, a physical significant shape of the dif- fusion front is ensured by the threshold concentration defi- ned by one particle at a grid point. This feature also saves memory and computing time by avoiding computations at lattice sites with irrelevant small concentrations, as in case when one gives up the particle indivisibility or in finite dif- ferences schemes [15]. In the following by GRW we shall denote only reduced fluctuations algorithms. For constant diffusion coefficient D the two- and three- dimensional algorithms can be simply built by repeating the same procedure for all space directions. Figures 1 and 2 illus- trate the two-dimensional GRW algorithm for constant D. The concentration plotted in Fig. 2 was computed at a given time t = k δt and at a lattice site (x 1 , x 2 ) = (i 1 δ x 1 , i 2 δ x 2 ), where δt is the time step, δ x 1 , and δ x 2 are the space steps, with the formula c(x 1 , x 2 , t ) = 1 N 1 δ x 1 δ x 2 i 1 +1 i 1 =i 1 −1 i 2 +1 i 2 =i 2 −1 n ( i 1 , i 2 , k δt ) . (6) Let us consider, more generally, an initial condition consis- ting of N particles uniformly distributed over N X 0 grid sites. 0.25 0.2 0.15 0.1 0.05 -2 0 2 4 6 8 10 12 14 x 1 -4 -3 -2 -1 0 1 2 3 4 x 2 0 0.05 0.1 0.15 0.2 0.25 0.3 c(x 1 ,x 2 ) Fig. 2 The resulting concentration field, at t = 0, t = 5δt and t = 10δt , for a uniform initial distribution of 10 10 particles in the square of sides δ x 1 and δ x 2 By n(i 1 , i 2 , k ; x 0,1 , x 0,2 ) we denote the distribution of par- ticles at the time step k given by the GRW procedure for a diffusion process staring at (x 0,1 , x 0,2 ). Since the distribution of the particles at time k can be written as n(i 1 , i 2 , k ) = x 0,1 ,x 0,2 n(i 1 , i 2 , k ; x 0,1 , x 0,2 ), it follows that 1 N i 1 ,i 2 n(i 1 , i 2 , k ) = 1 N X 0 x 0,1 ,x 0,2 N X 0 N i 1 ,i 2 n(i 1 , i 2 , k ; x 0,1 , x 0,2 ). Thus, the concentration (6), as well as its spatial moments, are averages over the trajectories of the diffusion process starting at given initial positions and over the distribution of the initial positions. It was shown that for Gaussian diffusion the numerical solution (6) converges to the Gaussian distribution 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 the finite diffe- rences scheme [15]. Therefore, GRW is “self-averaging”, in the sense that for large N the solution given by a single simu- lation obeys the same equation (the finite difference scheme) as the ensemble average of Eq. (6). The self-averaging pro- perty distinguishes GRW from the gas-lattice approaches, where the “exclusion principle” does not permit the pre- sence large numbers of identical particles at the same site and the solution of the problem is obtained by averaging over ensembles of simulations. 123
80 N. Suciu et al. The algorithm is obviously stable since the total number of particles N contained in the lattice is conserved. The condi- tion r ≤ 1, ensures that there is no numerical diffusion. In [15] it was also shown that GRW is orders of magnitude fas- ter than PT, in all its possible implementations, and is able to handle large numbers of particles that are prohibitive for the classical PT algorithms. 3 GRW algorithm for two-dimensional diffusion in non-homogeneous and anisotropic media In the following we consider a two-dimensional diffusion process described by the diagonal diffusion tensor D x (x , y ) 0 0 D y (x , y ) and the diffusion equation ∂ t c = ∂ x ( D x ∂ x c) + ∂ y ( D y ∂ y c). The GRW solution of the diffusion equation with spatially dependent diffusion coefficients requires correction terms in Eq. (2). These corrections have to account for a drift term equal to (∂ x D x ,∂ y D y ) (see e.g. [6]). However, for the pur- pose of illustration, we consider here only the particular case of slowly variable diffusion coefficients, i.e. ∂ x D x ≈ 0 and ∂ y D y ≈ 0, for which the diffusion equation can be approxi- mated by ∂ t c = D x ∂ 2 x c + D y ∂ 2 y c. (7) The Eq. (7) has no drift terms and the diffusion coefficients vary in space. In this case the GRW relation (2) is replaced by n(i , j , k ) = δn(i , j | i , j , k ) + δn(i − d x , j | i , j , k ) + δn(i + d x , j | i , j , k ) + δn(i , j −d y | i , j , k ) + δn(i , j + d y | i , j , k ) (8) where n(i , j , k ) is the number of particles at the site (x i , y j ) = (i δ x , j δ y ) at the time k δt . Unlike the two-dimensional algorithm presented in Fig. 1, where the one-dimensional procedure is applied successively for lon- gitudinal and transversal directions, the procedure (8) moves the particles simultaneously along the principal directions of the diffusion tensor. Figure 3 illustrates this two-dimensional algorithm, in the general case of a non vanishing drift. The anisotropy is taken into account by two different para- meters d x and d y which describe the diffusive jumps along the coordinates axes. The spatial variation of the diffusion coefficients is described through the variable parameters r x = 2 D x δt (d x δ x ) 2 , r y = 2 D y δt (d y δ y ) 2 . (9) 0 0.1 0.2 0.3 0.4 x 1 (m) -0.1 0 0.1 x 2 (m) 1e+09 1e+10 n(x 1 ,x 2 ) Fig. 3 The two-dimensional GRW algorithm for variable diffusion coefficients After the global change of state of the lattice by the pro- cedure (8) applied at every site containing particles, the new numbers of particles at sites are obtained similarly to Eq. (3) by summation over two spatial indices. The average over an ensemble of simulations of the terms in Eq. (8) are related by δn(i ± d x | j , k ) = 1 2 r x (i , j ) n(i , j , k ), δn(i | j ± d y , k ) = 1 2 r y (i , j ) n(i , j , k ), δn(i , j | i , j , k ) =[1 − r x (i , j ) − r y (i , j )] n( j , k ). These relations can be used to show that for slowly variable diffusion coefficients the two-dimensional GRW algorithm approximates the finite difference scheme. The lattice steps δ x and δ y are chosen accordingly to the desired resolution of the concentration field. The time step δt is inferred from the condition r x + r y ≤ 1, which states that the numbers of diffusing particles are limited by the numbers of particles at the lattice sites. From Eq. (9) it follows δt ≤ 2 D max x δt (d x δ x ) 2 + 2 D max y δt (d y δ y ) 2 , (10) where D max x and D max y are the upper bounds of the diffusion coefficients D x (x , y ) and D y (x , y ). 4 Diffusion in skin The human skin varies in texture, structure and thickness and is made up of three main layers: namely epidermis, dermis and hypodermis. The epidermis is the uppermost layer of the skin, measures less than 1 mm in thickness and can be also subdivided into several layers. The Stratum Corneum (SC) is located on the outer surface of the skin and its thickness is in the range 10–100 μm. It is made of keratinized flat, roughly hexagonally shaped, partly overlapping cells embedded in a lipid matrix. Its main constituents are proteins, lipids and 123
Global random walk modelling of transport in complex systems 81 water. The dermis is the second layer of the skin and acts as a supportive layer for the epidermis and the hypodermis, the third layer of the skin, acts as a shock absorber and an energy reserve. The skin is a complex heterogeneous medium that plays an essential part in body protection and regulation: its most outer layer, the SC, acts as a protective barrier between the body and the environment and its main role to prevent fluid loss and to exclude environmental toxins, is fundamental to life. The excellent diffusional resistance of the SC makes the transdermal delivery of drugs pretty difficult. Nevertheless, there remains an important pharmaceutical need to reliably predict the topical and/or transdermal bioavailability of cuta- neously applied drugs. In this sense, to measure percutaneous penetration of the externally applied chemicals into human skin is very important for understanding how drugs can be delivered through the skin. In most cases a constant diffusion coefficient for trans- ported material adequately describes the diffusion process. However, in real biological systems, because of the hetero- geneity of the skin, the diffusion coefficient behaves like a random variable. This is why adequate analytical and nume- rical models able to describe the diffusion process in com- plex structures must be developed. An example of such a model was developed by Johnson et al.[7], who focused on measuring and modeling the molecular diffusion coefficient in order to describe the macroscopic SC permeation via the interkeratinocyte lipid. The two-dimensional algorithm introduced in Sect. 3 can be used as a new approach able to simulate the diffusion through the human skin. A two-dimensional geometry with the x axis parallel and the y axis perpendicular to the surface of the skin was used. A thin film consisting of N = 10 20 water molecules with ∆x = 10 mm was considered to be applied on the surface of the skin. Since the skin structure is stratified, the diffusion can be described by the two-dimensional equa- tion (7). To describe the diffusion in the horizontal direction the lattice dimension on the x axis was 3∆x . An acceptable resolution was obtained with δ x = 0.1 mm, (300 nodes per lattice). A simplified two layers structure of the skin was considered, with thickness y 1 = 0.1 mm and y 2 = 0.5 mm, respectively. With a resolution of δ y = 0.01 mm the lattice extended over 10 nodes in the first layer and over 50 nodes in the second. Because of the nonhomogeneous structure of the skin the diffusion coefficients show spatial fluctuations about the mean value. At every lattice site the coefficients D x and D y were generated as normal random variables with the mean value equal to half the maximal values D max x and D max y and variance equal to a fraction p < 1 of the corresponding maxi- mal values. In the simulations presented in the following we chose p = 0.1. Since the cells are rather flat, here we conside- red D max x = 2 D max y = 5.8810 −7 m 2 /s in the first layer. In the second layer an isotropic diffusion coefficient ten times lar- ger than in the superior layer was introduced. Between the layers a transition zone of thickness 3δ y was placed, where diffusion coefficients vary linearly [16]. Two different boundary conditions were used. At the sur- face of the skin, the molecules which have been jumped in the exterior were blocked on the boundary. At the inferior boun- dary we imposed a “transmission boundary condition” in the first order of approximation [15]. This last condition has the property to not disturb the diffusion front at the boundary. The time step was chosen using the condition (10). Some simulation results are presented in Figs. 4, 5, 6 and 7. Fig. 4 The distribution of molecules after t = 5 min Fig. 5 The distribution of molecules after t = 200 min 123
82 N. Suciu et al. Fig. 6 Depth averaged probability density as a function of longitudinal coordinate for fixed times Fig. 7 Time behavior of longitudinal and transversal flux of molecules through the skin 5 Transport in heterogeneous porous media Some applications to transport in groundwater show that GRW is suitable to be used in large-scale simulations and in Monte Carlo approaches on the predictability of the sto- chastic model [10, 13]. We considered two-dimensional divergence-free velocity fields with constant mean V= U = (U, 0), U = 1 m/day, given by the Darcy law for exponential correlated normal log-hydraulic conductivity with correlation length λ = 1 m, and variances σ 2 = 0.1. The random velocity was nume- rically generated, as a first-order approximation in σ , with the procedure based on the Kraichnan routine, already used in Refs. [1–3, 9, 13, 14]. In Ref. [13] it was shown that such velocity fields fulfil the requirements of the limit theorem (on the existence of the Gaussian up-scaling) of Kesten and Papanicolaou [5]. In every realization of the velocity field the simulation of an isotropic diffusion ( D ii = D = 0.01 m 2 /day, D ij = 0 for i = j , i = 1, 2) was conducted for dimensionless times t / U λ corresponding to 3,700 correlation lengths, using the reduced fluctuations GRW algorithm illustrated by Figs. 1 and 2. In each realization of the velocity field, N = 10 10 particles were initially located at the origin of the grid. It was shown that increasing the number of particles between 10 10 and 10 24 , the values of the effective coefficients, computed at successive moments over 5,000 time steps, are identical in the limit of the double precision [10]. To assess the reliability of the procedure, the simulations were repeated for increasing numbers of periods used in the Kraichnan routine of Np = 64, 640, and 6,400. Ensembles of 1,024 realizations of the transport simulations were computed for each of the three Np values considered. The lattice size of 10 7 nodes was larger than the maxi- mum extension of the plume in all simulations. We used a uniform lattice with constant δ x = 0.1 m, a time step δt = 0.5 day and diffusive jumps in Eq. (2) d = 2. These parameters are related with r = 0.25 accordingly to Eq. (5). To reduce the “overshooting” errors in particles methods one impose that the mean displacement in a time step does not overtake a given threshold [8]. For the same overshooting, the error of GRW simulations is mainly influenced by the discretization of the velocity, as described by the parame- ter U δt /δ x . It was found that increasing this parameter from five (as in the present simulations) to ten results in differences smaller than 2% [10]. For U δt /δ x = 5 and the parameters above, the large-scale GRW simulations estimate the obser- vables of the transport in groundwater with errors which are one or two orders of magnitude smaller than the correspon- ding first-order approximations in σ 2 [11, 12]. The parallelization was trivial since the computation of the Kraichnan velocity field and the simulation of the transport in a given realization were performed on a single processor of an IBM Regatta p690+ computer (the JUMP system at Research Center Jülich, Germany). The costs were considerable, for instance of about 4,500 cpu hours for the ensemble of 1,024 realizations in the case Np = 64. Nevertheless, a comparison with other simulations (see e.g. [9]) shows that the costs are 4.4millions times smaller than they could have been for the PT method. The effective coefficients of transport in groundwater were estimated through the rate of increase of the second cen- tral moment of the plume D eff ii = σ 2 ii /(2t ), which in the large time limit should coincide with the coefficients of the up-scaled Gaussian process. For Np = 6,400 periodic modes in Kraichnan routine and for an initial solute plume with transverse dimension of 100λ, the root-mean-square deviation of the effective coefficients from the up-scaled coefficients, computed at distances larger than 200λ, was indeed found to be one order of magnitude smaller than the local coefficient D [13, Fig. 9], as predicted by theory [5, 13]. 123
Global random walk modelling of transport in complex systems 83 0 20 40 60 80 100 120 1 10 100 1000 (σ D eff 11 / <D eff 11 >) 2 (%) Ut/λ Np=64 Np=640 Np=6400 Fig. 8 Sample to sample fluctuations of the longitudinal effective coefficient as function of dimensionless time, number of realizations and number of periods Np, for 1,024 realizations of the transport simulation In the following we present some results concerning the predictive value of the stochastic approach for the contami- nation in actual aquifer systems, in the case of point instan- taneous sources. Following the approach of previous works [1, 13] we computed the sample to sample fluctuations of the effective coefficients σ D eff ii D eff ii 2 , where σ 2 D eff ii = D eff ii 2 − D eff ii 2 . The results presented in Fig. 8 indicate that at large times the fluctuations of the longitudinal effective coefficient tend to zero, i.e. the coefficient is “self-averaging”, only when the number of periods in the Kraichnan routine, Np, is of the order of 1,000s or larger. The self-averaging of the longitudi- nal coefficients has been also investigated by approximations of the transport equations up to the first order in velocity fluc- tuations, which were previously evaluated by comparisons with GRW simulations [14]. It was found that, in order to account for the self-averaging behavior of transport in Gaus- sian velocity fields, the number of periodic modes in the Kraichnan procedure must be at least as large as the total dimensionless time Ut /λ of the simulation [2]. The transverse coefficient, shown in Fig. 9, has much smaller fluctuations than the longitudinal coefficient and is self-averaging even for small Np. One can see that Np = 640 periodic modes yields practically the same result as Np = 6,400. The normalized concentration in given realizations c(x 1 , x 2 , t ) was computed, similarly with Eq. (6), as the num- ber of particles in a domain 1λ × 1λ centered at (x 1 , x 2 ) divided by the total number of particles N . Further, the point concentrations were averaged over a domain of dimensions 1λ × 160λ, oriented across the mean flow. This domain is larger than the transverse dimension of the plume and, 0 5 10 15 20 25 30 1 10 100 1000 (σ D eff 22 / <D eff 22 >) 2 (%) Ut/λ Np=64 Np=640 Np=6400 Fig. 9 Sample to sample fluctuations of the transverse effective coef- ficient as function of dimensionless time, number of realizations and number of periods Np, for 1,024 realizations of the transport simula- tion 0 5 10 15 20 25 30 35 1 10 100 1000 CV = σ C / <C> (%) Ut/λ Np=64 Np=640 Np=6400 Fig. 10 Coefficient of variation of the cross-section averaged concen- tration at the plume center of mass consequently the concentration field can be described by an one-dimensional problem for the transverse averaged concentration C (x 1 , t ). This quantity corresponds to the ideal sampling across reference planes aimed at in field experi- ments [17]. Further, We stored the mean concentration at the plume center of mass, C (x 1 , t ), and we computed the mean C (x 1 , t )over the ensemble of 1,024 realizations of the transport. Finally we computed the variance at the center of mass σ 2 c (x 1 , t ) = C 2 −C 2 (x 1 , t ) and the coefficient of variation CV = σ c / C . In Fig. 10 one remarks that for Np = 6,400 and after more than 2,000 dimensionless times the coefficient of variation becomes smaller than 10% and shows a monotonous decrease. The self-averaging behavior of the effective coefficients and cross-section concentration, shown in Figs. 8, 9 and 123
84 N. Suciu et al. 10, indicates that single realizations of the transport process behave in the large time limit as the ensemble averaged pro- cess. Since the latter tends to an up-scaled Gaussian pro- cess, as shown by both theoretical and numerical results, one concludes that the process is asymptotically “ergodic” in a large sense [13]. Concerning the predictive power of the sto- chastic up-scaled model, the present results complete and reinforce the conclusions drown by previous works [9, 1]: Even under rather ideal conditions (σ 2 = 0.1), in the case of small initial contaminant plumes the predictions of the up- scaled model are reliable for real groundwater systems only after travelled distances of 1,000s of heterogeneity scales. 6 Conclusions The simulated concentration field capture the effects of non- homogeneity and anisotropy of the skin structure. These first results indicate that the simulations can be used in inves- tigations on diffusion and other microscopic processes in biological structures. Comparisons between measured quan- tities (for instance fluxes of molecules) and GRW simulations could be helpful for the elaboration of more complex codes, able to account for less idealized skin structures. The stochastic simulations of the advection dominated transport in groundwater, based on GRW procedure, are in good agreement with the exact mathematical result on the existence of the up-scaled Gaussian process. Therefore, we expect that our method could be successfully used to analyze more realistic situations, with larger variance of the log-hydraulic conductivity, for which there are no exact mathematical results. The accuracy of the simulated concen- tration field and the possibility to describe large-scale beha- vior of the transport process make the GRW method appro- priate for simulations of field experiments, such as those pre- sented in [17]. The relatively small computational costs allow the implementation of GRW algorithm in inverse modeling procedures for the assessment of the aquifer’s properties. Though it can be considered as a particular cellular auto- maton, the GRW procedure differs from most of lattice gas models because no exclusion principle is necessary. As a consequence, the number of particles used in simulations can be as large as the number of molecules involved in a specific transport problem. This can be a great advantage in simulations for reactive and radioactive decay processes. The method presented in this paper does not completely remove the overshooting errors associated with the spatial variation of the velocity field. The overshooting can be com- pletely removed with a biased GRW code. The principle of the biased GRW consists in the simulation of the advec- tive displacement by a bias in the random walk jumps [11]. This procedure remains “global”, unlike similar approaches [4] where the updated lattice configurations are obtained by sequential random decisions for the jumps of the particles at a lattice site. Since it introduces new constraints, the bia- sed procedure requires larger computational domains. It is yet possible to use the biased GRW algorithm in large-scale problems owing to its cellular automaton property which permits domain-decomposition on parallel machines. When domain-decomposition does not increase the computational efficiency, as in case of advection-diffusion simulations for groundwater, biased GRW test simulations can be useful for the calibration and the evaluation of the less expensive un- biased GRW algorithm [14, 12]. Acknowledgments The work presented in this paper was supported by the Deutsche Forschungsgemeinschaft grant SU 415/1-2, the Fun- damental Programme No. 1 of the Romanian Academy, and the Roma- nian Ministry of Education and Research grants CEEX-M1-C2-2532 and 01-8-CPD-042- BIOTECH Program. The simulations of transport in groundwater were performed in the frame of the computing project JICG41 at ICG IV: Agrosphere, Research center Jülich, Germany. References 1. Eberhard, J.: Approximations for transport parameters and self- averaging properties for point-like injections in heterogeneous media. J. Phys. A Math. Gen. 37, 2549–2571 (2004) 2. Eberhard, J., Suciu N., Vamo¸ s: On the self-averaging of dispersion for transport in quasi-periodic random media. J. Phys. A: Math. Gen. 40, 597–610, doi:10.1088/1751-8113/40/4/002 (2007) 3. Jaekel, U., Vereecken, H.: Renormalization group analysis of macrodispersion in a directed random flow. Water Resour. Res. 33, 2287–2299 (1997) 4. Karapiperis, T., Blankleider, B.: Cellular automaton model of reaction-transport processes. Physica D 78, 30–64 (1991) 5. Kesten, H., Papanicolaou, G.C.: A limit theorem for turbulent dif- fusion. Commun. Math. Phys. 65, 97–128 (1979) 6. Kinzelbach, W., Uffink, G.: The random walk method and exten- sions in groundwater modelling. In: Bear, J., Corapcioglu, M.Y. (eds.) Transport Processes in Porous Media, pp. 761–787. Kluwer, Norwell (1991) 7. Johnson, M.E., Blankschtein, D., Langer, R.: Evaluation of solute permeation through the stratum corneum: lateral bilayer diffusion as the primary transport mechanism. J. Pharm. Sci. 86(10), 1162– 1172 (1997) 8. Roth, K., Hammel, K.: Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with steady state flow. Water Resour. Res. 32(6), 1653–1663 (1996) 9. Schwarze, H., Jaekel, U., Vereecken, H.: Estimation of macrodis- persivity by different approximation methods for flow and transport in randomly heterogeneous media. Transp. Porous Media 43, 265– 287 (2001) 10. Suciu, N., Vamo¸ s, C., Vanderborght, J., Hardelauf, H., Vereecken, H.: Numerical modeling of large scale transport of contminant solutes using the global random walk algorithm. Monte Carlo Methods Appl. 10(2), 153–177 (2004) 11. Suciu, N., Vamo¸ s, C., Knabner, P., Rüde, U.: Biased global ran- dom walk, a cellular automaton for diffusion. In: Hüsemann, F., Kowarschik, M., Rüde, U. (eds.) Simulationstechnique, 18th Sym- posium in Erlangen, September 2005., pp. 562–567. SCS Publi- shing House e. V., Erlangen (2005) 123
Global random walk modelling of transport in complex systems 85 12. Suciu, N., Vamo¸ s, C.: Evaluation of overshooting errors in particle methods for diffusion by biased global random walk. Rev. Anal. Num. Th. Approx. 35, 119–126 (2006) 13. Suciu, N., Vamo¸ s, C., Vanderborght, J., Hardelauf, H., Vereecken, H.: Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res. 42, W04409. doi:10.1029/2005WR004546 (2006) 14. Suciu, N., Vamo¸ s, C., Eberhard, J.: Evaluation of the first- order approximations for transport in heterogeneous media. Water Resour. Res. 42, W11504. doi:10.1029/2005WR004714 (2006) 15. Vamo¸ s, C., Suciu, N., Vereecken, H.: Generalized random walk algorithm for the numerical modeling of complex diffusion pro- cesses. J. Comput. Phys. 186(2), 527–544 (2003) 16. Vamo¸ s, C., Suciu, N., Turcu, I., Pop, C.V.L., Ciortea, L. I.: Programme BIOTECH, Project No. 01-8-CPD-042/19.10.2001. Research report (2004, in Romanian) 17. Vereecken, H., Döring, U., Hardelauf, H., Jaekel, U., Hashagen, U., Neuendorf, O., Schwarze, H., Seidemann, R.: Analysis of solute transport in a heterogeneous aquifer: the Krauthausen field expe- riment. J. Contam. Hydol. 45, 329–358 (2000) 18. Wilke, J., Pohl, T., Kowarschik, M., Rüde, U.: Cache performance optimizations for parallel lattice Boltzmann codes. In: Proceedings of the EuroPar-03 Conference. Lecture Notes in Computer Science, vol. 2790, pp. 441–450. Springer, Heidelberg (2003) 123
2007, 2009

Related Posts