Towards a filtered density function approach for reactive transport in groundwater

7 years ago

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

Abstract

Evolution equations for probability density functions (PDFs) and filtered density functions (FDFs) of random species concentrations weighted by conserved scalars are formulated as Fokker–Planck equations describing stochastically equivalent processes in concentration-position spaces. This approach provides consistent numerical PDF/FDF solutions, given by the density in the concentration-position space of an ensemble of computational particles governed by the associated Itô equations. The solutions are obtained by a global random walk (GRW) algorithm, which is stable, free of numerical diffusion, and practically insensitive to the increase of the number of particles. The general FDF approach and the GRW numerical solution are illustrated for a reduced complexity problem consisting of the transport of a single scalar in groundwater. Randomness is induced by the stochastic parameterization of the hydraulic conductivity, characterized by short range correlations and small variance. The objective is to infer the statistics of the random concentration sampled at the plume center of mass, integrated over the transverse dimension of a two-dimensional spatial domain. The PDF/FDF problem can therefore be formulated in a two-dimensional domain as well, a spatial dimension and one in the concentration space. The upscaled drift and diffusion coefficients describing the PDF transport in the physical space are estimated on single-trajectories of diffusion in velocity fields with short-range correlations, owing to their self-averaging property. The mixing coefficients describing the PDF transport in concentration spaces are parameterized by the trend and the noise inferred from the statistical analysis of an ensemble of simulated concentration time series, as well as by classical mixing models. A Gaussian spatial filter applied to a Kraichnan velocity field generator is used to construct coarse-grained simulations (CGS) for FDF problems. The purposes of the CGS simulations are two-fold: first to understand the significance of the FDF approach from a practical point of view and its relation to the PDF approach; second to investigate the limits of the mixing models considered here and the desirable features of the mixing models for groundwater systems.

Authors

N. Suciu
Tiberiu Popoviciu Institute of Numerical Analysis

L. Schüler

S. Attinger

P. Knabner

Keywords

PDF/FDF methods; Mixing; Global random walk; Groundwater

Cite this paper as:

N. Suciu, L. Schüler, S. Attinger, P. Knabner (2016), Towards a filtered density function approach for reactive transport in groundwater, Adv. Water Resour., 90, 83–98,
doi: 10.1016/j.advwatres.2016.02.016

PDF

ResearchGate.

About this paper

Journal

Advances in Water Resources

Publisher Name

Elsevier

DOI

10.1016/j.advwatres.2016.02.016

Print ISSN

0309-1708

Online ISSN

1573-1634

MR

ZBL

Google Scholar Profile

[1] Pope SB., The probability approach to the modelling of turbulent reacting ﬂows. Combust Flame 1976; 27:299–312.
CrossRef (DOI)

[2] Pope SB., PDF methods for turbulent reactive ﬂows. Prog Energy Combust Sci 1985;11(2):119–192.
CrossRef (DOI)

[3] Pope SB., Turbulent Flows. Cambridge: Cambridge University Press: 2000.

[4] Fox RO., Computational Models for Turbulent Reacting Flows. New York: Cambridge University Press; 2003.

[5] Haworth DC., Progress in probability density function methods for turbulent reacting ﬂows. Prog Energy Combust Sci 2010; 36:168-259.
CrossRef (DOI)

[6] Haworth DC, Pope SB., Transported probability density function methods for Reynolds-averaged and large-eddy simulations. In: EchekkiT, Mastorakos E, editors. Turbulent combustion modeling. Fluid mechanics and its applications, vol. 95. Dordrecht: Springer; 2011. p.119–42.
CrossRef (DOI)

[7] Colucci PJ, Jaberi FA, Givi P., Filtered density function for large eddy simulation of turbulent reacting ﬂows. PhysFluids 1998; 10(2):499–515.
CrossRef (DOI)

[8] Jaberi FA, Colucci PJ, James S, Givi P, Pope SB., Filtered mass density function for large-eddy simulation of turbulent reacting ﬂows. J. Fluid Mech. 1999; 401:85–121.
CrossRef (DOI)

[9] McDermott R, Pope SB., A particle formulation for treating differential diffusion in filtered density function methods. J Comput Phys 2007; 226:947–993.
CrossRef (DOI)

[10] Heinz S., Uniﬁed turbulence models for LES and RANS, FDF and PDF simulations. Theor Comput Fluid Dyn 2007;21:99118.
CrossRef (DOI)

[11] Dodoulas IA, Navarro-Martinez S., Large eddy simulation of premixed turbulent ﬂames unsing probability density approach. Flow Tur-bulence Combust 2013;90:645–678.
CrossRef (DOI)

[12] Schwede RL, Cirpka OA, Nowak W, Neuweiler I., Impact of sampling volume on the probability density function of steady state concentration. Water Resour Res 2008;44(12):W12433.
CrossRef (DOI)

[13] Sanchez-Vila X, Guadagnini A, Fernandez-Garcia D., Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers. Math Geosci 2009;41:32351.
CrossRef (DOI)

[14] Dentz M, Tartakovsky DM., Probability density functions for passive scalars dispersed in random velocity ﬁelds. Geophys Res Lett 2010;37:L24406.
CrossRef (DOI)

[15] Meyer DW, Jenny P, Tchelepi HA., A joint velocity-concentration PDF method for tracer ﬂow in heterogeneous porous media. Water Resour Res 2010;46:W12522.
CrossRef (DOI)

[16] Cirpka OA, de Barros FPJ, Chiogna G, Nowak W., Probability density function of steady state concentration in two-dimensional heterogeneous porous media. Water Resour Res 2011;47:W11523.
CrossRef (DOI)

[17] Venturi D, Tartakovsky DM, Tartakovsky AM, Karniadakis GE., Exact PDF equations and closure approximations for advective reactive transport. J Comput Phys 2013;243:32343.
CrossRef (DOI)

[18] Suciu N., Diﬀusion in random velocity ﬁelds with applications to contaminant transport in groundwater. Adv Water Resour 2014;69:114–133.
CrossRef (DOI)

[19] Suciu N, Radu FA, Attinger S, Schuler L, Knabner P., A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media. J Comput Appl Math 2015;289:241–252.
CrossRef (DOI)

[20] Suciu N, Schuler L, Attinger S, Vamos¸ C, Knabner P., Consistency issues in PDF methods. An Sti U Ovid Co-Mat 2015;23(3):187–208.
CrossRef (DOI)

[21] Beckie R, Aldama AA, Wood EF., Modeling the large-scale dynamics of saturated groundwater ﬂow using spatial ﬁltering theory: 1.Theoretical development. Water Resour Res 1996;32(5):1269–1280.
CrossRef (DOI)

[22] Beckie, R, Aldama AA, Wood EF., Modeling the large-scale dynamics of saturated groundwater ﬂow using spatial ﬁltering theory: 2.Numerical Evaluation. Water Resour Res 1996;32(5):1281–1288.
CrossRef (DOI)

[23] Efendiev YR, Durlofsky LJ, Lee SH., Modeling of subgrid eﬀects in coarse scale simulations of transport in heterogeneous porous media.Water Resour Res 2000;36:2031–2041.
CrossRef (DOI)

[24] Efendiev Y, Durlofsky LJ., A Generalized convection-diﬀusion model for subgrid transport in porous media. Multiscale Model Simul 2003;1(3):504–526.
CrossRef (DOI)

[25] Attinger S., Generalized Coarse Graining Procedures for Flow in Porous Media. Computational Geosciences 2003;7(4):253–273.
CrossRef (DOI)

[26] Heße F, Radu FA, Thullner M, Attinger S., Upscaling of the advection diﬀusion reaction equation with Monod reaction. Adv Water Resour2009;32:1336–1351.
CrossRef (DOI)

[27] Minier JP, Peirano E., The PDF approach to turbulent and polydispersed two-phase ﬂows. Phys Rep 2001;352:1–214.

[28] Klimenko AY, Bilger RW., Conditional moment closure for turbulent combustion. Progr Energ Combust Sci 1999; 25:595–687.
CrossRef (DOI)

[29] Morales-Casique E, Neuman SP, Gaudagnini A., Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: theoretical framework. Adv Water Resour 2006;29:1238-1255.
CrossRef (DOI)

[30] Morales-Casique E, Neuman SP, Gaudagnini A., Nonlocal and localized analyses of nonreactive solute transport inbounded randomly heterogeneous porous media: computational analysis. Adv Water Resour 2006;29:1399–1418.
CrossRef (DOI)

[31] Pope SB., A Monte Carlo method for the PDF equations of turbulent reactive ﬂow. Combustion Science and Technology 1981;25:159–174.
CrossRef (DOI)

[32] Mobus H, Gerlinger P, Brggemann D., Comparison of Eulerian and Lagrangian Monte Carlo PDF methods for turbulent diﬀusion ﬂames. Combust Flame 2001;124:519–534.
CrossRef (DOI)

[33] Jones WP, Marquis AJ, Prasad VN., LES of a turbulent premixed swirl burner using the Eulerian stochastic ﬁeld method. Combust Flame 2012;159:3079–3095.
CrossRef (DOI)

[34] Valino L., A Field Monte Carlo Formulation for Calculating the Probability Density Function of a Single Scalar in a Turbulent Flow. Flow Turb Combust 1998;60(2):157–172.
CrossRef (DOI)

[35] Sabelnikov V, Soulard O., Rapidly decorrelating velocity-ﬁeld model as a tool for solving one-point Fokker-Planck equations for probability density functions of turbulent reactive scalars. Phys Rev E 2005;72(1):016301.
CrossRef (DOI)

[36] Waclawczyk M, Pozorski J, Minier JP., New Molecular Transport Model for FDF/LES of Turbulence with Passive Scalar. Flow Turbulence Combust 2008;81:235–260.
CrossRef (DOI)

[37] Wang H, Popov PP, Pope SB., Weak second-order splitting schemes for Lagrangian Monte Carlo particle methods for the composition PDF/FDF transport equations. J Comput Phys 2010;229:1852–1878.
CrossRef (DOI)

[38] Kloeden PE, Platen E., Numerical solutions of stochastic diﬀerential equations. Berlin: Springer; 1999.

[39] Herz M., Mathematical modeling and analysis of electrolyte solutions. 2014; PhD thesis. http://www.mso.math.fau.de/ﬁleadmin/am1/projects/PhD Herz.pdf.

[40] Rubin Y, Sun A, Maxwell R, Bellin A., The concept of block eﬀective macrodispersivity and a uniﬁed approach for grid-scale and plume-scale-dependent transport. J Fluid Mech 1999;395:161–180.
CrossRef (DOI)

[41] de Barros FPJ, Rubin Y., Modelling of block-scale macrodispersion as a random function. J Fluid Mech 2011;676:514–545.
CrossRef (DOI)

[42] Bellin A, Tonina D., Probability density function of non-reactive solute concentration in heterogeneous porous formations, J. Contam. Hydrol. 2007;94:109–125.
CrossRef (DOI)

[43] Minier J-P, Chibbaro S, Pope SB., Guidelines for the formulation of Lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent ﬂows. Phys Fluids 2014;26:113303.
CrossRef (DOI)

[44] Klimenko AY., On simulating scalar transport by mixing between Lagrangian particles. Phys Fluids 2007;19:031702.
CrossRef (DOI)

[45] Vamos¸ C, Suciu N, Vereecken H., Generalized random walk algorithm for the numerical modeling of complex diﬀusion processes. J. Comput Phys 2003;186(2):52744.
CrossRef (DOI)

[46] Suciu N, Radu FA, Prechtel A, Brunner F, Knabner P., A coupled ﬁnite element-global random walk approach toadvection-dominated transport in porous media with random hydraulic conductivity. J Comput Appl Math 2013;24627–37.
CrossRef (DOI)

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

[48] Kraichnan RH., Diﬀusion by a random velocity ﬁeld. Phys Fluids 1970;13(1):2231.
CrossRef (DOI)

[49] Vamos¸ C, Craciun M., Separation of components from a scale mixture of Gaussian white noises. Phys Rev E 2010;81:051125.
CrossRef (DOI)

[50] Suciu N, Vamos C., Ergodic estimations of upscaled coeﬃcients for diﬀusion in random velocity ﬁelds. In: L’Ecuyer Pierre, Owen Art20 B, editors. Monte Carlo and quasi-Monte Carlo methods 2008. Berlin: Springer; 2009. p. 617-626.
CrossRef (DOI)

[51] Dagan G., Upscaling of dispersion coeﬃcients in transport through heterogeneous porous formations. In: Peters A et al., editors. Computational Methods in Water Resources X. Norwell, Mass: Kluwer Acad; 1994. p. 431–439.

[52] Schwarze H, Jaekel U, Vereecken H., Estimation of Macrodispersion by Diﬀerent Approximation Methods for Flow and Transport in Randomly Heterogeneous Media. Transport Porous Med 2001;43:265–287.
CrossRef (DOI)

[53] Dentz M, Kinzelbach H,Attinger S,Kinzelbach W., Temporal behavior of a solute cloud in a heterogeneous porous medium. 3. Numerical simulations. Water Resour Res 2002;38:1118.
CrossRef (DOI)

[54] Kurbanmuradov OA, Sabelfeld KK., Stochastic Flow Simulation and Particle Transport in a 2D Layer of Random Porous Medium Transp Porous Med 2010;85:347–373.
CrossRef (DOI)

