Predicting the transport of groundwater contaminations remains a demanding task, especially with respect to the heterogeneity of the subsurface, the large measurement uncertainties, and the increasing impact of human activities on groundwater systems [1]. Hence, a risk analysis also includes the quantification of the uncertainty in order to evaluate how accurate the predictions are.

Transport of dissolved contaminants in the subsurface is strongly influenced by the flow velocity of the groundwater, which in turn is determined by the properties of the surrounding geological structures. Their properties like the hydraulic conductivity are generally highly heterogeneous on many different scales. These heterogeneities range from the order of magnitude of individual grains to large geological structures like facies, fractures and sediment layers.

Spatial fluctuations of the aquifer properties make the transport of the contaminants highly heterogeneous too. Water parcels transporting a contaminant and travelling very closely together can be separated and follow different and distinct flow paths. As a consequence, an enhanced spreading of the plume is observed. The impact of the heterogeneities on the transport behaviour is well known in general [2,3]. But in order to predict the transport of contaminants accurately, it is necessary to know the heterogeneous aquifer properties influencing the transport everywhere. Monitoring the complete aquifer on all relevant scales is not feasible, but it is possible to retrieve partial knowledge by local measurements. These measurements can be used to generate a geostatistical representation of the aquifer. This way, the remaining uncertainty of aquifer properties and model parameters can be taken into account. For applying a stochastic framework, an ensemble of aquifer realisations is generated in accordance with the geostatistical representation. Hydraulic properties are modelled as spatial random functions, which in turn leads to the contaminant concentrations being modelled as spatial random functions too.

The moment approach uses transport equations of the concentration moments consistent with the geostatistical representation of the aquifer’s heterogeneity. If this heterogeneity is statistically homogeneous, the equation for the first moment, which is the mean concentration, has the following characteristics: The highly heterogeneous and spatially fluctuating groundwater velocity is replaced by an ensemble averaged velocity field and the effect of the fluctuating velocity on the transport is modelled by an enhanced dispersion called macrodispersion or ensemble dispersion [2]. This approach has the limitation that the ensemble averaged concentration only describes the mean plume behaviour. In general, the mean behaviour differs from that of a specific plume in a single aquifer. See Figure 1 for a comparison between the mean concentration and a concentration obtained from a simulation in a specific velocity field realisation. Only if such a single plume has sampled a representative part of the aquifer, it becomes ergodic and its transport behaviour can be modelled by the ensemble averaged behaviour, described above. In a first step, possible deviations from the mean behaviour can be quantified by the concentration variance. It is transported by the same processes as the mean concentration, thus it is advected by

the averaged velocity field and dispersed by an enhanced macrodispersion. But concentration variance is also generated by mean concentration gradients and simultaneously it is destroyed by dissipative processes, which are created by small-scale fluctuations in the velocity field. In order to calculate the influence of these small-scale fluctuations on the concentration variance, a so-called closure model is needed.

In the field of turbulence modelling, where very similar transport equations are used, different approaches exist for such closure models [e.g. 4]. Up to this point, the adoption of these approaches to groundwater transport modelling has been hampered by the vastly different flow conditions prevalent in both fields. Contrary to most other problems where turbulent flows are more challenging, the roles are reversed here. The strong mixing induced by turbulent flows causes this closure problem to be easier to tackle. The mixing induced by heterogeneities in the groundwater flow is slower and changes significantly in time and is therefore more difficult to model. As Dentz et al. [5] have shown, the mechanism which generates the physical mixing in a given aquifer realisation is more reliably described by the effective dispersion coefficients in comparison to the ensemble dispersion coefficients, which correspond to the turbulent diffusion coefficients. The effective dispersion is small at early times and increases only slowly with time. Therefore, concentration gradients at early times are steep and may remain steep for prolonged times, which in turn prevents the smoothing of concentration fluctuations and preserves concentration uncertainty. Andricevic [6] proposed a mixing mechanism based on a time variable effective length scale which, in principle, could be determined experimentally. Kapoor and Gelhar $[7,8]$ derived a transport equation for the concentration variance, including local
dispersivity and macrodispersive transport. By neglecting the local dispersivity, the results from Dagan [9] could be derived. But it was concluded that even very small local dispersivities create a qualitatively different behaviour compared to the zero local dispersivity case, as the local dispersivity is the only mechanism which can reduce the variance. They used an approach developed for turbulent flows to model the variance dissipation, created by the local dispersivity. Furthermore, analytical solutions for the long-time behaviour of the concentration variance were derived. These results were confirmed for globally integrated variances by numerical simulations [10].

