Solute transport in aquifers with evolving scale heterogeneity

Nicolae Suciu
(preprint), Numerical Modeling, paper

Abstract

Transport processes in groundwater systems with spatially heterogeneous properties often exhibit anomalous behavior.

Using first-order approximations in velocity fluctuations we show that anomalous superdiffusive behavior may result if velocity fields are modeled as superpositions of random space functions with correlation structures consisting of linear combinations of short-range correlations.

In particular, this corresponds to the superposition of independent random velocity fields with increasing integral scales proposed as model for evolving scale heterogeneity of natural porous media [Gelhar, L. W. Water Resour. Res. 22 (1986), 135S-145S].

Monte Carlo simulations of transport in such multi-scale fields support the theoretical results and demonstrate the approach to superdiffusive behavior as the number of superposed scales increases.

Authors

N. Suciu
-Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

S. Attinger

F. A.  Radu

C. Vamoș
-Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

J. Vanderborght

H. Vereecken

P. Knabner

Keywords

Porous media; Random fields; transport; random walk

Cite this paper as:

N. Suciu, S. Attinger, F.A. Radu, C. Vamoş, J. Vanderborght, H. Vereecken, P. Knabner, Solute transport in aquifers with evolving scale heterogeneity, Preprint No. 346, 2011, Institute of Applied Mathematics, Friedrich-Alexander University Erlangen-Nuremberg

References

see the expansion block below.

PDF

pr346.pdf

About this paper

Journal

Preprint

Publisher Name
DOI
Print ISSN
Online ISSN

MR

?

ZBL

?

[1] Attinger, S., Dentz, M., H. Kinzelbach, and W. Kinzelbach (1999), Temporal behavior of a solute cloud in a chemically heterogeneous porous medium, J. Fluid Mech., 386, 77-104.
CrossRef (DOI)

[2] Bellin, A., M. Pannone, A. Fiori, and A. Rinaldo (1996), On transport in porous formations characterized by heterogeneity of evolving scales, Water Resour. Res., 32, 3485-3496.
CrossRef (DOI)

[3] Cintoli, S., S. P. Neuman, and V. Di Federico (2005), Generating and scaling fractional Brownian motion on finite domains, Geophys. Res. Lett. 32, L08404,
CrossRef (DOI).

[4] Dagan, G. (1994), The significance of heterogeneity of evolving scales and of anomalous diffusion to transport in porous formations, Water Resour. Res., 30, 3327-3336, 1994.
CrossRef (DOI)

[5] Dagan, G. (1987), Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, 183-215.
CrossRef (DOI)

[6] Dagan, G. (1988), Time-dependent macrodispersion for solute transport in anisotropic heterogeneous aquifers, Water Resour. Res., 24, 1491-1500.
CrossRef (DOI)

[7] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000), Temporal behavior of a solute cloud in a heterogeneous porous medium 1. Point-like injection, Water Resour. Res., 36, 3591-3604.
CrossRef (DOI)

[8] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000), Temporal behavior of a solute cloud in a heterogeneous porous medium 2. Spatially extended injection, Water Resour. Res., 36, 3605-3614.
CrossRef (DOI)

[9] Di Federico, V., and S. P. Neuman (1997), Scaling of random fields by means of truncated power variograms and associated spectra, Water Resour. Res., 33, 1075-1085.
CrossRef (DOI)

[10] Fiori, A. (1996), Finite Peclet extensions of Dagan’s solutions to transport in anisotropic heterogeneous formations, Water Resour. Res., 32, 193-198.
CrossRef (DOI)

[11] Fiori, A. (2001), On the influence of local dispersion in solute transport through formations with evolving scales of heterogeneity, Water Resour. Res., 37, 235-242.
CrossRef (DOI)

[12] Fiori, A., and G. Dagan (2000), Concentration fluctuations in aquifer transport: A rigorous first-order solution and applications,  J. Contam. Hydrol., 45, 139-163.
CrossRef (DOI)

[13] Gelhar, L. W. (1986), Stochastic subsurface hydrology from theory to applications, Water Resour. Res., 22, 135S-145S.
CrossRef (DOI)

[14] Gelhar, L. W., and C. L. Axness (1983), Three-dimensional stochastic analysis of macrodispersion in aquifers, textit Water  Resour. Res., 19, 161-180.
CrossRef (DOI)

[15] Gradshteyn, I. S., and I. M. Ryzhik (2007), Table of Integrals, Series, and Products, Elsevier, Amsterdam.

[16] Jeon, J.-H., and R. Metzler (2010), Fractional Brownian motion and motion governed by the fractional Langevin equation in confined geometries, Phys. Rev. E 81, 021103,
CrossRef (DOI)

[17] McLaughlin, D., and F. Ruan (2001), Macrodispersivity and large-scale hydrogeologic variability, Transp. Porous Media, 42, 133-154.
CrossRef (DOI)

[18] Papoulis, A., and S. U. Pillai (2009), Probability, Random Variables and Stochastic Processes, McGraw-Hill, New York.

[19] Radu, F. A., N. Suciu, J. Hoffmann, A. Vogel, O. Kolditz, C-H. Park, and S. Attinger (2011), Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study, Adv. Water Resour. 34, 47–61.
CrossRef (DOI)

[20] Ross, K., and S. Attinger (2010), Temporal behaviour of a solute cloud in a fractal heterogeneous porous medium at different scales, paper presented at EGU General Assembly 2010, Vienna, Austria, 02-07 May 2010.

[21] Schwarze, H., U. Jaekel, and H. Vereecken (2001), Estimation of macrodispersivity by different approximation methods for flow and transport in randomly heterogeneous media, Transp. Porous Media, 43, 265-287.

[22] Suciu N. (2010), Spatially inhomogeneous transition probabilities as memory effects for diffusion in statistically homogeneous random velocity fields, Phys. Rev. E, 81, 056301,
CrossRef (DOI)

[23] Suciu, N. (2014), Diffusion in random velocity fields with applications to contaminant transport in groundwater, Adv. Water Resour. 69, 114–133.
CrossRef (DOI)

[24] Suciu, N., C. Vamo¸s, J. Vanderborght, H. Hardelauf, and H. Vereecken (2006), Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409,
CrossRef (DOI)

[25] Suciu N., C. Vamos, F. A. Radu, H. Vereecken, and P. Knabner (2009), Persistent memory of diffusing particles, Phys. Rev. E, 80, 061134,
CrossRef (DOI)