[55] Heße F, Prykhodko V, Schlter S, Attinger S., Generating random ﬁelds with a truncated power-law variogram: A comparison of several numerical methods. Environ Model Software 2014;55:32–48.
CrossRef (DOI)

[56] Bronstein I.N, Semendjajew KA, Musiol G, Muhlig H., Taschenbuch der Mathematik. Frankfurt am Main: Verlag Harri Deutsch; 2006.

Suciu_et_al

Towards a filtered density function approach for groundwater

N. Suciu $^{a,b,*}$ , L. Schüler $^{c,d}$ , S. Attinger $^{c,d}$ , 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}$ Institute of Geosciences, University of Jena Burgweg 11, 07749 Jena, Germany $^{d}$ Department Computational Hydrosystems, UFZ Centre for Environmental Research, Permoserstraße 15, 04318 Leipzig, Germany

Abstract

Evolution equations for probability density functions (PDFs) and filtered density functions (FDFs) of random concentrations weighted by conserved scalars are derived from Fokker-Planck equations describing stochastically equivalent processes in concentration-position spaces. This approach provides consistent numerical PDF/FDF solutions, given by the density in the concentration-position space of an ensemble of computational particles governed by the associated Itô equations. The solutions are obtained by a global random walk (GRW) algorithm, which is stable, free of numerical diffusion, and practically insensitive to the increase of the number of particles. The upscaled drift and diffusion coefficients describing the PDF transport in the physical space are estimated on single-trajectories of the advection-diffusion process governing the random concentration of the conserved scalar. The procedure uses the self-averaging property of the transport in velocity fields with short-range correlations. The drift and mixing coefficients describing the PDF transport in concentration spaces are parameterized by the trend and the noise inferred from analyses of simulated concentration time series as well as by classical mixing models. A Gaussian spatial filter, applied to a Kraichnan velocity field generator, is used to construct coarsegrained simulations (CGS) for FDF problems. The CGS approach helps to understand the significance of the FDF from a practical point of view and its relation to the PDF approach. The methodology to construct GRW-based PDF and FDF numerical solutions is illustrated for a problem of passive contaminant transport in heterogeneous groundwater systems.

Keywords: PDF/FDF methods, mixing, global random walk, groundwater
2010 MSC: 60J60, 60G60, 86A05, 65M75

1. Introduction

Geological formations are heterogeneous and their properties are normally not measurable everywhere. This lack of knowledge implies uncertainty in aquifer parameters like hydraulic conductivity. As a consequence, quantifying the transport of solutes through these formations is also uncertain. This uncertainty is dealt with by using a probabilistic description of the involved processes. Going beyond mean values and also considering the variance is a good starting point to take the uncertainty into account. But in the case of risk assessments, the need to predict extreme values of contaminant concentration becomes important. However, such occurrences cannot be described by the variance alone. Hence, the need to describe the solute concentrations by their complete probability density function (PDF) arises. Another major advantage of the PDF approach is the easy way to handle reactive solutes, for example solutes undergoing chemical reactions or decaying radioactively. Even highly nonlinear reactions can simply be included in the formulation without the need to approximate them.

PDF methods were originally developed in the context of modeling turbulent reacting flows as a powerful tool to close highly nonlinear terms arising from averaging chemical reaction rates [

1, 2, 3, 4, 5

]. The PDF approach is based on solving evolution equations for the one-dimensional (one-point one-time) Eulerian joint PDF of sets of variables describing the state of the system, for example velocity, chemical composition, turbulent

frequency, temperature or enthalpy. PDF equations are derived by different methods [1,2] from local balance equations governing the flow and the evolution of the thermochemical state of the system. The PDF equations are unclosed because they contain terms that cannot be determined by the one-dimensional PDF alone [2]. But, irrespective of the complexity of the set of state variables, the chemical source terms in the PDF equations are closed, namely, they have the same functional form as the unaveraged reaction rates in the local balance equations [6]. The unclosed terms which require modeling are those describing the turbulent frequency, the scalar mixing by molecular diffusion, and, when the velocity is not included among the state variables, the effects of the turbulent velocity fluctuations [4, 6].

Filtered density function (FDF) methods provide an alternative approach using spatially filtered quantities, that is, spatial averages, instead of stochastic averages on which the PDF methods are based [6, 7, 8, 9]. The meaning of the FDF is that of a PDF of state variables at scales smaller than the filter width [10,11]. FDF evolution equations can be derived by similar procedures and have the same structure, with source terms in closed form, as the PDF equations [7,5,6]. With rare exceptions, the only state variables considered in FDF methods are the scalars describing the thermochemical composition [6]. Unlike the PDF, which is a deterministic function, the FDF is a random quantity. Its expectation tends to the PDF only in the limit of a small filter width [6]. However, since, by invoking an ergodic theorem, the ensemble averaging may be seen as a filtering in space with a suitably chosen large filter width [10], the FDF should approach the PDF with increasing filter width. Even though this assertion has not been proved theoretically, in some cases the convergence may be demonstrated numerically, as we will see below in this article.

In modeling transport through highly heterogeneous natural groundwater systems, the randomness is introduced by the stochastic parameterization of the hydraulic conductivity, which accounts for the parameter uncertainty due to a lack of measurements. This stochasticity implies the randomness of the Darcy flow velocity and of the dependent variables of the transport equations [12]. While in earlier stochastic approaches the focus was mainly on the mean and in some cases the variance of the concentration, during the last decade the need to model the concentration PDF received an increased attention [12, 13, 14, 15, 16, 17, 18, 19]. Concentration PDFs of conserved scalars may be inferred without solving PDF evolution equations in case of small or moderate fluctuations of the hydraulic conductivity, modeled as a lognormal random function with finite correlation lengths. Then, a Gaussian shape of the concentration may be assumed or inferred, which is completely determined by its first two moments. The statistics of these moments, specified under various assumptions in a first-order perturbation approach, is finally "mapped", using the Gaussian functional shape of the concentration, onto the concentration PDF via numerical [12] or analytical techniques [14, 16]. Another favorable situation is that of stratified transport, when the Gaussian concentration can be expressed explicitly as a function of the hydraulic conductivity, with the only assumption of negligible transverse dispersion. This leads to an explicit functional dependence of the concentration PDF on the hydraulic conductivity [13]. This approach also provides the PDF of reacting chemical species if their concentrations are fully defined by monotonous functions of conserved scalars [13].

For advective-reactive transport, PDFs of reacting species can be computed by solving evolution equations similar to those used in turbulence [17]. Such PDF equations do not contain mixing terms, because molecular diffusion is neglected. The only closure problem concerns terms due to velocity fluctuations, which are modeled as effective, or upscaled, diffusion coefficients, leading to Fokker-Planck evolution equations [17, 19]. By considering the velocity among the state variables, velocity-concentration PDF equations similar to turbulence problems can be derived and no closure for velocity fluctuations is necessary. Mixing is modeled similarly to turbulence approaches and the concentration PDF is obtained by integrating the joint velocity-concentration PDF over the velocity state space [15]. Evolution equations of the concentration PDF weighted by a conserved scalar, which generalizes the mass density function used in turbulence [2], can be formulated as Fokker-Planck equations [19]. Closures are provided by stochastically upscaled diffusion coefficients [18, 19] and by mixing models, which are formulated as a diffusion in concentration space. The parameters for this diffusion process can be inferred from measured/simulated concentration time series [20].

While several PDF approaches are being developed [17, 15, 19], FDF methods were not yet used in modeling transport in groundwater systems. The purpose of this article is to assess the feasibility and the utility of the FDF approach in stochastic subsurface hydrology. To begin, we notice the similarities between large eddy simulations (LES) in turbulence [7,8,11] and some approaches to coarse-scale, or coarse-grained simulations (CGS) in porous media [21, 22, 23, 24, 25]. In both cases, spatial averages of dependent variables are used to coarsen the grid, whereas subgrid effects are modeled. The objective is to obtain results comparable to fine grid simulations at reduced computational costs. In case of reactive transport, the upscaled equations obtained by spatial filtering contain unclosed averaged reaction terms. In LES, the closure problem is solved by coupling the filtered equations to FDF evolution equations [11]. In subsurface hydrology, effective reaction rates needed to close the problem can only be determined for specific problems under simplifying assumptions [26]. FDF approaches can be used to
avoid the need to close the filtered reaction terms and CGS can be designed, based on numerical upscaling through volume averages [23,24] or, similarly to LES, by solving filtered equations [21,22].

Nevertheless, there are three major differences with respect to LES-FDF modeling for turbulent flows. The first one concerns the number of parameters. While only a few parameters are required to solve filtered LES equations [7, 8], upscaling flow and transport processes in groundwater, by either spatial or stochastic averaging, requires fields of parameters: the hydraulic conductivity [22,24] or the velocity field [26,18]. The second difference, related to the first one, is the origin of the randomness. Turbulent flows are governed by the deterministic Navier-Stokes equations but, for large Reynolds numbers, the flow velocity behaves as a random variable due the sensible dependence of the solution on initial and boundary conditions. This mathematical aspect corresponds to an experimental lack of reproducibility of the measurements in turbulent systems [2]. In groundwater systems, the spatial variability of the hydraulic conductivity cannot be completely described. Therefore, stochastic parameterizations by random space functions are used to account for this uncertainty. The flow equations are thus solved in a probabilistic sense [16]. In this case, randomness is caused by the uncertainty of the parameter fields propagated through the flow and transport equations, which have to be modeled as stochastic equations. The third difference is given by the available experimental data. In turbulence, detailed velocity, temperature, or concentration profiles are available from measurements. Consequently, in LES-FDF simulations the objective is to obtain a good agreement with the available measurements, at scales lying between that of the fine-scale simulations, which fully resolve the variability of the turbulent flows, and that of the ensemble averaged solutions of the Navier-Stokes and transport equations, for which the variability of the unresolved scales is modeled by the solution of the PDF equation [10,11]. In natural porous media the measurements are sparse and subject to uncertainty. Therefore, as long as the simulation relays on fields of stochastic parameters (hydraulic conductivity or velocity) and no detailed experimental data are available, the aim of the FDF approach should be a probabilistic quantification of the uncertainty for the whole hierarchy of scales. In other words, FDF approaches for modeling transport processes in groundwater systems, besides closing reaction terms, could alleviate the computational costs through estimations of the global PDF from filtered parameter fields and coarse-grained FDF simulations.

2. PDF/FDF equations

We consider a statistically homogeneous random velocity field

V (x)

with divergence-free samples. This field is a solution of the continuity and Darcy equations with a random parameterization of the hydraulic conductivity. Furthermore, we consider an isotropic diffusion process specified by a constant diffusion coefficient

D

, which describes the local dispersion in a saturated aquifer system. A system of reacting chemical species described by concentrations

C_{α} (x, t) \in Ω_{c}, x \in Ω_{x}, t \in R_{+}, α = 1, \dots, N_{α}

, is transported through the aquifer according to the system of balance equations

\begin{matrix} (1) & \frac{\partial C_{α}}{\partial t} + V_{i} \frac{\partial C_{α}}{\partial x_{i}} = D \frac{\partial^{2} C_{α}}{\partial x_{i} \partial x_{i}} + S_{α}, \end{matrix}

where

S_{α} (C)

denote the reaction rates. Since the velocity field

V

is a random function, the vector concentration

C (x, t)

is a random field as well.

The one-point one-time probability density function

f (c; x, t)

of the random concentration

C

solving (1) satisfies the PDF evolution equation

\begin{matrix} (2) & \frac{\partial f}{\partial t} + \frac{\partial}{\partial x_{i}} (V_{i} f) - \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} (D_{i j} f) = - \frac{\partial^{2}}{\partial c_{α} \partial c_{β}} (M_{α β} f) - \frac{\partial}{\partial c_{α}} (S_{α} f), \end{matrix}

where

V

and

D

are the stochastically upscaled drift vector and the diffusion tensor, respectively, and

M

is the conditional dissipation rate accounting for mixing by diffusion [19, 20]. As a common convention, a semicolon is used in writing the concentration PDF

f (c; x, t)

to emphasize the distinction between the value

c

taken by the random function

C (x, t)

in the state space

Ω_{c}

and the values

x

and

t

of its independent variables in

Ω_{x} \times R_{+}

[2,27].

The mostly used approaches to derive the PDF equation (2) are the delta function method and the test function method. The first one starts with the definition of the PDF given by the ensemble average of some delta function depending on random fields. For instance the concentration PDF is given by

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

, where the multidimensional

δ

function is defined by the product

δ (C (x, t) - c) = \prod_{α = 1}^{α = N_{α}} δ (C_{α} (x, t) - c_{α})

. Such singular space functions were shown to define consistent probability distributions [18] and formal calculus involving them corresponds to rigorous operations with Dirac functionals [28,19,20]. The PDF equation is obtained by evaluating

\partial f / \partial t

from formal derivatives of

δ

functions [1,13,5,20]. In the test function approach, the ensemble average of the operator

\partial / \partial t + V_{i} \partial / \partial x_{i}

applied to a test function of state variables

Q

, with suitable properties, is evaluated in two different ways: first, by interchanging derivatives and stochastic average and using the incompressibility of
the velocity field, and second, by multiplying the right hand side of (1) by

Q

and taking the ensemble average. The PDF equation follows by equating the two expressions for

⟨ \partial Q / \partial t + V_{i} \partial Q / \partial x_{i} ⟩