If the predictions made by a contaminant transport model are to be used for risk analysis, even more information than the mean concentration and the variance is needed. Risk thresholds, regulated e.g. by an environmental agency [1], can only be factored in by the so-called exceedance probability [e.g. 11]. It depends on the complete one-point probability density function (PDF) of the concentration. The concentration variance, as discussed above, can only be used to calculate a first estimation of the exceedance probability [12]. Even if such an estimation would be an acceptable approximation, rare events or extreme values cannot be mapped by the mean concentration and the concentration variance alone. This limitation stems from the fact that by using only the first two statistical moments, namely the mean and the variance, a Gaussian shape for the concentration PDF is implied. Such a distribution is short-tailed and therefore excludes the possibility of rare events. The first studies applying such a PDF framework to the transport in groundwater used a beta distribution, fully characterised by two parameters, to fit the concentration PDF in a non-Gaussian way [13]. But Srzic et al. [14] concluded that beta-shaped PDF’s only match the true PDF for low heterogeneities. As a consequence, a second approach - the PDF approach - is investigated in this study. This approach yields an equation for the whole PDF of the concentration and thus makes no assumptions about the shape of it. The crux of these PDF methods is finding a mixing model. We refer to Celis and Figueira da Silva [15] for a recent review of mixing models. Meyer et al. [16] simulated the concentration PDF for advection-dispersion processes tailored to the transport in groundwater without assuming a specific shape. The PDF transport equation was solved by particle simulations developed for turbulent flows. A need for better mixing models was identified $[17,18,16]$. In our recent publications, we derived consistency conditions in order to link PDF equations to Fokker-Planck equations for which efficient numerical schemes exist [19, 20]. We derived mixing models from simulated concentration time series. These mixing models performed well, but a trajectory needs to be prescribed on which the concentration is sampled. Prescribing such a trajectory is only possible in special cases. In addition, we used spatially filtered probability density functions (FDF) to further reduce the computational costs [21].

Another major advantage of the PDF approach is the possibility to include mass transfer, like chemical reactions or radioactive decay, even in case of nonlinear reactions. This intriguing property of the PDF approach is possible by assuming that the mass transfer solely depends on the concentration [17, 18, 20].

In this paper, we investigate the spatially resolved concentration variance and concentration PDF behaviour over a long time period and show that mixing models used before fail to reproduce the variance at all times. To that end, we first introduce a closed transport equation for the one-point concentration PDF in section 2. Using this equation as a starting point, we derive the transport equations of the first two statistical moments. We show that by prescribing a certain mixing model, both the PDF and the variance equations have the same closure problem and thus depend on the same mixing coefficients. The importance of this finding lies in the possibility of testing new mixing models with much simpler concentration variance simulations first and subsequently transferring them to PDF models. Next, we present analytical solutions of the moment equations and show the dependence of the analytical solution of the concentration variance on the mixing model. In section 3 we identify the need for a time dependency of the coefficients of the mixing model and we propose a new time dependent mixing model. This new model is explicit and in a closed form. It is then verified in section 4 by comparing the previously derived analytical solution of the concentration variance equation with the old and the new mixing model with reference Monte Carlo simulations. Afterwards, the new model is also used in the PDF framework and compared to reference Monte Carlo solutions. Finally, we conclude our work in section 5.

The non-reactive transport of a solute in groundwater can be described by

where $C(\mathbf{x},t)$ is the concentration of the solute which is transported by the statistically homogeneous random velocity field $\mathbf{V}(\mathbf{x})$ and the local dispersion $\mathbf{D}$. It is assumed to be diagonal with $D_{11}=D_{L}$ being the longitudinal component, $D_{ii}=D_{T}$ for $i>1$ being the transversal components, and $D_{ij}=0$ for $i\neq j$. Throughout this work, we will be using the Einstein summation convention. The stochastic partial differential equation (1) describes the time evolution of a plume of a dissolved substance in the groundwater. If an ensemble of statistically equivalent solutions of this equation is calculated, the mean behaviour can be calculated from the ensemble average $\langle C\rangle=\sum_{i=1}^{N}C_{i}/N$ over $N$ realisations. As a first measure, the variance $\sigma_{c}^{2}(\mathbf{x},t)$ can quantify how good the ensemble average approximates the behaviour in a specific realisation. The mean concentration $\langle C\rangle(\mathbf{x},t)$ and the concentration variance $\sigma_{c}^{2}(\mathbf{x},t)$ are the first and second statistical moments of the one-point one-time concentration PDF $P(c;\mathbf{x},t)$ :

Thus, if a transport equation for the PDF is derived, transport equations of the mean and variance can be derived too.

The IEM model describes the decrease of the concentration PDF too slow, as we have already pointed out $[19,21]$. It was developed for simulating turbulent flows. One major difference between turbulent flows and flows in porous media is the time scale on which mixing takes place. In the turbulent regime, it is often taken as a constant. And even there, a mixing time scale as a variable parameter has already been taken into account [27, 28, 29].

The original IEM model for turbulent flows approximates the conditional dissipation rate by equation (5). The variance decay coefficient $\chi$ is proportional to the inverse mixing time scale. In classical PDF approaches, the latter is usually assumed to be proportional to the turbulence time scale [17, 15]. In large eddy simulations (LES), the mixing time scale is often estimated as a velocity [29] or as a diffusion time scale [24]. Colucci et al. [24], for instance, used the subgrid length scale $\Delta_{l}$, which defines the transition from resolved to unresolved scales and the subgrid diffusion coefficient, which corresponds to an isotropic ensemble dispersion coefficient $D^{\text{ens }}$ in groundwater flows. Following this approach, we can formulate the variance decay coefficient as

The dimensionless model parameter $k_{\chi}$ is usually in the range of $1\leq k_{\chi}\leq 3$ [17, 24].

For groundwater systems, characterised by the anisotropic local dispersion coefficients $D_{ij}$, Kapoor and Gelhar [7] arrived at a very similar equation for the variance decay coefficient, by introducing the Taylor microscales $\Delta_{c_{i}}$, which characterise the gradients of the concentration fluctuations along the ith coordinate. The resulting variance decay coefficient is

But the Taylor microscales could only be fitted to measurements, as a closed formula was not given.

We should recall that the IEM model was developed to approximate the second derivative of the conditional dissipation rate with respect to $c(5)$. The conditional dissipation rate is defined by

Here, $\langle A\mid B\rangle=\langle AB\rangle/\langle B\rangle$ denotes the conditional expectation of $A$ given $B$. The conditional dissipation rate $\mathcal{M}$ depends on the squared concentration gradients, which clearly evolve in time. But the IEM model has no way of accounting for this evolution. It only takes the difference between the current concentration and the local mean concentration into account. As we have already observed [21], a more accurate mixing model would account for larger dissipation rates at early
times and smaller dissipation rates at later times, as the concentration gradients decrease. For turbulent reactive flows, a dependence of $\chi$ on the Reynolds number of the subgrid scale flow was already proposed [28, 29]. Furthermore, Sabel’nikov et al. [27] modelled the mixing frequency as a stochastic process in order to account for the entire range of time scales in the mixing process. They named their model extended interaction by exchange with the mean (EIEM). We elaborate the idea of using a variable variance decay coefficient $\chi(t)$ and propose a new time dependent extension of the model, adapted to the transport processes in groundwater.