[26] Suciu, N., S. Attinger, F.A. Radu, C. Vamos, J. Vanderborght, H. Vereecken, P. Knabner (2011), Solute transport in aquifers with evolving scale heterogeneity, Preprint No. 346, Mathematics Department, Friedrich-Alexander University Erlangen-Nuremberg
(http://fauams5.am.uni-erlan-gen.de/papers/pr346.pdf).
CrossRef (DOI)

[27] Suciu, N., F.A. Radu, A. Prechtel, F. Brunner, and P. Knabner (2013), A coupled finite element-global random walk approach to advectiondominated transport in porous media with random hydraulic conductivity, J. Comput. Appl. Math. 246, 27–37.
CrossRef (DOI)

[28] Suciu, N., F.A. Radu, S. Attinger, L. Schuler, Knabner (2014), A FokkerPlanck approach for probability distributions of species concentrations transported in heterogeneous media, J. Comput. Appl. Math., in press,
CrossRef (DOI)

[29] Vamos, C., N. Suciu, H. Vereecken, J. Vanderborght, and O. Nitzsche (2001), Path decomposition of discrete effective diffusion coefficient, Internal Report ICG-IV. 00501, Research Center Jülich.

[30] Vamos, C., N. Suciu, and H. Vereecken (2003), Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comput. Phys., 186, 527-544,
CrossRef (DOI)

soon

Paper (preprint) in HTML form

pr346

Solute transport in aquifers with evolving scale heterogeneity

Nicolae Suciu 1 , 3 1 , 3 ^(1,3)quad{ }^{1,3} \quad1,3 Sabine Attinger 2 2 ^(2)quad{ }^{2} \quad2 Florin A. Radu 2 2 ^(2)quad{ }^{2} \quad2 Călin Vamoş 3 3 ^(3){ }^{3}3Jan Vanderborght 4 4 ^(4){ }^{4}4 Harry Vereecken 4 4 ^(4){ }^{4}4 Peter Knabner 1 1 ^(1){ }^{1}1 1 1 ^(1){ }^{1}1 Chair for Applied Mathematics I, Friedrich-Alexander University Erlangen-Nuremberg, Germany, suciu@am.uni-erlangen.de, knabner@am.uni-erlangen.de 2 2 ^(2){ }^{2}2 Computational Hydrosystems, Helmholtz Center for Environmental Research - UFZ, Germany, sabine.attinger@ufz.de, florin.radu@ufz.de 3 3 ^(3){ }^{3}3 Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj Napoca, Romania, cvamos@ictp.acad.ro, nsuciu@ictp.acad.ro 4 4 ^(4){ }^{4}4 Agrosphere Institute ICG-IV, Research Center Jülich, Germany, h. vereecken@fz-juelich. de, J.vanderborght@fz-juelich.de

Abstract

A superposition of independent log-K fields with increasing integral scales has been proposed by many authors [ 13 , 4 , 9 , 17 , 3 ] [ 13 , 4 , 9 , 17 , 3 ] [13,4,9,17,3][13,4,9,17,3][13,4,9,17,3] as model for evolving scale heterogeneity of natural porous media. Here, we investigate the anomalous transport in multi-scale fields by first-order approximations to transport equations and numerical Monte Carlo simulations.

1 Lagrangian and Eulerian models of diffusion in random fields and their relationship

Let { X i , t = X i ( t ) , i = 1 , 2 , 3 , t ≥ 0 } X i , t = X i ( t ) , i = 1 , 2 , 3 , t ≥ 0 {X_(i,t)=X_(i)(t),i=1,2,3,t >= 0}\left\{X_{i, t}=X_{i}(t), i=1,2,3, t \geq 0\right\}{Xi,t=Xi(t),i=1,2,3,t≥0} be the advection-dispersion process with a space variable drift V ( x ) V ( x ) V(x)\mathbf{V}(\mathbf{x})V(x) which is a sample of a random velocity field, and a local diffusion coefficient D D DDD which is assumed constant, described by the Itô equation
(1) X i , t = x i , 0 + ∫ 0 t V i ( X t ′ ) d t ′ + W i , t (1) X i , t = x i , 0 + ∫ 0 t   V i X t ′ d t ′ + W i , t {:(1)X_(i,t)=x_(i,0)+int_(0)^(t)V_(i)(X_(t^(')))dt^(')+W_(i,t):}\begin{equation*} X_{i, t}=x_{i, 0}+\int_{0}^{t} V_{i}\left(\mathbf{X}_{t^{\prime}}\right) d t^{\prime}+W_{i, t} \tag{1} \end{equation*}(1)Xi,t=xi,0+∫0tVi(Xt′)dt′+Wi,t
where x i , 0 x i , 0 x_(i,0)x_{i, 0}xi,0 is a deterministic initial position and W i W i W_(i)W_{i}Wi are the components of a Wiener process of mean zero and variance ⟨ W i , t 2 ⟩ = 2 D t W i , t 2 = 2 D t (:W_(i,t)^(2):)=2Dt\left\langle W_{i, t}^{2}\right\rangle=2 D t⟨Wi,t2⟩=2Dt.
The solution g g ggg of the Fokker-Planck equation
(2) ∂ t g + u ∇ g = D ∇ 2 g (2) ∂ t g + u ∇ g = D ∇ 2 g {:(2)del_(t)g+ugrad g=Dgrad^(2)g:}\begin{equation*} \partial_{t} g+\mathbf{u} \nabla g=D \nabla^{2} g \tag{2} \end{equation*}(2)∂tg+u∇g=D∇2g
which corresponds to the Green's function in Eulerian approaches to dispersion in random fields, is the density of the transition probability P ( d x , t , x 0 ) = g ( x , t ∣ x 0 ) d x P d x , t , x 0 = g x , t ∣ x 0 d x P(dx,t,x_(0))=g(x,t∣x_(0))dxP\left(d \mathbf{x}, t, \mathbf{x}_{0}\right)=g\left(\mathbf{x}, t \mid \mathbf{x}_{0}\right) d \mathbf{x}P(dx,t,x0)=g(x,t∣x0)dx of the Itô process (1). Equation (2) is the starting point of the perturbation approach of [ 1 , 7 , 8 ] [ 1 , 7 , 8 ] [1,7,8][1,7,8][1,7,8].
The following three centered processes are useful to describe the effective and ensemble dispersion and the fluctuations of the center of mass: X i , t e f f = X i , t − ⟨ X i , t ⟩ X i , t e f f = X i , t − X i , t X_(i,t)^(eff)=X_(i,t)-(:X_(i,t):)X_{i, t}^{e f f}=X_{i, t}-\left\langle X_{i, t}\right\rangleXi,teff=Xi,t−⟨Xi,t⟩, the process centered about the mean ⟨ X i , t ⟩ X i , t (:X_(i,t):)\left\langle X_{i, t}\right\rangle⟨Xi,t⟩ of X i , t X i , t X_(i,t)X_{i, t}Xi,t in single realizations of the random velocity field, ⟨ X i , t e f f ⟩ = 0 X i , t e f f = 0 (:X_(i,t)^(eff):)=0\left\langle X_{i, t}^{e f f}\right\rangle=0⟨Xi,teff⟩=0, the process X i , t e n s = X i , t − ⟨ X i , t ⟩ ― X i , t e n s = X i , t − X i , t ¯ X_(i,t)^(ens)=X_(i,t)- bar((:X_(i,t):))X_{i, t}^{e n s}= X_{i, t}-\overline{\left\langle X_{i, t}\right\rangle}Xi,tens=Xi,t−⟨Xi,t⟩― centered about the average ⟨ X i , t ⟩ ― X i , t ¯ bar((:X_(i,t):))\overline{\left\langle X_{i, t}\right\rangle}⟨Xi,t⟩― of ⟨ X i , t ⟩ X i , t (:X_(i,t):)\left\langle X_{i, t}\right\rangle⟨Xi,t⟩ over the ensemble of realizations, ⟨ X i , t e n s ⟩ ― = 0 X i , t e n s ¯ = 0 bar((:X_(i,t)^(ens):))=0\overline{\left\langle X_{i, t}^{e n s}\right\rangle}=0⟨Xi,tens⟩―=0,
and the fluctuation of the center of mass X i , t c m = ⟨ X i , t ⟩ − ⟨ X i , t ⟩ ― X i , t c m = X i , t − X i , t ¯ X_(i,t)^(cm)=(:X_(i,t):)- bar((:X_(i,t):))X_{i, t}^{c m}=\left\langle X_{i, t}\right\rangle-\overline{\left\langle X_{i, t}\right\rangle}Xi,tcm=⟨Xi,t⟩−⟨Xi,t⟩―, with X i , t c m ― = 0 X i , t c m ¯ = 0 bar(X_(i,t)^(cm))=0\overline{X_{i, t}^{c m}}=0Xi,tcm―=0. Associated to these processes, one defines [1]:
\left.\left.\begin{array}{rlrl} S_{i i}(t) & =\overline{\left\langle\left(X_{t}^{e f f}\right)^{2}\right\rangle}, & & D_{i i}^{e f f}(t) \end{array}=\frac{1}{2} \frac{d S(t)}{d t}, \text { effective dispersion coefficients }, 75\right)=\frac{1}{2} \frac{d \Sigma(t)}{d t}, \text { ensemble dispersion coefficients }\right] \begin{aligned} \Sigma_{i i}(t)_{i i} & =\overline{\left\langle\left(X_{t}^{e n s}\right)^{2}\right\rangle}, \tag{3}\\ R_{i i}(t) & =\overline{\left\langle\left(X_{t}^{c m}\right)^{2}\right\rangle}, \end{aligned} \quad D_{i i}^{e n s}(t)=\frac{1}{2} \frac{d R(t)}{d t}=D^{e n s}-D^{e f f}, \text { center of mass dispersion coefficients }\tag not allowed in aligned environment
According to (1), for a process starting at x i , 0 = 0 x i , 0 = 0 x_(i,0)=0x_{i, 0}=0xi,0=0 we have
(4) X i , t = ∫ 0 t u i ( X t ′ ) d t ′ + W i , t (4) X i , t = ∫ 0 t   u i X t ′ d t ′ + W i , t {:(4)X_(i,t)=int_(0)^(t)u_(i)(X_(t^(')))dt^(')+W_(i,t):}\begin{equation*} X_{i, t}=\int_{0}^{t} u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right) d t^{\prime}+W_{i, t} \tag{4} \end{equation*}(4)Xi,t=∫0tui(Xt′)dt′+Wi,t
where u i = V i − ⟨ V i ⟩ u i = V i − V i u_(i)=V_(i)-(:V_(i):)u_{i}=V_{i}-\left\langle V_{i}\right\rangleui=Vi−⟨Vi⟩ for X i , t = X i , t e f f X i , t = X i , t e f f X_(i,t)=X_(i,t)^(eff)X_{i, t}=X_{i, t}^{e f f}Xi,t=Xi,teff and u i = V i − ⟨ V i ⟩ ― u i = V i − V i ¯ u_(i)=V_(i)- bar((:V_(i):))u_{i}=V_{i}-\overline{\left\langle V_{i}\right\rangle}ui=Vi−⟨Vi⟩― for X i , t = X i , t e n s X i , t = X i , t e n s X_(i,t)=X_(i,t)^(ens)X_{i, t}=X_{i, t}^{e n s}Xi,t=Xi,tens. From (4), the variance of X i , t = X i , t ens X i , t = X i , t ens  X_(i,t)=X_(i,t)^("ens ")X_{i, t}=X_{i, t}^{\text {ens }}Xi,t=Xi,tens  is obtained as follows,
Σ i i ( t ) i i = ∫ 0 t ∫ 0 t ⟨ u i ( X t ′ ) u i ( X t ′ ′ ) ⟩ ― d t ′ d t ′ ′ + 2 ∫ 0 t ⟨ u i ( X t ′ ) W i , t ⟩ ― d t ′ + ⟨ W i , t 2 ⟩ = ∫ 0 t ∫ 0 t ⟨ u i ( X t ′ ) u i ( X t ′ ′ ) ⟩ ― d t ′ d t ′ ′ + 2 ∫ 0 t ⟨ u i ( X t ′ ) ― W i , t ⟩ d t ′ + ⟨ W i , t 2 ⟩ = 2 D t + 2 ∫ 0 t d t ′ ∫ 0 t ′ ⟨ u i ( X t ′ ) u i ( X t ′ ′ ) ⟩ ― d t ′ ′ Σ i i ( t ) i i = ∫ 0 t   ∫ 0 t   u i X t ′ u i X t ′ ′ ¯ d t ′ d t ′ ′ + 2 ∫ 0 t   u i X t ′ W i , t ¯ d t ′ + W i , t 2 = ∫ 0 t   ∫ 0 t   u i X t ′ u i X t ′ ′ ¯ d t ′ d t ′ ′ + 2 ∫ 0 t   u i X t ′ ¯ W i , t d t ′ + W i , t 2 = 2 D t + 2 ∫ 0 t   d t ′ ∫ 0 t ′   u i X t ′ u i X t ′ ′ ¯ d t ′ ′ {:[Sigma_(ii)(t)_(ii)=int_(0)^(t)int_(0)^(t) bar((:u_(i)(X_(t^(')))u_(i)(X_(t^(''))):))dt^(')dt^('')+2int_(0)^(t) bar((:u_(i)(X_(t^(')))W_(i,t):))dt^(')+(:W_(i,t)^(2):)],[{:=int_(0)^(t)int_(0)^(t) bar((:u_(i)(X_(t^(')))u_(i)(X_(t^(''))):))dt^(')dt^('')+2int_(0)^(t) bar((:u_(i)(X_(t^('))))W_(i,t):)dt^(')+(:W_(i,t)^(2):)],[=2Dt+2int_(0)^(t)dt^(')int_(0)^(t^(')) bar((:u_(i)(X_(t^(')))u_(i)(X_(t^(''))):))dt^('')]:}\begin{aligned} \Sigma_{i i}(t)_{i i} & =\int_{0}^{t} \int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime \prime}}\right)\right\rangle} d t^{\prime} d t^{\prime \prime}+2 \int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right) W_{i, t}\right\rangle} d t^{\prime}+\left\langle W_{i, t}^{2}\right\rangle \\ & \left.=\int_{0}^{t} \int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime \prime}}\right)\right\rangle} d t^{\prime} d t^{\prime \prime}+2 \int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right)\right.} W_{i, t}\right\rangle d t^{\prime}+\left\langle W_{i, t}^{2}\right\rangle \\ & =2 D t+2 \int_{0}^{t} d t^{\prime} \int_{0}^{t^{\prime}} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime \prime}}\right)\right\rangle} d t^{\prime \prime} \end{aligned}Σii(t)ii=∫0t∫0t⟨ui(Xt′)ui(Xt′′)⟩―dt′dt′′+2∫0t⟨ui(Xt′)Wi,t⟩―dt′+⟨Wi,t2⟩=∫0t∫0t⟨ui(Xt′)ui(Xt′′)⟩―dt′dt′′+2∫0t⟨ui(Xt′)―Wi,t⟩dt′+⟨Wi,t2⟩=2Dt+2∫0tdt′∫0t′⟨ui(Xt′)ui(Xt′′)⟩―dt′′
where we used ⟨ u i ⟩ ― = 0 u i ¯ = 0 bar((:u_(i):))=0\overline{\left\langle u_{i}\right\rangle}=0⟨ui⟩―=0 and the independence of the Wiener process on velocity statistics, which, together, cancel the middle term. Further, the definitions (3) yield
(5) D i i e n s ( t ) = D + ∫ 0 t ⟨ u i ( X t ) u i ( X t ′ ) ⟩ ― d t ′ (6) D i i c m ( t ) = ∫ 0 t ⟨ u i ( X t ) ⟩ ⟨ u i ( X t ′ ) ⟩ ― d t ′ (7) D i i e f f ( t ) = D i , i e n s ( t ) − D i , i c m ( t ) (5) D i i e n s ( t ) = D + ∫ 0 t   u i X t u i X t ′ ¯ d t ′ (6) D i i c m ( t ) = ∫ 0 t   u i X t u i X t ′ ¯ d t ′ (7) D i i e f f ( t ) = D i , i e n s ( t ) − D i , i c m ( t ) {:[(5)D_(ii)^(ens)(t)=D+int_(0)^(t) bar((:u_(i)(X_(t))u_(i)(X_(t^('))):))dt^(')],[(6)D_(ii)^(cm)(t)=int_(0)^(t) bar((:u_(i)(X_(t)):)(:u_(i)(X_(t^('))):))dt^(')],[(7)D_(ii)^(eff)(t)=D_(i,i)^(ens)(t)-D_(i,i)^(cm)(t)]:}\begin{align*} D_{i i}^{e n s}(t) & =D+\int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right)\right\rangle} d t^{\prime} \tag{5}\\ D_{i i}^{c m}(t) & =\int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}}\right)\right\rangle\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right)\right\rangle} d t^{\prime} \tag{6}\\ D_{i i}^{e f f}(t) & =D_{i, i}^{e n s}(t)-D_{i, i}^{c m}(t) \tag{7} \end{align*}(5)Diiens(t)=D+∫0t⟨ui(Xt)ui(Xt′)⟩―dt′(6)Diicm(t)=∫0t⟨ui(Xt)⟩⟨ui(Xt′)⟩―dt′(7)Diieff(t)=Di,iens(t)−Di,icm(t)
The first iteration of (4) about the un-perturbed problem consisting of diffusion with constant coefficient D D DDD in the mean flow field U = ⟨ V i ⟩ ― U = V i ¯ U= bar((:V_(i):))U=\overline{\left\langle V_{i}\right\rangle}U=⟨Vi⟩―, with trajectory X i , t ( 0 ) = δ j , 1 U t + W i , t X i , t ( 0 ) = δ j , 1 U t + W i , t X_(i,t)^((0))=delta_(j,1)Ut+W_(i,t)X_{i, t}^{(0)}=\delta_{j, 1} U t+W_{i, t}Xi,t(0)=δj,1Ut+Wi,t, yields:
D i i e n s ( t ) − D = ∫ 0 t ⟨ u i ( X t ( 0 ) ) u i ( X t ′ ( 0 ) ) ⟩ ― d t ′ = ∫ 0 t d t ′ ∬ ⟨ u i ( x ) δ ( x − X t ( 0 ) ) u i ( x ′ ) δ ( x ′ − X t ′ ( 0 ) ) ⟩ ― d x d x ′ (8) = ∫ 0 t d t ′ ∬ u i ( x ) u i ( x ′ ) ― p ( x , t ; x ′ , t ′ ) d x d x ′ D i i e n s ( t ) − D = ∫ 0 t   u i X t ( 0 ) u i X t ′ ( 0 ) ¯ d t ′ = ∫ 0 t   d t ′ ∬ u i ( x ) δ x − X t ( 0 ) u i x ′ δ x ′ − X t ′ ( 0 ) ¯ d x d x ′ (8) = ∫ 0 t   d t ′ ∬ u i ( x ) u i x ′ ¯ p x , t ; x ′ , t ′ d x d x ′ {:[D_(ii)^(ens)(t)-D=int_(0)^(t) bar((:u_(i)(X_(t)^((0)))u_(i)(X_(t^('))^((0))):))dt^(')],[=int_(0)^(t)dt^(')∬ bar((:u_(i)(x)delta(x-X_(t)^((0)))u_(i)(x^('))delta(x^(')-X_(t^('))^((0))):))dxdx^(')],[(8)=int_(0)^(t)dt^(')∬ bar(u_(i)(x)u_(i)(x^(')))p(x,t;x^('),t^('))dxdx^(')]:}\begin{align*} D_{i i}^{e n s}(t)-D & =\int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}}^{(\mathbf{0})}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}^{(\mathbf{0})}\right)\right\rangle} d t^{\prime} \\ & =\int_{0}^{t} d t^{\prime} \iint \overline{\left\langle u_{i}(\mathbf{x}) \delta\left(\mathbf{x}-\mathbf{X}_{\mathbf{t}}^{(\mathbf{0})}\right) u_{i}\left(\mathbf{x}^{\prime}\right) \delta\left(\mathbf{x}^{\prime}-\mathbf{X}_{\mathbf{t}^{\prime}}^{(\mathbf{0})}\right)\right\rangle} d \mathbf{x} d \mathbf{x}^{\prime} \\ & =\int_{0}^{t} d t^{\prime} \iint \overline{u_{i}(\mathbf{x}) u_{i}\left(\mathbf{x}^{\prime}\right)} p\left(\mathbf{x}, t ; \mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x} d \mathbf{x}^{\prime} \tag{8} \end{align*}Diiens(t)−D=∫0t⟨ui(Xt(0))ui(Xt′(0))⟩―dt′=∫0tdt′∬⟨ui(x)δ(x−Xt(0))ui(x′)δ(x′−Xt′(0))⟩―dxdx′(8)=∫0tdt′∬ui(x)ui(x′)―p(x,t;x′,t′)dxdx′
where p ( x , t ; x ′ , t ′ ) = ⟨ δ ( x − X t ( 0 ) ) δ ( x ′ − X t ′ ( 0 ) ) ⟩ p x , t ; x ′ , t ′ = δ x − X t ( 0 ) δ x ′ − X t ′ ( 0 ) p(x,t;x^('),t^('))=(:delta(x-X_(t)^((0)))delta(x^(')-X_(t^('))^((0))):)p\left(\mathbf{x}, t ; \mathbf{x}^{\prime}, t^{\prime}\right)=\left\langle\delta\left(\mathbf{x}-\mathbf{X}_{\mathbf{t}}^{(\mathbf{0})}\right) \delta\left(\mathbf{x}^{\prime}-\mathbf{X}_{\mathbf{t}^{\prime}}^{(\mathbf{0})}\right)\right\ranglep(x,t;x′,t′)=⟨δ(x−Xt(0))δ(x′−Xt′(0))⟩ is the Gaussian joint probability density of the process X i , t ( 0 ) X i , t ( 0 ) X_(i,t)^((0))X_{i, t}^{(0)}Xi,t(0). Similarly,
(9) D i i c m = ∫ 0 t d t ′ ∬ u i ( x ) u i ( x ′ ) ― p ( x , t ) p ( x ′ , t ′ ) d x d x ′ (9) D i i c m = ∫ 0 t   d t ′ ∬ u i ( x ) u i x ′ ¯ p ( x , t ) p x ′ , t ′ d x d x ′ {:(9)D_(ii)^(cm)=int_(0)^(t)dt^(')∬ bar(u_(i)(x)u_(i)(x^(')))p(x","t)p(x^('),t^('))dxdx^('):}\begin{equation*} D_{i i}^{c m}=\int_{0}^{t} d t^{\prime} \iint \overline{u_{i}(\mathbf{x}) u_{i}\left(\mathbf{x}^{\prime}\right)} p(\mathbf{x}, t) p\left(\mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x} d \mathbf{x}^{\prime} \tag{9} \end{equation*}(9)Diicm=∫0tdt′∬ui(x)ui(x′)―p(x,t)p(x′,t′)dxdx′
Further, D i i e f f ( t ) D i i e f f ( t ) D_(ii)^(eff)(t)D_{i i}^{e f f}(t)Diieff(t) is obtained from (7). Together, (7-9) are the equivalent (stochastic) Lagrangian formulation of the first-order approximation which, in an Eulerian frame, has been presented in [1,7,8].
Remark 1. As follows from (8-9), if u i ( x ) u i ( x ) u_(i)(x)u_{i}(\mathbf{x})ui(x) is a superposition of independent velocity fields with different integral scales, then the dispersion coefficients estimated to the first-order are given by the sum of the coefficients associated to the terms of the velocity superposition.

2 Quantifying anomalous diffusion by memory terms of the ensemble, effective, and center of mass dispersion coefficients

In the following, by X t X t X_(t)X_{t}Xt we denote one of the components of X i , t X i , t X_(i,t)X_{i, t}Xi,t as well as those of the processes X i , t e f f , X i , t e n s X i , t e f f , X i , t e n s X_(i,t)^(eff),X_(i,t)^(ens)X_{i, t}^{e f f}, X_{i, t}^{e n s}Xi,teff,Xi,tens, and X i , t c m X i , t c m X_(i,t)^(cm)X_{i, t}^{c m}Xi,tcm. Further, we consider the uniform time partition 0 < Δ < ⋯ < k Δ < ⋯ < ( K − 1 ) Δ < K Δ 0 < Δ < ⋯ < k Δ < ⋯ < ( K − 1 ) Δ < K Δ 0 < Delta < cdots < k Delta < cdots < (K-1)Delta < K Delta0<\Delta<\cdots<k \Delta< \cdots<(K-1) \Delta<K \Delta0<Δ<⋯<kΔ<⋯<(K−1)Δ<KΔ and for k > 0 k > 0 k > 0k>0k>0 we define the increments X k , k − 1 = X t k − X t k − 1 X k , k − 1 = X t k − X t k − 1 X_(k,k-1)=X_(t_(k))-X_(t_(k-1))X_{k, k-1}=X_{t_{k}}-X_{t_{k-1}}Xk,k−1=Xtk−Xtk−1 and s k , k − 1 = ( X k , k − 1 ) 2 s k , k − 1 = X k , k − 1 2 s_(k,k-1)=(X_(k,k-1))^(2)s_{k, k-1}=\left(X_{k, k-1}\right)^{2}sk,k−1=(Xk,k−1)2. For a sum of two increments X k , 0 = X k − 1 , 0 + X k , k − 1 X k , 0 = X k − 1 , 0 + X k , k − 1 X_(k,0)=X_(k-1,0)+X_(k,k-1)X_{k, 0}=X_{k-1,0}+X_{k, k-1}Xk,0=Xk−1,0+Xk,k−1 we have
(10) s k , 0 = s k − 1 , 0 + s k , k − 1 + m k , k − 1 , 0 (10) s k , 0 = s k − 1 , 0 + s k , k − 1 + m k , k − 1 , 0 {:(10)s_(k,0)=s_(k-1,0)+s_(k,k-1)+m_(k,k-1,0):}\begin{equation*} s_{k, 0}=s_{k-1,0}+s_{k, k-1}+m_{k, k-1,0} \tag{10} \end{equation*}(10)sk,0=sk−1,0+sk,k−1+mk,k−1,0
where m k , k − 1 , 0 = 2 X k − 1 , 0 X k , k − 1 m k , k − 1 , 0 = 2 X k − 1 , 0 X k , k − 1 m_(k,k-1,0)=2X_(k-1,0)X_(k,k-1)m_{k, k-1,0}=2 X_{k-1,0} X_{k, k-1}mk,k−1,0=2Xk−1,0Xk,k−1. By iterating the binomial formula (10) one obtains:
s 1 , 0 = s 1 , 0 s 2 , 0 = s 1 , 0 + s 2 , 1 + m 2 , 1 , 0 … … … … … … … s k , 0 = s k − 1 , 0 + s k , k − 1 + m k , k − 1 , 0 … … … … … … … s K − 1 , 0 = s K − 2 , 0 + s K − 1 , K − 2 + m K − 1 , K − 2 , 0 s K , 0 = s K − 1 , 0 + s K , K − 1 + m K , K − 1 , 0 s 1 , 0 = s 1 , 0 s 2 , 0 = s 1 , 0 + s 2 , 1 + m 2 , 1 , 0 … … … … … … … s k , 0 = s k − 1 , 0 + s k , k − 1 + m k , k − 1 , 0 … … … … … … … s K − 1 , 0 = s K − 2 , 0 + s K − 1 , K − 2 + m K − 1 , K − 2 , 0 s K , 0 = s K − 1 , 0 + s K , K − 1 + m K , K − 1 , 0 {:[s_(1,0)=s_(1,0)],[s_(2,0)=s_(1,0)+s_(2,1)+m_(2,1,0)],[dots dots dots dots dots dots dots],[s_(k,0)=s_(k-1,0)+s_(k,k-1)+m_(k,k-1,0)],[dots dots dots dots dots dots dots],[s_(K-1,0)=s_(K-2,0)+s_(K-1,K-2)+m_(K-1,K-2,0)],[s_(K,0)=s_(K-1,0)+s_(K,K-1)+m_(K,K-1,0)]:}\begin{aligned} s_{1,0} & =s_{1,0} \\ s_{2,0} & =s_{1,0}+s_{2,1}+m_{2,1,0} \\ \ldots & \ldots \ldots \ldots \ldots \ldots \ldots \\ s_{k, 0} & =s_{k-1,0}+s_{k, k-1}+m_{k, k-1,0} \\ \ldots & \ldots \ldots \ldots \ldots \ldots \ldots \\ s_{K-1,0} & =s_{K-2,0}+s_{K-1, K-2}+m_{K-1, K-2,0} \\ s_{K, 0} & =s_{K-1,0}+s_{K, K-1}+m_{K, K-1,0} \end{aligned}s1,0=s1,0s2,0=s1,0+s2,1+m2,1,0…………………sk,0=sk−1,0+sk,k−1+mk,k−1,0…………………sK−1,0=sK−2,0+sK−1,K−2+mK−1,K−2,0sK,0=sK−1,0+sK,K−1+mK,K−1,0
Summing up these relations we get
(11) s K , 0 = ∑ k = 1 K s k , k − 1 + ∑ k = 1 K m k , k − 1 , 0 , (11) s K , 0 = ∑ k = 1 K   s k , k − 1 + ∑ k = 1 K   m k , k − 1 , 0 , {:(11)s_(K,0)=sum_(k=1)^(K)s_(k,k-1)+sum_(k=1)^(K)m_(k,k-1,0)",":}\begin{equation*} s_{K, 0}=\sum_{k=1}^{K} s_{k, k-1}+\sum_{k=1}^{K} m_{k, k-1,0}, \tag{11} \end{equation*}(11)sK,0=∑k=1Ksk,k−1+∑k=1Kmk,k−1,0,
which, using the relation
m k , k − 1 , 0 = 2 X k − 1 , 0 X k , k − 1 = 2 X k , k − 1 ∑ l = 1 k − 1 X l , l − 1 m k , k − 1 , 0 = 2 X k − 1 , 0 X k , k − 1 = 2 X k , k − 1 ∑ l = 1 k − 1   X l , l − 1 m_(k,k-1,0)=2X_(k-1,0)X_(k,k-1)=2X_(k,k-1)sum_(l=1)^(k-1)X_(l,l-1)m_{k, k-1,0}=2 X_{k-1,0} X_{k, k-1}=2 X_{k, k-1} \sum_{l=1}^{k-1} X_{l, l-1}mk,k−1,0=2Xk−1,0Xk,k−1=2Xk,k−1∑l=1k−1Xl,l−1
can be further expressed as
(12) s K , 0 = ∑ k = 1 K s k , k − 1 + 2 ∑ k = 1 K ∑ l = 1 k − 1 X k − 1 , 0 X l , l − 1 , (12) s K , 0 = ∑ k = 1 K   s k , k − 1 + 2 ∑ k = 1 K   ∑ l = 1 k − 1   X k − 1 , 0 X l , l − 1 , {:(12)s_(K,0)=sum_(k=1)^(K)s_(k,k-1)+2sum_(k=1)^(K)sum_(l=1)^(k-1)X_(k-1,0)X_(l,l-1)",":}\begin{equation*} s_{K, 0}=\sum_{k=1}^{K} s_{k, k-1}+2 \sum_{k=1}^{K} \sum_{l=1}^{k-1} X_{k-1,0} X_{l, l-1}, \tag{12} \end{equation*}(12)sK,0=∑k=1Ksk,k−1+2∑k=1K∑l=1k−1Xk−1,0Xl,l−1,
For X t = X t ens X t = X t ens  X_(t)=X_(t)^("ens ")X_{t}=X_{t}^{\text {ens }}Xt=Xtens , the expectation S k , 0 = ⟨ s k , 0 ⟩ ― S k , 0 = s k , 0 ¯ S_(k,0)= bar((:s_(k,0):))S_{k, 0}=\overline{\left\langle s_{k, 0}\right\rangle}Sk,0=⟨sk,0⟩― of (10) verifies
(13) Σ k , 0 = Σ k − 1 , 0 + Σ k , k − 1 + M k , k − 1 , 0 . (13) Σ k , 0 = Σ k − 1 , 0 + Σ k , k − 1 + M k , k − 1 , 0 . {:(13)Sigma_(k,0)=Sigma_(k-1,0)+Sigma_(k,k-1)+M_(k,k-1,0).:}\begin{equation*} \Sigma_{k, 0}=\Sigma_{k-1,0}+\Sigma_{k, k-1}+M_{k, k-1,0} . \tag{13} \end{equation*}(13)Σk,0=Σk−1,0+Σk,k−1+Mk,k−1,0.
For continuous time t = k Δ t = k Δ t=k Deltat=k \Deltat=kΔ, (13) becomes
(14) Σ t , 0 = Σ t − Δ , 0 + Σ t , t − Δ + M t , t − Δ , 0 . (14) Σ t , 0 = Σ t − Δ , 0 + Σ t , t − Δ + M t , t − Δ , 0 . {:(14)Sigma_(t,0)=Sigma_(t-Delta,0)+Sigma_(t,t-Delta)+M_(t,t-Delta,0).:}\begin{equation*} \Sigma_{t, 0}=\Sigma_{t-\Delta, 0}+\Sigma_{t, t-\Delta}+M_{t, t-\Delta, 0} . \tag{14} \end{equation*}(14)Σt,0=Σt−Δ,0+Σt,t−Δ+Mt,t−Δ,0.
According to the definition of the ensemble coefficient in (3), the half time derive of (13) yields
(15) D t , 0 e n s = D t , t − Δ e n s + M D t , t − Δ , 0 e n s (15) D t , 0 e n s = D t , t − Δ e n s + M D t , t − Δ , 0 e n s {:(15)D_(t,0)^(ens)=D_(t,t-Delta)^(ens)+MD_(t,t-Delta,0)^(ens):}\begin{equation*} D_{t, 0}^{e n s}=D_{t, t-\Delta}^{e n s}+M D_{t, t-\Delta, 0}^{e n s} \tag{15} \end{equation*}(15)Dt,0ens=Dt,t−Δens+MDt,t−Δ,0ens
The last term of (14) M t , t − Δ , 0 = ⟨ m t , t − Δ , 0 ⟩ ― M t , t − Δ , 0 = m t , t − Δ , 0 ¯ M_(t,t-Delta,0)= bar((:m_(t,t-Delta,0):))M_{t, t-\Delta, 0}=\overline{\left\langle m_{t, t-\Delta, 0}\right\rangle}Mt,t−Δ,0=⟨mt,t−Δ,0⟩― is a memory term quantifying the departure from normaldiffusion behavior of the process X t ens [ 23 , 21 ] X t ens  [ 23 , 21 ] X_(t)^("ens ")[23,21]X_{t}^{\text {ens }}[23,21]Xtens [23,21] and its half derivative M D t , t − Δ , 0 ens = ( 1 / 2 ) d M t , t − Δ , 0 / d t M D t , t − Δ , 0 ens  = ( 1 / 2 ) d M t , t − Δ , 0 / d t MD_(t,t-Delta,0)^("ens ")=(1//2)dM_(t,t-Delta,0)//dtM D_{t, t-\Delta, 0}^{\text {ens }}=(1 / 2) d M_{t, t-\Delta, 0} / d tMDt,t−Δ,0ens =(1/2)dMt,t−Δ,0/dt is a memory dispersion coefficient. Relations similar to (14-15) also hold for S , D e f f S , D e f f S,D^(eff)S, D^{e f f}S,Deff, and for R , D c m R , D c m R,D^(cm)R, D^{c m}R,Dcm. The expectation of the second term of Eq. (11),
(16) C M t = ∑ k = 1 t / Δ ⟨ m k Δ , ( k − 1 ) Δ , 0 ⟩ ― , (16) C M t = ∑ k = 1 t / Δ   m k Δ , ( k − 1 ) Δ , 0 ¯ , {:(16)CM_(t)=sum_(k=1)^(t//Delta) bar((:m_(k Delta,(k-1)Delta,0):))",":}\begin{equation*} C M_{t}=\sum_{k=1}^{t / \Delta} \overline{\left\langle m_{k \Delta,(k-1) \Delta, 0}\right\rangle}, \tag{16} \end{equation*}(16)CMt=∑k=1t/Δ⟨mkΔ,(k−1)Δ,0⟩―,
is a cumulated memory term which, according to Eq. (12), contains a hierarchy of correlations between increments along paths of the transport process. Such cumulated memory terms, on paths of length 3 Δ 3 Δ 3Delta3 \Delta3Δ were estimated from a "global random walk" (GRW) simulation of transport in a single realization of the velocity field by [24].
According to (15) and (5), the memory dispersion coefficient corresponding to the increment of the process X t X t X_(t)X_{t}Xt over the time interval τ τ tau\tauτ has the following Lagrangian expression
M D i i e n s ( t , τ , 0 ) = D i i e n s ( t , 0 ) − D i i e n s ( t , τ ) = ∫ 0 t ⟨ u i ( X t ) u i ( X t ′ ) ⟩ ― d t ′ − ∫ τ t ⟨ u i ( X t ) u i ( X t ′ ) ⟩ ― d t ′ (17) = ∫ 0 τ ⟨ u i ( X t ) u i ( X t ′ ) ⟩ ― d t ′ M D i i e n s ( t , τ , 0 ) = D i i e n s ( t , 0 ) − D i i e n s ( t , τ ) = ∫ 0 t   u i X t u i X t ′ ¯ d t ′ − ∫ τ t   u i X t u i X t ′ ¯ d t ′ (17) = ∫ 0 τ   u i X t u i X t ′ ¯ d t ′ {:[MD_(ii)^(ens)(t","tau","0)=D_(ii)^(ens)(t","0)-D_(ii)^(ens)(t","tau)],[=int_(0)^(t) bar((:u_(i)(X_(t))u_(i)(X_(t^('))):))dt^(')-int_(tau)^(t) bar((:u_(i)(X_(t))u_(i)(X_(t^('))):))dt^(')],[(17)=int_(0)^(tau) bar((:u_(i)(X_(t))u_(i)(X_(t^('))):))dt^(')]:}\begin{align*} M D_{i i}^{e n s}(t, \tau, 0) & =D_{i i}^{e n s}(t, 0)-D_{i i}^{e n s}(t, \tau) \\ & =\int_{0}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right)\right\rangle} d t^{\prime}-\int_{\tau}^{t} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right)\right\rangle} d t^{\prime} \\ & =\int_{0}^{\tau} \overline{\left\langle u_{i}\left(\mathbf{X}_{\mathbf{t}}\right) u_{i}\left(\mathbf{X}_{\mathbf{t}^{\prime}}\right)\right\rangle} d t^{\prime} \tag{17} \end{align*}MDiiens(t,τ,0)=Diiens(t,0)−Diiens(t,τ)=∫0t⟨ui(Xt)ui(Xt′)⟩―dt′−∫τt⟨ui(Xt)ui(Xt′)⟩―dt′(17)=∫0τ⟨ui(Xt)ui(Xt′)⟩―dt′
The first order approximation also yields
(18) M D i i e n s ( t , τ , 0 ) = ∫ 0 τ d t ′ ∬ u i ( x ) u i ( x ′ ) ― p ( x , t ; x ′ , t ′ ) d x d x ′ (18) M D i i e n s ( t , τ , 0 ) = ∫ 0 τ   d t ′ ∬ u i ( x ) u i x ′ ¯ p x , t ; x ′ , t ′ d x d x ′ {:(18)MD_(ii)^(ens)(t","tau","0)=int_(0)^(tau)dt^(')∬ bar(u_(i)(x)u_(i)(x^(')))p(x,t;x^('),t^('))dxdx^('):}\begin{equation*} M D_{i i}^{e n s}(t, \tau, 0)=\int_{0}^{\tau} d t^{\prime} \iint \overline{u_{i}(\mathbf{x}) u_{i}\left(\mathbf{x}^{\prime}\right)} p\left(\mathbf{x}, t ; \mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x} d \mathbf{x}^{\prime} \tag{18} \end{equation*}(18)MDiiens(t,τ,0)=∫0τdt′∬ui(x)ui(x′)―p(x,t;x′,t′)dxdx′
Memory coefficients (18) as well as their time integrals (i.e. memory terms) were computed in [23] for short range correlations of the velocity field and for diffusion in perfectly stratified fields (model of Matheron and de Marsily). In [21], the time behavior of the memory terms for fractional Gaussian noise has been also analyzed. Jeon and Metzler [16] used increments correlations to quantify the memory of two special anomalous diffusion processes (fractional Brownian motion and processes governed by fractional Langevin equation). Their correlations are in fact memory terms related with the variance of the process by Eq. (13). For instance, in case of fractional Brownian motion with variance Σ t 2 , t 1 = ( t 2 − t 1 ) 2 H Σ t 2 , t 1 = t 2 − t 1 2 H Sigma_(t_(2),t_(1))=(t_(2)-t_(1))^(2H)\Sigma_{t_{2}, t_{1}}= \left(t_{2}-t_{1}\right)^{2 H}Σt2,t1=(t2−t1)2H and Hurst coefficient 0 < H < 1 0 < H < 1 0 < H < 10<H<10<H<1, for two successive time intervals Δ Δ Delta\DeltaΔ, from (14) one obtains
(19) M Δ = 2 ⟨ X Δ , 0 X 2 Δ , Δ ⟩ ― = Σ 2 Δ , 0 − Σ Δ , 0 − Σ 2 Δ , Δ = ( 2 2 H − 2 ) Δ 2 H (19) M Δ = 2 X Δ , 0 X 2 Δ , Δ ¯ = Σ 2 Δ , 0 − Σ Δ , 0 − Σ 2 Δ , Δ = 2 2 H − 2 Δ 2 H {:(19)M_(Delta)=2 bar((:X_(Delta,0)X_(2Delta,Delta):))=Sigma_(2Delta,0)-Sigma_(Delta,0)-Sigma_(2Delta,Delta)=(2^(2H)-2)Delta^(2H):}\begin{equation*} M_{\Delta}=2 \overline{\left\langle X_{\Delta, 0} X_{2 \Delta, \Delta}\right\rangle}=\Sigma_{2 \Delta, 0}-\Sigma_{\Delta, 0}-\Sigma_{2 \Delta, \Delta}=\left(2^{2 H}-2\right) \Delta^{2 H} \tag{19} \end{equation*}(19)MΔ=2⟨XΔ,0X2Δ,Δ⟩―=Σ2Δ,0−ΣΔ,0−Σ2Δ,Δ=(22H−2)Δ2H
which is just the Eq. (4.4) of Jeon and Metzler [16].

3 Explicit first-order results for given statistics of the log-K field

The ensemble coefficient (8) and the memory coefficient (18) can be further expresses as function of the covariance C Y C Y C_(Y)C_{Y}CY of the log-hydraulic conductivity Y = ln ⁡ K Y = ln ⁡ K Y=ln KY=\ln KY=ln⁡K in the frame of approximation of flow and transport equations to the first-order in the variance σ l n K 2 σ l n K 2 sigma_(lnK)^(2)\sigma_{l n K}^{2}σlnK2. With the common convention for the Fourier transform
(20) u ^ j ( k ) = ∫ u j ( x ) exp ⁡ ( − i k ⋅ x ) d x (21) u j ( x ) = 1 ( 2 π ) d ∫ u ^ j ( k ) exp ⁡ ( i k ⋅ x ) d k , j = 1 ⋯ d , d = 2 , 3 (20) u ^ j ( k ) = ∫ u j ( x ) exp ⁡ ( − i k ⋅ x ) d x (21) u j ( x ) = 1 ( 2 π ) d ∫ u ^ j ( k ) exp ⁡ ( i k ⋅ x ) d k , j = 1 ⋯ d , d = 2 , 3 {:[(20) hat(u)_(j)(k)=intu_(j)(x)exp(-ik*x)dx],[(21)u_(j)(x)=(1)/((2pi)^(d))int hat(u)_(j)(k)exp(ik*x)dk","quad j=1cdots d","d=2","3]:}\begin{align*} \hat{u}_{j}(\mathbf{k}) & =\int u_{j}(\mathbf{x}) \exp (-i \mathbf{k} \cdot \mathbf{x}) d \mathbf{x} \tag{20}\\ u_{j}(\mathbf{x}) & =\frac{1}{(2 \pi)^{d}} \int \hat{u}_{j}(\mathbf{k}) \exp (i \mathbf{k} \cdot \mathbf{x}) d \mathbf{k}, \quad j=1 \cdots d, d=2,3 \tag{21} \end{align*}(20)u^j(k)=∫uj(x)exp⁡(−ik⋅x)dx(21)uj(x)=1(2π)d∫u^j(k)exp⁡(ik⋅x)dk,j=1⋯d,d=2,3
the first-order approximation yields (e.g., [20, 7, 8])
(22) u ^ j ( k ) u ^ j ( k ′ ) ― = ( 2 π ) d U 2 δ ( k + k ′ ) p j 2 ( k ) C ^ Y ( k ) , (22) u ^ j ( k ) u ^ j k ′ ¯ = ( 2 π ) d U 2 δ k + k ′ p j 2 ( k ) C ^ Y ( k ) , {:(22) bar(hat(u)_(j)(k) hat(u)_(j)(k^(')))=(2pi)^(d)U^(2)delta(k+k^('))p_(j)^(2)(k) hat(C)_(Y)(k)",":}\begin{equation*} \overline{\hat{u}_{j}(\mathbf{k}) \hat{u}_{j}\left(\mathbf{k}^{\prime}\right)}=(2 \pi)^{d} U^{2} \delta\left(\mathbf{k}+\mathbf{k}^{\prime}\right) p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}), \tag{22} \end{equation*}(22)u^j(k)u^j(k′)―=(2π)dU2δ(k+k′)pj2(k)C^Y(k),
where p j ( k ) = ( δ j , 1 − k 1 k j / k 2 ) p j ( k ) = δ j , 1 − k 1 k j / k 2 p_(j)(k)=(delta_(j,1)-k_(1)k_(j)//k^(2))p_{j}(\mathbf{k})=\left(\delta_{j, 1}-k_{1} k_{j} / k^{2}\right)pj(k)=(δj,1−k1kj/k2) are projectors which ensure the incompressibility of the flow.
Remark 2. The Dirac function δ ( k + k ′ ) δ k + k ′ delta(k+k^('))\delta\left(\mathbf{k}+\mathbf{k}^{\prime}\right)δ(k+k′) which ensures the statistical homogeneity in (22) has to be changed into δ ( k − k ′ ) δ k − k ′ delta(k-k^('))\delta\left(\mathbf{k}-\mathbf{k}^{\prime}\right)δ(k−k′) to obtain a similar first-order relation for u ^ ― j ( k ) u ^ j ∗ ( k ′ ) u ^ ¯ j ( k ) u ^ j ∗ k ′ bar(hat(u))_(j)(k) hat(u)_(j)^(**)(k^('))\overline{\hat{u}}_{j}(\mathbf{k}) \hat{u}_{j}^{*}\left(\mathbf{k}^{\prime}\right)u^―j(k)u^j∗(k′), where u ^ j ∗ u ^ j ∗ hat(u)_(j)^(**)\hat{u}_{j}^{*}u^j∗ denotes the complex conjugate [5,6,20]. Even though it leads to the same results, this relation render the computation a bit tedious.
Using the inverse Fourier transform (21) and (22), the Eulerian velocity correlation from (8) and (18) can be computed as follows
u j ( x ) u j ( x ′ ) ― = 1 ( 2 π ) 2 d ∬ u ^ j ( k ) u ^ j ( k ′ ) ― exp ⁡ ( i k ⋅ x ) exp ⁡ ( i k ′ ⋅ x ′ ) d k d k ′ = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ ( i k ⋅ x ) d k ∫ δ ( k + k ′ ) exp ⁡ ( i k ′ ⋅ x ′ ) d k ′ (23) = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k ⋅ ( x − x ′ ) ] d k u j ( x ) u j x ′ ¯ = 1 ( 2 π ) 2 d ∬ u ^ j ( k ) u ^ j k ′ ¯ exp ⁡ ( i k ⋅ x ) exp ⁡ i k ′ ⋅ x ′ d k d k ′ = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ ( i k ⋅ x ) d k ∫ δ k + k ′ exp ⁡ i k ′ ⋅ x ′ d k ′ (23) = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k ⋅ x − x ′ d k {:[ bar(u_(j)(x)u_(j)(x^(')))=(1)/((2pi)^(2d))∬ bar(hat(u)_(j)(k) hat(u)_(j)(k^(')))exp(ik*x)exp(ik^(')*x^('))dkdk^(')],[=(U^(2))/((2pi)^(d))intp_(j)^(2)(k) hat(C)_(Y)(k)exp(ik*x)dkint delta(k+k^('))exp(ik^(')*x^('))dk^(')],[(23)=(U^(2))/((2pi)^(d))intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik*(x-x^('))]dk]:}\begin{align*} \overline{u_{j}(\mathbf{x}) u_{j}\left(\mathbf{x}^{\prime}\right)} & =\frac{1}{(2 \pi)^{2 d}} \iint \overline{\hat{u}_{j}(\mathbf{k}) \hat{u}_{j}\left(\mathbf{k}^{\prime}\right)} \exp (i \mathbf{k} \cdot \mathbf{x}) \exp \left(i \mathbf{k}^{\prime} \cdot \mathbf{x}^{\prime}\right) d \mathbf{k} d \mathbf{k}^{\prime} \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp (i \mathbf{k} \cdot \mathbf{x}) d \mathbf{k} \int \delta\left(\mathbf{k}+\mathbf{k}^{\prime}\right) \exp \left(i \mathbf{k}^{\prime} \cdot \mathbf{x}^{\prime}\right) d \mathbf{k}^{\prime} \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i \mathbf{k} \cdot\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\right] d \mathbf{k} \tag{23} \end{align*}uj(x)uj(x′)―=1(2π)2d∬u^j(k)u^j(k′)―exp⁡(ik⋅x)exp⁡(ik′⋅x′)dkdk′=U2(2π)d∫pj2(k)C^Y(k)exp⁡(ik⋅x)dk∫δ(k+k′)exp⁡(ik′⋅x′)dk′(23)=U2(2π)d∫pj2(k)C^Y(k)exp⁡[ik⋅(x−x′)]dk
With (23), the integrand of (8) and (18) becomes
∬ u j ( x ) u j ( x ′ ) ― p ( x , t ; x ′ , t ′ ) d x d x ′ = = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) d k ∬ exp ⁡ [ i k ⋅ ( x − x ′ ) p ( x , t ; x ′ , t ′ ) ] d x d x ′ (24) = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U ( t − t ′ ) − k 2 D ( t − t ′ ) ] d k ∬ u j ( x ) u j x ′ ¯ p x , t ; x ′ , t ′ d x d x ′ = = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) d k ∬ exp ⁡ i k ⋅ x − x ′ p x , t ; x ′ , t ′ d x d x ′ (24) = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t − t ′ − k 2 D t − t ′ d k {:[∬ bar(u_(j)(x)u_(j)(x^(')))p(x,t;x^('),t^('))dxdx^(')=],[=(U^(2))/((2pi)^(d))intp_(j)^(2)(k) hat(C)_(Y)(k)dk∬exp[ik*(x-x^('))p(x,t;x^('),t^('))]dxdx^(')],[(24)=(U^(2))/((2pi)^(d))intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)U(t-t^('))-k^(2)D(t-t^('))]dk]:}\begin{align*} & \iint \overline{u_{j}(\mathbf{x}) u_{j}\left(\mathbf{x}^{\prime}\right)} p\left(\mathbf{x}, t ; \mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x} d \mathbf{x}^{\prime}= \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) d \mathbf{k} \iint \exp \left[i \mathbf{k} \cdot\left(\mathbf{x}-\mathbf{x}^{\prime}\right) p\left(\mathbf{x}, t ; \mathbf{x}^{\prime}, t^{\prime}\right)\right] d \mathbf{x} d \mathbf{x}^{\prime} \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t-t^{\prime}\right)\right] d \mathbf{k} \tag{24} \end{align*}∬uj(x)uj(x′)―p(x,t;x′,t′)dxdx′==U2(2π)d∫pj2(k)C^Y(k)dk∬exp⁡[ik⋅(x−x′)p(x,t;x′,t′)]dxdx′(24)=U2(2π)d∫pj2(k)C^Y(k)exp⁡[ik1U(t−t′)−k2D(t−t′)]dk
The passage to the last line in (24) was achieved through the relation
∬ exp ⁡ [ i k ⋅ ( x − x ′ ) ] p ( x , t ; x ′ , t ′ ) d x d x ′ = exp ⁡ [ i k 1 U ( t − t ′ ) − k 2 D ( t − t ′ ) ] ∬ exp ⁡ i k ⋅ x − x ′ p x , t ; x ′ , t ′ d x d x ′ = exp ⁡ i k 1 U t − t ′ − k 2 D t − t ′ ∬exp[ik*(x-x^('))]p(x,t;x^('),t^('))dxdx^(')=exp[ik_(1)U(t-t^('))-k^(2)D(t-t^('))]\iint \exp \left[i \mathbf{k} \cdot\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\right] p\left(\mathbf{x}, t ; \mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x} d \mathbf{x}^{\prime}=\exp \left[i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t-t^{\prime}\right)\right]∬exp⁡[ik⋅(x−x′)]p(x,t;x′,t′)dxdx′=exp⁡[ik1U(t−t′)−k2D(t−t′)]
for the characteristic function of the Gaussian increment ( X t ( 0 ) − X t ′ ( 0 ) ) X t ( 0 ) − X t ′ ( 0 ) (X_(t)^((0))-X_(t^('))^((0)))\left(\mathbf{X}_{t}^{(0)}-\mathbf{X}_{t^{\prime}}^{(0)}\right)(Xt(0)−Xt′(0)) of the un-perturbed process X j , t ( 0 ) = δ j , 1 U t + W j , t , j = 1 ⋯ d X j , t ( 0 ) = δ j , 1 U t + W j , t , j = 1 ⋯ d X_(j,t)^((0))=delta_(j,1)Ut+W_(j,t),j=1cdots dX_{j, t}^{(0)}=\delta_{j, 1} U t+W_{j, t}, j=1 \cdots dXj,t(0)=δj,1Ut+Wj,t,j=1⋯d (e.g., [18]).
Remark 3. Changing the sign convention of the Fourier transform,
u ^ j ( k ) = ∫ u j ( x ) exp ⁡ ( i k ⋅ x ) d x u ^ j ( k ) = ∫ u j ( x ) exp ⁡ ( i k ⋅ x ) d x hat(u)_(j)(k)=intu_(j)(x)exp(ik*x)dx\hat{u}_{j}(\mathbf{k})=\int u_{j}(\mathbf{x}) \exp (i \mathbf{k} \cdot \mathbf{x}) d \mathbf{x}u^j(k)=∫uj(x)exp⁡(ik⋅x)dx
changes the sign of the first term of the characteristic function without altering the sign of the second term, exp ⁡ [ − i k 1 U ( t − t ′ ) − k 2 D ( t − t ′ ) ] exp ⁡ − i k 1 U t − t ′ − k 2 D t − t ′ exp[-ik_(1)U(t-t^('))-k^(2)D(t-t^('))]\exp \left[-i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t-t^{\prime}\right)\right]exp⁡[−ik1U(t−t′)−k2D(t−t′)] (see detailed computation of characteristic function for Gaussian random variables of Papoulis and Pillai [18], Example 5-28, p. 153]. However, changing the sign of the imaginary exponent is equivalent to changing U U UUU into − U − U -U-U−U (this also changes the sign of p j p j p_(j)p_{j}pj (see Gelhar and Axness [14], Eq. (28a)), but (24) depends only on p j 2 p j 2 p_(j)^(2)p_{j}^{2}pj2 ). Since second moments and dispersion coefficients do not depend on the sign of the mean flow, the expression (26) gives the same result irrespective of the sign of the imaginary exponent.
Inserting (24) in (8) one obtains the explicit expression of the advective contribution to the ensemble dispersion coefficient,
δ { D i i e n s ( t ) } = D i i e n s ( t ) − D (25) = U 2 ( 2 π ) d ∫ 0 t d t ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U ( t − t ′ ) − k 2 D ( t − t ′ ) ] d k δ D i i e n s ( t ) = D i i e n s ( t ) − D (25) = U 2 ( 2 π ) d ∫ 0 t   d t ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t − t ′ − k 2 D t − t ′ d k {:[delta{D_(ii)^(ens)(t)}=D_(ii)^(ens)(t)-D],[(25)=(U^(2))/((2pi)^(d))int_(0)^(t)dt^(')intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)U(t-t^('))-k^(2)D(t-t^('))]dk]:}\begin{align*} \delta\left\{D_{i i}^{e n s}(t)\right\} & =D_{i i}^{e n s}(t)-D \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int_{0}^{t} d t^{\prime} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t-t^{\prime}\right)\right] d \mathbf{k} \tag{25} \end{align*}δ{Diiens(t)}=Diiens(t)−D(25)=U2(2π)d∫0tdt′∫pj2(k)C^Y(k)exp⁡[ik1U(t−t′)−k2D(t−t′)]dk
With the change of variable t ′ ′ = t − t ′ t ′ ′ = t − t ′ t^('')=t-t^(')t^{\prime \prime}=t-t^{\prime}t′′=t−t′, (25) becomes
δ { D i i e n s ( t ) } = − U 2 ( 2 π ) d ∫ t 0 d t ′ ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U t ′ ′ − k 2 D t ′ ′ ] d k (26) = U 2 ( 2 π ) d ∫ 0 t d t ′ ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U t ′ ′ − k 2 D t ′ ′ ] d k δ D i i e n s ( t ) = − U 2 ( 2 π ) d ∫ t 0   d t ′ ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t ′ ′ − k 2 D t ′ ′ d k (26) = U 2 ( 2 π ) d ∫ 0 t   d t ′ ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t ′ ′ − k 2 D t ′ ′ d k {:[delta{D_(ii)^(ens)(t)}=-(U^(2))/((2pi)^(d))int_(t)^(0)dt^('')intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)Ut^('')-k^(2)Dt^('')]dk],[(26)=(U^(2))/((2pi)^(d))int_(0)^(t)dt^('')intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)Ut^('')-k^(2)Dt^('')]dk]:}\begin{align*} \delta\left\{D_{i i}^{e n s}(t)\right\} & =-\frac{U^{2}}{(2 \pi)^{d}} \int_{t}^{0} d t^{\prime \prime} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U t^{\prime \prime}-k^{2} D t^{\prime \prime}\right] d \mathbf{k} \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int_{0}^{t} d t^{\prime \prime} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U t^{\prime \prime}-k^{2} D t^{\prime \prime}\right] d \mathbf{k} \tag{26} \end{align*}δ{Diiens(t)}=−U2(2π)d∫t0dt′′∫pj2(k)C^Y(k)exp⁡[ik1Ut′′−k2Dt′′]dk(26)=U2(2π)d∫0tdt′′∫pj2(k)C^Y(k)exp⁡[ik1Ut′′−k2Dt′′]dk
Remark 4.1. The time derivative of (26) is identical with the result of Dagan ([5], Eq. (3.14); [6], Eq. (17)), and Schwarze et al. ([20], Eq. (13)). Dagan [6] derived the result by using the changed sign in the homogeneity condition (Remark 2) and opposite Fourier transform convention (Remark 3).
Remark 4.2. The opposite Fourier transform convention and the same homogeneity condition as in (22) change i k 1 U t ′ ′ i k 1 U t ′ ′ ik_(1)Ut^('')i k_{1} U t^{\prime \prime}ik1Ut′′ into − i k 1 U t ′ ′ − i k 1 U t ′ ′ -ik_(1)Ut^('')-i k_{1} U t^{\prime \prime}−ik1Ut′′ and let − k 2 D t ′ ′ − k 2 D t ′ ′ -k^(2)Dt^('')-k^{2} D t^{\prime \prime}−k2Dt′′ unchanged (see Remark 3). This is the case in the derivation of the variance Σ j j Σ j j Sigma_(jj)\Sigma_{j j}Σjj by Fiori and Dagan [12] (Eq. (14)), which is exactly the integral of (25) from above, with i k 1 U t ′ ′ i k 1 U t ′ ′ ik_(1)Ut^('')i k_{1} U t^{\prime \prime}ik1Ut′′ replaced by − i k 1 U t ′ ′ − i k 1 U t ′ ′ -ik_(1)Ut^('')-i k_{1} U t^{\prime \prime}−ik1Ut′′. Dentz et al. [7] used the same conventions as Fiori and Dagan [12], (i.e. Eqs. (16) and (25) in [7]) and obtained the same result. Their result is given by Eq. (40) with g ~ 0 g ~ 0 tilde(g)_(0)\tilde{g}_{0}g~0 replaced by the characteristic function of the un-perturbed problem ([7], Eq. (25)), which is just the relation (26) from above, with exponent i k 1 U t ′ ′ i k 1 U t ′ ′ ik_(1)Ut^('')i k_{1} U t^{\prime \prime}ik1Ut′′ replaced by − i k 1 U t ′ ′ − i k 1 U t ′ ′ -ik_(1)Ut^('')-i k_{1} U t^{\prime \prime}−ik1Ut′′.
The corresponding memory coefficient is similarly obtained by inserting (24) into (18),
(27) M D i i e n s ( t , τ , 0 ) = U 2 ( 2 π ) d ∫ 0 τ d t ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U ( t − t ′ ) − k 2 D ( t − t ′ ) ] d k (27) M D i i e n s ( t , τ , 0 ) = U 2 ( 2 π ) d ∫ 0 τ   d t ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t − t ′ − k 2 D t − t ′ d k {:(27)MD_(ii)^(ens)(t","tau","0)=(U^(2))/((2pi)^(d))int_(0)^(tau)dt^(')intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)U(t-t^('))-k^(2)D(t-t^('))]dk:}\begin{equation*} M D_{i i}^{e n s}(t, \tau, 0)=\frac{U^{2}}{(2 \pi)^{d}} \int_{0}^{\tau} d t^{\prime} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t-t^{\prime}\right)\right] d \mathbf{k} \tag{27} \end{equation*}(27)MDiiens(t,τ,0)=U2(2π)d∫0τdt′∫pj2(k)C^Y(k)exp⁡[ik1U(t−t′)−k2D(t−t′)]dk
Making the change of variable t ′ ′ = t − t ′ t ′ ′ = t − t ′ t^('')=t-t^(')t^{\prime \prime}=t-t^{\prime}t′′=t−t′ in (27), we see that the ensemble memory coefficients can be computed, after a simple change of the integration limits in the time integral, with the explicit expression (26) of the advective contributions δ { D i i e n s } δ D i i e n s delta{D_(ii)^(ens)}\delta\left\{D_{i i}^{e n s}\right\}δ{Diiens} to the ensemble coefficients:
(28) M D i i e n s ( t , τ , 0 ) = U 2 ( 2 π ) d ∫ t − τ t d t ′ ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U t ′ ′ − k 2 D t ′ ′ ] d k (28) M D i i e n s ( t , τ , 0 ) = U 2 ( 2 π ) d ∫ t − τ t   d t ′ ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t ′ ′ − k 2 D t ′ ′ d k {:(28)MD_(ii)^(ens)(t","tau","0)=(U^(2))/((2pi)^(d))int_(t-tau)^(t)dt^('')intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)Ut^('')-k^(2)Dt^('')]dk:}\begin{equation*} M D_{i i}^{e n s}(t, \tau, 0)=\frac{U^{2}}{(2 \pi)^{d}} \int_{t-\tau}^{t} d t^{\prime \prime} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U t^{\prime \prime}-k^{2} D t^{\prime \prime}\right] d \mathbf{k} \tag{28} \end{equation*}(28)MDiiens(t,τ,0)=U2(2π)d∫t−τtdt′′∫pj2(k)C^Y(k)exp⁡[ik1Ut′′−k2Dt′′]dk
Remark 5. The memory coefficient (27) can also be obtained within the Eulerian approach of Dentz et al. [7] if in Eq. (18) from [8] the initial concentration ρ ^ ρ ^ hat(rho)\hat{\rho}ρ^ is replaced by the perturbation solution for point source at time τ τ tau\tauτ given in Eq. (26) from [7], i.e. ρ ^ ( k , τ ) = g ( k , τ ) ρ ^ ( k , τ ) = g ( k , τ ) hat(rho)(k,tau)=g(k,tau)\hat{\rho}(k, \tau)=g(k, \tau)ρ^(k,τ)=g(k,τ). Then Eq. (18) from [8] gives g ^ ( k , t − τ ) g ^ ( k , t − τ ) hat(g)(k,t-tau)\hat{g}(k, t-\tau)g^(k,t−τ). Further, inserting g ^ ( k , t − τ ) g ^ ( k , t − τ ) hat(g)(k,t-tau)\hat{g}(k, t-\tau)g^(k,t−τ) in Eq. (14) from [8] and "expanding the logarithm consistently up to second order in the fluctuations" (as described at the beginning of Sect. 3.2 of [8]), should reproduce the expression (27) of the memory coefficient M D ens ( t , τ , 0 ) M D ens  ( t , τ , 0 ) MD^("ens ")(t,tau,0)M D^{\text {ens }}(t, \tau, 0)MDens (t,τ,0).
To derive the explicit form of the effective coefficients (7) we first have to compute the center of mass coefficients (9). Similarly to (24), we obtain
∬ u j ( x ) u j ( x ′ ) ― p ( x , t ) p ( x ′ , t ′ ) d x d x ′ = = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) d k ∬ exp ⁡ [ i k ⋅ ( x − x ′ ) p ( x , t ) p ( x ′ , t ′ ) ] d x d x ′ (29) = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U ( t − t ′ ) − k 2 D ( t + t ′ ) ] d k ∬ u j ( x ) u j x ′ ¯ p ( x , t ) p x ′ , t ′ d x d x ′ = = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) d k ∬ exp ⁡ i k ⋅ x − x ′ p ( x , t ) p x ′ , t ′ d x d x ′ (29) = U 2 ( 2 π ) d ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t − t ′ − k 2 D t + t ′ d k {:[∬ bar(u_(j)(x)u_(j)(x^(')))p(x","t)p(x^('),t^('))dxdx^(')=],[=(U^(2))/((2pi)^(d))intp_(j)^(2)(k) hat(C)_(Y)(k)dk∬exp[ik*(x-x^('))p(x,t)p(x^('),t^('))]dxdx^(')],[(29)=(U^(2))/((2pi)^(d))intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)U(t-t^('))-k^(2)D(t+t^('))]dk]:}\begin{align*} & \iint \overline{u_{j}(\mathbf{x}) u_{j}\left(\mathbf{x}^{\prime}\right)} p(\mathbf{x}, t) p\left(\mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x} d \mathbf{x}^{\prime}= \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) d \mathbf{k} \iint \exp \left[i \mathbf{k} \cdot\left(\mathbf{x}-\mathbf{x}^{\prime}\right) p(\mathbf{x}, t) p\left(\mathbf{x}^{\prime}, t^{\prime}\right)\right] d \mathbf{x} d \mathbf{x}^{\prime} \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t+t^{\prime}\right)\right] d \mathbf{k} \tag{29} \end{align*}∬uj(x)uj(x′)―p(x,t)p(x′,t′)dxdx′==U2(2π)d∫pj2(k)C^Y(k)dk∬exp⁡[ik⋅(x−x′)p(x,t)p(x′,t′)]dxdx′(29)=U2(2π)d∫pj2(k)C^Y(k)exp⁡[ik1U(t−t′)−k2D(t+t′)]dk
The passage to the last line in (29) is now achieved through the relations
∬ exp ⁡ ( i k ⋅ x ) p ( x , t ) d x = exp ⁡ [ i k 1 U t − k 2 D t ] ∬ exp ⁡ ( − i k ⋅ x ′ ) p ( x ′ , t ′ ) d x ′ = exp ⁡ [ − i k 1 U t ′ − k 2 D t ′ ] ∬ exp ⁡ ( i k ⋅ x ) p ( x , t ) d x = exp ⁡ i k 1 U t − k 2 D t ∬ exp ⁡ − i k ⋅ x ′ p x ′ , t ′ d x ′ = exp ⁡ − i k 1 U t ′ − k 2 D t ′ {:[∬exp(ik*x)p(x","t)dx=exp[ik_(1)Ut-k^(2)Dt]],[∬exp(-ik*x^('))p(x^('),t^('))dx^(')=exp[-ik_(1)Ut^(')-k^(2)Dt^(')]]:}\begin{aligned} \iint \exp (i \mathbf{k} \cdot \mathbf{x}) p(\mathbf{x}, t) d \mathbf{x} & =\exp \left[i k_{1} U t-k^{2} D t\right] \\ \iint \exp \left(-i \mathbf{k} \cdot \mathbf{x}^{\prime}\right) p\left(\mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x}^{\prime} & =\exp \left[-i k_{1} U t^{\prime}-k^{2} D t^{\prime}\right] \end{aligned}∬exp⁡(ik⋅x)p(x,t)dx=exp⁡[ik1Ut−k2Dt]∬exp⁡(−ik⋅x′)p(x′,t′)dx′=exp⁡[−ik1Ut′−k2Dt′]
for the characteristic functions of the Gaussian un-perturbed process X j , t ( 0 ) = δ j , 1 U t + W j , t , j = 1 ⋯ d X j , t ( 0 ) = δ j , 1 U t + W j , t , j = 1 ⋯ d X_(j,t)^((0))=delta_(j,1)Ut+W_(j,t),j=1cdots dX_{j, t}^{(0)}=\delta_{j, 1} U t+W_{j, t}, j=1 \cdots dXj,t(0)=δj,1Ut+Wj,t,j=1⋯d [18], where we used the property pointed out in Remark 3. Inserting (29) into (9), we obtain the center of mass coefficient
(30) D i i c m ( t ) = U 2 ( 2 π ) d ∫ 0 t d t ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U ( t − t ′ ) − k 2 D ( t + t ′ ) ] d k (30) D i i c m ( t ) = U 2 ( 2 π ) d ∫ 0 t   d t ′ ∫ p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t − t ′ − k 2 D t + t ′ d k {:(30)D_(ii)^(cm)(t)=(U^(2))/((2pi)^(d))int_(0)^(t)dt^(')intp_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)U(t-t^('))-k^(2)D(t+t^('))]dk:}\begin{equation*} D_{i i}^{c m}(t)=\frac{U^{2}}{(2 \pi)^{d}} \int_{0}^{t} d t^{\prime} \int p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U\left(t-t^{\prime}\right)-k^{2} D\left(t+t^{\prime}\right)\right] d \mathbf{k} \tag{30} \end{equation*}(30)Diicm(t)=U2(2π)d∫0tdt′∫pj2(k)C^Y(k)exp⁡[ik1U(t−t′)−k2D(t+t′)]dk
With the change of variable t ′ ′ = t − t ′ t ′ ′ = t − t ′ t^('')=t-t^(')t^{\prime \prime}=t-t^{\prime}t′′=t−t′ (which implies t + t ′ = 2 t − t ′ ′ t + t ′ = 2 t − t ′ ′ t+t^(')=2t-t^('')t+t^{\prime}=2 t-t^{\prime \prime}t+t′=2t−t′′ ) the expression (30) of the center of mass coefficients becomes
D i i c m ( t ) = − U 2 ( 2 π ) d ∫ t 0 d t ′ ′ ∫ exp ⁡ ( − 2 k 2 D t ) p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U t ′ ′ + k 2 D t ′ ′ ] d k (31) = U 2 ( 2 π ) d ∫ 0 t d t ′ ′ ∫ exp ⁡ ( − 2 k 2 D t ) p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U t ′ ′ + k 2 D t ′ ′ ] d k D i i c m ( t ) = − U 2 ( 2 π ) d ∫ t 0   d t ′ ′ ∫ exp ⁡ − 2 k 2 D t p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t ′ ′ + k 2 D t ′ ′ d k (31) = U 2 ( 2 π ) d ∫ 0 t   d t ′ ′ ∫ exp ⁡ − 2 k 2 D t p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t ′ ′ + k 2 D t ′ ′ d k {:[D_(ii)^(cm)(t)=-(U^(2))/((2pi)^(d))int_(t)^(0)dt^('')int exp(-2k^(2)Dt)p_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)Ut^('')+k^(2)Dt^('')]dk],[(31)=(U^(2))/((2pi)^(d))int_(0)^(t)dt^('')int exp(-2k^(2)Dt)p_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)Ut^('')+k^(2)Dt^('')]dk]:}\begin{align*} D_{i i}^{c m}(t) & =-\frac{U^{2}}{(2 \pi)^{d}} \int_{t}^{0} d t^{\prime \prime} \int \exp \left(-2 k^{2} D t\right) p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U t^{\prime \prime}+k^{2} D t^{\prime \prime}\right] d \mathbf{k} \\ & =\frac{U^{2}}{(2 \pi)^{d}} \int_{0}^{t} d t^{\prime \prime} \int \exp \left(-2 k^{2} D t\right) p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U t^{\prime \prime}+k^{2} D t^{\prime \prime}\right] d \mathbf{k} \tag{31} \end{align*}Diicm(t)=−U2(2π)d∫t0dt′′∫exp⁡(−2k2Dt)pj2(k)C^Y(k)exp⁡[ik1Ut′′+k2Dt′′]dk(31)=U2(2π)d∫0tdt′′∫exp⁡(−2k2Dt)pj2(k)C^Y(k)exp⁡[ik1Ut′′+k2Dt′′]dk
The center of mass coefficients (30) and the ensemble contribution (26) determine the advective contributions to the effective coefficients according to (7):
(32) δ { D i i e f f ( t ) } = D i i e n s ( t ) − D i , i c m ( t ) − D = δ { D i i e n s ( t ) } − D i , i c m ( t ) (32) δ D i i e f f ( t ) = D i i e n s ( t ) − D i , i c m ( t ) − D = δ D i i e n s ( t ) − D i , i c m ( t ) {:(32)delta{D_(ii)^(eff)(t)}=D_(ii)^(ens)(t)-D_(i,i)^(cm)(t)-D=delta{D_(ii)^(ens)(t)}-D_(i,i)^(cm)(t):}\begin{equation*} \delta\left\{D_{i i}^{e f f}(t)\right\}=D_{i i}^{e n s}(t)-D_{i, i}^{c m}(t)-D=\delta\left\{D_{i i}^{e n s}(t)\right\}-D_{i, i}^{c m}(t) \tag{32} \end{equation*}(32)δ{Diieff(t)}=Diiens(t)−Di,icm(t)−D=δ{Diiens(t)}−Di,icm(t)
Remark 6. The ensemble contribution (26) and the center of mass coefficient (30) particularize to isotropic local diffusion coefficient the general relations M − M − M^(-)M^{-}M−and M + M + M^(+)M^{+}M+, respectively, derived by Dentz et al. ([7] Eq. (47)). We also note that (30) is the time derivative of the relation (15) of Fiori and Dagan [12], which, together with Remarks 4.1 and 4.2, proves the equivalence of Eulerian expressions of Dentz et al. [7, 8] and the Lagragian expressions derived by Fiori and Dagan [12] and by Vanderborght [26] for the ensemble and effective dispersion coefficients, respectively, variances. Nevertheless, the Eulerian approach has the important advantage of providing first-order approximations of the concentration fields (e.g., Eq. (26) in [7], on which the derivation of the dispersion coefficients is in actually based).
The memory coefficient of the center of mass is readily obtained, similarly to (28), by a change of integration limits in the expression (31) of the center of mass dispersion coefficients:
(33) M D i i c m ( t , τ , 0 ) = U 2 ( 2 π ) d ∫ t − τ t d t ′ ′ ∫ exp ⁡ ( − 2 k 2 D t ) p j 2 ( k ) C ^ Y ( k ) exp ⁡ [ i k 1 U t ′ ′ + k 2 D t ′ ′ ] d k (33) M D i i c m ( t , τ , 0 ) = U 2 ( 2 π ) d ∫ t − τ t   d t ′ ′ ∫ exp ⁡ − 2 k 2 D t p j 2 ( k ) C ^ Y ( k ) exp ⁡ i k 1 U t ′ ′ + k 2 D t ′ ′ d k {:(33)MD_(ii)^(cm)(t","tau","0)=(U^(2))/((2pi)^(d))int_(t-tau)^(t)dt^('')int exp(-2k^(2)Dt)p_(j)^(2)(k) hat(C)_(Y)(k)exp[ik_(1)Ut^('')+k^(2)Dt^('')]dk:}\begin{equation*} M D_{i i}^{c m}(t, \tau, 0)=\frac{U^{2}}{(2 \pi)^{d}} \int_{t-\tau}^{t} d t^{\prime \prime} \int \exp \left(-2 k^{2} D t\right) p_{j}^{2}(\mathbf{k}) \hat{C}_{Y}(\mathbf{k}) \exp \left[i k_{1} U t^{\prime \prime}+k^{2} D t^{\prime \prime}\right] d \mathbf{k} \tag{33} \end{equation*}(33)MDiicm(t,τ,0)=U2(2π)d∫t−τtdt′′∫exp⁡(−2k2Dt)pj2(k)C^Y(k)exp⁡[ik1Ut′′+k2Dt′′]dk
Finally, the memory coefficient of the effective dispersion is obtained from the relations (32), (31), and (28),
(34) M D i i e f f ( t ) = M D i i e n s ( t ) − M D i , i c m ( t ) . (34) M D i i e f f ( t ) = M D i i e n s ( t ) − M D i , i c m ( t ) . {:(34)MD_(ii)^(eff)(t)=MD_(ii)^(ens)(t)-MD_(i,i)^(cm)(t).:}\begin{equation*} M D_{i i}^{e f f}(t)=M D_{i i}^{e n s}(t)-M D_{i, i}^{c m}(t) . \tag{34} \end{equation*}(34)MDiieff(t)=MDiiens(t)−MDi,icm(t).
Remark 7. The memory coefficients (28), (33), and (34) can be computed with the explicit expressions of the dispersion coefficients derived in either Eulerian frames [[1,7,8] or by Lagrangian approaches [10,12,26]. To achieve that, we only have to replace the lower limit by t − τ t − τ t-taut-\taut−τ in the time integrals from these expressions.

4 Explicit first-order results for power law correlated log-K fields

The first-order results for log-K fields with finite correlation lengths can be further used to develop approaches for a superposition of scales [9, 19]. To proceed, let us consider the particular case of an isotropic power-law covariance of Y = ln ⁡ K Y = ln ⁡ K Y=ln KY=\ln KY=ln⁡K,
(35) C Y P L ( x ) = σ Y 2 z − β , where z = ( 1 + | x | 2 L 2 ) 1 2 , | x | 2 = ∑ j = 1 d x j 2 , and 0 ≤ β ≤ 2 (35) C Y P L ( x ) = σ Y 2 z − β ,  where  z = 1 + | x | 2 L 2 1 2 , | x | 2 = ∑ j = 1 d   x j 2 ,  and  0 ≤ β ≤ 2 {:(35)C_(Y)^(PL)(x)=sigma_(Y)^(2)z^(-beta)","" where "z=(1+(|x|^(2))/(L^(2)))^((1)/(2))","|x|^(2)=sum_(j=1)^(d)x_(j)^(2)","" and "0 <= beta <= 2:}\begin{equation*} C_{Y}^{P L}(\mathbf{x})=\sigma_{Y}^{2} z^{-\beta}, \text { where } z=\left(1+\frac{|\mathbf{x}|^{2}}{L^{2}}\right)^{\frac{1}{2}},|\mathbf{x}|^{2}=\sum_{j=1}^{d} x_{j}^{2}, \text { and } 0 \leq \beta \leq 2 \tag{35} \end{equation*}(35)CYPL(x)=σY2z−β, where z=(1+|x|2L2)12,|x|2=∑j=1dxj2, and 0≤β≤2
where L L LLL is some length scale associated with the spatial extension of the system. Power-law covariances can be constructed by integrating covariances with continuously increasing finite integral scales λ λ lambda\lambdaλ ([15], p. 370 , expression 3.478 (1.)) by the relation
(36) ∫ 0 ∞ λ β − 1 exp ⁡ ( − μ λ p ) d λ = 1 p μ − β p Γ ( β p , ) (36) ∫ 0 ∞   λ β − 1 exp ⁡ − μ λ p d λ = 1 p μ − β p Γ β p , {:(36)int_(0)^(oo)lambda^(beta-1)exp(-mulambda^(p))d lambda=(1)/(p)mu^(-(beta )/(p))Gamma((beta )/(p),):}\begin{equation*} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left(-\mu \lambda^{p}\right) d \lambda=\frac{1}{p} \mu^{-\frac{\beta}{p}} \Gamma\left(\frac{\beta}{p},\right) \tag{36} \end{equation*}(36)∫0∞λβ−1exp⁡(−μλp)dλ=1pμ−βpΓ(βp,)
where μ > 0 μ > 0 mu > 0\mu>0μ>0 and β > 0 β > 0 beta > 0\beta>0β>0. For μ = z 2 μ = z 2 mu=z^(2)\mu=z^{2}μ=z2 and p = 2 p = 2 p=2p=2p=2, (36) yields a superposition of Gaussian covariances:
C Y P L ( x ) / σ Y 2 = z − β = 2 Γ ( β 2 ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − z 2 λ 2 ) d λ = 2 Γ ( β 2 ) ∫ 0 ∞ λ β − 1 exp ⁡ [ − ( 1 + | x | 2 L 2 ) λ 2 ] d λ = 2 Γ ( β 2 ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − λ 2 ) exp ⁡ ( − | x | 2 L 2 λ 2 ) d λ = 2 Γ ( β 2 ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − λ 2 ) exp ⁡ ( − 1 2 | x | 2 1 2 L 2 λ 2 ) d λ (37) = ∫ 0 ∞ d Λ ( λ ) exp ⁡ ( − 1 2 | x | 2 l ( λ ) 2 ) C Y P L ( x ) / σ Y 2 = z − β = 2 Γ β 2 ∫ 0 ∞   λ β − 1 exp ⁡ − z 2 λ 2 d λ = 2 Γ β 2 ∫ 0 ∞   λ β − 1 exp ⁡ − 1 + | x | 2 L 2 λ 2 d λ = 2 Γ β 2 ∫ 0 ∞   λ β − 1 exp ⁡ − λ 2 exp ⁡ − | x | 2 L 2 λ 2 d λ = 2 Γ β 2 ∫ 0 ∞   λ β − 1 exp ⁡ − λ 2 exp ⁡ − 1 2 | x | 2 1 2 L 2 λ 2 d λ (37) = ∫ 0 ∞   d Λ ( λ ) exp ⁡ − 1 2 | x | 2 l ( λ ) 2 {:[C_(Y)^(PL)(x)//sigma_(Y)^(2)=z^(-beta)],[=(2)/(Gamma((beta)/(2)))int_(0)^(oo)lambda^(beta-1)exp(-z^(2)lambda^(2))d lambda],[=(2)/(Gamma((beta)/(2)))int_(0)^(oo)lambda^(beta-1)exp[-(1+(|x|^(2))/(L^(2)))lambda^(2)]d lambda],[=(2)/(Gamma((beta)/(2)))int_(0)^(oo)lambda^(beta-1)exp(-lambda^(2))exp(-(|x|^(2))/(L^(2))lambda^(2))d lambda],[=(2)/(Gamma((beta)/(2)))int_(0)^(oo)lambda^(beta-1)exp(-lambda^(2))exp(-(1)/(2)(|x|^(2))/((1)/(2)(L^(2))/(lambda^(2))))d lambda],[(37)=int_(0)^(oo)d Lambda(lambda)exp(-(1)/(2)(|x|^(2))/(l(lambda)^(2)))]:}\begin{align*} C_{Y}^{P L}(\mathbf{x}) / \sigma_{Y}^{2} & =z^{-\beta} \\ & =\frac{2}{\Gamma\left(\frac{\beta}{2}\right)} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left(-z^{2} \lambda^{2}\right) d \lambda \\ & =\frac{2}{\Gamma\left(\frac{\beta}{2}\right)} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left[-\left(1+\frac{|\mathbf{x}|^{2}}{L^{2}}\right) \lambda^{2}\right] d \lambda \\ & =\frac{2}{\Gamma\left(\frac{\beta}{2}\right)} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left(-\lambda^{2}\right) \exp \left(-\frac{|\mathbf{x}|^{2}}{L^{2}} \lambda^{2}\right) d \lambda \\ & =\frac{2}{\Gamma\left(\frac{\beta}{2}\right)} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left(-\lambda^{2}\right) \exp \left(-\frac{1}{2} \frac{|\mathbf{x}|^{2}}{\frac{1}{2} \frac{L^{2}}{\lambda^{2}}}\right) d \lambda \\ & =\int_{0}^{\infty} d \Lambda(\lambda) \exp \left(-\frac{1}{2} \frac{|\mathbf{x}|^{2}}{l(\lambda)^{2}}\right) \tag{37} \end{align*}CYPL(x)/σY2=z−β=2Γ(β2)∫0∞λβ−1exp⁡(−z2λ2)dλ=2Γ(β2)∫0∞λβ−1exp⁡[−(1+|x|2L2)λ2]dλ=2Γ(β2)∫0∞λβ−1exp⁡(−λ2)exp⁡(−|x|2L2λ2)dλ=2Γ(β2)∫0∞λβ−1exp⁡(−λ2)exp⁡(−12|x|212L2λ2)dλ(37)=∫0∞dΛ(λ)exp⁡(−12|x|2l(λ)2)
where
(38) ∫ 0 ∞ d Λ ( λ ) = 2 Γ ( β 2 ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − λ 2 ) d λ and l ( λ ) = L λ 2 (38) ∫ 0 ∞   d Λ ( λ ) = 2 Γ β 2 ∫ 0 ∞   λ β − 1 exp ⁡ − λ 2 d λ  and  l ( λ ) = L λ 2 {:(38)int_(0)^(oo)d Lambda(lambda)=(2)/(Gamma((beta)/(2)))int_(0)^(oo)lambda^(beta-1)exp(-lambda^(2))d lambdaquad" and "quad l(lambda)=(L)/(lambdasqrt2):}\begin{equation*} \int_{0}^{\infty} d \Lambda(\lambda)=\frac{2}{\Gamma\left(\frac{\beta}{2}\right)} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left(-\lambda^{2}\right) d \lambda \quad \text { and } \quad l(\lambda)=\frac{L}{\lambda \sqrt{2}} \tag{38} \end{equation*}(38)∫0∞dΛ(λ)=2Γ(β2)∫0∞λβ−1exp⁡(−λ2)dλ and l(λ)=Lλ2
The integrand in (37) has the from of the Gaussian C Y C Y C_(Y)C_{Y}CY considered by Dentz et al. [7,8].
Redefining z z zzz in (35) as z = ( 1 + | x | / L ) z = ( 1 + | x | / L ) z=(1+|x|//L)z=(1+|\mathbf{x}| / L)z=(1+|x|/L) and choosing μ = z μ = z mu=z\mu=zμ=z and p = 1 p = 1 p=1p=1p=1, one constructs a power-law covariance via (36) as a superposition of exponential covariances:
C Y P L ( x ) / σ Y 2 = z − β = 1 Γ ( β ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − z λ ) d λ = 1 Γ ( β ) ∫ 0 ∞ λ β − 1 exp ⁡ [ − ( 1 + | x | L ) λ ] d λ = 1 Γ ( β ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − λ ) exp ⁡ ( − | x | L λ ) d λ (39) = ∫ 0 ∞ d Λ ( λ ) exp ⁡ ( − | x | l ( λ ) ) C Y P L ( x ) / σ Y 2 = z − β = 1 Γ ( β ) ∫ 0 ∞   λ β − 1 exp ⁡ ( − z λ ) d λ = 1 Γ ( β ) ∫ 0 ∞   λ β − 1 exp ⁡ − 1 + | x | L λ d λ = 1 Γ ( β ) ∫ 0 ∞   λ β − 1 exp ⁡ ( − λ ) exp ⁡ − | x | L λ d λ (39) = ∫ 0 ∞   d Λ ( λ ) exp ⁡ − | x | l ( λ ) {:[C_(Y)^(PL)(x)//sigma_(Y)^(2)=z^(-beta)],[=(1)/(Gamma(beta))int_(0)^(oo)lambda^(beta-1)exp(-z lambda)d lambda],[=(1)/(Gamma(beta))int_(0)^(oo)lambda^(beta-1)exp[-(1+(|x|)/(L))lambda]d lambda],[=(1)/(Gamma(beta))int_(0)^(oo)lambda^(beta-1)exp(-lambda)exp(-(|x|)/(L)lambda)d lambda],[(39)=int_(0)^(oo)d Lambda(lambda)exp(-(|x|)/(l(lambda)))]:}\begin{align*} C_{Y}^{P L}(\mathbf{x}) / \sigma_{Y}^{2} & =z^{-\beta} \\ & =\frac{1}{\Gamma(\beta)} \int_{0}^{\infty} \lambda^{\beta-1} \exp (-z \lambda) d \lambda \\ & =\frac{1}{\Gamma(\beta)} \int_{0}^{\infty} \lambda^{\beta-1} \exp \left[-\left(1+\frac{|\mathbf{x}|}{L}\right) \lambda\right] d \lambda \\ & =\frac{1}{\Gamma(\beta)} \int_{0}^{\infty} \lambda^{\beta-1} \exp (-\lambda) \exp \left(-\frac{|\mathbf{x}|}{L} \lambda\right) d \lambda \\ & =\int_{0}^{\infty} d \Lambda(\lambda) \exp \left(-\frac{|\mathbf{x}|}{l(\lambda)}\right) \tag{39} \end{align*}CYPL(x)/σY2=z−β=1Γ(β)∫0∞λβ−1exp⁡(−zλ)dλ=1Γ(β)∫0∞λβ−1exp⁡[−(1+|x|L)λ]dλ=1Γ(β)∫0∞λβ−1exp⁡(−λ)exp⁡(−|x|Lλ)dλ(39)=∫0∞dΛ(λ)exp⁡(−|x|l(λ))
where
(40) ∫ 0 ∞ d Λ ( λ ) = 1 Γ ( β ) ∫ 0 ∞ λ β − 1 exp ⁡ ( − λ ) d λ and l ( λ ) = L λ (40) ∫ 0 ∞   d Λ ( λ ) = 1 Γ ( β ) ∫ 0 ∞   λ β − 1 exp ⁡ ( − λ ) d λ  and  l ( λ ) = L λ {:(40)int_(0)^(oo)d Lambda(lambda)=(1)/(Gamma(beta))int_(0)^(oo)lambda^(beta-1)exp(-lambda)d lambdaquad" and "quad l(lambda)=(L)/( lambda):}\begin{equation*} \int_{0}^{\infty} d \Lambda(\lambda)=\frac{1}{\Gamma(\beta)} \int_{0}^{\infty} \lambda^{\beta-1} \exp (-\lambda) d \lambda \quad \text { and } \quad l(\lambda)=\frac{L}{\lambda} \tag{40} \end{equation*}(40)∫0∞dΛ(λ)=1Γ(β)∫0∞λβ−1exp⁡(−λ)dλ and l(λ)=Lλ
The integrand in (39) is the isotropic version of the exponential C Y C Y C_(Y)C_{Y}CY considered by Fiori [10] and Vanderborght [26].
For β = 1 β = 1 beta=1\beta=1β=1 in (37) (or in (39)) one obtains the covariance of the " 1 / z 1 / z 1//z1 / z1/z noise" as an integral of Gaussian (or exponential) covariances with continuous scale parameter λ λ lambda\lambdaλ. Further, by defining z = | x | / L z = | x | / L z=|x|//Lz= |\mathbf{x}| / Lz=|x|/L in (39) one obtains d Λ ( λ ) = d λ d Λ ( λ ) = d λ d Lambda(lambda)=d lambdad \Lambda(\lambda)=d \lambdadΛ(λ)=dλ and (37), respectively (39), gives the covariance of the 1 / z 1 / z 1//z1 / z1/z noise as a simple integral of Gaussian, respectively exponential, modes. This corresponds to the limit of
infinite superposition of log-K fields with finite correlation scales, like those considered in the numerical investigations presented in the following.
By the linearity of the superposition relations (37), respectively (39), we have the general relations
(41) C Y P L ( x ) = ∫ 0 ∞ d Λ ( λ ) C Y ( x , λ ) , C Y P L ^ ( k ) = ∫ 0 ∞ d Λ ( λ ) C Y ^ ( k , λ ) (41) C Y P L ( x ) = ∫ 0 ∞   d Λ ( λ ) C Y ( x , λ ) , C Y P L ^ ( k ) = ∫ 0 ∞   d Λ ( λ ) C Y ^ ( k , λ ) {:(41)C_(Y)^(PL)(x)=int_(0)^(oo)d Lambda(lambda)C_(Y)(x","lambda)","quad widehat(C_(Y)^(PL))(k)=int_(0)^(oo)d Lambda(lambda) widehat(C_(Y))(k","lambda):}\begin{equation*} C_{Y}^{P L}(\mathbf{x})=\int_{0}^{\infty} d \Lambda(\lambda) C_{Y}(\mathbf{x}, \lambda), \quad \widehat{C_{Y}^{P L}}(\mathbf{k})=\int_{0}^{\infty} d \Lambda(\lambda) \widehat{C_{Y}}(\mathbf{k}, \lambda) \tag{41} \end{equation*}(41)CYPL(x)=∫0∞dΛ(λ)CY(x,λ),CYPL^(k)=∫0∞dΛ(λ)CY^(k,λ)
By inserting (41) in (24) and (5), we finally obtain the superposition relations
(42) δ { D i i P L , e n s } ( t ) = ∫ 0 ∞ d Λ ( λ ) { D i i e n s } ( t , λ ) , δ { D i i P L , e f f } ( t ) = ∫ 0 ∞ d Λ ( λ ) { D i i e f f } ( t , λ ) (42) δ D i i P L , e n s ( t ) = ∫ 0 ∞   d Λ ( λ ) D i i e n s ( t , λ ) , δ D i i P L , e f f ( t ) = ∫ 0 ∞   d Λ ( λ ) D i i e f f ( t , λ ) {:(42)delta{D_(ii)^(PL,ens)}(t)=int_(0)^(oo)d Lambda(lambda){D_(ii)^(ens)}(t","lambda)","quad delta{D_(ii)^(PL,eff)}(t)=int_(0)^(oo)d Lambda(lambda){D_(ii)^(eff)}(t","lambda):}\begin{equation*} \delta\left\{D_{i i}^{P L, e n s}\right\}(t)=\int_{0}^{\infty} d \Lambda(\lambda)\left\{D_{i i}^{e n s}\right\}(t, \lambda), \quad \delta\left\{D_{i i}^{P L, e f f}\right\}(t)=\int_{0}^{\infty} d \Lambda(\lambda)\left\{D_{i i}^{e f f}\right\}(t, \lambda) \tag{42} \end{equation*}(42)δ{DiiPL,ens}(t)=∫0∞dΛ(λ){Diiens}(t,λ),δ{DiiPL,eff}(t)=∫0∞dΛ(λ){Diieff}(t,λ)
According to (42), the memory coefficients (28), (33), and (34) are given by integrals with respect to d Λ ( λ ) d Λ ( λ ) d Lambda(lambda)d \Lambda(\lambda)dΛ(λ) of those for given λ λ lambda\lambdaλ, computed in the previous Section.

5 Numerical simulations of transport in aquifers with multiscale hydraulic conductivity

To investigate numerically the behavior of several transport observables we computed an ensemble of 256 GRW simulations of isotropic diffusion, with coefficient D = 0.01 m / d D = 0.01 m / d D=0.01m//dD=0.01 \mathrm{~m} / \mathrm{d}D=0.01 m/d which is a typical value for local dispersion coefficient, in multiscale velocity fields. The latter, generated with the Kraichnan method, are first-order approximations of the flow equations for a ln ⁡ K ln ⁡ K ln K\ln \mathrm{K}ln⁡K field consisting of a superposition of statistically independent ln ⁡ K ln ⁡ K ln K\ln \mathrm{K}ln⁡K with isotropic exponential correlation and constant variance 0.1 . The integral scales (defined as integrals of the correlation function, not divided by the variance) of the ln ⁡ K ln ⁡ K ln K\ln \mathrm{K}ln⁡K fields were increased with a constant step of 1 m , so that the superposition of 2 fields has the variance 0.2 and integral scale λ = 3 m λ = 3 m lambda=3m\lambda=3 \mathrm{~m}λ=3 m and that of seven fields (the largest one investigated here) has the variance 0.7 and λ = 28 m λ = 28 m lambda=28m\lambda=28 \mathrm{~m}λ=28 m. In our GRW numerical simulations it is more convenient to estimate dispersion coefficients for processes starting from instantaneous point injections by the mean half-slope, Σ i i ( t ) / ( 2 t ) , S i i ( t ) / ( 2 t ) , R i i ( t ) / ( 2 t ) Σ i i ( t ) / ( 2 t ) , S i i ( t ) / ( 2 t ) , R i i ( t ) / ( 2 t ) Sigma_(ii)(t)//(2t),S_(ii)(t)//(2t),R_(ii)(t)//(2t)\Sigma_{i i}(t) /(2 t), S_{i i}(t) /(2 t), R_{i i}(t) /(2 t)Σii(t)/(2t),Sii(t)/(2t),Rii(t)/(2t). The GRW algorithm is presented in detail by Vamoş et al. [25] and its implementation for transport in saturated aquifers is described in reference [22]. The effective, ensemble, and center of mass coefficients and the corresponding one-step memory dispersion coefficients ( M i i ( t ) / ( 2 t ) M i i ( t ) / ( 2 t ) (M_(ii)(t)//(2t):}\left(M_{i i}(t) /(2 t)\right.(Mii(t)/(2t), where M i i ( t ) = M i i ( t , t − Δ , 0 ) M i i ( t ) = M i i ( t , t − Δ , 0 ) M_(ii)(t)=M_(ii)(t,t-Delta,0)M_{i i}(t)=M_{i i}(t, t-\Delta, 0)Mii(t)=Mii(t,t−Δ,0) ), normalized by the local dispersion coefficient, are presented in Figs. 1-12. The extinction of the memory coefficients (Figs. 7-12) coincides with the approach to a normal diffusive behavior, in agreement with (15) and (17).
Figures 13 and 14 show the cumulated memory coefficients C M ( t ) / ( 2 t ) C M ( t ) / ( 2 t ) CM(t)//(2t)C M(t) /(2 t)CM(t)/(2t), with C M C M CMC MCM given by (16), compared with the corresponding ensemble, effective and center of mass coefficients, in case of the seven-scales velocity field. The cumulated memory coefficients are close to the corresponding dispersion coefficients, the small difference between the two coefficients being just the half slope of the average of the first term in (11). The latter is a dispersion coefficient describing the displacement of solute particles over one step Δ Δ Delta\DeltaΔ (which in our case is also the simulation time step). In case of ensemble dispersion Σ i i Σ i i Sigma_(ii)\Sigma_{i i}Σii, we can see that the average of (11) is a discrete Taylor formula approximating the ensemble coefficient (5). Thus, ensemble coefficients are mainly determined by the increment correlations or by correlations of the Lagrangian velocity at different moments. This makes the difference from genuine diffusion processes (Brownian motions) where the increments are uncorrelated and the one-step dispersion defines the constant diffusion coefficient.
Figures 15 and 16 show the ensemble dispersion coefficients for one time-step increments of the process X i , t e n s X i , t e n s X_(i,t)^(ens)X_{i, t}^{e n s}Xi,tens, normalized by the local dispersion coefficient as well. According to (5), Σ i i ( t , t − Δ ) Σ i i ( t , t − Δ ) Sigma_(ii)(t,t-Delta)\Sigma_{i i}(t, t-\Delta)Σii(t,t−Δ) are given by the sum between the local dispersion coefficient and the integral of the Lagrangian velocity over a time step Δ Δ Delta\DeltaΔ and they increase with the integral scale.
The asymmetry of the solute plume at 1000 days, presented in Fig. 17, is more pronounced than that recorded at 2000 days (Fig. 18). The asymmetry of the plume indicates that many solute particles move slower than the center of mass and remain behind the plume's front. The explanation is, very likely, the persistence of negative velocity fluctuations, rather than the increase of fluctuations' magnitude (note that the variance of the ln ⁡ K ln ⁡ K ln K\ln \mathrm{K}ln⁡K field is only 0.7 ). This agrees with the anomalous behavior of the dispersion coefficients over the first 1000 days (Figs. 1 and 3) and with expected approach to a Fickian regime in the limit of very large times, ensured by the finiteness of the integral scale [23].
The breakthrough curves (Figs. 19-20) are also strongly asymmetric at early stages and become more and more symmetrical as the process approaches the Fickian regime. Ensembles of simulated breakthrough curves are available and can be used to compute means and variances of the travel times (Figs. 21-22). Complete probability distributions of breakthrough concentration or of the arrival times at reference planes can be readily obtained from the ensemble of 256 transport simulations.
The concentration at the center of mass, integrated over the transverse dimension of the plume decreases, and the corresponding coefficient of variation increases, with the increase of the integral scale (Figs. 23-24). One remarks a strong non-Gaussian behavior in seven-scales case (compare black and red lines). More information on effect of multiscale fields provide the spatial profiles of the cross-section mean concentration and representations based on them presented in Figs. 25-26. Here one observes small asymmetry of the concentration profiles, and mainly that of their standard deviation, for multiscale fields in the pre-asymptotic regime, which diminishes at larger travel distances. More evident is the asymmetry shown by Figs. 27-28. There we represented the asymmetry of the longitudinal dispersion coefficient computed as follows. First we decomposed the effective coefficient in single realizations as a sum of cross-section contributions
D 11 e f f ( t ) = ∬ ( x 1 − ⟨ x 1 ⟩ ) 2 p ( x 1 , x 2 , t ) d x 1 d x 2 = ∫ d x 1 ∫ ( x 1 − ⟨ x 1 ⟩ ) 2 p ( x 1 , x 2 , t ) d x 2 (43) = ∫ d 11 e f f ( x 1 , t ) d x 1 D 11 e f f ( t ) = ∬ x 1 − x 1 2 p x 1 , x 2 , t d x 1 d x 2 = ∫ d x 1 ∫ x 1 − x 1 2 p x 1 , x 2 , t d x 2 (43) = ∫ d 11 e f f x 1 , t d x 1 {:[D_(11)^(eff)(t)=∬(x_(1)-(:x_(1):))^(2)p(x_(1),x_(2),t)dx_(1)dx_(2)],[=int dx_(1)int(x_(1)-(:x_(1):))^(2)p(x_(1),x_(2),t)dx_(2)],[(43)=intd_(11)^(eff)(x_(1),t)dx_(1)]:}\begin{gather*} D_{11}^{e f f}(t)=\iint\left(x_{1}-\left\langle x_{1}\right\rangle\right)^{2} p\left(x_{1}, x_{2}, t\right) d x_{1} d x_{2} \\ =\int d x_{1} \int\left(x_{1}-\left\langle x_{1}\right\rangle\right)^{2} p\left(x_{1}, x_{2}, t\right) d x_{2} \\ =\int d_{11}^{e f f}\left(x_{1}, t\right) d x_{1} \tag{43} \end{gather*}D11eff(t)=∬(x1−⟨x1⟩)2p(x1,x2,t)dx1dx2=∫dx1∫(x1−⟨x1⟩)2p(x1,x2,t)dx2(43)=∫d11eff(x1,t)dx1
The ensemble averages d 11 e f f ( x 1 , t ) ― d 11 e f f x 1 , t ¯ bar(d_(11)^(eff)(x_(1),t))\overline{d_{11}^{e f f}\left(x_{1}, t\right)}d11eff(x1,t)― and the corresponding standard deviations are presented in Fig. 27 and 28, respectively. Here, again, one can see that asymmetry persists in pre-asymptotic regime and increases with the integral scale of the velocity field. Further, we quantified the asymmetry by
(44) ∫ x 1 < ⟨ x 1 ⟩ d 11 e f f ― ( x 1 , t ) d x 1 / ∫ x 1 > ⟨ x 1 ⟩ d 11 e f f ― ( x 1 , t ) d x 1 (44) ∫ x 1 < x 1   d 11 e f f ¯ x 1 , t d x 1 / ∫ x 1 > x 1   d 11 e f f ¯ x 1 , t d x 1 {:(44)int_(x_(1) < (:x_(1):)) bar(d_(11)^(eff))(x_(1),t)dx_(1)//int_(x_(1) > (:x_(1):)) bar(d_(11)^(eff))(x_(1),t)dx_(1):}\begin{equation*} \int_{x_{1}<\left\langle x_{1}\right\rangle} \overline{d_{11}^{e f f}}\left(x_{1}, t\right) d x_{1} / \int_{x_{1}>\left\langle x_{1}\right\rangle} \overline{d_{11}^{e f f}}\left(x_{1}, t\right) d x_{1} \tag{44} \end{equation*}(44)∫x1<⟨x1⟩d11eff―(x1,t)dx1/∫x1>⟨x1⟩d11eff―(x1,t)dx1
The asymmetry coefficient (44) is presented in Fig. 29 and the corresponding effective coefficient (43) (averaged over realizations) in Fig. 30. Fig. 29 supplies another evidence of the relation between integral scale and asymmetry, which manifests in pre-asymptotic regime. (The large positive asymmetry at 2000 days in seven-scales case is very likely an artifact due to numerical trapping of particles in twodimensional GRW simulations [22]).
To further quantify the influence of multiscale velocity, we fitted dispersion and memory coefficients with power-law functions. The corresponding lines in log-log plots from Figs. 31-32 fit the behavior of the dispersion coefficients, variances, and memory terms in time windows proportional to the integral scale, of 15 , 25 , 35 , 55 15 , 25 , 35 , 55 15,25,35,5515,25,35,5515,25,35,55, and 300 days for single-, two-, three-, five-, and seven-scale fields, respectively. The increase of time window with integral scale is consistent with theoretical and numerical results obtained in the past [9, 3]. Figures 33 and 34 show that ensemble and effective coefficients and variances have different time behavior, with larger slope of the effective quantities. This seems to be at variance with the theoretical result for power-law correlated fields obtained by Fiori [11].
There is a long standing attempt to introduce power law correlations and anomalous dispersion in modeling groundwater transport by a rather abstract mathematical approach [4, 2, 9, 11, 3]. Nevertheless, a superposition of a finite number of scales seems to be a more natural assumption in many contamination scenarios for evolving scale groundwater systems [13, 17]. Numerical simulations of
multiscale transport may be a bridge to theoretical modeling by fractional Brownian motion ( fBm ) and fractional Gaussian noise (fGn) models. For instance, the super-diffusive behavior shown in Figs. 3134 can be modeled as a fBm and the fGn generating it can be readily obtained by two time derivatives, which reduce the power-law exponents by 2 [ 23 , 21 , 16 ] 2 [ 23 , 21 , 16 ] 2[23,21,16]2[23,21,16]2[23,21,16].
As we discussed in the previous Section, power-law covariances can be constructed by integrating covariances with continuously increasing integral scales [15]. An infinite sum of equal-weight exponential correlations results in an 1 / x 1 / x 1//x1 / x1/x-type correlation (according to [15], p. 370, expression 3.478 (1.); see also the discussion after Eq. (40) in Section 4). Such correlation leads to a ( t ln ⁡ t − t t ln ⁡ t − t t ln t-tt \ln t-ttln⁡t−t ) behavior of the ensemble variance Σ i i Σ i i Sigma_(ii)\Sigma_{i i}Σii [21]. Since this would be the theoretical infinite limit of our simulations, we also fitted dispersion coefficients and variances with ln ⁡ t ln ⁡ t ln t\ln tln⁡t and t ln ⁡ t − t t ln ⁡ t − t t ln t-tt \ln t-ttln⁡t−t respectively (Figures 35-36). At a first sight the fitting results are also acceptable. Nevertheless, a more detailed analysis is required to decide which time-behavior, the power-law one from Figs. 33-34 or the lnt-type behavior from Figs. 35-36 provides the better fit with the results of our numerical experiment. An independent selection criterium would be the numerical evaluation of dispersion and memory terms derived in Section 3 for an 1 / x 1 / x 1//x1 / x1/x-correlated log-K field, using the integral representation (39-40) derived in Section 4.

References

[1] Attinger, S., Dentz, M., H. Kinzelbach, and W. Kinzelbach (1999), Temporal behavior of a solute cloud in a chemically heterogeneous porous medium, J. Fluid Mech., 386, 77-104.
[2] Bellin, A., M. Pannone, A. Fiori, and A. Rinaldo (1996), On transport in porous formations characterized by heterogeneity of evolving scales, Water Resour. Res., 32, 3485-3496.
[3] Cintoli, S., S. P. Neuman, and V. Di Federico (2005), Generating and scaling fractional Brownian motion on finite domains, Geophys. Res. Lett. 32, L08404, doi:10.1029/2005GL022608.
[4] Dagan, G. (1994), The significance of heterogeneity of evolving scales and of anomalous diffusion to transport in porous formations, Water Resour. Res., 30, 33273336, 1994.
[5] Dagan, G. (1987), Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, 183215.
[6] Dagan, G. (1988), Time-dependent macrodispersion for solute transport in anisotropic heterogeneous aquifers, Water Resour. Res., 24, 1491-1500.
[7] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000), Temporal behavior of a solute cloud in a heterogeneous porous medium 1. Point-like injection, Water Resour. Res., 36, 35913604.
[8] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000), Temporal behavior of a solute cloud in a heterogeneous porous medium 2. Spatially extended injection, Water Resour. Res., 36, 3605-3614.
[9] Di Federico, V., and S. P. Neuman (1997), Scaling of random fields by means of truncated power variograms and associated spectra, Water Resour. Res., 33, 1075-1085.
[10] Fiori, A. (1996), Finite Peclet extensions of Dagan's solutions to transport in anisotropic heterogeneous formations, Water Resour. Res., 32, 193198.
[11] Fiori, A. (2001), On the influence of local dispersion in solute transport through formations with evolving scales of heterogeneity, Water Resour. Res., 37, 235242.
[12] Fiori, A., and G. Dagan (2000), Concentration fluctuations in aquifer transport: A rigorous firstorder solution and applications, J. Contam. Hydrol., 45, 139163.
[13] Gelhar, L. W. (1986), Stochastic subsurface hydrology from theory to applications, Water Resour. Res., 22, 135S-145S.
[14] Gelhar, L. W., C. L. and Axness (1993), Three-dimensional stochastic analysis of macrodispersion in aquifers, textitWater Resour. Res., 19, 161-180.
[15] Gradshteyn, I. S., and I. M. Ryzhik (2007) Table of Integrals, Series, and Products, Elsevier, Amsterdam.
[16] Jeon, J.-H., and R. Metzler (2010), Fractional Brownian motion and motion governed by the fractional Langevin equation in confined geometries, Phys. Rev. E 81, 021103, doi:10.1103/PhysRevE.81.021103.
[17] McLaughlin, D., and F. Ruan (2001), Macrodispersivity and large-scale hydrogeologic variability, Transp. Porous Media, 42, 133154.
[18] Papoulis, A., and S. U. Pillai (2009), Probability, Random Variables and Stochastic Processes, McGraw-Hill, New York.
[19] Ross, K., and S. Attinger (2010), Temporal behaviour of a solute cloud in a fractal heterogeneous porous medium at different scales, paper presented at EGU General Assembly 2010, Vienna, Austria, 02-07 May 2010.
[20] Schwarze, H., U. Jaekel, and H. Vereecken (2001), Estimation of macrodispersivity by different approximation methods for flow and transport in randomly heterogeneous media, Transp. Porous Media, 43, 265287.
[21] Suciu N. (2010), Spatially inhomogeneous transition probabilities as memory effects for diffusion in statistically homogeneous random velocity fields, Phys. Rev. E, 81, 056301, doi:10.1103/PhysRevE.81.056301.
[22] Suciu, N., C. Vamoş, J. Vanderborght, H. Hardelauf, and H. Vereecken (2006), Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi:10.1029/2005WR004546.
[23] Suciu N., C. Vamos, F. A. Radu, H. Vereecken, and P. Knabner (2009), Persistent memory of diffusing particles,Phys. Rev. E, 80, 061134, doi:10.1103/PhysRevE.80.061134.
[24] Vamoş, C., N. Suciu, H. Vereecken, J. Vanderborht, and O. Nitzsche (2001), Path decomposition of discrete effective diffusion coecient, Internal Report ICG-IV. 00501, Research Center Jülich.
[25] Vamoş, C., N. Suciu, and H. Vereecken (2003), Generalized random walk algorithm for the numerical modeling of complex diusion processes, J. Comput. Phys., 186, 527-544, doi:10.1016/S0021-9991(03)00073-1.
[26] Vanderborght, J. (2001), Concentration variance and spatial covariance in second-order stationary heterogeneous conductivity fields, Water Resour. Res., 37, 1893-1912.
2011

Related Posts