, performing integration by parts, and considering the vanishing of

Q

at the boundaries of

Ω_{c} [2, 4, 5, 19]

In LES approaches, spatial filters are used to separate the dynamics of the scales larger than the filter width from sub-filter effects. The former are obtained as solutions of filtered equations and the latter, corresponding to unresolved scales, are modeled

[10, 5, 6]

. The filtered value of a physical quantity

Q

is given by the spatial average

\begin{matrix} (3) & ⟨ Q ⟩_{λ} (x, t) = \int_{Ω_{x}} Q (x^{'}, t) G (x^{'} - x) d x^{'} \end{matrix}

with

λ

being the filter width, which is implicitly defined by the filter function

G (x^{'} - x)

. The filter

G (x^{'} - x

is spatially invariant, non-negative, and normalized to unity,

\int_{Ω_{a}} G (x^{'} - x) d x^{'} = 1

. Furthermore, the filtering operation commutes with differentiation [36]. Under these conditions, the FDF can be defined by

f_{λ} (c; x, t) = ⟨ δ (C (x, t) - c) ⟩_{λ}

and FDF equations can be derived, analogous to PDF equations, by the delta function method [5]. The FDF equation has the form of PDF equation (2), but its coefficients are now obtained by spatial averages (3) instead of ensemble averages. Note that, technically, neither the derivation of the PDF equation nor the derivation of the FDF equation requires the statistical homogeneity of the random velocity field [19,20]. Nevertheless, there is an important difference between PDF and FDF approaches. While in the case of the PDF approach to transport in groundwater the statistical homogeneity of the random velocity is essential for the existence of the stochastically upscaled coefficients [29], spatially filtered upscaled coefficients to be used in FDF equations can be derived by methods free of homogeneity assumptions [23,24]. This opens the perspective for FDF methods applicable to realistic situations, such as transport simulations conditional on measurements of hydraulic conductivity [30].

Because of the high dimensionality of the PDF/FDF equations (time and space dimensions and the

N_{α}

dimensions of the concentration space

Ω_{c}

) solutions by standard discretization methods (finite-differences or finiteelements) are impracticable [2,5]. Therefore, numerical solutions are usually obtained by Monte Carlo methods. "Eulerian particle methods" simulate the finite-difference solution of the PDF equation by locating an ensemble of

N

notional particles at each point of an Eulerian grid. These particles have assigned representative values of the state variable (e.g. concentration), initially distributed according to the initial PDF. The particles move on the grid following rules consistent with the finite-difference scheme and, as

N

tends to infinity, averaging over ensembles of particles converges to the expectation estimated from the finite-difference scheme [31]. Even if they are computationally simpler than other particle methods in use, Eulerian particle methods are numerically dissipative and have low spatial accuracy [32, 5]. Another Monte Carlo approach is the "stochastic Eulerian field method" which uses a representation of the PDF by an ensemble of stochastically equivalent space-time random fields, with the same one-point one-time PDF as the solution of the PDF equation [33, 11]. Stochastic fields are governed by partial differential equations with a linear multiplicative noise term, interpreted in either Itô [34] or Stratonovich [35] form. The computational effort of solving a partial differential equation for each field renders the method of stochastic Eulerian field less competitive than other Monte Carlo methods when large numbers of fields are required [5]. "Lagrangian particle methods" use systems of stochastic particles moving in continuous space, according to grid-free particle tracking procedures. Eventually, grids are used to compute averages and to interpolate the output of averaged or filtered transport equations to the particles positions. Lagrangian methods are now the dominant numerical approach for PDF equations. The performance of the above Monte Carlo approaches, with several variants, is analyzed in [5].

3. The Fokker-Planck approach

3.1. Direct Fokker-Planck approach

Lagrangian particle methods are based on the similarity between PDF/FDF equations of type (2) and FokkerPlanck equations. It is common to compare [7, 8, 15], or even to assimilate [36] PDF/FDF equations to FokkerPlanck equations and to use the associated Itô equations as a model for Lagrangian particles. Further, numerical solutions to equation (2) are constructed by imposing a uniform distribution of Lagrangian particles during the simulations. To understand this constraint, we consider the concentration PDF problem and note that the corresponding PDF fulfills the normalization condition

\int_{Ω_{c}} f (c; x, t) d c = 1

. On the other side, if (2) were a Fokker-Planck equation, its solution would be a PDF

p (c, x, t)

defined on the concentration-position state space

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

[19], which yields

\int_{Ω_{c}} p (c, x, t) d c = p_{x} (x, t)

by integration over

Ω_{c}

, that is, the position PDF of the Lagrangian particles. A uniform particle distribution,

p_{x} (x, t) \equiv

const, would suffice to make

p (c, x, t)

proportional to

f (c; x, t)

, which allows estimating the concentration PDF from the solution of the Fokker-Planck equation. As will be shown in
the following, a more general consistent relation between the two PDFs

f (c; x, t)

and

p (c, x, t)

may be established if a suitable weighting function

Θ

exists such that

\int_{Ω_{c}} Θ (c) f (c; x, t) d c = p_{x} (x, t)

In operator splitting schemes the transport of the PDF in physical space is treated in separate advection and diffusion steps [31]. These two steps solve equation (2) with the right-hand side set to zero [2]. This equation has the form of a Fokker-Planck equation describing the position PDF of a passive scalar. The corresponding Itô equation, used to simulate the transport step in Lagrangian particle methods [37], has the form

\begin{matrix} (4) & d X_{i} (t) = V_{i} (X, t) d t + d {\tilde{W}}_{i} (X, t) \end{matrix}

where

{X_{i}, i = 1, 2, 3}

are trajectories of an Itô diffusion process and

{\tilde{W}}_{i}

is a Wiener process with expectation

E {{\tilde{W}}_{i} (X, t)} = 0

and variance

E {{\tilde{W}}_{i}^{2} (X, t)} = 2 \int_{0}^{t} D_{i i} (X, t^{'}) d t^{'}

. The other fractional steps of the Lagrangian approach [2,37] correspond to the transport in concentration space and may be formulated in a general way as

\begin{matrix} (5) & d C_{α} (t) = M_{α} (C_{α} (t)) d t + S_{α} (C (t)) d t \end{matrix}

where

C_{α} (t) = C_{α} (X (t))

and the coefficients

M_{α}

are provided by mixing models for the term containing the dissipation rate

M_{α β}

in equation (2) [19].

To design a Monte Carlo method based on equations (4) and (5), we first have to establish a correspondence between the two mathematical objects involved, namely: the stochastic process

{C (t), X (t)}

, indexed by a single index,

t

, and the multi index random function

C (x, t)

. To do so, we remark that the Itô equation (4) and the corresponding Fokker-Planck equation [38] for the position PDF

p_{x}

\begin{matrix} (6) & \frac{\partial p_{x}}{\partial t} + \frac{\partial}{\partial x_{i}} (V_{i} p_{x}) = \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} (D_{i j} p_{x}) \end{matrix}

do not depend on concentrations (i.e. on the process

C_{α} (t)

or on the state space variable

c

) and may be used to describe any conserved scalar transported in the same system (with the same parameters

V_{i}

and

D_{i j}

), under the same initial conditions. Equation (6) also coincides with the equation satisfied by the ensemble averaged [15] or by the filtered scalar [36], which can be derived by multiplying the PDF/FDF equation of (2) without the reaction term by the scalar and by taking the ensemble average, or the filter average. It follows that for any conserved scalar

Θ (x, t)

solving equation (1) without reaction terms, the ensemble averaged

⟨ Θ ⟩

(or the filtered scalar

⟨ Θ ⟩_{λ}

) solves the Fokker-Planck equation (6). Thus, we have,

\begin{matrix} (7) & ⟨ Θ ⟩ (x, t) / Θ^{*} = p_{x} (x, t), where Θ^{*} = \int_{Ω_{x}} ⟨ Θ ⟩ (x, t) d x \end{matrix}

Equation (7) also holds for

⟨ Θ ⟩

replaced by

⟨ Θ ⟩_{λ}

.
Next, we consider a conserved scalar depending on

(x, t)

through concentrations of chemical species,

Θ (x, t) = Θ (C (x, t))

. Making use of the definition of the PDF via delta functions, we rewrite equation (7) as

\begin{matrix} (8) & ⟨ Θ (C (x, t)) ⟩ / Θ^{*} = \frac{1}{Θ^{*}} ⟨ \int_{Ω_{c}} Θ (c) δ (C (x, t) - c) d c ⟩ = \frac{1}{Θ^{*}} \int_{Ω_{c}} Θ (c) f (c; x, t) d c = \int_{Ω_{c}} p (c, x, t) d c \end{matrix}

The last equality in (8) is ensured by choosing

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

Relation (9) provides a correspondence between the one-point statistics of the random concentration

C (x, t)

and that of the stochastic process

{C (t), X (t)}

. Hence, the scalar

Θ

is the weighting function we were looking for in the beginning. Such a correspondence has been previously introduced in [19] and further analyzed in [20]. It is worth mentioning that the PDF

f (c; x, t)

of the random field cannot be uniquely determined in this way by the PDF

p (c, x, t)

of the stochastic process. Instead, based on (9), the system of Itô equations (4) and (5) can be used to compute the concentration PDF

f

weighted by a conserved scalar.

The conserved scalar

Θ

can be chosen as the sum of all species concentrations composing the reaction system,

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

. This sum is conserved in closed systems, as a consequence of mass conservation of the total amount of chemical elements contained in the reacting species molecules [20]. In particular, if the transport problem is formulated in terms of mass concentrations,

C_{α} = ρ_{α}

, and all the components of the fluid system are included among the

N_{α}

species, we have

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

, where

ρ

is the fluid density. Then, according to (7),

Θ^{*} = M

is the total mass of fluid in

Ω_{x}, ⟨ ρ ⟩ (x, t) = M p_{x} (x, t)

, and the correspondence relation (9) takes the form

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

which relates the mass density function

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

to the PDF

p (c, x, t)

of the system of stochastic particles used in Lagrangian Monte Carlo solution algorithms [2].

To facilitate the derivation of the FDF evolution equation, the filtered mass density function is usually written as

F_{λ} (c, x; t) = \int_{Ω_{x}} ρ (x^{'}, t) δ (C (x^{'}, t) - c) G (x^{'} - x) d x^{'} [8, 5]

. With

ρ (x, t) = ρ (C (x, t))

, using the FDF definition by the filtering operation (3) applied to

δ (C (x, t) - c)

and the "shifting property"

ρ (c) δ (C (x, t) - c) = ρ (C (x, t)) δ (C (x, t) - c)

[20], it is easy to recognize that

F_{λ} (c, x; t) = ρ (c) f_{λ} (c; x, t)

. Further, (8) still holds when the ensemble average is replaced by filtering and the FDF version of the correspondence relation is given by (10), with

f

replaced by

f_{λ}

If the weighting function

Θ

obeys the continuity equation

\begin{matrix} (11) & \frac{\partial Θ}{\partial t} + V_{i} \frac{\partial Θ}{\partial x_{i}} = 0 \end{matrix}

then, by using the correspondence (9), the Fokker-Planck equation governing the concentration-position PDF

p (c, x, t)

of the Itô process with trajectories given by (4) and (5) may be derived directly from the transport equations (1). To do so, either the test function [19] approach or the delta function [20] approach can be used. This transport equation, also satisfied by weighted concentration PDFs (in particular by

F

and

F_{λ}

), is

\begin{matrix} (12) & \frac{\partial p (c, x, t)}{\partial t} + ⟨ V_{i} ⟩ \frac{\partial}{\partial x_{i}} p (c, x, t) = - \frac{\partial}{\partial x_{i}} [⟨ U_{i} ∣ c ⟩ p (c, x, t)] - \frac{\partial}{\partial c_{α}} {[⟨ D \frac{\partial^{2} C_{α}}{\partial x_{i} \partial x_{i}} | c ⟩ + S_{α} (c)] p (c, x, t)} . \end{matrix}

Equation (12) contains the conditional averages given the concentration of the velocity fluctuation

U = V - ⟨ V ⟩

and of the diffusion flux

D \frac{\partial^{2} C_{α}}{\partial x_{i} \partial x_{i}}

. The first one is traditionally closed by a gradient-diffusion hypothesis

⟨ U_{i} ∣ c ⟩ f = - D_{i, j}^{*} \partial f / \partial x_{j} [2, 4]

. The second one can be related to the conditional dissipation rate, but only if the weighting factor

Θ

is a constant [19]. Hence, for

Θ =

const, the two closures introduce the coefficients [19]

\begin{matrix} (13) & V_{i} = ⟨ V_{i} ⟩ + \frac{\partial}{\partial x_{j}} D_{i, j}^{*}, D_{i j} = D + D_{i, j}^{*}, M_{α β} = ⟨ D \frac{\partial C_{α}}{\partial x_{i}} \frac{\partial C_{β}}{\partial x_{i}} | c ⟩ \end{matrix}

which, when inserted into (12), leads to the PDF equation (2) with

f (c; x, t) = p (c, x, t) Θ^{*} / Θ

, according to (9). In this case, (7) also implies a uniform position PDF equal to the inverse of the volume of the physical domain,

Θ / Θ^{*} = p_{x} (x, t) \equiv 1 / \int_{Ω_{x}} d x

. With

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

being the conditional PDF of the concentration given the position, it follows that

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

In particular, the framework of

PDF / FDF

approaches for constant density flows is obtained for

Θ = ρ =

const. Then, the coefficients of equation (2) are given by (13) and relation (14) estimates the concentration PDF by the solution of the system of Itô equations (4) and (5) [4, 15, 7, 36].

Within the direct Fokker-Planck approach, the coefficients of equation (12) are determined by those of the transport equations (1), through unclosed conditional averages. The joint concentration-position PDF solving this equation can be numerically approximated by the associated Itô equations (4) and (5) and interpreted as a weighted concentration PDF by using the correspondence relation (9). Relations between joint concentration-position PDFs of Itô processes and weighted concentration PDFs equivalent to (10) are the core of the Lagrangian particle Monte Carlo approaches. These approaches where introduced by somewhat involved arguments using the concept of "fluid particles" (rather questionable for systems undergoing diffusion [19]) [2, 27, 4].

3.2. Reverse Fokker-Planck approach

When the balance equation for the solvent is considered together with those for

N_{α - 1}

reacting species concentrations, then, summing up these equations, not only the sum of reaction rates but also that of diffusion fluxes vanishes. Then, it follows that the solution

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

of the continuity equation (11) is precisely the density

ρ

of the fluid system [39]. Considering the balance equations for all the components of the fluid is a natural choice in PDF approaches for reacting gas mixtures [2], where different components may have comparable weights. Including the carrying fluid among species components of a dilute solution transported in groundwater may cause difficulties in solving PDF/FDF problems. Even if the complicated balance equation for the solvent may be avoided by using the

N_{α - 1}

equations (1) and the continuity equation (11) [39], the numerical solution of the system of Itô equations (4) and (5) would be complicated by the need to consider the initial condition for the carrying fluid. If one neglects the variations of the solvent concentration, which is tantamount to considering constant density flows, in a direct Fokker-Planck approach, then equation (7) imposes a uniform position PDF. In turn, a space-time constant solution to equation (6) requires some relations between the spatial derivatives of the drift and
diffusion coefficients [19], which may not be fulfilled in case of statistically inhomogeneous velocity fields [29]. Such issues can be avoided in a Fokker planck approach which does not require the fulfillment of condition (11).

Such an alternative Fokker-Planck approach has already been suggested in [18]. Instead of specifying the system of Itô equations by the (modeled) coefficients of the Fokker-Planck equation (12), modeled Itô equations may be used to derive a Fokker-Planck equation. In our case, given the system of Itô equations (4) and (5), the corresponding Fokker-Planck equation is [38]

\begin{matrix} (15) & \frac{\partial p}{\partial t} + \frac{\partial}{\partial x_{i}} (V_{i} p) + \frac{\partial}{\partial c_{α}} (M_{α} p) = \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} (D_{i j} p) - \frac{\partial}{\partial c_{α}} (S_{α} p) . \end{matrix}

According to (9), the solution

p

of the Fokker-Planck equation (15) coincides with the concentration PDF

f

weighted by a conserved scalar

Θ

solving (1) without reaction terms but not necessarily solving the continuity equation (11). This "reverse" Fokker-Planck approach could be a valid option in modeling groundwater systems, for which dispersion coefficients

D_{i j}

are provided by various theoretical

[40, 16, 41]

or numerical

[23, 18]

methods and mixing models may be inferred from measurements [20]. Similar reverse approaches have been used in a hydrological context to obtain the Fokker-Planck equation for the concentration PDF of a process generated by a given mixing model for fixed positions in physical space [42], and, in a general setting, to formulate the evolution equation for the velocity-concentration PDF [15]. Also similar to this approach is the derivation of the evolution equation for the mass density function from a diffusion model for the velocity of discrete particles dispersed in turbulent flows presented in [43]. The novelty of our reverse Fokker-Planck approach is that it does not require that the conserved scalar

Θ

be a solution of the continuity equation.

3.3. Numerical solutions

A straightforward solution to the system of Itô equations (4) and (5) is given by a grid-free particle tracking in the concentration-position state space

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

. A particle follows a trajectory which starts at an initial position (

c_{0}, x_{0}

) drawn from the initial PDF

p (c, x, 0)

N

particles initialized in the same way are used to form a statistical ensemble. The approximation of

p (c, x, t)

t > 0

is given by the particle distribution in computational cells around (

c, x

). As an example, let us consider a two-dimensional passive transport in groundwater with initial concentration

C (x, 0)

set to

C = 1

in a rectangular initial support in the (

x_{1}, x_{2}

) plane. The initial PDF

p (c, x, 0)

is approximated by a uniform distribution of particles in the same rectangle lifted in a plane parallel to the (

x_{1}, x_{2}

) plane which intersects the concentration axis at

c = 1

. The actual joint concentration-position PDF

p (c, x, t)

is then approximated at the center of the cubic cells of a regular lattice in the (

c, x_{1}, x_{2}

) space by the number

n

of particles in each cell, through the number density

n / N

In Lagrangian particle approaches, the solution algorithm differs from the straightforward solution to the system of Itô equations. A "notional particle" carries a "composition"

c = {c_{1}, \dots, c_{N_{α}}}

of species concentrations in physical space. Equation (4) is solved for

N

notional particles and equation (5), solved for each particle, updates its composition [2,4]. For a fixed initial position

x_{0}

of the notional particle, the composition

c_{0}

consistent with the initial PDF

p (c, x, 0)

has to be extracted from the conditional PDF

p (c ∣ x, 0)

. It follows that the initial distribution of notional particles has to approximate the position PDF

p_{x} (x, 0)

, to ensure that the joint event "extracting

c_{0}

and

x_{0}

" belongs to the ensemble with the PDF

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

. The algorithm is mainly useful when one wants to impose a uniform position PDF

p_{x}

which allows representing the concentration PDF through conditional PDFs of the Itô process by using equation (14). As an illustration, for the case of two-dimensional passive transport discussed above, the notional particles are initially uniformly distributed and their concentrations are set to unity for particles inside the rectangular support of the initial concentration and to zero outside [15]. The concentration PDF

p (c; x, t)

may be approximated by histograms obtained from ensembles of particles in cells, further smoothed by spline functions [2], as taking spatially constant values in cells [15], or in a weak sense, through estimations of cell averaged quantities [27].

Both particle tracking and Lagrangian particle methods suffer from the increase of the computational costs with the number of particles and from numerical diffusion generated by interpolation of cell averages to the particles’ positions [44, 19]. Such inconveniences are overcome by using a global random walk algorithm (GRW) equivalent to a superposition of many weak Euler schemes for systems of Itô equations projected on regular lattices [45]. Unlike in sequential particle tracking procedures, in GRW algorithms all the

N

particles used to construct a numerical estimate of the probability density solving the Fokker-Planck equation are distributed on the lattice, according to the initial PDF, from the beginning of the simulation. The particles from each lattice site are then globally spread over new positions on the lattice determined by drift and diffusion coefficients, according to binomial distributions. This numerical procedure is by construction free of numerical diffusion. Moreover, the GRW algorithm is practically insensitive to the increase of the number of particles [18]. Since mean values are
defined at lattice sites, the GRW solutions are not affected by the artificial diffusion caused by cell averaging and interpolation in particle tracking and Lagrangian particle approaches. Details on implementation and convergence properties of the GRW algorithms can be found in [18, 45, 46]. For the PDF/FDF simulations presented in the following we will use the GRW algorithm presented in [19].

4. Coarse-grained FDF simulations

4.1. A $P D F / F D F$ problem for passive transport in aquifers

A two-dimensional problem for passive transport in saturated aquifers, previously considered for investigations on ergodicity, memory effects, and asymptotic behavior [47, 18], is chosen here to illustrate the principles of the CGS and the utility of the FDF approach in subsurface hydrology. The same problem has been used in our recent publications to design a GRW numerical solution for PDF equations [18, 19, 20]. Following the idea announced in the Introduction, we use spatial filtering to smooth the velocity field and to infer upscaled dispersion coefficients, needed to parameterize FDF equations for this two-dimensional transport problem, in order to test the convergence of the CGS-FDF solutions towards reference results obtained by PDF solutions and Monte Carlo simulations.

Samples of a log-normally distributed hydraulic conductivity

K

with

\ln K

of variance

σ_{\ln K}^{2} = 0.1

and Gaussian correlations, with correlation length set to

λ_{\ln K} = 1 m

and corresponding integral scale

I_{\ln K} = 0.09 m

, are numerically generated with the Kraichnan routine [48]. Further, samples of the random velocity field with ensemble mean components

⟨ V_{1} ⟩ = 1 m /

day and

⟨ V_{2} ⟩ = 0 m /

day are generated using the Kraichnan routine by the usual approximation for small variances of

\ln K

(see Appendix A). The passive transport is governed by equation (1) with the right hand side set to zero and

N_{α} = 1

. The local dispersion is described as an isotropic diffusion with constant coefficient

D = 0.01 m^{2} /

day. With these parameters, the longitudinal ensemble dispersion coefficient [47] takes the long time asymptotic value of

D_{11}^{ens} = D + ⟨ V_{1} ⟩ I_{\ln K} = 0.1 m^{2} /

day , so that

D_{11}^{ens} / D = 10

. The simulations were conducted in a two-dimensional computational grid larger than the largest extension of the plume (to mimic an infinite medium) with the GRW algorithm. The total number of particles released in every GRW simulation was

N = 10^{10}

. The initial condition was, for every velocity realization, a transverse slab of

1 m \times 100 m

. Statistical ensembles were obtained by repeating the simulations for

R = 256

R = 3072

realizations of the

K

field. As an observable quantity affected by uncertainty we consider the cross-section concentration at the plume center of mass obtained by summing the number of particles

n (x, y, t)

over transverse slabs

Δ x \times L_{y}

, where

Δ x = 1 m

and

L_{y}

is the the transverse dimension of the two-dimensional domain. This concentration,

C (x, t) = \frac{1}{Δ x L_{y}} \int_{0}^{L_{y}} \int_{x - Δ x / 2}^{x + Δ x / 2} n (x^{'}, y, t) d x^{'} d y

was estimated on the trajectory

x = ⟨ V_{1} ⟩ t

of the ensemble averaged flow velocity at one day intervals. Reference Monte Carlo estimates for concentrations and dispersion coefficients were obtained from the ensembles of transport simulations presented in [47], as well as from new simulations performed during the present study.

Within the frame of the reverse Fokker-Planck approach presented in Section (3.2), the PDF equation is derived from a particular form of the system of Itô equations (4) and (5),

\begin{matrix} (16) & d X (t) = V d t + \sqrt{2 D} d W_{1} (t) \\ (17) & d C (t) = V_{c} d t + \sqrt{2 D_{c}} d W_{2} (t) \end{matrix}

where

W_{1} (t)

and

W_{2} (t)

are two independent standard Wiener processes [19]. The Fokker-Planck equation (15) derived from (16) and (17) is

\begin{matrix} (18) & \frac{\partial p}{\partial t} + 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 transport of the cross-section averaged concentration

C

is essentially one-dimensional. For now, we assert that the upscaled diffusion coefficient

D

is the ensemble dispersion coefficient

D_{11}^{ens}

[18] in case of PDF simulations and, for FDF simulations, it is given by coarse-grained coefficients estimated in Section 4.3 below. The correctness of this assertion will be demonstrated numerically by the results presented in Section 5. Since mixing models inspired from turbulence studies have shown a poor performance in PDF simulations for groundwater [19], in our previous studies

[49, 19, 19]

we have chosen to model the term

M d t

from equation (5) as a diffusion in the concentration space described by the Itô equation (17). This model, as well as a classical mixing model and a linear combination of both models will be introduced in Section 4.2. Finally, the convergence with the filter width

λ

of the FDF solution based on the filtering operation (3) to the PDF solution will be investigated numerically in Section 5.

The PDF

f (c; x, t)

is estimated from the solution

p (c, x, t)

of the Fokker-Planck equation (18) via (9) with

Θ (c) = c

. The comparisons presented in figures 1 and 2 show that there is no significant difference between the concentration PDFs

f (c; x, t)

and the weighted PDFs

c f (c; x, t) / ⟨ C ⟩

, denoted by "pdf" and "w-pdf", respectively. This is a consequence of the chosen large transverse slab source, for which the transport is almost ergodic [47, 18], i.e.

C \approx ⟨ C ⟩

. Therefore, relation (9) may be approximated by (14). In Section 5, the concentration PDF

f (c; x, t)

will be estimated by the conditional PDF

p (c ∣ x, t)

, which is a solution of the Fokker-Planck equation (18). The conditional PDF

p (c ∣ x, t)

is obtained by a GRW simulation consisting of a superposition of about

10^{25}

particle tracking solutions of the Itô equations (16) and (17) on a regular lattice. The detailed implementation of the numerical scheme is presented in [19].

Figure 1: Concentration PDFs at

x = ⟨ V_{1} ⟩ t

, sampled at 2 days intervals, from

t = 0

t = 100

days, inferred from 3072 GRW simulations of the two dimensional transport problem (from right to left).

Figure 3: Time series of cross-section concentrations at the plume center of mass given by Monte Carlo simulations.

Figure 2: Cumulative distribution functions (cdf), corresponding to the PDFs presented in figure 1.

Figure 4: Increments

d C (t)

for the sample time series presented in figure 3.

4.2. Mixing models

The coefficients

V_{c}

and

D_{c}

occurring in equation (17) are estimated in the following by a stochastic time series analysis. An ensemble of 500 time series

C (t) = C (x = ⟨ V_{1} ⟩ t, t)

of simulated concentrations, previously described in [20], is presented in figure 3. As common in time series analysis,

t

is an integer and represents a dimensionless time. Since the velocity field is statistically homogeneous, these time series correspond to the concentration

C

sampled at the plume center of mass [47]. The corresponding increments

d C (t) = C (t + 1) - C (t)

, approximating the slope

d C (t) / d t

of a continuous time series, are shown in figure 4. The residual noise in the time series increments,

ξ (t) = d C (t) - ⟨ d C (t) ⟩

, shown in figure 5 has an approximately exponentially decaying amplitude.

Figure 5: Noisy components

ξ (t)

of the concentration increments

d C (t)

from figure 4.

Figure 7: GRW solutions for the PDF

f (c; x, t)

for different mixing models, compared to Monte Carlo estimates (MC), sampled at

x = V t

at 2 days intervals, from

t = 0

t = 100

days (from right to left).

Figure 6: Noisy components

ξ (t)

of the concentration increments normalized by their maximum amplitude

‖ ξ ‖

Figure 8: Cumulative distribution functions estimated from GRW-PDF simulations, compared to MC estimates, sampled at

x = V t

, for

t = 0, 10, 20, 30, 40

, and 100 days (from right to left).

After normalizing each

ξ

-sample by its maximum amplitude,

‖ ξ ‖ = max | ξ |

, the noise standardized in this way,

ξ / ‖ ξ ‖

, may be approximated by a white noise (figure 6). In turn, it follows that

ξ (t)

may be approximated by a "heteroskedastic process" [49] consisting of a white noise with rapidly decaying amplitude. The drift coefficient

V_{c}

in equation (17) was inferred from the slope of the ensemble mean concentration

⟨ C (t) ⟩

, represented by a thick line in figure 4. The maximum value of the diffusion coefficient

D_{c}

is specified by an initial amplitude of diffusive jumps in the GRW algorithm of

5 δ c

, where

δ c = 0.001

is the space step in concentration space. This value is chosen as close as possible to the standard deviation 0.0053 of the maximum amplitude

‖ ξ ‖

of the noise shown in figure 5. The time variation of the diffusion coefficient

D_{c}

which mimics the behavior of the noise in figure 5 is optimized through a crude calibration of the PDF simulations.

The time series (TS) mixing model, described above, accurately describes the transport of the PDF in concentration space but it fails to reproduce the narrowing of the PDF at large times and small values of the concentration shown in figure

1 [18, 19, 20]

. Such a narrowing corresponds to an "attraction" towards the mean of the concentration realizations which is well described by the Interaction by Exchange with the Mean (IEM) model [2, 7],

V_{c} = - C_{c} (C - ⟨ C ⟩) / (2 τ)

. For the purpose of the present illustration of the IEM model, we set the constant

C_{c}

to the mean value

C_{c} = 2

suggested by turbulence studies [2]. As for the the characteristic time

τ

we chose the time scale defined by the effective dispersion coefficient and a characteristic length scale

L

. Since the transport process is almost ergodic, the effective and the ensemble coefficients are almost identical [18], so that

τ = L / D

. In PDF simulations, the scale

L

is obviously given by the correlation length

λ_{\ln K}

of the log-hydraulic conductivity field. For FDF simulations the filter width should be accounted for, so, we propose the length scale

L = λ_{\ln K} + λ

. With
these parameters, the IEM model reads

\begin{matrix} (19) & V_{c}^{I E M} = - \frac{D}{λ_{\ln K} + λ} (C - ⟨ C ⟩) \end{matrix}

Since, as we will see, the IEM model is not satisfactory at all, we also consider a linear combination of TS and IEM models (hereafter referred to as the TS-IEM model),

\begin{matrix} (20) & V_{c}^{T S - I E M} (t) = \frac{T - t}{T} V_{c}^{T S} (t) + \frac{t}{T} V_{c}^{I E M} (t) \end{matrix}

where

T

denotes the total simulation time. For both IEM and TS-IEM models of the drift

V_{c}

, we also consider the diffusion term in equation (17), with the diffusion coefficient

D_{c}

provided by the TS model. In figures 7 and 8 we compare PDF-GRW simulations (

λ = 0

) carried out with the three mixing models. Similar FDF-GRW simulations are deferred to Section 5. As we see in figure 7, both the IEM model (19) and the TS-IEM model (20) yield a narrow PDF at very small concentrations but they behave differently at larger values of

c

. The IEM model produces smaller PDFs (figure 7) and fails to reproduce the transport in concentration space of the cumulative distribution functions (figure 8). The TS-IEM model behaves satisfactorily for small to intermediate times (concentration range from 0.4 to 0.7 in figure 7 and from 0.15 to 0.4 in figure 8). Since the TS model does not depend on the filter width

λ

, we shall use the TS-IEM model (20) to illustrate the effect of filtering on FDF simulations presented below in Section 5.

4.3. Filtered velocity fields and upscaled diffusion coefficients

The computation of the filtered velocity

⟨ V ⟩_{λ}

through (3) is described in detail in Appendix A. The final expression (A.4), used in our field generator code, reduces to the usual Kraichnan field generator when

λ = 0

. Therefore, the cases with

λ = 0

presented in the following correspond to the fine-grained, non-filtered velocity field. The effect of filtering upon the components of the Kraichnan velocity field is illustrated in figure 9, where one can see that with increasing

λ

the longitudinal and transverse components of the filtered velocity,

{⟨ V_{1} ⟩}_{λ}

and

{⟨ V_{2} ⟩}_{λ}

, approach their ensemble averages

⟨ V_{1} ⟩ = 1 m /

day and

⟨ V_{2} ⟩ = 0 m /

day, respectively.

Thanks to their self-averaging property, upscaled diffusion coefficients can be efficiently estimated from computations on single trajectories of the advection-diffusion process [50,18] specified in Section 4.1. Such estimations were validated numerically for diffusion in velocity fields with finite correlation lengths [50]. One expects these estimations also to be valid in the case of power law correlated velocities as long as their covariance sampled on trajectories is ergodic [18]. The method of self-averaging estimation is presented in Appendix B.

The longitudinal ensemble dispersion coefficient estimated by (B.3), as well as those given by Monte Carlo GRW simulations of two-dimensional transport [47], are defined as one-half of the mean slope of the time dependent variance of the ensemble process. The diffusion coefficient

D

from (16), needed in PDF/FDF simulations, is defined as one-half of the local variance, which, since the process has finite first and second spatial moments at finite times, is tantamount to one-half of the derivative of the global variance [18].

D

is approximated in the GRW-PDF scheme from the input parameter numerically estimated by (B.3) as

\begin{matrix} (21) & D (t) = \frac{(t + δ t) \tilde{D} (t + δ t) - t \tilde{D} (t)}{δ t} \end{matrix}

where

δ t

is the simulation time step.

4.4. Coarse-grained simulations of transport

Ensemble dispersion coefficients computed for diffusive transport in filtered velocity fields only account for the variability of the velocity at scales larger than the filter width

λ

. Therefore, as seen in figure 10, the coefficients

D_{11}^{ens} (t, λ)

are smaller than the coefficients obtained for fine-grained, unfiltered velocity fields (the case

λ = 0

) and they approach the local dispersion coefficient

D

λ

increases and the velocity fluctuations are smoothed out. The sub-filter variability is modeled by fictitious diffusion coefficients. When used in a coarse-grained description, such coefficients should retrieve the statistics of the spatial moments, in stochastic approaches [51, 40, 41], or they should produce results close to those from fine scale simulations, in numerical upscaling approaches [23,24].

Following this paradigm, we conducted GRW-CGS of the two dimensional transport problem formulated in Section 4.1 by using filtered velocities (

{⟨ V_{1} ⟩}_{λ}, {⟨ V_{2} ⟩}_{λ}

) and coarse-grained longitudinal diffusion coefficients obtained by adding the difference between fine-grained and filtered ensemble dispersion coefficients to the local dispersion coefficient

D

\begin{matrix} (22) & D_{11}^{c g} (t, λ) = D + δ D (t, λ), δ D (t, λ) = D_{11}^{e n s} (t) - D_{11}^{e n s} (t, λ) \end{matrix}

Figure 9: Samples of filtered velocity fields (

{⟨ V_{1} ⟩}_{λ}, {⟨ V_{2} ⟩}_{λ}

) computed by applying a Gaussian spatial filter of width

λ

to finegrained velocity fields (

V_{1}, V_{2}

) obtained with a Kraichnan random field generator.

Figure 11: Single realization CGS estimations of the ensemble dispersion coefficient

D_{11}^{ens}

for increasing filter width

λ

, performed with single realization coarse-grained diffusion coefficients

D_{11}^{c g} (t, λ)

Figure 10: Ensemble dispersion coefficients

D_{11}^{ens} (t, λ)

derived from GRW simulations for single realizations of the filtered velocity field compared to their ensemble average (smooth curves).

Figure 12: Single realization CGS estimations of the finegrained ensemble dispersion coefficient

D_{11}^{ens}

for increasing filter width

λ

, performed with ensemble averaged coarse-grained diffusion coefficients

D_{11}^{c g} (t, λ)

Because the corrections to the transverse coefficient were found to be negligible, we set

D_{22}^{c g} = D

. By using the representation of the dispersion coefficients in terms of trajectories of Itô processes, it is easy to see that the corrections

δ D

in (22) are entirely due to sub-filter effects only in absence of correlations between the filtered and the sub-filter velocity fields (Appendix C). Such correlations are generally non-vanishing in upscaling by spatial averages of the flow equations [21]. In analytical estimations by first order approximations in

σ_{\ln K}^{2}

the condition of vanishing correlations is ensured by integrating over disjoint domains in the Fourier space [40]. In our numerical approach the lack of correlations is ensured by using independent sets of random numbers to generate filtered and fine-grained velocity fields.

The longitudinal ensemble dispersion coefficients were computed in GRW simulations from the variance of the ensemble process

X_{1}^{ens} = X_{1} - ⟨ X_{1} ⟩

D_{11}^{ens} = ⟨ {(X_{1}^{ens})}^{2} ⟩ / (2 t)

, where

⟨ \cdot ⟩

is the average over the

N = 10^{10}

particles (realizations of the local diffusion process) and over the

R = 256

realizations of the velocity field. Single realization estimations of

D_{11}^{ens}

were computed by averaging

{(X_{1}^{ens})}^{2}

for

R = 1

over the

N

particles, while keeping the full averaged mean

⟨ X_{1} ⟩

. In figure 11 it is shown that CGS, performed with filtered velocity fields and coarse-grained diffusion coefficients (22) determined from single realization coefficients

D_{11}^{ens} (t, λ)

, recover the fine-grained ensemble coefficient computed in a single realization (corresponding to

λ = 0

) for large enough filter widths. Figure 12 shows that for the same single realization CGS, but with coarse-grained coefficients (22) computed from

D_{11}^{ens} (t, λ)

averaged over

R = 256

realizations, for a filter width

λ = 4 m

one recovers the ensemble averaged

D_{11}^{ens}

given by fine-grained simulations (

λ = 0

). This demonstrates the usefulness of CGS, which reduces

Figure 13: Ensemble dispersion coefficients

D_{11}^{ens} (t, λ)

for filtered velocity fields obtained by self-averaging estimations and corresponding standard deviations (represented by error bars). Solid black lines represent ensemble averaged coefficients

D_{11}^{ens} (t, λ)

obtained from 256 GRW simulations.

Figure 15: Cross section concentration

C (x_{1}, t)

recorded at

t =

10,50 , and 100 days and the concentration at the expected center of mass

C (x_{1} = ⟨ V_{1} ⟩ t, t)

(monotonous curves) given by single realizations GRW simulations with filtered velocities.

Figure 14: Estimations of the fine-grained ensemble averaged coefficient

D_{11}^{ens} (t)

for increasing filter width

λ

, obtained from 256 GRW simulations with filtered velocities and selfaveraging estimates of the coarse-grained diffusion coefficients

D_{11}^{c g} (t, λ)

Figure 16: Fine-grained ensemble averages corresponding to the concentrations shown in figure 15, recovered from 256 GRW simulations with filtered velocities and self-averaging coarse-grained diffusion coefficients

D_{11}^{c g} (t, λ)

the computational costs by a factor of 256 , provided that the coarse-grained coefficient

D_{11}^{c g} (t, λ)

can be obtained without carrying out ensembles of GRW simulations.

In the present numerical setup, the low variance

σ_{\ln K}^{2} = 0.1

allows us to generate velocity fields at given points in space (Appendix A) and to construct trajectories on which self-averaging estimations of the coefficients

D_{11}^{ens} (t, λ)

for filtered velocity fields are readily available (Appendix B). These self-averaging estimations, shown in figure 13, were used to infer the coarse-grained coefficients (22). With coefficients

D_{11}^{c g} (t, λ)

obtained in this way, we found that GRW simulations, for

R = 256

velocity realizations, yield an acceptable recovery of the fine-grained ensemble coefficient

D_{11}^{ens}

(figure 14). The suitability of the self-averaging estimations of the coarse-grained coefficients for ensemble average estimates could significantly alleviate the costs in Monte Carlo simulations by using discretization elements several times coarser than in fine-grained simulations.

Figure 15 shows the cross-section concentration

C (x_{1}, t)

recorded at three simulation times, together with the concentration sampled in time at the expected position of the center of mass

⟨ V_{1} ⟩ t

, computed from a single GRW simulation with filtered velocities and with an unchanged local diffusion coefficient

D

. An ensemble of

R = 256

GRW simulations, using self-averaging estimations of the coarse-grained coefficients

D_{11}^{c g} (t, λ)

and filtered velocities, recover the fine-grained ensemble averaged cross-section concentration with a good precision (figure 16). This result indicates that the self-averaging estimations of the coarse-grained coefficients

D_{11}^{c g} (t, λ)

are accurate enough for the purpose of simulating the transport of the PDF/FDF in physical space, which is governed by the
same equation (16) as the mean concentration.

5. PDF/FDF solutions

GRW-FDF solutions for the problem of passive scalar transport formulated in Section 4.1 were computed with drift coefficients for transport in physical space

V

(see equation (16)) given by filtered velocity fields

{⟨ V_{1} ⟩}_{λ}

, generated in a first-order approximation by the modified Kraichnan expression (A.4). The corresponding diffusion coefficients,

D

in equation (16), were computed by (21) from coarse-grained longitudinal ensemble dispersion coefficients

D_{11}^{c g} (t, λ)

. The coefficients

D_{11}^{c g} (t, λ)

were obtained by (22) from the mean values of the self-averaging ensemble dispersion coefficients presented in figure 13. The GRW-PDF solutions were obtained similarly, with

V

given by the ensemble average

⟨ V_{1} ⟩ = 1 m /

day of the fine-grained velocity field and

D

given by the mean value of the self averaging ensemble dispersion coefficient

D_{11}^{ens} (t, 0)

(case

λ = 0

in figure 13). The coefficients describing the transport in concentration space,

V_{c}

and

D_{c}

in equation (17), were those of the TS mixing model in PDF simulations and those of the combined TS-IEM model (20) in FDF simulations.

In the present reverse Fokker Planck approach to PDF/FDF solutions, the weighting factor is the scalar concentration itself,

Θ (c) = c

. Hence, as follows from (7), the mean concentration equals the position PDF,

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

, which is simply obtained in the GRW-PDF/FDF simulation by summing up the joint concentrationposition PDF/FDF

p (c, x, t)

over the c-axis. We stress here the advantage of using the reverse Fokker-Planck approach. In a direct approach, the scalar concentration

c

cannot serve as a weighting factor

Θ

, because it undergoes diffusion and does not verify the continuity equation (11).

Figure 17:

⟨ C ⟩_{λ} (x, t)

recorded at fixed times

t = 10

, 50, and 100 days and

⟨ C ⟩_{λ} (t)

recorded at

x = ⟨ V_{1} ⟩

t. TS-IEM results for increasing

λ

are compared to Monte Carlo estimates (MC), and PDF solutions using the TS model (PDF).

Figure 18: Cumulative distribution functions (cdf) for the same cases as in figure 17, estimated at

x = ⟨ V_{1} ⟩ t

, for

t = 10, 20

, and 30 days (from right to left).

Results for the mean concentration

⟨ C ⟩_{λ}

obtained from GRW-FDF simulations using the linear combination of mixing models TS-IEM (20) for increasing filter width

λ

are compared to Monte Carlo results and to GRWPDF simulations using the TS mixing model in figure 17. We can see that the accuracy of the GRW-FDF solutions increases with

λ

and for

λ \geq 5 m

the results obtained by the three methods practically coincide. We also remark that the FDF solution for small

λ

is less smooth than similar results obtained from single realization GRW simulations of the two dimensional transport shown in figure 15. In both cases, one can see a misfit between

⟨ C ⟩ (t)

, respectively

⟨ C ⟩_{λ} (t)

, both recorded at the expected center of mass

x = ⟨ V_{1} ⟩ t

, and the peaks in the spatial distribution of the concentration located at the position of the actual center of mass, which is random and generally does not coincide with its expectation.

Cumulative distribution functions sampled at several times on the same trajectory

x = ⟨ V_{1} ⟩ t

are shown in figure 18. The GRW-FDF results with the TS-IEM model are less accurate at extreme values of the filter width. The differences with respect to the Monte Carlo and GRW-PDF results are mainly significant for large filter widths (

λ = 10 m

), when the decrease of the term (19) reduces the weight of the IEM drift in the TS-IEM linear combination (20). The inaccurate results for the smallest filter width

λ = 1 m

have to be related to the failure of the filtered velocity

{⟨ V_{1} ⟩}_{λ}

and coarse-grained dispersion coefficient

D_{11}^{c g} (t, λ)

to simulate the behavior of the mean concentration (see figure 17). The TS-IEM model provides acceptable results for intermediate values of the filter width, from

λ = 3 m

λ = 5 m

. At early times, the performance of the TS-IEM model could be improved

Figure 19: GRW-FDF solutions

f (c; x, t)

for the TS-IEM mixing model, compared to Monte Carlo estimates and PDF solutions for the TS model, sampled at

x = ⟨ V_{1} ⟩ t

at 2 days intervals, from

t = 0

t = 100

days (from right to left).

Figure 20: Cumulative distribution functions (cdf) corresponding to PDFs/FDFs from figure 19, estimated at

x = ⟨ V_{1} ⟩ t

, for

t = 0, 10, 20, 30

, and 50 days (from right to left).

by increasing the left-drift, towards small concentrations, of the cumulative distributions from figure 18. Since according to (20) the IEM weight is small at early times, the increase of

V_{c}^{T S - I E M}

can be achieved by augmenting

V_{c}^{I E M}

. Also, the larger the filter width

λ

, the stronger should the increase of

V_{c}^{I E M}

be.

For the favorable cases

λ = 3 m

and

λ = 5 m

, the FDFs presented in figure 19 show some agreement with the Monte Carlo results only at small times. The corresponding cumulative distribution functions (figure 20) are also close to Monte Carlo and PDF curves at moderately large times smaller than

t = 50

days, which is half of the simulation time

T = 100

days. At larger times, when both

V_{c}^{I E M}

and

V_{c}^{T S}

contribute to

V_{c}^{T S - I E M}

, the TS-IEM model produces a too large left-drift of the FDF. Comparing to the case

λ = 0

shown in figure 7, we see that the excessive drift increases with

λ

. Summarizing,

V_{c}^{I E M}

should be larger at small times and smaller at large times. That is, its magnitude should depend on time similarly to that of

V_{c}^{T S}

, shown in figure 4. Hence, in a TS-IEM model,

V_{c}^{I E M}

and

V_{c}^{T S}

should rather be proportional than linearly combined.

A combination of TS and IEM models is desirable because it improves the overall bad performance of the IEM model (see figure 8) and could produce the increase of the maximum PDF at large times shown by the Monte Carlo simulations (figures 9 and 19). The tests with the TS-IEM model used in the FDF simulations presented here (figures 18-20) suggest that the optimal combination between the two models should be nonlinear. Moreover, the results presented in figures 18 and 19 also indicate that the dependence on

λ

could be more complex than in classical IEM models.

6. Conclusions

PDF/FDF evolution equations are formulated as Fokker-Planck equations describing processes in concentrationposition spaces stochastically equivalent to concentration random fields. Within this approach, which we call a "reverse Fokker-Planck approach", the PDF (or the FDF) weighted by a conserved scalar is approximated numerically by the one-time PDF of the Itô process in the concentration-position space described by the Fokker-Planck equation. The conserved scalar which weights the one-point one-time PDF/FDF to render it equivalent to the one-time PDF of the concentration-position process is not involved in chemical reactions but, unlike in classical approaches, it does not necessarily solve a continuity equation. For systems of mobile reacting species in groundwater with species-independent diffusion coefficients, the conserved scalar can be chosen as the sum of molecular species concentrations, governed by an advection-diffusion equation.

The numerical solution is provided by the density on a lattice imbedded in the concentration-position space of an ensemble of computational particles governed by the GRW algorithm, which is stable, free of numerical diffusion, and practically insensitive to the increase of the number of particles. The transport of the system of particles in the physical space is described by drift and diffusion coefficients obtained through stochastic upscaling, in case of PDF simulations, or by spatial coarse-graining, for FDF simulations. For small variances of the loghydraulic conductivity, the upscaled coefficients can be determined through self-averaging estimations using a single trajectory of diffusion in the fine-grained velocity field, for PDF, or in the spatially filtered velocity field, for FDF simulations. Some tests done in the present study show that the FDF converges to the PDF as the filter
width increases up to a few correlation lengths of the log-hydraulic conductivity. This demonstrates the utility of the coarse-grained FDF simulation as a mean to estimate the global PDF, describing the whole variability of the velocity field, at significantly reduced computational costs.

One of the most stringent problems in PDF/FDF approaches is to design the appropriate mixing model which describes the transport in concentration spaces of the system of computational particles. For the problem considered here, we found that an advection-diffusion process in concentration space, inferred from time series analyses of simulated concentrations, yields a correct solution for the transport of the cumulative distribution towards low concentration values. However, this time series model cannot simulate the narrowing and the raising peak of the maximum PDF/FDF values at large times and small concentrations. A combination between the TS model and the classical IEM model could recover this large time behavior. The linear combination TS-IEM used in our CGS-FDF approach yields an acceptable accuracy only for moderately large times and filter widths. A suggestion would be to test nonlinear combinations of the two mixing models. However, the TS model is limited to cases where the path in physical space, on which the concentration of the conserved scalar has to be recorded, is given. For the problem considered here, the cross-section concentration was sampled on the mean-flow trajectory. A promising solution for more complex situations seems to be that of an IEM model with time variable coefficients, achieving a time-decaying magnitude of the drift in concentration space similar to that of the TS model.

Volume averaging or multi-grid numerical upscaling procedures can be used to extend the CGS-FDF strategy illustrated here to problems of reactive transport in conditions of highly variable and statistically inhomogeneous hydraulic conductivity of the groundwater system. Nonetheless, the reverse Fokker-Planck approach could also be used to design robust numerical schemes based on the GRW algorithm to solve PDF/FDF problems for other fluid systems, such as those studied in turbulence and combustion theory.

Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft-Germany under Grant SU 415/2-1, KN 229/15-1, SU 415/2-1, by the Jülich Supercomputing Centre-Germany as part of the Project JICG41, and by the Bundesministerium für Wirtschaft und Energie-Germany as part of the H-DuR Project (funding number: 02 E 11062 B).

Appendix A. Filtered Kraichnan fields

A Gaussian random velocity field, corresponding to a log-normal hydraulic conductivity, can be generated in a first order approximation by the Kraichnan method [48] as a randomized spectral representation given by the well known formula

\begin{matrix} (A.1) & V_{i} (x) = ⟨ V_{1} ⟩ δ_{i 1} + σ_{\ln K} ⟨ V_{1} ⟩ \sqrt{\frac{2}{N}} \sum_{j = 1}^{N} p_{i} (k^{(j)}) \cos (k^{(j)} \cdot x + ϕ^{(j)}), \end{matrix}

where

i = 1, \dots, d

and

d

is the spatial dimension of the problem [52,53]. The simple representation (A.1) is obtained by choosing the probability density of the random vector

k^{(j)} = (k_{1}^{(j)}, \dots, k_{d}^{(j)})

as given by the Fourier transform of the correlation function of the statistically homogeneous Gaussian random field

\ln K

divided by the variance

σ_{\ln K}^{2}

(see e.g. [54]). In the present study, we consider a Gaussian correlation of the

\ln K

field which, in the general anisotropic case, is given by

\begin{matrix} (A.2) & ⟨ \ln K (x) \ln K (x + r) ⟩ = σ_{\ln K}^{2} \exp (- \sum_{i = 1}^{d} \frac{r_{i}^{2}}{λ_{\ln K, i}^{2}}) . \end{matrix}

The random vectors

k^{(j)}

are thus distributed according to the PDF

f (k^{(j)}) = \prod_{i = 1}^{d} π^{1 / 2} λ_{\ln K, i} \exp [- π^{2} {(k_{i}^{(j)} λ_{\ln K, i})}^{2}] = \prod_{i = 1}^{d} f (k_{i}^{(j)}),

that is, they have mutually independent components. With the change of variables

k_{i}^{(j)} = η / (\sqrt{2} π λ_{\ln K, i})

, where

η

is a normally distributed random variable of zero mean and unit variance, the corresponding probability measures of the components become Gaussian,

f (k^{(j)}) d k_{i}^{(j)} = (2 π)^{- 1 / 2} \exp (η^{2} / 2) d η

. The random components

k_{i}^{(j)}

are therefore extracted from Gaussian distributions with vanishing mean and variances

1 / (2 π^{2} λ_{\ln K, i}^{2})

The phases

ϕ^{(j)}

, also being independent random numbers, are extracted from a uniform distribution in the interval

[0, 2 π]

. Randomized spectral representations like (A.1), are particular cases of more general randomization formulas using sums of cosine and sine functions [54, 55]. In our case, the general representation is obtained when

\sqrt{2} \cos (k^{(j)} \cdot x + ϕ^{(j)})

in (A.1) is replaced by

ξ^{(j)} \cos (2 π k^{(j)} \cdot x) + ζ^{(j)} \sin (2 π k^{(j)} \cdot x)

, where

ξ^{(j)}

and

ζ^{(j)}

are normal variables, mutually independent and independent of

k_{i}^{(j)}

. Our tests have shown that the two representations ar equally accurate, with the only difference that the general formula is about two times slower.

The projector

p (k) = e_{1} - \frac{k_{1} k}{k^{2}},

with

e_{1}

being the unit vector to which the mean flow is aligned, ensures the incompressibility of the flow. We note that when dropping the first term, as well as

⟨ V ⟩

and

p_{i} (k^{(j)})

from the second term, (A.1) reduces to the randomization formula for the fluctuations of the

\ln K

field.

We choose a Gaussian filter function

\begin{matrix} (A.3) & G (x^{'} - x) = {(\frac{b}{π})}^{d / 2} \exp [- b {(x^{'} - x)}^{2}] \end{matrix}

where

b = a / λ^{2}, a \in R^{+}

, and

λ

is the filter width. Using (A.3) in the filtering operation (3) applied to equation (A.1), one obtains

\begin{aligned} {⟨ V_{i} ⟩}_{λ} (x) & = \int_{- \infty}^{\infty} d x^{'} {(\frac{b}{π})}^{d / 2} [⟨ V_{1} ⟩ δ_{i 1} + σ_{\ln K} ⟨ V_{1} ⟩ \sqrt{\frac{2}{N}} \sum_{j = 1}^{N} p_{i} (k^{(j)}) \cos (k^{(j)} \cdot x + ϕ^{(j)})] \exp [- b {(x^{'} - x)}^{2}] \\ = ⟨ V_{1} ⟩ δ_{i 1} + σ_{\ln K} ⟨ V_{1} ⟩ \sqrt{2 / N} \sum_{j = 1}^{N} p_{i} (k^{(j)}) \underset{= I}{\underset{⏟}{\int_{- \infty}^{\infty} d x^{'} \cos (k^{(j)} \cdot x^{'} + ϕ^{(j)}) {(\frac{b}{π})}^{d / 2} \exp [- b {(x^{'} - x)}^{2}]}}, \end{aligned}

where

\int_{- \infty}^{\infty} d x^{'}

stands for

\int_{- \infty}^{\infty} \dots \int_{- \infty}^{\infty} d x_{1}^{'} \dots d x_{d}^{'}

.
The integral

I

can be simplified by the substitution

q = x^{'} - x

, which results in

I = {(\frac{b}{π})}^{d / 2} \int_{- \infty}^{\infty} d q \cos (k^{(j)} \cdot q + k^{(j)} \cdot x + ϕ^{(j)}) \exp (- b q^{2}) .

The cosine function has an argument with

2 d + 1

terms, of which

d

depend on a component of the vector q. The goal now is to separate these

d

terms from the rest. This separation can be achieved by a trigonometric identity. If we take a look at this identity for fewer terms

\begin{aligned} \cos (u + v) = & \cos (u) \cos (v) - \sin (u) \sin (v) \\ \cos (u_{1} + u_{2} + v) = & \cos (u_{1}) \cos (u_{2}) \cos (v) - \sin (u_{1}) \sin (u_{2}) \cos (v) \\ - \sin (u_{1}) \cos (u_{2}) \sin (v) - \cos (u_{1}) \sin (u_{2}) \sin (v) \\ \dots, \end{aligned}

we see that there is always only one term which has no sine functions. As sine functions are antisymmetric and the exponential function with which they are multiplied under the integral is symmetric, the resulting term is antisymmetric and the evaluated integral over all these terms vanishes. Hence, only the single term with no sine functions gives a non-zero value. Thus, by using the result

\int_{- \infty}^{\infty} d q \cos (u q) \exp (- w^{2} q^{2}) = \sqrt{π} \exp (- u^{2} / 4 w^{2}) / w (

see [56]), the integral gives

\begin{aligned} I & = {(\frac{b}{π})}^{d / 2} \int_{- \infty}^{\infty} d q [\prod_{α = 1}^{d} \cos (k_{α}^{(j)} q_{α}) \cos (k^{(j)} \cdot x + ϕ^{(j)}) + \dots] \exp (- b q^{2}) \\ = \cos (k^{(j)} \cdot x + ϕ^{(j)}) \exp [- \frac{{(k^{(j)})}^{2}}{4 b}] \end{aligned}

The final expression for the Gaussian filtered velocity field is

\begin{matrix} (A.4) & {⟨ V_{i} ⟩}_{λ} (x) = ⟨ V_{1} ⟩ δ_{i 1} + σ_{\ln K} ⟨ V_{1} ⟩ \sqrt{2 / N} \sum_{j = 1}^{N} p_{i} (k^{(j)}) \cos (k^{(j)} \cdot x + ϕ^{(j)}) \exp [- \frac{{(k^{(j)} λ)}^{2}}{4 a}] . \end{matrix}

For the present computations of filtered velocities,

a \in R^{+}

was chosen to be

a = 2

[25]. The only difference compared to the unfiltered velocity equation (A.1) is one factor for each mode, which can be precalculated and saved before the actual simulations start.

The number of Fourier modes was set to

N = 6400

for all the simulations performed in the present study. We considered the isotropic version of the correlation function (A.2), specified by the correlation length

λ_{\ln K} = 1

m and by the variance

σ_{\ln K}^{2} = 0.1

, which ensures the accuracy of the first-order approximation. With these parameters, integrating (A.2) from zero to infinity along any straight line through the origin of the coordinate system one obtains the integral scale

I_{\ln K} = \sqrt{π} σ_{\ln K}^{2} λ_{\ln K} / 2 \approx 0.09 m

Appendix B. Self-averaging estimations of diffusion coefficients

We consider a component of the ensemble process

{X_{t}, t \geq 0}

of mean zero starting from

X_{0} = 0

, obtained by subtracting the ensemble mean

⟨ X_{i, t} ⟩

from one of the components

X_{i, t}, i = 1, 2

, of the transport process. If the continuous time interval

[0, t]

is partitioned into

S

subintervals of equal length

τ

, so that

t = S τ

, a position on the trajectory

X_{t}

and its square

X_{t}^{2}

can be expressed in terms of increments

δ X_{s} = X_{s τ} - X_{(s - 1) τ}

X_{t} = \sum_{s = 1}^{S} δ X_{s}, X_{t}^{2} = \sum_{s = 1}^{S} {(δ X_{s})}^{2} + 2 \sum_{r = 1}^{S - 1} \sum_{s = 1}^{S - r} δ X_{s} δ X_{s + r} .

Further, introducing the time averages

\begin{matrix} (B.1) & ρ (r) = \frac{1}{S - r} \sum_{s = 1}^{S - r} δ X_{s} δ X_{s + r} \end{matrix}

X_{t}^{2}

can be rewritten as

\begin{matrix} (B.2) & X_{t}^{2} = S ρ (0) + 2 \sum_{r = 1}^{S - 1} (S - r) ρ (r) . \end{matrix}

Since the averages

ρ

defined by (B.1) are similar to stochastic correlations, (B.2) is a single-realization form of the Taylor formula for processes with stationary variance [18].

The process of the increments

δ X_{s}

is determined through advective displacements by the velocity field sampled on trajectories [50], which, in a first order approximation, has ergodic covariances [18]. Therefore, one expects that the ergodic estimation (B.1) converges to the correlation of

δ X_{s}

, so that (B.2) yields a good approximation of the variance of

X_{t}

. Then, the diffusion coefficient can be estimated on a single sample of

X_{t}

\tilde{D} = X_{t}^{2} / (2 S τ)

. With this, (B.2) yields

\begin{matrix} (B.3) & \tilde{D} = \frac{1}{2 τ} ρ (0) + \frac{1}{τ} \sum_{r = 1}^{S - 1} (1 - \frac{r}{S}) ρ (r) . \end{matrix}

Since, as

r

gets closer to

S

the time averages

ρ

are poor ergodic estimates of the correlations, a rigorous ergodicity statement is not possible in this case. Nevertheless, the reliability of the self-averaging estimator (B.3) has been demonstrated by numerical experiments [50]. In practice, acceptable results are obtained if the estimation time

S τ

is larger than the total simulation time

T

. For the PDF simulations presented here, with

T = 100

days, we have chosen

S τ = 400

days. For an increased accuracy, the ergodic estimations (B.1) of the correlations

ρ

were further averaged over

10^{6}

paths

S

sampled on the same trajectory. With these parameters, the computation time for the estimation (B.3) of the diffusion coefficient is of about 4 minutes. For comparison, on the same computing platform, the computation of an ensemble of GRW simulations which provides the time variation of the ensemble dispersion coefficients over 100 days requires about 10 minutes and 256 processors.

Appendix C. Coarse-grained diffusion coefficients

Decomposing the velocity field into ensemble mean and fluctuation,

V = ⟨ V ⟩ + U

, and taking its statistical homogeneity into account, we see that it is related to the filtered velocity by

\begin{matrix} (C.1) & V = ⟨ V ⟩ + ⟨ U ⟩_{λ} + U_{λ} \end{matrix}

where

⟨ U ⟩_{λ}

is the filtered fluctuation of the fine-grained velocity and

U_{λ}

is the residual sub-filter velocity fluctuation about the filtered velocity,

⟨ V ⟩_{λ} = ⟨ V ⟩ + ⟨ U ⟩_{λ}

. The longitudinal component of the transport process starting from

X = 0

is described by the Itô equation

\begin{matrix} (С.2) & X_{1} (t) = \int_{0}^{t} V_{1} (t^{'}) d t^{'} + W (t) \end{matrix}

where

V_{1} (t^{'}) = V_{1} (X (t^{'}))

and

W (t)

is a Wiener process of mean zero and variance equal to

2 D t

. For statistically homogeneous velocity fields, the ensemble average of (C.2) is

⟨ X_{1} ⟩ (t) = ⟨ V_{1} ⟩ t

. The longitudinal ensemble process

X_{1}^{ens} (t) = X_{1} (t) - ⟨ X_{1} (t) ⟩

verifies the same equation (C.2), with the velocity component

V_{1}

replaced by its fluctuation

U_{1}

, and the half derivative of its variance gives the ensemble dispersion coefficient in the following form [18]:

\begin{matrix} (С.3) & D_{11}^{e n s} (t) = D + \int_{0}^{t} ⟨ U_{1} (t) U_{1} (t^{'}) ⟩ d t^{'} \end{matrix}

From (C.1) and (C.2) one obtains the equivalent expression

\begin{aligned} D_{11}^{e n s} (t) = D & + \int_{0}^{t} ⟨ {⟨ U_{1} ⟩}_{λ} (t) {⟨ U_{1} (t^{'}) ⟩}_{λ} ⟩ d t^{'} \\ + \int_{0}^{t} ⟨ U_{λ, 1} (t) U_{λ, 1} (t^{'}) ⟩ d t^{'} \\ (С.4) & + \int_{0}^{t} ⟨ {⟨ U_{1} ⟩}_{λ} (t) U_{λ, 1} (t^{'}) ⟩ d t^{'} + \int_{0}^{t} ⟨ {⟨ U_{1} ⟩}_{λ} (t^{'}) U_{λ, 1} (t) ⟩ d t^{'} \end{aligned}

The first line in (C.4) corresponds to (C.3) with

U_{1}

replaced by

{⟨ U_{1} ⟩}_{λ}

, that is to the coefficient

D_{11}^{ens} (t, λ)

derived for a filtered velocity field. The second line is given by the integral from (C.3) written for the correlation of the sub-filter velocity fluctuations. As for the last line in (C.4), it describes a contribution to the ensemble coefficient produced by correlations between filtered and sub-filter velocity fluctuations. If this last contribution vanishes, then the contribution from the second line of (C.4) corresponds to the correction

δ D_{11}^{ens} (t, λ)

of the coarse-grained diffusion coefficient defined by (22).
[1] Pope SB. The probability approach to the modelling of turbulent reacting flows. Combust Flame 1976;27:299-312. http://dx.doi.org/10.1016/0010-2180(76)90035-3
[2] Pope SB. PDF methods for turbulent reactive flows. Prog Energy Combust Sci 1985;11(2):119-192. http://dx.doi.org/10.1016/0360-1285(85)90002-4.
[3] Pope SB. Turbulent Flows. Cambridge: Cambridge University Press: 2000.
[4] Fox RO. Computational Models for Turbulent Reacting Flows. New York: Cambridge University Press; 2003.
[5] Haworth DC. Progress in probability density function methods for turbulent reacting flows. Prog Energy Combust Sci 2010;36:168-259. http://dx.doi.org/10.1016/j.pecs.2009.09.003.
[6] Haworth DC, Pope SB. Transported probability density function methods for Reynolds-averaged and large-eddy simulations. In: Echekki T, Mastorakos E, editors. Turbulent combustion modeling. Fluid mechanics and its applications, vol. 95. Dordrecht: Springer; 2011. p. 119-42. http://dx.doi.org/10.1007/978-94-007-0412-1.6.
[7] Colucci PJ, Jaberi FA, Givi P. Filtered density function for large eddy simulation of turbulent reacting flows. Phys Fluids 1998;10(2):499515. http://dx.doi.org/10.1063/1.869537.
[8] Jaberi FA, Colucci PJ, James S, Givi P, Pope SB. Filtered mass density function for large-eddy simulation of turbulent reacting flows. J Fluid Mech 1999;401:85-121. http://dx.doi.org/10.1017/S0022112099006643.
[9] McDermott R, Pope SB. A particle formulation for treating differential diffusion in filtered density models. J Comput Phys 2007;226:947993. http://dx.doi.org/10.1016/j.jcp.2007.05.006.
[10] Heinz S. Unified turbulence models for LES and RANS, FDF and PDF simulations. Theor Comput Fluid Dyn 2007;21:99118. http://dx.doi.org/10.1007/s00162-006-0036-8.
[11] Dodoulas IA, Navarro-Martinez s. Large eddy simulation of premixed turbulent flames unsing probability density approach. Flow Turbulence Combust 2013;90:645-678. http://dx.doi.org/10.1007/s10494-013-9446-z.
[12] Schwede RL, Cirpka OA, Nowak W, Neuweiler I. Impact of sampling volume on the probability density function of steady state concentration. Water Resour Res 2008;44(12):W12433. http://dx.doi.org/10.1029/2007WR006668.
[13] Sanchez-Vila X, Guadagnini A, Fernndez-Garcia D. Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers. Math Geosci 2009;41:32351. http://dx.doi.org/10.1007/s11004-008-9204-2.
[14] Dentz M, Tartakovsky DM. Probability density functions for passive scalars dispersed in random velocity fields. Geophys Res Lett 2010;37:L24406. http://dx.doi.org/10.1029/2010GL045748.
[15] Meyer DW, Jenny P, Tchelepi HA. A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media. Water Resour Res 2010;46:W12522. http://dx.doi.org/10.1029/2010WR009450.
[16] Cirpka OA, de Barros FPJ, Chiogna G, Nowak W. Probability density function of steady state concentration in two-dimensional heterogeneous porous media. Water Resour Res 2011;47:W11523. http://dx.doi.org/10.1029/2011WR010750.
[17] Venturi D, Tartakovsky DM, Tartakovsky AM, Karniadakis GE. Exact PDF equations and closure approximations for advectivereactive transport. J Comput Phys 2013;243:32343. http://dx.doi.org/10.1016/j.jcp.2013.03.001.
[18] Suciu N. Diffusion in random velocity fields with applications to contaminant transport in groundwater. Adv Water Resour 2014;69:114133. http://dx.doi.org/10.1016/j.advwatres.2014.04.002.
[19] Suciu N, Radu FA, Attinger S, Schüler L, Knabner P. A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media. J Comput Appl Math 2015;289:241-252. http://dx.doi.org/10.1016/j.cam.2015.01.030.
[20] Suciu N, Schüler L, Attinger S, Vamoş C, Knabner P. Consistency issues in PDF methods. An Sti U Ovid Co-Mat 2015;23(3):187-208. http://dx.doi.org/10.1515/auom-2015-0055.
[21] Beckie R, Aldama AA, Wood EF. Modeling the large-scale dynamics of saturated groundwater flow using spatial filtering theory: 1. Theoreticald evelopment. Water Resour Res 1996;32(5):1269-1280. http://dx.doi.org/10.1029/96WR00276.
[22] Beckie, R, Aldama AA, Wood EF. Modeling the large-scale dynamics of saturated groundwater flow using spatial filtering theory: 2. Numerical Evaluation. Water Resour Res 1996;32(5):1281-1288. http://dx.doi.org/10.1029/96WR00277.
[23] Efendiev YR, Durlofsky LJ, Lee SH. Modeling of subgrid effects in coarse scale simulations of transport in heterogeneous porous media. Water Resour Res 2000;36:2031-2041. http://dx.doi.org/10.1029/2000WR900141
[24] Efendiev Y, Durlofsky LJ. A Generalized convection-diffusion model for subgrid transport in porous media. Multiscale Model Simul 2003;1(3):504-526. http://dx.doi.org/10.1137/S1540345902413693.
[25] Attinger S. Generalized Coarse Graining Procedures for Flow in Porous Media. 2003;7(4):253-273. http://dx.doi.org/10.1023/B:COMG.0000005243.73381.e3.
[26] Heße F, Radu FA, Thullner M, Attinger S. Upscaling of the advectiondiffusionreaction equation with Monod reaction. Adv Water Resour 2009;32:1336-1351. http://dx.doi.org/10.1016/j.advwatres.2009.05.009.
[27] Minier JP, Peirano E. The PDF approach to turbulent and polydispersed two-phase flows. Phys Rep 2001;352:1-214.
[28] Klimenko AY, Bilger RW. Conditional moment closure for turbulent combustion. Progr Energ Combust Sci 1999; 25:595-687. http://dx.doi.org/10.1016/S0360-1285(99)00006-4.
[29] Morales-Casique E, Neuman SP, Gaudagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: theoretical framework. Adv Water Resour 2006;29:1238-1255. http://dx.doi.org/10.1016/j.advwatres.2005.10.002.
[30] Morales-Casique E, Neuman SP, Gaudagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: computational analysis. Adv Water Resour 2006;29:1399-1418. http://dx.doi.org/10.1016/j.advwatres.2005.10.014.
[31] Pope SB. A Monte Carlo method for the PDF equations of turbulent reactive flow. Combustion Science and Technology 1981;25:159174. http://dx.doi.org/10.1080/00102208108547500.
[32] Möbus H, Gerlinger P, Brggemann D. Comparison of Eulerian and Lagrangian Monte Carlo PDF methods for turbulent diffusion flames. Combust Flame 2001;124:519-534. http://dx.doi.org/10.1016/S0010-2180(00)00207-8.
[33] Jones WP, Marquis AJ, Prasad VN. LES of a turbulent premixed swirl burner using the Eulerian stochastic field method. Combust Flame 2012;159:3079-3095. http://dx.doi.org/10.1016/j.combustflame.2012.04.008.
[34] Valiño L. A Field Monte Carlo Formulation for Calculating the Probability Density Function of a Single Scalar in a Turbulent Flow. Flow Turb Combust 1998;60(2):157-172. http://dx.doi.org/10.1023/A:1009968902446.
[35] Sabelnikov V, Soulard O. Rapidly decorrelating velocity-field model as a tool for solving one-point Fokker-Planck equations for probability density functions of turbulent reactive scalars. Phys Rev E 2005;72(1):016301. http://dx.doi.org/10.1103/PhysRevE.72.016301.
[36] Waclawczyk M, Pozorski J, Minier JP. New Molecular Transport Model for FDF/LES of Turbulence with Passive Scalar. Flow Turbulence Combust 2008;81:235-260. http://dx.doi.org/10.1007/s10494-007-9112-4.
[37] Wang H, Popov PP, Pope SB. Weak second-order splitting schemes for Lagrangian Monte Carlo particle methods for the composition PDF/FDF transport equations. J Comput Phys 2010;229:1852-1878. http://dx.doi.org/10.1016/j.jcp.2009.11.012.
[38] Kloeden PE, Platen E. Numerical solutions of stochastic differential equations. Berlin: Springer; 1999.
[39] Herz M. Mathematical modeling and analysis of electrolyte solutions. 2014; PhD thesis. http://www.mso.math.fau.de/fileadmin/am1/projects/PhD_Herz.pdf.
[40] Rubin Y, Sun A, Maxwell R, Bellin A. The concept of block effective macrodispersivity and a unified approach for grid-scale and plume-scale-dependent transport. J Fluid Mech 1999;395:161-180. http://dx.doi.org/10.1017/S0022112099005868.
[41] de Barros FPJ, Rubin Y. Modelling of block-scale macrodispersion as a random function. J Fluid Mech 2011;676:514-545. http://dx.doi.org/10.1017/jfm.2011.65.
[42] Bellin A, Tonina D. Probability density function of non-reactive solute concentration in heterogeneous porous formations, J. Contam. Hydrol. 2007;94:109-125. http://dx.doi.org/10.1016/j.jconhyd.2007.05.005.
[43] Minier J-P, Chibbaro S, Pope SB. Guidelines for the formulation of Lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent flows. Phys Fluids 2014;26:113303. http://dx.doi.org/10.1063/1.4901315.
[44] Klimenko AY. On simulating scalar transport by mixing between Lagrangian particles. Phys Fluids 2007;19:031702. http://dx.doi.org/10.1063/1.2711233
[45] Vamoş C, Suciu N, Vereecken H. Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J Comput Phys 2003;186(2):52744. http://dx.doi.org/10.1016/S0021-9991(03)00073-1.
[46] Suciu N, Radu FA, Prechtel A, Brunner F, Knabner P. A coupled finite element-global random walk approach to advection-dominated transport in porous media with random hydraulic conductivity. J Comput Appl Math 2013;24627-37. http://dx.doi.org/10.1016/j.cam.2012.06.027.
[47] Suciu N, Vamoş C, Vanderborght J, Hardelauf H, Vereecken H. Numerical investigations on ergodicity of solute transport in heterogeneous aquifers. Water Resour Res 2006;42:W04409. http://dx.doi.org/10.1029/2005WR004546.
[48] Kraichnan RH. Diffusion by a random velocity field. Phys Fluids 1970;13(1):2231. http://dx.doi.org/10.1063/1.1692799.
[49] Vamoş C, Crăciun M. Separation of components from a scale mixture of Gaussian white noises. Phys Rev E 2010;81:051125. http://dx.doi.org/10.1103/PhysRevE.81.051125.
[50] Suciu N, Vamos C. Ergodic estimations of upscaled coefficients for diffusion in random velocity fields. In: LEcuyér Pierre, Owen Art

B, editors. Monte Carlo and quasi-Monte Carlo methods 2008. Berlin: Springer; 2009. p. 61726. http://dx.doi.org/10.1007/978-3-642-04107-5_40.
[51] Dagan G. Upscaling of dispersion coefficients in transport through heterogeneous porous formations. In: Peters A et al., editors. Computational Methods in Water Resources X. Norwell, Mass: Kluwer Acad; 1994. p. 431-439.
[52] Schwarze H, Jaekel U, Vereecken H. Estimation of Macrodispersion by Different Approximation Methods for Flow and Transport in Randomly Heterogeneous Media. Transport Porous Med 2001;43:265-287. http://dx.doi.org/10.1023/A:1010771123844.
[53] Dentz M, Kinzelbach H, Attinger S, Kinzelbach W. Temporal behavior of a solute cloud in a heterogeneous porous medium. 3. Numerical simulations. Water Resour Res 2002;38:1118. http://dx.doi.org/10.1029/2001WR000436.
[54] Kurbanmuradov OA, Sabelfeld KK. Stochastic Flow Simulation and Particle Transport in a 2D Layer of Random Porous Medium Transp Porous Med 2010;85:347-373. http://dx.doi.org/10.1007/s11242-010-9567-y.
[55] Heße F, Prykhodko V, Schlter S, Attinger S. Generating random fields with a truncated power-law variogram: A comparison of several numerical methods. Environ Model Software 2014;55:32-48. http://dx.doi.org/10.1016/j.envsoft.2014.01.013.
[56] Bronstein IN, Semendjajew KA, Musiol G, Mühlig H. Taschenbuch der Mathematik. Frankfurt am Main: Verlag Harri Deutsch; 2006.

*Corresponding author.
Email address: suciu@math.fau.de (N. Suciu)

2016

Julia Robinson and Hilbert’s tenth problem

10 Jul at 7:51 am

Existence, localization and multiplicity of positive solutions to φ-Laplace equations and systems

April 20, 2022

On a robust Aitken-Newton method based on the Hermite polynomial

February 8, 2018

Towards a filtered density function approach for reactive transport in groundwater

Abstract

Authors

Keywords

Cite this paper as:

PDF

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

MR

ZBL

Google Scholar Profile

References

Paper in html format

Towards a filtered density function approach for groundwater

Abstract

1. Introduction

2. PDF/FDF equations

3. The Fokker-Planck approach

3.1. Direct Fokker-Planck approach

3.2. Reverse Fokker-Planck approach

3.3. Numerical solutions

4. Coarse-grained FDF simulations

4.1. A $P D F / F D F$ problem for passive transport in aquifers

4.2. Mixing models

4.3. Filtered velocity fields and upscaled diffusion coefficients

4.4. Coarse-grained simulations of transport

5. PDF/FDF solutions

6. Conclusions

Acknowledgements

Appendix A. Filtered Kraichnan fields

Appendix B. Self-averaging estimations of diffusion coefficients

Appendix C. Coarse-grained diffusion coefficients

Related Posts

Towards a filtered density function approach for reactive transport in groundwater

Abstract

Authors

Keywords

Cite this paper as:

PDF

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

MR

ZBL

Google Scholar Profile

References

Paper in html format

Towards a filtered density function approach for groundwater

Abstract

1. Introduction

2. PDF/FDF equations

3. The Fokker-Planck approach

3.1. Direct Fokker-Planck approach

3.2. Reverse Fokker-Planck approach

3.3. Numerical solutions

4. Coarse-grained FDF simulations

4.1. A P D F / F D F P D F / F D F PDF//FDFP D F / F D FPDF/FDF problem for passive transport in aquifers

4.2. Mixing models

4.3. Filtered velocity fields and upscaled diffusion coefficients

4.4. Coarse-grained simulations of transport

5. PDF/FDF solutions

6. Conclusions

Acknowledgements

Appendix A. Filtered Kraichnan fields

Appendix B. Self-averaging estimations of diffusion coefficients

Appendix C. Coarse-grained diffusion coefficients

Related Posts

4.1. A $P D F / F D F$ problem for passive transport in aquifers