Next Article in Journal
Development and Evaluation of the Combined Machine Learning Models for the Prediction of Dam Inflow
Previous Article in Journal
Environmental Benefits and Economical Sustainability of Urban Wastewater Reuse for Irrigation—A Cost-Benefit Analysis of an Existing Reuse Project in Puglia, Italy
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Theoretical Study about Ergodicity Issues in Predicting Contaminant Plume Evolution in Aquifers

School of Engineering, University of Basilicata, 85100 Potenza, Italy
Water 2020, 12(10), 2929; https://doi.org/10.3390/w12102929
Submission received: 31 August 2020 / Revised: 4 October 2020 / Accepted: 17 October 2020 / Published: 20 October 2020
(This article belongs to the Section Hydrology)

Abstract

:
A large-time Eulerian–Lagrangian stochastic approach is employed to: (1) estimate centroid position uncertainty of contaminant plumes that originate from instantaneous point sources in statistically stationary and isotropic porous formations; (2) assess the time needed for achieving ergodic conditions, which would allow for the evaluation of local concentration values based on the only ensemble mean distribution; (3) derive the concentration coefficient of variation (CV) as a function of asymptotic macro-dispersion coefficients and centroid trajectory variances. The results indicate that the decay time of plume position uncertainty is so large that there is practically no chance for effective ergodicity. The concentration coefficient of variation is zero at the centroid but rapidly increases when moving away from it. The dissipative effect of local dispersion in the presence of relatively high Péclet numbers is considerably exalted by marked flow field heterogeneity, which confirms the previously postulated synergic, non-additive effect of advection and local dispersion in passive solute dilution. A further result from this study is the derivation of the power law that relates dimensionless concentration micro-scale to dimensionless local dispersive area. The exponent of this power law is the same that appears in the relationship between dimensionless Kolmogorov turbulent micro-scale and flow Reynolds number.

1. Introduction

Protecting groundwater resources and attempting to limit the damages deriving from their sometimes unavoidable deterioration is one of the main goals in the field of environmental engineering. For that reason, many scientists in the last few decades have addressed, using different approaches, all the issues related to water flow and solute transport in heterogeneous porous formations (among the benchmark treatises: [1,2,3]). As a matter of fact, the marked heterogeneity of groundwater flow fields considerably complicates the theoretical analysis, tying it to the choice of a well-defined reference space–time scale and to the introduction of a number of simplifying assumptions. When affording the study of porous media by the classic continuum theory (which averages the microscopic balance equations over a representative elementary volume), one implicitly admits the existence of a lower-limit natural scale. As a consequence, the unknown microscopic characteristics of the medium are incorporated into macroscopic parameters that can more easily be estimated. However, as for the complex pore network, at a field scale the aquifer can be so heterogeneous that a detailed description of macro-units distribution is impossible (or at least impractical). Provided that it is not possible to analyze the related flow and transport processes in a mathematically exact form, an acceptable compromise is resorting to stochastic models that are able to consistently encompass their intrinsic uncertainty (e.g., [2,3,4,5,6]). The geometrical organization of the porous formation will then be considered only one of the infinite possible scenarios, all of them belonging to the same statistical population. The characterizing physical properties will be interpreted as space–time random functions, operatively described by ensemble mean and expected deviation. Nevertheless, the statistical methods are only a tool for handling phenomena that, while obeying basically deterministic laws, cannot be described at an acceptable level of detail. Calculating ensemble mean and standard deviation of a random variable provides an estimate of its most probable value and corresponding degree of reliability. The feasibility and the usefulness of such an operation depends on the occurrence of the conditions for the validity of the ergodic theorem (e.g., [7]). In simple words, the theorem postulates the possibility of assuming the equivalence of spatial or temporal means of a given variable calculated in a single realization (in the present case, the real aquifer) and theoretical point/instantaneous ensemble means of the same variable referring to the whole statistical population. For this to be possible, it is required that the spatial domain or the time interval used to compute the spatial/temporal means is much larger or longer than the spatial domain or the time interval over which the values of the given variable are correlated (or, equivalently, kept very close to each other). In the present study, the estimation of position and dimensions of contaminant plumes and corresponding point concentrations in heterogeneous and statistically isotropic saturated porous formations is addressed. In terms of position and dimensions, the discussion is based on centroid and central inertia moments, as well as on the conditions under which they can be considered as equivalent to same-order single-particle statistical moments; in terms of point concentrations, the discussion is based on concentration ensemble mean, variance and coefficient of variation (which in turn involve the single-particle moments as well as the statistics of the barycenter of mass) and the conditions under which the Gaussian ensemble mean for point instantaneous mass injection can be representative of the actual distribution. The mathematical treatment, which also hinges on previous author’s results, makes use of both Lagrangian (e.g., [8,9,10,11,12,13]) and Eulerian (e.g., [14,15,16]) framework for tracer transport in heterogeneous flow fields.

2. Methods

