Global random walk solutions to PDF evolution equations

Nicolae Suciu
(preprint), html, paper

Abstract

Authors

Nicolae Suciu
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy
Friedrich-Alexander University of Erlangen-Nuremberg
Mathematics Department

Călin Vamoș
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

Sabine Attinger
UFZ-Helmholtz Center for Environmental Research
Division of Computational Environmental Systems

Peter Knabner
Friedrich-Alexander University of Erlangen-Nuremberg
Mathematics Department

Keywords

Paper coordinates

N. Suciu, C. Vamoş, S. Attinger, P. Knabner, Global random walk solutions to PDF evolution equations, Proceedings of XIX International Conference Computational Methods in Water Resources CMWR 2012, University of Illinois at Urbana-Champaign, June 17-21, 2012

PDF

Suciu.Nicolae.pdf

About this paper

Journal

 Proceedings of XIX International Conference Computational Methods in Water Resources CMWR 2012

Publisher Name
DOI
Print ISSN
Online ISSN

google scholar link

[1] S. B. Pope. PDF methods for turbulent reactive flows. Prog. Energy Combust. Sci., 11(2), 119–192, (1885).
[2] D. C. Haworth. Progress in probability density function methods for turbulent reacting flows. Prog. Energy Combust. Sci., 36, 168–259, doi:10.1016/j.pecs.2009.09.003, (2010).
[3] S. Viswanathan, H. Wang, and S. B. Pope. S. Numerical implementation of mixing and molecular transport in LES/PDF studies of turbulent reacting flows. J. Comput. Phys., 230, 6916–6957, doi:10.1016/j.jcp.2011.05.020, (2011).
[4] C. Vamo¸s, N. Suciu and H. Vereecken. Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys., 186(2), 527
544, doi:10.1016/S0021-9991(03)00073-1, (2003).
[5] G. Dagan. Flow and Transport in Porous Formations, Springer, Berlin (1989).
[6] N. Suciu. Spatially inhomogeneous transition probabilities as memory effects for dif fusion in statistically homogeneous random velocity fields. Phys. Rev. E, 81, 056301, (2010).
[7] P. E. Kloeden, and E. Platen. Numerical Solutions of Stochastic Differential Equa tions, Springer, Berlin, (1999).
[8] D. W. Meyer, P. Jenny and H. A. Tchelepi. A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media. Water Resour. Res., 46,
W12522, (2010). doi:10.1029/2010WR009450.
[9] S. Attinger, M. Dentz, H. Kinzelbach and W. Kinzelbach. Temporal behavior of asolute cloud in a chemically heterogeneous porous medium. J. Fluid Mech., 386,
77-104, (1999).
[10] N. Suciu, C. Vamos, J. Vanderborght, H. Hardelauf and H. Vereecken. Numericalinvestigations on ergodicity of solute transport in heterogeneous aquifers. Water Re
sour. Res., 42, W04409, doi:10.1029/2005WR004546, (2006) .
[11] N. Suciu, C. Vamos, F. A. Radu, H. Vereecken and P. Knabner. Persisten memoryof diffusing particles. Phys. Rev. E, 80, 061134, doi:10.1103/PhysRevE.80.061134, (2009).
[12] S. B. Pope. Lagrangian PDF methods for turbulent flows. Annu. Rev. Fluid Mech., 26, 23-63, (1994)

Paper (preprint) in HTML form

Suciu.Nicolae

GLOBAL RANDOM WALK SOLUTIONS TO PDF EVOLUTION EQUATIONS

Nicolae Suciu ∗ , † ∗ , † ^(**,†){ }^{*, \dagger}∗,†, Călin Vamoş † † ^(†){ }^{\dagger}†, Sabine Attinger † † † † ^(††){ }^{\dagger \dagger}†† and Peter Knabner**Friedrich-Alexander University of Erlangen-NurembergMathematics DepartmentMartensstrasse 3, 91058 Erlangen, Germanye-mail: {suciu,knabner}@am.uni-erlangen.de, web page:http://www.mso.math.uni-erlangen.de/ † † ^(†){ }^{\dagger}† Romanian AcademyTiberiu Popoviciu Institute of Numerical AnalysisFântanele 57, 400320 Cluj-Napoca, Romaniae-mail: cvamos@ictp.acad.ro, web page: http://www.ictp.acad.ro/ † † † † ^(††){ }^{\dagger \dagger}†† UFZ-Helmholtz Center for Environmental ResearchDivision of Computational Environmental SystemsPermoserstrasse 15, 04318 Leipzig, Germanye-mail: sabine.attinger@ufz.de, web page: http://www.ufz.de/

Key words: Advanced numerical methods, Random walk, PDF evolution equations
Summary. We introduce a global random walk method to solve modeled equations for the evolution of the probability density function of the random concentration in passive transport through heterogeneous aquifers. The algorithm consists of a superposition of many weak Euler schemes for the model Itô equations. The new method avoids the limitations concerning the number of particles in Lagrangian approaches, completely removes the numerical diffusion and speeds up the computation by orders of magnitude.

1 INTRODUCTION

Modeled Probability Density Function (PDF) evolution equations are in most cases solved by methods based on a representation of the PDF by an ensemble of notional particles, each carrying information on scalar values (e.g. species concentrations) in the physical space. The particles move in physical space along stochastic-Lagrangian trajectories governed by Itô equations, with drift coefficients given by the local values of the resolved velocity field and diffusion coefficients obtained by upscaling procedures 1 1 ^(1){ }^{1}1. The evolution of their scalar composition is described by various mixing models, a general one consisting of a system of Itô equations solving for trajectories in the composition space 2 2 ^(2){ }^{2}2.
The Lagrangian PDF approach is much more efficient than classical methods for partial differential equations, mainly in case of high dimensional multi-species reactive transport problems 2 2 ^(2){ }^{2}2. In spite of their efficiency, Lagrangian approaches suffer from two severe
limitations. Since the particle trajectories are constructed sequentially, the demanded computing resources increase linearly with the number of particles 3 3 ^(3){ }^{3}3. Moreover, the need to gather particles at the center of computational cells to perform the mixing step and to estimate statistical parameters as well as the interpolation of various terms to particles positions inevitably produce numerical diffusion in either particle-mesh or grid-free particles methods 1 , 3 1 , 3 ^(1,3){ }^{1,3}1,3.
Our idea to overcome these limitations is to construct weak solutions, which give the probability distributions without solving for individual trajectories. The sequential Lagrangian approach is essentially looking for strong solutions to the Itô equations consisting of path-wise unique trajectories in physical and concentration spaces of the notional particles describing the evolution of the PDF. Instead of Gaussian noises, on which the construction of strong solutions is based, we use surrogate-noise random variables of constant amplitude, determined by the diffusion coefficients, and unbiased random orientation. With this choice, the movement of the ensemble of notional particles is described by a superposition of weak Euler schemes on a regular lattice and the particles from a lattice site can be spread simultaneously according to a binomial distribution. This is the principle of the Global Random Walk algorithm (GRW), which is practically insensitive to the increase of the number of particles 4 4 ^(4){ }^{4}4. Because mixing in concentration space takes place only at lattice sites, the GRW approach also avoids the occurrence of the artificial diffusion during the mixing step.
The GRW solution to PDF evolution equations will be illustrated for passive transport in aquifers. The method will be validated by comparisons with Monte Carlo simulations of passive advection-dispersion transport in saturated aquifers with randomly distributed hydraulic conductivity.

2 MODELED PDF EQUATIONS

A very simple example of modeled PDF equations for transport in groundwater, which illustrates the main features of the PDF methods, can be constructed for the "macrodispersion model" of Dagan 5 5 ^(5){ }^{5}5. Starting with the stochastic partial differential equation
(1) ∂ t c + V ∇ c = D ∇ 2 c (1) ∂ t c + V ∇ c = D ∇ 2 c {:(1)del_(t)c+Vgrad c=Dgrad^(2)c:}\begin{equation*} \partial_{t} c+\mathbf{V} \nabla c=D \nabla^{2} c \tag{1} \end{equation*}(1)∂tc+V∇c=D∇2c
with random coefficients V V V\mathbf{V}V (because of the randomness of the hydraulic conductivity field K K KKK ), the macrodispersion model is derived in from of a deterministic advection-diffusion equation,
(2) ∂ t c ~ + V ∇ c ~ = D i j ∂ x i ∂ x j c ~ (2) ∂ t c ~ + V ∇ c ~ = D i j ∂ x i ∂ x j c ~ {:(2)del_(t) tilde(c)+Vgrad tilde(c)=D_(ij)del_(x_(i))del_(x_(j)) tilde(c):}\begin{equation*} \partial_{t} \tilde{c}+\mathcal{V} \nabla \tilde{c}=\mathcal{D}_{i j} \partial_{x_{i}} \partial_{x_{j}} \tilde{c} \tag{2} \end{equation*}(2)∂tc~+V∇c~=Dij∂xi∂xjc~
If the ln ⁡ K ln ⁡ K ln K\ln Kln⁡K field is a statistically homogeneous random field with small variance, the consistent first-order approximation in the variance of ln ⁡ K ln ⁡ K ln K\ln Kln⁡K consists of sampling the random velocity on the mean-flow trajectory 6 6 ^(6){ }^{6}6 and the coefficients of (2) become (at the most) time functions
(3) V = ⟨ V ⟩ and D i j ( t ) = D + ∫ 0 t Cov V , i j ( t ′ ) d t ′ (3) V = ⟨ V ⟩  and  D i j ( t ) = D + ∫ 0 t   Cov V , i j t ′ d t ′ {:(3)V=(:V:)quad" and "quadD_(ij)(t)=D+int_(0)^(t)Cov_(V,ij)(t^('))dt^('):}\begin{equation*} \mathcal{V}=\langle\mathbf{V}\rangle \quad \text { and } \quad \mathcal{D}_{i j}(t)=D+\int_{0}^{t} \operatorname{Cov}_{V, i j}\left(t^{\prime}\right) d t^{\prime} \tag{3} \end{equation*}(3)V=⟨V⟩ and Dij(t)=D+∫0tCovV,ij(t′)dt′
Equation (2) describes the evolution of the ensemble mean concentration, c ~ = ⟨ c ⟩ c ~ = ⟨ c ⟩ tilde(c)=(:c:)\tilde{c}=\langle c\ranglec~=⟨c⟩, at finite times as well as the asymptotic long-time behavior, provided that the integral of the Lagrangian velocity covariance in the second equation (3) is finite 6 6 ^(6){ }^{6}6. A more complete solution of the stochastic partial differential equation (1) is the one-point concentration PDF , p ( c ; x , t ) PDF , p ( c ; x , t ) PDF,p(c;x,t)\mathrm{PDF}, p(c ; \mathbf{x}, t)PDF,p(c;x,t). Its evolution equation can be inferred as follows. First, by Itô - FokkerPlanck equivalence 7 7 ^(7){ }^{7}7, one associates to (2) "fluid particles" of coordinates X i , i = 1 , 2 , 3 X i , i = 1 , 2 , 3 X_(i),i=1,2,3X_{i}, i=1,2,3Xi,i=1,2,3, and (stochastic) characteristics, or (stochastic) Lagrangian trajectories
(4) d X i ( t ) = V i ( t ) d t + d W ~ i ( t ) (4) d X i ( t ) = V i ( t ) d t + d W ~ i ( t ) {:(4)dX_(i)(t)=V_(i)(t)dt+d tilde(W)_(i)(t):}\begin{equation*} d X_{i}(t)=\mathcal{V}_{i}(t) d t+d \tilde{W}_{i}(t) \tag{4} \end{equation*}(4)dXi(t)=Vi(t)dt+dW~i(t)
where W ~ i W ~ i tilde(W)_(i)\tilde{W}_{i}W~i is a diffusion process of mean zero and covariance 2 ∫ 0 t D i j ( t ′ ) d t ′ 2 ∫ 0 t   D i j t ′ d t ′ 2int_(0)^(t)D_(ij)(t^('))dt^(')2 \int_{0}^{t} \mathcal{D}_{i j}\left(t^{\prime}\right) d t^{\prime}2∫0tDij(t′)dt′. The random concentration C ( x , t ) C ( x , t ) C(x,t)C(\mathbf{x}, t)C(x,t) at a fixed space point x x x\mathbf{x}x also can be modeled by an Itô equation 1 , 2 , 8 1 , 2 , 8 ^(1,2,8){ }^{1,2,8}1,2,8
(5) d C ( t ) = A ( t ) d t + B ( t ) d W ( t ) (5) d C ( t ) = A ( t ) d t + B ( t ) d W ( t ) {:(5)dC(t)=A(t)dt+B(t)dW(t):}\begin{equation*} d C(t)=A(t) d t+B(t) d W(t) \tag{5} \end{equation*}(5)dC(t)=A(t)dt+B(t)dW(t)
where W ( t ) W ( t ) W(t)W(t)W(t) is a standard Wiener process. The coefficients A A AAA and B B BBB are subject of modeling, perhaps the most difficult task in modeled PDF methods 1 , 2 1 , 2 ^(1,2){ }^{1,2}1,2.
By Itô - Fokker-Planck equivalence again, one associates to the positions ( X i X i X_(i)X_{i}Xi ) and concentrations ( C ) ( C ) (C)(C)(C) Itô equations (4-5) a Fokker-Planck equation describing the evolution of the joint PDF for concentration and position p ( c , x , t ) p ( c , x , t ) p(c,x,t)p(c, \mathbf{x}, t)p(c,x,t),
(6) ∂ t p + V ∇ p + ∂ c ( A p ) = D i j ∂ x i ∂ x j p + ∂ c 2 ( B 2 p ) (6) ∂ t p + V ∇ p + ∂ c ( A p ) = D i j ∂ x i ∂ x j p + ∂ c 2 B 2 p {:(6)del_(t)p+Vgrad p+del_(c)(Ap)=D_(ij)del_(x_(i))del_(x_(j))p+del_(c)^(2)(B^(2)p):}\begin{equation*} \partial_{t} p+\mathcal{V} \nabla p+\partial_{c}(A p)=\mathcal{D}_{i j} \partial_{x_{i}} \partial_{x_{j}} p+\partial_{c}^{2}\left(B^{2} p\right) \tag{6} \end{equation*}(6)∂tp+V∇p+∂c(Ap)=Dij∂xi∂xjp+∂c2(B2p)
For D = 0 D = 0 D=0D=0D=0, there is no mixing and equation (4) with coefficients approximated to the first-order by (3) describes the movement of the erratic center of mass 9 9 ^(9){ }^{9}9. Since for infinite Peclet number equation (5) reduces to d C ( t ) = 0 d C ( t ) = 0 dC(t)=0d C(t)=0dC(t)=0 the initial concentration is transported without dilution along the characteristics 8 8 ^(8){ }^{8}8 and (6) becomes identical to equation (2) for the mean. The consistency requirement that marginal probabilities are obtained by integrating the joint PDF 10 PDF 10 PDF^(10)\mathrm{PDF}^{10}PDF10 is also fulfilled for D ≥ 0 D ≥ 0 D >= 0D \geq 0D≥0 when (6) is integrated over the concentration space.
For reactive transport and vector concentrations c = ( c 1 , ⋯ , c S ) c = c 1 , ⋯ , c S c=(c_(1),cdots,c_(S))\mathbf{c}=\left(c_{1}, \cdots, c_{S}\right)c=(c1,⋯,cS), where S S SSS is the number of molecular species, the equation (1) includes source (reaction) terms R s ( c ) , s = 1 , ⋯ , S R s ( c ) , s = 1 , ⋯ , S R_(s)(c),s=1,cdots,SR_{s}(\mathbf{c}), s= 1, \cdots, SRs(c),s=1,⋯,S, which also occur as supplementary drift terms in the system of Itô equations 1 , 2 1 , 2 ^(1,2){ }^{1,2}1,2. Unlike in turbulence applications, immobile species c ~ c ~ tilde(c)\tilde{\mathbf{c}}c~ may also be present in groundwater systems. The PDF of c ~ c ~ tilde(c)\tilde{\mathbf{c}}c~ is no longer transported by (4) and only depends on space through drift terms R ( c , c ~ ) R ( c , c ~ ) R(c, tilde(c))R(\mathbf{c}, \tilde{\mathbf{c}})R(c,c~) in (5).

3 GRW ALGORITHM

We consider particles trajectories governed by the one-dimensional version of the Itô equation (4) with constant diffusion coefficient D D DDD. The one-dimensional GRW algorithm describes the evolution of n ( j , k ) n ( j , k ) n(j,k)n(j, k)n(j,k) particles over the sites i i iii of a regular lattice, at discrete times indexed by k k kkk, by spreading the particles according to the rule
(7) n ( j , k ) = δ n ( j + v j , j , k ) + δ n ( j + v j − d , j , k ) + δ n ( j + v j + d , j , k ) (7) n ( j , k ) = δ n j + v j , j , k + δ n j + v j − d , j , k + δ n j + v j + d , j , k {:(7)n(j","k)=delta n(j+v_(j),j,k)+delta n(j+v_(j)-d,j,k)+delta n(j+v_(j)+d,j,k):}\begin{equation*} n(j, k)=\delta n\left(j+v_{j}, j, k\right)+\delta n\left(j+v_{j}-d, j, k\right)+\delta n\left(j+v_{j}+d, j, k\right) \tag{7} \end{equation*}(7)n(j,k)=δn(j+vj,j,k)+δn(j+vj−d,j,k)+δn(j+vj+d,j,k)
where v j = [ V j δ t / δ x ] , [ ⋅ ] v j = V j δ t / δ x , [ ⋅ ] v_(j)=[V_(j)delta t//delta x],[*]v_{j}=\left[V_{j} \delta t / \delta x\right],[\cdot]vj=[Vjδt/δx],[⋅] denoting the integer part, are discrete displacements produced by the velocity field and d d ddd describes the diffusive jumps. The quantities δ n ( j + v j ± d , j , k ) δ n j + v j ± d , j , k delta n(j+v_(j)+-d,j,k)\delta n\left(j+v_{j} \pm d, j, k\right)δn(j+vj±d,j,k) in (7) are Bernoulli random variables associated to the number of particles jumping to the left and to the right from the advected position j + v j j + v j j+v_(j)j+v_{j}j+vj and δ n ( j + v j , j , k ) δ n j + v j , j , k delta n(j+v_(j),j,k)\delta n\left(j+v_{j}, j, k\right)δn(j+vj,j,k) is the number of particles that are allowed to stay at the lattice site reached after an advective displacement during the current time step. The distribution of the particles at the next time ( k + 1 k + 1 k+1k+1k+1 ) is given by
(8) n ( i , k + 1 ) = δ n ( i , i , k ) + ∑ j ≠ i δ n ( i , j , k ) (8) n ( i , k + 1 ) = δ n ( i , i , k ) + ∑ j ≠ i   δ n ( i , j , k ) {:(8)n(i","k+1)=delta n(i","i","k)+sum_(j!=i)delta n(i","j","k):}\begin{equation*} n(i, k+1)=\delta n(i, i, k)+\sum_{j \neq i} \delta n(i, j, k) \tag{8} \end{equation*}(8)n(i,k+1)=δn(i,i,k)+∑j≠iδn(i,j,k)
The averages over repeated simulations of the number of particles undergoing diffusive jumps and of the number of particles remaining at the same node after the displacement v j v j v_(j)v_{j}vj are given by the relations
(9) δ n ( j + v j ± d , j , k ) ― = 1 2 r n ( j , k ) ― , δ n ( j , j + v j , k ) ― = ( 1 − r ) n ( j , k ) ― (9) δ n j + v j ± d , j , k ¯ = 1 2 r n ( j , k ) ¯ , δ n j , j + v j , k ¯ = ( 1 − r ) n ( j , k ) ¯ {:(9) bar(delta n(j+v_(j)+-d,j,k))=(1)/(2)r bar(n(j,k))","quad bar(delta n(j,j+v_(j),k))=(1-r) bar(n(j,k)):}\begin{equation*} \overline{\delta n\left(j+v_{j} \pm d, j, k\right)}=\frac{1}{2} r \overline{n(j, k)}, \quad \overline{\delta n\left(j, j+v_{j}, k\right)}=(1-r) \overline{n(j, k)} \tag{9} \end{equation*}(9)δn(j+vj±d,j,k)―=12rn(j,k)―,δn(j,j+vj,k)―=(1−r)n(j,k)―
where 0 ≤ r ≤ 1 0 ≤ r ≤ 1 0 <= r <= 10 \leq r \leq 10≤r≤1 is a rational number. The space step δ x δ x delta x\delta xδx and the time step δ t δ t delta t\delta tδt are related to the diffusion coefficient D D DDD by
(10) D = r ( d δ x ) 2 2 δ t (10) D = r ( d δ x ) 2 2 δ t {:(10)D=r((d delta x)^(2))/(2delta t):}\begin{equation*} D=r \frac{(d \delta x)^{2}}{2 \delta t} \tag{10} \end{equation*}(10)D=r(dδx)22δt
Using (10), for given δ x δ x delta x\delta xδx and δ t δ t delta t\delta tδt and adjusting the parameters d d ddd and r r rrr, any possible value of the diffusion coefficient D D DDD can be approximated 4 4 ^(4){ }^{4}4.
According to (9) r r rrr is the particles' jump probability. Each particle can either stay at the position reached after the advection step with probability 1 − r 1 − r 1-r1-r1−r or can execute unbiased jumps of length d δ x d δ x d delta xd \delta xdδx with probability r r rrr. The random variable defined by the two probabilities is a surrogate Wiener process which fulfills the conditions ensuring the weak convergence of order one of the weak Euler scheme solving the Itô equation for individual particles 7 7 ^(7){ }^{7}7. It follows that the GRW algorithm is a superposition of weak Euler schemes. We also note that r ( d δ x ) 2 r ( d δ x ) 2 r(d delta x)^(2)r(d \delta x)^{2}r(dδx)2 is the local variance over the time δ t δ t delta t\delta tδt of the diffusion process starting at a given lattice site. Thus, the relation (10) is the general definition of the diffusion coefficient 7 7 ^(7){ }^{7}7 discretized on the GRW lattice. This ensures that the GRW algorithm is unconditionally free of numerical diffusion.
The simplest way to construct the random variables δ n δ n delta n\delta nδn is to define
(11) δ n ( j + v j − d , j , k ) = { n / 2 if n is even [ n / 2 ] + θ if n is odd (11) δ n j + v j − d , j , k = n / 2  if  n  is even  [ n / 2 ] + θ  if  n  is odd  {:(11)delta n(j+v_(j)-d,j,k)={[n//2," if "n" is even "],[[n//2]+theta," if "n" is odd "]:}:}\delta n\left(j+v_{j}-d, j, k\right)=\left\{\begin{array}{cl} n / 2 & \text { if } n \text { is even } \tag{11}\\ {[n / 2]+\theta} & \text { if } n \text { is odd } \end{array}\right.(11)δn(j+vj−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 ] n = n ( j , k ) − δ n j + v j , j , k , [ n / 2 ] n=n(j,k)-delta n(j+v_(j),j,k),[n//2]n=n(j, k)-\delta n\left(j+v_{j}, j, k\right),[n / 2]n=n(j,k)−δn(j+vj,j,k),[n/2] is the integer part of n / 2 n / 2 n//2n / 2n/2 and θ θ theta\thetaθ is a random variable taking the values 0 and 1 with probability 1 / 2 1 / 2 1//21 / 21/2. Further, the number of particles jumping in the opposite direction, δ n ( j + v j + d , j , k ) δ n j + v j + d , j , k delta n(j+v_(j)+d,j,k)\delta n\left(j+v_{j}+d, j, k\right)δn(j+vj+d,j,k), is determined by (7). We designed an efficient implementation of (11) which avoids using random number generators and
considerably increases the computational efficiency. It is obtained by redistributing the rests R j = n ( j , k ) − [ ( 1 − r ) n ( j , k ) ] R j = n ( j , k ) − [ ( 1 − r ) n ( j , k ) ] R_(j)=n(j,k)-[(1-r)n(j,k)]R_{j}=n(j, k)-[(1-r) n(j, k)]Rj=n(j,k)−[(1−r)n(j,k)], appearing in the computation of δ n ( j + v j , j , k ) = [ ( 1 − r ) n ( j , k ) ] δ n j + v j , j , k = [ ( 1 − r ) n ( j , k ) ] delta n(j+v_(j),j,k)=[(1-r)n(j,k)]\delta n\left(j+v_{j}, j, k\right)= [(1-r) n(j, k)]δn(j+vj,j,k)=[(1−r)n(j,k)], and the rests R j + v j = n / 2 − [ n / 2 ] R j + v j = n / 2 − [ n / 2 ] R_(j+v_(j))=n//2-[n//2]R_{j+v_{j}}=n / 2-[n / 2]Rj+vj=n/2−[n/2] of the division by 2 of the number of jumping particles to calculate δ n ( j + v j − d , j , k ) = [ n / 2 + R j + R j + v j ] δ n j + v j − d , j , k = n / 2 + R j + R j + v j delta n(j+v_(j)-d,j,k)=[n//2+R_(j)+R_(j+v_(j))]\delta n\left(j+v_{j}-d, j, k\right)=\left[n / 2+R_{j}+R_{j+v_{j}}\right]δn(j+vj−d,j,k)=[n/2+Rj+Rj+vj]. To ensure the strict conservation of the number of particles, besides the condition requiring that ( 1 − r ) N ( 1 − r ) N (1-r)N(1-r) N(1−r)N be an integer, related to the use of the parameter r r rrr, the total number of particles N N NNN should also be a power of 2 . Nevertheless, we have checked that if N ≥ 10 8 N ≥ 10 8 N >= 10^(8)N \geq 10^{8}N≥108 the truncation errors are negligible and further increasing N N NNN over 10 10 10 10 10^(10)10^{10}1010 renders the results of the computations presented in the next section indistinguishable in the limit of double precision 4 4 ^(4){ }^{4}4.
It has been checked numerically that GRW converges as O ( δ x 2 ) + O ( 1 / N ) O δ x 2 + O ( 1 / N ) O(deltax^(2))+O(1//sqrtN)\mathcal{O}\left(\delta x^{2}\right)+\mathcal{O}(1 / \sqrt{N})O(δx2)+O(1/N) and when the condition 1 / N = O ( δ x 2 ) 1 / N = O δ x 2 1//sqrtN=O(deltax^(2))1 / \sqrt{N}=\mathcal{O}\left(\delta x^{2}\right)1/N=O(δx2) is met, the convergence is of the order O ( δ x 2 ) O δ x 2 O(deltax^(2))\mathcal{O}\left(\delta x^{2}\right)O(δx2) as for the finite differences scheme 4 4 ^(4){ }^{4}4. While sequential particle tracking cannot practical handle more than N = 10 9 N = 10 9 N=10^(9)N=10^{9}N=109 particles even in simple diffusion problems, the GRW algorithm based on (11) is insensitive to the increase of the number of particles.
Figure 1: Two-dimensional GRW algorithm for variable diffusion coefficients (left) and plume contours at t = 0 , 10 , 100 t = 0 , 10 , 100 t=0,10,100\mathrm{t}=0,10,100t=0,10,100, for a small Gaussian pulse as initial condition (right).
Figure 1: Two-dimensional GRW algorithm for variable diffusion coefficients (left) and plume contours at t = 0 , 10 , 100 t = 0 , 10 , 100 t=0,10,100\mathrm{t}=0,10,100t=0,10,100, for a small Gaussian pulse as initial condition (right).
The two-dimensional GRW uses two different parameters r x r x r_(x)r_{x}rx and r y r y r_(y)r_{y}ry, which also can be space-time variable, and are constrained by the relation r x + r y ≤ 1 r x + r y ≤ 1 r_(x)+r_(y) <= 1r_{x}+r_{y} \leq 1rx+ry≤1, which ensures the conservation of the number of particles and the stability of the numerical solutions. Accordingly, two jump amplitudes are chosen, d x d x d_(x)d_{x}dx and d y d y d_(y)d_{y}dy, such that for given space steps δ x δ x delta x\delta xδx and δ y δ y delta y\delta yδy the time step verifies the inequality
(12) δ t ≤ [ 2 D x max ( d x δ x ) 2 + 2 D y max ( d y δ y ) 2 ] − 1 , (12) δ t ≤ 2 D x max d x δ x 2 + 2 D y max d y δ y 2 − 1 , {:(12)delta t <= [(2D_(x)^(max))/((d_(x)delta x)^(2))+(2D_(y)^(max))/((d_(y)delta y)^(2))]^(-1)",":}\begin{equation*} \delta t \leq\left[\frac{2 D_{x}^{\max }}{\left(d_{x} \delta x\right)^{2}}+\frac{2 D_{y}^{\max }}{\left(d_{y} \delta y\right)^{2}}\right]^{-1}, \tag{12} \end{equation*}(12)δt≤[2Dxmax(dxδx)2+2Dymax(dyδy)2]−1,
where D x max D x max  D_(x)^("max ")D_{x}^{\text {max }}Dxmax  and D y max D y max  D_(y)^("max ")D_{y}^{\text {max }}Dymax  are the upper bounds of the diffusion coefficients D x ( x , y , t ) D x ( x , y , t ) D_(x)(x,y,t)D_{x}(x, y, t)Dx(x,y,t) and D y ( x , y , t ) D y ( x , y , t ) D_(y)(x,y,t)D_{y}(x, y, t)Dy(x,y,t).
In the following, the two-dimensional GRW is used to solve the modeled PDF equations (4-5). On the y y yyy-axis is now the concentration c , 0 < c < 1 c , 0 < c < 1 c,0 < c < 1c, 0<c<1c,0<c<1. GRW provides the number of particles as function of time at lattice sites ( i δ x , j δ c i δ x , j δ c i delta x,j delta ci \delta x, j \delta ciδx,jδc ). One obtains thus a numerical approximation of the joint PDF p ( c , x , t ) p ( c , x , t ) p(c,x,t)p(c, x, t)p(c,x,t), solution of the evolution equation (6). This is the main difference from the Lagrangian approaches in turbulence 1 , 2 , 3 , 10 1 , 2 , 3 , 10 ^(1,2,3,10){ }^{1,2,3,10}1,2,3,10 or hydrology 8 8 ^(8){ }^{8}8, where a grid-free particle method is used to compute trajectories sequentially and a mesh is used to estimate in a post-processing step the full PDF 8 PDF 8 PDF^(8)\mathrm{PDF}^{8}PDF8 or, in most cases, only its moments 3 , 10 3 , 10 ^(3,10){ }^{3,10}3,10. A schematic of the two-dimensional GRW rule, together with the contours of the "solute plume" in the ( x , c ) ( x , c ) (x,c)(x, c)(x,c) domain are shown in Fig. 1.

4 NUMERICAL TESTS

The new method, referred to in the following as GRW-PDF, was tested for the one dimensional transport problem describing the vertically integrated concentration c ( x , t ) c ( x , t ) c(x,t)c(x, t)c(x,t) in a two-dimensional heterogeneous aquifer. We considered the numerical setup of our Monte Carlo simulations done in the last years 11 , 12 11 , 12 ^(11,12){ }^{11,12}11,12, which also provided the upscaled diffusion coefficient D x ( t ) D x ( t ) D_(x)(t)\mathcal{D}_{x}(t)Dx(t) from (6), the upscaled constant velocity V = 1 m / V = 1 m / V=1m//\mathcal{V}=1 \mathrm{~m} /V=1 m/ day , and the initial PDF. The latter has been mutiplied by Abogadro's number N = 10 24 N = 10 24 N=10^(24)N=10^{24}N=1024 to obtain the initial distribution of particles. For a fixed discretization δ x = 0.1 m , δ c = 0.001 δ x = 0.1 m , δ c = 0.001 delta x=0.1m,delta c=0.001\delta x=0.1 \mathrm{~m}, \delta c=0.001δx=0.1 m,δc=0.001, we chose a time step δ t = 0.5 δ t = 0.5 delta t=0.5\delta t=0.5δt=0.5 days which fulfills the condition (12), with d x = 4 d x = 4 d_(x)=4d_{x}=4dx=4 and d y = 5 d y = 5 d_(y)=5d_{y}=5dy=5. In absence of a model for the parameter B B BBB in equation (6), which cannot be simply transfered from turbulence models to groundwater, we considered a constant diffusion coefficient in concentration space D c = 10 − 6 D c = 10 − 6 D_(c)=10^(-6)D_{c}=10^{-6}Dc=10−6, which is the coefficient of the Wiener process defined by (10) for the chosen discretization. The coefficient A A AAA was modeled as a combination of mixing through "interaction by exchange with the mean" (the first term in the equation (13) below) and a drift term Δ ⟨ C ⟩ Δ ⟨ C ⟩ Delta(:C:)\Delta\langle C\rangleΔ⟨C⟩ which incorporates effects of local scale transport 3 3 ^(3){ }^{3}3,
(13) A ( t ) = − a D x ( t ) λ [ C ( t ) − ⟨ C ⟩ ] + Δ ⟨ C ⟩ , (13) A ( t ) = − a D x ( t ) λ [ C ( t ) − ⟨ C ⟩ ] + Δ ⟨ C ⟩ , {:(13)A(t)=-a(D_(x)(t))/(lambda)[C(t)-(:C:)]+Delta(:C:)",":}\begin{equation*} A(t)=-a \frac{\mathcal{D}_{x}(t)}{\lambda}[C(t)-\langle C\rangle]+\Delta\langle C\rangle, \tag{13} \end{equation*}(13)A(t)=−aDx(t)λ[C(t)−⟨C⟩]+Δ⟨C⟩,
where λ λ lambda\lambdaλ is a characteristic length scale. A mean drift Δ ⟨ C ⟩ = − 0.004 ⟨ C ⟩ Δ ⟨ C ⟩ = − 0.004 ⟨ C ⟩ Delta(:C:)=-0.004(:C:)\Delta\langle C\rangle=-0.004\langle C\rangleΔ⟨C⟩=−0.004⟨C⟩, describing the attenuation of the mean concentration due to local scale diffusion, was inferred from a Gaussian distribution of the mean concentration at t = 20 t = 20 t=20t=20t=20 days, for the local diffusion coefficient D = 0.01 m / D = 0.01 m / D=0.01m//D=0.01 \mathrm{~m} /D=0.01 m/ day. A constant a = 1 a = 1 a=1a=1a=1 and λ = 1 λ = 1 lambda=1\lambda=1λ=1, the correlation scale of the hydraulic conductivity, have been found acceptable for the first tests of the algorithm and enabled us to produce the simulations for illustration purposes presented in the following. The runs take about 12 seconds on a laptop with 3 GB RAM and 2 GHz .
We checked first that in absence of mixing ( A = B = 0 A = B = 0 A=B=0A=B=0A=B=0 ) the initial PDF is passively transported without any alteration of its shape. Then, we summed up the number of particles over the c c ccc lattice coordinates to compute the marginal position-probability density, which is in fact the ensemble mean concentration. The excellent agreement with the results derived from an ensemble of Monte Carlo simulations 11 , 12 11 , 12 ^(11,12){ }^{11,12}11,12 shown in Fig. 2 proves that GRW-PDF fulfills the consistency requirement 10 10 ^(10){ }^{10}10, mentioned at the end of Section 2.
Figure 2: Vertically integrated concentration as function of x x xxx at times 10 and 100 (left) and as function of time at points on the mean flow trajectory x = V t x = V t x=Vtx=\mathcal{V} tx=Vt (right).
We also computed the "reconstructed" mean concentration (dot line in Fig. 2) by summing up the c c ccc lattice coordinates weighted by the corresponding PDF. The reconstructed mean underestimates the Monte Carlo results, indicating the lack of precision of the PDF approximation. This is not surprising, since no calibration of the GRW-PDF method has been done. Nevertheless, the comparison of cumulative distribution functions presented in fig. 3 demonstrates the ability of the GRW-PDF approach to produce dense representations of the PDF and to describe its evolution in time and space.
Figure 3: Cumulative distribution function at x = V t x = V t x=Vtx=\mathcal{V} tx=Vt transported in time by the GRW-PDF algorithm (left) compared with the Monte Carlo results (right).

REFERENCES

[1] S. B. Pope. PDF methods for turbulent reactive flows. Prog. Energy Combust. Sci., 11(2), 119-192, (1885).
[2] D. C. Haworth. Progress in probability density function methods for turbulent reacting flows. Prog. Energy Combust. Sci., 36, 168-259, doi:10.1016/j.pecs.2009.09.003, (2010).
[3] S. Viswanathan, H. Wang, and S. B. Pope. S. Numerical implementation of mixing and molecular transport in LES/PDF studies of turbulent reacting flows. J. Comput. Phys., 230, 6916-6957, doi:10.1016/j.jcp.2011.05.020, (2011).
[4] C. Vamoş, N. Suciu and H. Vereecken. Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys., 186(2), 527544, doi:10.1016/S0021-9991(03)00073-1, (2003).
[5] G. Dagan. Flow and Transport in Porous Formations, Springer, Berlin (1989).
[6] N. Suciu. Spatially inhomogeneous transition probabilities as memory effects for diffusion in statistically homogeneous random velocity fields. Phys. Rev. E, 81, 056301, (2010).
[7] P. E. Kloeden, and E. Platen. Numerical Solutions of Stochastic Differential Equations, Springer, Berlin, (1999).
[8] D. W. Meyer, P. Jenny and H. A. Tchelepi. A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media. Water Resour. Res., 46, W12522, (2010). doi:10.1029/2010WR009450.
[9] S. Attinger, M. Dentz, H. Kinzelbach and W. Kinzelbach. Temporal behavior of a solute cloud in a chemically heterogeneous porous medium. J. Fluid Mech., 386, 77-104, (1999).
[10] N. Suciu, C. Vamos, J. Vanderborght, H. Hardelauf and H. Vereecken. Numerical investigations on ergodicity of solute transport in heterogeneous aquifers. Water Resour. Res., 42, W04409, doi:10.1029/2005WR004546, (2006) .
[11] N. Suciu, C. Vamos, F. A. Radu, H. Vereecken and P. Knabner. Persisten memory of diffusing particles. Phys. Rev. E, 80, 061134, doi:10.1103/PhysRevE.80.061134, (2009).
[12] S. B. Pope. Lagrangian PDF methods for turbulent flows. Annu. Rev. Fluid Mech., 26, 23-63, (1994).
2012

Related Posts