Theoretical investigation of nonergodic solute dispersion in natural porous formations characterized by persistent and antipersistent power-law log-conductivity correlations

Abstract

Tracer dispersion within log-conductivity fields represented by power-law semivariograms is investigated by an analytical first-order Lagrangian approach that, in treating subsurface flow and transport, resorts to the superposition principle of an infinite double hierarchy of mutually independent scales of heterogeneity. The results of the investigation, which are corroborated by a preliminary field validation, and also interpreted in terms of probabilistic risk analysis, say that transport anomaly is intrinsically associated with evolving-scale heterogeneous porous formations, regardless of their semivariogram scaling exponent b. In contrast with what was found by previous studies that dealt with this subject in a Lagrangian framework, it is demonstrated that: (1) the magnitude of nonergodic dispersion is nonmonotonically related to b; (2) consistently assuming a characteristic advective-heterogeneity length-scale leads to a universal (and quadratic) dependence of the dimensionless macrodispersion coefficients on the dimensionless time. Additionally, it is demonstrated that, in the presence of fractal heterogeneity, and unlike what happens for short-range correlations, diffusion acts as an antagonist mechanism in terms of Fickian dispersion achievement. Finally, the reinterpretation of antipersistent and persistent correlations as a double hierarchy of oscillatory nonperiodic and periodic fields, respectively, besides allowing for a technical explanation of all the detected trends, envisions a possible alternative methodology for their numerical generation.

Résumé

La dispersion des traceurs dans les champs de conductivité logarithmique représentés par des semivariogrammes en loi de puissance est étudiée par une approche lagrangienne analytique du premier ordre qui, en traitant l’écoulement et le transport souterrain, recourt au principe de superposition d’une double hiérarchie infinie d’échelles d’hétérogénéité indépendantes les unes des autres. Les résultats de l’étude, qui sont corroborés par une validation préliminaire sur le terrain et également interprétés en termes d’analyse de risque probabiliste, indiquent que l’anomalie de transport est intrinsèquement associée aux formations poreuses hétérogènes à échelle évolutive, quel que soit l’exposant d’échelle b de leur semivariogramme. Contrairement à ce qui a été constaté dans les études précédentes qui ont traité ce sujet dans un cadre lagrangien, il est démontré que: (1) l’amplitude de la dispersion non ergodique est liée de manière non monotone à b; (2) l’hypothèse cohérente d’une échelle de longueur caractéristique de l’advection et de l’hétérogénéité conduit à une dépendance universelle (et quadratique) des coefficients de macro-dispersion adimensionnels sur le temps adimensionnel. En outre, il est démontré qu’en présence d’hétérogénéité fractale, et contrairement à ce qui se passe pour les corrélations à courte portée, la diffusion agit comme un mécanisme antagoniste en termes de réalisation de la dispersion de Fick. Enfin, la réinterprétation des corrélations antipersistantes et persistantes comme une double hiérarchie de champs oscillatoires non périodiques et périodiques, respectivement, en plus de permettre une explication technique de toutes les tendances détectées, envisage une méthodologie alternative possible pour leur génération numérique.

Resumen

La dispersión del trazador dentro de campos de conductividad logarítmica representados por semivariogramas de ley de potencias se investiga mediante un enfoque lagrangiano analítico de primer orden que, al tratar el flujo y el transporte subsuperficiales, recurre al principio de superposición de una doble jerarquía infinita de escalas de heterogeneidad mutuamente independientes. Los resultados de la investigación, corroborados por una validación preliminar sobre el terreno e interpretados también en términos de análisis probabilístico de riesgos, afirman que la anomalía de transporte está intrínsecamente asociada a formaciones porosas heterogéneas de escala evolutiva, independientemente de su exponente de escala b de semivariograma. En contraste con lo hallado por estudios anteriores que trataban este tema en un marco lagrangiano, se demuestra que: (1) la magnitud de la dispersión no ergódica está relacionada de forma no monotónica con b; (2) asumir de forma consistente una escala de longitud advectivo-heterogénea característica conduce a una dependencia universal (y cuadrática) de los coeficientes adimensionales de macrodispersión con respecto al tiempo adimensional. Adicionalmente, se demuestra que, en presencia de heterogeneidad fractal, y a diferencia de lo que ocurre para correlaciones de corto alcance, la difusión actúa como un mecanismo antagonista en términos de logro de dispersión Fickiana. Finalmente, la reinterpretación de las correlaciones antipersistentes y persistentes como una doble jerarquía de campos oscilatorios no periódicos y periódicos, respectivamente, además de permitir una explicación técnica de todas las tendencias detectadas, vislumbra una posible metodología alternativa para su generación numérica.

摘要

对具有幂律半变异函数表示的对数传导性场内的示踪扩散进行了解析的一阶拉格朗日方法研究，该方法在处理地下流动和运移时，采用了无限双层次异质性互相独立尺度的叠加原理。研究结果通过初步的现场验证，并通过概率风险分析进行解释，得出结论:无论其半变异函数标度指数b的值取何大小，运移异常与不断演化尺度异质多孔介质形成本质上存在联系。与以前在拉格朗日框架中处理这个问题的研究所发现的相反，它证明了以下几点:(1)非遍历性扩散的幅度与b的关系非单调；(2)一致地假设一个特征对流异质性长度尺度导致了无量纲宏观弥散系数与无量纲时间之间的普遍(且二次)依赖关系。此外，它证明了在存在分形异质性的情况下，与短程相关不同，扩散在Fickian弥散的实现方面起到了对抗性机制的作用。最后，将非连续性和连续性相关解释为一种双层次振荡非周期和周期场的方法，除了能够技术上解释所有检测到的趋势外，还展望了一种可能的替代方法来生成它们的结果。

Resumo

A dispersão do traçador dentro dos campos de log-condutividade representados por semivariogramas de lei de potência é investigada por uma abordagem Lagrangiana analítica de primeira ordem que, ao tratar o fluxo e transporte subterrâneo, recorre ao princípio da superposição de uma hierarquia dupla infinita de escalas de heterogeneidade mutuamente independentes. Os resultados da investigação, que são corroborados por uma validação preliminar de campo, e também interpretados em termos de análise de risco probabilística, dizem que a anomalia de transporte está intrinsecamente associada a formações porosas heterogêneas de escala evolutiva, independentemente de seu expoente de escala de semivariograma b. Em contraste com o que foi encontrado por estudos anteriores que trataram deste assunto em uma estrutura Lagrangiana, demonstra-se que: (1) a magnitude da dispersão não ergódica é não monotonicamente relacionada a b; (2) assumir consistentemente uma escala de comprimento de heterogeneidade advectiva característica leva a uma dependência universal (e quadrática) dos coeficientes de macrodispersão adimensional no tempo adimensional. Adicionalmente, demonstra-se que, na presença de heterogeneidade fractal, e ao contrário do que ocorre para correlações de curto alcance, a difusão atua como um mecanismo antagonista em termos de alcance da dispersão Fickiana. Finalmente, a reinterpretação das correlações antipersistentes e persistentes como uma dupla hierarquia de campos oscilatórios não periódicos e periódicos, respectivamente, além de permitir uma explicação técnica de todas as tendências detectadas, vislumbra uma possível metodologia alternativa para sua geração numérica.

Introduction

Due to the great economic and environmental value of groundwater resources, many scientists in the last few decades (among the benchmark classic studies: Bear 1972; Dagan 1989; Rubin 2003) have focused their efforts on addressing the issues related to water flow and solute transport in natural porous formations. The marked structural and functional heterogeneity of these formations typically ties the mathematical treatment of the topic to the choice of a well-defined reference space-time scale, and to the introduction of several more or less drastic simplifications. The classic continuum theory (which averages the microscopic balance equations over a representative elementary volume) implies the existence of a lower-limit natural scale, and the possibility to incorporate the smaller-range irregularities into macroscopic parameters (sort of constitutive variables) that can more easily be estimated. Unfortunately, as for the complex pore network, the aquifer can be so heterogeneous at field scale that an exact description of the distribution of the geologic macro-units is impossible (or at least impractical). An acceptable compromise between the need for faithfully describing reality and attempts to limit the costs and complexity of groundwater management is to resort to stochastic models able to consistently encompass this intrinsic uncertainty (e.g., Bear 1972; Dagan 1989; Kapoor and Gelhar 1994a, b; Pannone and Kitanidis 1999; Rubin 2003). Classic porous media stochastic theories assume that the log-conductivity Y(x) = ln K(x) (where K is the local hydraulic conductivity and x indicates the vector position in a Cartesian reference frame), is a statistically stationary and normally distributed random function, completely characterized by a constant mean and fast-decaying correlation function. As well known, statistical stationarity implies the existence of a single representative scale of heterogeneity, i.e., the so-called integral scale. However, in the case of regional aquifers, solute transport is typically influenced by several scales of structural variability, which progressively come into play as the travel distance increases. In these conditions, the semivariogram of Y (i.e., the variance of log-K spatial increments) may mathematically be represented by a simple-scaling, power-law continuous model according to Neuman (1990): γ_Y(r) = ar^b, where a is a dimensional constant and one-half of the scaling exponent H = b/2 is known as the Hurst coefficient (e.g., Feder 1988).

The spreading of solutes in heterogeneous velocity fields is classically investigated in terms of the longitudinal macrodispersion coefficient, be it the so-called ergodic macrodispersion coefficient, given by the (time) rate of change of the single-particle position covariance, and incorporates the uncertainty related to plume centroid location, or the effective, nonergodic macrodispersion coefficient, which consists in the rate of change of the expected longitudinal central inertia moment and only accounts for scales of heterogeneity smaller than the actual plume dimension. The difference between them and the specificity of their use were widely discussed by Fischer et al. (1979) in the context of turbulent mixing, by Kitanidis (1988), Dagan (1990) and Rajaram and Gelhar (1993) in the context of transport in single-scale porous formations, and by Glimm et al. (1993) in the context of transport in evolving-scale porous formations. The quantitative implications of the ergodic and nonergodic assumption on transport in evolving-scale formations for purely advective regimes were then afforded by Dagan (1994), who obtained first-order analytical solutions for the longitudinal macrodispersion coefficient in formations described by power-law semivariograms. The main result of Dagan’s study was that the ergodic assumption, given for granted in the formulation of most stochastic models, and invariably leading to super-diffusive transport (single-particle covariance more than linearly increasing in time even at large times), cannot be assumed as universally valid. On the other hand, rather than showing an anomalous continuous growth, the physically more relevant expected central inertia moment would reach a Fickian asymptotic limit in evolving-scale formations characterized by b < 1, with the anomalous case that would only be recovered for b > 1.

