A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media

7 years ago

(original), global random walk, ISI/JCR, Numerical Modeling, paper

Abstract

We identify sufficient conditions under which evolution equations for probability density functions (PDF) of random concentrations are equivalent to Fokker–Planck equations. The novelty of our approach is that it allows consistent PDF approximations by densities of computational particles governed by Itô processes in concentration–position spaces. Accurate numerical solutions are obtained with a global random walk (GRW) algorithm, stable, free of numerical diffusion, and insensitive to the increase of the total number of computational particles. The system of Itô equations is specified by drift and diffusion coefficients describing the PDF transport in the physical space, provided by up-scaling procedures, as well as by drift and mixing coefficients describing the PDF transport in concentration spaces. Mixing models can be obtained similarly to classical PDF approaches or, alternatively, from measured or simulated concentration time series. We compare their performance for a GRW-PDF numerical solution to a problem of contaminant transport in heterogeneous groundwater systems.

Authors

N. Suciu
-Mathematics Department, Friedrich-Alexander University of Erlangen-Nuremberg, Cauerstraße 11, 91058 Erlangen, Germany
-Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

F.A. Radu
-Department of Mathematics, University of Bergen, Allegaten 41, 5008 Bergen, Norway

S. Attinger
-Faculty for Chemistry and Earth sciences, University of Jena, Burgweg 11, 07749 Jena, Germany
-Department Computational Hydrosystems, UFZ Centre for Environmental Research, Permoserstraße 15, 04318 Leipzig, Germany

L. Schüler
-Faculty for Chemistry and Earth sciences, University of Jena, Burgweg 11, 07749 Jena, Germany
-Department Computational Hydrosystems, UFZ Centre for Environmental Research, Permoserstraße 15, 04318 Leipzig, Germany

P. Knabner
Mathematics Department, Friedrich-Alexander University of Erlangen-Nuremberg, Cauerstraße 11, 91058 Erlangen, Germany

Keywords

PDF methods; Mixing; Random walk; Porous media.

Cite this paper as:

N. Suciu, F.A. Radu, S. Attinger, L. Schüler, P. Knabner, A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media, J. Comput. Appl. Math., 289 (2015), 241-252,
doi:10.1016/j.cam.2015.01.030

PDF

A Fokker–Planck approach for probability distributions of species concentrations transported in heterogeneous media

About this paper

Journal

Journal of Computational and Applied Mathematics

Publisher Name

Elsevier

DOI

10.1016/j.cam.2015.01.030

Print ISSN

0377-0427

Online ISSN

Not available yet.

Google Scholar Profile

[1] N. Suciu, Diffusion in random velocity fields with applications to contaminant transport in groundwater, Advances in Water Resources, 69 (2014), pp. 114-133.
CrossRef (DOI)

[2] C. Vamoş, M. Crăciun, Separation of components from a scale mixture of Gaussian white noises, Physical Review E, 81 (2010) no. 5
CrossRef (DOI)

[3] C. Vamoş, M. Crăciun, Automatic trend estimation, Jan 2013, Springer, Dordrecht, 2012, pp. 131, ISBN: 978-94-007-4824-8 (Print), 978-94-007-4825-5 (Online),
CrossRef (DOI)

[4] N. Suciu, F.A. Radu, A. Prechtel, F. Brunner, P. Knabner, A coupled finite element–global random walk approach to advection-dominated transport in porous media with random hydraulic conductivity, Journal of Computational and Applied Mathematics, 246 (2013), pp. 27-37
CrossRef (DOI)

[5] F.A. Radu, N. Suciu, J. Hoffmann, A. Vogel, S. Attinger, Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: A comparative study, Advances in Water Resources, 34 (2011) no. 1, pp. 47-61.
CrossRef (DOI)

[6] R.W. Bilger, The Structure of Diffusion Flames, 1976, Combust Sci Technol
CrossRef (DOI)

[7] S.B. Pope, The Statistical Theory of Turbulent Flames, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 291 (1979) no. 1384, pp. 529-568
CrossRef (DOI)

[8] D.W. Meyer, P. Jenny, H.A. Tchelepi, A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media, Water Resources Research, 46 (2010) no. 12
CrossRef (DOI)

[9] F.A Jaberi, P. J. Colucci, S. James, P. Givi, S. B. Pope, Filtered Mass Density Function for Large-Eddy Simulation of Turbulent Reacting Flows, J Fluid Mech (1999)
CrossRef (DOI)

[10] Haworth D.C., Progress in probability density function methods for turbulent reacting flows, Prog. Energy Combust. Sci., 36 (2010), 168-259
CrossRef (DOI)

[11] N. Suciu, C. Vamoş, J. Vanderborght, H. Hardelauf, H. Vereecken, Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, 2006.
CrossRef (DOI)

[12] S. B. Pope, PDF Methods for Turbulent Reacting Flows, Prog Energ Combust, (1985)
CrossRef (DOI)

[13] R. McDermott, S. B. Pope, A particle formulation for treating differential diffusion in filtered density function methods, Journal of Computational Physics, 226 (2007) no 1, pp. 947-993
CrossRef (DOI)

[14] C. Vamoş, N. Suciu, H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, Journal of Computational Physics, 186 (2003), 527-544
CrossRef (DOI)

[15] Haifeng Wang, Pavel P. Popov, Stephen B. Pope, Weak second-order splitting schemes for Lagrangian Monte Carlo particle methods for the composition PDF/FDF transport equations, Journal of Computational Physics, 229 (2010) issue 5, pp. 1852-1878
CrossRef (DOI)

[16] Venkatramanan Raman, Heinz Pitsch, Rodney O., Fox, Hybrid large-eddy simulation/Lagrangian filtered-density-function approach for simulating turbulent combustion, Combust Flame (2005)
CrossRef (DOI)

[17] S. B. Pope, The probability approach to the modelling of turbulent reacting flows, Combustion and Flame, 27 (1976), pp. 299-312
CrossRef (DOI)

[18] R. O. Fox, Computational Models for Turbulent Reacting Flow, 2003

[19] P.E. Kloeden, E. Platen, Numerical Solutions of Stochastic Differential Equations, Springer, Berlin, 1999.

[20] A.Y. Klimenko, On simulating scalar transport by mixing between Lagrangian particles, Physics of Fluids, 19 (2007) no. 3, p 031702
CrossRef (DOI)

[21] S. B. Pope, Advances in PDF Methods for Turbulent Reactive Flows
CrossRef? yr? J? bk? ISBN

Paper in HTML form

1-s2.0-S0377042715000424-main

A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media

N. Suciu $^{a, b, *}$ , F.A. Radu $^{c}$ , S. Attinger $^{d, e}$ , L. Schüler $^{d, e}$ , P. Knabner $^{a}$ $^{a}$ Mathematics Department, Friedrich-Alexander University of Erlangen-Nuremberg, Cauerstraße 11, 91058 Erlangen, Germany $^{b}$ Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Fantanele 57, 400320 Cluj-Napoca, Romania $^{c}$ Department of Mathematics, University of Bergen, Allegaten 41, 5008 Bergen, Norway $^{d}$ Faculty for Chemistry and Earth sciences, University of Jena, Burgweg 11, 07749 Jena, Germany $^{e}$ Department Computational Hydrosystems, UFZ Centre for Environmental Research, Permoserstraße 15, 04318 Leipzig, Germany

ARTICLE INFO

Article history:
Received 16 August 2014
Received in revised form 5 November 2014

MSC:
60J60
60G60
86A05
Keywords:
PDF methods
Mixing
Random walk
Porous media

Abstract

We identify sufficient conditions under which evolution equations for probability density functions (PDF) of random concentrations are equivalent to Fokker-Planck equations. The novelty of our approach is that it allows consistent PDF approximations by densities of computational particles governed by Itô processes in concentration-position spaces. Accurate numerical solutions are obtained with a global random walk (GRW) algorithm, stable, free of numerical diffusion, and insensitive to the increase of the total number of computational particles. The system of Itô equations is specified by drift and diffusion coefficients describing the PDF transport in the physical space, provided by up-scaling procedures, as well as by drift and mixing coefficients describing the PDF transport in concentration spaces. Mixing models can be obtained similarly to classical PDF approaches or, alternatively, from measured or simulated concentration time series. We compare their performance for a GRW-PDF numerical solution to a problem of contaminant transport in heterogeneous groundwater systems.

