Predicting the transport of groundwater contaminants 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 [51]. Hence, a risk analysis also includes the quantification of the uncertainty in order to evaluate how accurate the predictions are.

It is well known from papers published in the last decades that a major source of uncertainty associated with predicting contaminant concentrations is the lack of detailed information about the spatial heterogeneity of the hydraulic conductivity in the subsurface (see e.g. [23, 4]). A long standing approach to deal with this uncertainty is the stochastic parameterisation of the hydraulic conductivity through random space functions with statistics inferred from field and laboratory data. Via flow and transport equations, the contaminant concentrations being modelled are random functions too. Their statistics may be inferred from Monte Carlo ensembles of transport simulations, done for realisations of the random hydraulic conductivity, by stochastic perturbation approaches, or by the probability density function (PDF) method, which received an increased attention during the last decade [42].

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 [23]. But 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 the hydraulic conductivity has finite correlation lengths and the 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. 47]. 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. [12] 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 [2] proposed a mixing mechanism based on a time variable effective length scale which, in principle, could be determined experimentally. Kapoor and Gelhar [28, 27] 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 [29].

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 [51], can only be factored in by the so-called exceedance probability [e.g. 1]. 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 [10]. 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. Yee and Chan [52] analysed a large set of experimental data from tracer tests in the turbulent atmosphere. From this analysis they identified a collapse of higher-order concentration moments, which means that the higher-order moments can be expressed through lower-order moments. Srzic et al. [38] applied such an analysis to the transport of solutes in groundwater and also found a collapse. This might be a promising way of computing the concentration PDF from the mean concentration and the concentration variance. 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 [20, 21]. Cirpka et al. [7] extended the assumption of a $\beta$-shaped PDF to some special cases of reactive transport by mapping the statistics of conservative transport to those of mixing-controlled reactions. But Srzic et al. [39] concluded that $\beta$-shaped PDF’s only match the true PDF for low heterogeneities. By strictly assuming a multivariate Gaussian random velocity field and by assuming that the PDF of the centre of mass of the solute is also Gaussian, Dentz and Tartakovsky [14] derived a formula for the concentration PDF without need to solve a PDF evolution equation. Their approach consists of mapping the PDF of the random centre of mass of the plume, assumed to be Gaussian, onto the concentration PDF. A similar approach, using a stream function coordinate system, was further developed by Cirpka et al. [6]. A “Lagrangian concentration” framework for calculating the integrated PDF, the cumulative distribution function (CDF), was presented by de Barros and Fiori [11]. They derived a semi-analytical equation for the concentration CDF by assuming a small plume size and a normally distributed and statistically stationary conductivity field with low to mild heterogeneity. Compared to the $\beta$-CDF, both models perform similar. A review on assumed $\beta$-PDFs and mapping random variable approaches is given by Fiori et al. [19].

A more flexible approach - the PDF method - 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 which describes the transport of the PDF in the concentration space. 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 [33, 22, 43]. We refer to Celis and Figueira da Silva [5] for a recent review of mixing models.

The term “mixing” will be used from now on with a precise meaning, as in turbulence literature [5], which, although related, is different from its meaning in stochastic sub-surface hydrology, where it is associated to the effective dispersion coefficient [12]. The latter is the diffusion coefficient of the stochastic process modelling the transport, centred on the actual plume center of mass. It differs from the diffusion coefficient of the process centred on the centre of mass averaged over the ensemble of velocity realisations, i.e. the ensemble dispersion coefficient describing the spreading of the solute plume [12], by the diffusion coefficient of the centre of mass process [45, equation (3.3)]. The local dispersion coefficient, consisting of the molecular diffusion coefficient and a term describing the hydrodynamic dispersion by unresolved velocity fluctuations, enters additively into both the effective and the ensemble dispersion coefficients. All these are processes are in physical space. Instead, a mixing model generically consists of an advection and a diffusion in concentration space. The mixing models provide closures for the conditional average of the diffusion flux, determined by the local dispersion coefficients [42] and are related in this way to the transport processes in physical space.

Beyond approaches which overcome the need of mixing models, for instance in case of advective transport [46] or under simplifying hypotheses for stratified aquifers [37] and simplified solutions which only consider PDF transport in concentration space [3], Meyer et al. [31] proposed a solution to the full complexity of the PDF problem, in line with classical approaches in turbulence literature. In our recent publications, we derived consistency conditions in order to link PDF equations to Fokker-Planck equations for which efficient numerical schemes exist [40]. 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. Our findings also highlighted the inadequacy of classical mixing models used in turbulence for groundwater systems [42].

In our recent publications, we derived consistency conditions in order to link PDF equations to Fokker-Planck equations for which efficient numerical schemes exist [43, 40]. 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 [42].

In this paper, we investigate the 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 [42, 40]. 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 [36, 26, 15].

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 [33, 5]. In large eddy simulations (LES), the mixing time scale is often estimated as a velocity [15] or as a diffusion time scale [8]. Colucci et al. [8], 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^{\mathrm{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$ [33, 8].

For groundwater systems, characterised by the anisotropic local dispersion coefficients $D_{ij}$, Kapoor and Gelhar [28] 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|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 [42], 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 [26, 15]. Furthermore, Sabel’nikov et al. [36] 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 [2], 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}^{\mathrm{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}^{\mathrm{eff}}$ excludes this artificial dispersion and converges to $\mathbf{D}^{\mathrm{ens}}$ in the long-time limit for velocity fields with short range correlations [12, 45]. 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 [33, 8]. 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}^{\mathrm{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 [28], 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\to\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\text{\,}\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 [25, 32]. 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 [32], 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.

The TIEM model explicitly takes into account the dimensionality of the transport problem trough the expression (14) of the variance decay coefficient. This shows that by increasing the dimension of the physical space the variance decay coefficient increases, which results in enhanced mixing. The TIEM model also depends on the Péclet number implicitly, through the dependence on the effective diffusion coefficient given in (14). The limit of an infinite Péclet number causes a singularity in (14). However, for vanishing local dispersion there is no mixing at all and the PDF equation takes the form of a Fokker-Planck equation [49]. In case of passive transport as considered here, this Fokker-Planck equation reduces to equation (4) with the right hand side set to zero, which describes the PDF transport in physical space.

Unlike the approaches based on mapping random variables [14, 6, 11], the derivation of the PDF equation and the closure by IEM or TIEM mixing models are free of the low heterogeneity variance assumption. Another important difference is the way the influence of the sampling volume is taken into account. In mapping approaches [11], as well as in Monte Carlo simulations [39] the sampling volume is explicitly considered in the computation of the concentration. A sampling volume, associated to the spatial scale of the measurement, can be accounted for by a spatial filtering of the transport equations, similarly to the LES approach in turbulence [42]. The filtered density function (FDF), which describes the unresolved concentration fluctuations, verifies the PDF equation (4), with coefficients defined by spatial filtering and is solved by the GRW algorithm described in Section 4.2. The GRW solution also provides the filtered concentration, which corresponds to a spatial average over the sampling volume in this approach. In the limit of large filter widths, the filtered concentration tends to its ensemble average. For finite filter widths it is a random quantity. Its statistics can be inferred from a Monte Carlo ensemble of GRW-FDF solutions, obtained at low computational costs (e.g. computing time of the order of seconds [42]). The TIEM model for FDF simulations is easily obtained by adding the filter width at the denominators of terms in expression (14).

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).