We consider an array of species concentrations $\mathbf{C}(\mathbf{x},t)=\left\{C_{\alpha}(\mathbf{x},t)\right\}$ related by reaction terms $\mathbf{S}(\mathbf{C})=\left\{S_{\alpha}(\mathbf{C})\right\},\alpha=1,\ldots,N_{\alpha}$, where $N_{\alpha}$ is the number of chemical species. For constant diffusion coefficient $D$, reaction-diffusion processes in statistically homogeneous random velocity fields $\mathbf{V}$ with divergencefree samples, are governed by the system of $N_{\alpha}$ stochastic balance equations

The one-point one-time $\operatorname{PDF}f(\mathbf{c};\mathbf{x},t)$ of the random concentrations $\mathbf{C}$ solving (1) satisfies

where $\mathcal{V}$ and $\mathcal{D}$ are upscaled drift and diffusion coefficients and $\mathcal{M}$ is the tensor of the average dissipation rate conditional on concentration, which accounts for mixing by molecular diffusion. The most important feature of the PDF methods is that the reaction term $\mathbf{S}$ in (2) is in a closed form, the same as in the balance equation (1) [1].

The numerical solution of the PDF equation (2) is usually obtained by solving a system of Itô equations describing the evolution of an ensemble of computational particles in physical and concentration spaces [2],

where $\mathbf{C}(t)=\mathbf{C}(\mathbf{X}(t),t),\mathbf{W}$ is a Wiener process with $E\{\mathbf{W}(\mathbf{X}(t),t)\}=0$ and $E\left\{\mathbf{W}^{2}(\mathbf{X}(t),t)\right\}=2\int_{0}^{t}\mathcal{D}\left(\mathbf{X}(t),t^{\prime}\right)dt^{\prime}$, and $\mathbf{M}$ is a mixing model for the term containing the dissipation rate $\mathcal{M}_{\alpha\beta}$ in equation (2) [3]. The equations (3) and (4) describe a stochastic process $\left\{X_{i}(t),C_{\alpha}(t)\right\},C_{\alpha}(t)=C_{\alpha}(\mathbf{X}(t),t),t\geq 0$, $i=1,2,3,\alpha=1,\ldots,N_{\alpha}$. Throughout the paper, we follow the convention which denotes random functions by capital letters and their values in the state space by corresponding lowercase letters. At given times $t$, the process takes fixed values $\mathbf{x}=\mathbf{X}(t)$ and $\mathbf{c}=\mathbf{C}(t)$ in the Cartesian product state space $\Omega=\Omega_{x}\times\Omega_{c}$, where $\Omega_{x}$ is the three-dimensional physical space and $\Omega_{c}$ is the $N_{\alpha}$-dimensional concentration space.

The joint concentration-position PDF $p(\mathbf{c},\mathbf{x},t)$ of the process described by the system of Itô equations (3) and (4) solves a Fokker-Planck equation which is in general different from the PDF equation (2). In this paper we introduce consistency conditions relating the statistics of the random field $\mathbf{C}(x,t)$ to that of the stochastic process $\{\mathbf{X}(t),\mathbf{C}(t)\}$ and derive the corresponding FokkerPlanck equation. Further, we propose a new approach to approximate the concentration PDF by the solution of the Fokker-Planck equation, based on a global random walk (GRW) algorithm. Finally, we illustrate the GRW-PDF approach for non-reactive transport in saturated aquifers.

In studies on turbulent reacting flows probability densities are usually defined by expectations of Dirac- $\delta$ functions $[1,2,4,5]$. For instance, the concentration

PDF is given by

where the angular bracket indicates stochastic average with respect to the probability measure of the random concentration vector $\mathbf{C}(\mathbf{x},t)$ and the multidimensional $\delta$ function is defined by the product

The consistency of (5) with stochastic averaging is ensured by the definition of the $\delta$ functional. For instance, the integral weighted by the PDF (5) of a continuous function $Q(\mathbf{c})$ over the concentration space $\Omega_{c}$ yields the expectation of $Q$ as a function of the random concentration $\mathbf{C}(\mathbf{x},t)$ :