Let us consider steady flow and passive solute transport in statistically stationary and isotropic porous formations. After [17,18], the following relationship governs the evolution of the expected longitudinal inertia moment I11 of a plume that originates from an instantaneous point solute source at x = 0:
I 11 t = X 11 t S 11 t
where the angle brackets indicate the ensemble mean operator, t is the time, X11 the longitudinal single-particle trajectory variance, S11 the longitudinal centroid trajectory variance:
S 11 t = S 1 2 t
and S 1 the longitudinal deviatory centroid trajectory. For unit normalized total mass (M/n = 1, where M indicates solute mass and n medium porosity):
S i t = x i c x , t d x   i = 1 ,   2 ,   3
where x is the position vector and c x , t the actual concentration distribution. From the mass balance equation:
L A D c x , t = c t + v c D 2 c = 0
where LAD indicates the advection–dispersion operator, D the local dispersion coefficient and v the zero-divergence velocity field, one obtains:
d S 1 d t = d S 1 d t + d S 1 d t = V 1 + v 1 x c x , t d x
In Equation (5) and in what follows, V1 is the mean longitudinal velocity and v 1 the corresponding deviatory component:
v 1 = v 1 V 1
The underlying assumption for Equation (5) to be valid is that the concentration and its spatial derivatives vanish at large distances from origin:
c ,   c x i 0   x i ±   i = 1 ,   2 ,   3
Ensemble averaging leads to:
d S 1 d t = V 1
and, as a consequence:
d S 1 d t = v 1 x c x , t d x
or:
S 1 t = 0 t v 1 x c x , s d x d s
The centroid trajectory variance is then obtained by:
d S 1 2 d t = 2 S 1 d S 1 d t
and:
d S 1 2 d t = 2 0 t v 1 x v 1 y c x , t c y , s d x d y d s
Following the standard linearized transport formulation (e.g., [2]), the velocity field sampling trajectory consists of a mean-velocity straight line (a + Vt), with a indicating the generic particle starting position, only perturbed by the local dispersive Brownian displacement (XB). In these conditions:
c x , t = Ω 0 C 0 a δ x a V t X B t d a
where Ω0 indicates the injection volume, δ the Dirac Delta, and C0(a) the initial concentration distribution. As mentioned when introducing Equation (1), in the present study it is assumed that the initial concentration distribution is a point-like one: C 0 a δ a . Therefore:
d S 1 2 d t = 2 0 t v 1 x v 1 y δ x V t X B t δ y V s Y B s f X B Y B X B , Y B ; t , s , D d x d y d X B d Y B d s
where f X B Y B indicates the bivariate Normal probability density function of two different and statistically independent Brownian trajectories (which are also not statistically related to the Darcian flow field):
f X B Y B X B , Y B ; t , s , D = f X B X B , t , D f Y B Y B , s , D = 1 4 π D 3 t s 3 / 2 i = 1 3 exp X B i 2 4 D t exp Y B i 2 4 D s
From Equation (14), with C v 11 x y = v 1 x v 1 y indicating the longitudinal velocity stationary covariance, one obtains:
d S 1 2 d t = 2 0 t C v 11 x y f X B Y B x V t , y V s ; t , s , D d x d y d s
or:
S 1 2 = 0 t 0 t C v 11 x y f X B Y B x V t 1 , y V t 2 ; t 1 , t 2 , D d x d y d t 1 d t 2
Resorting to the Fourier transformation, one can write:
C v 11 x y = Z v 11 k exp j 2 π k x y d k
where Z v 11 k is the co-called spectrum of the longitudinal velocity and j = 1 . Substitution of Equation (18) into Equation (17) yields:
S 1 2 = 0 t 0 t Z v 11 k exp j 2 π k x y f X B Y B x V t 1 , y V t 2 ; t 1 , t 2 , D d x d y d k d t 1 d t 2
Finally, the Cartesian space–time double integrations performed accounting for Equation (15) lead to:
S 1 2 = Z v 11 k Γ k , t 2 π V k 2 + 16 π 4 D 2 k 4 d k
with:
Γ k , t = 1 + exp 8 π 2 D k 2 t 2 cos 2 π V k t exp 4 π 2 D k 2 t
and k = |k|.
In the present work, we derive S11 for Gaussian isotropic log-conductivity covariance function:
C Y r = σ Y 2 exp π r 2 4 I Y 2
where r = |r| is the scalar distance between two generic points of the domain, σ Y 2 is the log-conductivity variance and IY the log-conductivity integral scale. The corresponding spectrum is:
Z Y k = C Y r exp j 2 π k r d r = 8 σ Y 2 I Y 3 exp 4 π k 2 I Y 2
Note that considering a 3-D (three-dimensional) isotropic conductivity distribution makes sense in the presence of relatively reduced horizontal extension of the geologic units, so that the corresponding integral scale can be assumed of the same order of magnitude of the vertical one, which in turn is controlled by the stratigraphic structure.
The classic first-order theory of steady flow in stationary porous media (e.g., [2]) straightforwardly relates log-conductivity and velocity spectra by:
Z u 11 k = p = 1 3 q = 1 3 V p V q δ p 1 k 1 k p k 2 δ q 1 k 1 k q k 2 Z Y k
where δij indicates Kronecker Delta. Assuming that mean flow is directed along x1, with V = (V,0,0), λi = kiIY, τ = tV/IY, λ = |λ| and:
F 11 λ 1 , λ 2 , λ 3 = 1 λ 1 2 λ 2 2
one obtains:
S 1 2 σ Y 2 I Y 2 = 2 π 2 F 11 λ 1 , λ 2 , λ 3 exp 4 π λ 2 λ 1 2 + 4 π 2 P e 2 λ 4 1 + exp 8 π 2 P e λ 2 τ 2 cos 2 π λ 1 τ exp 4 π 2 P e λ 2 τ d λ
where Pe = VIY/D indicates the Péclet number. The triple integration in Equation (26) can more easily be afforded by switching to spherical coordinates:
λ 1 = λ sin θ cos φ λ 2 = λ sin θ sin φ λ 3 = λ cos θ
Thus:
S 1 2 σ Y 2 I Y 2 = 2 π 2 0 0 π 0 2 π 1 sin 2 θ cos 2 φ 2 exp 4 π λ 2 sin θ sin 2 θ cos 2 φ + 4 π 2 P e 2 λ 2 1 + exp 8 π 2 P e λ 2 τ 2 cos 2 π λ sin θ cos φ τ exp 4 π 2 P e λ 2 τ d φ d θ d λ
and, in terms of dimensionless time derivative:
d S 11 * d τ = 8 0 0 π 0 2 π 1 sin 2 θ cos 2 φ 2 exp 4 π λ 2 sin θ λ 2 4 π 2 λ 2 sin 2 θ cos 2 φ + 16 π 4 P e 2 λ 4 d Γ d τ d φ d θ d λ
where:
d Γ d τ = 8 π 2 P e λ 2 exp 8 π 2 P e λ 2 τ + 2 exp 8 π 2 P e λ 2 τ 2 π λ sin θ cos φ sin 2 π λ sin θ cos φ τ exp 4 π 2 P e λ 2 τ + 4 π 2 P e λ 2 cos 2 π λ sin θ cos φ τ exp 4 π 2 P e λ 2 τ
and S 11 * = S 11 / σ Y 2 I Y 2 . The evaluation of the integrals in Equation (29) is performed after some mathematical manipulation that transforms it into:
d S 11 * d τ = 16 0 0 π 0 2 π 1 sin 2   θ cos 2   φ 2 sin θ λ 2 0 τ exp 4 π 2 P e 2 τ τ + 4 π λ 2 cos 2 π λ sin θ   cos φ τ d τ d φ d θ d λ
Note that the equivalence of Equation (31) and Equations (29) and (30) can easily be verified by performing the integration over τ . At large τ, taking the integral principal value, one obtains:
0 τ exp 4 π 2 P e 2 τ τ + 4 π λ 2 cos 2 π λ sin θ cos φ τ d τ exp 4 π 2 P e τ + 4 π λ 2 0 τ cos 2 π λ sin θ cos φ τ d τ = exp 4 π 2 P e τ + 4 π λ 2 sin 2 π λ sin θ cos φ τ 2 π λ sin θ cos φ
and
d S 11 * d τ 8 π 0 0 π 0 2 π 1 sin 2 θ cos 2 φ 2 cos φ exp 4 π 2 P e τ + 4 π λ 2 sin 2 π λ sin θ cos φ τ λ d φ d θ d λ
Starting from Equation (33), we will examine two different cases: the behavior of the centroid large-time variance for highly advective regimes (Pe→∞) and for advective-diffusive regimes (finite Pe). In the first case, with τ / P e 1 :
d S 11 * d τ = 16 π τ 4 π 2 P e τ + 4 π 3 / 2 0 π 0 2 π 1 sin 2 θ cos 2 φ 2 sin θ exp π 2 sin 2 θ cos 2 φ 4 π 2 P e τ + 4 π τ 2 d φ d θ
Since, for τ→∞ and τ / P e 1 , when performing the integration over φ one can assume:
exp π 2 sin 2 θ cos 2 φ 4 π 2 P e τ + 4 π τ 2 π δ cos φ π sin θ τ 4 π 2 P e τ + 4 π 1 / 2 ,
the solution of Equation (34) turns out to be:
d S 11 * d τ 8 1 + π τ P e 2
and S 11 t 2 σ Y 2 V I Y t . Returning to Equation (1), with X 11 t = 2 σ Y 2 V I Y t (e.g., [2]):
I 11 t 2 σ Y 2 I Y V t 2 σ Y 2 I Y V t = 0
Thus, for zero local dispersion, an initial point injection is expected to remain a point solute body, randomly moving within the flow domain with high position uncertainty (S11X11). In the more realistic case that Pe has a finite value, at large times (4π2τ/Pe>>4π) one gets:
d S 11 * d τ = 2 P e 3 / 2 π 5 / 2 τ 1 / 2 0 π 0 2 π 1 sin 2 θ cos 2 φ 2 sin θ exp sin 2 θ cos 2 φ 4 P e τ d φ d θ
and:
S 11 * τ = 4 P e π 2 0 π 0 2 π 1 sin 2 θ cos 2 φ 2 cos φ Φ sin θ cos φ 2 P e τ d φ d θ const
with Φ indicating Error Function. Note that the subtractive constant in Equation (39) comes from the time integration in Equation (38) between zero and a time large enough for the asymptotic approximation of the integrand function to be valid. The large-τ integration over φ, with Φ→1, would invariably give zero, except for cosφ = 0 (which happens twice: for φ = π/2 and φ = 3π/2). In this case, the Error Function tends to 2 / π times its argument. The final result is therefore:
S 11 * τ 16 P e 3 / 2 π 5 / 2 τ const
and:
I 11 t 2 σ Y 2 I Y V t + 2 D t 16 π 5 / 2 σ Y 2 I Y 3 V 2 D 3 / 2 t + const
where 2 D t = X B 1 2 represents the contribution of the Brownian motion to X11. Incidentally, it can be shown that at short times, even for Pe ≠ ∞, S11 does not depend on it, is quadratic in τ and tends to X11. Equation (41) says that, in asymptotic conditions and for non-negligible local dispersion, a solute plume that originates from an instantaneous point source would see its expected longitudinal dimension to increase in time due to two different mechanisms: the first, represented by the term 2Dt, comes from pure local dispersion; the second, much more important, represented by a subtractive term that increases at a time rate lower than X11 (whose advection-related part, on the other hand, is not affected by local dispersion—see [13]), comes from the fundamental interplay of local dispersion and advection. Thus, in scalar transport and mixing in heterogeneous flow fields, the simple rule of effect superposition does not apply. Note that, in the case of anisotropic porous formations (represented by anisotropic Gaussian log-conductivity covariance function), Equation (40) would become:
S 11 * τ 16 P e L 3 / 2 e π 5 / 2 ε 3 / 2 τ const
with e = I Y v / I Y h indicating the anisotropy ratio (i.e., the ratio of vertical to horizontal log-conductivity integral scale), ε = D T / D L indicating the local dispersion anisotropy ratio (i.e., the ratio of transverse to longitudinal local dispersion coefficient) and P e L = U I Y h / D L the longitudinal Peclét number. Provided that both e and ε are typically < 1, the modified asymptotic two-particle covariance (Equation (42)) says that the more anisotropic the porous formation (i.e., the more layered-like), the less persistent the solute particle correlation. This makes perfect sense because, in the presence of non-negligible local dispersion, a pronounced stratification accelerates the sampling of markedly different velocity values. On the other hand, the more anisotropic the local dispersion (with the critical transverse coefficient that is much smaller than the longitudinal one), the weaker the dissipative effect and the slower the loss of particle correlation.
It will now be shown that the centroid variance coincides with the two-particle covariance (see also [15] for steady stream-flow transport):
S i i Θ i i = X i Y i
with X i = X i X i = X i V i t and Y i = Y i Y i = Y i V i t representing the deviatory longitudinal trajectory of two different particles. Indeed, from Equation (3):
S i 2 t = x i y i c x , t c y , t d x d y
where, at large times, for unit normalized total mass and instantaneous point injection at x = 0 (see [19]):
c x , t c y , t = 2 π 3 P   W W   P 1 / 2 exp 1 2 x V t T y V t T P   W W T   P 1 x V t y V t
with P = [Xii] and W = [Θii] indicating the diagonal tensors of one- and two-particle covariances respectively, T vector-matrix transpose and |A| = det A. Therefore:
S i i t = S i 2 t = S i 2 t S i t 2 = x i y i c x , t c y , t d x d y V i 2 t 2 = Θ i i t
Furthermore, in the regime of equilibrium between concentration large-scale irregularity production and small-scale irregularity dissipation, the following relationship is valid for the concentration variance, an important measure of concentration uncertainty:
σ c 2 x , t = c 2 x , t = i = 1 3 Θ i i t c x , t x i 2
Equation (47) can be obtained from the following concentration variance approximation (see [6]), which is valid for Gaussian and jointly Gaussian particle trajectories:
σ c 2 x , t = c x , t 2 c x , t 2 = 2 π 3 P   W W   P 1 / 2 exp 1 2 x V t T x V t T P   W W T   P 1 x V t x V t 2 π 3 / 2 P 1 / 2 exp 1 2 x V t T P 1 x V t 2 ,
by equating to zero L A D σ c 2 x , t with Θ i i / t 0 . The structure of Equation (47) is clearly similar to that involving turbulent velocity correlation and corresponding mean gradient after Prandtl:
v i v j = l i j 2 v i x j 2
where lij represents mixing length in the ij plane. Thus, one can conclude that the role of S i i Θ i i is that of solute transport (squared) mixing length in steady heterogeneous flow fields (that is, a sort of tracer (squared) micro-scale). Interestingly enough, while Kolmogorov universal equilibrium theory (e.g., [20]) gives the following order-of-magnitude relationship between the ratio of turbulence micro-scale η to flow domain length L and flow Reynolds number:
η L Re 3 / 4
the same power (−3/4) is found in the present study when analyzing the relationship between the ratio of square root of Θ 11 to average distance actually covered by the plume (Vt) and a dimensionless number that is obtained as the ratio of squared local dispersive length to squared log-conductivity integral scale. Indeed, from Equation (40):
Θ 11 V t 4 σ Y π 5 / 4 D t I Y 2 3 / 4
As a matter of fact, high Reynolds numbers pertain to flows that are increasingly chaotic and prone to the loss of correlation among fluid trajectories, exactly like subsurface tracer transport processes that are characterized by relatively high values of the ratio of local dispersive (diffusive-like) area to log-conductivity high-correlation area.
In order to complete the mathematical treatment, and to derive the concentration coefficient of variation, the transverse centroid variances/two-particle covariances (S22 and S33) as well as the transverse macro-dispersion coefficients (Dm22 = (1/2)dX22/dt and Dm33 = (1/2)dX33/dt) are derived as follows. First, rewriting Equation (26) for τ→∞ and i = 2,3 in spherical coordinates, one obtains:
S i i σ Y 2 I Y 2 = 2 π 2 0 0 π 0 2 π F i i λ , θ , φ exp 4 π λ 2 sin θ sin 2 θ cos 2 φ + 4 π 2 P e 2 λ 2 d φ d θ d λ
where:
F 22 λ , θ , φ = sin 4 θ sin 2 φ cos 2 φ
and:
F 33 λ , θ , φ = sin 2 θ cos 2 θ cos 2 φ
The integration over λ in Equation (52) is performed by taking the principal value:
0 exp 4 π λ 2 sin 2 θ cos 2 φ + 4 π 2 P e 2 λ 2 d λ 1 sin 2 θ cos 2 φ 0 exp 4 π λ 2 d λ = 1 4 sin 2 θ cos 2 φ
The final result is:
S 22 σ Y 2 I Y 2 = S 33 σ Y 2 I Y 2 = 2 3 π
The asymptotic transverse macro-dispersion coefficients are derived after [21] obtaining:
D m i i D = σ Y 2 F i i λ 1 , λ 2 , λ 3 exp 4 π λ 2 λ 2 λ 1 2 + 4 π 2 P e 2 λ 4 d λ + 1
where:
F 22 λ 1 , λ 2 , λ 3 = λ 1 2 λ 2 2 λ 4
and:
F 33 λ 1 , λ 2 , λ 3 = λ 1 2 λ 3 2 λ 4
Switching to spherical coordinates and performing the triple integration, one gets:
D m 22 D = D m 33 D = σ Y 2 3 + 1
The asymptotic longitudinal macro-dispersion coefficient for Gaussian log-conductivity covariance and finite Péclet was already computed by [13]: as mentioned above when talking about X 11 = 2 0 t D m 11 τ d τ and Equation (41), its advective part does not differ from the high-Péclet case. Thus:
D m 11 D = σ Y 2 P e + 1
Finally, the concentration coefficient of variation:
C V x , t = σ c x , t c x , t
with σ c x , t = σ c 2 x , t , is obtained at large times from Equation (47) and the Gaussian ensemble mean concentration, which is the solution of the constant–coefficient advection–dispersion equation for instantaneous point source at x = 0:
c x , t = i = 1 3 1 4 π D m i i t exp x i V δ i 1 t 2 4 D m i i t
Accounting for Equations (40), (56), (60) and (61), one gets:
C V x ˜ , τ 4 σ Y 2 P e 3 / 2 π 5 / 2 σ Y 2 + 1 P e 2 τ 3 / 2 x ˜ 1 τ 2 + σ Y 2 P e 2 x ˜ 2 2 + x ˜ 3 2 3 π σ Y 2 3 + 1 2 τ 2
with x ˜ i = x i / I Y .