Like Andricevic [6], we approximate the concentration gradients using an evolving effective spatial scale $\lambda_{c}(t)$. This assumption implies that the squared concentration gradients evolve inversely proportional to that squared characteristic length scale $(\nabla C)^{2}\sim\lambda_{c}(t)^{-2}$ as the plume spreads and its fringes and fluctuations smooth out. Since the concentration fluctuations are smoothed out by the local dispersion with the characteristic scale $\lambda_{c}(t)=\sqrt{2Dt}$, we assume a decay of the conditional dissipation rate (12) to be proportional to $\mathcal{M}\sim\lambda_{c}(t)^{-2}$. In order to improve the IEM model, we include this dependency on $\lambda_{c}(t)$ into the variance decay coefficient (10).

Furthermore, the ensemble dispersion coefficient $\mathbf{D}^{\text{ens }}$ accounts for an artificial dispersion which is caused by centre of mass fluctuations of the solute plume from realisation to realisation. The effective dispersion coefficient $\mathbf{D}^{\text{eff }}$ excludes this artificial dispersion and converges to $\mathbf{D}^{\text{ens }}$ in the long-time limit for velocity fields with short range correlations [30, 5]. Because the mixing in turbulent flows is so much faster than it is in groundwater flows, the difference does not matter for turbulent flows. Therefore, the mathematically simpler to handle ensemble dispersion coefficient is used in studies concerning flows in the turbulent regime [17, 24]. But in groundwater flows the difference is significant and because the centre of mass fluctuations do not influence the dissipation, the effective dispersion coefficient $\mathbf{D}^{\text{eff }}$ describes the correct behaviour for the mixing model. With these physical arguments and choosing $k_{\chi}=2$ from the middle of the interval of reported values, we propose following time-dependent variance decay coefficient:

In order to show the similarities between this newly proposed variance decay coefficient and the previous ones, we assume an isotropic correlation length of the $\log$ conductivity field $\lambda$ and an isotropic local dispersion coefficient $D$, which gives us the dispersive time scale $\tau_{D}=\lambda^{2}/D$. With this relationship, we can transform the variance decay coefficient to

Equation (14) generalises equation (10) to a time-variable characteristic length

scale and to anisotropic effective dispersion coefficients. Compared to the coefficient (11) introduced by Kapoor and Gelhar [7], the new variance decay coefficient (14) depends on the effective dispersion coefficients instead of the local dispersion coefficients. Furthermore, the unclosed Taylor microscale was replaced by the correlation length and a dimensionless time factor $\tau_{D}/t$ was added.

As shown in figure 2, this variance decay coefficient has larger values than the constant one at early times, which causes a stronger dissipation. But then it drops below the constant variance decay coefficient and approaches $\lim_{t\rightarrow\infty}\chi(t)=0$. In order to distinguish this model from other extensions of the IEM model, we name it the time dependent interaction by exchange with the mean model (TIEM).

With this extension, the simplicity and low computational costs of the IEM model are preserved, while at the same time, it incorporates the time dependent physical processes causing the dissipation.

This paper presents a new and time-dependent mixing model: an extended IEM model for groundwater, named TIEM. We showed that the same mixing model is used for both the concentration variance evolution equation (7) and the concentration PDF evolution equation (4). This link was used to verify the new TIEM model (13) with the much simpler to handle variance equation. The verification was done by comparing an analytical solution of the variance equation (9), which depends on a mixing closure model, to two independent numerical models. The TIEM model shows a strong improvement over the IEM model. A significant deviation from the reference simulations can only be observed at times $t<15\mathrm{~d}$. And even for these very short times, the new model is a significant improvement.

Based on these promising results, the model was transferred to the PDF framework. The results obtained from the PDF simulations with the TIEM model are not quite as satisfying as the results from the variance simulations
mentioned above. Although there are mismatches at early and also at later times, the new model is still a major improvement over the classical IEM model. One possible way of further improving the IEM model is to derive a partial differential equation as a dynamic model for the variance decay coefficient [39, 40]. Such a model could include the actual and instantaneous length scales of the processes destroying the variance. This feature would make it possible to also apply the model to statistically non-homogeneous conductivity fields [40], as needed if the fields are to be conditioned on measurements.

The GRW-simulations, together with the TIEM model, can easily be extended to three-dimensional problems, to anisotropic dispersion coefficients, and to reactive transport. Especially the latter point is worth highlighting. The reaction terms can simply be plugged into the PDF equations, which makes the PDF framework the method of choice for modelling complex reactive transport in groundwater.

We kindly thank Veronika Zeiner, Malte Schüler, Jannis Weimar, and Falk Heße for their helpful suggestions and proof reading and David Schäfer for his help with parts of the numerical code. 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 as part of the H-DuR Project (funding number: 02 E 11062 B).

The transport equation for the mean concentration can be derived from the PDF transport equation (4) by multiplying it with $c$ and integrating over the entire concentration space. Doing so yields

The order of integration and derivation is swapped on the left hand side and on the right hand side the product rule is applied:

On the left hand side, the definition of the mean concentration $\langle C\rangle:=\int cP\mathrm{~d}c$ can already be inserted and on the right hand side, the integral is evaluated:

Both terms on the right hand side vanish, the second one is obvious, but the first one needs further comment. The case $c=0$ is clear, but for $c=1$, the term could potentially result in a non-zero value, if all concentration is gathered at one singular point as a Dirac delta function. But this case is excluded as it is not relevant if studying natural systems. Thus, the transport equation for the mean concentration is

It can be seen that the mixing term does not influence the mean concentration, as it cancels itself out.

The variance is defined by (3). Thus, as with the derivation of the mean concentration, we start again from the PDF transport equation (4), but now we multiply it with $c^{2}$ and integrate over the whole concentration space. The order of integration and derivation is swapped and the product rule is applied to the right hand side:

The first term on the right hand side vanishes for the same reason as in the derivation of the mean concentration. The definition of the mean concentration (2) can be inserted into the second term on the right hand side:

The term inside the squared brackets on the right hand side could already be replaced by the concentration variance (3), but to do this for every term, the transport equation of $\langle C\rangle^{2}$ needs to be subtracted from equation (A.6). Therefore, the equation of the squared mean concentration needs to be derived first. This is done by multiplying equation (A.4) by $\langle C\rangle$, yielding

By making extensive use of the product rule we arrive at

The dispersion term is further modified by using the product rule:

Hence, equation (A.8) is transformed into

If we compare the second line of equation (A.10) with equation (A.7), we see that it vanishes and the transport equation of $\langle C\rangle^{2}$ is

Now, equation (A.11) can be subtracted from equation (A.6):

Finally, the definition of the concentration variance (3) is inserted, which yields the transport equation for the variance:

An analytical solution of the variance transport equation (7) will now be derived. As this equation is a linear inhomogeneous partial differential equation, a fundamental solution is derived and convolved with the inhomogeneity of equation (7) in order to derive the general solution. This way, an analytical solution can be found without making any approximations or assumptions. This derivation is similar to the one presented by Kapoor and Gelhar [8], but we believe the derivation presented here is easier to comprehend and we include a time dependency of the variance decay coefficient $\chi(t)$.

If we define the differential operator $L(\mathbf{x},t)$ by

with the inhomogeneity

then we can rewrite the concentration variance transport equation (7) as

The fundamental solution (or Green’s function) $G\left(\mathbf{x}-\mathbf{x}^{\prime},t,t^{\prime}\right)$ is defined as the solution of the differential operator $L(\mathbf{x},t)$ with delta functions as the inhomogeneity:

The fundamental solution is translation invariant in space, because the operator $L$ has constant coefficients with respect to $\mathbf{x}$. The general solution of equation (7) is given by

where $\sigma_{c_{h}}^{2}(\mathbf{x},t)$ is the solution of equation (7) without the inhomogeneity. But because we assume that the initial condition is known exactly, the variance at time $t=0$ is $\sigma_{c}^{2}(\mathbf{x},t=0)=0$. Thus, without the inhomogeneity, which acts as the only source term, the solution of the homogeneous partial differential equation is $\sigma_{c_{h}}^{2}(\mathbf{x},t)=0$ for all times. Therefore, the solution of the homogeneous equation can be dropped.

If Green’s function is known, the solution of equation (7) can be calculated from equation (B.5), which is a convolution of Green’s function and the inhomogeneity in physical space:

where $\mathcal{F}^{-1}$ denotes the inverse Fourier transform and a tilde denotes the Fourier transform of a function. Hence, $\tilde{G}$ and $\tilde{g}$ need to be calculated in order to obtain the solution $\sigma_{c}^{2}$.

Fourier transforming both sides of equation (B.4) gives an inhomogeneous ordinary differential equation in the frequency domain:

with $I$ being the imaginary unit. This ordinary differential equation can be solved by separation of variables for a solution of the homogeneous equation. But with the initial condition of $\sigma_{c}^{2}(\mathbf{x},t=0)=0$ as explained above and thus $G\left(\mathbf{x}-\mathbf{x}^{\prime},t=0,t^{\prime}\right)=0$, the homogeneous solution always stays zero, as the inhomogeneity is the only source term. Nevertheless, the homogeneous solution can be used as a starting point to guess a particular solution of the inhomogeneous solution, yielding

where $\Theta$ is the Heaviside step function. In order to transform the inhomogeneity (B.2), the transformed mean concentration $\langle C\rangle$ from equation (8) needs to be plugged in:

At this point, the time shift $t_{0}$ is needed. Otherwise, a singularity for $t=0$ would cause problems, as the Gaussian distribution would tend to a Dirac delta function for small times. The Fourier transformed mean concentration is

With this solution, the Fourier transformed inhomogeneity can be calculated:

Finally, the transformed Green’s function (B.8) and the transformed inhomogeneity (B.11) are inserted into equation (B.6):

By completing the square for the variable $\mathbf{k}$, the Fourier integrand is lead back to a Gaussian function which can be integrated. Because the ensemble dispersion tensor is diagonal, we can simplify the expression by only considering the diagonal elements. Now, only a final time integral remains to be calculated:

This integral can either be evaluated analytically by using a long-time approximation or by applying numerical methods.

[1] WWAP, The United Nations World Water Development Report 4: Managing Water under Uncertainty and Risk, no. Vol. 1 in World Water Assessment Programme, Unesco, Paris, 2012.
[2] L. W. Gelhar, C. L. Axness, Three Dimensional Stochastic Analysis of Macrodispersion in Aquifers, Water Resour. Res. 19 (1) (1983) 161-180. doi:10.1029/WR019i001p00161.
[3] D. T. Burr, E. A. Sudicky, R. L. Naff, Nonreactive and reactive solute transport in three-dimensional heterogeneous porous media: Mean displacement, plume spreading, and uncertainty, Water Resour. Res. 30 (3) (1994) 791815. doi:10.1029/93WR02946.
[4] H. Tennekes, J. L. Lumley, A First Course in Turbulence, MIT Press, Cambridge, Mass., 1972.
[5] M. Dentz, H. Kinzelbach, S. Attinger, W. Kinzelbach, Temporal behavior of a solute cloud in a heterogeneous porous medium 1 .Point-like injection, Water Resour. Res. 36 (12) (2000) 3591-3604. doi:10.1029/ 2000WR900162.
[6] R. Andričević, Effects of local dispersion and sampling volume on the evolution of concentration fluctuations in aquifers, Water Resour. Res. 34 (5) (1998) 1115-1129. doi:10.1029/98WR00260.
[7] V. Kapoor, L. W. Gelhar, Transport in three-dimensionally heterogeneous aquifers: 1. Dynamics of concentration fluctuations, Water Resour. Res. 30 (6) (1994) 1775-1788. doi:10.1029/94WR00076.
[8] V. Kapoor, L. W. Gelhar, Transport in three-dimensionally heterogeneous aquifers 2. Predictions and observations of concentration fluctuations, Water Resour. Res. 30 (6) (1994) 1789-1801. doi:10.1029/94WR00075.
[9] G. Dagan, Stochastic Modeling of Groundwater Flow by Unconditional and Conditional Probabilities 1.Conditional Simulation and the Direct Problem, Water Resour. Res. 18 (4) (1982) 813-833. doi:10.1029/ WR018i004p00835.
[10] V. Kapoor, P. K. Kitanidis, Advection-diffusion in spatially random flows: Formulation of concentration covariance, Stoch. Hydrol. Hydraul. 11 (5) (1997) 397-422. doi:10.1007/BF02427926.
[11] R. Andričević, V. Cvetković, Evaluation of Risk from Contaminants Migrating by Groundwater, Water Resour. Res. 32 (3) (1996) 611-621. doi:10.1029/95WR03530.
[12] F. P. J. de Barros, A. Fiori, A. Bellin, A simple closed-form solution for assessing concentration uncertainty, Water Resour. Res. 47 (12) (2011) 1-5. doi:10.1029/2011WR011107.
[13] V. Fiorotto, E. Caroni, Solute concentration statistics in heterogeneous aquifers for finite Peclet values, Transp. Porous Media 48 (3) (2002) 331351. doi:10.1023/A:1015744421033.
[14] V. Srzic, V. Cvetkovic, R. Andricevic, H. Gotovac, Impact of aquifer heterogeneity structure and local-scale dispersion on solute concentration uncertainty: Impact of Aquifer Heterogeneity on Concentration Uncertainty, Water Resour. Res. 49 (6) (2013) 3712-3728. doi:10.1002/wrcr. 20314.
[15] C. Celis, L. F. Figueira da Silva, Lagrangian Mixing Models for Turbulent Combustion: Review and Prospects, Flow, Turbul. Combust. 94 (3) (2015) 643-689. doi:10.1007/s10494-015-9597-1.
[16] D. W. Meyer, P. Jenny, H. A. Tchelepi, A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media, Water Resour. Res. 46 (12) (2010) 1-17. doi:10.1029/2010WR009450.
[17] S. B. Pope, PDF Methods for Turbulent Reactive Flows, Prog. Energy Combust. Sci. 11 (2) (1985) 119-192. doi:10.1016/0360-1285(85) 90002-4.
[18] R. O. Fox, Computational Models for Turbulent Reacting Flows, Cambridge Series in Chemical Engineering, Cambridge University Press, New York, 2003.
[19] N. Suciu, F. A. Radu, S. Attinger, L. Schüler, P. Knabner, A Fokker-Planck approach for probability distributions of species concentrations transported in heterogeneous media, J. Comput. Appl. Math. 289 (2015) 241-252. doi: 10.1016/j.cam.2015.01.030.
[20] N. Suciu, L. Schüler, S. Attinger, C. Vamos, P. Knabner, Consistency issues in PDF methods, An. St. Univ. Ovidius Constanta, Ser. Mat. 23 (3) (2015) 187-208. doi:10.1515/auom-2015-0055.
[21] N. Suciu, L. Schüler, S. Attinger, P. Knabner, Towards a filtered density function approach for reactive transport in groundwater, Adv. Water Resour.Accepted.
[22] J. Villermaux, J. C. Devillon, Représentation de la coalescence et de la redispersion des domaines de ségrégation dans un fluide par un modèle d’interaction phénoménologique., in: Proceedings of the 2nd International symposium on chemical reaction engineering, Elsevier New York, 1972, pp. 1-13.
[23] C. Dopazo, E. E. O’Brien, An approach to autoignition of a turbulent mixture, Acta Astronaut. 1 (1974) 1239-1266. doi:10.1016/0094-5765(74) 90050-2.
[24] P. J. Colucci, F. A. Jaberi, P. Givi, S. B. Pope, Filtered density function for large eddy simulation of turbulent reacting flows, Phys. Fluids 10 (2) (1998) 499-515. doi:10.1063/1.869537.
[25] V. Raman, H. Pitsch, A consistent LES/filtered-density function formulation for the simulation of turbulent flames with detailed chemistry, Proc. Combust. Inst. 31 (2) (2007) 1711-1719. doi:10.1016/j.proci.2006.07. 152.
[26] P. P. Popov, S. B. Pope, Implicit and explicit schemes for mass consistency preservation in hybrid particle/finite-volume algorithms for turbulent reactive flows, J. Comput. Phys. 257 (2014) 352-373. doi:10.1016/j.jcp. 2013.09.005.
[27] V. Sabel’nikov, M. Gorokhovski, N. Baricault, The extended IEM mixing model in the framework of the composition PDF approach: applications to diesel spray combustion, Combust. Theory Modell. 10 (1) (2006) 155-169. doi:10.1080/13647830500348109.
[28] W. Jones, A. Marquis, V. Prasad, LES of a turbulent premixed swirl burner using the Eulerian stochastic field method, Combust. Flame 159 (10) (2012) 3079-3095. doi:10.1016/j.combustflame.2012.04.008.
[29] I. A. Dodoulas, S. Navarro-Martinez, Large Eddy Simulation of Premixed Turbulent Flames Using the Probability Density Function Approach, Flow, Turbul. Combust. 90 (3) (2013) 645-678. doi:10.1007/ s10494-013-9446-z.
[30] N. Suciu, Diffusion in random velocity fields with applications to contaminant transport in groundwater, Adv. Water Resour. 69 (2014) 114-133. doi:10.1016/j.advwatres.2014.04.002.
[31] M. Dentz, H. Kinzelbach, S. Attinger, W. Kinzelbach, Temporal behavior of a solute cloud in a heterogeneous porous medium 3. Numerical simulations, Water Resour. Res. 38 (7) (2002) 23-1-23-13. doi: $10.1029/2001$ WR 000436.
[32] C. Vamoş, N. Suciu, H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comput. Phys. 186 (2) (2003) 527-544. doi:10.1016/S0021-9991(03)00073-1.
[33] R. H. Kraichnan, Diffusion by a Random Velocity Field, Phys. Fluids 13 (1) (1970) 22-31. doi:10.1063/1.1692799.
[34] F. Heße, V. Prykhod’ko, S. Schlüter, S. Attinger, Generating random fields with a truncated power-law variogram. A comparison of several numerical methods with respect to accurary and reproduction of structural features., Environ. Model. Softw. 55 (2014) 32-48. doi:10.1016/j.envsoft.2014. 01.013.
[35] J. P. Eberhard, N. Suciu, C. Vamoş, On the self-averaging of dispersion for transport in quasi-periodic random media, J. Phys. A: Math. Gen. 40 (4) (2007) 597. doi:10.1088/1751-8113/40/4/002. URL http://iopscience.iop.org/1751-8121/40/4/002
[36] I. T. Drummond, S. Duane, R. R. Horgan, Scalar diffusion in simulated helical turbulence with molecular diffusivity, J. Fluid Mech. 138 (1984) 75-91. doi:10.1017/S0022112084000045.
[37] N. Suciu, C. Vamoş, J. Vanderborght, H. Hardelauf, H. Vereecken, Numerical Investigations on Ergodicity of Solute Transport in Heterogeneous Aquifers, Water Resour. Res. 42 (4) (2006) 1-17. doi:10.1029/ 2005WR004546.
[38] N. Suciu, F. A. Radu, A. Prechtel, F. Brunner, P. Knabner, A coupled finite element-global random walk approach to advection-dominated transport in porous media with random hydraulic conductivity, J. Comput. Appl. Math. 246 (2013) 27-37. doi:10.1016/j.cam.2012.06.027.
[39] H. G. Im, T. S. Lund, J. H. Ferziger, Large eddy simulation of turbulent front propagation with dynamic subgrid models, Phys. Fluids 9 (12) (1997) 3826-3833. doi:10.1063/1.869517.
URL http://scitation.aip.org/content/aip/journal/pof2/9/12/ 10.1063/1.869517
[40] C. D. Pierce, P. Moin, A dynamic model for subgrid-scale variance and dissipation rate of a conserved scalar, Phys. Fluids 10 (12) (1998) 3041-3044. doi:10.1063/1.869832.
URL http://scitation.aip.org/content/aip/journal/pof2/10/12/ 10.1063/1.869832