As a consequence of the existence of multiple scales of medium heterogeneity, anomalous transport typically manifests itself in long-tail spatial and/or temporal distributions of solute concentration at given times and locations. It should be emphasized that the strong heterogeneity associated with evolving-scale antipersistently (b < 1) and, mainly, persistently (b > 1) correlated porous structures (e.g., Painter and Mahinthakumar 1999) determines the presence of very low-conductivity transition zones and preferential flow paths (Edery et al. 2014), as also demonstrated by the numerical generation of single-realization samples in Bellin et al. (1996). Therefore, it can easily be understood how the emergence of non-Fickian behavior is a highly probable event also in the case of transport through fractures (Frampton and Cvetkovic 2009; Kang et al. 2016; Wang and Cardenas 2017; Shahkarami et al. 2019; Hu et al. 2020; Liu et al. 2022), where high-velocity channels delimit almost impervious matrix portions. For such types of geologic formations, the evolving-scale persistent log-K correlation may therefore represent a suitable modelling tool. An important consequence of this sort of “pipeline effect” in markedly heterogeneous aquifers subject to chemical/biological contamination is represented by high concentrations with very long residence times. The implications of the scaling behavior in the context of the probabilistic risk analysis were first discussed, among other things, by Moslehi et al. (2016) and Moslehi and de Barros (2017). At the first-order in the velocity fluctuations, the study by Suciu et al. (2015) showed how anomalous super-diffusive behavior may result from the linear combination of independent random fields characterized by short-range correlation functions and increasing variance and integral scales (e.g., Gelhar 1986; Neuman 1990). According to the same type of linear decomposition of the log-conductivity fluctuations, a preliminary study by the author (Pannone 2017) proposed a new theoretical approach for the determination of the ergodic macrodispersion coefficient as a function of initial plume size and related Peclét number magnitude. The results of the investigation invariably predicted asymptotically Fickian macrodispersion for purely advective regimes (Pe→∞), and super-diffusive transport in the presence of non-negligible local dispersion. Note that, in the second case, starting from the asymptotic expression of the longitudinal macrodispersion coefficient for high Peclét numbers and short-range log-K correlations, the effect of diffusion was “a posteriori” superposed by a local dispersion-related, mobile-extreme integration over the hierarchy of scales. The essential role of local dispersion (sometimes simply named “diffusion”) in solute macrodispersion and dilution had already been explored in the context of subsurface flow and transport by Kapoor and Gelhar (1994a, b), Pannone and Kitanidis (1999) and, afterwards, in the context of river-flow and transport, by Pannone (2010, 2012, 2014). Overall, these studies showed that, in the case of heterogeneous flow fields characterized by short-range correlations, macro-dispersion and dilution are singularly driven by the interplay of advective heterogeneities and diffusive-like mechanisms. Accounting for this interplay, the two-dimensional (2D) Lagrangian study by Fiori (2001) derived the nonergodic macrodispersion coefficients for evolving-scale structures directly from the corresponding stationary-increments semivariogram. The direct use of the semivariogram determined the exigence of adopting invariably finite initial dimensions, at which the definition of Peclét was tied, of the consequently arbitrary-shape and already partially diluted solute body. Finite initial dimensions were indeed necessary in order to filter out the imposed wave-length lower cutoff. The present 3-D investigation explicitly resorts to the superposition principle—an “emanation” of the underlying linear formulation, which is ubiquitous in all the analytical approaches to the topic—for point-like solute pulses. The superposition principle allows for relating the macrodispersion coefficients to a “fractal” log-K covariance obtained by a (double) integration of periodic or pseudo-periodic correlation functions that belong to the double hierarchy of translationally-invariant log-K fields into which the original fractal field is decomposed, with no need for lower or upper cutoff. Furthermore, by reinterpreting antipersistent (b < 1) and persistent (b > 1) evolving-scale correlations as a double hierarchy of oscillatory nonperiodic and periodic fields, respectively, the present study envisions a possible alternative methodology for their numerical generation.

Materials and methods

The aim here is to begin by showing, via simple analytical identities, that power law log-conductivity semivariograms with exponent b ranging between 0 and 2 can indeed be obtained by the superposition of a double hierarchy of isotropic stationary fields (i.e., Y = log K stationary fields for which the covariance function R_Y and the semivariogram ${\gamma}_Y={\sigma}_Y^2-{R}_Y$ (with ${\sigma}_Y^2$ indicating the log-K variance) just depend on the norm of the distance vector r = |r|). Such a double hierarchy has to be intended as belonging to the space of wave-number λ (the inverse of the primary-hierarchy integral scales) and μ (the inverse of the secondary-hierarchy periods):

$${\displaystyle \begin{array}{l}{\gamma}_Y(r)=a{r}^b={\int}_0^{\infty }{\int}_0^{\infty}\frac{\partial^2{\gamma}_Y}{\partial \lambda \partial \mu } d\mu d\lambda =\\ {}{\int}_0^{\infty }{\int}_0^{\infty}\frac{\phi }{\lambda^{1+b}}\left\{\sqrt{\frac{1}{\lambda \pi}}\exp \left(-\frac{\mu^2}{4\lambda}\right)\left[1-\cos \left(\sqrt{r}\mu \right)\right]\right\} d\mu d\lambda\ 0 < b < 1\\ {}\kern13.75em a=-\phi \Gamma \left(-b\right)\end{array}}$$

(1)

and

$${\displaystyle \begin{array}{l}{\gamma}_Y(r)=a{r}^b={\int}_0^{\infty }{\int}_0^{\infty}\frac{\partial^2{\gamma}_Y}{\partial \lambda \partial \mu } d\mu d\lambda =\\ {}{\int}_0^{\infty }{\int}_0^{\infty}\frac{\phi }{\lambda^{1+b}}\left\{\frac{2}{\lambda \pi}\exp \left(-\frac{\mu^2}{\pi {\lambda}^2}\right)\left[1-\cos \left( r\mu \right)\right]\right\} d\mu d\lambda \kern4.75em 1\le b<2\\ {}\kern14.25em a=-\frac{\phi }{2}\Gamma \left(-\frac{b}{2}\right){\left(\frac{\sqrt{\pi }}{2}\right)}^b\end{array}}$$

(2)

where $r=\sqrt{x_1^2+{x}_2^2+{x}_3^2}$ indicates the norm of the distance, Γ the Gamma function (Gradshteyn and Ryzhik 1994), ϕ is a dimensional proportionality constant, and

$$\frac{\phi }{\lambda^{1+b}}=\frac{d{\sigma}_Y^2}{d\lambda}$$

(3)

identifies the variance density associated with the primary-hierarchy frequency domain. Specifically, the integration over μ leads to semivariograms given by finite variance and exponential (0 < b < 1) or Gaussian (1 ≤b < 2) covariance; the further integration over λ respectively leads to the power law for the two ranges of scaling exponent. That means that the evolving-scale log-conductivity fluctuation ${\overset{\sim }{Y}}^{\prime}\left(\textbf{x}\right)=Y\left(\textbf{x}\right)-\left\langle Y\right\rangle$ (with 〈Y〉 indicating the log-conductivity ensemble mean, that is, the log-conductivity mean computed over the whole ensemble of possible realizations) can formally be interpreted as a double series of mutually independent components:

$${\overset{\sim }{Y}}^{\prime}\left(\textbf{x}\right)=\sum\limits_{\lambda =0}^{\infty }\sum\limits_{\mu =0}^{\infty }{Y}^{\prime}\left(\textbf{x};\lambda, \mu \right)$$

(4)

In the hypothesis of reasonably small ϕ, even the larger heterogeneous subunits (μ→0, λ→0) will be characterized by a reasonably small variance:

$$d{\sigma}_Y^2=\frac{\partial^2{\sigma}_Y^2}{\partial \lambda \partial \mu } d\mu d\lambda =\frac{\phi }{\lambda^{3/2+b}{\pi}^{1/2}}\ \exp \left(-\frac{\mu^2}{4\lambda}\right) d\mu d\lambda \kern3.75em 0<b<1$$

(5)

and

$$d{\sigma}_Y^2=\frac{\partial^2{\sigma}_Y^2}{\partial \lambda \partial \mu } d\mu d\lambda =\frac{2\phi }{\lambda^{2+b}\pi }\ \exp \left(-\frac{\mu^2}{\pi {\lambda}^2}\right) d\mu d\lambda \kern4em 1\le b<2$$

(6)

respectively, allowing for flow and transport first-order treatment (e.g., Dagan 1989; Rubin, 2003). From Darcy’s law, v =  –K(x) ∇ h(x)/n, where v indicates fluid velocity vector, h hydraulic head and n medium porosity, and continuity for incompressible fluids in incompressible solid matrices ∇ ∙ v = 0, the following steady flow equation is derived (e.g., Bear 1972; Dagan 1989):

$${\displaystyle \begin{array}{c}{\nabla}^2h\left(\textbf{x}\right)+\boldsymbol{\nabla}h\left(\textbf{x}\right)\cdot \boldsymbol{\nabla}Y\left(\textbf{x}\right)=\\ {}={\nabla}^2{h}^{\prime}\left(\textbf{x}\right)-\textbf{J}\cdot \boldsymbol{\nabla}{Y}^{\prime}\left(\textbf{x}\right)+\boldsymbol{\nabla}{h}^{\prime}\left(\textbf{x}\right)\cdot \boldsymbol{\nabla}{Y}^{\prime}\left(\textbf{x}\right)=0\end{array}}$$

(7)

In Eq. (7), <h(x)> is the hydraulic head ensemble mean, h′ the corresponding fluctuation, and

$$\textbf{J}=-\boldsymbol{\nabla}<h\left(\textbf{x}\right)>$$

(8)

is the constant head mean gradient vector. In the Fourier domain (wave-number vector k), Eq. (7) becomes:

$${\left|\textbf{k}\right|}^2\hat{h}\left(\textbf{k}\right)+{\int}_{\bar {\textbf{k}}}\bar {\textbf{k}}\cdot \left(\textbf{k}-\bar {\textbf{k}}\right)\hat{Y}\left(\textbf{k}-\bar {\textbf{k}}\right)\hat{h}\left(\bar {\textbf{k}}\right)d\bar {\textbf{k}}=j\textbf{J}\cdot \textbf{k}\hat{Y}\left(\textbf{k}\right)$$

(9)

where $j=\sqrt{-1}$, $\overline{\textbf{k}}$ is the integration variable, and the circumflex accent indicates Fourier transforms:

$$\hat{h}\left(\textbf{k}\right)={\int}_{\textbf{x}}{h}^{\prime}\left(\textbf{x}\right)\exp \left(-j2\pi \textbf{k}\cdot \textbf{x}\right)d\textbf{x}$$

(10)

$$\hat{Y}\left(\textbf{k}\right)={\int}_{\textbf{x}}{Y}^{\prime}\left(\textbf{x}\right)\exp \left(-j2\pi \textbf{k}\cdot \textbf{x}\right)d\textbf{x}$$

(11)

Provided that flow and transport linear theory applies, the convolution term in Eq. (9) can be neglected (it is a function of the product of two small deviatoric terms) leading to the Fourier transform of a Poisson equation:

$${\left|\textbf{k}\right|}^2\hat{h}\left(\textbf{k}\right)=j\textbf{J}\cdot \textbf{k}\hat{Y}\left(\textbf{k}\right)$$