1. Introduction

PDF methods have been developed in turbulence and combustion theory, motivated by the need to close the Reynolds averaged transport equations. Beyond achieving the closure (e.g. for convection and reaction terms) these approaches may be valuable tools in modeling transport in random environments in that they provide the one-dimensional (one-point onetime) PDF of transported quantities, such as turbulent velocity fluctuations and scalars (concentrations, temperature, enthalpy) [1-3]. In this paper, we only consider scalar PDF methods, with an illustration for the concentration of chemical species transported in saturated groundwater formations. The development of the Fokker-Planck approach and the numerical illustration presented in the following are restricted to the case of incompressible flows and constant diffusion coefficients.

Thus, we consider an array of species concentrations

C (x, t) = {C_{α} (x, t)}

linked through reaction terms

S (C) = {S_{α} (C)}

α = 1, \dots, N_{α}

, where

N_{α}

is the number of chemical species. The diffusion coefficient

D

is constant and species independent

and the chemical components are transported in a statistically homogeneous random velocity field

V

with divergence-free realizations according to a system of

N_{α}

stochastic balance equations, written in a compact form as

\begin{matrix} (1) & \partial_{t} C + V \cdot \nabla C = D \nabla^{2} C + S (C) \end{matrix}

The Eulerian

PDF f (c; x, t)

of the random concentrations

C

solving (1) satisfies

\begin{matrix} (2) & \partial_{t} f + \nabla \cdot (V f) = \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} (D_{i j} f) - \frac{\partial^{2}}{\partial c_{α} \partial c_{β}} (M_{α β} f) - \nabla_{c} \cdot (S f) \end{matrix}

where

V = ⟨ V ⟩ + \nabla \cdot D

is the upscaled drift and

⟨ V ⟩

the mean velocity,

D = D + D^{*}

is the diffusion coefficient tensor upscaled through the gradient-diffusion closure

f ⟨ V - ⟨ V ⟩ ∣ C (x, t) = c ⟩ = - D^{*} \nabla f

of the mean velocity fluctuations conditional on concentration,

\nabla_{c}

is the gradient in the concentration space, and

M = ⟨ D \nabla C \nabla C ∣ C (x, t) = c ⟩

is the conditional dissipation rate tensor which accounts for molecular mixing [2, Section 12.7.4]. The reaction term

S

occurring in (2) is in a closed form, the same as in the concentration balance equation (1), which is the essential advantage of the PDF methods [3].

The Eulerian PDF

f (c; x, t)

is usually obtained from the numerical solution of a system of Itô equations describing the evolution of an ensemble of computational particles in physical and concentration spaces [3, Section 6.7]. The system of Itô equations can be written in a general form as

\begin{aligned} (3) & d X (t) & = V (X (t), t) d t + d W (X (t), t) \\ (4) & d C (t) & = M (C (t), X (t), t) d t + S (C (t)) d t \end{aligned}

where

C (t) = C (X (t), t), W

is a Wiener process with

E {W (X (t), t)} = 0

and