To derive the above relation we used the continuity of $Q(\mathbf{c})$, which allows the permutation of the integral with the stochastic average, and the obvious prop$\operatorname{erty}\int_{\Omega_{c}}Q(\mathbf{c})\delta(\mathbf{C}(\mathbf{x},t)-\mathbf{c})d\mathbf{c}=\int_{\mathbb{R}^{N_{\alpha}}}I_{\Omega_{c}}(\mathbf{c})Q(\mathbf{c})\delta(\mathbf{C}(\mathbf{x},t)-\mathbf{c})d\mathbf{c}=Q(\mathbf{C}(\mathbf{x},t))$, where $I_{\Omega_{c}}$ is the indicator function of the set $\Omega_{c}$. The definition (5) can be generalized, by using products of $\delta$ functions, to obtain a consistent hierarchy of multi-point probability distributions which completely define the random concentration $\mathbf{C}(\mathbf{x},t)$ as a random function [3, Appendix A.1].

A straightforward derivation of the evolution equation for the PDF $f(\mathbf{c};\mathbf{x},t)$ is based on the evaluation of the time derivative of its definition (5) through formal manipulations of derivatives of $\delta$ functions [6]. In terms of Dirac distributions, the derivative of the $\delta$ function is defined by the relation

for any smooth function $\varphi$ with compact support in $\mathbb{R}$. The usual notation for the distributional derivative is $\delta^{\prime}[\varphi]=-\delta\left[\varphi^{\prime}\right]=-\varphi^{\prime}\left(y_{0}\right)$. Further, (6) can be generalized for the case where $y_{0}$ is a given value of a composite function of some variable $x,y_{0}=g(x)$. Then, the Dirac functional applied to the test function $\varphi$ reads

The derivative of (7) with respect to $x$ follows as

where the second equality is implied by (7) and the third equality by (6). In the common compact notation for distributions we get

which is often written in applications to turbulent flows as $[1,4,6]$

The correctness of the results obtained with formal manipulations of the relation (8) can be checked by multiplying them with test functions, integrating over $y$ and using (6) (see e.g. [6, Sect. 2.2.1]).

Considering multidimensional $\delta$ functions and the vectorial function $\mathbf{g}=\mathbf{C}(\mathbf{x},t)$ one obtains, similarly to (8), the formal expression of the partial derivative with respect to the time variable,

and the partial derivatives with respect to the spatial variables,

Using (9), the time derivative of the PDF (5) is computed as follows:

Performing the stochastic average on the right hand side of (11) requires more information on the statistics of the random concentration $\mathbf{C}(\mathbf{x},t)$ than the knowledge of the one-point $\operatorname{PDF}f(\mathbf{C};\mathbf{x},t)$. To see that, note that the time derivative is the limit of the increment of $\mathbf{C}$ to the corresponding time increment, which is a two-point quantity. Thus, a two-point PDF is needed if one wants to perform the average before taking the limit. Alternatively, the average can be performed if one knows the joint statistics of $\partial\mathbf{C}/\partial t$ and $\mathbf{C}$. To obtain meaningful expressions of these kind of stochastic averages, we follow the approach of Fox [2] and consider a generic random function $\mathbf{Z}(\mathbf{x},t)$ which is not described by the one-point $\operatorname{PDF}f(\mathbf{C};\mathbf{x},t)$. Similarly to the average occurring in (11), we have, in general and for the random function $F(\mathbf{Z}(\mathbf{x},t)$, the average

where $f(\mathbf{c},\mathbf{z};\mathbf{x},t)=\langle\delta(\mathbf{Z}(\mathbf{x},t)-\mathbf{z})\delta(\mathbf{C}(\mathbf{x},t)-\mathbf{c})\rangle$ defines, similarly to (5), the joint PDF in the ( $\mathbf{c},\mathbf{z}$ ) space and $f(\mathbf{z}\mid\mathbf{c};\mathbf{x},t)=f(\mathbf{c},\mathbf{z};\mathbf{x},t)/f(\mathbf{c};\mathbf{x},t)$ is a conditional PDF. Thus, we finally obtain

which is just the expectation of the random function $F(\mathbf{Z}(\mathbf{x},t))$ conditional on a fixed value of the concentration vector $\mathbf{c}$ multiplied by the one-point PDF $f(\mathbf{c};\mathbf{x},t)$. In the following, we use the shorthand notation $\langle\cdot\mid\mathbf{c}\rangle$ to denote conditional averages $\langle\cdot\mid\mathbf{C}(\mathbf{x},t)=\mathbf{c}\rangle$.

Now, substituting the time derivative $\partial C_{\alpha}(\mathbf{x},t)/\partial t$ from (1) into (11) and using (12), we obtain the evolution equation for the PDF $f(\mathbf{c};\mathbf{x},t)$ :

We note that the last term in (13) is in a closed form because the reaction term $S_{\alpha}$ in (1) is completely determined by the random concentration $\mathbf{C}$, which implies that its conditional expectation is just the value of $S_{\alpha}$ evaluated for
the sample space value $\mathbf{c}$,

The first two terms on the right hand side of (13) are unclosed and require modeling.

The first unclosed term can be transformed as follows:

where the first equality follows from (12) and to obtain the final result we used (10) and the incompressibility condition $\partial V_{i}/\partial x_{i}=0$. Defining the velocity fluctuation $\mathbf{U}=\mathbf{V}-\langle\mathbf{V}\rangle$ and using the gradient diffusion closure $\langle\mathbf{U}\mid\mathbf{c}\rangle f(\mathbf{c};\mathbf{x},t)=-\mathbf{D}^{*}\nabla f(\mathbf{c};\mathbf{x},t)[1,2,5]$, (14) becomes

The upscaled diffusion tensor $D_{i,j}^{*}$ is provided by turbulence models [4] or by stochastic upscaling of diffusion in random velocity fields [3,10].

The second unclosed term in (13) is the conditional expectation of the molecular diffusion term in equation (1). We will show that this term is related to the divergence of the diffusive flux in physical space of the PDF $f(\mathbf{c};\mathbf{x},t)$, with the same diffusion coefficient $D$ as in equation (1). Since by (5) $f(\mathbf{c};\mathbf{x},t)=\langle\delta(\mathbf{C}(\mathbf{x},t)-\mathbf{c})\rangle$, the diffusive flux of $f$ is determined by the expectation of (10) and its divergence can be expressed as follows:

Using (15) and (16) to define the coefficients we get

and (13) takes the form of the PDF equation (2). The same result has been obtained by equating the stochastic average of the operator from the left hand side of (1) applied to a test function with that corresponding to the right hand side of (1) and by integrating by parts [1, 2, 5, 7].

The system of Itô equations (3) and (4), used to numerically approximate the solution of the PDF evolution equation (2) by computational particles evolving in the state space $\Omega=\Omega_{x}\times\Omega_{c}$, describes a stochastic process $\{\mathbf{C}(t),\mathbf{X}(t)\}$. As a mathematical object, this process is a random function indexed by time, with one-time statistics completely described by the joint concentration-position $\mathrm{PDF}p(\mathbf{c},\mathbf{x},t)$. On the other side, the random concentration $\mathbf{C}(\mathbf{x},t)$ is a random function with four indices, the three spatial coordinates and the time, and its one-point (in time and space) PDF is the solution $f(\mathbf{c};\mathbf{x},t)$ of the evolution
equation (2). Consequently, the two PDFs verify different normalization conditions, $\int_{\Omega}p(\mathbf{c},\mathbf{x},t)d\mathbf{c}d\mathbf{x}=1$ and $\int_{\Omega_{c}}f(\mathbf{c};\mathbf{x},t)d\mathbf{c}=1$, respectively. While the parameter $\mathbf{x}$ of the random concentration is a point in the physical space, $\mathbf{X}(t)$ is a stochastic process described by (3), with the one-time PDF $p_{x}(\mathbf{x},t)$, which is a solution of the associated Fokker-Planck equation (equation (2) with the right hand side set to zero). Since this equation is equivalent to an advectiondiffusion equation describing the transport of a passive scalar [3], the position PDF $p_{x}(\mathbf{x},t)$ has the meaning of a normalized concentration. Since (3) describes an upscaled transport process, it follows that $p_{x}(\mathbf{x},t)$ represents the corresponding mean concentration. A consistency condition, relating $p_{x}(\mathbf{x},t)$ to the statistics of the random concentration $\mathbf{C}(\mathbf{x},t)$, is the starting point in designing numerical solutions to the PDF equation (2) based on the solution of the system of Itô equations (3) and (4).

In case of a single chemical species $c$ which is conserved, the statistics of the Itô process is obviously consistent with that of the random concentration if and only if

Rewritten with the aid of the corresponding PDFs, the relation (20) becomes

A sufficient condition for (21) is given by the equality of the integrants,

We note that when particle representations of the concentration PDF $f(\mathbf{c};\mathbf{x},t)$ are used to compute concentration statistics, the mean concentration is, according to (20), the position PDF $p_{x}(\mathbf{x},t)$ of the system of computational particles. Furthermore, the average with respect to $p(c,\mathbf{x},t)$ of the state space variable $c$ yields, according to (21), the expectation of the squared concentration. Thus, subtracting from this expectation the squared position PDF, one obtains the variance of the random concentration,

In case of multicomponent reactive transport, the position PDF $p_{x}(\mathbf{x},t)$ can be expressed as

that is, as an arithmetic average of position PDFs $p_{x}^{(\alpha)}=\int_{\Omega_{c_{\alpha}}}p_{c_{\alpha},x}\left(c_{\alpha},\mathbf{x},t\right)dc_{\alpha}$ associated to each chemical species, where $p_{c_{\alpha},x}\left(c_{\alpha},\mathbf{x},t\right)$ are marginal PDFs of $p(\mathbf{c},\mathbf{x},t)$ [7]. The natural conjecture is that the position PDF $p_{x}(\mathbf{x},t)$ corresponds to the expectation of the arithmetic mean of the species concentrations,

where $C^{*}$ is a normalization constant. This conjecture is thus formulated as

Analogous to (22), the choice of (25) implies the following relation between the statistics of the one-dimensional PDFs of the random concentration and that of the Itô process:

The conjecture (25) is strictly valid if the arithmetic average (24) is a conserved quantity solving the balance equation (1) without reaction term,

Then, the expectation $\langle\Theta\rangle$ of a conserved scalar verifies the equation (2) with the right hand terms set to zero [8]. Since this equation coincides with the Fokker-Planck equation associated to the Itô equation (3), verified by $p_{x}(\mathbf{x},t)$, the equality (25) holds true. That this is indeed the case can be proved by a slight extension of the method used by Bilger [9] to construct conserved scalars as concentrations of chemical elements. Let $r_{\alpha k}$ be the weight (e.g. the mass fraction) of the chemical element indexed by $k$ in the composition of the molecules $\alpha$ and let $C_{k}$ be the total concentration of the element $k$. Obviously, the elemental mass sum to unity, $\sum_{k=1}^{N_{k}}r_{\alpha k}=1$, and $C_{k}=\sum_{\alpha=1}^{N_{\alpha}}r_{\alpha k}C_{\alpha}$. It follows that

that is, the sum of elemental concentrations equals the sum of species concentrations. Since elemental concentrations are conserved under chemical reactions, the sum of species concentrations is a conserved scalar. Further, summing up the equations (1) with species independent coefficients, one obtains
the relation $\sum_{\alpha=1}^{N_{\alpha}}S_{\alpha}=0$, which, when $C_{\alpha}$ are mass concentrations, expresses the conservation of mass of the reacting system.

The conjecture (25) is thus true for $\Theta$ defined by the sum of species concentrations or by their arithmetic mean (24), as well as for any other conserved scalar properly normalized. For instance, if the problem is formulated in terms of mass concentrations, $C_{\alpha}=\rho_{\alpha}$, we can choose $\Theta=\sum_{\alpha=1}^{N_{\alpha}}\rho_{\alpha}=\rho$, where $\rho$ is the fluid density. Then, (26) becomes $\rho f(\mathbf{c};\mathbf{x},t)=p(\mathbf{c},\mathbf{x},t)$ and (25) takes the form

where $M=\int_{\Omega}\langle\rho\rangle(\mathbf{x},t)d\mathbf{x}$ is the total mass of fluid in $\Omega_{x}$.
If the species concentrations are defined as mass fractions, $\rho_{\alpha}/\rho$, and their sum is used to define $\Theta=\sum_{\alpha=1}^{N_{\alpha}}\rho_{\alpha}/\rho=1$, the relation (26) becomes $f(\mathbf{c};\mathbf{x},t)=p(\mathbf{c},\mathbf{x},t)$ and (25) implies $p_{x}(\mathbf{x},t)=1$, thus $p(\mathbf{c},\mathbf{x},t)=p(\mathbf{c}\mid\mathbf{x},t)p_{x}(\mathbf{x},t)=p(\mathbf{c}\mid\mathbf{x},t)$ and the concentration PDF equals the conditional PDF of the system of Itô equations (3)-(4), $f(\mathbf{c};\mathbf{x},t)=p(\mathbf{c}\mid\mathbf{x},t)$. This choice, which implies a constant position PDF $p_{x}(\mathbf{x},t)$, is however consistent only if the drift and diffusion coefficients in (2) are constants or if they are subject to some constraints, for instance $\frac{\partial}{\partial x_{i}}\mathcal{V}_{i}=\frac{\partial^{2}}{\partial x_{i}\partial x_{i}}\mathcal{D}$ in case of isotropic $\mathcal{D}$ [4]. Uniform position PDF is also implied by (28) in case of constant density flows, with $p_{x}(\mathbf{x},t)=\rho/M$.

Remark 1. In combustion theory $[1,5]$, one uses a mixture of the two choices presented above: concentrations defined as mass fractions and the relation (28), based on a conserved scalar constructed as a sum of mass concentrations, to connect the statistics of the random concentration to that of the computational particles. This approach is nevertheless consistent and may be used to construct a system of computational particles stochastically equivalent to the PDF evolution equation [1, Sect. 3.4]. The relation (27) presumes that the ensemble of $N_{\alpha}$ species contains both the solvent and the solutes, otherwise the problem is not closed because the density, which is the sum of mass concentrations, is not determined. The density weighted PDF defined by (27) with $\Theta=\rho$ is adequate for combustion problems but it could be inappropriate for dilute solutions, like those transported in groundwater. The scalar $\Theta$ defined by the sum or the arithmetic mean of the species concentration, irrespective of the particular definition chosen for concentrations, would have the advantage that it avoids including the carrying fluid into the reaction system and does not force a uniform position PDF for constant density fluids.

It has been shown [7] that if the weighting function $\Theta$ in (26) obeys the continuity equation

the evolution of the joint concentration-position PDF $p(\mathbf{c},\mathbf{x},t)$ is described by the following Fokker-Planck equation associated with the Itô process (3)-(4)

To derive equation (30) in the $\delta$-function approach from Section 2 we make use of the "shifting" property

which can be readily checked by integrating both sides of (31) with respect to c. To derive the equation for $p(\mathbf{c},\mathbf{x},t)=\Theta(\mathbf{c})f(\mathbf{c};\mathbf{x},t)=\Theta(\mathbf{c})\langle\delta(\mathbf{C}(\mathbf{x},t)-\mathbf{c})\rangle$ we start, as for (11), by computing its time derivative:

In the first equality of (32) we introduced $\Theta(\mathbf{c})$ under average, because it is a non-random function of the state space variable $\mathbf{c}$, the second equality follows from the shifting property (31), and for the first term in the last line we substituted the partial spatial derivative (10) of the $\delta$-function. With the time
derivative of $C_{\alpha}$ from the transport equation (1), this latter term becomes

We used the shifting property (31) and the conditional average (12) to obtain the second term on the right hand side of (33), which is just the last term of the Fokker-Planck equation (30). The first term on the right hand side of (33) can be rewritten by using (10) as

We note that due to the incompressibility condition $\partial V_{i}/\partial x_{i}=0$ and because $\Theta$ obeys the continuity equation (29), the first term in the final expression (34) cancels the last term of (32). Further, using (31), (12), (26), the incompressibility condition, and writing the velocity as a sum of mean and fluctuations, the second term of (34) becomes

The Fokker-Planck equation (30) is finally obtained by recursively substituting (35), (34), and (33) into (32).

Remark 2. A conserved combination of species concentrations $\Theta$ defined by (24) which does not include all the constituents of the fluid satisfies equation (27). The latter reduces to (29) only if $D=0$. Then, the conditional diffusion term on the right hand side of (30) drops out and the resulting equation describes the evolution of the PDF in advection-reaction problems or, as an approximation, in advection-dominated reactive transport. Equation (30) is strictly verified when both the solutes and the solvents are considered and $\Theta$ is defined as the total mass density of the fluid, $\rho$, which verifies the continuity equation (29). Then, according to (26), the solution of the Fokker-Planck
equation (30) is $p(\mathbf{c},\mathbf{x},t)=\rho(\mathbf{c})f(\mathbf{c};\mathbf{x},t)=\mathcal{F}(\mathbf{c},\mathbf{x},t)$, which defines the mass density function (e.g. [1, Eq. (3.59)]).

Remark 3. It is also readily to check, by using (16), that the FokkerPlanck equation (30) takes the form of the PDF equation (2) only if $\Theta$ is a constant (see e.g. [7]), which implies a uniform position PDF $p_{x}(\mathbf{x},t)=\langle\Theta\rangle(\mathbf{x},t)$.

Expressing $p(\mathbf{c},\mathbf{x},t)$ as a product of conditional PDF $p(\mathbf{c}\mid\mathbf{x},t)$ and position $\operatorname{PDF}p_{x}(\mathbf{x},t),p(\mathbf{c},\mathbf{x},t)=p(\mathbf{c}\mid\mathbf{x},t)p_{x}(\mathbf{x},t)$, and using (25) and (26) one obtains

If $\Theta=\rho$, then the conditional PDF of the system of computational particles is a discrete representation of the density weighted PDF, $p(\mathbf{c}\mid\mathbf{x},t)=\rho(\mathbf{c})f(\mathbf{c};\mathbf{x},t)/\langle\rho\rangle$.