(12)

From Eq. (12), and from first-order Darcy’s law expressed in terms of fluctuating (or deviatoric) quantities:

$${\textbf{v}}^{\prime}\left(\textbf{x}\right)\cong \frac{\exp \left\langle Y\right\rangle }{n}\left[{Y}^{\prime}\left(\textbf{x}\right)\textbf{J}-\boldsymbol{\nabla}{h}^{\prime}\left(\textbf{x}\right)\right]$$

(13)

one obtains in the Fourier domain (Dagan 1989):

$${\hat{v}}_i\left(\textbf{k}\right)=\hat{Y}\left(\textbf{k}\right)\left({U}_i-\sum\limits_{l=1}^3\frac{k_i{k}_l}{{\left|\textbf{k}\right|}^2}{U}_l\right)$$

(14)

where ${\hat{v}}_i$ is the Fourier transform of the ith component of velocity fluctuation and U = (U₁, U₂, U₃) = exp < Y > J/n the corresponding ensemble mean. What Eqs. (12) and (14) say is that, at the first-order, each Fourier component of the log-conductivity field corresponds to a single component of hydraulic head and to a single component of velocity. Therefore, the same properties of linear superposition holding for Y (see Eq. 4) apply to h and v as well. It should in any case be noted that the superposition principle could be applied irrespective of the linearization of the flow and transport problem: as a matter of fact, taking the ensemble mean of Eq. (9), with the ensemble mean of the deviatoric $\hat{Y}$ and $\hat{h}$ equal to zero by definition, one obtains, for $\textbf{k}\ne \bar {\textbf{k}}$, $\left\langle \hat{Y}\left(\textbf{k}-\bar {\textbf{k}}\right)\hat{h}\left(\bar {\textbf{k}}\right)\right\rangle =0$, meaning that the Fourier amplitude of log-conductivity and hydraulic head are statistically uncorrelated at different scales.

Based on the classic Lagrangian theory (e.g., Dagan 1989), the trajectory of the generic solute particle belonging to a plume that originates from an instantaneous point injection at x = 0 can be expressed as:

$${\displaystyle \begin{array}{c}\textbf{X}(t)={\int}_0^t\textbf{v}\left[\textbf{X}\left({t}_1\right)\right]d{t}_1+{\textbf{X}}_{\mathrm{B}}(t)={\int}_0^t\left\langle \textbf{v}\left[\textbf{X}\left({t}_1\right)\right]\right\rangle d{t}_1+{\int}_0^t{\textbf{v}}^{\prime}\left[\textbf{X}\left({t}_1\right)\right]d{t}_1+{\textbf{X}}_{\mathrm{B}}(t)\\ {}=\textbf{U}t+{\textbf{X}}^{\prime }(t)+{\textbf{X}}_{\mathrm{B}}(t)\end{array}}$$

(15)

where t is the time and X_B the zero-mean Brownian component. Additionally,

$$\left\langle \textbf{X}\left(\mathrm{t}\right)\right\rangle =\textbf{U}t={\int}_0^t\left\langle \textbf{v}\left[\textbf{X}\left({t}_1\right)\right]\right\rangle d{t}_1$$

(16)

indicates the ensemble mean particle position vector,

$${\textbf{X}}^{\prime }(t)={\int}_0^t{\textbf{v}}^{\prime}\left[\textbf{X}\left({t}_1\right)\right]d{t}_1={\int}_0^t{\textbf{u}}^{\prime}\left({t}_1\right)d{t}_1$$

(17)

its deviatoric (fluctuating) component, and u′ the Lagrangian velocity fluctuation. Notice that, based on subsurface flow and transport Lagrangian linear theory for finite Peclét numbers (e.g., Pannone and Kitanidis 1999), solute particles sample the velocity distribution along the ensemble mean trajectory (Ut) only perturbed by the local-dispersive contribution represented by X_B, leading to Eq. (18), which is the first-order, linearized version of Eq. (17):

$${\textbf{X}}^{\prime }(t)\cong {\int}_0^t{\textbf{v}}^{\prime}\left[\textbf{U}{t}_1+{\textbf{X}}_{\mathrm{B}}\left({t}_1\right)\right]d{t}_1$$

(18)

When applied to multiple-scale (or hierarchical) heterogeneity according to Eqs. (18), (4) and (14), the superposition principle allows for extending to the hierarchical total trajectory fluctuation the double linear decomposition rule:

$${\overset{\sim }{\textbf{X}}}^{\prime }(t)={\sum}_{\lambda =0}^{\infty }{\sum}_{\mu =0}^{\infty }{\textbf{X}}^{\prime}\left(t;\lambda, \mu \right)$$

(19)

The next step is to proceed with the derivation of one- (X_ii) and two- (Θ_ii) particle moments rate of change from Eqs. (18) and (19), based on the formulation by Pannone and Kitanidis (1999) for single-scale heterogeneity, and in the plume principal reference frame (no off-diagonal components of one- and two-particle moment tensors):

$$\frac{d{X}_{ii}}{dt}=\frac{d\left\langle {\overset{\sim }{X}}_i^{\prime 2}\right\rangle }{dt}=\left\langle 2{\overset{\sim }{X}}_i^{\prime}\frac{d{\overset{\sim }{X}}_i^{\prime }}{dt}\right\rangle =2{\int}_0^t\left\langle {\tilde{v}}_i^{\prime}\left[\textbf{U}t+{\textbf{X}}_{\mathrm{B}}(t)\right]{\tilde{v}}_i^{\prime}\left[\textbf{U}{t}_1+{\textbf{X}}_{\mathrm{B}}\left({t}_1\right)\right]\right\rangle d{t}_1+2D$$

(20)

$${\displaystyle \begin{array}{c}\frac{d{\Theta}_{ii}}{dt}=\frac{d\left\langle {\overset{\sim }{X}}_i^{\prime }{\overset{\sim }{Z}}_i^{\prime}\right\rangle }{dt}=\left\langle {\overset{\sim }{X}}_i^{\prime}\frac{d{\overset{\sim }{Z}}_i^{\prime }}{dt}\right\rangle +\left\langle {\overset{\sim }{Z}}_i^{\prime}\frac{d{\overset{\sim }{X}}_i^{\prime }}{dt}\right\rangle =\\ {}{\int}_0^t\left\langle {\tilde{v}}_i^{\prime}\left[\textbf{U}t+{\textbf{Z}}_{\mathrm{B}}(t)\right]{\tilde{v}}_i^{\prime}\left[\textbf{U}{t}_1+{\textbf{X}}_{\mathrm{B}}\left({t}_1\right)\right]\right\rangle d{t}_1\\ {}\kern15.5em +{\int}_0^t\left\langle {\tilde{v}}_i^{\prime}\left[\textbf{U}t+{\textbf{X}}_{\mathrm{B}}(t)\right]{\tilde{v}}_i^{\prime}\left[\textbf{U}{t}_1+{\textbf{Z}}_{\mathrm{B}}\left({t}_1\right)\right]\right\rangle d{t}_1\end{array}}$$

(21)

where $D=\left\langle {X}_{\mathrm{B}i}^2\right\rangle /2t$ indicates the isotropic local dispersion coefficient and $\overset{\sim }{\boldsymbol{Z}}(t)$ the second particle trajectory. Note that the local dispersion coefficient D appears in the first Fick’s law, expressing the specific mass flux vector q from higher to lower concentration zones in any isotropic and purely diffusive transport process: q =  –D $\nabla$ c (e.g., Fischer et al. 1979). The reason why it is accounted for in Eq. (20) and not in Eq. (21) is that the purely diffusive displacements of two distinct solute particles are uncorrelated by definition. The ensemble means indicated by angle brackets in Eqs. (20) and (21) can be obtained by a conditional ensemble averaging operation over the whole space of possible particle positions for fixed Brownian displacements, followed by a further ensemble averaging operation based on the probability density function of these displacements. Due to the stationarity of each component of the double-hierarchy, and invoking the superposition principle, with the hierarchical velocity covariance function that just depends on the generic distance vector x-z, one obtains (adapted from the single-scale expression by Pannone and Kitanidis 1999):

$$\frac{d{X}_{ii}}{dt}=2{\int}_0^t{\int}_{\textbf{x}-\textbf{z}}{\overset{\sim }{R}}_{vii}\left(\textbf{x}-\textbf{z}|{\boldsymbol{X}}_{\mathrm{B}}\right){f}_{\boldsymbol{X}\mathrm{B}}\left[\left(\textbf{x}-\textbf{z}\right)-\textbf{U}\left(t-{t}_1\right);t-{t}_1,D\right]d\left(\textbf{x}-\textbf{z}\right)d\left(t-{t}_1\right)+2D$$

(22)

and

$$\frac{d{\Theta}_{ii}}{dt}=2{\int}_0^t{\int}_{\textbf{x}}{\int}_{\textbf{z}}{\overset{\sim }{R}}_{vii}\left(\textbf{x}-\textbf{z}|{\boldsymbol{X}}_{\mathrm{B}},{\boldsymbol{Z}}_{\mathrm{B}}\right){f}_{\textbf{X}\mathrm{B}\textbf{Z}\mathrm{B}}\left(\textbf{x}-\textbf{U}t,\textbf{z}-\textbf{U}{t}_1;t,{t}_1,D\right)d\textbf{x}d\textbf{z}d{t}_1$$

(23)

where f_XB indicates the Gaussian probability density function of the single Brownian displacement:

$${f}_{\textbf{X}\mathrm{B}}\left(\textbf{r}-\boldsymbol{\upeta}; s,D\right)=\prod\limits_{i=1}^3\frac{1}{\sqrt{4\pi Ds}}\exp \left[-\frac{{\left({r}_i-{\eta}_i\right)}^2}{4 Ds}\right]$$

(24)

f_XBZB is the joint probability density function of the two independent Brownian displacements:

$${f}_{\textbf{X}\mathrm{B}\textbf{Z}\mathrm{B}}\left(\textbf{x}-\boldsymbol{\upeta}, \textbf{z}-\boldsymbol{\upupsilon}; s,w,D\right)={f}_{\textbf{X}\mathrm{B}}\left(\textbf{x}-\boldsymbol{\upeta}; s,D\right){f}_{\textbf{Y}\mathrm{B}}\left(\textbf{z}-\boldsymbol{\upupsilon}; w,D\right)$$

(25)

${\overset{\sim }{R}}_{vii}\left(\textbf{x}-\textbf{z}\right)$ is the hierarchical covariance of the ith velocity component:

$${\overset{\sim }{R}}_{vii}\left(\textbf{x}-\textbf{z}\right)=\left\langle {\tilde{v}}_i^{\prime}\left(\textbf{x}\right){\tilde{v}}_i^{\prime}\left(\textbf{z}\right)\right\rangle ={\int}_{\boldsymbol{k}}{\overset{\sim }{S}}_{vii}\left(\textbf{k}\right)\exp \left[j2\pi \textbf{k}\cdot \left(\textbf{x}-\textbf{z}\right)\right]d\textbf{k}$$

(26)

${\overset{\sim }{S}}_{vii}\left(\textbf{k}\right)$ is the corresponding Fourier transform and the vertical slash indicates the statistical conditioning. The substitution of Eqs. (24), (25) and (26) into Eqs. (22) and (23), and the subsequent space-time integrations yield, for U = (U,0,0):

$$\frac{d{X}_{ii}}{dt}=2{\int}_{\textbf{k}}\frac{{\overset{\sim }{S}}_{vii}\left(\textbf{k}\right){\Psi}_{\textbf{X}}\left(\textbf{k},t\right)}{{\left(2\pi U{k}_1\right)}^2+16{\pi}^4{D}^2{k}^4}d\textbf{k}+2D$$

(27)

and

$$\frac{d{\Theta}_{ii}}{dt}=2{\int}_{\textbf{k}}\frac{{\overset{\sim }{S}}_{vii}\left(\textbf{k}\right){\Psi}_{\boldsymbol{\Theta}}\left(\textbf{k},t\right)}{{\left(2\pi U{k}_1\right)}^2+16{\pi}^4{D}^2{k}^4}d\textbf{k}$$

(28)

with:

$${\displaystyle \begin{array}{c}{\Psi}_{\textbf{X}}\left(\textbf{k},t\right)=4{\pi}^2D{k}^2+\kern19.25em \\ {}\exp \left[-4{\pi}^2D{k}^2t\right]\left[2\pi U{k}_1\sin \left(2\pi U{k}_1t\right)-4{\pi}^2D{k}^2\cos \left(2\pi U{k}_1t\right)\right]\end{array}}$$

(29)

$${\displaystyle \begin{array}{c}{\Psi}_{\boldsymbol{\Theta}}\left(\textbf{k},t\right)=-4{\pi}^2D{k}^2\exp \left[-8{\pi}^2D{k}^2t\right]+\exp \left[-4{\pi}^2D{k}^2t\right]\\ {}\left[2\pi U{k}_1\sin \left(2\pi U{k}_1t\right)+4{\pi}^2D{k}^2\cos \left(2\pi U{k}_1t\right)\right]\end{array}}$$

(30)

and $k=\sqrt{k_1^2+{k}_2^2+{k}_3^2}$.

The conceptual interpretation of quantities Ψ_X(k, t) and Ψ_Θ(k, t) is that of two different frequency-time kernels that characterize and crucially distinguish the relationship between (1) the hierarchical Fourier transform of the covariance of the velocity components and (2) the one- and two-particle moment rates of change. The hierarchical Fourier transform of covariance of the velocity components, in turn, is obtained by superposition from the classical linear theory (e.g., Dagan 1989):

$${\overset{\sim }{S}}_{vii}\left(\textbf{k}\right)={U}^2{\left({\delta}_{i1}-\frac{k_1{k}_i}{\left|\textbf{k}\right|}\right)}^2{\overset{\sim }{S}}_Y\left(\textbf{k}\right)$$

(31)

where ${\overset{\sim }{S}}_Y\left(\textbf{k}\right)$ indicates the hierarchical Fourier transform of the log-conductivity covariance and δ_ij is the Kronecker delta function. It must be stressed that, in its original version, Eq. (31) holds for the single stationary component of the hierarchy. The superposition principle, and the statistical independence of all the hierarchical components, allow for applying it to the evolving-scale field by a simple μ-λ double integration. Following this procedure, the hierarchical Fourier transform of the log-conductivity covariance in Eq. (31)) is obtained as:

$${\overset{\sim }{S}}_{Y\mathrm{E}}\left(\textbf{k}\right)=4\pi \phi \mathrm{B}\left(\frac{3+b}{2},\frac{1-b}{2}\right){\left[4{\pi}^2\left({k}_1^2+{k}_2^2+{k}_3^2\right)\right]}^{-\frac{3+b}{2}}\kern2.75em 0<b<1$$

(32)

and

$${\overset{\sim }{S}}_{Y\mathrm{G}}\left(\textbf{k}\right)=4\phi \Gamma \left(\frac{3+b}{2}\right){\left[4\pi \left({k}_1^2+{k}_2^2+{k}_3^2\right)\right]}^{-\frac{3+b}{2}}\kern3.5em 1\le b<2$$

(33)

with B here indicating Beta function (Gradshteyn and Ryzhik 1994), and subscripts E and G respectively referring to the exponential and Gaussian primary-hierarchy sequence. Finally, as demonstrated by Pannone and Kitanidis (2004) for finite Peclét and point solute pulses in single-scale heterogeneous formations, one- and two-particle covariance combine to give the expected central inertia moment:

$$\left\langle {I}_{ii}(t)\right\rangle ={X}_{ii}(t)-{\Theta}_{ii}(t)$$

(34)

and, therefore, the nonergodic macrodispersion coefficient:

$${D}_{mii}(t)=\frac{1}{2}\frac{d\left\langle {I}_{ii}(t)\right\rangle }{dt}=\frac{1}{2}\left(\frac{d{X}_{ii}}{dt}-\frac{d{\Theta}_{ii}}{dt}\right)$$

(35)

which is a measure of the expected dispersion about the expected center of mass, that is, an average measure of the aquifer volume affected by the presence of the solute plume, regardless of the effective point concentration values. See the Appendix for a mathematical discussion about the statistical representativeness of the ensemble mean central inertia moments in the context of evolving-scale log-K distributions. See also Dentz and Barros (2013) for a discussion of the time-behavior of nonergodic dispersion variance in short-range log-K distributions.

Analytical derivation of the non-ergodic macrodispersion coefficients

Equation (35) will be here applied to evolving-scale porous media, accounting for the superposition principle, at large times and for large, though finite, Peclét numbers:

$$\mathrm{Pe}=\frac{U}{\phi^{1/b}D}>>1$$

(36)

It must be stressed that 1/ϕ^1/b, related to the semivariogram scaling coefficient a after Eqs. (1) and (2), is the only possible spatial scale that can be defined in such a context. The use of a Peclét number related to solute body initial dimension as in previous similar works would not be possible in the present study, since it deals with highly concentrated, point-like injections. On the other hand, the opinion of the author is that such a choice would not in any case make sense, provided that Pe should by definition represent the ratio of the involved spatial range of advective to diffusive transport mechanisms. The identification of the initial (transverse) dimension with a well-defined and fixed advective correlation range is only consistent with a first-order formulation for purely advective regimes, that is, in the presence of negligible transverse displacements (no chance to sample transverse heterogeneity scales larger than it). In a spherical coordinate system, one obtains from Eqs. (27) to (31), and for the quantitatively dominant advection-related parts of the macrodispersion coefficients:

$${D}_{mii}(t)=\frac{\phi C{U}^2}{16{\pi}^4{D}^2}{\int}_0^{\infty }{\int}_0^{\pi }{\int}_0^{2\pi}\frac{F_{ii}\left(\theta, \varphi \right){\Psi}_{\mathrm{X}{\Theta }}\left(k,\theta, \varphi, t\right)\sin \theta }{\left({k}^2+\frac{U^2}{4{\pi}^{2}D^{2}}{\sin}^2\theta {\cos}^2\varphi \right){k}^{3+b}} d\varphi d\theta dk$$

(37)

where:

$${\displaystyle \begin{array}{l}{F}_{11}\left(\theta, \varphi \right)={\left(1-{\sin}^2\theta {\cos}^2\varphi \right)}^2\\ {}{F}_{22}\left(\theta, \varphi \right)={\sin}^4\theta {\cos}^2\varphi {\sin}^2\varphi \\ {}{F}_{33}\left(\theta, \varphi \right)={\sin}^2\theta {\cos}^2\theta {\cos}^2\varphi \end{array}}$$

(38)

$$\displaystyle \begin{aligned}{\Psi}_{\mathrm{X}{\Theta}}\left(k,\theta, \varphi, t\right)&={\Psi}_{\mathrm{X}}\left(k,\theta, \varphi, t\right)-{\Psi}_{\Theta}\left(k,\theta, \varphi, t\right)\\ & = 4{\pi}^{2}D{k}^{2}\left[1+\exp \left(-8{\pi}^2D{k}^2t\right)-\right. \\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\left.2\exp \left(-4{\pi}^{2} D{k}^2t\right)\cos \left(2\pi k\!\sin \theta \cos\!\varphi Ut\right)\right]\end{aligned}$$

(39)

$$C=\left\{\begin{array}{c}\frac{\mathrm{B}\left(\frac{3+b}{2},\frac{1-b}{2}\right)}{4^{\left(1+b\right)/2}{\pi}^{2+b}}\kern0.5em 0<b<1\\ {}\frac{\Gamma \left(\frac{3+b}{2}\right)}{4^{\left(1+b\right)/2}{\pi}^{\left(3+b\right)/2}}\kern0.75em 1\le b<2\end{array}\right.$$

(40)

and

$$\left\{\begin{array}{l}{k}_1=k\sin \theta \cos \varphi \\ {}{k}_2=k\sin \theta \sin \varphi \\ {}{k}_3=k\cos \theta \\ {}d{k}_1d{k}_2d{k}_3={k}^2\sin \theta d\varphi d\theta dk\end{array}\right.$$

(41)

Switching to dimensionless coordinates (ϕ^–1/b as a length-scale and ϕ^–1/b/U as a time-scale), and considering that the integration over φ from 0 to 2π can be transformed into 4 times the integration over φ from 0 to π/2 due to the properties of sin²φ, cos²φ, and the even cosφ-function cos(2πk sin θ cos φUt), Eq. (34) becomes:

$${\displaystyle \begin{array}{l}{D}_{mii}^{\ast}\left(\tau \right)=\frac{C\mathrm{Pe}}{\pi^2}.\\ {}{\int}_0^{\infty }{\int}_0^{\pi }{\int}_0^{\pi /2}{F}_{ii}\left(\theta, \varphi \right)\sin \theta\cdot \\ {}\frac{\left\{{k}^{\ast 2}\left[1+\exp \left(-\frac{8{\pi}^2{k}^{\ast 2}}{\mathrm{P}\mathrm{e}}\tau \right)-2\exp \left(-\frac{4{\pi}^2{k}^{\ast 2}}{\mathrm{P}\mathrm{e}}\tau \right)\cos \left(2\pi {k}^{\ast}\sin \theta \cos \varphi \tau \right)\right]\right\}}{\left({k^{\ast}}^2+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi \right){k}^{\ast 3+b}} d\varphi d\theta d{k}^{\ast}\end{array}}$$

(42)

with:

$${D}_{mii}^{\ast }=\frac{D_{mii}}{U}{\phi}^{1/b}$$

(43)

$${k}^{\ast }=k{\phi}^{1/b}$$

(44)

and

$$\tau = tU{\phi}^{1/b}$$

(45)

The integration over the dimensionless wave-number k^* can be split as follows:

$${\int}_0^{\infty }f\left({k}^{\ast}\right)g\left({k}^{\ast}\right)d{k}^{\ast }={\int}_0^1f\left({k}^{\ast}\right)g\left({k}^{\ast}\right)d{k}^{\ast }+{\int}_1^{\infty }f\left({k}^{\ast}\right)g\left({k}^{\ast}\right)d{k}^{\ast }$$

(46)

with:

$$f\left({k}^{\ast}\right)=\left\{{k}^{\ast 2}\left[1+\exp \left(-\frac{8{\pi}^2{k}^{\ast 2}}{\mathrm{Pe}}\tau \right)-2\exp \left(-\frac{4{\pi}^2{k}^{\ast 2}}{\mathrm{Pe}}\tau \right)\cos \left(2\pi {k}^{\ast}\sin \theta \cos \varphi \tau \right)\right]\right\}$$

(47)

and:

$$g\left({k}^{\ast}\right)=\frac{1}{\left({k^{\ast}}^2+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi \right){k}^{\ast 3+b}}$$

(48)

In the case of the large-scale contribution (0 < k^* < 1) at large Peclét, provided that 0 < sin θ < 1, 0 < cos φ < 1, the Taylor expansion of f around k^* = 0 yields:

$$f\left({k}^{\ast}\right)\to 4{\pi}^4{k}^{\ast 4}{\sin}^2\theta {\cos}^2\varphi {\tau}^2$$

(49)

Thus, for 4π²k^∗2/Pe² << 1:

$${D}_{mii\mathrm{A}}^{\ast}\left(\tau \right)=\frac{16{\pi}^2C}{\left(2-b\right)\mathrm{Pe}}{\int}_0^{\pi }{\int}_0^{\pi /2}{F}_{ii}\left(\theta, \varphi \right)\sin \theta d\varphi d\theta$$

(50)

and, respectively, for 0 < b < 1 and 1 ≤ b < 2:

$${\displaystyle \begin{array}{l}{D}_{m11\mathrm{A}}^{\ast}\left(\tau \right)={D}_{m\mathrm{LA}}^{\ast}\left(\tau \right)=\frac{224}{15\mathrm{Pe}}\frac{\pi^{1-b}}{\left(2-b\right)}\frac{\mathrm{B}\left(\frac{3+b}{2},\frac{1-b}{2}\right)}{2^b}{\tau}^2\\ {}{D}_{m22\mathrm{A}}^{\ast}\left(\tau \right)={D}_{m33\mathrm{A}}^{\ast}\left(\tau \right)={D}_{mT\mathrm{A}}^{\ast}\left(\tau \right)=\frac{8}{15\mathrm{Pe}}\frac{\pi^{1-b}}{\left(2-b\right)}\frac{\mathrm{B}\left(\frac{3+b}{2},\frac{1-b}{2}\right)}{2^b}{\tau}^2\end{array}}$$

(51)

$${\displaystyle \begin{array}{l}{D}_{m11\mathrm{A}}^{\ast}\left(\tau \right)={D}_{m\mathrm{LA}}^{\ast}\left(\tau \right)=\frac{224}{15\mathrm{Pe}}\frac{\pi^{3/2}}{\left(2-b\right)}\frac{\Gamma \left(\frac{3+b}{2}\right)}{2^b{\pi}^{b/2}}{\tau}^2\\ {}{D}_{m22\mathrm{A}}^{\ast}\left(\tau \right)={D}_{m33\mathrm{A}}^{\ast}\left(\tau \right)={D}_{m\mathrm{TA}}^{\ast}\left(\tau \right)=\frac{8}{15\mathrm{Pe}}\frac{\pi^{3/2}}{\left(2-b\right)}\frac{\Gamma \left(\frac{3+b}{2}\right)}{2^b{\pi}^{b/2}}{\tau}^2\end{array}}$$

(52)

where subscript A indicates the contribution of dimensionless wave-number between 0 and 1. In the case of the smaller-scale contribution (1 < k^* < ∞, subscript B) at τ >> 1, f(k^*)→k^*2 and:

$${D}_{mii\mathrm{B}}^{\ast}\left(\tau \right)=\frac{C\mathrm{Pe}}{\pi^2}{\int}_1^{\infty }{\int}_0^{\pi }{\int}_0^{\pi /2}\frac{F_{ii}\left(\theta, \varphi \right)\sin \theta {k}^{\ast 2}}{\left({k^{\ast}}^2+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi \right){k}^{\ast 3+b}} d\varphi d\theta d{k}^{\ast }$$

(53)

The application of the second mean value theorem for infinite integrals (e.g., Gradshteyn and Ryzhik 1994) allows for writing:

$${\displaystyle \begin{array}{l}{\int}_1^{\infty}\frac{k^{\ast -\left(1+b\right)}}{k^{\ast 2}+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi }d{k}^{\ast }=\frac{1}{1+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi }{\int}_1^{\xi }{k}^{\ast -\left(1+b\right)}d{k}^{\ast }+\\ {}\frac{1}{\left(\infty +\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi \right)}{\int}_{\xi}^{\infty }{k}^{\ast -\left(1+b\right)}d{k}^{\ast}\\ {}\kern8em =\frac{1}{\left(1+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi \right)}{\int}_1^{\xi }{k}^{\ast -\left(1+b\right)}d{k}^{\ast}\end{array}}$$

(54)

with 1 < ξ <∞. Since, for any finite ξ:

$${\int}_1^{\xi }{k}^{\ast -\left(1+b\right)}d{k}^{\ast}\cong {\int}_1^{\infty }{k}^{\ast -\left(1+b\right)}d{k}^{\ast }=\frac{1}{2+b}$$

(55)

one obtains:

$${D}_{mii\mathrm{B}}^{\ast}\left(\tau \right)=\frac{C\mathrm{Pe}}{\pi^2\left(2+b\right)}{\int}_0^{\pi }{\int}_0^{\pi /2}\frac{F_{ii}\left(\theta, \varphi \right)\sin \theta }{\left(1+\frac{\mathrm{P}{\mathrm{e}}^2}{4{\pi}^2}{\sin}^2\theta {\cos}^2\varphi \right)} d\varphi d\theta$$

(56)

Finally, for large Peclét, and, respectively for 0 < b < 1 and 1 ≤ b < 2:

$${\displaystyle \begin{array}{l}{D}_{m11\mathrm{B}}^{\ast}\left(\tau \right)={D}_{m\mathrm{LB}}^{\ast}\left(\tau \right)=\frac{\mathrm{B}\left(\frac{3+b}{2},\frac{1-b}{2}\right)}{\left(2+b\right){4}^{\left(1+b\right)/2}{\pi}^{1+b}}\left(\frac{1}{2}-\frac{20}{3\mathrm{Pe}}\right)\\ {}{D}_{m22\mathrm{B}}^{\ast}\left(\tau \right)={D}_{m33\mathrm{B}}^{\ast}\left(\tau \right)={D}_{m\mathrm{TB}}^{\ast}\left(\tau \right)=\frac{\mathrm{B}\left(\frac{3+b}{2},\frac{1-b}{2}\right)}{3\left(2+b\right){4}^{\left(1+b\right)/2}{\pi}^{1+b}\mathrm{Pe}}\end{array}}$$

(57)

$${\displaystyle \begin{array}{l}{D}_{m11\mathrm{B}}^{\ast}\left(\tau \right)={D}_{m\mathrm{LB}}^{\ast}\left(\tau \right)=\frac{\Gamma \left(\frac{3+b}{2}\right)}{\left(2+b\right){4}^{\left(1+b\right)/2}{\pi}^{\left(1+b\right)/2}}\left(\frac{1}{2}-\frac{20}{3\mathrm{Pe}}\right)\\ {}{D}_{m22\mathrm{B}}^{\ast}\left(\tau \right)={D}_{m33\mathrm{B}}^{\ast}\left(\tau \right)={D}_{m\mathrm{TB}}^{\ast}\left(\tau \right)=\frac{\Gamma \left(\frac{3+b}{2}\right)}{3\left(2+b\right){4}^{\left(1+b\right)/2}{\pi}^{\left(1+b\right)/2}\mathrm{Pe}}\end{array}}$$

(58)

Results

Figures 1 and 2 show ${D}_{m\mathrm{L},\mathrm{T}}^{\ast}\left(\tau \right)$ for b = 0.5, 1.0 and 1.5 at Pe = 10⁴, with subscript L indicating x₁ (longitudinal) direction and subscript T indicating x₂ and x₃ (transverse) directions. As one can see, due to the invariably quadratic dependence on time of the large-scale components, the shorter-scale (Fickian) contributions turn out to be comparatively negligible, except for t→0. For the tested values of the scaling exponent b, the non-Fickian macrodispersion components appreciably decrease with it, or, equivalently, appreciably increase with the fractal dimension d = 2 – b/2. Additionally, the non-Fickian macrodispersion components decrease with Pe and, in the case of the longitudinal components, much faster than the Fickian counterparts (see Figs. 3 and 4, showing ${D}_{m\mathrm{LA},\mathrm{TA}}^{\ast}\left(\tau \right)$ and ${D}_{m\mathrm{LB},\mathrm{TB}}^{\ast}\left(\tau \right)$ at Pe = 10⁷). It should in any case be emphasized that, even for this considerably high value of Pe, the temporal transient in which Fickian and non-Fickian longitudinal components are quantitatively comparable is rather short. As auxiliary interpretation tools, Figs. 5, 6, 7 and 8 show, with ${D}_{mii}^{\ast}\left(\tau \right)={D}_{mii\mathrm{A}}^{\ast}\left(\tau \right)+{D}_{mii\mathrm{B}}^{\ast}\left(\tau \right)={\alpha}_{ii}{\tau}^2+{\beta}_{ii}$, the dependence of non-Fickian (α_ii) and Fickian (β_ii) coefficients on b at Pe = 10⁴. For both α_ii and β_ii, a minimum is detected in both antipersistent and persistent range. The reason for this behavior is that the smaller scales dominate at relatively smaller b (relatively larger fractal dimensions), while the larger scales dominate at relatively larger b (relatively smaller fractal dimensions). Indeed, the fractal dimension may synthetically be defined as a statistical measure of geometrical complexity: the larger the fractal dimension, the smaller the minimum scale at which the (self-similar) geometrical variability is replied. In the central part of the two ranges, Fickian and non-Fickian coefficients reach the respective minimum because of the more balanced contribution from all the involved scales. As a consequence of the wider range of scales of heterogeneity affecting it (1 < k^* < ∞), this minimum is more pronounced for the Fickian component, particularly in the antipersistent range. As a matter of fact, Figs. 9 and 10, which display an example of covariance functions densities at λ = 1:

$$\frac{\partial^2{\overset{\sim }{R}}_Y}{\partial \lambda \partial \mu }=\frac{\partial^2{\gamma}_Y}{\partial \lambda \partial \mu }-\frac{\partial^2{\overset{\sim }{\sigma}}_Y^2}{\partial \lambda \partial \mu }=\frac{\phi }{\lambda^{3/2+b}{\pi}^{1/2}}\ \exp \left(-\frac{\mu^2}{4\lambda}\right)\cos \left(\sqrt{r}\mu \right)\kern1.5em 0<b<1$$

(59)

and

$$\frac{\partial^2{\overset{\sim }{R}}_Y}{\partial \lambda \partial \mu }=\frac{\partial^2{\gamma}_Y}{\partial \lambda \partial \mu }-\frac{\partial^2{\overset{\sim }{\sigma}}_Y^2}{\partial \lambda \partial \mu }=\frac{2\phi }{\lambda^{2+b}\pi }\ \exp \left(-\frac{\mu^2}{\pi {\lambda}^2}\right)\cos \left( r\mu \right)\kern2em 1\le b<2,$$

(60)

show that the antipersistent range is characterized by additional non-negligible contributions from μ-scales appreciably different from zero. Specifically, from these figures it is clearly seen that, whereas the antipersistent correlation structures (Fig. 9) are characterized by a more distributed and slow variability (oscillations) associated with the presence of the pseudo-periodic $\cos \left(\sqrt{r}\mu \right)$, for the persistent counterparts, and λ being the same, the variability is essentially concentrated near μ = 0 (Fig. 10). Note that this behavior respectively translates into larger and smaller fractal dimensions (or greater and lower smaller-scale complexity), according to the preceding explanation.

Fig. 1

Fickian and non-Fickian longitudinal dispersion coefficients at Pe = 10⁴. Equations (51), (52), (53) and (58). Subscript A refers to the non-Fickian part of the dimensionless longitudinal macrodispersion coefficient ${D}_{m\mathrm{LA}}^{\ast }$; subscript B refers to the dimensionless Fickian part ${D}_{m\mathrm{LB}}^{\ast }$

Fig. 2

Fickian and non-Fickian transverse dispersion coefficients at Pe = 10⁴. Equations (51), (52), (53) and (58). Subscript A refers to the non-Fickian part of the dimensionless transverse macrodispersion coefficient ${D}_{m\mathrm{TA}}^{\ast }$; subscript B refers to the dimensionless Fickian part ${D}_{m\mathrm{TB}}^{\ast }$

Fig. 3

Fickian and non-Fickian longitudinal dispersion coefficients at Pe = 10⁷. Equations (51), (52), (53) and (58). Subscript A refers to the non-Fickian part of the dimensionless longitudinal macrodispersion coefficient ${D}_{m\mathrm{LA}}^{\ast }$; subscript B refers to the dimensionless Fickian part ${D}_{m\mathrm{LB}}^{\ast }$

Fig. 4

Fickian and non-Fickian transverse dispersion coefficients at Pe = 10⁷. Equations (51), (52), (53) and (58). Subscript A refers to the non-Fickian part of the dimensionless transverse macrodispersion coefficient ${D}_{m\mathrm{TA}}^{\ast }$; subscript B refers to the dimensionless Fickian part ${D}_{m\mathrm{TB}}^{\ast }$

Fig. 5

Non-Fickian longitudinal coefficient α₁₁ over the whole range of scaling exponent (Pe = 10⁴). Black full line identifies the antipersistent range (0 < b < 1); gray full line identifies the persistent range (1 ≤ b < 2)

Fig. 6

Non-Fickian transverse coefficient α₂₂ ≡ α₃₃ over the whole range of scaling exponent (Pe = 10⁴). Black full line identifies the antipersistent range (0 < b < 1); gray full line identifies the persistent range (1 ≤ b < 2)

Fig. 7

Fickian longitudinal coefficient β₁₁ over the whole range of scaling exponent (Pe = 10⁴). Black full line identifies the antipersistent range (0 < b < 1); gray full line identifies the persistent range (1 ≤ b < 2)

Fig. 8

Fickian transverse coefficient β₂₂ ≡ β₃₃ over the whole range of scaling exponent (Pe = 10⁴). Black full line identifies the antipersistent range (0 < b < 1); gray full line identifies the persistent range (1 ≤ b < 2)

Fig. 9

Hierarchical density of the log-conductivity covariance function for antipersistent structures (Eq. (59) with λ = 1 and ϕ = 1). The hierarchical density map in the figure represents the contribution to the evolving-scale log-K covariance from the combination of primary-hierarchy frequency λ = 1 and a finite interval of secondary-hierarchy frequencies 0 ≤ μ ≤ 10

Fig. 10

Hierarchical density of the log-conductivity covariance function for persistent structures (Eq. 60 with λ = 1 and ϕ = 1). The hierarchical density map in the figure represents the contribution to the evolving-scale log-K covariance from the combination of primary-hierarchy frequency λ = 1 and a finite interval of secondary-hierarchy frequencies 0 ≤ μ ≤ 10

Interpretation in terms of risk analysis

As an example of a practical application of the nonergodic macrodispersion coefficients derived in the present study, one may for instance refer to the estimation of the residence time of a nonreactive contaminant (originating from a point pulse) upstream of the area of an urban agglomeration affected by well-withdrawal operations, within an extremely heterogeneous aquifer utilized for drinking supply. Provided that the meaning of the nonergodic macrodispersion coefficient is that of half-time rate of change of the expected central inertia moment,

$${D}_{mii}=\frac{1}{2}\frac{d\left\langle {I}_{ii}\right\rangle }{dt}$$

(61)

and that the related ensemble-mean advection-dispersion equation

$$\frac{\partial \left\langle c\right\rangle }{\partial t}+U\frac{\partial \left\langle c\right\rangle }{\partial {x}_1}=\sum_{i=1}^3\frac{1}{2}\frac{d\left\langle {I}_{ii}\right\rangle }{dt}\frac{\partial^2\left\langle c\right\rangle }{{{\partial x}_i}^2}$$

(62)

is solved, for a point solute pulse of mass M at x = 0 and t = 0, by

$$\left\langle c\left(\textbf{x},t\right)\right\rangle =\frac{M}{n}\prod_{i=1}^3\frac{1}{\sqrt{2\pi \left\langle {I}_{ii}\right\rangle }}\exp \left[-\frac{{\left({x}_i-U{\delta}_{1i}t\right)}^2}{2\left\langle {I}_{ii}\right\rangle}\right]$$

(63)

one can straightforwardly build the (complementary) cumulative probability distribution function of the single-particle position between –∞ and $\overline{l}$ (with $\overline{l}$ indicating the distance of the downstream boundary of the area affected by well-withdrawal operations from the contaminant source). This would require resorting to the well-known statistical definition

$$\varPhi \left(\overline{l},t\right)=1-\frac{n}{M}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\overline{l}}\left\langle c\left({x}_1,{x}_2,{x}_3,t\right)\right\rangle d{x}_1d{x}_2d{x}_3$$

(64)

Note that the probability density function of the single-particle trajectory coincides, for a point solute source, with the normalized (M/n = 1) ensemble-mean concentration (e.g., Pannone and Kitanidis 1999). In dimensionless terms, one obtains:

$$\varPhi \left(l,\tau \right)=1-\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{l-\tau }{\sqrt{\frac{2}{3}k{\tau}^3+\frac{4\tau }{Pe}}}\right)\right]$$

(65)

where erf indicates the error function (Gradshteyn and Ryzhik 1994). Additionally, from the first equations of Eqs. (51) and (52), which allow for the computation of the expected longitudinal central inertia moment 〈I₁₁〉 in the antipersistent and persistent range, respectively:

$$k=\frac{224}{15 Pe}\frac{\pi^{1-b}}{\left(2-b\right)}\frac{\mathrm{B}\left[\frac{\left(3+b\right)}{2},\frac{\left(1-b\right)}{2}\right]}{2^b}\kern2.5em 0<b<1$$

(66)

and

$$k=\frac{224}{15 Pe}\frac{\pi^{3/2}}{\left(2-b\right)}\frac{\Gamma \left[\frac{\left(3+b\right)}{2}\right]}{2^b{\pi}^{b/2}}\kern2.5em 1\le b<2$$

(67)

with

$$l=\overline{l}{\phi}^{1/b}$$

(68)

Practically speaking, Eq. (65) provides, at each dimensionless time, the percentage of mass that, at that time, is already beyond the critical downstream boundary. Small values of Φ (close to zero) indicate that most of the contaminant is still inside the area affected by well-withdrawal, with high water contamination risk. Conversely, high values of Φ (close to 1) indicate that most of the contaminant, thanks to the combined action of advection and dispersion, has already moved out of the supplying area, and that the withdrawal operations would practically be safe. Figures 11, 12 and 13 respectively show, for Pe = 10³, 10⁴ and 10⁵, the time-behavior of function Φ for b = 0.5 and 1.5, and for l = 20, 50 and 100. As one can see, relatively smaller values of Peclét (which trigger a more intense interplay of advective heterogeneity and diffusive mechanisms, with consequent more intense macrodispersion) are responsible for an anomalous phenomenon, singularly related to evolving-scale log-K distributions. As Figure 11 clearly makes evident, for b = 0.5 at a relatively short distance from the contaminant source (l = 20) (and, to a gradually lesser extent, for the other parameter combinations), solute dispersion about the center of mass can progressively become so intense that function Φ, after getting close to 0.9, decreases towards asymptotic values considerably smaller than 1. This means that, despite the advective average downstream translation at rate U, plume backward dispersion can result in the dominant process leading to the entrapment of the contaminant upstream of the safety limit l (highlighted by the infinite right tail in Fig. 11) and to the permanent risk of withdrawal-area slow contamination. As Pe increases (Figs. 12 and 13), the downstream advection takes over and the crossing of boundary l becomes faster and faster. Overall, it is observed that the smaller scaling exponent (b = 0.5) or, equivalently, the higher fractal dimension, determines a more marked sensitivity to increasing Pe, with a faster achievement of the theoretical asymptotic Φ = 1. Finally, as one expects, larger distances from the contaminant source are associated with a crossing time dilatation, as well as with a higher probability of incomplete crossing, due to the longer time allowed for (mainly backward) dispersion.

Fig. 11

Breakthrough-curve analysis: complementary cumulative probability distribution function of the single-particle position for Pe = 10³, b = 0.5 and b = 1.5 (Eq. 65)

Fig. 12

Breakthrough-curve analysis: Complementary cumulative probability distribution function of the single-particle position for Pe = 10⁴, b = 0.5 and b = 1.5 (Eq. 65)

Fig. 13

Breakthrough-curve analysis: Complementary cumulative probability distribution function of the single-particle position for Pe = 10⁵, b = 0.5 and b = 1.5 (Eq. 65)

The use of curves like those appearing in Figs. 11, 12 and 13, for a fast preliminary risk assessment related to a contamination event, requires: (1) knowledge of the log-K statistical structure in terms of experimentally reconstructed ensemble mean 〈Y〉 and semivariogram γ_Y (with definition of scaling coefficient a and exponent b), as well as of the average head gradient/Darcy velocity U; (2) selection of the percentage of acceptable residual resident mass (1-Φ) as a function of the contaminant toxicity level; (3) based on the input datum $\overline{l}$, determination of the corresponding achievement time t = τ/(Uϕ^1/b) where:

$$\phi =-\frac{a}{\Gamma \left(-b\right)}\kern2em 0<b<1$$

(69)

or

$$\phi =-\frac{2a}{\Gamma \left(-\frac{b}{2}\right){\left(\frac{\sqrt{\pi }}{2}\right)}^b}\kern3.75em 1\le b<2$$

(70)

Field validation (Cape Cod tracer test data)

As reported by Leblanc et al. (1991), the aquifer of the Cape Cod (Massachusetts, USA) tracer test was composed of about 100 m of unconsolidated sediments overlying a practically impermeable crystalline bedrock. The upper 30 m of the aquifer consisted of a stratified, sand-gravel outwash; during the tracer test, the bromide tracer plume kept always well inside the sand-gravel aquifer. All the statistics related to the main hydrogeological and hydraulic parameters were reconstructed by taking the needed measurements within that geological unit, which was assumed to be represented by a single-scale, short-range correlation structure: ${\sigma}_Y^2=0.24$, I_Yh = 2.6 m, I_Yv = 0.19 m, 〈K〉 = 110 m/day, U = 0.42 m/day and D = 0.000063 m²/day (with I_Yh and I_Yv, respectively, indicating the horizontal and the vertical log-K integral scale). Beginning on 28 July 1985 and ending on 29 July 1985, a total mass M = 4,900 g of bromide (Br^–) was injected into three 5.08-cm wells located 0.9 m apart along a line perpendicular to the groundwater flow at ~1.2–2.4 m below the water table. A total of 19 local concentration sampling rounds were completed between July 1985 and June 1988 via a multilevel sampler array. After December 1986 (511 days from injection), the edge of the bromide plume had moved out of the array, and the related sampling was stopped. Among other things, local concentration data (from a minimum of 597 samples 13 days after injection to a maximum of 2,270 samples 241 days after injection) allowed for the reconstruction of average trajectory and plume central inertia moments (Garabedian et al. 1991). From then on, many statistical transport theories, whose common denominator was represented by the single-scale, short-range correlation hypothesis, were tested based on these data. However, neither the Eulerian approach by Gelhar and Axness (1983), nor the Lagrangian approach by Dagan (1989), seemed to be able to closely reconstruct the longitudinal dispersivity A_L (i.e., the ratio of longitudinal macrodispersion coefficient to mean velocity) obtained from the field data as the ratio of the average (time) rate of change of the actual longitudinal central inertia moment I_ii to mean velocity: A_L = 0.96 m. Note that, were the Cape Cod aquifer really represented by a single-scale, short-range correlation structure, at large times after injection, ergodic and nonergodic macrodispersion coefficients should tend to coincide, with

$${A}_{\mathrm{L}}=\frac{1}{2U}\frac{d{X}_{ii}}{dt}\cong \frac{1}{2U}\frac{d\left\langle {I}_{ii}\right\rangle }{dt}\cong \frac{1}{2U}\frac{d{I}_{ii}}{dt}$$

(71)

Nevertheless, both the ergodic Eulerian and the Lagrangian theoretical dispersivity considerably underestimate the field value. Specifically, the result of the Eulerian approach by Gelhar and Axness (1983) as reported by Garabedian et al. (1991):

$${A}_{\mathrm{LEul}}=\frac{\sigma_Y^2{I}_{Y\mathrm{h}}e}{\gamma^2\zeta }$$

(72)

yields a longitudinal dispersivity equal to 0.5 m. In Eq. (72), e indicates the anisotropy ratio: e = I_Yv/I_Yh,

$$\gamma =\frac{q_{\mathrm{D}}}{K_{\mathrm{G}}J}$$

(73)

q _D = |q_D| = |nv| indicates Darcy specific discharge, K_G = exp < Y > is the log-K geometric mean, and ζ is a correction factor that accounts for the small angle of flow to the bedding θ = 5°:

$$\zeta =\sqrt{\sin^2\theta +{e}^2{\cos}^2\theta }$$

(74)

The result of the Lagrangian approach by Dagan (1989), expressed by:

$${A}_{\mathrm{LLag}}={\sigma}_Y^2{I}_{Y\mathrm{h}}$$

(75)

leads to a longitudinal dispersivity equal to 0.624 m.

The underestimation of the field value by the theoretical ergodic macrodispersivity is rather surprising since, being inclusive of the uncertainty related to the location of the centroid, it should at most overestimate the real quantity. Note that the application of the present Lagrangian theory in the single-scale version (Pannone and Kitanidis 1999; Pannone and Kitanidis 2004), and in terms of large-time nonergodic macrodispersivity for Gaussian anisotropic covariance,

$${R}_Y\left(\textbf{r}\right)={\sigma}_Y^2\exp \left[-\frac{\pi }{4{I}_{Y\mathrm{h}}^2}\left({r}_1^2+{r}_2^2+\frac{r_3^2}{e^2}\right)\right]$$

(76)

leads to:

$${A}_{\mathrm{L}}=\underset{t\to \infty }{\lim}\frac{1}{2U}\frac{d\left\langle {I}_{ii}\right\rangle }{dt}\cong 0.496$$

(77)

This value, which is consistently very close to the outcome of the ergodic formulation by Gelhar and Axness (1983) as expressed by Eq. (72), was specifically obtained for the present study via the time numerical integrations of Eqs. (27) and (28), which required, for the convergence, a very large (but still computationally possible) number of harmonics. Incidentally, a similar numerical integration was not possible for the investigated case of the evolving-scale uncut log-K spectrum within the limits of an acceptable computational effort.

The idea underlying the present attempt of validation consists in thinking about the outwash sand-gravel layer of the Cape Cod aquifer as an evolving-scale geologic formation with a 30-m cutoff represented by its thickness. In this case, the suitable equation to be used was Eq. (57) (or, equivalently, Fig. 7, in the antipersistent range 0 < b < 1), which expresses the Fickian part of the dimensionless macrodispersion coefficient ${D}_{m\mathrm{LB}}^{\ast }$, with the characteristic length scale being represented by the cutoff: L_max = ϕ^–1/b = 30 m. For b = 0.75, one gets ${D}_{m\mathrm{LB}}^{\ast }$= 0.031 and A_L = ${D}_{m\mathrm{LB}}^{\ast }$ϕ^–1/b = 0.93 m. Based on the theory proposed by the present study, and on a suitably modified version of Eq. (3), the subsequent estimation of the variance that corresponds to the single-scale heterogeneity detected in 1991 at the Cape Cod tracer test site yielded:

$${\left.\frac{d{\sigma}_Y^2}{d\left(\frac{\lambda }{e}\right)}\right|}_{\lambda_{\mathrm{f}}=\frac{1}{I_{Y\mathrm{h}}}}=\frac{\phi }{{\lambda_{\mathrm{f}}}^{1+b}}$$

(78)

and, therefore:

$$d{\sigma}_Y^2\left({\lambda}_{\mathrm{f}}=\frac{1}{I_{Y\mathrm{h}}}\right)=\phi {I}_{Y\mathrm{h}}^{1+b}\frac{d\lambda}{e}\cong \frac{\phi }{e}{I}_{Y\mathrm{h}}^{1+b}\left(\frac{1}{L_{\mathrm{max}}}\right)=0.2$$

(79)

This value is very close to the detected 0.24. Note that in Eq. (78), the application of the present isotropic model, which implies compression of the domain in the horizontal directions as a result of the real-formation anisotropy, imposes a corresponding dilatation of the interval of frequencies contributing to the single-scale hierarchical variance $d{\sigma}_Y^2$ (dλ→ dλ/e). Figure 14 shows the comparison between observed longitudinal central inertia moments and large-time Fickian predictions after Eq. (57).

Fig. 14

Observed vs theoretically estimated (Eq. (57) of the present study with b = 0.75) central longitudinal inertia moments at the Cape Cod tracer test site; t indicates, in days, the time elapsed from tracer injection

As one can see, the evolving-scale-with-cutoff model produces a quite satisfactory agreement between field (single-realization) observations and theoretical (ensemble-mean) predictions. Only at the beginning of the sampling period does the present (large-time) model slightly and systematically overestimate the effective plume central inertia moments.

Discussion

Many authors have debated the existence of geologic power-law structures (e.g., Zech et al. 2015). A recent study by Brunetti et al. (2022) has experimentally shown evidence of sequential ranges of scaling behavior of hydraulic conductivity at the mesoscale, by resorting to pumping tests within highly heterogeneous materials packed into a laboratory tank. Such experiments have also enabled identification of the proper law allowing for the transition from the laboratory to the field scale. As a preliminary field validation, the present study has proposed a comparison between observed longitudinal central inertia moments at the Cape Cod tracer test site (Garabedian et al. 1991) and the theoretical predictions obtained by assuming that the involved aquifer was represented by an evolving-scale structure with vertical cutoff. As a matter of fact, both Eulerian and Lagrangian theories (e.g., Gelhar and Axness 1983; Dagan 1989; Pannone and Kitanidis 1999), in terms of ergodic as well as nonergodic (the present study) estimates for single-scale short-range correlations, fail in predicting the field values, which are in all cases more or less considerably underestimated. The adoption of the evolving-scale-with-cutoff model, in resorting to the theoretical approach proposed by the present study to estimate the expected Cape Cod central inertia moments, practically translated into the use of the Fickian component of the fractal macrodispersion coefficient, which was indeed obtained by integration from a nonzero wave-number—the inverse of the characteristic length-scale—to infinity. The suitable characteristic length scale could in this case be identified with the spatial cutoff, that is, with the sand-gravel outwash thickness. The scaling exponent that better fitted to the observed plume spatial moments time-series was b = 0.75 (antipersistent range), leading to a very good agreement with the theoretical predictions. The single-scale variance detected in the field was reconstructed starting from the theoretical expression of the primary-hierarchy variance density, with a suitable modification represented by the dilatation of the horizontal frequency domain as a consequence of the anisotropy-related compression of the Cartesian horizontal coordinates. The field-detected variance was ${\sigma}_Y^2=0.24$; the theoretical estimate was equal to ~0.2.

As widely explained, the key point of the proposed analytical approach was represented by the decomposition of the evolving-scale log-conductivity distribution into a double hierarchy of stationary fields (see Eqs. 59 and 60) characterized by variable-amplitude, pseudo-periodic or periodic covariances. Previous studies had resorted to a single-hierarchy decomposition, involving exponential or Gaussian log-K covariance. The double hierarchy allows one to go into detail in the interpretation of antipersistent and persistent correlations, and their role in determining the dispersive power of the related flow field. In summary, each field of the primary hierarchy (identified by wave-number λ) is given by the superposition of: (1) cosinusoidal slower-than-periodic fields whose variability is distributed over a wider range of secondary-hierarchy wave-number μ for 0 < b < 1 (antipersistent structures) and (2) cosinusoidal periodic fields whose variability is distributed over a narrower range (with values relatively closer to zero) of secondary-hierarchy wave-number μ for 1 ≤ b < 2 (persistent structures). For that reason, and far from the extremes of b-value in both ranges, where the large amplitude of the oscillations dominates due to the strong contribution of smallest and largest scales of heterogeneity, respectively, the antipersistent structures are associated with a more marked dispersion about the center of mass, thanks to the greater high/medium frequency content. As already inferred by Pannone (2017), diffusion acts as an antagonist mechanism for the achievement of Fickian dispersion conditions. Indeed, as shown by the present study, being the longitudinal Fickian (constant) part of the large-time macrodispersion coefficient practically independent of the Peclét number, when this last becomes very large (almost purely advective regime), and the coefficient of the non-Fickian time-dependent component tends to zero, the process should in principle tend to become globally diffusive. Nevertheless, the value of Pe for which this would happen is always unrealistically high; therefore, Fickian conditions are practically never achieved. Finally, another important point that needs to be emphasized is that the use made by the present investigation of a length-scale (ϕ^–1/b) related to the specific log-conductivity/velocity correlation range, in contrast with what was done by previous studies focusing on the same topic, consistently leads to a universal (quadratic) dependence of the nonergodic macrodispersion coefficients on the dimensionless time.

Conclusions

The protection of groundwater resources and attempts to limit or slow-down their unavoidable deterioration is one of the crucial objectives in the field of environmental engineering. In order to construct efficient prediction and monitoring tools, many researchers in the last few decades have addressed, via different types of deterministic or stochastic approaches, water flow and solute transport in geologic formations. The use of stochastic methodologies is often suggested by the typically marked heterogeneity of groundwater flow fields, which considerably complicates the whole theoretical analysis, tying it to the definition of a specific space-time reference scale, and to the introduction of several simplifying assumptions. Classic porous media stochastic theories assume that the log-conductivity Y(x) = ln K(x), where K is the local hydraulic conductivity and x is the vector of spatial coordinates, is a statistically stationary and normally distributed random function, completely characterized by a constant mean and fast-decaying correlation function. Statistical stationarity implies the existence of a single representative scale of heterogeneity, i.e., the so-called integral scale. In the case of regional aquifers, solute transport is unavoidably influenced by several scales of structural variability, which progressively come into play as the travel distance increases. In these conditions, the semivariogram of Y (i.e., the variance of log-K spatial increments) may suitably be represented by a simple-scaling power-law model. The present study has analyzed, by a first-order Lagrangian approach and for instantaneous point injections, the effect of evolving-scale log-conductivity distributions on nonergodic tracer dispersion in the presence of realistically weak but nonzero diffusion. In contrast with what was found by previous Lagrangian studies that focused on finite-dimension solute injections, the main result of the investigation, which provides large-time closed-form solutions for the nonergodic macrodispersion coefficients, was that, overall, antipersistent correlations (scaling exponent between 0 and 1) lead to a more intense dispersion about the center of mass, mathematically represented by a universal (quadratic) dependence on the dimensionless time built with a characteristic heterogeneity length. Specifically, the nonergodic macrodispersion coefficients resulted in the sum of a quadratically time-dependent, non-Fickian component (coming from the largest scales of the heterogeneity) and a constant Fickian component (obtained from large to infinitely small scales). The Fickian component exhibited a clearly nonmonotonic behavior as a function of the scaling exponent, with a well-defined minimum, asymmetrically located in the central part of the antipersistent range, and a practically vanishing contribution in the persistent one. Conversely, the non-Fickian component coefficient (the multiplier of the square dimensionless time) slightly decreases around the lower extremes and considerably increases around the higher extremes of both ranges. The predictive ability of the Fickian component, which can in itself be thought of as referred to an evolving-scale-with-cutoff log-K distribution, was successfully tested against the experimental observations at the Cape Cod tracer test site. Additionally, the interpretation of the derived macrodispersion coefficients in terms of breakthrough curves (BTCs) allowed for highlighting a very peculiar behavior of the evolving-scale log-conductivity structures. This very peculiar behavior consists of the association of relatively-low Peclét number processes with practically infinite right tails of BTCs, and, due to dominant backward dispersion, with the permanent risk of slow contamination of the aquifer area upstream of the selected breakthrough point. There is no doubt that the role played by flow and transport nonlinearity in determining solute “fractal” dispersion deserves a specific and detailed investigation, although the author’s feeling is that it may likely only accelerate and intensify the trends and the dynamics that have emerged from the present study. This may happen by a self-feeding process in which the particle sampling of larger and larger scales of heterogeneity is obviously faster.

Appendix: On the representativeness of the ensemble mean central inertia moments

Here, the generic single-realization central inertia moment (i.e., the central inertia moment obtained by spatial averages), is considered:

$${I}_{ii}(t)=\frac{1}{V}{\int}_V\left[{\overset{\sim }{X}}_i\left(t,\textbf{a}\right)-\frac{1}{V}{\int}_V{\overset{\sim }{X}}_i\left(t,{\textbf{a}}^{\prime}\right)d{\textbf{a}}^{\prime}\right]\left[{\overset{\sim }{X}}_i\left(t,\textbf{a}\right)-\frac{1}{V}{\int}_V{\overset{\sim }{X}}_i\left(t,{\textbf{a}}^{\prime \prime}\right)d{\textbf{a}}^{\prime \prime}\right]d\textbf{a}$$

(80)

or

$${I}_{ii}(t)={I}_{ii}\left(t- dt\right)+\frac{1}{V}{\int}_V\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)d\textbf{a}-\frac{1}{V^2}{\int}_V{\int}_V\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)\overset{\sim }{X}^{\prime}_i\left(t,{\textbf{a}}^{\prime}\right)d\textbf{a}d{\textbf{a}}^{\prime}\kern0.5em$$

(81)

where V = V(t) is the instantaneous plume volume, a, a′ and a′′ indicate the particles’ starting position at time t – dt,

$${I}_{ii}\left(t- dt\right)=\frac{1}{V}{\int}_V{\left({a}_i-\frac{1}{V}{\int}_V{a}_i^{\prime }d{\textbf{a}}^{\prime}\right)}^2d\textbf{a}$$

(82)

is the corresponding central inertia moment,

$${\overset{\sim }{X}}_i\left(t,\textbf{a}\right)=\left\langle {\overset{\sim }{X}}_i\left(t,\textbf{a}\right)\right\rangle +\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)={a}_i+U{\delta}_{1i} dt+\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)$$

(83)

identifies particle position at time t, and

$${S}_i(t)=\frac{1}{V}{\int}_V{\overset{\sim }{X}}_i\left(t,\textbf{a}\right)d\textbf{a}$$

(84)

indicates the single-realization centroid. The single-realization Eq. (81) can thus be read as:

$${I}_{ii}(t)={I}_{ii}\left(t- dt\right)+{X}_{ii}^r(t)-{S}_{ii}^r(t)$$

(85)

with:

$${X}_{ii}^r(t)=\frac{1}{V}{\int}_V\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)d\textbf{a}$$

(86)

indicating the single-realization one-particle variance and

$${S}_{ii}^r(t)=\frac{1}{V^2}{\int}_V{\int}_V\overset{\sim }{X}^{\prime}_i\left(t,\textbf{a}\right)\overset{\sim }{X}^{\prime}_i\left(t,{\textbf{a}}^{\prime}\right)d\textbf{a}d{\textbf{a}}^{\prime }$$

(87)

indicating the single-realization centroid variance. The corresponding ensemble mean and deviation are, respectively:

$$\left\langle {I}_{ii}(t)\right\rangle =\left\langle {I}_{ii}\left(t- dt\right)\right\rangle +\left\langle {X}_{ii}^r(t)\right\rangle -\left\langle {S}_{ii}^r(t)\right\rangle =\left\langle {I}_{ii}\left(t- dt\right)\right\rangle +{X}_{ii}(t)-{S}_{ii}(t)$$

(88)

and:

$${I}_{ii}^{\prime }(t)={I}_{ii}(t)-\left\langle {I}_{ii}(t)\right\rangle ={I}_{ii}^{\prime}\left(t- dt\right)+{X}_{ii}^{\prime }(t)-{S}_{ii}^{\prime }(t)$$

(89)

Note that the deviation is essentially due to the failure of the ergodic hypothesis (plume volume at finite times is never large enough to sample the whole heterogeneous log-K distribution and to make spatial and statistical mean equivalent). The failure of the ergodic hypothesis, in turn, is determined by the largest scales of the heterogeneity (larger than plume dimensions). After Eq. (19), one can write:

$${\overset{\sim }{X}}_i^{\prime}\left(t,\textbf{a}\right)=\sum\limits_{\lambda =0}^{\infty }\sum\limits_{\mu =0}^{\infty }{X}_i^{\prime}\left(t,\textbf{a};\lambda, \mu \right)=\sum\limits_{\lambda =0}^{\overline{\lambda}(t)}\sum\limits_{\mu =0}^{\infty }{X}_i^{\prime}\left(t,\textbf{a};\lambda, \mu \right)+\sum\limits_{\lambda =\overline{\lambda}(t)}^{\infty }\sum\limits_{\mu =0}^{\infty }{X}_i^{\prime}\left(t,\textbf{a};\lambda, \mu \right)={X}_I^{\prime}\left(t,\textbf{a}\right)+{X}_{II}^{\prime}\left(t,\textbf{a}\right)$$

(90)

with ${X}_I^{\prime }$ which is associated to such large scales ($\overline{\lambda}$ is the wave-length that corresponds to the instantaneous representative plume dimension). Due to the related large range of correlation, the ${X}_I^{\prime }$-field “sees” the particles starting at different locations a and a′ within volume V (t – dt), as they practically were the same particle. In these conditions:

$${S}_{ii}^{\prime }(t)=\frac{1}{V^2}{\int}_V{\int}_V{X}_{iI}\left(t,\textbf{a}\right){X}_{iI}\left(t,{\textbf{a}}^{\prime}\right)d\textbf{a}d{\textbf{a}}^{\prime}\cong \frac{1}{V^2}{\int}_V{\int}_V{X}_{iI}\left(t,\textbf{a}\right){X}_{iI}\left(t,\textbf{a}\right)d\textbf{a}d{\textbf{a}}^{\prime }=\frac{1}{V}{\int}_V{X}_{iI}\left(t,\textbf{a}\right){X}_{iI}\left(t,\textbf{a}\right)d\textbf{a}={X}_{ii}^{\prime }(t)$$

(91)

that is:

$${X}_{ii}^{\prime }(t)\cong {S}_{ii}^{\prime }(t)$$

(92)

Since similar considerations apply at any t (and, therefore, also at t – dt):

$${I}_{ii}^{\prime}\left(t- dt\right)\cong 0\kern0.5em \mathrm{and}\kern1em {I}_{ii}^{\prime }(t)\cong 0$$

(93)

Therefore, the single realization inertia moments, in an evolving-scale log-K field such as those considered in the present study, are almost fully represented by the corresponding ensemble mean.