E {W^{2} (X (t), t)} = 2 \int_{0}^{t} D (X (t), t^{'}) d t^{'}

[4, Section 7.1], and

M

is a mixing model for the conditional dissipation rate

M

[1-3].

In this paper, we identify consistency conditions relating the statistics of the random field

C (x, t)

, solution of the system of stochastic equations (1), to that of the stochastic process

{X (t), C (t)}

, described by the Itô equations (3) and (4), which reveals the natural relationship between the Eulerian

PDF f (c; x, t)

and the joint concentration-position

PDF p (c, x, t)

. Further, we establish sufficient conditions under which

p (c, x, t)

solves a closed-form equation and, in particular, an equation of the form (2). Such equations are thus Fokker-Planck equations associated to systems of Itô equations (3) and (4). By equivalence between Itô and Fokker-Planck representations of a diffusion process, ensured by regularity conditions on coefficients [5] which are assumed fulfilled,

p (c, x, t)

is approximated numerically by the density of a large ensemble of computational particles which evolve in the (

x, c

) space accordingly to Eqs. (3) and (4) [4].

The classical approach to solve the PDF equation (2) is discussed in Section 2. The consistency conditions and the correspondence between concentration random fields and stochastic processes in concentration-position spaces, as well as the derivation of the Fokker-Planck equation for the joint concentration-position PDF, are presented in Section 3. Section 4 introduces the new solution approach for the concentration PDF using the GRW algorithm. An application to non-reactive transport in saturated aquifers is presented in Section 5 and some conclusions are formulated in Section 6.

2. Classical PDF methods

The solution of the system of Itô equations (3) and (4) is the stochastic process

{X_{i} (t), C_{α} (t)}, C_{α} (t) = C_{α} (X (t), t), t \geq 0

i = 1, 2, 3, α = 1, \dots, N_{α}

. At given time

t

, the process takes values (

x = X (t), c = C (t)

). Here and throughout the paper, we use the convention which denotes random functions by capital letters and their values in the state space by corresponding lowercase letters. The state space

Ω = Ω_{x} \times Ω_{c}

is the Cartesian product of the three-dimensional physical space

Ω_{x}

and the

N_{α}

-dimensional concentration space

Ω_{c}

. The one-dimensional (one-time) PDF of this process is the joint PDF of concentrations and positions

p (c, x, t)

. Classical PDF methods approximate the PDF

f (c; x, t)

by solving a system of stochastic differential equations stochastically equivalent to Eqs. (3) and (4) with coefficients defined by averages conditional on concentration, which correspond to the coefficients of the PDF equation ( 2 ) [

1, 3, 6

For constant density flows, the one-dimensional Eulerian

PDF f (c; x, t)

of the random field

C (x, t)

is assimilated to the conditional PDF of the stochastic process

{X_{i} (t), C_{α} (t)}, p (c ∣ x, t) = p (c, x, t) / p_{x} (x, t)

, where

p_{x} (x, t)

is the position PDF. A necessary ingredient to render this approach feasible is to force a uniform position PDF

p_{x} (x, t)

, so that

f (c; x, t)

becomes proportional to the solution

p (c, x, t)

of the Fokker-Planck equation associated to the system of Itô equations (3) and (4). The position PDF

p_{x} (x, t)

is the solution of the Fokker-Planck equation associated to the position Itô equation (3), which is Eq. (2) with the last two terms set to zero. The spatially uniform position PDF

p_{x} (x, t) =

const

= 1 / \int_{Ω_{x}} d x

is also time independent (

Ω_{x}

is a fixed domain). This is a trivial solution of the Fokker-Planck equation in case of space-constant velocity and diffusion coefficients [7] or, in the more general case of space-variable drift and isotropic diffusion coefficient, when the condition

\nabla \cdot V = \nabla^{2} D

is fulfilled [2, Section 12.7].

The general solution algorithm for variable density flows is based on a discrete representation of the PDF which requires the equality between the joint concentration-position PDF multiplied by the total mass

M

of fluid in

Ω_{x}, M p (c, x, t)

, and the "mass density function"

ρ (c) f (c; x, t)

, where

ρ (c)

is the density of the fluid which, for constant temperature, depends only on species concentrations [1, Eqs. (3.78) and (3.79)]. For the case of scalar transport considered here, one assumes that the carrying fluid is also an element of the chemical composition described by the concentration vector

C

. Equating the integrals
over the concentration space

Ω_{c}

of the two densities one finds that

M p_{x} (x, t) = \int_{Ω_{x}} ρ (c) f (c; x, t) d c = ⟨ ρ ⟩ (x, t)

. Since for constant density flows

⟨ ρ ⟩

is independent of

x

, the position

PDF p_{x} (x, t)

has to be uniform, irrespective of the properties of the coefficients

V

and

D

. We note that while a uniform

p_{x} (x, t)

is physically meaningful in case of homogeneous turbulence with space-independent statistics [1, p. 140], it is not consistent with a general mechanism of the transport process. Even for divergence-free upscaled velocity

V

, the condition

\nabla^{2} D = 0

still has to be fulfilled for the existence of a uniform solution

p_{x} (x, t)

of the position Fokker-Planck equation.

The implication of a uniform position PDF for the solution algorithm is that the expected particle number density is uniform in physical space at all times. Instead of representing a point in the (

3 + N_{α}

)-dimensional space

Ω

which describes the state of the process

{X_{i} (t), C_{α} (t)}

at a given time, the "notional particles" introduced by Pope [1] carry a "composition"

{C_{α} (t)}, α = 1, \dots, N_{α}

, of the chemical system in the physical space

Ω_{x}

. Concentration evolution equations similar to Eq. (4) are solved for each particle to update its composition and, in a subsequent step, the system of particles is transported in the physical space. Finally, mean values and PDFs are estimated by counting the number of particles with different compositions in cells of the physical space. In this way, one avoids the curse of dimensionality for multi-component reactive transport problems. The drawback is the linear increase of computational costs with the number of particles and the numerical diffusion produced by interpolation to particle positions of mean values estimated by averaging over cells [8].

The equivalence between the system of notional particles and the physical system is established by requiring the equality between the joint PDF

p (c, x, t)

of the system of notional particles and that of a system of fluid particles with the same composition

{C_{α} (t)}

(see e.g. [1, Section 4.2] and [3, Section 6.7]). The mathematical definition of the fluid particle uses Lagrangian equations for the fluid flow which allow a consistent derivation of the evolution equation for the Eulerian velocity PDF [2, Section 12.6]. An intermediate situation is that of reaction-diffusion systems in turbulent flows described by fluid particles with trajectories specified by the upscaled flow velocity, through Eq. (3) with the diffusion term set to zero, and compositions modeled by Eq. (4). The latter accounts for the diffusion process through a mixing term which models the conditional average of the diffusion flux according to the transport balance equation

(1) [1, 3]

. This formulation of the problem already introduces a disagreement between the Lagrangian description and the diffusion process, which destroys the individuality of the fluid particles. This is even better seen in the general form of the equations governing the system of notional particles, where the particle's position is an Itô diffusion process of type (3) [3,6,7]. Because Eq. (3) describes a diffusion process, there is an infinity of trajectories starting from the same initial state [5] and a fluid particle uniquely specified by its initial position cannot be defined. Nevertheless, the heuristic approach based on the equivalence between notional and fluid particles carrying compositions in physical space succeeded in designing efficient solution algorithms in a broad field of applications of the PDF methods [9].

3. The Fokker-Planck approach

To avoid the above conceptual inconsistency and to overcome the numerical issues of the classical PDF approach, we derive the Fokker-Planck equation for the joint concentration-position PDF,

p (c, x, t)

, consistent with the transport equations (1). Further, we estimate

p (c, x, t)

from a direct GRW solution of the system of Itô equations (3) and (4) with the initial condition consisting of a distribution of computational particles in the state space

Ω

which approximates the initial joint concentration-position PDF.

3.1. Correspondence between random fields and stochastic processes

The implementation of particle methods for solving PDF equations requires relationships between the Eulerian PDF

f (c; x, t)

which solves (2) and the joint concentration-position PDF

p (c, x, t)

associated to the system of Itô equations (3) and (4). A difficulty, not yet explicitly formulated in the literature, is caused by the use of two quite different stochastic models: random fields and stochastic processes. On one side, the Eulerian

PDF f (c; x, t)

is the one-point, one-time PDF of the random concentration

C (x, t)

, which is a random field indexed by three spatial parameters and time. On the other side,

p (c, x, t)

is the one-time PDF of the stochastic process

{C (t), X (t)}

described by the system of Itô equations (3) and (4), with a single parameter, the time

t

. While in the Eulerian picture

x

is a parameter of the random concentration,

X (t)

is a stochastic process with the one-time

PDF p_{x} (x, t)

. The position

PDF p_{x} (x, t)

is the normalized concentration at a given time

t

of the points

X (t)

on trajectories of the Itô process (3) [4]. Since (3) describes an upscaled transport process, it follows that

p_{x} (x, t)

represents an ensemble mean concentration of points

X (t)

, irrespective of the species concentrations associated with them via Eq. (4).

Using relations between joint and marginal probabilities,

p_{x} (x, t)

can be expressed as

\begin{aligned} p_{x} (x, t) & = \int_{Ω_{c}} p (c, x, t) d c \\ = \int_{Ω_{c_{1}}} d c_{1} \int_{Ω_{c} ∖ Ω_{c_{1}}} p (c, x, t) d c_{2} \dots d c_{α} = \int_{Ω_{c_{1}}} p_{c_{1}, x} (c_{1}, x, t) d c_{1} = \dots = \int_{Ω_{c_{N_{α}}}} p_{c_{N_{α}}, x} (c_{N_{α}}, x, t) d c_{N_{α}} \\ = \frac{1}{N_{α}} \sum_{α = 1}^{N_{α}} \int_{Ω_{c_{α}}} p_{c_{α}, x} (c_{α}, x, t) d c_{α} \end{aligned}

that is, as an arithmetic average of position PDFs associated to each species,

p_{x}^{(α)} = \int_{Ω_{c_{α}}} p_{c_{α}, x} (c_{α}, x, t) d c_{α}

, where

p_{c_{α}, x} (c_{α}

x, t)

are marginal PDFs of

p (c, x, t)

. Since

p_{x}^{(α)}

have the meaning of ensemble average concentrations of species

α

, we conjecture that the Eulerian counterpart of the above expression is given in terms of the random field

C (x, t)

by the stochastic average of the arithmetic mean of the species concentrations,

\begin{matrix} (5) & Θ (C (x, t)) = \frac{1}{C^{*} N_{α}} \sum_{α = 1}^{N_{α}} C_{α} (x, t) \end{matrix}

where

C^{*}

is a normalization constant. This conjecture is true if

Θ

is a conserved scalar which verifies Eq. (1) without reaction term. Then, its average solves Eq. (2) with the last two terms set to zero (e.g. [7]), which is a Fokker-Planck equation associated to the Itô equation (3). This is always the case if the components of the concentration vector

C

include all the chemical species of the reaction system. Indeed, since each

C_{α}

is a linear combination of chemical elements [10] and the sum of chemical elements in the reaction system is conserved, the sum of species concentrations

\sum_{α = 1}^{N_{α}} C_{α}

is conserved as well.

Equating the position PDF

p_{x} (x, t)

with the mean

⟨ Θ ⟩

, expressed by the usual change of variables relating integrals over probability spaces and state spaces [5] as average over the state space

Ω_{c}

weighted by the Eulerian PDF, one obtains the relations

\begin{aligned} \int_{Ω_{c}} Θ (c) f (c; x, t) d c & = ⟨ Θ ⟩ (x, t) \\ = p_{x} (x, t) \\ (6) & = \int_{Ω_{c}} p (c, x, t) d c \end{aligned}

The second equality in (6) allows the determination of the normalization constant

C^{*}

. These relations are obviously satisfied if one chooses

\begin{matrix} (7) & Θ (c) f (c; x, t) = p (c, x, t) \end{matrix}

The relation (7) establishes a consistent correspondence between the one-dimensional statistics of the random field

C (x, t)

and the process

{C (t), X (t)}

. By virtue of (7),

p (c, x, t)

also solves the evolution equation for

F (c; x, t) = Θ (c) f (c; x, t)

As we have seen above, the weighting function

Θ

in (7) must be a conserved scalar, but its choice is not unique. It may be defined by the sum or by the arithmetic mean of the species concentrations, as well as by any other conserved combination of concentrations, for example the mixture fraction used in diffusion flames approaches [10]. For a given choice of

Θ

, (7) establishes a particular relation between

f (c; x, t)

and

p (c, x, t)

. For instance, in combustion theory the weighting factor is chosen as

Θ = \sum_{α = 1}^{N_{α}} ρ_{α} = ρ

, where

ρ_{α}

represents mass concentrations of species

α

and

ρ

is the fluid density [1,9]. In this case, both the solutes and the solvent have to be included in the ensemble of species

α

to close the problem. The relation (7) takes the form

\begin{matrix} (8) & ρ (c) f (c; x, t) = M p (c, x, t) \end{matrix}

which by integration over the concentration state space variable

c

yields

\begin{matrix} (9) & ⟨ ρ ⟩ (x, t) = M p_{x} (x, t), \end{matrix}

where

M = \int_{Ω} ⟨ ρ ⟩ (x, t) d x

is the total mass of fluid in

Ω_{x}

(see e.g. [1, Section 3.4] and [9, Section 4]). If the relation (8) is used instead of (7), then

F = ρ (c) f (c; x, t)

defines the mass density function (e.g. [1, Eq. (3.59)]). For constant density flows, (9) implies a uniform position PDF

p_{x} (x, t)

, which, as seen in Section 2, may be inconsistent with the spatial variability of the coefficients occurring in Eq. (1). The advantage of using (7) with a weighting function

Θ

different from

ρ

is that it does not force a uniform position PDF for constant density and is adequate for the case of dilute solutions, when the solvent is not included among the species concentrations.

3.2. Derivation of the Fokker-Planck equation

In studies of turbulent reacting flows [1,3], evolution equations for

f (c; x, t)

or for the weighted concentration PDF

F (c, x, t)

are derived by integrating over the velocity space evolution equations for joint velocity-concentration PDFs. In the following, we derive evolution equations for

F (c, x, t)

in the form of Fokker-Planck equations which, according to (7), also govern the evolution of the joint concentration-position PDF

p (c, x, t)

of the Itô process (3)-(4).

Proposition 1. Assume that the velocity field

V

is statistically homogeneous with divergence-free samples, the diffusion coefficient is constant and species independent, the system of reacting species is governed by the system of equations (1), and the weighting function

Θ

verifies the continuity equation

\begin{matrix} (10) & \partial_{t} Θ + V \cdot \nabla Θ = 0 \end{matrix}

Then, the weighted concentration PDF

F (c, x, t)

satisfies the evolution equation

\begin{matrix} (11) & \partial_{t} F + ⟨ V ⟩ \cdot \nabla F = - \nabla \cdot [⟨ U ∣ C = c ⟩ F] - \nabla_{c} \cdot [⟨ D \nabla^{2} C ∣ C = c ⟩ F] - \nabla_{c} \cdot (S F) . \end{matrix}

Proof. Let us consider the differential operator

\begin{matrix} (12) & A = \partial_{t} + V \cdot \nabla \end{matrix}

applied to an arbitrary differentiable function

Q

which depends on space-time variables only through the random concentration vector,

Q (x, t) = Q (C (x, t))

, with compact support in

Ω_{c}^{0}

(the interior of

Ω_{c}

). Similarly to the Pope derivation of the velocity-composition PDF equation [1], we derive Eq. (11) by equating two independent expressions for the ensemble average of the product between

Θ

and the operator (12) applied to

Q

The first expression is given by

\begin{matrix} (13) & ⟨ Θ A Q ⟩ = ⟨ Θ \partial_{t} Q + Θ V \cdot \nabla Q ⟩ = ⟨ \partial_{t} (Θ Q) + \nabla \cdot (V Θ Q) ⟩ = \partial_{t} ⟨ Θ Q ⟩ + ⟨ V ⟩ \cdot \nabla ⟨ Θ Q ⟩ + \nabla \cdot ⟨ U Θ Q ⟩, \end{matrix}

where we used the continuity equation (10), the incompressibility condition

\nabla \cdot V = 0

, and by

U = V - ⟨ V ⟩

we defined the fluctuation about the constant mean

⟨ V ⟩

of the statistical homogeneous velocity field. To compute the ensemble average (13) we also need information not included in the one-point

PDF (c; x, t)

, in this case, the statistics of the velocity fluctuations. Following Fox [3, Section 6.2], we lump the unknown statistics in a random vector

Z

and consider the joint Eulerian

PDF f (c, z; x, t)

. The ensemble average of a function

F (C, Z) = F_{1} (C) F_{2} (Z)

will then be computed as

\begin{matrix} (14) & ⟨ F (C (x, t), Z (x, t)) ⟩ = \int_{Ω_{c}} \int_{Ω_{z}} F_{1} (c) F_{2} (z) f (c, z; x, t) d c d z = \int_{Ω_{c}} F_{1} (c) ⟨ F_{2} ∣ C = c ⟩ f (c; x, t) d c \end{matrix}

where the conditional average is performed with respect to the conditional

PDF f (z ∣ c; x, t) = f (c, z; x, t) / f (c; x, t)

\begin{matrix} (15) & ⟨ F_{2} ∣ C = c ⟩ = \int_{Ω_{z}} F_{2} (z) f (z ∣ c; x, t) d z \end{matrix}

With (14) and (15), the ensemble average (13) of

Θ A Q

becomes

\begin{matrix} (16) & ⟨ Θ A Q ⟩ = \int_{Ω_{c}} Q (c) Θ (c) {\partial_{t} f (c; x, t) + ⟨ V ⟩ \cdot \nabla f (c; x, t) + \nabla \cdot [⟨ U ∣ C = c ⟩ f (c; x, t)]} d c \end{matrix}

The second expression of

⟨ Θ A Q ⟩

follows from the fact that

Q

depends on time-space variables through the random concentration

C (x, t)

and

A C

is given by the right hand side of the transport equation (1):

\begin{aligned} ⟨ Θ A Q ⟩ = ⟨ Θ \nabla_{c} Q \cdot A C ⟩ & = ⟨ Θ \nabla_{c} Q \cdot [D \nabla^{2} C + S (C)] ⟩ \\ (17) & = \int_{Ω_{c}} Θ (c) \nabla_{c} Q (c) \cdot [⟨ D \nabla^{2} C ∣ C = c ⟩ + S (c)] f (c; x, t) d c \end{aligned}

To obtain (17) we used the chain rule of calculus and the averaging procedure (14). Integration by parts yields

\begin{aligned} ⟨ Θ A Q ⟩ = & \int_{\partial Ω_{c}} Q (c) {Θ (c) [⟨ D \nabla^{2} C ∣ C = c ⟩ + S (c)] f (c; x, t)} \cdot n d Γ_{c} \\ - \int_{Ω_{c}} Q (c) \nabla_{c} \cdot {Θ (c) [⟨ D \nabla^{2} C ∣ C = c, ⟩ + S (c)] f (c; x, t)} d c \end{aligned}

where

n

is the outward unit normal to the boundary

\partial Ω_{c}

of the concentration space

Ω_{c}

and

d Γ_{c}

is the surface element. Since the function

Q (c)

has compact support in

Ω_{c}

it vanishes on

\partial Ω_{c}

and cancels the first integral, so that we finally obtain

\begin{matrix} (18) & ⟨ Θ A Q ⟩ = - \int_{Ω_{c}} Q (c) \nabla_{c} \cdot {Θ (c) [⟨ D \nabla^{2} C ∣ C = c, X = x, t ⟩ + S (c)] f (c; x, t)} d c \end{matrix}

Since the expressions (16) and (18) should give the same result for any function

Q

with compact support, equating them one obtains the evolution equation (11) for

F (c; x, t) = Θ (c) f (c; x, t)

Eq. (11) is not yet in a closed form and the conditional average of the velocity fluctuations and that of the diffusive flux from the right hand side of (11) require modeling. Often used in turbulence [

1, 3, 9

] (with some caution in case of variable density flows) is a gradient-diffusion closure of the first conditional average in (11),

\begin{matrix} (19) & ⟨ U ∣ C = c ⟩ f (c; x, t) = - D^{*} \nabla f (c; x, t) . \end{matrix}

Introduced in (11), the closure (19) gives the diffusion term

\frac{\partial}{\partial x_{i}} D_{i, j}^{*} \frac{\partial}{\partial x_{j}} F

. To put it in a Fokker-Planck form we add

\frac{\partial}{\partial x_{i}} F \frac{\partial}{\partial x_{j}} D_{i, j}^{*}

in both sides of Eq. (11), which also implies defining the new drift coefficients

\begin{matrix} (20) & V_{i} = ⟨ V_{i} ⟩ + \frac{\partial}{\partial x_{j}} D_{i, j}^{*} \end{matrix}

Proposition 1 ensures the existence of the evolution equation (11) only if the weighting factor verifies the continuity equation (10). This is the case when

Θ = ρ

and Eq. (11) solves for the mass density function

F (c; x, t) = ρ (c) f (c; x, t)

. In
case of dilute solutions, when the solvent is not included in the reaction system, an evolution equation for

F

can still be derived for advection-reaction problems. Then,

D = 0

, all conserved scalars

Θ

verify (10), and we have the following corollary.

Corollary 1. Assume

D = 0

and the closure (19). Then the weighted concentration PDF

F (c, x, t)

satisfies the evolution equation

\begin{matrix} (21) & \partial_{t} F + \nabla \cdot (V f) = \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} (D_{i j}^{*} F) - \nabla_{c} \cdot (S F) \end{matrix}

According to Eq. (7),

F (c; x, t) = p (c, x, t)

. Hence, (11) is the general form of the Fokker-Planck equation associated to the system of Itô equations (3) and (4). Expressing

p (c, x, t)

as a product of conditional PDF

p (c ∣ x, t)

and position PDF

p_{x} (x, t), p (c, x, t) = p (c ∣ x, t) p_{x} (x, t)

, and using (6) and (7) one obtains

\begin{matrix} (22) & p (c ∣ x, t) = Θ (c) f (c; x, t) / ⟨ Θ ⟩ (x, t) . \end{matrix}

Eq. (22) shows that the concentration PDF conditional on position, approximated by numerical solutions of the Itô equations (3) and (4), is a numerical solution for the one-dimensional Eulerian concentration PDF weighted by

Θ / ⟨ Θ ⟩

. In classical PDF methods

Θ = ρ

, Eq. (11) describes the evolution of the mass density function

F (c; x, t) = ρ (c) f (c; x, t)

, and the conditional PDF

p (c ∣ x, t)

of the system of computational particles approximates, according to (8) and (9), the density weighted Eulerian PDF,

\begin{matrix} (23) & p (c ∣ x, t) = ρ (c) f (c; x, t) / ⟨ ρ ⟩ . \end{matrix}

The solution of the classical "stochastic particles algorithm" of Pope is based on Eq. (23), which allows approximating density weighted averages by ensemble averages over the system of computational particles (see e.g. [1, Eqs. (3.84) and (3.89)]). But, even if in this case (11) is a Fokker-Planck equation, the Pope particles algorithm, shortly described in Section 2, does not solve the system of Itô equations (3) and (4).

The numerical solution for the weighted PDF based on Eq. (22) which does not require including the solvent among the ensemble of species concentrations has the advantage that it does not force a uniform position PDF for constant density flows. As shown by Corollary 1, this approach is valid for advection-reaction problems, which also may be thought as a consistent first-order approximation of the full, advection-diffusion-reaction problem [4]. Another situation when a weighting function

Θ \neq ρ

could be used is when the approximation

Θ ≃ ⟨ Θ ⟩ (x, t)

may be adopted, as in case of "ergodic behavior" of transport in groundwater

[4, 11]

. Then the Eulerian

PDF (c; x, t)

can be approximated by the conditional

PDF p (c ∣ x, t)

, even if

p_{x} (x, t) = ⟨ Θ ⟩ (x, t)

is not uniform. Within this ergodic approximation, Eq. (2) is the Fokker-Planck equation describing the evolution of

p (c ∣ x, t)

3.3. The treatment of the conditional diffusion term in Eq. (11)

The conditional averages in (11) depend not only on independent state-space variables

c

and

x

but also on space-time dependent variables describing the random field

C (x, t)

. A way to relate the two types of variables is to use probability densities defined by averages of Dirac-

δ

functions

[1, 9]

. Even though such densities are singular, their integration with respect to the state space variables yields well-defined statistics [4, Appendix A.1].

In this frame, a useful expression for the conditional average (15) is obtained by using the definition of the joint PDF in the (

c, z

) space,

f (c, z; x, t) = ⟨ δ (Z (x, t) - z) δ (C (x, t) - c) ⟩

, as follows:

\begin{aligned} ⟨ F ∣ C = c ⟩ (x, t) & = \int_{Ω_{z}} F (z) f (z ∣ c; x, t) d z = \int_{Ω_{z}} F (z) \frac{f (c, z; x, t)}{f (c; x, t)} d z \\ = \frac{1}{f (c; x, t)} \int_{Ω_{z}} F (z) ⟨ δ (Z (x, t) - z) δ (C (x, t) - c) ⟩ d z \\ = \frac{1}{f (c; x, t)} ⟨ δ (C (x, t) - c) \int_{Ω_{z}} F (z) δ (Z (x, t) - z) d z ⟩ \\ (24) & = \frac{1}{f (c; x, t)} ⟨ F (Z (x, t)) δ (C (x, t) - c) ⟩ \end{aligned}

In the above calculations we used multi-dimensional

δ

functions,

δ (y_{0} - y) = \prod_{i = 1}^{i = n} δ (y_{0, i} - y_{i})

, and the obvious property

\int_{Ω} f (y) δ (y_{0} - y) d y = \int_{R^{n}} I_{Ω} (y) f (y) δ (y_{0} - y) d y = f (y_{0})

for all

y \in Ω \subset R^{n}

, where

I_{Ω}

is the indicator function of

Ω

. Expressions of conditional averages similar to (24) are often used in classical PDF methods [1,2].

Another ingredient to work out the conditional averages in (11) is the derivative of the

δ

function. This is properly defined in terms of Dirac distributions by

\begin{matrix} (25) & \int_{- \infty}^{\infty} δ^{'} (y_{0} - y) φ (y) d y = - \int_{- \infty}^{\infty} δ (y_{0} - y) φ^{'} (y) d y \end{matrix}

for any smooth function

φ

with compact support in

R

. The usual notation for the distributional derivative is

δ^{'} [φ] = - δ [φ^{'}] = - φ^{'} (y_{0})

. Using (25) for a composite function

φ (g (x))

we have

\begin{aligned} \frac{d}{d x} \int_{- \infty}^{\infty} δ (g (x) - y) φ (y) d y & = \frac{d}{d x} φ (g (x)) = φ^{'} (g (x)) \frac{d g (x)}{d x} = \frac{d g (x)}{d x} \int_{- \infty}^{\infty} δ (g (x) - y) φ^{'} (y) d y \\ = - \frac{d g (x)}{d x} \int_{- \infty}^{\infty} δ^{'} (g (x) - y) φ (y) d y \end{aligned}

or, in compact notation,

\begin{matrix} (26) & \frac{d}{d x} δ (g (x) - y) [φ] = - \frac{d g (x)}{d x} δ^{'} (g (x) - y) [φ] . \end{matrix}

Bearing in mind its interpretation as distributional derivative, the relation (26) is formally written in applications as [1,2]

\begin{matrix} (27) & \frac{d}{d x} δ (g (x) - y) = - \frac{d g (x)}{d x} \frac{d}{d y} δ (g (x) - y) . \end{matrix}

For the vectorial function

g = C (x, t)

, (27) can be generalized to the following expression of the partial derivative with respect to the space variable,

\begin{matrix} (28) & \frac{\partial}{\partial x_{i}} δ (C (x, t) - c) = - \frac{\partial C_{α}}{\partial x_{i}} (x, t) \frac{\partial}{\partial c_{α}} δ (C (x, t) - c) . \end{matrix}

Proposition 2. If the weighting factor

Θ \equiv 1

the conditional diffusion term in Eq. (11) can be decomposed as sum of diffusion terms in physical and concentration spaces,

\begin{matrix} (29) & - \frac{\partial}{\partial c_{α}} [⟨ D \frac{\partial^{2} C_{α}}{\partial x_{i} \partial x_{i}} | c ⟩ f (c; x, t)] = D \frac{\partial^{2}}{\partial x_{i} \partial x_{i}} f (c; x, t) - \frac{\partial^{2}}{\partial c_{α} \partial c_{β}} [⟨ D \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial C_{β}}{\partial x_{i}} | c ⟩ f (c; x, t)] . \end{matrix}

Proof.

Θ \equiv 1

implies

F (c; x, t) \equiv f (c; x, t)

. Using (24), (28), and the definition

⟨ δ (C (x, t) - c) ⟩ = f (c; x, t)

for the onedimensional Eulerian PDF, the second conditional average from (11) can be computed as follows:

\begin{aligned} - \frac{\partial}{\partial c_{α}} [⟨ D \frac{\partial^{2} C_{α}}{\partial x_{i} \partial x_{i}} | c ⟩ f (c; x, t)] & = - \frac{\partial}{\partial c_{α}} ⟨ D \frac{\partial^{2} C_{α}}{\partial x_{i} \partial x_{i}} δ (C (x, t) - c) ⟩ \\ = - D \frac{\partial}{\partial c_{α}} ⟨ \frac{\partial}{\partial x_{i}} [\frac{\partial C_{α}}{\partial x_{i}} δ (C (x, t) - c)] - \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial}{\partial x_{i}} δ (C (x, t) - c) ⟩ \\ = - D ⟨ \frac{\partial}{\partial x_{i}} [\frac{\partial C_{α}}{\partial x_{i}} \frac{\partial}{\partial c_{α}} δ (C (x, t) - c)] - \frac{\partial}{\partial c_{α}} \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial}{\partial x_{i}} δ (C (x, t) - c) ⟩ \\ = D \frac{\partial^{2}}{\partial x_{i} \partial x_{i}} ⟨ δ (C (x, t) - c) ⟩ - \frac{\partial}{\partial c_{α}} ⟨ D \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial C_{β}}{\partial x_{i}} \frac{\partial}{\partial c_{β}} δ (C (x, t) - c) ⟩ \\ = D \frac{\partial^{2}}{\partial x_{i} \partial x_{i}} f (c; x, t) - \frac{\partial^{2}}{\partial c_{α} \partial c_{β}} [⟨ D \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial C_{β}}{\partial x_{i}} | c ⟩ f (c; x, t)] . \end{aligned}