3. Results

One of the main issues when dealing with the estimation of location, extension and peak concentration values of contaminant plumes in heterogeneous subsurface flows concerns the predictive ability of the ergodic spatial moments ( X i t and X i i t ) and the corresponding Gaussian concentration ensemble mean. As a matter of fact, when the conditions for the validity of the ergodic theorem were achieved, and the statistics of the single solute particle could be substituted by the same-order spatial moments referred to the single-realization plume, one would have: S i i / S i 0 ,   S D I i i / I i i 0 ,   S i S i X i , I i i I i i X i i and the space–time distribution of the concentration in case of instantaneous point source at x = 0 could be estimated by Equation (63) at a good level of approximation. Therefore, before analyzing the large-time behavior of CV at a given short distance from the centroid and at a fixed point in space, it is useful to evaluate the time needed for S 11   (the macroscopic flow direction is obviously the more critical one) to decrease to some small percentage ω of X 11 = 2 σ Y 2 V I Y t + 2 D t :
τ ω = t ω V I Y = 0.209 P e 3 ω 2 1 + 1 σ Y 2 P e 2
Table 1 and Table 2 respectively list, for σ Y 2 = 0.1 and σ Y 2 = 0.5 , τω computed from Equation (65) when ω = 0.1, 0.05 and 0.01, and Pe = 10, 100 and 500.
As one can see, except for Pe = 10 (a rather uncommon value of Péclet deriving from particularly intense local dispersion, slow flow and reduced log-conductivity integral scale), τω is characterized by very high values, even when one only requires that ω = 0.1 and the porous formation is mildly heterogeneous ( σ Y 2 < 1 ). Consider, for example, τω for σ Y 2 = 0.1 and Pe = 100: 1.727 × 107. That means that, based on the large-time approximation for both X11 and S11, 1.727 × 107 times the ratio I Y / V is the time needed for a 90% reduction of the ratio S11/X11 (which still cannot be considered enough to represent really ergodic conditions). Note that IY is usually in the order of meters and V is frequently in the order of 10−6 to 10−5 meters/second. Therefore, one can conclude that there is practically no chance to achieve the theoretical ergodic domain where c x , t = c x , t . The best one can do is therefore to evaluate CV and eventually to estimate c x , t as c x , t plus a suitable multiple of the standard deviation σ c x , t . By way of example, Figure 1, Figure 2 and Figure 3 respectively show, for σ Y 2 = 1, 0.5 and 0.1, the large-time behavior of CV at a relatively short distance ( r = r = 3 I Y ) from the expected moving centroid and peak of the Gaussian distribution (63) S t = V t , 0 , 0 , when Pe = 10, 100 and 500.
Note that based on the large-time concentration variance in Equation (47), the concentration coefficient of variation is zero at S t and increases with increasing distance from it. That means that the reliability of the concentration estimates deteriorates as their value decreases. As one can see, for given log-conductivity variance, CV is an increasing function of Pe. On the other hand, the three figures reveal that, whereas for diffusive or mildly advective regimes (Pe = 10 and Pe = 100), the level of heterogeneity expressed by σ Y 2 only slightly affects CV, in the case of markedly advective regimes (Pe = 500), a more heterogeneous log-conductivity distribution ( σ Y 2 = 1 ) leads to faster CV reduction and, therefore, faster solute plume dilution. From a quantitative standpoint, Figure 3 (the worst scenario) says that a weakly heterogeneous formation in a markedly advective transport regime could require a travel distance of up to 10,000 integral scales to allow a standard deviation equal to 20% of the (relatively close to the peak) expected concentration value. For IY = 1 m, that would mean 10 km: that is, a travel distance hardly compatible with the typical dimensions of a hydraulically continuous aquifer and a residence time much larger than the typical times of technical interest. Finally, Figure 4, Figure 5 and Figure 6 show the large-time behavior of CV at a fixed point in space (x1 = 5000IY, x2 = 1IY, x3 = 1IY). Such a type of estimate (a sort of “Eulerian” CV) could be useful, as an example, when choosing the best location of agricultural or urban supply wells in the presence of possible upstream groundwater pollution sources.
As one can see, the concentration coefficient of variation experiences its minimum in correspondence of the centroid passage (τ = 5000). Overall, the values of CV at different times are always very high, except for very small Péclet. Obviously, it should be observed that the expected concentration quickly goes to zero far from the peak and the evaluation of the corresponding CV makes sense only if the toxicity threshold is very low. As for CV at a fixed distance from the moving centroid (the “Lagrangian” CV), we see that advective regimes are characterized by lower uncertainty when coupled with higher heterogeneity.
Thus, one can conclude that in transport processes that are dominated by advective mechanisms, the dissipative effect of local dispersion is exalted by flow field heterogeneity. This fully confirms what was postulated by [6] in terms of synergetic non-additive effect of advection and local dispersion in tracer dilution.

