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

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.


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 Published in the special issue "Geostatistics and hydrogeology" * Marilena Pannone marilena.pannone@unibas.it 1 School of Engineering -University of Basilicata, 10 Ateneo Lucano Avenue, 85100 Potenza, Italy 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 superdiffusive 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 (2010Pannone ( , 2012Pannone ( , 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 translationallyinvariant 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 Y = 2 Y − R Y (with 2 Y 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 secondaryhierarchy periods): and where r = √ x 2 1 + x 2 2 + x 2 3 indicates the norm of the distance, Γ the Gamma function (Gradshteyn and Ryzhik 1994), ϕ is a dimensional proportionality constant, and (1) identifies the variance density associated with the primaryhierarchy 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 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: In the hypothesis of reasonably small ϕ, even the larger heterogeneous subunits (μ→0, λ→0) will be characterized by a reasonably small variance: and 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): In Eq. (7), <h(x)> is the hydraulic head ensemble mean, h′ the corresponding fluctuation, and is the constant head mean gradient vector. In the Fourier domain (wave-number vector k), Eq. (7) becomes: where j = √ −1 , k is the integration variable, and the circumflex accent indicates Fourier transforms: 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: From Eq. (12), and from first-order Darcy's law expressed in terms of fluctuating (or deviatoric) quantities: one obtains in the Fourier domain (Dagan 1989): where v i is the Fourier transform of the ith component of velocity fluctuation and U = (U 1 , U 2 , U 3 ) = 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 Ŷ and ĥ equal to zero by definition, one obtains, for k ≠k , ⟨Ŷ k −k ĥ k ⟩ = 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: where t is the time and X B the zero-mean Brownian component. Additionally, indicates the ensemble mean particle position vector, 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): 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: 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): indicates the isotropic local dispersion coefficient and ∼ 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 c (e.g., Fischer et al. 1979). The reason why it is accounted for in Eq. (20) and (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 singlescale expression by Pannone and Kitanidis 1999): where f XB indicates the Gaussian probability density function of the single Brownian displacement: f XBZB is the joint probability density function of the two independent Brownian displacements: is the hierarchical covariance of the ith velocity component: 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): and with: and k = √ k 2 1 + k 2 2 + k 2 3 . 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): 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: and 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: and, therefore, the nonergodic macrodispersion coefficient: 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: 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: where: and F 11 ( , ) = 1 − sin 2 cos 2 2 F 22 ( , ) = sin 4 cos 2 sin 2 F 33 ( , ) = sin 2 cos 2 cos 2 (39) = 4 2 Dk 2 1 + exp −8 2 Dk 2 t − 2 exp −4 2 Dk 2 t cos (2 ksin cos Ut) 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 2 φ, cos 2 φ, and the even cosφ-function cos(2πk sin θ cos φUt), Eq. (34) becomes: with: Pe cos (2 k * sin cos ) k * 2 + Pe 2 4 2 sin 2 cos 2 k * 3+b d d dk *

Results
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 ( ) = D * miiA ( ) + D * miiB ( ) = ii 2 + ii , the dependence of non-Fickian (α ii ) and Fickian (β ii ) coefficients on b at Pe = 10 4 . 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: 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 � √ r � , 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.

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 halftime rate of change of the expected central inertia moment, and that the related ensemble-mean advection-dispersion equation is solved, for a point solute pulse of mass M at x = 0 and t = 0, by one can straightforwardly build the (complementary) cumulative probability distribution function of the single-particle Note that the probability density function of the singleparticle 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: 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 11 〉 in the antipersistent and persistent range, respectively: 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 3 , 10 4 and 10 5 , 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 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.
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 l , determination of the corresponding achievement time t = τ/ (Uϕ 1/b ) where: Fig. 8 Fickian transverse coefficient β 22 ≡ β 33 over the whole range of scaling exponent (Pe = 10 4 ). 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 secondaryhierarchy 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 primaryhierarchy frequency λ = 1 and a finite interval of secondary-hierarchy frequencies 0 ≤ μ ≤ 10

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: 2 Y = 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 2 /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 . 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 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): yields a longitudinal dispersivity equal to 0.5 m. In Eq. (72), e indicates the anisotropy ratio: e = I Yv /I Yh , 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°: The result of the Lagrangian approach by Dagan (1989), expressed by: 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, leads to: 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 evolvingscale 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, variance that corresponds to the single-scale heterogeneity detected in 1991 at the Cape Cod tracer test site yielded: and, therefore: 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 2 Y (dλ→ dλ/e). Figure 14 shows the comparison between observed longitudinal central inertia moments and large-time Fickian predictions after Eq. (57).
As one can see, the evolving-scale-with-cutoff model produces a quite satisfactory agreement between field (singlerealization) 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.  ) and the theoretical predictions obtained by assuming that the involved aquifer was represented by an evolvingscale 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 anisotropyrelated compression of the Cartesian horizontal coordinates. The field-detected variance was 2 Y = 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 wavenumber μ 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 welldefined 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 evolvingscale-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 relativelylow 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.