Inserting (29) into (11) and using the gradient-diffusion closure (19) and the definition (20) of the drift coefficients one obtains the evolution equation (2) for the Eulerian

PDF f (c; x, t)

, with upscaled diffusion coefficients defined by

\begin{matrix} (30) & D_{i j} = D + D_{i, j}^{*} \end{matrix}

and conditional dissipation tensor defined, according to the second term of the decomposition (29), by

\begin{matrix} (31) & M_{α β} = ⟨ D \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial C_{β}}{\partial x_{i}} | c ⟩ . \end{matrix}

We have thus the following corollary of Propositions 1 and 2.
Corollary 2. Under the assumptions of Proposition 1 and for weighting factor

Θ \equiv 1

, the Eulerian PDF

f (c; x, t)

verifies Eq. (2).

The PDF evolution equation (2) has the form of a Fokker-Planck equation and, according to (22), it is verified by the conditional PDF

p (c ∣ x, t)

of the system of Itô equations (3)-(4) for uniform position PDF

p_{x} (x, t)

. Corollary 2 also applies if as weighting factor

Θ

one chooses

ρ =

const , which corresponds to the case of homogeneous turbulence in classical PDF methods [1, p. 140].

For variable

Θ

and

F (c; x, t) = Θ (c) f (c; x, t)

, one obtains a Fokker-Planck equation of form (2) but with additional source terms. A computation similar to the proof of Proposition 2 shows that these terms depend on partial derivatives of

Θ

with respect to

c_{α}

. The equation is thus unclosed and it is no longer equivalent to the system of Itô equations (3) and (4). However, in classical PDF methods Eq. (2), with diffusion coefficients (30) and mixing coefficients (31), is assumed valid in case of variable density, without further justification [12,13], or in some cases, e.g. [14], by referring to its derivation for constant density flows [15].

Summarizing this section, we note that the exact Fokker-Planck equation has the from (11) and has been derived for weighting functions

Θ

which verify the continuity equation (10). This condition can be relaxed for advection-reaction problems when, by Corollary 1, one obtains the Fokker-Planck equation (21), without mixing term. With the gradient-diffusion closure (19), diffusion in physical space may be introduced in both equations, yet with incomplete upscaled diffusion coefficient

D

, which does not contain the molecular diffusion coefficient

D

. The general form of the mixing coefficients in (11) is that given by the conditional average of the diffusion flux from the transport equation (1). The evolution equation (2) for the Eulerian PDF

f

(

c; x, t

), with diffusion in physical space described by the full upscaled coefficients (30) and mixing modeled as a diffusion in concentration space with coefficients (31) given by the conditional dissipation rate, can be derived only in case of constant

Θ

, e.g., constant density

ρ

as in turbulence studies, or as an approximation under ergodicity assumptions.

4. GRW solutions to PDF equations

The solutions of the PDF evolution equations, defined in the state space

Ω = Ω_{x} \times Ω_{c}

, depend on

3 + N_{α}

independent variables plus the time variable. Numerical methods for PDF equations are thus confronted with the curse of dimensionality. Therefore, instead of using finite difference/element schemes, solution methods for PDF equations are based on equivalent representations as systems of computational particles. Particle methods also avoid the increase of computational errors by numerical diffusion [16], unavoidable in classical schemes for PDF equations, which usually result in strongly advectiondominated problems. For either Pope's particles approach or particle tracking solutions to the system of Itô equations (3)-(4), the computational cost increases linearly with the number of particles. The GRW algorithm [17] described in the following is insensitive to the increase of the number of particles. For advection-diffusion problems with constant coefficients GRW is equivalent to a finite difference scheme [4] but in case of variable coefficients it is faster and free of numerical diffusion. For instance, when solving flow and transport problems in highly heterogeneous media by finite element schemes, if the transport solver is replaced by a GRW scheme the numerical diffusion is completely removed and the total computation time is about 20 times smaller [18]. As compared with sequential procedures, we note that for a grid-free particle tracking scheme the computing time is of the order of the number of particles and for a GRW solution to the same problem and the same number of time steps it is of the order of grid points occupied by particles. For the solution of the PDF problem solved in the next section, using

\sim 10^{24}

particles which move on a lattice with

\sim 10^{5}

points, the GRW computing time is about 0.5 s while a sequential particle tracking would require a computing time

10^{19}

times larger.

To illustrate the GRW-PDF approach we consider a two dimensional PDF problem for joint concentration-position PDF

p (c, x, t)

, solution of the particular form of (2)

\begin{matrix} (32) & \partial_{t} p + V \frac{\partial p}{\partial x} + V_{c} \frac{\partial p}{\partial c} = D \frac{\partial^{2} p}{\partial x^{2}} + D_{c} \frac{\partial^{2} p}{\partial c^{2}} \end{matrix}

The solution of the Fokker-Planck equation (32) is approximated by the point-density at lattice sites of a large number of computational particles evolving according to Itô equations of the form (3)-(4). At a given time step, the computational particles from a lattice site (blue column in Fig. 1) are globally scattered in groups of particles which remain at the position determined by the drift coefficients and particles undergoing diffusive jumps (red columns). The numbers of particles in each group are binomial random variables with parameters determined by the coefficients of the Itô equations, the time step, and the lattice constants. The GRW algorithm is thus a superposition of many weak Euler schemes for systems of Itô equations projected on a regular lattice [4,17].

The two-dimensional GRW illustrated in Fig. 1 is defined by the relations

\begin{aligned} n (i, j, k) = & δ n (i + v_{i}, j + v_{j} ∣ i, j, k) + δ n (i + v_{i} + d_{i}, j + v_{j} ∣ i, j, k) + δ n (i + v_{i} - d_{i}, j + v_{j} ∣ i, j, k) \\ (33) & + δ n (i + v_{i}, j + v_{j} + d_{j} ∣ i, j, k) + δ n (i + v_{i}, j + v_{j} - d_{j} ∣ i, j, k) \\ (34) & n (l, m, k + 1) = & δ n (l, m, k) + \sum_{i \neq l, j \neq m} δ n (l, m ∣ i, j, k) \end{aligned}

where

n (i, j, k)

is the number of particles at the site (

x, c

) = (i

δ x, j δ c

) at the time

t = k δ t

and the

δ n

are binomial random variables describing the spread of

n (i, j, k)

. To the drift and diffusion coefficients of the transport problem,

V (i δ x, k δ t)

D (i δ x, k δ t), V_{c} (i δ x, j δ c, k δ t)

, and

D_{c} (i δ x, j δ c, k δ t)

, one associates dimensionless parameters

\begin{matrix} (35) & v_{i} = V \frac{δ t}{δ x}, v_{j} = V_{c} \frac{δ t}{δ c}, r_{i} = D \frac{2 δ t}{{(d_{i} δ x)}^{2}}, r_{j} = D_{c} \frac{2 δ t}{{(d_{j} δ c)}^{2}} \end{matrix}

where

d_{i}

and

d_{j}

are integers describing the amplitude of the diffusive jumps.

Fig. 1. Two-dimensional GRW algorithm in (x, c) space. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

The average over GRW runs of the terms in (33) are related by

\begin{aligned} \overset{―}{δ n (i + v_{i}, j + v_{j} ∣ i, j, k)} = (1 - r_{i} - r_{j}) \overset{―}{n (i, j, k)}, \\ \overset{―}{δ n (i + v_{i} \pm d_{i}, j + v_{j} ∣ i, j, k)} = \frac{r_{i}}{2} \overset{―}{n (i, j, k)}, \\ (36) & \overset{―}{δ n (i + v_{i}, j + v_{j} \pm d_{j} ∣ i, j, k)} = \frac{r_{j}}{2} \overset{―}{n (i, j, k)} . \end{aligned}

The major feature of the GRW algorithm is its self-averaging property [17, Fig. 5]. That means that for large enough numbers of computational particles the averages over runs are well approximated by the corresponding quantities computed in a single GRW run, e.g.

\bar{n} ≃ n

and

\overset{―}{δ n} ≃ δ n

An efficient implementation of the GRW algorithm is obtained when the binomial random variables

δ n

are approximated by

\begin{matrix} (37) & δ n (i + v_{i} - d, j + v_{j} ∣ i, j, k) = {\begin{array}{cc} n / 2 & if n is even \\ [n / 2] + θ & if n is odd \end{array} \end{matrix}

where

n = r_{i} \overset{―}{n (i, j, k)}

and

θ

is a random variable taking the values 0 and 1 with probability

1 / 2

. The number of particles jumping in the opposite direction is simply given by

δ n (i + v_{i} + d, j + v_{j} ∣ i, j, k) = n - δ n (i + v_{i} - d, j + v_{j} ∣ i, j, k)

. The same procedure is further applied for the jumps on the

j

-axis. In practice, (37) is implemented by summing up reminders of division by 2 and multiplication of

n (i, j, k)

r_{i}

and

r_{j}

and by assigning a particle to the lattice site where the sum of reminders reaches the unity. In this way, one avoids the need to use a random number generator to compute

θ

[18].

The GRW algorithm is free of numerical diffusion, because the diffusive jumps and the nominal diffusion coefficients are related by the last two equations (35), and, when implemented by (37), it is practically insensitive to the increase of the number of particles (see

[4, 17, 18]

for implementation details and convergence tests). Since the particles move between lattice sites on which mean values are also defined through particle densities, the GRW-PDF solution presented in the next section avoids the artificial diffusion generated in classical PDF methods by interpolation to particle positions of the means computed by averaging over computational cells.

5. Application to transport in groundwater

Monte Carlo simulations of two-dimensional passive transport of a single chemical species in saturated aquifers were used to estimate the

PDF f (c; x, t)

of the cross-section space-average concentration

C (x, t) = \int_{0}^{L_{y}} c (x, y, t) d y

, where

L_{y}

is the transverse dimension of the computational domain, estimated at the

x

-coordinate of the plume center of mass,

x = ⟨ V ⟩ t [11]

The upscaled drift coefficient

V

in Eq. (32) is the ensemble mean velocity, equal to the velocity of the center of mass

⟨ V ⟩ = 1 m /

day. The upscaled diffusion coefficient

D

is the longitudinal component of the time dependent "ensemble dispersion coefficient". The latter is a self-averaging quantity for transport in random velocity fields with finite correlation range considered here. Hence,

D (t)

was efficiently determined by using a single trajectory of diffusion in a single realization of the random velocity field [4, Section 7.2].

Two different mixing models

M

were considered in the concentration Itô equation (4). The first one is a general mixing model written as a sum of the interaction by exchange with the mean model (IEM) [1], a drift term given by the attenuation of the mean concentration due to the local diffusion [19], and an additional noise term modeled as a Wiener process, needed to control the shape of the PDF [3], written on the right hand side of Eq. (38),

\begin{matrix} (38) & d C (t) = - a (t) [C (t) - ⟨ C (t) ⟩] d t + Δ ⟨ C (t) ⟩ + b d W (t) . \end{matrix}

The coefficient

a (t)

was estimated, similarly to turbulence problems [1,12], by the inverse of diffusion time scale

D (t) / λ^{2}

, with

λ = 1 m

, the characteristic correlation scale of the transport problem [11].

Δ ⟨ C (t) ⟩

was estimated from a fractional

Fig. 2. Distribution of particles in the (x, c) GRW lattice at successive times,

n (x, c; t = 0, 10, 50, 100)

Fig. 3.

⟨ C ⟩ (x)

at fixed

t = 10, 50

, and 100 days (peaks), and

⟨ C ⟩ (t)

at the plume center of mass

x = V t

(monotone curves) compared with Monte Carlo results. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

step, performed at each time step of the GRW-PDF simulation, consisting of a GRW simulation of the one-dimensional diffusion with the constant local diffusion coefficient

D = 0.01 m /

day considered in the Monte Carlo simulations. The amplitude

b \approx 10^{- 6}

day

^{- 1}

of the Wiener process

W (t)

was adjusted by comparisons with the Monte Carlo inferred Eulerian

PDF f (c; x, t)

. The results obtained by this extended IEM model are represented by green curves in Figs. 3, 4, and 5.

The second model consisted of a drift term equal to the rate of decay of the mean concentration at the center of mass,

⟨ C (x, t) ⟩

, determined from the ensemble of Monte Carlo simulations [11], and a small diffusion in the concentration space. The corresponding diffusion coefficient starts with an initial value adjusted in the same way as the noise term in (38) and decays exponentially in time, as suggested by a preliminary analysis of concentration time series generated during the Monte Carlo simulations, carried out with an automatic algorithm [20].

The Eulerian concentration

PDF f (c; x, t)

was simulated by the GRW algorithm and compared with the one obtained by Monte Carlo simulations. Since the cross-section concentration

C (x, t)

is ergodic with a good approximation [11],

C (x, t) ≃ ⟨ C (x, t) ⟩

and

f (c; x, t)

is well approximated by the conditional PDF

p (c ∣ x, t)

, according to (22). The latter was determined by the concentration-position PDF

p (c, x, t)

, approximated by the GRW solution of the Fokker-Planck equation (32), divided by the position

PDF p_{x} (x, t)

. The initial distribution of particles in the (

x, c

) plane was approximated by multiplying the Monte Carlo PDF at

t = 1

day by

10^{24}

particles. Fig. 2 shows the evolution of the computational particles during the GRW simulation.

The position PDF

p_{x} (x, t)

, obtained by integrating over the

c

-coordinate the joint PDF

p (c, x, t)

, defines the ensemble average concentration

⟨ C ⟩ (x, t)

, according to (6). Fig. 3 shows a good agreement, for both mixing models, between Monte Carlo results and GRW-PDF simulations. This result is already expected because

⟨ C ⟩ (x, t)

is the probability density of the computational particles governed by the Itô equation (3), which is independent of the concentration Itô equation (4) and of the mixing model. The comparison presented in Figs. 4 and 5 shows that the mixing model

M

based on the rate of decay of "measured" (i.e. simulated) concentration resembles quite well the Monte Carlo results. The extended IEM model (38) instead fails to capture the transport of the PDF in physical and concentration spaces. The discrepancy may be attributed to structural differences between groundwater and turbulent flows. Such models are better suited for homogeneous systems [1,19], the IEM model being exact in case of statistically homogeneous Gaussian concentrations [2, p. 551].

Fig. 4. Transported

PDF f (c; x, t)

at the plume center of mass

x = V t

. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Fig. 5. GRW and Monte Carlo cumulative distribution functions

c d f (c; x, t), x = V t, t = 0, 10, 30, 50

, and 100 days (from right to left). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

6. Conclusions

We establish relations between Eulerian random fields describing concentrations of chemical species transported in random environments and Itô processes in physical and concentration spaces by equating the PDF of the position Itô process to the ensemble average of a conserved scalar specified by a combination species concentrations. Based on this correspondence we show that the Fokker-Planck equation for the concentration-position PDF of the Itô process coincides with the Eulerian evolution equation of the concentration PDF weighted by the conserved scalar, provided that the latter obeys a continuity equation. For a constant weighting factor or under the ergodic assumption of small randomness of the conserved scalar, the Fokker-Planck equation takes the form of a diffusion in the (

x, c

) space.

In classical PDF methods for turbulent reacting flows the weighting function is the fluid density. This implies a uniform position PDF for constant density flows which may be inconsistent with the physical problem. To overcome this difficulty, numerical solutions to PDF equations are obtained by tracking in the physical space notional particles with an associated composition of species concentrations, which is updated by solving concentration evolution equations for each particle. This approach, mainly useful in problems for large reaction systems, is affected by the increase of the computational cost with the number of particles and by the numerical diffusion induced by interpolations of cell-averages to particle positions.

The alternative GRW-PDF approach proposed in this paper provides solutions to Fokker-Planck equations for realistic initial conditions represented by particles distributions in the (

x, c

) space. The Eulerian PDF weighted by a conserved scalar is estimated by dividing the joint concentration-position PDF by the position PDF of the system of computational particles. The procedure is insensitive to the increase of the number of particles and is strictly free of numerical diffusion because it does not involve cell-averages and interpolation procedures. The computational effort in PDF simulations for multi-component reactive transport could be reduced by parallelization with respect to the number of reactants.

Our numerical illustration for transport in groundwater shows that the mixing model consisting of the rate of decay of the mean concentration and an additive noise performs much better than the IEM model used in turbulence studies. Some tests done with mixing models inferred from single-realizations of Monte Carlo simulated concentrations also yield fairly good GRW-PDF solutions. Therefore, detrending time series of measured concentrations by efficient methods [21,20] seems to be a promising approach to model mixing in practical applications.

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft-Germany under Grants AT 102/7-1, KN 229/15-1, SU 415/2-1 and by the Jülich Supercomputing Centre-Germany in the frame of the Project JICG41.

References

[1] S.B. Pope, PDF methods for turbulent reactive flows, Prog. Energy Combust. Sci. 11 (2) (1985) 119-192.
[2] S.B. Pope, Turbulent Flows, Cambridge University Press, Cambridge, 2000.
[3] R.O. Fox, Computational Models for Turbulent Reacting Flows, Cambridge University Press, New York, 2003.
[4] N. Suciu, Diffusion in random velocity fields with applications to contaminant transport in groundwater, Adv. Water Resour. 69 (2014) 114-133.
[5] P.E. Kloeden, E. Platen, Numerical Solutions of Stochastic Differential Equations, Springer, Berlin, 1999.
[6] H. Wang, P.P. Popov, S.B. Pope, Weak second-order splitting schemes for Lagrangian Monte Carlo particle methods for the composition PDF/FDF transport equations, J. Comput. Phys. 229 (2010) 1852-1878.
[7] D.W. Meyer, P. Jenny, H.A. Tchelepi, A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media, Water Resour. Res. 46 (2010) W12522.
[8] A.Y. Klimenko, On simulating scalar transport by mixing between Lagrangian particles, Phys. Fluids 19 (2007) 031702.
[9] D.C. Haworth, Progress in probability density function methods for turbulent reacting flows, Prog. Energy Combust. Sci. 36 (2010) 168-259.
[10] R.W. Bilger, The structure of diffusion flames, Combust. Sci. Technol. 13 (1976) 155-170.
[11] N. Suciu, C. Vamos, J. Vanderborght, H. Hardelauf, H. Vereecken, Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res. 42 (2006) W04409.
[12] F.A. Jaberi, P.J. Colucci, S. James, P. Givi, S.B. Pope, Filtered mass density function for large-eddy simulation of turbulent reacting flows, J. Fluid Mech. 401 (1999) 85-121.
[13] V. Raman, H. Pitsch, R.O. Fox, Hybrid large-eddy simulation/Lagrangian filtered-density-function approach for simulating turbulent combustion, Combust. Flame 143 (2005) 56-78.
[14] S.B. Pope, The statistical theory of turbulent flames, Philos. Trans. R. Soc. Lond. Ser. A 291 (1979) 529-568.
[15] S.B. Pope, The probability approach to the modelling of turbulent reacting flows, Combust. Flame 27 (1976) 299-312.
[16] F.A. Radu, N. Suciu, J. Hoffmann, A. Vogel, O. Kolditz, C.-H. Park, S. Attinger, Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study, Adv. Water Resour. 34 (2011) 47-61.
[17] C. Vamoş, N. Suciu, H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comput. Phys. 186 (2003) 527-544.
[18] N. Suciu, F.A. Radu, A. Prechtel, F. Brunner, P. Knabner, A coupled finite element-global random walk approach to advection-dominated transport in porous media with random hydraulic conductivity, J. Comput. Appl. Math. 246 (2013) 27-37.
[19] R. McDermott, S.B. Pope, A particle formulation for treating differential diffusion in filtered density models, J. Comput. Phys. 226 (2007) 947-993.
[20] C. Vamoş, M. Crăciun, Automatic Trend Estimation, Springer, Dortrecht, 2012.
[21] C. Vamoş, M. Crăciun, Separation of components from a scale mixture of Gaussian white noises, Phys. Rev. E 81 (2010) 051125.

- Corresponding author at: Mathematics Department, Friedrich-Alexander University of Erlangen-Nuremberg, Cauerstraße 11, 91058 Erlangen, Germany.
E-mail address: suciu@math.fau.de (N. Suciu).

2015

Some properties of solutions of the homogeneous nonlinear second order differential equations

January 16, 2019

Conjugate free convection in a vertical channel filled with nanofluid

December 3, 2018

Sandwich theorems for radiant functions

April 27, 2023

A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media

Abstract

Authors

Keywords

Cite this paper as:

PDF

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

Google Scholar Profile

References

Paper in HTML form

A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media

ARTICLE INFO

Abstract

1. Introduction

2. Classical PDF methods

3. The Fokker-Planck approach

3.1. Correspondence between random fields and stochastic processes

3.2. Derivation of the Fokker-Planck equation

3.3. The treatment of the conditional diffusion term in Eq. (11)

4. GRW solutions to PDF equations

5. Application to transport in groundwater

6. Conclusions

Acknowledgments

References

Related Posts