4. Discussion and Conclusions

The protection of groundwater resources, and the prediction of the time evolution of their accidental contamination, is one of the priorities in the field of environmental engineering. Many researchers in the last few decades have addressed theoretical, computational and experimental issues related to water flow and solute transport in heterogeneous porous formations. The importance of theoretical studies in this field (against the costs of extensive field surveys and the computational burden of realistic numerical simulations) resides in the possibility to design better targeted, smaller-scale interventions. Note that groundwater transport theories are typically applied to the saturated zone of aquifers. As a matter of fact, within the vadose zone, the lack of hydraulic continuity implies the absence of a longitudinal gradient. Fluid displacement only occurs along the vertical direction due to gravity and suction and usually ends (rather quickly) at the aquifer phreatic surface (or at the impervious boundary in the presence of semi-dry layers). One could say that infiltration across the unsaturated zone can eventually be viewed as the pre-initial condition for saturated transport. Moreover, the ubiquitous asymptotic character of the analytical formulations would not be consistent with the reduced temporal horizon of contaminant motion driven by gravity and capillary forces. Finally, the dangerousness of pollution events is mostly related to downstream propagation and its typical space–time uncertainty.
The present study makes use of a combined Eulerian–Lagrangian stochastic approach to evaluate centroid position uncertainty in the case of passive contaminant plumes originating from an instantaneous point source in statistically stationary and isotropic saturated porous formations, to assess the time needed for achieving the so-called ergodic conditions (which would allow for the evaluation of point concentrations based on the only Gaussian ensemble mean distribution), and to derive the large-time concentration coefficient of variation as a function of asymptotic macro-dispersion coefficients and centroid trajectory variances. The Eulerian part of the present formulation consists of making use of the advection–dispersion equation to derive the governing equation of the plume centroid second-order statistics and the equilibrium concentration variance approximation. The Lagrangian part consists of expressing the actual concentration distribution and, as a consequence, all its moments in terms of solute particle actual trajectories and related statistics. The results indicate that, even in mildly heterogeneous aquifers and mildly advective regimes, there is practically no chance to achieve the theoretical ergodic domain. Therefore, estimating the concentration coefficient of variation is essential for the prediction of point levels of exposure to toxic non-reactive substances accidentally carried out by groundwater, even at large times after contamination. The large-time coefficient of variation is zero at the moving centroid. At (short) fixed distances from it, CV is an increasing function of Péclet for given log-conductivity variance, and decays in time at a rate that is higher for more heterogeneous log-conductivity distributions. The concentration coefficient of variation evaluated at a fixed point in space experiences its minimum in correspondence of the centroid passage, but rapidly increases when the peak of the distribution moves away. Its values are always very large, except for unrealistically small Péclet. Obviously, it must be stressed that the expected concentration quickly goes to zero far from the peak of the distribution and the evaluation of the corresponding CV makes sense if the toxicity threshold is very low (which is the case for several contaminant chemicals). For both CV at a fixed distance from the moving centroid (the “Lagrangian” CV) and CV at a fixed point in space (the “Eulerian” CV), we see that the advective regimes (high Péclet) are characterized by lower uncertainty when coupled with higher heterogeneity. Thus, one can conclude that in the presence of relatively high Péclet numbers (transport processes dominated by advection), the dissipative effect of a relatively weak local dispersion is enhanced by large values of log-conductivity (and velocity) variance. As a matter of fact, large log-conductivity (and velocity) variance imply highly irregular concentration distributions and, as a consequence, large mass fluxes based on Fick’s law and more efficient homogenization processes. This fully confirms what was previously postulated by the author in terms of synergetic non-additive effect of advection and local dispersion in passive solute dilution. In the present study, it has also been shown that the ratio of the square root of the two-particle covariance (a sort of passive scalar micro-scale) to the average distance actually covered by the plume is proportional to the power −3/4 of the ratio of squared local dispersive length to squared log-conductivity integral scale. This is the same type of power law that governs the relationship between dimensionless turbulent Kolmogorov’s micro-scale and Reynolds number. As a matter of fact, momentum and mass transfer by chaotic flows are related by Reynolds’ analogy. The square root of the two-particle covariance can be considered as the tracer transport microscale and, therefore, the equivalent of Kolmogorov microscale η. The ratio of Kolmogorov microscale to the representative flow domain dimension L scales with the power −3/4 of Reynolds number. Reynolds number measures the relative magnitude of inertial and viscous effects, which can respectively be associated with the tendency of the flow to chaos and fast-decaying trajectories correlation, and to order and long-range trajectories correlation. As mass transport is an intrinsically unsteady process, the equivalent of L must be represented by the travel distance Vt, which measures the portion of the domain actually sampled by the solute body. Finally, the role of Reynolds number must be played by the ratio of two quantities that respectively symbolize uncorrelated and correlated particle displacements: the transverse area of the local dispersive spot about the average longitudinal trajectory (proportional to Dt) and the correlated area (proportional to IY2).
Finally, a further issue when dealing with real porous formations may be represented by their confinement by natural or artificial boundaries. In mathematical terms, confinement is equivalent to the exclusion of the contribution from the wave-number region around k = (0,0,0) (that is, the contribution from the largest scales of heterogeneity). As a consequence, in this case, log-conductivity and velocity could be considered as periodic functions to be expanded in Fourier series with a finite maximum period. With Fvii(k) indicating the discrete spectrum of the pair (vi, vi), one obtains at t→∞, from Equation (20) rewritten for the generic i:
S i 2 = Θ i i = Z v i i k Γ k , t 2 π V k 2 + 16 π 4 D 2 k 4 d k
the following asymptotic two-particle covariance:
S i 2 = Θ i i = k 0 F v i i k 2 π V k 2 + 16 π 4 D 2 k 4 = const
Asymptotically constant two-particle covariance is synonym of large-time complete dilution, which can be associated with maximum system entropy as for any confined domain (see [6] for subsurface flow transport and [22] for stream-flow suspended load transport).

Funding

This research received no external funding.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Bear, J. Dynamics of Fluids in Porous Media; Dover Publications Inc.: New York, NY, USA, 1972; 764p, ISBN 0486656756. [Google Scholar]
  2. Dagan, G. Flow and Transport in Porous Formations; Springer: New York, NY, USA, 1989; 465p, ISBN 3540510982. [Google Scholar]
  3. Rubin, Y. Applied Stochastic hydrogeology; Oxford University Press: London, UK, 2003; 416p, ISBN 9780198031543. [Google Scholar]
  4. Kapoor, V.; Gelhar, L.W. Transport in three-dimensionally heterogeneous aquifers: 1. Dynamics of concentration fluctuations. Water Resour. Res. 1994, 30, 1775–1788. [Google Scholar] [CrossRef]
  5. Kapoor, V.; Gelhar, L.W. Transport in three-dimensionally heterogeneous aquifers: 2. Predictions and observations of concentration fluctuations. Water Resour. Res. 1994, 30, 1789–1801. [Google Scholar] [CrossRef]
  6. Pannone, M.; Kitanidis, P.K. Large-time behavior of concentration variance and dilution in heterogeneous formations. Water Resour. Res. 1999, 35, 623–634. [Google Scholar] [CrossRef]
  7. Walters, P. An Introduction to Ergodic Theory; Springer: New York, NY, USA, 1982; 250p, ISBN 9780387951522. [Google Scholar]
  8. Pannone, M. Transient Hydrodynamic Dispersion in Rough Open Channels: Theoretical Analysis of Bed-Form Effects. J. Hydraul. Eng. ASCE 2010, 136, 155–164. [Google Scholar] [CrossRef]
  9. Pannone, M. Effect of nonlocal transverse mixing on river flows dispersion: A numerical study. Water Resour. Res. 2010, 46, W08534. [Google Scholar] [CrossRef]
  10. Pannone, M. Longitudinal Dispersion in River Flows Characterized by Random Large-Scale Bed Irregularities: First-Order Analytical Solution. J. Hydraul. Eng. ASCE 2012, 138, 400–411. [Google Scholar] [CrossRef]
  11. Pannone, M.; De Vincenzo, A.; Brancati, F. A Mathematical Model for the Flow Resistance and the Related Hydrodynamic Dispersion Induced by River Dunes. J. Appl. Math. 2013, 432610. [Google Scholar] [CrossRef] [Green Version]
  12. Pannone, M.; De Vincenzo, A. Stochastic numerical analysis of anomalous longitudinal dispersion and dilution in shallow decelerating stream flows. Stoch. Environ. Res. Risk Assess. 2015, 29, 2087–2100. [Google Scholar] [CrossRef]
  13. Pannone, M. An Analytical Model of Fickian and Non-Fickian Dispersion in Evolving-Scale Log-Conductivity Distributions. Water 2017, 9, 751. [Google Scholar] [CrossRef] [Green Version]
  14. Pannone, M. On the exact analytical solution for the spatial moments of the cross-sectional average concentration in open channel flows. Water Resour. Res. 2012, 48, W08511. [Google Scholar] [CrossRef]
  15. Pannone, M. Predictability of tracer dilution in large open channel flows: Analytical solution for the coefficient of variation of the depth-averaged concentration. Water Resour. Res. 2014, 50, 2617–2635. [Google Scholar] [CrossRef]
  16. Pannone, M.; Mirauda, D.; De Vincenzo, A.; Molino, B. Longitudinal Dispersion in Straight Open Channels: Anomalous Breakthrough Curves and First-Order Analytical Solution for the Depth-Averaged Concentration. Water 2018, 10, 478. [Google Scholar] [CrossRef] [Green Version]
  17. Kitanidis, P.K. Prediction by the method of moments of transport in heterogeneous formation. J. Hydrol. 1988, 102, 453–473. [Google Scholar] [CrossRef]
  18. Dagan, G. Transport in heterogeneous porous formations: Spatial moments, ergodicity, and effective dispersion. Water Resour. Res. 1990, 26, 1281–1290. [Google Scholar] [CrossRef]
  19. Pannone, M.; Kitanidis, P.K. Large-time spatial covariance of concentration of conservative solute and application to the Cape Cod tracer test. Transp. Porous Media 2001, 42, 109–132. [Google Scholar] [CrossRef]
  20. Tennekes, H.; Lumley, J.L. A First Course in Turbulence; MIT Press: Boston, MA, USA, 1972; 320p, ISBN 9780262200196. [Google Scholar]
  21. Gelhar, L.W.; Axness, C.L. Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resour. Res. 1983, 19, 161–180. [Google Scholar] [CrossRef]
  22. Mirauda, D.; De Vincenzo, A.; Pannone, M. Simplified entropic model for the evaluation of suspended load concentration. Water 2018, 10, 378. [Google Scholar] [CrossRef] [Green Version]
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure 1. Concentration coefficient of variation at a given distance from the moving centroid σ Y 2 = 1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Figure 1. Concentration coefficient of variation at a given distance from the moving centroid σ Y 2 = 1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Water 12 02929 g001
Figure 2. Concentration coefficient of variation at a given distance from the moving centroid σ Y 2 = 0.5 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Figure 2. Concentration coefficient of variation at a given distance from the moving centroid σ Y 2 = 0.5 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Water 12 02929 g002
Figure 3. Concentration coefficient of variation at a given distance from the moving centroid σ Y 2 = 0.1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Figure 3. Concentration coefficient of variation at a given distance from the moving centroid σ Y 2 = 0.1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Water 12 02929 g003
Figure 4. Concentration coefficient of variation at a fixed point in space σ Y 2 = 1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Figure 4. Concentration coefficient of variation at a fixed point in space σ Y 2 = 1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Water 12 02929 g004
Figure 5. Concentration coefficient of variation at a fixed point in space σ Y 2 = 0.5 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Figure 5. Concentration coefficient of variation at a fixed point in space σ Y 2 = 0.5 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Water 12 02929 g005
Figure 6. Concentration coefficient of variation at a fixed point in space σ Y 2 = 0.1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Figure 6. Concentration coefficient of variation at a fixed point in space σ Y 2 = 0.1 . CV is the concentration coefficient of variation; τ is the dimensionless time.
Water 12 02929 g006
Table 1. Decay time of the ratio of two-particle to one-particle covariance. Pe: Péclet number; ω: decay percentage.
Table 1. Decay time of the ratio of two-particle to one-particle covariance. Pe: Péclet number; ω: decay percentage.
ωPe10100500
0.15.225 × 1031.727 × 1072.511 × 109
0.052.09 × 1046.909 × 1071.004 × 1010
0.015.225 × 1051.727 × 1092.511 × 1011
Table 2. Decay time of the ratio of two-particle to one-particle covariance. Pe: Péclet number; ω: decay percentage.
Table 2. Decay time of the ratio of two-particle to one-particle covariance. Pe: Péclet number; ω: decay percentage.
ωPe10100500
0.11.451 × 1042.009 × 1072.592 × 109
0.055.806 × 1048.035 × 1071.037 × 1010
0.011.451 × 1062.009 × 1092.592 × 1011

Share and Cite

MDPI and ACS Style

Pannone, M. A Theoretical Study about Ergodicity Issues in Predicting Contaminant Plume Evolution in Aquifers. Water 2020, 12, 2929. https://doi.org/10.3390/w12102929

AMA Style

Pannone M. A Theoretical Study about Ergodicity Issues in Predicting Contaminant Plume Evolution in Aquifers. Water. 2020; 12(10):2929. https://doi.org/10.3390/w12102929

Chicago/Turabian Style

Pannone, Marilena. 2020. "A Theoretical Study about Ergodicity Issues in Predicting Contaminant Plume Evolution in Aquifers" Water 12, no. 10: 2929. https://doi.org/10.3390/w12102929

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop