Dependence on initial conditions, memory effects, and ergodicity of transport in heterogeneous media

Nicolae Suciu
html, paper

Abstract

Authors

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

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

K. Sabelfeld

H. Vereecken

P. Knabner

Keywords

Cite this paper as:

N. Suciu, C. Vamoş, H. Vereecken, K. Sabelfeld, P. Knabner, Dependence on initial conditions, memory effects, and ergodicity of transport in heterogeneous media, Preprint No. 324, 2008, Institute of Applied Mathematics, Friedrich-Alexander University Erlangen-Nuremberg

References

see the expansion block below.

PDF

pr324.pdf

About this paper

Journal

Preprint

Publisher Name
DOI
Print ISSN
Online ISSN

MR

?

ZBL

?

[1] Balescu, R. (2000), Memory e ects in plasma transport theory, Plasma Phys. Control. Fusion 42, B1B13.
[2] Balescu, R., H.-D. Wang, and J. H. Misguich (1994), Langevin equation versus kinetic equation:Subdiffusive usive behavior of charged particles in a stochastic magnetic field, Phys. Plasmas 1(12), 38263842.
[3] Berkowitz, B., and H. Scher (2001), Probabilistic approaches to transport in heterogeneous media, Transport in Porous Media, 42, 241263.
[4] Bhattacharya, R. N., and V. K. Gupta (1983), A Theoretical Explanation of Solute Dispersion in Saturated Porous Media at the Darcy Scale, Water Resour. Res., 19, 934944.
[5] Bouchaud, J.-P., and A. Georges (1990), Anomalous di usion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep., 195, 127 293.
[6] Brouwers, J. J. H. (2002), On di usion theory in turbulence, J. Eng. Mat., 44, 277-295.
[7] Chiles, J. P., and P. Del¯ner (1999), Geostatisctics: Modeling Spatial Uncertainty, John Wiley& Sons, New York.
[8] Clincy, M., and H. Kinzelbach (2001), Strati¯ed disordered media: exact solutions for transport parameters and their self-averaging properties, J. Phys. A: Math. Gen. 34, 71427152.
[9] Compte, A., and M. O. Caceres (1998), Fractional dynamics in random velocity ¯elds, Phys. Rev. Lett. 81(15),3140-3143.
[10] Cortis, A., H. Scher, and B. Berkowitz (2004), Numerical simulation of non-Fickian transport in geological formations with multiple scale heterogeneity, Water Resour. Res, 40, W04209, doi:10.1029/2003WRR002750.
[11] Dagan, G. (1987), Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, 183215.
[12] Dagan, G. (1990), Transport in heterogeneous porous formations: Spatial moments, ergodicity and e ective dispersion, Water Resour. Res., 26, 12811290.
[13] Dagan, G. (1991), Dispersion of a solute in non-ergodic transport in steady velocity ¯elds in
heterogeneous formations, J. Fluid Mech., 233, 197 210.
[14] Dagan, G. (2004), Comment on `Exact averaging of stochastic equations for transport in random velocity field, Transport in porous media 50, 223 241, 2003, and Probability density functions for solute transport in random field Transport in porous media 50, 243266, 2003, Transport in porous media, 55, 113 116. 23
[15] Demmy, G., S. Berglund, and W. Graham (1999), Injection mode implications for solute transport in porous media: Analysis in a stochastic Lagrangian framework, Water Resour. Res., 35(7), 1965 1973.
[16] 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.
[17] Dentz, M., and J. Carrera (2005), E ective solute transport in temporally uctuating ow through heterogeneous media, Water Resour. Res., 41, W08414, doi:10.1029/2004WR003571.
[18] Dentz, M., and J. Carrera (2007), Mixing and spreading in strati¯ed ow, Phys. Fluids 19, 017107, doi:10.1063/1.2427089.
[19] Doob, J. L. (1953), Stochastic Processes, John Wiley & Sons, London.
[20] Eberhard J., N. Suciu, and C. Vamos (2007), On the self-averaging of dispersion for transport in quasi-periodic random media, J. Phys. A: Math. Theor., 40, 597-610, doi: 10.1088/1751 8113/40/4/002.
[21] Fannjiang, A., and T. Komorowski (1999), Turbulent di usion in Markovian ows, Ann. Appl. Probab., 9(3), 591-610.
[22] Fiori, A., and G. Dagan (2000), Concentration fluctuations in aquifer transport: a rigorous first-order solution and applications, J.Contam. Hydrol. 45, 139 163.
[23] Fiori, A., and I. Jancovic (2005), Can we determine the transverse macrodis persivity by using the method of moments?, Adv. Water Resour., 28, 589599,
doi:10.1016/j.advwaters.2004.09.909.
[24] Gardiner, C. W. (1985), Handbook of Stochastic Methods (for Physics, Chemistry and Natural Science), Springer, New York.
[25] Gleeson, J. P. (2002), Comment on Di usion in biased turbulence , Phys. Rev. E, 66, 038301.
[26] Jaekel, U., and H. Vereecken (1997), Renormalization group analysis of macrodispersion in a directed random ow, Water Resour. Res., 33, 2287 2299.
[27] Kitanidis, P. K. (1988), Prediction by the method of moments of transport in a heterogeneous formation, J. Hydrol., 102, 453 473.
[28] Kloeden, P. E. and E, Platen (1995), Numerical Solutions of Stochastic Di erential Equations, Springer, Berlin.
[29] Komorowski, T., and G. Papanicolaou (1997), Motion in a Gaussian incompressible ow, Ann. Appl. Probab., 7(1), 229-264.
[30] Kramer, P. R., and A. J. Majda (2007), Lectures on Turbulent Di usion, Springer (in preparation).
[31] LaBolle, E. M., J. Quastel, G. E. Fogg, and J. Gravner (2000), Di usion processes in composite media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coe cients, Water Resour. Res., 36(3), 651662.
[32] Le Doussal, P., and J. Machta (1989), Annealed versus quenched di usion coe cient in randommedia, Phys. Rev. B, 40(12), 9427 9430.24
[33] Lichtner, P. C., S. Kelkar, and B. Robinson (2002), New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking, Water Resour. Res., 38(8), 10.1029/2000WR000100.
[34] Lumley, J. L. (1962a), An approach to the Eulerian-Lagrangian problem, J. Math. Phys., 3(2),309-312.
[35] Lumley, J. L. (1962b), The mathematical nature of the problem of relating Lagrangian and Eulerian statistical functions in turbulence, pp , 17-26 in Mecanique de la Turbulence (Coll.Intern. du CNRS a Marseille), Ed. CNRS, Paris.
[36] Lundgren, T. S. (1981), Turbulent pair dispersion and scalar di usion, J. Fluid Mech. 111, 27-57.
[37] Matheron, G., and G. de Marsily (1980), Is transport in porous media always di usive?, Water Resour. Res., 16, 901 917.
[38] Monin, A. S., and A. M. Yaglom (1975) Statistical Fluid Mechanics: Mechanics of Turbulence, MIT Press, Cambridge, M A.
[39] Morales-Casique, E., S.P. Neuman, and A. Gaudagnini (2006), Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: Theoretical framework, Adv. Water Resour. 29(8), 1238-1255, doi:10.1016/j.advwatres.2005.10.002.
[40] Morgado, R., F. A. Oliveira, G. G. Batrouni, and A. Hansen (2002), Relation between anomalous and normal difusion in systems with memory, Phys. Rev. Lett., 89(10), 100601.
[41] Mokshin, A. V., Yulmetyev, R. M., and P. Hnggi (2005), Simple measure of memory for dynamical processes described by a generalized Langevin equation, Phys. Rev. Lett., 95(20), 200601.
[42] Na , R. L., D. F. Haley, and E. A. Sudicky (1998), High-resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media 2. Transport results, Water Resour. Res., 34(4), 679 697.
[43] Port, S. C., and C. J. Stone (1976), Random measures and their application to motion in an incompressible uid. J. Appl. Prob., 13, 498506.
[44] Sposito, G. (2001), Topological groundwater hydrodynamics, Adv. Water Resour., 24,, 793 801.
[45] Sposito, G. and G. Dagan (1994), Predicting solute plume evolution in heterogeneous porous formations, Water Resour. Res., 30(2), 585589.
[46] Suciu N., C. Vamos, J. Vanderborght, H. Hardelauf, and H. Vereecken (2006a), Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi:10.1029/2005WR004546.
[47] Suciu N., C. Vamos, and J. Eberhard (2006b), Evaluation of the ¯rst-order approx imations for transport in heterogeneous media, Water Resour. Res. 42, W11504, doi: 10.1029/2005WR004714.
[48] Suciu N., and C. Vamos (2007), Comment on Nonstationary ow and nonergodic transport in random porous media by G. Darvini and P. Salandin, Water Resour. Res., 43, W12601, doi:10.1029/2007WR005946. 25
[49] Suciu N., C. Vamos, and H. Vereecken (2007a), Multiple meanings of ergodicity in real life problems, pp.196-211 in Proceedings of the 5th Workshop on Mathematical Modelling of Environmental and Life Sciences Problems, Ed. G. Marinoschi, S. Ion and C. Popa, Rom. Acad.
Publishing House, Bucharest.
[50] Suciu N., C. Vamos, and K. Sabelfeld (2007b), Ergodic simulations of diffusion in random velocity fields, pp. 659-668 in Monte Carlo and Quasi-Monte Carlo Methods 2006, Ed. A. Keller, S. Heinrich, and H. Niederreiter, Springer Verlag, Heidelberg.
[51] Suciu N., C. Vamos, K. Sabelfeld, and. C. Andronache (2007c), Memory effects and ergodicity for diffsion in spatially correlated velocity fields, Proc. Appl. Math. Mech., 7, 20100152010016, doi:10.1002/pamm.20070057.
[52] Suciu N., F. A. Radu, and C. Vamos (2008) Lagrangian-Eulerian statistics and numerical modeling of transport in homogeneous random media, in Proceedings of the 6th Workshop on Mathematical Modelling of Environmental and Life Sciences Problems, Ed. G. Marinoschi, S.Ion, and C. Popa, Rom. Acad. Publishing House, Bucharest (to appear).
[53] Trefry, M. G., F. P. Ruan, and D. McLaughlin (2003), Numerical simulations of preasymptotic transport in heterogeneous porous media: Departures from the Gaussian limit, Water Resour. Res., 39(3), 1063, doi:10.1029/2001WRR001101.
[54] Vamos, C., N. Suciu, and H. Vereecken (2003), Generalized random walk algorithm for the numerical modeling of complex di usion processes, J. Comp. Phys., 186(2), 527 544, doi:10.1016/S0021-9991(03)00073-1.
[55] Vamos, C., N. Suciu, H. Vereecken, J. Vanderborht, and O. Nitzsche (2001), Path decomposi tion of discrete e ective di usion coe cient, Internal Report ICG-IV.00501, Research Center Julich.
[56] van Kampen, N. G. (1981), Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam.
[57] Wentzell, A. D. (1981), A Course in the Theory of Stochastic Processes, McGraw-Hill, New York.
[58] Yaglom, A. M. (1987) Correlation Theory of Stationary and Related Random Functions, Vol ume I: Basic Results, Springer-Verlag, New York.
[59] Zirbel, C., L. (2001), Lagrangian observations of homogeneous random environments, Adv. Appl. Prob., 33, 810-835.
[60] Zhang, Y.-K., and B.-m. Seo (2004), Stochastic analyses and Monte-Carlo simulations of non ergodic solute transport in three-dimensional heterogeneous statistically anisotropic aquifers, Water Resour. Res., 40, W05103, doi:10.1029/2003WR002871.
[61] Zwanzig, R. (1961), Memory e ects and irreversible thermodynamics, Phys. Rev. 124(4), 983 992.

soon

Paper (preprint) in HTML form

pr324

Dependence on Initial Conditions, Memory Effects, and Ergodicity of Transport in Heterogeneous Media

N. Suciu 1 1 ^(1){ }^{1}1, C. Vamoş 2 2 ^(2){ }^{2}2, H. Vereecken 3 3 ^(3){ }^{3}3, K. Sabelfeld 4 , 5 4 , 5 ^(4,5){ }^{4,5}4,5, and P. Knabner 1 1 ^(1){ }^{1}1 1 1 ^(1){ }^{1}1 Friedrich-Alexander University of Erlangen-Nuremberg, Institute of Applied Mathematics, Martensstrasse 3, 91058 Erlangen, Germany. 2 2 ^(2){ }^{2}2 Romanian Academy, Tiberiu Popoviciu Institute of Numerical Analysis, 400320 Cluj-Napoca, P. O. Box 68-1, Romania. 3 3 ^(3){ }^{3}3 Research Center Jülich, ICG-IV: Agrosphere Institute, 52425 Jülich, Germany. 4 4 ^(4){ }^{4}4 Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstr., 39, 10117 Berlin, Germany. 5 5 ^(5){ }^{5}5 Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch of Russian Academy of Sciences, Prospect Lavrentjeva 6, 630090 Novosibirsk, Russia

Abstract

Memory effects in the pre-asymptotic regime of transport in heterogeneous media, indicated by reliable numerical experiments, are explained by a dependence on initial conditions which is specific for systems with space variable properties such as aquifers, turbulent atmosphere or ionized plasmas. For statistically homogeneous random velocity fields, with finite correlation range, incompressible, and continuously differentiable, the one-particle dispersion (average over velocity realizations of the actual dispersion for fixed initial positions of the solute molecules) is a "memory-free" quantity, solely determined by the velocity correlation and the local dispersion coefficient and independent of initial conditions. Nonergodic behavior of large initial plumes, as manifest in deviations of actual solute dispersion from one-particle dispersion, is associated with a "memory term" consisting of correlations between initial positions and displacements of solute molecules. Increasing the source dimensions has two opposite effects: it reduces the uncertainty related to the randomness of center of mass, but, at the same time, it yields large memory terms. The memory effects increase with the source dimension and depend on its shape and orientation. Large narrow sources oriented transverse to the mean flow direction yield ergodic behavior with respect to the one-particle dispersion of the longitudinal dispersion and nonergodic behavior of the transverse dispersion, whereas for large longitudinal sources the longitudinal dispersion behaves nonergodically and the transverse dispersion ergodically. Such memory effects are significantly large over hundreds of heterogeneity scales and should therefore be considered in practical applications, as for instance calibration of model parameters, forecasting, identification of the contaminant source.

KEYWORDS: memory effects, ergodicity, macrodispersion, Itô equation, global random walk.

1 Introduction

The transport in heterogeneous media is conveniently characterized by the second central spatial moment tensor of the concentration field s l m ( l , m = 1 , 2 , 3 ) s l m ( l , m = 1 , 2 , 3 ) s_(lm)(l,m=1,2,3)s_{l m}(l, m=1,2,3)slm(l,m=1,2,3). As well as providing a measure for the spatial extension of the solute plume, the dependence on time t t ttt of the second moment is commonly used to investigate whether the transport is diffusive, i.e. s l m ∼ t s l m ∼ t s_(lm)∼ts_{l m} \sim tslm∼t [Sposito and Dagan, 1994]. Since it can be estimated, by either analytical approximations or numerical simulations, without solving the transport equations [Suciu et al., 2006b; Eberhard et al., 2007] this quantity is particularly useful in investigations on pre-asymptotic transport regime, for which generally there are no close form solutions [see e.g. Morales-Casique et al., 2006].
A frequently used approach considers a local dispersion process and models the heterogeneity of the velocity at larger scales by a space random function. Let X l X l X_(l)X_{l}Xl be the l l lll component of the trajectory of a solute molecule, governed by an advection-dispersion process. To simplify matters, we consider only the diagonal components l l l l lll lll of the various second moments. For a fixed realization of the velocity, the second moment s l l s l l s_(ll)s_{l l}sll of the actual concentration is given by the dispersion, or the
mean-square displacement of the molecules at a given time
(1) s l l = ⟨ [ X l − ⟨ X l ⟩ D X 0 ] 2 ⟩ D X 0 (1) s l l = X l − X l D X 0 2 D X 0 {:(1)s_(ll)=(:[X_(l)-(:X_(l):)_(DX_(0))]^(2):)_(DX_(0)):}\begin{equation*} s_{l l}=\left\langle\left[X_{l}-\left\langle X_{l}\right\rangle_{D X_{0}}\right]^{2}\right\rangle_{D X_{0}} \tag{1} \end{equation*}(1)sll=⟨[Xl−⟨Xl⟩DX0]2⟩DX0
The subscripts D D DDD and X 0 X 0 X_(0)X_{0}X0 in (1) denote respectively the average over the realizations of the local dispersion and the space average with respect to the initial distribution of molecules.
Various equivalent forms of (1) can be derived to emphasize the role of different factors which influence the dispersion. The use in (1) of the displacement X ~ l = X l − X 0 l X ~ l = X l − X 0 l widetilde(X)_(l)=X_(l)-X_(0l)\widetilde{X}_{l}=X_{l}-X_{0 l}X~l=Xl−X0l relative to the initial position X 0 l X 0 l X_(0l)X_{0 l}X0l yields
(2) s l l = S l l ( 0 ) + ⟨ [ X ~ l − ⟨ X ~ l ⟩ D X 0 ] 2 ⟩ D X 0 + m l l , (2) s l l = S l l ( 0 ) + X ~ l − X ~ l D X 0 2 D X 0 + m l l , {:(2)s_(ll)=S_(ll)(0)+(:[ widetilde(X)_(l)-(: widetilde(X)_(l):)_(DX_(0))]^(2):)_(DX_(0))+m_(ll)",":}\begin{equation*} s_{l l}=S_{l l}(0)+\left\langle\left[\widetilde{X}_{l}-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0}}\right]^{2}\right\rangle_{D X_{0}}+m_{l l}, \tag{2} \end{equation*}(2)sll=Sll(0)+⟨[X~l−⟨X~l⟩DX0]2⟩DX0+mll,
where S l l ( 0 ) = ⟨ [ X 0 l − ⟨ X 0 l ⟩ X 0 ] 2 ⟩ X 0 S l l ( 0 ) = X 0 l − X 0 l X 0 2 X 0 S_(ll)(0)=(:[X_(0l)-(:X_(0l):)_(X_(0))]^(2):)_(X_(0))S_{l l}(0)=\left\langle\left[X_{0 l}-\left\langle X_{0 l}\right\rangle_{X_{0}}\right]^{2}\right\rangle_{X_{0}}Sll(0)=⟨[X0l−⟨X0l⟩X0]2⟩X0 is the deterministic initial second moment. The last term in (2),
(3) m l l = 2 ⟨ [ X 0 l − ⟨ X 0 l ⟩ X 0 ] X ~ l ⟩ D X 0 (3) m l l = 2 X 0 l − X 0 l X 0 X ~ l D X 0 {:(3)m_(ll)=2(:[X_(0l)-(:X_(0l):)_(X_(0))] widetilde(X)_(l):)_(DX_(0)):}\begin{equation*} m_{l l}=2\left\langle\left[X_{0 l}-\left\langle X_{0 l}\right\rangle_{X_{0}}\right] \widetilde{X}_{l}\right\rangle_{D X_{0}} \tag{3} \end{equation*}(3)mll=2⟨[X0l−⟨X0l⟩X0]X~l⟩DX0
describes spatial correlations between relative displacements on trajectories and initial positions and, therefore, it has been called "memory term" [Suciu and Vamoş, 2007; Suciu et al., 2007c, 2008]. The implications of such terms for predicting solute dispersion were first investigated by Sposito and Dagan [1994], in the context of deterministic advective transport. A decade later, Fiori and Jancović [2005] simulated the time behavior of m l l m l l m_(ll)m_{l l}mll, for advective transport as well, and concluded that a representative transverse dispersivity cannot be inferred from experiments done for large transverse sources. Recent numerical investigations which considered both advection and local dispersion [Suciu et al, 2006a] also indicated that the transverse dispersion for large transverse sources significantly differs from theoretical "ergodic" results. The goal of this paper is to study the role of memory terms on the nonergodic behavior of the pre-asymptotic dispersion. This study is mainly based on "global random walk" (GRW) numerical simulations [Vamoş et al., 2003], for different initial conditions, which complete previous ones presented in [Suciu et al, 2006a].
Averages over ensembles of realizations of the random velocity field, noted in the following by a subscript V V VVV, yield idealized descriptions of solute dispersion. Because in practice we have a single velocity field, the ergodicity issue has to be considered in stochastic modeling. In other words, the relevance of the stochastic description for the dispersion observable in a real case has to be assessed. To proceed, we add and subtract ⟨ X ~ l ⟩ D X 0 V 2 X ~ l D X 0 V 2 (: widetilde(X)_(l):)_(DX_(0)V)^(2)\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}^{2}⟨X~l⟩DX0V2 in (2) and express the actual dispersion as
(4) s l l = S l l ( 0 ) + x ~ l l + m l l − r l l (5) x ~ l l = ⟨ [ X ~ l − ⟨ X ~ l ⟩ D X 0 V ] 2 ⟩ D X 0 (6) r l l = [ ⟨ X ~ l ⟩ D X 0 − ⟨ X ~ l ⟩ D X 0 V ] 2 (4) s l l = S l l ( 0 ) + x ~ l l + m l l − r l l (5) x ~ l l = X ~ l − X ~ l D X 0 V 2 D X 0 (6) r l l = X ~ l D X 0 − X ~ l D X 0 V 2 {:[(4)s_(ll)=S_(ll)(0)+ widetilde(x)_(ll)+m_(ll)-r_(ll)],[(5) widetilde(x)_(ll)=(:[ widetilde(X)_(l)-(: widetilde(X)_(l):)_(DX_(0)V)]^(2):)_(DX_(0))],[(6)r_(ll)=[(: widetilde(X)_(l):)_(DX_(0))-(: widetilde(X)_(l):)_(DX_(0)V)]^(2)]:}\begin{align*} s_{l l} & =S_{l l}(0)+\widetilde{x}_{l l}+m_{l l}-r_{l l} \tag{4}\\ \widetilde{x}_{l l} & =\left\langle\left[\widetilde{X}_{l}-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\right]^{2}\right\rangle_{D X_{0}} \tag{5}\\ r_{l l} & =\left[\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0}}-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\right]^{2} \tag{6} \end{align*}(4)sll=Sll(0)+x~ll+mll−rll(5)x~ll=⟨[X~l−⟨X~l⟩DX0V]2⟩DX0(6)rll=[⟨X~l⟩DX0−⟨X~l⟩DX0V]2
The average of (4) over the ensemble of velocity realizations gives the expected second moment S l l = ⟨ s l l ⟩ V S l l = s l l V S_(ll)=(:s_(ll):)_(V)S_{l l}=\left\langle s_{l l}\right\rangle_{V}Sll=⟨sll⟩V,
(7) S l l = S l l ( 0 ) + X ~ l l + M l l − R l l (7) S l l = S l l ( 0 ) + X ~ l l + M l l − R l l {:(7)S_(ll)=S_(ll)(0)+ widetilde(X)_(ll)+M_(ll)-R_(ll):}\begin{equation*} S_{l l}=S_{l l}(0)+\widetilde{X}_{l l}+M_{l l}-R_{l l} \tag{7} \end{equation*}(7)Sll=Sll(0)+X~ll+Mll−Rll
where X ~ l l = ⟨ x ~ l l ⟩ V , M l l = ⟨ m l l ⟩ V X ~ l l = x ~ l l V , M l l = m l l V widetilde(X)_(ll)=(: widetilde(x)_(ll):)_(V),M_(ll)=(:m_(ll):)_(V)\widetilde{X}_{l l}=\left\langle\widetilde{x}_{l l}\right\rangle_{V}, M_{l l}=\left\langle m_{l l}\right\rangle_{V}X~ll=⟨x~ll⟩V,Mll=⟨mll⟩V, and R l l = ⟨ r l l ⟩ V R l l = r l l V R_(ll)=(:r_(ll):)_(V)R_{l l}=\left\langle r_{l l}\right\rangle_{V}Rll=⟨rll⟩V. While the last two terms of (7) have the straightforward meaning of mean memory term ( M l l M l l M_(ll)M_{l l}Mll ) and that of variance of actual center of mass ( R l l ) R l l (R_(ll))\left(R_{l l}\right)(Rll), the term X ~ l l X ~ l l widetilde(X)_(ll)\widetilde{X}_{l l}X~ll still needs to be related to other meaningful ensemble quantities.
Hereinafter, we assume all necessary joint measurability conditions which allow permutations of averages and integrals [Wentzell, 1981, p.22; Zirbel, 2001, p.814]. Then, X ~ l l X ~ l l widetilde(X)_(ll)\widetilde{X}_{l l}X~ll can be calculated by ensemble averaging (5) as follows:
(8) X ~ l l = ⟨ X l l ⟩ X 0 + Q l l , (9) X l l = ⟨ [ X ~ l − ⟨ X ~ l ⟩ D V ] 2 ⟩ D V , (10) Q l l = ⟨ [ ⟨ X ~ l ⟩ D V − ⟨ X ~ l ⟩ D X 0 V ] 2 ⟩ X 0 . (8) X ~ l l = X l l X 0 + Q l l , (9) X l l = X ~ l − X ~ l D V 2 D V , (10) Q l l = X ~ l D V − X ~ l D X 0 V 2 X 0 . {:[(8) widetilde(X)_(ll)=(:X_(ll):)_(X_(0))+Q_(ll)","],[(9)X_(ll)=(:[ widetilde(X)_(l)-(: widetilde(X)_(l):)_(DV)]^(2):)_(DV)","],[(10)Q_(ll)=(:[(: widetilde(X)_(l):)_(DV)-(: widetilde(X)_(l):)_(DX_(0)V)]^(2):)_(X_(0)).]:}\begin{gather*} \widetilde{X}_{l l}=\left\langle X_{l l}\right\rangle_{X_{0}}+Q_{l l}, \tag{8}\\ X_{l l}=\left\langle\left[\widetilde{X}_{l}-\left\langle\widetilde{X}_{l}\right\rangle_{D V}\right]^{2}\right\rangle_{D V}, \tag{9}\\ Q_{l l}=\left\langle\left[\left\langle\widetilde{X}_{l}\right\rangle_{D V}-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\right]^{2}\right\rangle_{X_{0}} . \tag{10} \end{gather*}(8)X~ll=⟨Xll⟩X0+Qll,(9)Xll=⟨[X~l−⟨X~l⟩DV]2⟩DV,(10)Qll=⟨[⟨X~l⟩DV−⟨X~l⟩DX0V]2⟩X0.
Now, X ~ l l X ~ l l widetilde(X)_(ll)\widetilde{X}_{l l}X~ll is expressed by (8) as sum of the space mean of one-particle displacements variance X l l X l l X_(ll)X_{l l}Xll (defined by averaging with respect to D D DDD and V V VVV for a fixed initial position) and the spatial variance (computed by averages over X 0 X 0 X_(0)X_{0}X0 ) of the one-particle center of mass ⟨ X ~ l ⟩ D V X ~ l D V (: widetilde(X)_(l):)_(DV)\left\langle\widetilde{X}_{l}\right\rangle_{D V}⟨X~l⟩DV. From (7) and (8) one obtains
(11) S l l = S l l ( 0 ) + ⟨ X l l ⟩ X 0 + M l l − R l l + Q l l . (11) S l l = S l l ( 0 ) + X l l X 0 + M l l − R l l + Q l l . {:(11)S_(ll)=S_(ll)(0)+(:X_(ll):)_(X_(0))+M_(ll)-R_(ll)+Q_(ll).:}\begin{equation*} S_{l l}=S_{l l}(0)+\left\langle X_{l l}\right\rangle_{X_{0}}+M_{l l}-R_{l l}+Q_{l l} . \tag{11} \end{equation*}(11)Sll=Sll(0)+⟨Xll⟩X0+Mll−Rll+Qll.
The general expression (11) of the expected dispersion applies in the purely advective case as well, with the single difference that the average over D D DDD is missing in the definitions of its terms [ [ [[[ Suciu and Vamoş, 2007].
A candidate idealized description of the average second moment is the one-particle dispersion X l l X l l X_(ll)X_{l l}Xll, provided that it does not depend on initial conditions,
(12) S l l ≈ S l l ( 0 ) + X l l (12) S l l ≈ S l l ( 0 ) + X l l {:(12)S_(ll)~~S_(ll)(0)+X_(ll):}\begin{equation*} S_{l l} \approx S_{l l}(0)+X_{l l} \tag{12} \end{equation*}(12)Sll≈Sll(0)+Xll
One expects that (12) applies for plumes with large initial dimensions, often called "ergodic plumes", [e.g. Dagan, 1990, 1991, 2004; Zhang and Seo, 2004]. Even if (12) can approximate the expectation, it might not be accurate for the actual dispersion. Indeed, (3) shows that for large sources the memory terms can be important. So, at best, the actual dispersion can be approximated as
(13) s l l ≈ S l l ( 0 ) + X l l + m l l (13) s l l ≈ S l l ( 0 ) + X l l + m l l {:(13)s_(ll)~~S_(ll)(0)+X_(ll)+m_(ll):}\begin{equation*} s_{l l} \approx S_{l l}(0)+X_{l l}+m_{l l} \tag{13} \end{equation*}(13)sll≈Sll(0)+Xll+mll
As shown below, the approximations (12-13) can be adopted for suitable properties of the random velocity field and the ergodicity issue can be answered by estimating the memory terms.
The statistics of terms in (11) depends on that of the "Lagrangian velocity process", which consists of observations of the random velocity field at random locations on the trajectory of the advection-dispersion process [Zirbel, 2001]. If the local dispersion is described by a Wiener process and the "Eulerian velocity" (observed at fixed space points) is modeled as a statistically homogeneous random field, incompressible, and continuously differentiable, then the Lagrangian and Eulerian one-point statistics coincide [Port and Stone, 1976; Zirbel, 2001; Kramer and Majda, 2007]. This result implies the independence on X 0 X 0 X_(0)X_{0}X0 of ⟨ X ~ l ⟩ D V X ~ l D V (: widetilde(X)_(l):)_(DV)\left\langle\widetilde{X}_{l}\right\rangle_{D V}⟨X~l⟩DV and X l l X l l X_(ll)X_{l l}Xll [Suciu et al., 2007c, 2008]. Since under the conditions stated above the one-particle dispersion X l l X l l X_(ll)X_{l l}Xll no longer depends on X 0 X 0 X_(0)X_{0}X0, it is a memory-free quantity. The independence on X 0 X 0 X_(0)X_{0}X0 of ⟨ X ~ l ⟩ D V X ~ l D V (: widetilde(X)_(l):)_(DV)\left\langle\widetilde{X}_{l}\right\rangle_{D V}⟨X~l⟩DV leads to the cancelation of M l l M l l M_(ll)M_{l l}Mll and Q l l Q l l Q_(ll)Q_{l l}Qll. However, there are no special consequences for R l l R l l R_(ll)R_{l l}Rll, which depends on two-points Lagrangian statistics [Suciu et al., 2008]. The same conclusions can be drawn in the case of purely advective transport from the earlier result of Lumley [1962b]. In these conditions, (11) reduces to the simpler form
(14) S l l = S l l ( 0 ) + X l l − R l l (14) S l l = S l l ( 0 ) + X l l − R l l {:(14)S_(ll)=S_(ll)(0)+X_(ll)-R_(ll):}\begin{equation*} S_{l l}=S_{l l}(0)+X_{l l}-R_{l l} \tag{14} \end{equation*}(14)Sll=Sll(0)+Xll−Rll
This relation was obtained in the hydrological context for the first time by Dagan [1990, equation (11)], under the hypothesis of "Lagrangian stationarity" for advective transport, and it is retrieved by different approximations of the plume moments [Dagan, 1991; Suciu and Vamoş, 2007] or of the equivalent dispersion coefficients [Dentz et al., 2000; Suciu et al., 2006a, b].
Equation (14) shows that the expectation of s l l − S l l ( 0 ) s l l − S l l ( 0 ) s_(ll)-S_(ll)(0)s_{l l}-S_{l l}(0)sll−Sll(0) differs from X l l X l l X_(ll)X_{l l}Xll by − R l l − R l l -R_(ll)-R_{l l}−Rll, i.e. the actual dispersion is a biased estimator of the one-particle dispersion. Therefore, the ergodicity of the actual dispersion with respect to the one-particle dispersion is quantified in the mean square norm by both the bias of its expected value and by its standard deviation [Suciu et al., 2006a, relations (5)]. This ergodic property in a large sense differs from the common definition of ergodicity [Suciu et al., 2007a, b]. The latter characterizes, for instance, homogeneous random fields with finite correlation range, for which space averages of the velocity are unbiased estimators strongly convergent to the ensemble mean velocity [Chilès and Delfiner, 1999]. Ergodicity of the random velocity field can be helpful in the present context if we note that (5) and the center of mass ⟨ X ~ l ⟩ D X 0 X ~ l D X 0 (: widetilde(X)_(l):)_(DX_(0))\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0}}⟨X~l⟩DX0 in (6) are expressed as space averages with respect to X 0 X 0 X_(0)X_{0}X0. Arguments built on the assumption that the velocity field is ergodic and Gaussian suggest that, for large sources, these space averages can be replaced by ensemble averages [Suciu et al., 2007c]. Thus, as follows from (6), r l l ≈ 0 r l l ≈ 0 r_(ll)~~0r_{l l} \approx 0rll≈0, the term
(5) approximates the one-particle dispersion (9), x ~ l l ≈ X l l x ~ l l ≈ X l l widetilde(x)_(ll)~~X_(ll)\widetilde{x}_{l l} \approx X_{l l}x~ll≈Xll, and from (14) and (4) one obtains, respectively, the approximations (12) and (13).
In numerical approaches based on simulated concentration fields [Naff et al., 1998; Trefry et al., 2003; Suciu et al., 2006a], as well as in experiments, it is not possible to use trajectories of individual solute molecules to compute one-particle quantities. Therefore, global descriptions of the solute dispersion have to be used. These are obtained by gathering together the first three terms of (4),
(15) σ l l = S l l ( 0 ) + x ~ l l + m l l = ⟨ [ X l − ⟨ X l ⟩ D X 0 V ] 2 ⟩ D X 0 , (15) σ l l = S l l ( 0 ) + x ~ l l + m l l = X l − X l D X 0 V 2 D X 0 , {:(15)sigma_(ll)=S_(ll)(0)+ widetilde(x)_(ll)+m_(ll)=(:[X_(l)-(:X_(l):)_(DX_(0)V)]^(2):)_(DX_(0))",":}\begin{equation*} \sigma_{l l}=S_{l l}(0)+\widetilde{x}_{l l}+m_{l l}=\left\langle\left[X_{l}-\left\langle X_{l}\right\rangle_{D X_{0} V}\right]^{2}\right\rangle_{D X_{0}}, \tag{15} \end{equation*}(15)σll=Sll(0)+x~ll+mll=⟨[Xl−⟨Xl⟩DX0V]2⟩DX0,
which leads to the following expressions for the actual and ensemble average dispersion:
(16) s l l = σ l l − r l l (17) S l l = Σ l l − R l l (16) s l l = σ l l − r l l (17) S l l = Σ l l − R l l {:[(16)s_(ll)=sigma_(ll)-r_(ll)],[(17)S_(ll)=Sigma_(ll)-R_(ll)]:}\begin{align*} s_{l l} & =\sigma_{l l}-r_{l l} \tag{16}\\ S_{l l} & =\Sigma_{l l}-R_{l l} \tag{17} \end{align*}(16)sll=σll−rll(17)Sll=Σll−Rll
The ensemble average of (15), Σ l l = ⟨ σ l l ⟩ V Σ l l = σ l l V Sigma_(ll)=(:sigma_(ll):)_(V)\Sigma_{l l}=\left\langle\sigma_{l l}\right\rangle_{V}Σll=⟨σll⟩V, is the second moment of the mean concentration, or the "annealed second moment" of Le Doussal and Machta [1989], which is used to define "ensemble dispersion coefficients" 1 2 d Σ l l / d t 1 2 d Σ l l / d t (1)/(2)dSigma_(ll)//dt\frac{1}{2} d \Sigma_{l l} / d t12dΣll/dt [Dentz et al., 2000; Dentz and Carrera, 2005]. From (17) and (11) one obtains the relation between Σ l l Σ l l Sigma_(ll)\Sigma_{l l}Σll and the one-particle quantities X l l X l l X_(ll)X_{l l}Xll and Q l l Q l l Q_(ll)Q_{l l}Qll,
(18) Σ l l = S l l ( 0 ) + ⟨ X l l ⟩ X 0 + M l l + Q l l . (18) Σ l l = S l l ( 0 ) + X l l X 0 + M l l + Q l l . {:(18)Sigma_(ll)=S_(ll)(0)+(:X_(ll):)_(X_(0))+M_(ll)+Q_(ll).:}\begin{equation*} \Sigma_{l l}=S_{l l}(0)+\left\langle X_{l l}\right\rangle_{X_{0}}+M_{l l}+Q_{l l} . \tag{18} \end{equation*}(18)Σll=Sll(0)+⟨Xll⟩X0+Mll+Qll.
With the assumptions from above on the velocity field, (18) reduces to
(19) Σ l l = S l l ( 0 ) + X l l (19) Σ l l = S l l ( 0 ) + X l l {:(19)Sigma_(ll)=S_(ll)(0)+X_(ll):}\begin{equation*} \Sigma_{l l}=S_{l l}(0)+X_{l l} \tag{19} \end{equation*}(19)Σll=Sll(0)+Xll
and using (17) one retrieves (14). For sufficiently large plumes, x ~ l l ≈ X l l x ~ l l ≈ X l l widetilde(x)_(ll)~~X_(ll)\widetilde{x}_{l l} \approx X_{l l}x~ll≈Xll and from (13) and (15) it follows that
(20) σ l l ≈ s l l ≈ S l l ( 0 ) + X l l + m l l . (20) σ l l ≈ s l l ≈ S l l ( 0 ) + X l l + m l l . {:(20)sigma_(ll)~~s_(ll)~~S_(ll)(0)+X_(ll)+m_(ll).:}\begin{equation*} \sigma_{l l} \approx s_{l l} \approx S_{l l}(0)+X_{l l}+m_{l l} . \tag{20} \end{equation*}(20)σll≈sll≈Sll(0)+Xll+mll.
Here is the place to observe that (16) and (17) are identities which can be simply derived by adding and subtracting the mean center of mass ⟨ X l ⟩ D X 0 V X l D X 0 V (:X_(l):)_(DX_(0)V)\left\langle X_{l}\right\rangle_{D X_{0} V}⟨Xl⟩DX0V in (1). They are also strictly equivalent to (4) and (7), respectively, and do not require any conditions for permutation of averages and assumptions on the random velocity field. Such identities naturally occur whenever at least two different stochastic averages are involved [Monin and Yaglom 1971, p. 594; Kitanidis, 1988; Le Doussal and Machta, 1989; Naff et al., 1998; Dentz and Carrera, 2005; Suciu et al., 2006a]. It is important to avoid confusion between the identity (17) and the relation (14) which is based on strong assumptions on the random velocity field. This happens, for instance, when these relations are cited together as "the fundamental relation of Kitanidis [1988] and Dagan [1990]".
We shall use the relations (16-17) in Section 2 to compute expectations and standard deviations which provide a numerical evidence of memory effects and of their relation with the ergodicity issue. Section 3 shortly introduces the Itô formalism used in Section 4 to relate the dispersion quantities discussed in Introduction to velocity statistics and in Section 5 to derive some explicit first-order approximations. These results will be utilized in Section 6, together with the results of simulations presented Section 2, to get quantitative estimates of memory effects for large sources. Section 7 presents some conclusive remarks on the relevance of memory effects for theoretical investigations and practical applications.

2 Nonergodic Behavior of Large Plumes and Memory Effects

That solute plumes originating from large sources can have nonergodic behavior has already been shown by several numerical investigations [Naff et al., 1998; Trefry et al., 2003; Suciu et al., 2006a, Suciu and Vamoş, 2007]. Here, by ergodicity we mean a central issue in stochastic modeling, as formulated for instance by Trefry et al. [2003]: to determine whether theories of solute transport adequately predict the behavior of individual subsurface plumes over time and length scales
Figure 1: Variance of the plume center of mass and corresponding standard deviation (thin lines) for different shapes and extensions of the initial plume.
relevant to practical purposes. The approach to a quasi-ergodic state can be quantified by an "ergodicity range", defined as root mean square deviation of observable quantities from stochastic model predictions [Suciu et al., 2006a]. Since the results presented in the papers quoted above indicate that the ergodicity range does not decrease monotonically and uniformly in time with increasing dimension of initial plumes, strictly speaking, large initial plumes do not ensure ergodic behavior of transport. In the following we present a tableau of two-dimensional GRW numerical simulations, by completing the results for transverse initial slabs [Suciu et al., 2006a] with the cases of longitudinal slabs and square sources of increasing dimensions.
We considered an isotropic two-dimensional aquifer system, characterized by log-hydraulic conductivity with small variance equal to 0.1 and exponentially decaying isotropic correlation with correlation length λ = 1 m λ = 1 m lambda=1m\lambda=1 \mathrm{~m}λ=1 m. Incompressible velocity fields were approximated numerically by Gaussian fields. For fixed mean flow velocity U = 1 m / d U = 1 m / d U=1m//dU=1 \mathrm{~m} / \mathrm{d}U=1 m/d and isotropic local dispersion with constant coefficient D = 0.01 m 2 / d D = 0.01 m 2 / d D=0.01m^(2)//dD=0.01 \mathrm{~m}^{2} / \mathrm{d}D=0.01 m2/d, the Péclet number got a typical value P e = U λ / D = 100 P e = U λ / D = 100 Pe=U lambda//D=100P e=U \lambda / D=100Pe=Uλ/D=100. We obtained ensembles of transport simulations with the GRW algorithm, namely by tracking simultaneously in every velocity realization 10 10 10 10 10^(10)10^{10}1010 particles that were initially uniformly distributed in rectangular domains L 1 λ × L 2 λ L 1 λ × L 2 λ L_(1)lambda xxL_(2)lambdaL_{1} \lambda \times L_{2} \lambdaL1λ×L2λ or released from the origin of the computational grid (for details on the implementation of the numerical method see [Suciu et al., 2006a, Appendix A]). We simulated 1024 realizations for all cases investigated here and we verified that the statistical oscillations for both expectations and standard deviations were bounded uniformly in time by the value D t D t DtD tDt, i.e. by half the local dispersion (see Appendix A). For the purpose of the comparisons presented in Section 6 below, statistical estimations for P e = ∞ P e = ∞ Pe=ooP e=\inftyPe=∞ were also done by simulating advective transport for the same ensemble of velocity realizations. To avoid trapping phenomena [Suciu et al., 2006a, Appendix B1], these simulations were limited to a period equal to 100 U t / λ 100 U t / λ 100 Ut//lambda100 U t / \lambda100Ut/λ.
The accurate GRW simulations on which we base our analyzes can be considered as an ideal experiment in which the concentrations are uniformly sampled throughout the plumes, over large distances of 2000 correlation lengths λ λ lambda\lambdaλ in the mean flow direction, and for initial concentrations uniformly distributed in domains of different shapes and orientations with dimensions up to 100 λ 100 λ 100 lambda100 \lambda100λ. The moments related by (17), S l l , Σ l l S l l , Σ l l S_(ll),Sigma_(ll)S_{l l}, \Sigma_{l l}Sll,Σll, and R l l , l = 1 , 2 R l l , l = 1 , 2 R_(ll),l=1,2R_{l l}, l=1,2Rll,l=1,2, were computed respectively by averaging
Figure 2: Ensemble average second moment of the actual concentration for different shapes and extensions of the initial plume.
over realizations the second moments centered at the center of mass of the actual plumes, s l l s l l s_(ll)s_{l l}sll (1), the second moments with respect to the ensemble average center of mass, σ l l ( 15 ) σ l l ( 15 ) sigma_(ll)(15)\sigma_{l l}(15)σll(15), and the squared differences between the centers of mass of actual plumes and the ensemble average center of mass, r l l ( 6 ) r l l ( 6 ) r_(ll)(6)r_{l l}(6)rll(6). The standard deviations ( S D S D SDS DSD ) were also computed for all quantities and in all cases investigated.
Figure 1 shows that the dispersion of the center of mass, R l l R l l R_(ll)R_{l l}Rll, decays monotonically with the increasing dimensions of the source, irrespective of its shape and orientation. The standard deviations S D [ r l l ] S D r l l SD[r_(ll)]S D\left[r_{l l}\right]SD[rll] were found to be of the same order of magnitude. They are significantly larger than R l l R l l R_(ll)R_{l l}Rll only for small sources. This provides a numerical support for the ergodic arguments which suggest that for large plumes R l l ≈ 0 R l l ≈ 0 R_(ll)~~0R_{l l} \approx 0Rll≈0 and, according to ( 17 ) , S l l ≈ Σ l l ( 17 ) , S l l ≈ Σ l l (17),S_(ll)~~Sigma_(ll)(17), S_{l l} \approx \Sigma_{l l}(17),Sll≈Σll, i.e the expectation of the second moment of the actual concentration can be approximated by the second moment of the ensemble averaged concentration. Another confirmation is given by the visual comparison between S l l − S l l ( 0 ) S l l − S l l ( 0 ) S_(ll)-S_(ll)(0)S_{l l}-S_{l l}(0)Sll−Sll(0) (Figure 2) and Σ l l − S l l ( 0 ) Σ l l − S l l ( 0 ) Sigma_(ll)-S_(ll)(0)\Sigma_{l l}-S_{l l}(0)Σll−Sll(0) (Figure 3), which shows that for large plumes (i.e. L ≥ 50 L ≥ 50 L >= 50L \geq 50L≥50 ) the two moments behave closely at all times. The case of longitudinal dispersion for longitudinal slabs, which shows slower decrease of R 11 R 11 R_(11)R_{11}R11 with increasing L L LLL, indicates that an accurate estimation of the mean dispersion by Σ l l Σ l l Sigma_(ll)\Sigma_{l l}Σll requires transverse dimensions of the initial plume larger than λ λ lambda\lambdaλ.
Ergodicity with respect to the upscaled macrodispersion process (Gaussian with constant coefficients D l l ∗ ) D l l ∗ {:D_(ll)^(**))\left.D_{l l}^{*}\right)Dll∗) is quantified by the deviation of the mean moments from the macrodispersion moments, S l l − S l l ( 0 ) − 2 D l l ∗ t = η 1 S l l − S l l ( 0 ) − 2 D l l ∗ t = η 1 S_(ll)-S_(ll)(0)-2D_(ll)^(**)t=eta_(1)S_{l l}-S_{l l}(0)-2 D_{l l}^{*} t=\eta_{1}Sll−Sll(0)−2Dll∗t=η1, and by their standard deviation S D [ s l l ] = η 2 S D s l l = η 2 SD[s_(ll)]=eta_(2)S D\left[s_{l l}\right]=\eta_{2}SD[sll]=η2, which define the ergodicity range η [ s l l − S l l ( 0 ) ] = [ η 1 2 + η 2 2 ] 1 / 2 η s l l − S l l ( 0 ) = η 1 2 + η 2 2 1 / 2 eta[s_(ll)-S_(ll)(0)]=[eta_(1)^(2)+eta_(2)^(2)]^(1//2)\eta\left[s_{l l}-S_{l l}(0)\right]=\left[\eta_{1}^{2}+\eta_{2}^{2}\right]^{1 / 2}η[sll−Sll(0)]=[η12+η22]1/2 [Suciu et al., 2006a]. The macrodispersion coefficients for the transport problem considered here have the values D 11 ∗ = 0.11 m 2 / d D 11 ∗ = 0.11 m 2 / d D_(11)^(**)=0.11m^(2)//dD_{11}^{*}=0.11 \mathrm{~m}^{2} / \mathrm{d}D11∗=0.11 m2/d and D 22 ∗ = 0.01 m 2 / d D 22 ∗ = 0.01 m 2 / d D_(22)^(**)=0.01m^(2)//dD_{22}^{*}=0.01 \mathrm{~m}^{2} / \mathrm{d}D22∗=0.01 m2/d. The ergodicity ranges η η eta\etaη, normalized with the moments of the macrodispersion process 2 D l l ∗ t 2 D l l ∗ t 2D_(ll)^(**)t2 D_{l l}^{*} t2Dll∗t, computed for all the initial conditions considered in our GRW simulations are presented in Figure 4. Contradicting the common belief that large plumes have ergodic behavior, the ergodicity range does not decrease monotonically and uniformly in time with increasing source dimension L L LLL. On the contrary, η η eta\etaη may considerably increase with L L LLL, may show non-monotonous time behavior and depends strongly on the shape and orientation of the source. In spite of its intricate behavior, the general decreasing trend of η η eta\etaη indicates the mean-square convergence [ s l l ( t ) − S l l ( 0 ) ] / ( 2 t ) → D l l ∗ s l l ( t ) − S l l ( 0 ) / ( 2 t ) → D l l ∗ [s_(ll)(t)-S_(ll)(0)]//(2t)rarrD_(ll)^(**)\left[s_{l l}(t)-S_{l l}(0)\right] /(2 t) \rightarrow D_{l l}^{*}[sll(t)−Sll(0)]/(2t)→Dll∗ for t → ∞ t → ∞ t rarr oot \rightarrow \inftyt→∞,
Figure 3: Second moment of the ensemble average concentration for different shapes and extensions of the initial plume.
i.e. the "asymptotic ergodicity" with respect to macrodispersion moments.
At finite times of interest for predictions of groundwater contamination the macrodispersion model is useless, as shown by the large values of the ergodicity range η η eta\etaη in Figure 4 . It is therefore desirable to have maximum estimations of dispersion through second moments independent of the initial conditions and typical of a given aquifer system. Figure 3 shows that reliable memory-free moments can be defined by Σ l l − S l l ( 0 ) Σ l l − S l l ( 0 ) Sigma_(ll)-S_(ll)(0)\Sigma_{l l}-S_{l l}(0)Σll−Sll(0) for slab sources perpendicular to l l lll-axis: longitudinal moment for transverse sources, and transverse moment for longitudinal sources. Thus, we can conclude that for such situations the second moment of the ensemble mean concentration approximates the oneparticle dispersion, Σ l l − S l l ( 0 ) ≈ X l l Σ l l − S l l ( 0 ) ≈ X l l Sigma_(ll)-S_(ll)(0)~~X_(ll)\Sigma_{l l}-S_{l l}(0) \approx X_{l l}Σll−Sll(0)≈Xll, in agreement with (19). However, for sources with large extensions on the l l lll-axis, Σ l l Σ l l Sigma_(ll)\Sigma_{l l}Σll "remembers" the shape and the orientation of the initial plume. As indicated by Figure 3, this memory effect is enhanced by the dimension of the source. Since the spatial variance of the one-particle center of mass, Q l l Q l l Q_(ll)Q_{l l}Qll, was found to be much smaller than the local dispersion [Suciu and Vamoş, 2007, Suciu et al., 2008], the memory effect shown by Σ l l Σ l l Sigma_(ll)\Sigma_{l l}Σll is mainly caused by non-vanishing mean memory term M l l M l l M_(ll)M_{l l}Mll, in keeping with (18). For large sources, the expected second moment S l l S l l S_(ll)S_{l l}Sll shown in Figure 2 follows closely the behavior of Σ l l Σ l l Sigma_(ll)\Sigma_{l l}Σll and undergoes the same memory effect. This is an unexpected result, since theoretical considerations presented in Introduction indicate that the ensemble average memory term should vanish irrespective of initial conditions. A discussion on the discrepancy between simulations and theory is deferred to Section 6 below.
The memory effects are much more significant for single realizations. In Figure 5 one can see that for large extensions of the source on the l l lll-axis the standard deviations S D [ σ l l ] S D σ l l SD[sigma_(ll)]S D\left[\sigma_{l l}\right]SD[σll] are one or two orders of magnitude larger than for a point source. S D [ s l l ] S D s l l SD[s_(ll)]S D\left[s_{l l}\right]SD[sll] is smaller than S D [ σ l l ] S D σ l l SD[sigma_(ll)]S D\left[\sigma_{l l}\right]SD[σll] in cases with noticeable randomness of the center of mass (shown in Figure 1). For large sources ( L ≥ 50 L ≥ 50 L >= 50L \geq 50L≥50 ), regardless their shape and orientation, S D [ σ l l ] S D σ l l SD[sigma_(ll)]S D\left[\sigma_{l l}\right]SD[σll] practically coincide with S D [ s l l ] S D s l l SD[s_(ll)]S D\left[s_{l l}\right]SD[sll]. We have thus a numerical validation of the relation (20) which, for sufficiently large sources, allows the estimation of the standard deviation of either σ l l σ l l sigma_(ll)\sigma_{l l}σll or s l l s l l s_(ll)s_{l l}sll by the standard deviation of the memory terms m l l m l l m_(ll)m_{l l}mll. Since the mean memory effects indicated by Figure 3 are much smaller than the corresponding standard deviations in Figure 5, from (20) we conclude that S D [ m l l ] S D m l l SD[m_(ll)]S D\left[m_{l l}\right]SD[mll] is also a measure of the root
Figure 4: Ergodicity ranges η η eta\etaη with respect to macrodispersion moments of the second moment of the actual plume for different shapes and extensions of the initial plume.
mean square deviation of s l l − S l l ( 0 ) s l l − S l l ( 0 ) s_(ll)-S_(ll)(0)s_{l l}-S_{l l}(0)sll−Sll(0) from X l l X l l X_(ll)X_{l l}Xll, i.e. of the ergodicity range of the actual second moment with respect to the one-particle dispersion.
The unexpected behavior for large transverse and square sources of ergodicity ranges, with respect to both the macrodispersion (Figure 4) and one-particle dispersion (Figure 5), deserves a particular attention. Increasing the transverse dimension of the source delays the monotonic decrease of ergodicity range for the longitudinal dispersion and increases the late time ergodicity range of the transverse dispersion. This feature has already been reported in [Suciu et al., 2006a], where it was also shown that the cross-section space average concentration at the plume center of mass behaves similarly to the longitudinal dispersion and it was conjectured that the timescale of asymptotic ergodicity increases with the square of transverse dimension of the source. As for the pre-asymptotic regime, we remark two opposite effects of the source dimensions. While for large dimensions in both longitudinal and transverse directions the randomness of the plume center of mass, described by R l l R l l R_(ll)R_{l l}Rll and S D [ r l l ] S D r l l SD[r_(ll)]S D\left[r_{l l}\right]SD[rll], becomes negligible (Figure 1), the standard deviation of the actual dispersion S D [ s l l ] S D s l l SD[s_(ll)]S D\left[s_{l l}\right]SD[sll] increases for sources with large extension on l l lll-axis (Figure 5). Ergodic behavior with respect to one-particle dispersion can be only expected for narrow sources transverse to l l lll.

3 Itô Formalism

When local dispersion is described as a diffusion-like process, the Itô equation provides mathematically consistent and physically relevant descriptions of the transport in random media. For instance, the validity of the Itô equation as a model for the motion of the solute molecules at Darcy scale was inferred from a microscopic description of the transport in saturated porous media [Bhattacharya and Gupta, 1983]. For a given realization V ( x ) V ( x ) V(x)\mathbf{V}(\mathbf{x})V(x) of a time independent random velocity field and a constant local dispersion coefficient D D DDD, the trajectories of the local advection-dispersion process starting at t = 0 t = 0 t=0t=0t=0 from X 0 X 0 X_(0)\mathbf{X}_{0}X0 are the solutions of the integral Itô equation
(21) X l ( t ) = X 0 l + ∫ 0 t V l ( X ( t ′ ) ) d t ′ + W l ( t ) (21) X l ( t ) = X 0 l + ∫ 0 t   V l X t ′ d t ′ + W l ( t ) {:(21)X_(l)(t)=X_(0l)+int_(0)^(t)V_(l)(X(t^(')))dt^(')+W_(l)(t):}\begin{equation*} X_{l}(t)=X_{0 l}+\int_{0}^{t} V_{l}\left(\mathbf{X}\left(t^{\prime}\right)\right) d t^{\prime}+W_{l}(t) \tag{21} \end{equation*}(21)Xl(t)=X0l+∫0tVl(X(t′))dt′+Wl(t)
Figure 5: Standard deviation of the second moment with respect to ensemble averaged center of mass and standard deviation of the second moment of the actual plume (thin lines) for different shapes and extensions of the initial plume.
where W l ( t ) = ∫ 0 t d W l ( t ′ ) W l ( t ) = ∫ 0 t   d W l t ′ W_(l)(t)=int_(0)^(t)dW_(l)(t^('))W_{l}(t)=\int_{0}^{t} d W_{l}\left(t^{\prime}\right)Wl(t)=∫0tdWl(t′) is a stochastic integral which defines a Wiener process with the properties
(22) ⟨ W l ⟩ D = 0 , ⟨ W l 2 ⟩ D = 2 D t . (22) W l D = 0 , W l 2 D = 2 D t . {:(22)(:W_(l):)_(D)=0","(:W_(l)^(2):)_(D)=2Dt.:}\begin{equation*} \left\langle W_{l}\right\rangle_{D}=0,\left\langle W_{l}^{2}\right\rangle_{D}=2 D t . \tag{22} \end{equation*}(22)⟨Wl⟩D=0,⟨Wl2⟩D=2Dt.
The particular construction of the partial sums, which evaluates the integrand at the left of the discretization interval, and the independence of the increments d W l d W l dW_(l)d W_{l}dWl of the Wiener process lead to the cancelation of the average of the Itô integral [Gardiner, 1985, Kloeden and Platen, 1995]. For instance,
(23) ⟨ ∫ 0 t X l ( t ′ ) d W l ( t ′ ) ⟩ D = 0 . (23) ∫ 0 t   X l t ′ d W l t ′ D = 0 . {:(23)(:int_(0)^(t)X_(l)(t^('))dW_(l)(t^(')):)_(D)=0.:}\begin{equation*} \left\langle\int_{0}^{t} X_{l}\left(t^{\prime}\right) d W_{l}\left(t^{\prime}\right)\right\rangle_{D}=0 . \tag{23} \end{equation*}(23)⟨∫0tXl(t′)dWl(t′)⟩D=0.
Taking into account that ⟨ d W l ⟩ D 2 = 2 D d t d W l D 2 = 2 D d t (:dW_(l):)_(D)^(2)=2Ddt\left\langle d W_{l}\right\rangle_{D}^{2}=2 D d t⟨dWl⟩D2=2Ddt, as follows from (22), and collecting terms of first-order in d t d t dtd tdt and d W d W dWd WdW, the differential of a function of time and of the l l lll component of the Itô process (21), f ( t , x ) = f ( t , X l ( t ) ) f ( t , x ) = f t , X l ( t ) f(t,x)=f(t,X_(l)(t))f(t, x)=f\left(t, X_{l}(t)\right)f(t,x)=f(t,Xl(t)), is computed as
d f = ∂ f ∂ t d t + V l ∂ f ∂ x d t + D d 2 f ∂ x 2 d t + ∂ f ∂ x d W , d f = ∂ f ∂ t d t + V l ∂ f ∂ x d t + D d 2 f ∂ x 2 d t + ∂ f ∂ x d W , df=(del f)/(del t)dt+V_(l)(del f)/(del x)dt+D(d^(2)f)/(delx^(2))dt+(del f)/(del x)dW,d f=\frac{\partial f}{\partial t} d t+V_{l} \frac{\partial f}{\partial x} d t+D \frac{d^{2} f}{\partial x^{2}} d t+\frac{\partial f}{\partial x} d W,df=∂f∂tdt+Vl∂f∂xdt+Dd2f∂x2dt+∂f∂xdW,
where the derivatives are evaluated at x = X l ( t ) x = X l ( t ) x=X_(l)(t)x=X_{l}(t)x=Xl(t). This is a particular form of the stochastic chain rule known as the Itô formula [e.g. Kloeden and Platen, 1995, Section 3.3], which is mainly used in its integral form
f ( t , X l ( t ) ) − f ( 0 , X l ( 0 ) ) = ∫ 0 t ∂ f ( t ′ , X l ( t ′ ) ) ∂ t ′ d t ′ + ∫ 0 t [ V l ( X l ( t ′ ) ) ∂ f ( t ′ , X l ( t ′ ) ) d x + D ∂ 2 f ( t ′ , X l ( t ′ ) ) ∂ x 2 ] d t ′ (24) + ∫ 0 t ∂ f ( t ′ , X l ( t ′ ) ) ∂ x d W ( t ′ ) f t , X l ( t ) − f 0 , X l ( 0 ) = ∫ 0 t   ∂ f t ′ , X l t ′ ∂ t ′ d t ′ + ∫ 0 t   V l X l t ′ ∂ f t ′ , X l t ′ d x + D ∂ 2 f t ′ , X l t ′ ∂ x 2 d t ′ (24) + ∫ 0 t   ∂ f t ′ , X l t ′ ∂ x d W t ′ {:[f(t,X_(l)(t))-f(0,X_(l)(0))=int_(0)^(t)(del f(t^('),X_(l)(t^('))))/(delt^('))dt^(')],[+int_(0)^(t)[V_(l)(X_(l)(t^(')))(del f(t^('),X_(l)(t^('))))/(dx)+D(del^(2)f(t^('),X_(l)(t^('))))/(delx^(2))]dt^(')],[(24)+int_(0)^(t)(del f(t^('),X_(l)(t^('))))/(del x)dW(t^('))]:}\begin{align*} & f\left(t, X_{l}(t)\right)-f\left(0, X_{l}(0)\right)=\int_{0}^{t} \frac{\partial f\left(t^{\prime}, X_{l}\left(t^{\prime}\right)\right)}{\partial t^{\prime}} d t^{\prime} \\ & +\int_{0}^{t}\left[V_{l}\left(X_{l}\left(t^{\prime}\right)\right) \frac{\partial f\left(t^{\prime}, X_{l}\left(t^{\prime}\right)\right)}{d x}+D \frac{\partial^{2} f\left(t^{\prime}, X_{l}\left(t^{\prime}\right)\right)}{\partial x^{2}}\right] d t^{\prime} \\ & +\int_{0}^{t} \frac{\partial f\left(t^{\prime}, X_{l}\left(t^{\prime}\right)\right)}{\partial x} d W\left(t^{\prime}\right) \tag{24} \end{align*}f(t,Xl(t))−f(0,Xl(0))=∫0t∂f(t′,Xl(t′))∂t′dt′+∫0t[Vl(Xl(t′))∂f(t′,Xl(t′))dx+D∂2f(t′,Xl(t′))∂x2]dt′(24)+∫0t∂f(t′,Xl(t′))∂xdW(t′)
The joint n n nnn-times probability densities [van Kampen, 1981] of the Itô process (21),
(25) p ( x 1 , t 1 ; x 2 , t 2 ; ⋯ ; x n , t n ) = ⟨ δ [ x 1 − X ( t 1 ) ] δ [ x 2 − X ( t 2 ) ] ⋯ δ [ x n − X ( t n ) ] ⟩ D X 0 (25) p x 1 , t 1 ; x 2 , t 2 ; ⋯ ; x n , t n = δ x 1 − X t 1 δ x 2 − X t 2 ⋯ δ x n − X t n D X 0 {:[(25)p(x_(1),t_(1);x_(2),t_(2);cdots;x_(n),t_(n))=],[(:delta[x_(1)-X(t_(1))]delta[x_(2)-X(t_(2))]cdots delta[x_(n)-X(t_(n))]:)_(DX_(0))]:}\begin{align*} & p\left(\mathbf{x}_{1}, t_{1} ; \mathbf{x}_{2}, t_{2} ; \cdots ; \mathbf{x}_{n}, t_{n}\right)= \tag{25}\\ & \left\langle\delta\left[\mathbf{x}_{1}-\mathbf{X}\left(t_{1}\right)\right] \delta\left[\mathbf{x}_{2}-\mathbf{X}\left(t_{2}\right)\right] \cdots \delta\left[\mathbf{x}_{n}-\mathbf{X}\left(t_{n}\right)\right]\right\rangle_{D X_{0}} \end{align*}(25)p(x1,t1;x2,t2;⋯;xn,tn)=⟨δ[x1−X(t1)]δ[x2−X(t2)]⋯δ[xn−X(tn)]⟩DX0
are projections of the probability measure of the process starting from initial positions X 0 X 0 X_(0)\mathbf{X}_{\mathbf{0}}X0 on the cylindrical sets { X ( t 1 ) = x 1 ; X ( t 2 ) = x 2 ; ⋯ X ( t n ) = x n } X t 1 = x 1 ; X t 2 = x 2 ; ⋯ X t n = x n {X(t_(1))=x_(1);X(t_(2))=x_(2);cdotsX(t_(n))=x_(n)}\left\{\mathbf{X}\left(t_{1}\right)=\mathbf{x}_{1} ; \mathbf{X}\left(t_{2}\right)=\mathbf{x}_{2} ; \cdots \mathbf{X}\left(t_{n}\right)=\mathbf{x}_{n}\right\}{X(t1)=x1;X(t2)=x2;⋯X(tn)=xn} [see e.g. Doob, 1953; Zwanzig, 1961; Gardiner, 1985; Kloeden and Platen, 1995]. In particular, the normalized concentration is defined by the one-time probability, c ( x , t ) = p ( x , t ) c ( x , t ) = p ( x , t ) c(x,t)=p(x,t)c(\mathbf{x}, t)=p(\mathbf{x}, t)c(x,t)=p(x,t).
Doob [1953, chap. VI] has proved that with the Itô interpretation of the stochastic integral the description in terms of trajectories (21) is equivalent to the description given by the Fokker-Planck equation for the transition probabilities of the process. The normalized concentration c c ccc verifies the same Fokker-Planck equation. In the case of constant coefficient D D DDD which we consider, this is simply the advection-dispersion equation
(26) ∂ t c + ∇ ( V c ) = D ∇ 2 c (26) ∂ t c + ∇ ( V c ) = D ∇ 2 c {:(26)del_(t)c+grad(Vc)=Dgrad^(2)c:}\begin{equation*} \partial_{t} c+\nabla(\mathbf{V} c)=D \nabla^{2} c \tag{26} \end{equation*}(26)∂tc+∇(Vc)=D∇2c
The connection with the Stratonovich equation and Fick's law, for smoothly variable D D DDD, can be established by using the Itô formula [Gardiner, 1985, section 4.3.6]. Itô formalism can also be used to account for variable porosity and discontinuous D D DDD [LaBolle et al., 2000]. The equivalence between (21) and (26) is the basis for numerical methods such as the particle tracking procedure in its Itô implementation, i.e. the forward-time Euler scheme for equation (21) [Kloeden and Platen, 1995; Lichtner et al., 2002], or the GRW algorithm [Vamoş et al., 2003] used in simulations presented in Section 2.
We note here the difference between the Itô formalism and a somewhat similar "Lagrangian" approach. In the latter, the concentration is defined by the relation (25) for n = 1 n = 1 n=1n=1n=1, without average over the realizations of the local dispersion process. This Lagrangian concentration is the solution of a purely advective equation, like that obtained by removing the dispersive term on the right hand side of (26) [Gardiner, 1985, p. 54]. When the Lagrangian approach is extended to advection-dispersion processes, the velocity field is of the form V + d W / d t V + d W / d t V+dW//dt\mathbf{V}+d \mathbf{W} / d tV+dW/dt [see e.g. Dagan, 1987; Fiori and Dagan, 2000]. Since the derivative of the Wiener process does not exist [Gardiner, 1985, p. 69], the transformation from X 0 X 0 X_(0)\mathbf{X}_{\mathbf{0}}X0 to X X X\mathbf{X}X is not one-to-one and has no unique inverse. Or, these are necessary conditions for the definition of the "fluid particle" and for the rigorous Lagrangian approach [Lundgren, 1981]. However, the smoothness and boundedness of the time derivative of the trajectory X l ( t ) X l ( t ) X_(l)(t)X_{l}(t)Xl(t) which permit the use of a Lagrangian description [Sposito and Dagan, 1994] can be ensured by modeling the local dispersion as a spatially correlated noise [Zirbel, 2001, Section 7]. Alternatively, a Fokker-Planck equation for the probability density of the particle position can be derived by applying a rigorous asymptotic analysis for small velocity fluctuations to a Lagrangian trajectory equation [Brouwers, 2002]. Even if the conceptual inconsistencies can be avoided in this way and equation (21) can be retrieved as a limit case, the Lagrangian approach still differs from Itô formalism by two important features. Firstly, joint probabilities of advective and dispersive displacements have to be inferred from hypotheses and approximations [Fiori and Dagan, 2000]. Secondly, by using a Fourier representation for trajectories, the dependence of the second moment on the local dispersion coefficient is introduced via the relation between the characteristic function and the variance of the Gaussian local dispersion process [Dagan, 1987, equation (3.11); Fiori and Dagan, 2000, equation (14)]. For our purposes, the Itô formalism has the advantage of simplicity and clarity. It relates dispersion and velocity statistics and yields first-order approximations in real time-space domain, by using explicit averaging over local dispersion, initial concentration distribution, and velocity realizations.

4 Second Moments and Memory Terms

Assuming the existence and the uniqueness of the solutions of the Itô equation, the terms of the actual dispersion (4) can be computed as follows. From (21) and (22) one obtains the average of
the relative displacement X ~ l = X l − X 0 l X ~ l = X l − X 0 l widetilde(X)_(l)=X_(l)-X_(0l)\widetilde{X}_{l}=X_{l}-X_{0 l}X~l=Xl−X0l
(27) ⟨ X ~ l ⟩ D X 0 V ( t ) = ∫ 0 t ⟨ V l ( X ( t ′ ) ) ⟩ D X 0 V d t ′ (27) X ~ l D X 0 V ( t ) = ∫ 0 t   V l X t ′ D X 0 V d t ′ {:(27)(: widetilde(X)_(l):)_(DX_(0)V)(t)=int_(0)^(t)(:V_(l)(X(t^('))):)_(DX_(0)V)dt^('):}\begin{equation*} \left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}(t)=\int_{0}^{t}\left\langle V_{l}\left(\mathbf{X}\left(t^{\prime}\right)\right)\right\rangle_{D X_{0} V} d t^{\prime} \tag{27} \end{equation*}(27)⟨X~l⟩DX0V(t)=∫0t⟨Vl(X(t′))⟩DX0Vdt′
By applying the Itô formula (24) to the function f ( t , X l ( t ) ) = [ X ~ l ( t ) − ⟨ X ~ l ⟩ D X 0 V ( t ) ] 2 , f ( 0 , X l ( 0 ) ) = f t , X l ( t ) = X ~ l ( t ) − X ~ l D X 0 V ( t ) 2 , f 0 , X l ( 0 ) = f(t,X_(l)(t))=[ widetilde(X)_(l)(t)-(: widetilde(X)_(l):)_(DX_(0)V)(t)]^(2),f(0,X_(l)(0))=f\left(t, X_{l}(t)\right)=\left[\widetilde{X}_{l}(t)-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}(t)\right]^{2}, f\left(0, X_{l}(0)\right)=f(t,Xl(t))=[X~l(t)−⟨X~l⟩DX0V(t)]2,f(0,Xl(0))= 0 , one obtains
[ X ~ l ( t ) − ⟨ X ~ l ⟩ D X 0 V ( t ) ] 2 = − 2 ∫ 0 t [ X ~ l ( t ′ ) − ⟨ X ~ l ⟩ D X 0 V ( t ′ ) ] ⟨ V l ( X ( t ′ ) ) ⟩ D X 0 V d t ′ + 2 ∫ 0 t [ X ~ l ( t ′ ) − ⟨ X ~ l ⟩ D X 0 V ( t ′ ) ] V l ( X l ( t ′ ) ) d t ′ + 2 D t (28) + 2 ∫ 0 t [ X ~ l ( t ′ ) − ⟨ X ~ l ⟩ D X 0 V ( t ′ ) ] d W ( t ′ ) X ~ l ( t ) − X ~ l D X 0 V ( t ) 2 = − 2 ∫ 0 t   X ~ l t ′ − X ~ l D X 0 V t ′ V l X t ′ D X 0 V d t ′ + 2 ∫ 0 t   X ~ l t ′ − X ~ l D X 0 V t ′ V l X l t ′ d t ′ + 2 D t (28) + 2 ∫ 0 t   X ~ l t ′ − X ~ l D X 0 V t ′ d W t ′ {:[[ widetilde(X)_(l)(t)-(: widetilde(X)_(l):)_(DX_(0)V)(t)]^(2)=],[-2int_(0)^(t)[ widetilde(X)_(l)(t^('))-(: widetilde(X)_(l):)_(DX_(0)V)(t^('))](:V_(l)(X(t^('))):)_(DX_(0)V)dt^(')],[+2int_(0)^(t)[ widetilde(X)_(l)(t^('))-(: widetilde(X)_(l):)_(DX_(0)V)(t^('))]V_(l)(X_(l)(t^(')))dt^(')+2Dt],[(28)+2int_(0)^(t)[ widetilde(X)_(l)(t^('))-(: widetilde(X)_(l):)_(DX_(0)V)(t^('))]dW(t^('))]:}\begin{align*} & {\left[\widetilde{X}_{l}(t)-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}(t)\right]^{2}=} \\ & -2 \int_{0}^{t}\left[\widetilde{X}_{l}\left(t^{\prime}\right)-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\left(t^{\prime}\right)\right]\left\langle V_{l}\left(\mathbf{X}\left(t^{\prime}\right)\right)\right\rangle_{D X_{0} V} d t^{\prime} \\ & +2 \int_{0}^{t}\left[\widetilde{X}_{l}\left(t^{\prime}\right)-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\left(t^{\prime}\right)\right] V_{l}\left(X_{l}\left(t^{\prime}\right)\right) d t^{\prime}+2 D t \\ & +2 \int_{0}^{t}\left[\widetilde{X}_{l}\left(t^{\prime}\right)-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\left(t^{\prime}\right)\right] d W\left(t^{\prime}\right) \tag{28} \end{align*}[X~l(t)−⟨X~l⟩DX0V(t)]2=−2∫0t[X~l(t′)−⟨X~l⟩DX0V(t′)]⟨Vl(X(t′))⟩DX0Vdt′+2∫0t[X~l(t′)−⟨X~l⟩DX0V(t′)]Vl(Xl(t′))dt′+2Dt(28)+2∫0t[X~l(t′)−⟨X~l⟩DX0V(t′)]dW(t′)
After averaging (28) with respect to D D DDD and X 0 X 0 X_(0)X_{0}X0 and using the properties (22-23) of the Itô integral, one obtains the explicit form of (5),
x ~ l l ( t ) = 2 D t + (29) 2 ∫ 0 t ⟨ [ X ~ l ( t ′ ) − ⟨ X ~ l ⟩ D X 0 V ( t ′ ) ] u l ( X ( t ′ ) , t ′ ) ⟩ D X 0 d t ′ x ~ l l ( t ) = 2 D t + (29) 2 ∫ 0 t   X ~ l t ′ − X ~ l D X 0 V t ′ u l X t ′ , t ′ D X 0 d t ′ {:[ widetilde(x)_(ll)(t)=2Dt+],[(29)2int_(0)^(t)(:[ widetilde(X)_(l)(t^('))-(: widetilde(X)_(l):)_(DX_(0)V)(t^('))]u_(l)(X(t^(')),t^(')):)_(DX_(0))dt^(')]:}\begin{align*} & \widetilde{x}_{l l}(t)=2 D t+ \\ & 2 \int_{0}^{t}\left\langle\left[\widetilde{X}_{l}\left(t^{\prime}\right)-\left\langle\widetilde{X}_{l}\right\rangle_{D X_{0} V}\left(t^{\prime}\right)\right] u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right)\right\rangle_{D X_{0}} d t^{\prime} \tag{29} \end{align*}x~ll(t)=2Dt+(29)2∫0t⟨[X~l(t′)−⟨X~l⟩DX0V(t′)]ul(X(t′),t′)⟩DX0dt′
where u l ( X ( t ) , t ) = V l ( X ( t ) ) − ⟨ V l ( X ( t ) ) ⟩ D X 0 V u l ( X ( t ) , t ) = V l ( X ( t ) ) − V l ( X ( t ) ) D X 0 V u_(l)(X(t),t)=V_(l)(X(t))-(:V_(l)(X(t)):)_(DX_(0)V)u_{l}(\mathbf{X}(t), t)=V_{l}(\mathbf{X}(t))-\left\langle V_{l}(\mathbf{X}(t))\right\rangle_{D X_{0} V}ul(X(t),t)=Vl(X(t))−⟨Vl(X(t))⟩DX0V. The fluctuations u l u l u_(l)u_{l}ul have a non-autonomous dependence on time because the total average of the velocity sampled on trajectories ⟨ V l ( X ( t ) ) ⟩ D X 0 V V l ( X ( t ) ) D X 0 V (:V_(l)(X(t)):)_(DX_(0)V)\left\langle V_{l}(\mathbf{X}(t))\right\rangle_{D X_{0} V}⟨Vl(X(t))⟩DX0V is, in general, a time function. The memory term (3) and the fluctuation of the center of mass (6) are obtained using (21) and (27) as
(30) m l l ( t ) = 2 ∫ 0 t ⟨ [ X 0 l − ⟨ X 0 l ⟩ X 0 ] u l ( X ( t ′ ) , t ′ ) ⟩ D X 0 d t ′ (31) r l l ( t ) = ∫ 0 t ∫ 0 t ⟨ u l ( X ( t ′ ) , t ′ ) ⟩ D X 0 ⟨ u l ( X ( t ′ ′ ) , t ′ ′ ) ⟩ D X 0 d t ′ d t ′ ′ (30) m l l ( t ) = 2 ∫ 0 t   X 0 l − X 0 l X 0 u l X t ′ , t ′ D X 0 d t ′ (31) r l l ( t ) = ∫ 0 t   ∫ 0 t   u l X t ′ , t ′ D X 0 u l X t ′ ′ , t ′ ′ D X 0 d t ′ d t ′ ′ {:[(30)m_(ll)(t)=2int_(0)^(t)(:[X_(0l)-(:X_(0l):)_(X_(0))]u_(l)(X(t^(')),t^(')):)_(DX_(0))dt^(')],[(31)r_(ll)(t)=int_(0)^(t)int_(0)^(t)(:u_(l)(X(t^(')),t^(')):)_(DX_(0))(:u_(l)(X(t^('')),t^('')):)_(DX_(0))dt^(')dt^('')]:}\begin{gather*} m_{l l}(t)=2 \int_{0}^{t}\left\langle\left[X_{0 l}-\left\langle X_{0 l}\right\rangle_{X_{0}}\right] u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right)\right\rangle_{D X_{0}} d t^{\prime} \tag{30}\\ r_{l l}(t)=\int_{0}^{t} \int_{0}^{t}\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right)\right\rangle_{D X_{0}}\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime \prime}\right), t^{\prime \prime}\right)\right\rangle_{D X_{0}} d t^{\prime} d t^{\prime \prime} \tag{31} \end{gather*}(30)mll(t)=2∫0t⟨[X0l−⟨X0l⟩X0]ul(X(t′),t′)⟩DX0dt′(31)rll(t)=∫0t∫0t⟨ul(X(t′),t′)⟩DX0⟨ul(X(t′′),t′′)⟩DX0dt′dt′′
Now, using the explicit forms (29-31) of the terms composing the actual dispersion (4), we are in a position to discuss the implications of the special properties of the random velocity field. We recall that the Eulerian velocity field was assumed to be statistically homogeneous, incompressible, continuously differentiable, and for the existence of unique pathwise strong solutions of Itô equation (21) it should also satisfy the Lipschitz continuity and linear growth bound conditions [Kloeden and Platen, 1995, Theorem 4.5.3]. Inspired by the terminology of [Port and Stone, 1976] we use the short names "regular fields", to denote this set of conditions, and "ergodic field" when the random velocity has finite correlation range (which implies ergodicity [Yaglom, 1987; Chilès and Delfiner, 1999]).
First, we consider a fixed realization of the Wiener process (21) which allows us to define a random function of space and time called Lagrangian velocity field V L ( X 0 , t ) = V ( X ( t ) ) V L X 0 , t = V ( X ( t ) ) V^(L)(X_(0),t)=V(X(t))\mathbf{V}^{L}\left(\mathbf{X}_{0}, t\right)=\mathbf{V}(\mathbf{X}(t))VL(X0,t)=V(X(t)) [Port and Stone, 1976]. We mentioned in Introduction that for regular fields the one-point Lagrangian and Eulerian statistics are the same. As a consequence, the ensemble average of the Lagrangian velocity is constant and equals the Eulerian mean velocity, ⟨ V L ( X 0 , t ) ⟩ V = ⟨ V ( X 0 ) ⟩ V = U V L X 0 , t V = V X 0 V = U (:V^(L)(X_(0),t):)_(V)=(:V(X_(0)):)_(V)=U\left\langle\mathbf{V}^{L}\left(\mathbf{X}_{0}, t\right)\right\rangle_{V}=\left\langle\mathbf{V}\left(\mathbf{X}_{0}\right)\right\rangle_{V}=\mathbf{U}⟨VL(X0,t)⟩V=⟨V(X0)⟩V=U, and the velocity fluctuation has the mean zero. The independence on time and the statistical homogeneity of the Eulerian velocity also imply the stationarity and the independence on X 0 X 0 X_(0)\mathbf{X}_{0}X0 of the Lagrangian
correlation [Lumley 1962b; Kramer and Majda, 2007, Proposition 4.2]. Thus, the Lagrangian fluctuations of a regular field satisfy
(32) ⟨ u l ( X ( t ) ) ⟩ V = ⟨ V l ( X ( t ) ) − U l ⟩ V = 0 (33) ⟨ u l ( X ( t ′ ) ) u l ( X ( t ′ ′ ) ) ⟩ V = u l l L ( t ′ − t ′ ′ ; D ) (32) u l ( X ( t ) ) V = V l ( X ( t ) ) − U l V = 0 (33) u l X t ′ u l X t ′ ′ V = u l l L t ′ − t ′ ′ ; D {:[(32)(:u_(l)(X(t)):)_(V)=(:V_(l)(X(t))-U_(l):)_(V)=0],[(33)(:u_(l)(X(t^(')))u_(l)(X(t^(''))):)_(V)=u_(ll)^(L)(t^(')-t^('');D)]:}\begin{gather*} \left\langle u_{l}(\mathbf{X}(t))\right\rangle_{V}=\left\langle V_{l}(\mathbf{X}(t))-U_{l}\right\rangle_{V}=0 \tag{32}\\ \left\langle u_{l}\left(\mathbf{X}\left(t^{\prime}\right)\right) u_{l}\left(\mathbf{X}\left(t^{\prime \prime}\right)\right)\right\rangle_{V}=u_{l l}^{L}\left(t^{\prime}-t^{\prime \prime} ; D\right) \tag{33} \end{gather*}(32)⟨ul(X(t))⟩V=⟨Vl(X(t))−Ul⟩V=0(33)⟨ul(X(t′))ul(X(t′′))⟩V=ullL(t′−t′′;D)
The argument D D DDD in the right side of (33) indicates the dependence of the Lagrangian correlation u l l L u l l L u_(ll)^(L)u_{l l}^{L}ullL on the fixed realization of the Wiener process.
Using (21) and (27), (29) can be expressed as
(34) x ~ l l ( t ) = 2 D t + ∫ 0 t ∫ 0 t ⟨ u l ( X ( t ′ ) , t ′ ) u l ( X ( t ′ ′ ) , t ′ ′ ) ⟩ D X 0 d t ′ d t ′ ′ + 2 ∫ 0 t ⟨ W ( t ′ ) u l ( X ( t ′ ) , t ′ ) ⟩ D X 0 d t ′ (34) x ~ l l ( t ) = 2 D t + ∫ 0 t   ∫ 0 t   u l X t ′ , t ′ u l X t ′ ′ , t ′ ′ D X 0 d t ′ d t ′ ′ + 2 ∫ 0 t   W t ′ u l X t ′ , t ′ D X 0 d t ′ {:(34) widetilde(x)_(ll)(t)=2Dt+int_(0)^(t)int_(0)^(t)(:u_(l)(X(t^(')),t^('))u_(l)(X(t^('')),t^('')):)_(DX_(0))dt^(')dt^('')+2int_(0)^(t)(:W(t^('))u_(l)(X(t^(')),t^(')):)_(DX_(0))dt^('):}\begin{equation*} \widetilde{x}_{l l}(t)=2 D t+\int_{0}^{t} \int_{0}^{t}\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right) u_{l}\left(\mathbf{X}\left(t^{\prime \prime}\right), t^{\prime \prime}\right)\right\rangle_{D X_{0}} d t^{\prime} d t^{\prime \prime}+2 \int_{0}^{t}\left\langle W\left(t^{\prime}\right) u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right)\right\rangle_{D X_{0}} d t^{\prime} \tag{34} \end{equation*}(34)x~ll(t)=2Dt+∫0t∫0t⟨ul(X(t′),t′)ul(X(t′′),t′′)⟩DX0dt′dt′′+2∫0t⟨W(t′)ul(X(t′),t′)⟩DX0dt′
Owing to the permutation of averages, (32) implies the cancelation of the ensemble average of the last term of (34) and from (33) it follows that for regular fields X ~ l l = ⟨ x ~ l l ⟩ V X ~ l l = x ~ l l V widetilde(X)_(ll)=(: widetilde(x)_(ll):)_(V)\widetilde{X}_{l l}=\left\langle\widetilde{x}_{l l}\right\rangle_{V}X~ll=⟨x~ll⟩V coincides with the memory-free one-particle dispersion
(35) X l l ( t ) = 2 D t + ∫ 0 t ∫ 0 t ⟨ u l l L ( t ′ − t ′ ′ ; D ) ⟩ D d t ′ d t ′ ′ (35) X l l ( t ) = 2 D t + ∫ 0 t   ∫ 0 t   u l l L t ′ − t ′ ′ ; D D d t ′ d t ′ ′ {:(35)X_(ll)(t)=2Dt+int_(0)^(t)int_(0)^(t)(:u_(ll)^(L)(t^(')-t^('');D):)_(D)dt^(')dt^(''):}\begin{equation*} X_{l l}(t)=2 D t+\int_{0}^{t} \int_{0}^{t}\left\langle u_{l l}^{L}\left(t^{\prime}-t^{\prime \prime} ; D\right)\right\rangle_{D} d t^{\prime} d t^{\prime \prime} \tag{35} \end{equation*}(35)Xll(t)=2Dt+∫0t∫0t⟨ullL(t′−t′′;D)⟩Ddt′dt′′
The equality X ~ l l = X l l X ~ l l = X l l widetilde(X)_(ll)=X_(ll)\widetilde{X}_{l l}=X_{l l}X~ll=Xll and (8) implies Q l l = 0 Q l l = 0 Q_(ll)=0Q_{l l}=0Qll=0. This also results directly from (32), after having used (10) and (21) to express Q l l Q l l Q_(ll)Q_{l l}Qll as a double time integral of the spatial variance of the ensemble mean velocity [see also Suciu and Vamoş, 2007]. Since as a consequence of (32) the ensemble average of the memory term (30) also vanishes, ⟨ m l l ⟩ V = M l l = 0 m l l V = M l l = 0 (:m_(ll):)_(V)=M_(ll)=0\left\langle m_{l l}\right\rangle_{V}=M_{l l}=0⟨mll⟩V=Mll=0, from (11) we get the result that for regular velocity fields relation (14) holds true.
If the regular Eulerian field is ergodic and Gaussian, then the spatial mean and correlation tend in mean square limit to the corresponding ensemble quantities [Yaglom, 1987, p. 234]. The space mean velocity ⟨ u l ( X ( t ) ) ⟩ X 0 u l ( X ( t ) ) X 0 (:u_(l)(X(t)):)_(X_(0))\left\langle u_{l}(\mathbf{X}(t))\right\rangle_{X_{0}}⟨ul(X(t))⟩X0 and correlation ⟨ u l ( X ( t ′ ) ) u l ( X ( t ′ ′ ) ) ⟩ X 0 u l X t ′ u l X t ′ ′ X 0 (:u_(l)(X(t^(')))u_(l)(X(t^(''))):)_(X_(0))\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime}\right)\right) u_{l}\left(\mathbf{X}\left(t^{\prime \prime}\right)\right)\right\rangle_{X_{0}}⟨ul(X(t′))ul(X(t′′))⟩X0 of the Lagrangian field V L V L V^(L)\mathbf{V}^{L}VL also are unbiased estimators of the ensemble averages (32-33). To prove the ergodicity of the Lagrangian velocity field one has yet to show that the standard deviations of the estimates go to zero for large spatial averaging domains. This requires the knowledge of two-points distributions. But the equality of one-dimensional marginal probability distributions is as far as our knowledge on regular fields goes [Zirbel, 2001]. Mathematical proofs that Lagrangian velocity inherits the entire probability distribution or at least the rapid decorrelation of Eulerian velocity (both implying ergodicity), are extremely difficult and can only be done for constraints stronger than assumed here [e.g. Fannjiang and Komorowski, 1999; Komorowski and Papanicolaou, 1997]. Nevertheless, there is a numerical evidence that space averages with respect to the location of the source of the concentration moments resulted from simulations of diffusion in Gaussian fields converge to their ensemble averages [Suciu et al., 2007b]. With the Itô representation of trajectory (21), the convergence of the first two moments implies the convergence of the Lagrangian estimators for mean and correlation. This argument suggests that the ergodicity of the Lagrangian field may be assumed, provided that the regular Eulerian field is ergodic and Gaussian. With this ergodic hypothesis, by interchanging space and ensemble averages, from (31) and (34) one obtains, respectively,
(36) r l l ( t ) = 0 , x ~ l l ( t ) = X l l ( t ) (36) r l l ( t ) = 0 , x ~ l l ( t ) = X l l ( t ) {:(36)r_(ll)(t)=0"," widetilde(x)_(ll)(t)=X_(ll)(t):}\begin{equation*} r_{l l}(t)=0, \widetilde{x}_{l l}(t)=X_{l l}(t) \tag{36} \end{equation*}(36)rll(t)=0,x~ll(t)=Xll(t)
which, as seen in Introduction, imply (12), (13), and (20). While the assumption of regular field seems to be a necessary ingredient for the ergodicity of the Lagrangian field, the analysis of numerical simulations, presented in Section 6 below, indicates that for Gaussian ergodic fields the behavior described by relations (36) can occur for large sources even if regularity conditions are not met.
Although mathematically consistent, the Lagrangian process which drives all the solute molecules using a single trajectory of the local dispersion and all the realizations of the velocity field is not a
physically satisfying model [Zirbel, 2001, p. 824]. Such a model does not describe the real dispersion for a representative of the random environment, for which the natural description is rather the ensemble of trajectories of an Itô process with drift coefficient given by a realization of the Eulerian velocity field. Therefore, we used this abstract Lagrangian process only to derive the memory-free dispersion (35) and to state the ergodic hypothesis which led to (36).
To write the dispersion terms (29-31) as moments of the joint probabilities (25) of the Itô process, we define the projection of the continuous Eulerian velocity V l ( x ) V l ( x ) V_(l)(x)V_{l}(\mathbf{x})Vl(x) on the trajectory X ( t ) X ( t ) X(t)\mathbf{X}(t)X(t) by using the Dirac measure δ δ delta\deltaδ,
V l ( X ( t ) ) = ∫ V l ( x ) δ [ x − X ( t ) ] d x V l ( X ( t ) ) = ∫ V l ( x ) δ [ x − X ( t ) ] d x V_(l)(X(t))=intV_(l)(x)delta[x-X(t)]dxV_{l}(\mathbf{X}(t))=\int V_{l}(\mathbf{x}) \delta[\mathbf{x}-\mathbf{X}(t)] d \mathbf{x}Vl(X(t))=∫Vl(x)δ[x−X(t)]dx
With c ( x 0 ) = c ( x , 0 ) c x 0 = c ( x , 0 ) c(x_(0))=c(x,0)c\left(\mathbf{x}_{0}\right)=c(\mathbf{x}, 0)c(x0)=c(x,0), the initial normalized concentration distribution given by (25) for n = 1 n = 1 n=1n=1n=1 and t = 0 t = 0 t=0t=0t=0, the two-times joint probability densities written as transition probabilities multiplied by c ( x 0 ) , p ( x 0 ; x , t ) = p ( x , t ∣ x 0 ) c ( x 0 ) c x 0 , p x 0 ; x , t = p x , t ∣ x 0 c x 0 c(x_(0)),p(x_(0);x,t)=p(x,t∣x_(0))c(x_(0))c\left(\mathbf{x}_{0}\right), p\left(\mathbf{x}_{0} ; \mathbf{x}, t\right)=p\left(\mathbf{x}, t \mid \mathbf{x}_{0}\right) c\left(\mathbf{x}_{0}\right)c(x0),p(x0;x,t)=p(x,t∣x0)c(x0), and the actual concentration c ( x , t ) = ∫ 0 t p ( x 0 ; x ′ , t ′ ) d x 0 c ( x , t ) = ∫ 0 t   p x 0 ; x ′ , t ′ d x 0 c(x,t)=int_(0)^(t)p(x_(0);x^('),t^('))dx_(0)c(\mathbf{x}, t)=\int_{0}^{t} p\left(\mathbf{x}_{0} ; \mathbf{x}^{\prime}, t^{\prime}\right) d \mathbf{x}_{0}c(x,t)=∫0tp(x0;x′,t′)dx0 [Gardiner, 1985], (29-31) become, respectively,
x ~ l l ( t ) = 2 D t + 2 ∫ 0 t d t ′ ∫ c ( x 0 ) d x 0 ∫ u l ( x ′ , t ′ ) [ x l ′ − x 0 l (37) − ∫ 0 t d t ′ ′ ∫ ⟨ V l ( x ′ ′ ) c ( x ′ ′ , t ′ ′ ) ⟩ V d x ′ ′ ] p ( x ′ , t ′ ∣ x 0 ) d x ′ m l l ( t ) = 2 ∫ 0 t d t ′ ∫ c ( x 0 ) d x 0 [ x 0 l − ∫ x 0 l ′ c ( x 0 ′ ) d x 0 ′ ] (38) ∫ u l ( x ′ , t ′ ) p ( x ′ , t ′ ∣ x 0 ) d x ′ r l l ( t ) = ∫ 0 t ∫ 0 t d t ′ d t ′ ′ ∬ c ( x 01 ) c ( x 02 ) d x 01 d x 02 (39) ∬ u l ( x ′ , t ′ ) u l ( x ′ ′ , t ′ ′ ) p ( x ′ , t ′ ∣ x 01 ) p ( x ′ ′ , t ′ ′ ∣ x 02 ) d x ′ d x ′ ′ x ~ l l ( t ) = 2 D t + 2 ∫ 0 t   d t ′ ∫ c x 0 d x 0 ∫ u l x ′ , t ′ x l ′ − x 0 l (37) − ∫ 0 t   d t ′ ′ ∫ V l x ′ ′ c x ′ ′ , t ′ ′ V d x ′ ′ p x ′ , t ′ ∣ x 0 d x ′ m l l ( t ) = 2 ∫ 0 t   d t ′ ∫ c x 0 d x 0 x 0 l − ∫ x 0 l ′ c x 0 ′ d x 0 ′ (38) ∫ u l x ′ , t ′ p x ′ , t ′ ∣ x 0 d x ′ r l l ( t ) = ∫ 0 t   ∫ 0 t   d t ′ d t ′ ′ ∬ c x 01 c x 02 d x 01 d x 02 (39) ∬ u l x ′ , t ′ u l x ′ ′ , t ′ ′ p x ′ , t ′ ∣ x 01 p x ′ ′ , t ′ ′ ∣ x 02 d x ′ d x ′ ′ {:[ widetilde(x)_(ll)(t)=2Dt+2int_(0)^(t)dt^(')int c(x_(0))dx_(0)intu_(l)(x^('),t^('))[x_(l)^(')-x_(0l):}],[(37){:-int_(0)^(t)dt^('')int(:V_(l)(x^(''))c(x^(''),t^('')):)_(V)dx^('')]p(x^('),t^(')∣x_(0))dx^(')],[m_(ll)(t)=2int_(0)^(t)dt^(')int c(x_(0))dx_(0)[x_(0l)-intx_(0l)^(')c(x_(0)^('))dx_(0)^(')]],[(38) intu_(l)(x^('),t^('))p(x^('),t^(')∣x_(0))dx^(')],[r_(ll)(t)=int_(0)^(t)int_(0)^(t)dt^(')dt^('')∬c(x_(01))c(x_(02))dx_(01)dx_(02)],[(39)∬u_(l)(x^('),t^('))u_(l)(x^(''),t^(''))p(x^('),t^(')∣x_(01))p(x^(''),t^('')∣x_(02))dx^(')dx^('')]:}\begin{align*} & \widetilde{x}_{l l}(t)=2 D t+2 \int_{0}^{t} d t^{\prime} \int c\left(\mathbf{x}_{0}\right) d \mathbf{x}_{0} \int u_{l}\left(\mathbf{x}^{\prime}, t^{\prime}\right)\left[x_{l}^{\prime}-x_{0 l}\right. \\ & \left.-\int_{0}^{t} d t^{\prime \prime} \int\left\langle V_{l}\left(\mathbf{x}^{\prime \prime}\right) c\left(\mathbf{x}^{\prime \prime}, t^{\prime \prime}\right)\right\rangle_{V} d \mathbf{x}^{\prime \prime}\right] p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}_{0}\right) d \mathbf{x}^{\prime} \tag{37}\\ & m_{l l}(t)=2 \int_{0}^{t} d t^{\prime} \int c\left(\mathbf{x}_{0}\right) d \mathbf{x}_{0}\left[x_{0 l}-\int x_{0 l}^{\prime} c\left(\mathbf{x}_{0}^{\prime}\right) d \mathbf{x}_{0}^{\prime}\right] \\ & \int u_{l}\left(\mathbf{x}^{\prime}, t^{\prime}\right) p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}_{0}\right) d \mathbf{x}^{\prime} \tag{38}\\ & r_{l l}(t)=\int_{0}^{t} \int_{0}^{t} d t^{\prime} d t^{\prime \prime} \iint c\left(\mathbf{x}_{01}\right) c\left(\mathbf{x}_{02}\right) d \mathbf{x}_{01} d \mathbf{x}_{02} \\ & \iint u_{l}\left(\mathbf{x}^{\prime}, t^{\prime}\right) u_{l}\left(\mathbf{x}^{\prime \prime}, t^{\prime \prime}\right) p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}_{01}\right) p\left(\mathbf{x}^{\prime \prime}, t^{\prime \prime} \mid \mathbf{x}_{02}\right) d \mathbf{x}^{\prime} d \mathbf{x}^{\prime \prime} \tag{39} \end{align*}x~ll(t)=2Dt+2∫0tdt′∫c(x0)dx0∫ul(x′,t′)[xl′−x0l(37)−∫0tdt′′∫⟨Vl(x′′)c(x′′,t′′)⟩Vdx′′]p(x′,t′∣x0)dx′mll(t)=2∫0tdt′∫c(x0)dx0[x0l−∫x0l′c(x0′)dx0′](38)∫ul(x′,t′)p(x′,t′∣x0)dx′rll(t)=∫0t∫0tdt′dt′′∬c(x01)c(x02)dx01dx02(39)∬ul(x′,t′)ul(x′′,t′′)p(x′,t′∣x01)p(x′′,t′′∣x02)dx′dx′′
These terms still depend on the realization of the velocity field. By taking a formal ensemble average of (37-39) and using (4), one obtains a more explicit form for the identity (17). For a point-like source, the memory term (38) vanishes and our result coincides with that obtained in an Eulerian frame by Kitanidis [1988, equations (29) and (33)], with the only difference that we considered a constant dispersion coefficient. The decomposition of the actual dispersion (4) in terms of form (37-39), for the general case of space-time dependent Eulerian velocity and local dispersion, also results from the definition of the drift and diffusion coefficients of the Fokker-Planck equation and the assumption that the first two moments of the concentration field are finite at finite times (Suciu et al., manuscript, 2005).
The terms of identity (17) can be also computed by ensemble averaging the first two moments of the concentration field, ⟨ X l α ⟩ D X 0 , α = 1 , 2 X l α D X 0 , α = 1 , 2 (:X_(l)^(alpha):)_(DX_(0)),alpha=1,2\left\langle X_{l}^{\alpha}\right\rangle_{D X_{0}}, \alpha=1,2⟨Xlα⟩DX0,α=1,2. Using (25) one obtains
⟨ ⟨ X l α ⟩ D X 0 ⟩ V = ∫ c ( x 0 ) d x 0 ∫ x l α ⟨ p ( x ′ , t ′ ∣ x 0 ) ⟩ V d x ′ X l α D X 0 V = ∫ c x 0 d x 0 ∫ x l α p x ′ , t ′ ∣ x 0 V d x ′ (:(:X_(l)^(alpha):)_(DX_(0)):)_(V)=int c(x_(0))dx_(0)intx_(l)^(alpha)(:p(x^('),t^(')∣x_(0)):)_(V)dx^(')\left\langle\left\langle X_{l}^{\alpha}\right\rangle_{D X_{0}}\right\rangle_{V}=\int c\left(\mathbf{x}_{0}\right) d \mathbf{x}_{0} \int x_{l}^{\alpha}\left\langle p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}_{0}\right)\right\rangle_{V} d \mathbf{x}^{\prime}⟨⟨Xlα⟩DX0⟩V=∫c(x0)dx0∫xlα⟨p(x′,t′∣x0)⟩Vdx′
If conditions weaker than those for regular fields but still sufficient for the homogeneity of the Lagrangian velocity field are met, then the ensemble average transition probability is translation invariant, as shown in Appendix B, and ⟨ ⟨ X l α ⟩ D X 0 ⟩ V X l α D X 0 V (:(:X_(l)^(alpha):)_(DX_(0)):)_(V)\left\langle\left\langle X_{l}^{\alpha}\right\rangle_{D X_{0}}\right\rangle_{V}⟨⟨Xlα⟩DX0⟩V become independent of initial conditions. Then M l l M l l M_(ll)M_{l l}Mll and Q l l Q l l Q_(ll)Q_{l l}Qll in (18) vanish, the difference between the second central moment of the ensemble average concentration and the initial second moment of the plume gives the one-particle dispersion, Σ l l − S l l ( 0 ) = X l l Σ l l − S l l ( 0 ) = X l l Sigma_(ll)-S_(ll)(0)=X_(ll)\Sigma_{l l}-S_{l l}(0)=X_{l l}Σll−Sll(0)=Xll, and (17) takes on the special form (14).

5 First-Order Results

With the exception of some solvable problems, as for instance in the case of transport in perfectly stratified flows [Balescu et al., 1994; Compte and Cáceres, 1998; Clincy and Kinzelbach, 2001], the ensemble average of (37-39) cannot be computed exactly. The difficulty arises from the highly non-linear dependence of probability densities in these relations on velocity V l ( X ( t ) ) V l ( X ( t ) ) V_(l)(X(t))V_{l}(\mathbf{X}(t))Vl(X(t)), which is sampled on trajectories which in their turn are determined by the velocity field [Balescu et al., 1994; Balescu, 2000]. Exact calculations of the ensemble averages require the knowledge of the entire infinite-dimensional probability distribution in a function space of all the realizations of the velocity field [Monin and Yaglom, 1971, p. 571]. Even in the advective case ( D = 0 D = 0 D=0D=0D=0 ) when this functional probability can be constructed, the functional integration can be carried out only in particular circumstances for which the results can be obtained directly [Lumley, 1962a]. In principle, the problem of ensemble averaging can be handled by iterating indefinitely the transport equation around an un-perturbed solution independent of velocity realization [Bouchaud and Georges, 1990; Balescu et al., 1994; Jaekel and Vereecken, 1997].
Iterations of the advection-dispersion equation (26) are usually done in Laplace-Fourier space, but one expects to obtain an exact ensemble average dispersion only when the velocity field is modeled as a white-noise in space [Compte and Cáceres, 1998]. Alternatively, the solution X ( t ) X ( t ) X(t)\mathbf{X}(t)X(t) of the Itô equation (21) can be approximated by the k k kkk-th iteration, starting with the unperturbed trajectory X ( 0 ) ( t ) X ( 0 ) ( t ) X^((0))(t)\mathbf{X}^{(0)}(t)X(0)(t),
(40) X l ( k ) ( t ) = X 0 l + ∫ 0 t V l ( X ( k − 1 ) ( t ′ ) ) d t ′ + W l ( t ) (40) X l ( k ) ( t ) = X 0 l + ∫ 0 t   V l X ( k − 1 ) t ′ d t ′ + W l ( t ) {:(40)X_(l)^((k))(t)=X_(0l)+int_(0)^(t)V_(l)(X^((k-1))(t^(')))dt^(')+W_(l)(t):}\begin{equation*} X_{l}^{(k)}(t)=X_{0 l}+\int_{0}^{t} V_{l}\left(\mathbf{X}^{(k-1)}\left(t^{\prime}\right)\right) d t^{\prime}+W_{l}(t) \tag{40} \end{equation*}(40)Xl(k)(t)=X0l+∫0tVl(X(k−1)(t′))dt′+Wl(t)
The second iteration is already highly complex. In this approximation, the argument of u l u l u_(l)u_{l}ul in (2931) has to be replaced by X ( 1 ) X ( 1 ) X^((1))\mathbf{X}^{(1)}X(1). As shown by (40), X ( 1 ) X ( 1 ) X^((1))\mathbf{X}^{(1)}X(1) depends on the values of the velocity V l V l V_(l)V_{l}Vl at X ( 0 ) ( t ′ ) X ( 0 ) t ′ X^((0))(t^('))\mathbf{X}^{(0)}\left(t^{\prime}\right)X(0)(t′) for all times between 0 and t t ttt. Hence, the infinite-dimensional probability distribution is needed to compute ensemble averages. Or, as we have seen, our knowledge about LagrangianEulerian relationships does not go beyond the one point statistics.
The approximation by the first iteration of the Itô equation, X l ( t ) ≈ X l ( 1 ) ( t ) X l ( t ) ≈ X l ( 1 ) ( t ) X_(l)(t)~~X_(l)^((1))(t)X_{l}(t) \approx X_{l}^{(1)}(t)Xl(t)≈Xl(1)(t), which is equivalent to the most-used first-order approximations in velocity fluctuations [Suciu et al., 2006b], simplifies matters considerably. Now, as follows from (40), the argument of V l V l V_(l)V_{l}Vl is replaced by X ( 0 ) X ( 0 ) X^((0))\mathbf{X}^{(0)}X(0), independent of the realization of the velocity field. The probability densities (25) computed from the unperturbed trajectories X ( 0 ) X ( 0 ) X^((0))\mathbf{X}^{(0)}X(0) become independent of velocity statistics as well. Since the nonlinearity of the initial problem is removed, useful ensemble averages can be computed under the only assumptions that the velocity field is statistically homogeneous and satisfies conditions for existence and uniqueness of the trajectories of the Itô process. For instance, the mean velocity on trajectory equals the Eulerian mean and the fluctuations vanish, accordingly to (32). It follows that the mean M l l M l l M_(ll)M_{l l}Mll of the memory term (30) vanishes. For the same reason, the ensemble mean of the last term in (34) is also zero. Further, by expressing (34) in terms of probability densities, similarly to (37), and taking the ensemble average one obtains the first-order approximation
(41) X ~ l l ( t ) = 2 D t + ∫ 0 t ∫ 0 t d t ′ d t ′ ′ ∫ c ( x 0 ) d x 0 ∬ p ( x ′ , t ′ ∣ x ′ ′ , t ′ ′ ) p ( x ′ ′ , t ′ ′ ∣ x 0 ) u l l ( x ′ , x ′ ′ ) d x ′ d x ′ ′ (41) X ~ l l ( t ) = 2 D t + ∫ 0 t   ∫ 0 t   d t ′ d t ′ ′ ∫ c x 0 d x 0 ∬ p x ′ , t ′ ∣ x ′ ′ , t ′ ′ p x ′ ′ , t ′ ′ ∣ x 0 u l l x ′ , x ′ ′ d x ′ d x ′ ′ {:(41) widetilde(X)_(ll)(t)=2Dt+int_(0)^(t)int_(0)^(t)dt^(')dt^('')int c(x_(0))dx_(0)∬p(x^('),t^(')∣x^(''),t^(''))p(x^(''),t^('')∣x_(0))u_(ll)(x^('),x^(''))dx^(')dx^(''):}\begin{equation*} \widetilde{X}_{l l}(t)=2 D t+\int_{0}^{t} \int_{0}^{t} d t^{\prime} d t^{\prime \prime} \int c\left(\mathbf{x}_{0}\right) d \mathbf{x}_{0} \iint p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}^{\prime \prime}, t^{\prime \prime}\right) p\left(\mathbf{x}^{\prime \prime}, t^{\prime \prime} \mid \mathbf{x}_{0}\right) u_{l l}\left(\mathbf{x}^{\prime}, \mathbf{x}^{\prime \prime}\right) d \mathbf{x}^{\prime} d \mathbf{x}^{\prime \prime} \tag{41} \end{equation*}(41)X~ll(t)=2Dt+∫0t∫0tdt′dt′′∫c(x0)dx0∬p(x′,t′∣x′′,t′′)p(x′′,t′′∣x0)ull(x′,x′′)dx′dx′′
where u l l ( x ′ , x ′ ′ ) = ⟨ u l ( x ′ ) u l ( x ′ ′ ) ⟩ V u l l x ′ , x ′ ′ = u l x ′ u l x ′ ′ V u_(ll)(x^('),x^(''))=(:u_(l)(x^('))u_(l)(x^('')):)_(V)u_{l l}\left(\mathbf{x}^{\prime}, \mathbf{x}^{\prime \prime}\right)=\left\langle u_{l}\left(\mathbf{x}^{\prime}\right) u_{l}\left(\mathbf{x}^{\prime \prime}\right)\right\rangle_{V}ull(x′,x′′)=⟨ul(x′)ul(x′′)⟩V is the Eulerian velocity correlation. The approximated variance of the center of mass is obtained directly by ensemble averaging (39),
R l l ( t ) = ∫ 0 t ∫ 0 t d t ′ d t ′ ′ ∬ c ( x 01 ) c ( x 02 ) d x 01 d x 02 ∬ p ( x ′ , t ′ ∣ x 01 ) p ( x ′ ′ , t ′ ′ ∣ x 02 ) u l l ( x ′ , x ′ ′ ) d x ′ d x ′ ′ R l l ( t ) = ∫ 0 t   ∫ 0 t   d t ′ d t ′ ′ ∬ c x 01 c x 02 d x 01 d x 02 ∬ p x ′ , t ′ ∣ x 01 p x ′ ′ , t ′ ′ ∣ x 02 u l l x ′ , x ′ ′ d x ′ d x ′ ′ R_(ll)(t)=int_(0)^(t)int_(0)^(t)dt^(')dt^('')∬c(x_(01))c(x_(02))dx_(01)dx_(02)∬p(x^('),t^(')∣x_(01))p(x^(''),t^('')∣x_(02))u_(ll)(x^('),x^(''))dx^(')dx^('')R_{l l}(t)=\int_{0}^{t} \int_{0}^{t} d t^{\prime} d t^{\prime \prime} \iint c\left(\mathbf{x}_{01}\right) c\left(\mathbf{x}_{02}\right) d \mathbf{x}_{01} d \mathbf{x}_{02} \iint p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}_{01}\right) p\left(\mathbf{x}^{\prime \prime}, t^{\prime \prime} \mid \mathbf{x}_{02}\right) u_{l l}\left(\mathbf{x}^{\prime}, \mathbf{x}^{\prime \prime}\right) d \mathbf{x}^{\prime} d \mathbf{x}^{\prime \prime}Rll(t)=∫0t∫0tdt′dt′′∬c(x01)c(x02)dx01dx02∬p(x′,t′∣x01)p(x′′,t′′∣x02)ull(x′,x′′)dx′dx′′.
If the unperturbed solution X ( 0 ) X ( 0 ) X^((0))\mathbf{X}^{(0)}X(0) is a diffusion in the ensemble mean velocity field U U U\mathbf{U}U [see e.g. Dentz et al., 2000; Fiori and Dagan, 2000] the probability densities in (41-42) are those of a Gaussian process of mean x 0 + U t x 0 + U t x_(0)+Ut\mathbf{x}_{0}+\mathbf{U} tx0+Ut and dispersion 2 D t 2 D t 2Dt2 D t2Dt. Since the Gaussian transition probability
and the Eulerian correlation of the homogeneous velocity field are invariant to space translations, (41) is independent of x 0 x 0 x_(0)\mathbf{x}_{0}x0. In the previous Section we have seen that Q l l = 0 Q l l = 0 Q_(ll)=0Q_{l l}=0Qll=0 as a consequence of independence on initial conditions of the velocity fluctuations (32). Hence, it follows from (8) that (41) is a first-order approximation of the one-particle dispersion X l l X l l X_(ll)X_{l l}Xll. The mean M l l M l l M_(ll)M_{l l}Mll also vanishes, and from (11) we have an explicit form of relation (14) which approximates the mean dispersion S l l S l l S_(ll)S_{l l}Sll with the aid of (41-42). From (19), using (41), one can see that the ensemble dispersion coefficient 1 2 d Σ l l / d t 1 2 d Σ l l / d t (1)/(2)dSigma_(ll)//dt\frac{1}{2} d \Sigma_{l l} / d t12dΣll/dt is also memory-free. This result was obtained by Dentz et al. [2000], with a perturbation approach using a Fourier representation for u l l u l l u_(ll)u_{l l}ull and for the concentration field of the unperturbed problem. Due to the simplicity of this approximation which uses Gaussian probability densities, the averages (41-42) and the corresponding standard deviations can be computed, at least by numerical integration, for given two- and four-points Eulerian correlations. We note that such computations can be achieved without using Fourier representations (which remain, however, essential tools for higher order approximations [Jaekel and Vereecken, 1997]).
For the advection-dominated transport problem considered in simulations presented in Section 2 the velocity fluctuations are of the order of P e − 1 / 2 P e − 1 / 2 Pe^(-1//2)P e^{-1 / 2}Pe−1/2 and consistent asymptotic expansions of (21) are obtained for the unperturbed solution X ( 0 ) X ( 0 ) X^((0))\mathbf{X}^{(0)}X(0) given by the deterministic trajectory X 0 + U t X 0 + U t X_(0)+Ut\mathbf{X}_{0}+\mathbf{U} tX0+Ut of the ensemble mean velocity [Suciu et al., 2006b]. Then, as follows from (25), the conditional probabilities from (41-42) degenerate to Dirac delta measures, p ( x ′ , t ′ ∣ x 0 ) = δ [ x ′ − ( x 0 + U t ′ ) ] p x ′ , t ′ ∣ x 0 = δ x ′ − x 0 + U t ′ p(x^('),t^(')∣x_(0))=delta[x^(')-(x_(0)+Ut^('))]p\left(\mathbf{x}^{\prime}, t^{\prime} \mid \mathbf{x}_{0}\right)=\delta\left[\mathbf{x}^{\prime}-\left(\mathbf{x}_{0}+\mathbf{U} t^{\prime}\right)\right]p(x′,t′∣x0)=δ[x′−(x0+Ut′)] and p ( x ′ , t ′ ; x ′ ′ , t ′ ′ ∣ x 0 ) = δ [ x ′ − ( x 0 + U t ′ ) ] δ [ x ′ ′ − ( x 0 + U t ′ ′ ) ] p x ′ , t ′ ; x ′ ′ , t ′ ′ ∣ x 0 = δ x ′ − x 0 + U t ′ δ x ′ ′ − x 0 + U t ′ ′ p(x^('),t^(');x^(''),t^('')∣x_(0))=delta[x^(')-(x_(0)+Ut^('))]delta[x^('')-(x_(0)+Ut^(''))]p\left(\mathbf{x}^{\prime}, t^{\prime} ; \mathbf{x}^{\prime \prime}, t^{\prime \prime} \mid \mathbf{x}_{0}\right)=\delta\left[\mathbf{x}^{\prime}-\left(\mathbf{x}_{0}+\mathbf{U} t^{\prime}\right)\right] \delta\left[\mathbf{x}^{\prime \prime}-\left(\mathbf{x}_{0}+\mathbf{U} t^{\prime \prime}\right)\right]p(x′,t′;x′′,t′′∣x0)=δ[x′−(x0+Ut′)]δ[x′′−(x0+Ut′′)] and one obtains
(43) X ~ l l ( t ) = 2 D t + ∫ 0 t ∫ 0 t d t ′ d t ′ ′ ∫ c ( x 0 ) u l l ( x 0 + U t ′ ; x 0 + U t ′ ′ ) d x 0 (44) R l l ( t ) = ∫ 0 t ∫ 0 t d t ′ d t ′ ′ ∬ c ( x 01 ) c ( x 02 ) u l l ( x 01 + U t ′ ; x 02 + U t ′ ′ ) d x 01 d x 02 (43) X ~ l l ( t ) = 2 D t + ∫ 0 t   ∫ 0 t   d t ′ d t ′ ′ ∫ c x 0 u l l x 0 + U t ′ ; x 0 + U t ′ ′ d x 0 (44) R l l ( t ) = ∫ 0 t   ∫ 0 t   d t ′ d t ′ ′ ∬ c x 01 c x 02 u l l x 01 + U t ′ ; x 02 + U t ′ ′ d x 01 d x 02 {:[(43) widetilde(X)_(ll)(t)=2Dt+int_(0)^(t)int_(0)^(t)dt^(')dt^('')int c(x_(0))u_(ll)(x_(0)+Ut^(');x_(0)+Ut^(''))dx_(0)],[(44)R_(ll)(t)=int_(0)^(t)int_(0)^(t)dt^(')dt^('')∬c(x_(01))c(x_(02))u_(ll)(x_(01)+Ut^(');x_(02)+Ut^(''))dx_(01)dx_(02)]:}\begin{gather*} \widetilde{X}_{l l}(t)=2 D t+\int_{0}^{t} \int_{0}^{t} d t^{\prime} d t^{\prime \prime} \int c\left(\mathbf{x}_{0}\right) u_{l l}\left(\mathbf{x}_{0}+\mathbf{U} t^{\prime} ; \mathbf{x}_{0}+\mathbf{U} t^{\prime \prime}\right) d \mathbf{x}_{0} \tag{43}\\ R_{l l}(t)=\int_{0}^{t} \int_{0}^{t} d t^{\prime} d t^{\prime \prime} \iint c\left(\mathbf{x}_{01}\right) c\left(\mathbf{x}_{02}\right) u_{l l}\left(\mathbf{x}_{01}+\mathbf{U} t^{\prime} ; \mathbf{x}_{02}+\mathbf{U} t^{\prime \prime}\right) d \mathbf{x}_{01} d \mathbf{x}_{02} \tag{44} \end{gather*}(43)X~ll(t)=2Dt+∫0t∫0tdt′dt′′∫c(x0)ull(x0+Ut′;x0+Ut′′)dx0(44)Rll(t)=∫0t∫0tdt′dt′′∬c(x01)c(x02)ull(x01+Ut′;x02+Ut′′)dx01dx02
These are well known expressions [e.g. Dagan, 1991] which relate the ensemble average dispersion (14) with the one- and two-particles velocity correlations. Since the velocity field is assumed to be statistically homogeneous, the one-particle correlation depends only on t ′ ′ − t ′ t ′ ′ − t ′ t^('')-t^(')t^{\prime \prime}-t^{\prime}t′′−t′. Obviously, for regular velocity fields (43) is a first-order approximation of the memory-free dispersion (35). Unlike in (42), the two-particle correlation from (44) depends only on t ′ ′ − t ′ t ′ ′ − t ′ t^('')-t^(')t^{\prime \prime}-t^{\prime}t′′−t′ and the separation lag x 02 − x 01 x 02 − x 01 x_(02)-x_(01)\mathbf{x}_{02}-\mathbf{x}_{01}x02−x01 [see also Fiori and Dagan, 2000]. This proves that the "Lagrangian stationarity" assumed by Dagan [1990] holds true in this consistent first-order approximation, provided that the Itô equation has unique solutions and the velocity field is statistically homogeneous. We note that incompressibility and continuity of the velocity derivatives are not necessary to derive these first-order results.
For the transport problem considered in Section 2, the approximations (41) and (43) of the memory-free one-particle dispersion are very close to each other at all times [Suciu et al., 2006b, Figure 2]. The first-order approximations (41-42) also give reasonably good estimations for the ensemble average S l l S l l S_(ll)S_{l l}Sll but they yield considerable overestimation, by 40 % 40 % 40%40 \%40% to 80 % 80 % 80%80 \%80%, of the fluctuations of the actual dispersion at early times [Suciu et al., 2006b, Figure 4]. Therefore, such approximations might not be accurate enough for investigations of memory effects which manifest themselves in the pre-asymptotic regime.

6 Memory Effects and Ergodicity

The alternative to analytical computations is to derive the mean values and the standard deviations of all terms in (4) by direct simulations followed by ensemble averages. For the term x ~ l l x ~ l l widetilde(x)_(ll)\widetilde{x}_{l l}x~ll, this approach is equivalent to individually tracking large numbers of particles and to compute the velocity-displacements correlation under the integral in (29). Since the statistical stability of second moments s l l s l l s_(ll)s_{l l}sll was ensured by 10 10 10 10 10^(10)10^{10}1010 particles [Suciu et al., 2006a], a particle-tracking simulation with the same precision is practically impossible [see e.g. Vamoş et al., 2003]. For memory terms m l l m l l m_(ll)m_{l l}mll, however, direct simulations for a single realization were done as follows. We considered a
number of 121 initial positions X 0 X 0 X_(0)\mathbf{X}_{0}X0 uniformly distributed in domains L 1 λ × L 2 λ L 1 λ × L 2 λ L_(1)lambda xxL_(2)lambdaL_{1} \lambda \times L_{2} \lambdaL1λ×L2λ. By releasing 10 10 10 10 10^(10)10^{10}1010 particles from each initial position, we performed GRW simulations to compute the velocity fluctuations along the trajectories of the processes starting from these positions. Finally, the correlations velocity-initial position from (30) were computed by averages over the initial positions. The results presented in Figures 6 and 7 are consistent with memory effects shown, indirectly, in Figures 3 and 5 by the behavior of the mean and standard deviations of the moments of the ensemble averaged concentrations: m l l m l l m_(ll)m_{l l}mll increase with the dimension of the source in l l lll-direction and are mainly significant for asymmetric sources. The increase with the transverse dimension of the source of the transverse memory terms in single realizations of transport is also indicated by the results of Fiori and Jancović [2005, Figure 5], who computed a quantity which corresponds to our (30) when local dispersion is neglected.
Repeating the above numerical experiment for a large ensemble of velocity realizations is very costly in terms of computing resources. Therefore, to evaluate means and standard deviations of the memory terms, we consider only the following cases of large sources: slabs ( λ , 100 λ ) , ( 100 λ , λ ) ( λ , 100 λ ) , ( 100 λ , λ ) (lambda,100 lambda),(100 lambda,lambda)(\lambda, 100 \lambda),(100 \lambda, \lambda)(λ,100λ),(100λ,λ), and square ( 100 λ , 100 λ ) ( 100 λ , 100 λ ) (100 lambda,100 lambda)(100 \lambda, 100 \lambda)(100λ,100λ). As shown below, for these cases the memory-free moments X l l X l l X_(ll)X_{l l}Xll and the memory terms can be estimated from the simulation results presented in Section 2.
Figures 8 and 9 show the mean moments Σ l l ( t ) − S l l ( 0 ) Σ l l ( t ) − S l l ( 0 ) Sigma_(ll)(t)-S_(ll)(0)\Sigma_{l l}(t)-S_{l l}(0)Σll(t)−Sll(0) for point source and for slabs of 100 λ 100 λ 100 lambda100 \lambda100λ perpendicular to l l lll-direction, for which the memory terms are negligible (see Figures 6 and 7). The good agreement of mean values for point and slab sources indicates that, according to (19), Σ l l ( t ) − S l l ( 0 ) Σ l l ( t ) − S l l ( 0 ) Sigma_(ll)(t)-S_(ll)(0)\Sigma_{l l}(t)-S_{l l}(0)Σll(t)−Sll(0) approximates the memory-free dispersion X l l X l l X_(ll)X_{l l}Xll. In the following we shall use the estimations
X 11 ( t ) ≈ Σ 11 ( t ; 1 , 100 ) − S 11 ( 0 ; 1 , 100 ) (45) X 22 ( t ) ≈ Σ 22 ( t ; 100 , 1 ) − S 22 ( 0 ; 100 , 1 ) X 11 ( t ) ≈ Σ 11 ( t ; 1 , 100 ) − S 11 ( 0 ; 1 , 100 ) (45) X 22 ( t ) ≈ Σ 22 ( t ; 100 , 1 ) − S 22 ( 0 ; 100 , 1 ) {:[X_(11)(t)~~Sigma_(11)(t;1","100)-S_(11)(0;1","100)],[(45)X_(22)(t)~~Sigma_(22)(t;100","1)-S_(22)(0;100","1)]:}\begin{align*} & X_{11}(t) \approx \Sigma_{11}(t ; 1,100)-S_{11}(0 ; 1,100) \\ & X_{22}(t) \approx \Sigma_{22}(t ; 100,1)-S_{22}(0 ; 100,1) \tag{45} \end{align*}X11(t)≈Σ11(t;1,100)−S11(0;1,100)(45)X22(t)≈Σ22(t;100,1)−S22(0;100,1)
We also represented in Figures 8 and 9 the first-order consistent approximations of memory-free moments computed according to (43) from the Eulerian correlation of the numerically generated velocity. Since their deviations from simulation results is much smaller than the memory effects indicated by Figure 5, the approximated memory-free moments can also be used for practical purposes as estimates of X l l X l l X_(ll)X_{l l}Xll [Suciu et al., 2007a]. The standard deviations S D [ σ l l ] S D σ l l SD[sigma_(ll)]S D\left[\sigma_{l l}\right]SD[σll] (represented by thin lines in Figures 8 and 9) for slab sources are much smaller than for point source (and even smaller when compared with the other cases from Figure 5), so that we can approximate σ l l − S l l ( 0 ) ≈ X l l , l = 1 , 2 σ l l − S l l ( 0 ) ≈ X l l , l = 1 , 2 sigma_(ll)-S_(ll)(0)~~X_(ll),l=1,2\sigma_{l l}-S_{l l}(0) \approx X_{l l}, l=1,2σll−Sll(0)≈Xll,l=1,2. But for negligible memory terms σ l l − S l l ( 0 ) ≈ x ~ l l σ l l − S l l ( 0 ) ≈ x ~ l l sigma_(ll)-S_(ll)(0)~~ widetilde(x)_(ll)\sigma_{l l}-S_{l l}(0) \approx \widetilde{x}_{l l}σll−Sll(0)≈x~ll, as follows from (15), thus x ~ l l ≈ X l l x ~ l l ≈ X l l widetilde(x)_(ll)~~X_(ll)\widetilde{x}_{l l} \approx X_{l l}x~ll≈Xll. Since the variance of center of mass (Figure 1) is also negligible, we can adopt the ergodic hypothesis (36) which allows us to estimate the memory terms by the approximate relation (20) as
(46) m l l ( t ) = s l l ( t ) − S l l ( 0 ) − X l l ( t ) (46) m l l ( t ) = s l l ( t ) − S l l ( 0 ) − X l l ( t ) {:(46)m_(ll)(t)=s_(ll)(t)-S_(ll)(0)-X_(ll)(t):}\begin{equation*} m_{l l}(t)=s_{l l}(t)-S_{l l}(0)-X_{l l}(t) \tag{46} \end{equation*}(46)mll(t)=sll(t)−Sll(0)−Xll(t)
with X l l X l l X_(ll)X_{l l}Xll given by (45).
The results obtained from GRW simulations by using (45-46) for longitudinal and transverse memory terms are presented in Figures 10 and 11, respectively. The thin lines in these figures correspond to M l l / ( 2 D t ) ± S D [ m l l / ( 2 D t ) ] / R 1 / 2 M l l / ( 2 D t ) ± S D m l l / ( 2 D t ) / R 1 / 2 M_(ll)//(2Dt)+-SD[m_(ll)//(2Dt)]//R^(1//2)M_{l l} /(2 D t) \pm S D\left[m_{l l} /(2 D t)\right] / R^{1 / 2}Mll/(2Dt)±SD[mll/(2Dt)]/R1/2, where S D [ m l l ] = S D [ s l l ] S D m l l = S D s l l SD[m_(ll)]=SD[s_(ll)]S D\left[m_{l l}\right]=S D\left[s_{l l}\right]SD[mll]=SD[sll], according to (46), and R = 1024 R = 1024 R=1024R=1024R=1024 is the number of realizations. The mean memory terms M 11 M 11 M_(11)M_{11}M11 for longitudinal slabs (Figure 10) and M 22 M 22 M_(22)M_{22}M22 for transverse slabs (Figure 11) have significant non-vanishing values at finite times. The comparison with the case P e = ∞ P e = ∞ Pe=ooP e=\inftyPe=∞ shows that there are no important differences with respect to the purely advective transport. For isotropic sources, ( 100 λ , 100 λ ) ( 100 λ , 100 λ ) (100 lambda,100 lambda)(100 \lambda, 100 \lambda)(100λ,100λ), the averages M l l M l l M_(ll)M_{l l}Mll are smaller than the local dispersion 2 D t 2 D t 2Dt2 D t2Dt. In the pre-asymptotic regime the standard deviations (see Figure 5) are considerably large: S D [ m l l ( 100 , 100 ) ] / ( 2 D t ) ∼ 8 S D m l l ( 100 , 100 ) / ( 2 D t ) ∼ 8 SD[m_(ll)(100,100)]//(2Dt)∼8S D\left[m_{l l}(100,100)\right] /(2 D t) \sim 8SD[mll(100,100)]/(2Dt)∼8 for l = 1 , 2 , S D [ m 22 ( 1 , 100 ) ] / ( 2 D t ) ∼ 20 l = 1 , 2 , S D m 22 ( 1 , 100 ) / ( 2 D t ) ∼ 20 l=1,2,SD[m_(22)(1,100)]//(2Dt)∼20l=1,2, S D\left[m_{22}(1,100)\right] /(2 D t) \sim 20l=1,2,SD[m22(1,100)]/(2Dt)∼20, and S D [ m 11 ( 100 , 1 ) ] / ( 2 D t ) ∼ 100 S D m 11 ( 100 , 1 ) / ( 2 D t ) ∼ 100 SD[m_(11)(100,1)]//(2Dt)∼100S D\left[m_{11}(100,1)\right] /(2 D t) \sim 100SD[m11(100,1)]/(2Dt)∼100. The memory terms show a remarkable persistence, over hundreds of dimensionless times. However, the decrease of both mean values and standard deviations indicates that the single realization memory terms m l l m l l m_(ll)m_{l l}mll converge to zero for t → ∞ t → ∞ t rarr oot \rightarrow \inftyt→∞, in the mean square limit.
The non-vanishing mean memory terms disagree with the theoretical result for regular fields presented in Section 4. To check whether this discrepancy is a numerical artifact, we shortly comment here possible sources of errors. The number of 10 10 10 10 10^(10)10^{10}1010 particles guarantees highly accurate single
Figure 6: Longitudinal memory terms for different sources of dimensions ( L 1 , L 2 L 1 , L 2 L_(1),L_(2)L_{1}, L_{2}L1,L2 ); direct GRW simulations for a given velocity field.
Figure 8: Longitudinal memory-free moments estimated from GRW simulations for large transverse slabs.
Figure 7: Transverse memory terms for different sources of dimensions ( L 1 , L 2 L 1 , L 2 L_(1),L_(2)L_{1}, L_{2}L1,L2 ); direct GRW simulations for a given velocity field.
Figure 9: Transverse memory-free moments estimated from GRW simulations for large longitudinal slabs.
realization transport simulations [Suciu et al., 2006a]. Averaging over 1024 velocity realizations also ensures reliable estimations, at least for the transverse mean memory terms, as shown by the confidence intervals represented by thin lines in Figure 11. As for the longitudinal mean memory terms in Figure 10, where the value 0 lies inside, or is at the limit of the confidence intervals, we have checked that increasing the number of realizations to 2048 (results not presented here) one obtains a trustworthy estimation, like that shown in Figure 11. Since the non-vanishing mean memory terms are not a consequence of the finite size of the statistical ensembles, they may be attributed to the bad behavior of the sample velocity fields. It is well known that, excepting some particular cases, exponentially correlated fields are not differentiable [see e.g. Yaglom, 1987, p. 67 and references therein]. Hence, there is no guarantee that results for regular fields or those derived from assumed translation invariance of the mean transition probability apply. By increasing the number of harmonics in the Kraichnan routine, the numerically generated velocity field approximates a Gaussian field with given correlation. It was shown that 6400 periodic modes, as in the present simulations, accurately reproduces the behavior of transport in Gaussian fields over thousands of heterogeneity scales [Eberhard et al., 2007]. We also note that for a coarser approximation, by a superposition of only 640 random sine modes with exponential correlation, the transverse memory terms were found to be smaller than half the local dispersion uniformly in time, hence they are negligible as compared to the resolution of the simulations [Suciu et al., 2007c]. One expects that for Gaussian correlation, which yields smooth differentiable fields [Trefry et al., 2003], the mean memory terms can be neglected as well. Of coarse, it is also possible that
Figure 10: Longitudinal memory terms for large longitudinal slab and square sources.
Figure 11: Transverse memory terms for large transverse slab and square sources.
mean memory terms arise as a consequence of the discrete nature of numerical simulations, which always reproduce the nominal values of the velocity statistics with some finite precision [see e.g. Gleeson, 2002, Figure 3; Trefry et al., 2003, Table 1; Suciu et al., 2006a, Table B1]. Relatively weak dependence of Lagrangian velocity on initial positions, caused by either non-smooth sample velocities or finite precision of simulations, [see e.g. Suciu and Vamoş, 2007, Figure 2; Suciu et al., 2008, Figures 1 and 2], for which the variance of the ensemble mean velocity and, consequently, the term Q l l ( 10 ) Q l l ( 10 ) Q_(ll)(10)Q_{l l}(10)Qll(10) are negligible, can yet produce significant mean memory terms when multiplied by large [ X 0 l − ⟨ X 0 l ⟩ X 0 ] X 0 l − X 0 l X 0 [X_(0l)-(:X_(0l):)_(X_(0))]\left[X_{0 l}-\left\langle X_{0 l}\right\rangle_{X_{0}}\right][X0l−⟨X0l⟩X0], as shown by (30).
For the cases presented in Figures 10 and 11, the memory terms (46) measure the deviation of the dispersions s l l − S l l ( 0 ) s l l − S l l ( 0 ) s_(ll)-S_(ll)(0)s_{l l}-S_{l l}(0)sll−Sll(0) from the memory-free moments X l l X l l X_(ll)X_{l l}Xll. Since the mean memory terms are small as compared with their fluctuations, the ergodicity ranges are practically determined by the standard deviations S D [ m l l ] = S D [ s l l ] S D m l l = S D s l l SD[m_(ll)]=SD[s_(ll)]S D\left[m_{l l}\right]=S D\left[s_{l l}\right]SD[mll]=SD[sll] presented in Figure 5. But this nice feature breaks down for moderately large plumes, for which the center of mass fluctuations r l l r l l r_(ll)r_{l l}rll cannot be neglected (see Figure 1). Then, the ergodicity range η [ s l l − S l l ( 0 ) ] = [ η 1 2 + η 2 2 ] 1 / 2 η s l l − S l l ( 0 ) = η 1 2 + η 2 2 1 / 2 eta[s_(ll)-S_(ll)(0)]=[eta_(1)^(2)+eta_(2)^(2)]^(1//2)\eta\left[s_{l l}-S_{l l}(0)\right]=\left[\eta_{1}^{2}+\eta_{2}^{2}\right]^{1 / 2}η[sll−Sll(0)]=[η12+η22]1/2 is given by the deviation of the mean S l l − S l l ( 0 ) − X l l = η 1 S l l − S l l ( 0 ) − X l l = η 1 S_(ll)-S_(ll)(0)-X_(ll)=eta_(1)S_{l l}-S_{l l}(0)-X_{l l}=\eta_{1}Sll−Sll(0)−Xll=η1 and the standard deviation S D [ s l l ] = η 2 S D s l l = η 2 SD[s_(ll)]=eta_(2)S D\left[s_{l l}\right]=\eta_{2}SD[sll]=η2. The results for all the initial conditions of the GRW simulations presented in Section 2 are given in Figure 12. One remarks that for times smaller than 100 U t λ 100 U t λ 100 Ut lambda100 U t \lambda100Utλ the ergodicity range is of the order of the local dispersion 2 D t 2 D t 2Dt2 D t2Dt or smaller, in the case of large slab sources perpendicular to l l lll-axis, but it is very large for slab sources parallel with l l lll and for symmetric sources. In the latter case, η η eta\etaη is practically the same as the ergodicity range with respect to the macrodispersion moments presented in Figure 4. The ergodicity ranges with respect to the consistent first-order approximation (43) of X l l [ X l l [ X_(ll)[X_{l l}[Xll[ Suciu et al., 2007a] were found to be in good agreement with the ergodicity ranges with respect to the numerically inferred memory-free moments presented in Figure 12.

7 Discussion of Results and Conclusions

The memory terms describe the persistent influence of the shape and dimension of the source on solute dispersion and are defined, according to (30) or (38), by correlations between initial positions and velocity fluctuations along trajectories integrated in time. Using the Itô formalism we have shown that the memory term m l l m l l m_(ll)m_{l l}mll naturally occurs as a component of the second moment of the concentration field, whenever the velocity varies in space and the source has a non-zero extension on the l l lll-direction. For smooth, incompressible, and statistically homogeneous velocity fields (i.e. regular fields), well known results on the equality of Lagrangian and Eulerian one-point probability distributions imply that the average over the ensemble of velocity realizations of m l l m l l m_(ll)m_{l l}mll is zero. Under the same conditions, the one-particle dispersion X l l X l l X_(ll)X_{l l}Xll is shown to be independent of initial conditions. First-order approximations lead to the same results under weakened assumptions, without requiring
Figure 12: Ergodicity ranges η η eta\etaη with respect to memory-free moments of the second moment of the actual plume for different shapes and extensions of the initial plume.
continuity of the velocity derivatives and incompressibility.
Numerical GRW simulations for transport in Gaussian fields have shown that, contrary to theoretical prediction for regular fields, non-vanishing mean memory terms of the order of local dispersion can occur for large sources. This discrepancy can be tracked back to the non-differentiability of the Gaussian fields with exponential correlations or to the discrete nature of the numerical simulations. Nevertheless, this numerical finding might have serious implications for real situations, for which the model of regular field is always an idealization [see also Suciu and Vamoş, 2007]. But the memory terms are mainly important for single realizations, which correspond to observable transport process. Our simulations indicate that the larger the initial plumes are, the stronger is the dependence on the initial state as described by memory terms. For a typical situation of contaminant transport in aquifers, considered in simulations, and for sources of tens of heterogeneity scales, the memory terms can be tens to hundreds times larger than the local dispersion. (Remember that the value of the upscaled longitudinal macrodispersion, in the same conditions, is of about ten local dispersions).
The memory terms for P e = 100 P e = 100 Pe=100P e=100Pe=100 and P e = ∞ P e = ∞ Pe=ooP e=\inftyPe=∞ (Figures 10 and 11) are quite close to each other and both decay to zero in the large time limit. With regard to the advective transport ( P e = ∞ P e = ∞ Pe=ooP e=\inftyPe=∞ ) in aquifers, however, this result is not conclusive. The numerically generated Gaussian field with finite correlation range is a first-order approximation of the solution of Darcy and continuity equations. For deterministic Darcy flows it was shown by Sposito [2001] that trajectories are confined on invariant subsets of the flow domain. It is therefore possible that solute molecules, when driven by Darcy flows and in absence of local dispersion, never lose the memory of initial position and memory terms persist indefinitely, as suggested by numerical investigations of Fiori and Jancović [2005].
The other extreme case is when the relative displacement X ~ l X ~ l widetilde(X)_(l)\widetilde{X}_{l}X~l is independent of initial position X 0 l X 0 l X_(0l)X_{0 l}X0l and the memory term m l l m l l m_(ll)m_{l l}mll defined by (3) vanishes. This happens, for instance, when X ~ l X ~ l widetilde(X)_(l)\widetilde{X}_{l}X~l is the superposition of a diffusion process and a uniform movement with constant velocity, as in the case of perfectly stratified aquifer model of Matheron and de Marsily [1980]. The initial conditions can yet influence the dispersion in case of confined stratified flows, as in a single fracture in geological media,
when the transverse dimension of the source governs the interplay between the local dispersion and the coherent cross-section velocity profile [see e.g. Dentz and Carrera, 2007]. In this case, the significant relation for the longitudinal dispersion is obtained from (2) by adding and subtracting ⟨ ⟨ X ~ 1 ⟩ D 2 ⟩ X 0 X ~ 1 D 2 X 0 (:(: widetilde(X)_(1):)_(D)^(2):)_(X_(0))\left\langle\left\langle\widetilde{X}_{1}\right\rangle_{D}^{2}\right\rangle_{X_{0}}⟨⟨X~1⟩D2⟩X0, where ⟨ X ~ 1 ⟩ D X ~ 1 D (: widetilde(X)_(1):)_(D)\left\langle\widetilde{X}_{1}\right\rangle_{D}⟨X~1⟩D is the center of mass of a solute plume originating from a point source,
s 11 = S 11 ( 0 ) + ⟨ ⟨ [ X ~ 1 − ⟨ X ~ 1 ⟩ D ] 2 ⟩ D ⟩ X 0 + ⟨ [ ⟨ X ~ 1 ⟩ D − ⟨ X ~ 1 ⟩ D X 0 ] 2 ⟩ X 0 s 11 = S 11 ( 0 ) + X ~ 1 − X ~ 1 D 2 D X 0 + X ~ 1 D − X ~ 1 D X 0 2 X 0 s_(11)=S_(11)(0)+(:(:[ widetilde(X)_(1)-(: widetilde(X)_(1):)_(D)]^(2):)_(D):)_(X_(0))+(:[(: widetilde(X)_(1):)_(D)-(: widetilde(X)_(1):)_(DX_(0))]^(2):)_(X_(0))s_{11}=S_{11}(0)+\left\langle\left\langle\left[\widetilde{X}_{1}-\left\langle\widetilde{X}_{1}\right\rangle_{D}\right]^{2}\right\rangle_{D}\right\rangle_{X_{0}}+\left\langle\left[\left\langle\widetilde{X}_{1}\right\rangle_{D}-\left\langle\widetilde{X}_{1}\right\rangle_{D X_{0}}\right]^{2}\right\rangle_{X_{0}}s11=S11(0)+⟨⟨[X~1−⟨X~1⟩D]2⟩D⟩X0+⟨[⟨X~1⟩D−⟨X~1⟩DX0]2⟩X0
The last term in the above relation is a spatial variance of the point-source center of mass which carries the memory of initial conditions: when it becomes negligible, the reduced dispersion s 11 − S 11 ( 0 ) s 11 − S 11 ( 0 ) s_(11)-S_(11)(0)s_{11}- S_{11}(0)s11−S11(0) behaves as a superposition of point-source dispersions, ⟨ ⟨ [ X ~ 1 − ⟨ X ~ 1 ⟩ D ] 2 ⟩ D ⟩ X 0 X ~ 1 − X ~ 1 D 2 D X 0 (:(:[ widetilde(X)_(1)-(: widetilde(X)_(1):)_(D)]^(2):)_(D):)_(X_(0))\left\langle\left\langle\left[\widetilde{X}_{1}-\left\langle\widetilde{X}_{1}\right\rangle_{D}\right]^{2}\right\rangle_{D}\right\rangle_{X_{0}}⟨⟨[X~1−⟨X~1⟩D]2⟩D⟩X0.
In this study we limited ourselves to the analysis of memory effects induced by deterministic uniform initial concentrations. One expects that for non-uniform initial solute concentration proportional to velocity [Demmy et al., 1999], or in the case of a random initial state which is the outcome of a postinjection stage of the transport [Sposito and Dagan, 1994], memory effects can show new and relevant features. A postinjection initial position, X τ l = X l ( τ ) X τ l = X l ( τ ) X_(tau l)=X_(l)(tau)X_{\tau l}=X_{l}(\tau)Xτl=Xl(τ), is the solution of the Itô equation (21) at t = τ t = τ t=taut=\taut=τ and the evolution of the solute molecule is described by (21) rewritten as
(47) X l ( t − τ ) = X τ l + ∫ τ t V l ( X ( t ′ ) ) d t ′ + W l ( t − τ ) (47) X l ( t − τ ) = X τ l + ∫ τ t   V l X t ′ d t ′ + W l ( t − τ ) {:(47)X_(l)(t-tau)=X_(tau l)+int_(tau)^(t)V_(l)(X(t^(')))dt^(')+W_(l)(t-tau):}\begin{equation*} X_{l}(t-\tau)=X_{\tau l}+\int_{\tau}^{t} V_{l}\left(\mathbf{X}\left(t^{\prime}\right)\right) d t^{\prime}+W_{l}(t-\tau) \tag{47} \end{equation*}(47)Xl(t−τ)=Xτl+∫τtVl(X(t′))dt′+Wl(t−τ)
The dispersion (4) can be expressed as
s l l ( t − τ ) = s l l ( τ ) + x ~ l l ( t − τ ) + m l l ( t − τ ) − r l l ( t − τ ) (48) = S l l ( 0 ) + x ~ l l ( τ ) + m l l ( τ ) − r l l ( τ ) + x ~ l l ( t − τ ) + m l l ( t − τ ) − r l l ( t − τ ) s l l ( t − τ ) = s l l ( τ ) + x ~ l l ( t − τ ) + m l l ( t − τ ) − r l l ( t − τ ) (48) = S l l ( 0 ) + x ~ l l ( τ ) + m l l ( τ ) − r l l ( τ ) + x ~ l l ( t − τ ) + m l l ( t − τ ) − r l l ( t − τ ) {:[s_(ll)(t-tau)=s_(ll)(tau)+ widetilde(x)_(ll)(t-tau)+m_(ll)(t-tau)-r_(ll)(t-tau)],[(48)=S_(ll)(0)+ widetilde(x)_(ll)(tau)+m_(ll)(tau)-r_(ll)(tau)+ widetilde(x)_(ll)(t-tau)+m_(ll)(t-tau)-r_(ll)(t-tau)]:}\begin{align*} s_{l l}(t-\tau) & =s_{l l}(\tau)+\widetilde{x}_{l l}(t-\tau)+m_{l l}(t-\tau)-r_{l l}(t-\tau) \\ & =S_{l l}(0)+\widetilde{x}_{l l}(\tau)+m_{l l}(\tau)-r_{l l}(\tau)+\widetilde{x}_{l l}(t-\tau)+m_{l l}(t-\tau)-r_{l l}(t-\tau) \tag{48} \end{align*}sll(t−τ)=sll(τ)+x~ll(t−τ)+mll(t−τ)−rll(t−τ)(48)=Sll(0)+x~ll(τ)+mll(τ)−rll(τ)+x~ll(t−τ)+mll(t−τ)−rll(t−τ)
Now, let us assume that at time t t ttt the plume is large enough to render the fluctuations of the center of mass negligible, i.e r l l ≈ 0 r l l ≈ 0 r_(ll)~~0r_{l l} \approx 0rll≈0. Disregarding the memory terms and equating (48) and (4) one obtains
(49) x ~ l l ( t ) = x ~ l l ( τ ) + x ~ l l ( t − τ ) (49) x ~ l l ( t ) = x ~ l l ( τ ) + x ~ l l ( t − τ ) {:(49) widetilde(x)_(ll)(t)= widetilde(x)_(ll)(tau)+ widetilde(x)_(ll)(t-tau):}\begin{equation*} \widetilde{x}_{l l}(t)=\widetilde{x}_{l l}(\tau)+\widetilde{x}_{l l}(t-\tau) \tag{49} \end{equation*}(49)x~ll(t)=x~ll(τ)+x~ll(t−τ)
Equation (49) shows that x ~ l l x ~ l l widetilde(x)_(ll)\widetilde{x}_{l l}x~ll is a linear function of time. Evidently, this is true only if x ~ l l x ~ l l widetilde(x)_(ll)\widetilde{x}_{l l}x~ll, which for large plumes approximates the one-particle dispersion X l l X l l X_(ll)X_{l l}Xll, has reached the macrodispersion limit behavior [Sposito and Dagan, 1994, page 587]. However, one cannot conclude that memory terms never vanish in the pre-asymptotic regime. In fact, for point or line sources perpendicular to l l lll, m l l ( τ ) = 0 m l l ( τ ) = 0 m_(ll)(tau)=0m_{l l}(\tau)=0mll(τ)=0 at all τ τ tau\tauτ. To get more insight, we analyze the memory term m l l ( t − τ ) m l l ( t − τ ) m_(ll)(t-tau)m_{l l}(t-\tau)mll(t−τ), computed from (30) by using (21) and (47),
(50) m l l ( t − τ ) = m 1 , l l ( t − τ ) + m 2 , l l ( t − τ ) m 1 , l l ( t − τ ) = 2 ∫ τ t ⟨ [ X 0 l − ⟨ X 0 l ⟩ X 0 ] u l ( X ( t ′ ) , t ′ ) ⟩ D X 0 d t ′ m 2 , l l ( t − τ ) = 2 ∫ τ t ∫ 0 τ [ ⟨ u l ( X ( t ′ ) , t ′ ) u l ( X ( t ′ ′ ) , t ′ ′ ) ⟩ D X 0 − ⟨ u l ( X ( t ′ ) , t ′ ) ⟩ D X 0 ⟨ u l ( X ( t ′ ′ ) , t ′ ′ ) ⟩ D X 0 ] d t ′ d t ′ ′ (50) m l l ( t − τ ) = m 1 , l l ( t − τ ) + m 2 , l l ( t − τ ) m 1 , l l ( t − τ ) = 2 ∫ τ t   X 0 l − X 0 l X 0 u l X t ′ , t ′ D X 0 d t ′ m 2 , l l ( t − τ ) = 2 ∫ τ t   ∫ 0 τ   u l X t ′ , t ′ u l X t ′ ′ , t ′ ′ D X 0 − u l X t ′ , t ′ D X 0 u l X t ′ ′ , t ′ ′ D X 0 d t ′ d t ′ ′ {:[(50)m_(ll)(t-tau)=m_(1,ll)(t-tau)+m_(2,ll)(t-tau)],[m_(1,ll)(t-tau)=2int_(tau)^(t)(:[X_(0l)-(:X_(0l):)_(X_(0))]u_(l)(X(t^(')),t^(')):)_(DX_(0))dt^(')],[m_(2,ll)(t-tau)=2int_(tau)^(t)int_(0)^(tau)[(:u_(l)(X(t^(')),t^('))u_(l)(X(t^('')),t^('')):)_(DX_(0)):}],[{:-(:u_(l)(X(t^(')),t^(')):)_(DX_(0))(:u_(l)(X(t^('')),t^('')):)_(DX_(0))]dt^(')dt^('')]:}\begin{align*} & m_{l l}(t-\tau)=m_{1, l l}(t-\tau)+m_{2, l l}(t-\tau) \tag{50}\\ & m_{1, l l}(t-\tau)=2 \int_{\tau}^{t}\left\langle\left[X_{0 l}-\left\langle X_{0 l}\right\rangle_{X_{0}}\right] u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right)\right\rangle_{D X_{0}} d t^{\prime} \\ & m_{2, l l}(t-\tau)=2 \int_{\tau}^{t} \int_{0}^{\tau}\left[\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right) u_{l}\left(\mathbf{X}\left(t^{\prime \prime}\right), t^{\prime \prime}\right)\right\rangle_{D X_{0}}\right. \\ & \left.-\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime}\right), t^{\prime}\right)\right\rangle_{D X_{0}}\left\langle u_{l}\left(\mathbf{X}\left(t^{\prime \prime}\right), t^{\prime \prime}\right)\right\rangle_{D X_{0}}\right] d t^{\prime} d t^{\prime \prime} \end{align*}(50)mll(t−τ)=m1,ll(t−τ)+m2,ll(t−τ)m1,ll(t−τ)=2∫τt⟨[X0l−⟨X0l⟩X0]ul(X(t′),t′)⟩DX0dt′m2,ll(t−τ)=2∫τt∫0τ[⟨ul(X(t′),t′)ul(X(t′′),t′′)⟩DX0−⟨ul(X(t′),t′)⟩DX0⟨ul(X(t′′),t′′)⟩DX0]dt′dt′′
The term m 1 , l l ( t − τ ) m 1 , l l ( t − τ ) m_(1,ll)(t-tau)m_{1, l l}(t-\tau)m1,ll(t−τ) describes the memory effect induced by deterministic initial conditions and vanishes as well if X 0 l = ⟨ X 0 l ⟩ X 0 X 0 l = X 0 l X 0 X_(0l)=(:X_(0l):)_(X_(0))X_{0 l}=\left\langle X_{0 l}\right\rangle_{X_{0}}X0l=⟨X0l⟩X0. The second term of (50), m 2 , l l ( t − τ ) m 2 , l l ( t − τ ) m_(2,ll)(t-tau)m_{2, l l}(t-\tau)m2,ll(t−τ), describes the memory effect produced by correlations between velocities sampled on trajectories of the process starting at τ τ tau\tauτ and velocities on segments of trajectories in the interval 0 ≤ t ≤ τ 0 ≤ t ≤ τ 0 <= t <= tau0 \leq t \leq \tau0≤t≤τ. Since m 1 , l l ( t − τ ) + m l l ( τ ) = m l l ( t ) m 1 , l l ( t − τ ) + m l l ( τ ) = m l l ( t ) m_(1,ll)(t-tau)+m_(ll)(tau)=m_(ll)(t)m_{1, l l}(t-\tau)+m_{l l}(\tau)=m_{l l}(t)m1,ll(t−τ)+mll(τ)=mll(t) whatever the initial conditions are, it follows that, if the fluctuations of the center of mass r l l r l l r_(ll)r_{l l}rll can be neglected, a sufficient condition for linear dispersion (49) is a vanishing m 2 , l l ( t − τ ) m 2 , l l ( t − τ ) m_(2,ll)(t-tau)m_{2, l l}(t-\tau)m2,ll(t−τ). Hence, terms of type m 2 , l l m 2 , l l m_(2,ll)m_{2, l l}m2,ll describe memory effects as persistence of the pre-asymptotic regime. We also note that relations (47-48), written for successive times τ τ tau\tauτ and for displacements of the solute molecules
described by the GRW algorithm, result in a numerical procedure useful in analyses of dispersion in terms of correlation structure of the velocity field [Vamoş et al., 2001].
Statistical mechanics approaches for spatially homogeneous systems [Zwanzig, 1961; Mokshin et al., 2005] or for transport in inhomogeneous and rapidly fluctuating velocity fields, as in plasma physics [Balescu, 2000], are concerned with memory effects which characterize the departure of the process from diffusive behavior, like those described by m 2 , l l m 2 , l l m_(2,ll)m_{2, l l}m2,ll above. Such effects are related to the occurrence of anomalous transport [Balescu, 2000; Morgado et al., 2002] and are usually described by non-local equations, containing "memory kernels", governing the velocity autocorrelation function or the ensemble average concentration. Similar effects in groundwater systems are captured by memory kernels in models for small scale [Berkowitz and Scher, 2001] as well as in models for integrated larger scales [Cortis et al., 2004; Morales-Casique et al., 2006].
Whereas the memory terms m 2 , l l m 2 , l l m_(2,ll)m_{2, l l}m2,ll are relevant for ergodic behavior with respect to the macrodispersion process, the terms of type m 1 , l l m 1 , l l m_(1,ll)m_{1, l l}m1,ll investigated here have a bearing on ergodicity with respect to the memory-free one-particle dispersion X l l X l l X_(ll)X_{l l}Xll, which characterizes the pre-asymptotic regime. Such ergodic properties are neither the usual ergodicity nor the self-averaging property, because the reference is not just the ensemble average of the observable (actual dispersion) but a stochastic model of it (macrodispersion or one-particle dispersion) [Suciu et al., 2006a]. However, the ergodicity of the velocity field with finite correlation range, in the common sense that unbiased space average estimators converge to ensemble averages, is intimately related to the pre-asymptotic behavior of dispersion. Increasing the source dimensions in all space directions reduces the uncertainty caused by the randomness of the center of mass and the ensemble average of the actual dispersion approaches the second moment of the ensemble average concentration. In the past this has led to the thinking that large plumes are ergodic plumes. But, theoretical considerations on ergodic Gaussian regular fields and numerical simulations also show that for large sources the uncertainty is mainly quantified by memory terms. We have found that, for a given direction, the memory term increases with the increase of the source dimension in the same direction. The dependence of the memory effects on the source extension and anisotropy should be incorporated into the prediction process and in procedures for identification of the contaminant source from available measurement data. For small sources, a best fit of memory-free dispersion X l l X l l X_(ll)X_{l l}Xll and measured second moments can underestimate the variance and the correlation length of the hydraulic conductivity [Suciu et al., 2006b]. Erroneous estimations will also be obtained for sources with large extension on the l l lll direction. Instead, for large narrow sources perpendicular to l l lll, the ergodic behavior of the actual dispersion with respect to X l l X l l X_(ll)X_{l l}Xll can be used in practice to improve the parameter identification from field experiments.

Appendix A: Monte Carlo Convergence

To assess the Monte Carlo convergence with the number of velocity realizations n n nnn of a statistical estimate f f fff we used an estimation error defined as
E ( f ( n ) ) = | f ( n + Δ n ) − f ( n ) | E ( f ( n ) ) = | f ( n + Δ n ) − f ( n ) | E(f(n))=|f(n+Delta n)-f(n)|E(f(n))=|f(n+\Delta n)-f(n)|E(f(n))=|f(n+Δn)−f(n)|
Figures 13 and 14 show the decrease of errors with n n nnn for the second moments of the ensemble average concentration Σ l l Σ l l Sigma_(ll)\Sigma_{l l}Σll and for the corresponding standard deviations S D ( σ l l ) , l = 1 , 2 S D σ l l , l = 1 , 2 SD(sigma_(ll)),l=1,2S D\left(\sigma_{l l}\right), l=1,2SD(σll),l=1,2. The results correspond to initial conditions consisting of slab sources λ × 100 λ λ × 100 λ lambda xx100 lambda\lambda \times 100 \lambdaλ×100λ oriented across the l l lll direction, which required the larger number of realization to stabilize the Monte Carlo oscillations of the ensemble average estimations. For longitudinal quantities ( l = 1 ) ( l = 1 ) (l=1)(l=1)(l=1) the errors were estimated at dimensionless times of 1 and 20 (where the ensemble averages were most sensitive to n n nnn ) while for transverse ones ( l = 2 l = 2 l=2l=2l=2 ) the estimation times were 1 and 10 . The increment of the number of realizations was set to Δ n = 32 Δ n = 32 Delta n=32\Delta n=32Δn=32 in both cases. We found that longitudinal quantities E E EEE require more realizations to stabilize (Figure 13) but for n ≥ 1000 n ≥ 1000 n >= 1000n \geq 1000n≥1000 realizations E E EEE falls under the value D t D t DtD tDt,
Figure 13: Monte Carlo convergence of mean and standard deviation of longitudinal second moments.
Figure 14: Monte Carlo convergence of mean and standard deviation of transverse second moments.
which corresponds to half the local dispersion.

Appendix B: Translation Invariant Mean Transition Probability

The density of the transition probability p ( x , t ∣ x 0 , 0 ) p x , t ∣ x 0 , 0 p(x,t∣x_(0),0)p\left(\mathbf{x}, t \mid \mathbf{x}_{0}, 0\right)p(x,t∣x0,0) for the Itô process (21) is the solution of the Fokker-Planck equation (26) for initial condition p ( x , 0 ∣ x 0 , 0 ) = δ ( x − x 0 ) p x , 0 ∣ x 0 , 0 = δ x − x 0 p(x,0∣x_(0),0)=delta(x-x_(0))p\left(\mathbf{x}, 0 \mid \mathbf{x}_{0}, 0\right)=\delta\left(\mathbf{x}-\mathbf{x}_{0}\right)p(x,0∣x0,0)=δ(x−x0).
In the following we show that for translation-invariant statistics of the Eulerian velocity, and under conditions which ensure the statistical homogeneity if the Lagrangian velocity field, the ensemble average transition probability is invariant to space translations as well, i.e.
(B1) ⟨ p ( x , t ∣ x 0 , 0 ) ⟩ V = p ¯ ( x − x 0 , t ) . (B1) p x , t ∣ x 0 , 0 V = p ¯ x − x 0 , t . {:(B1)(:p(x,t∣x_(0),0):)_(V)= bar(p)(x-x_(0),t).:}\begin{equation*} \left\langle p\left(x, t \mid x_{0}, 0\right)\right\rangle_{V}=\bar{p}\left(x-x_{0}, t\right) . \tag{B1} \end{equation*}(B1)⟨p(x,t∣x0,0)⟩V=p¯(x−x0,t).
(To avoid cumbersome notations, we consider the one-dimensional problem only. The extension to more dimensions is straightforward.)
To proceed, let us consider the Fourier transform of p p ppp (i.e. the characteristic function) expressed as a series of moments μ k = ⟨ X ( t ) k ⟩ D μ k = X ( t ) k D mu_(k)=(:X(t)^(k):)_(D)\mu_{k}=\left\langle X(t)^{k}\right\rangle_{D}μk=⟨X(t)k⟩D,
(B2) ϕ ( s ) = ∫ − ∞ ∞ e i s x p ( x ) d x = ∑ k = 0 ∞ ( i s ) k k ! μ k . (B2) ϕ ( s ) = ∫ − ∞ ∞   e i s x p ( x ) d x = ∑ k = 0 ∞   ( i s ) k k ! μ k . {:(B2)phi(s)=int_(-oo)^(oo)e^(isx)p(x)dx=sum_(k=0)^(oo)((is)^(k))/(k!)mu_(k).:}\begin{equation*} \phi(s)=\int_{-\infty}^{\infty} e^{i s x} p(x) d x=\sum_{k=0}^{\infty} \frac{(i s)^{k}}{k!} \mu_{k} . \tag{B2} \end{equation*}(B2)ϕ(s)=∫−∞∞eisxp(x)dx=∑k=0∞(is)kk!μk.
The transition probability p ( x ) p ( x ) p(x)p(x)p(x) is uniquely specified by its characteristic function ϕ ( s ) ϕ ( s ) phi(s)\phi(s)ϕ(s). Because of the linearity of the expansion (B2), the ensemble average p ¯ p ¯ bar(p)\bar{p}p¯ is also determined by ϕ ¯ ϕ ¯ bar(phi)\bar{\phi}ϕ¯ and μ k ― μ k ¯ bar(mu_(k))\overline{\mu_{k}}μk―. The ensemble average moments μ k ― μ k ¯ bar(mu_(k))\overline{\mu_{k}}μk― computed from Itô equation (21) depend on velocity statistics only through the k k kkk-order moments of the Lagrangian velocity, evaluated on trajectories starting from the same initial position x 0 x 0 x_(0)x_{0}x0,
(B3) ∫ 0 t ⋯ ∫ 0 t ⟨ V ( X ( t 1 ) ) ⋯ V ( X ( t k ) ) ⟩ V d t 1 ⋯ d t k (B3) ∫ 0 t   ⋯ ∫ 0 t   V X t 1 ⋯ V X t k V d t 1 ⋯ d t k {:(B3)int_(0)^(t)cdotsint_(0)^(t)(:V(X(t_(1)))cdots V(X(t_(k))):)_(V)dt_(1)cdotsdt_(k):}\begin{equation*} \int_{0}^{t} \cdots \int_{0}^{t}\left\langle V\left(X\left(t_{1}\right)\right) \cdots V\left(X\left(t_{k}\right)\right)\right\rangle_{V} \mathrm{~d} t_{1} \cdots \mathrm{~d} t_{k} \tag{B3} \end{equation*}(B3)∫0t⋯∫0t⟨V(X(t1))⋯V(X(tk))⟩V dt1⋯ dtk
If the Eulerian velocity field is statistically homogeneous and the necessary conditions for the existence and the uniqueness of the solutions of the Itô equation and for their measurability with respect to the probability space of the Eulerian velocity are fulfilled, then the Lagrangian velocity is statistically homogeneous as well [Zirbel, 2001, Remark 6.1 and Proposition 7.1]. Consequently the moments of the Lagrangian velocity from (B3) are independent of x 0 x 0 x_(0)x_{0}x0. From this the conclusion (B1) follows.
We note that the translation invariance of the ensemble average transition probability is ensured by conditions milder than those for equality of one-point probability distributions of Eulerian and
Lagrangian velocity, which also require the continuity of the partial derivatives and the incompressibility of the Eulerian random velocity field.
Acknowledgments The research reported in this paper was supported by Deutsche Forschungsgemeinschaft grant SU 415/1-2, Project JICG41 at Jülich Supercomputing Centre, Fundamental Program No. 1 of the Romanian Academy, Romanian Ministry of Education and Research grant 2-CEx06-11-96, NATO Collaborative Linkage Grant ESP.NR.CLG 981426, RFBR Grant 06-01-00498. The authors thank C. Andronache, G. Grün, A. Millet, G. Keller, P. Kramer, O. Kurbanmuradov, F. Radu, U. Rüde, M. Trefry, and W. Wagner for stimulating discussions.

References

[1] Balescu, R. (2000), Memory effects in plasma transport theory, Plasma Phys. Control. Fusion 42, B1-B13.
[2] Balescu, R., H.-D. Wang, and J. H. Misguich (1994), Langevin equation versus kinetic equation: Subdiffusive behavior of charged particles in a stochastic magnetic field, Phys. Plasmas 1(12), 3826-3842.
[3] Berkowitz, B., and H. Scher (2001), Probabilistic approaches to transport in heterogeneous media, Transport in Porous Media, 42, 241-263.
[4] Bhattacharya, R. N., and V. K. Gupta (1983), A Theoretical Explanation of Solute Dispersion in Saturated Porous Media at the Darcy Scale, Water Resour. Res., 19, 934-944.
[5] Bouchaud, J.-P., and A. Georges (1990), Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep., 195, 127-293.
[6] Brouwers, J. J. H. (2002), On diffusion theory in turbulence, J. Eng. Mat., 44, 277-295.
[7] Chilès, J. P., and P. Delfiner (1999), Geostatisctics: Modeling Spatial Uncertainty, John Wiley & Sons, New York.
[8] Clincy, M., and H. Kinzelbach (2001), Stratified disordered media: exact solutions for transport parameters and their self-averaging properties, J. Phys. A: Math. Gen. 34, 7142-7152.
[9] Compte, A., and M. O. Cáceres (1998), Fractional dynamics in random velocity fields, Phys. Rev. Lett. 81 (15), 3140-3143.
[10] Cortis, A., H. Scher, and B. Berkowitz (2004), Numerical simulation of non-Fickian transport in geological formations with multiple scale heterogeneity, Water Resour. Res, 40, W04209, doi:10.1029/2003WRR002750.
[11] Dagan, G. (1987), Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, 183-215.
[12] Dagan, G. (1990), Transport in heterogeneous porous formations: Spatial moments, ergodicity and effective dispersion, Water Resour. Res., 26, 1281-1290.
[13] Dagan, G. (1991), Dispersion of a solute in non-ergodic transport in steady velocity fields in heterogeneous formations, J. Fluid Mech., 233, 197-210.
[14] Dagan, G. (2004), Comment on 'Exact averaging of stochastic equations for transport in random velocity field', Transport in porous media 50, 223-241, 2003, and 'Probability density functions for solute transport in random field' Transport in porous media 50, 243-266, 2003, Transport in porous media, 55, 113-116.
[15] Demmy, G., S. Berglund, and W. Graham (1999), Injection mode implications for solute transport in porous media: Analysis in a stochastic Lagrangian framework, Water Resour. Res., 35(7), 1965-1973.
[16] 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.
[17] Dentz, M., and J. Carrera (2005), Effective solute transport in temporally fluctuating flow through heterogeneous media, Water Resour. Res., 41, W08414, doi:10.1029/2004WR003571.
[18] Dentz, M., and J. Carrera (2007), Mixing and spreading in stratified flow, Phys. Fluids 19, 017107, doi:10.1063/1.2427089.
[19] Doob, J. L. (1953), Stochastic Processes, John Wiley & Sons, London.
[20] Eberhard J., N. Suciu, and C. Vamoş (2007), On the self-averaging of dispersion for transport in quasi-periodic random media, J. Phys. A: Math. Theor., 40, 597-610, doi: 10.1088/17518113/40/4/002.
[21] Fannjiang, A., and T. Komorowski (1999), Turbulent diffusion in Markovian flows, Ann. Appl. Probab., 9(3), 591-610.
[22] Fiori, A., and G. Dagan (2000), Concentration fluctuations in aquifer transport: a rigorous first-order solution and applications, J.Contam. Hydrol. 45, 139-163.
[23] Fiori, A., and I. Jancović (2005), Can we determine the transverse macrodispersivity by using the method of moments?, Adv. Water Resour., 28, 589-599, doi:10.1016/j.advwaters.2004.09.909.
[24] Gardiner, C. W. (1985), Handbook of Stochastic Methods (for Physics, Chemistry and Natural Science), Springer, New York.
[25] Gleeson, J. P. (2002), Comment on "Diffusion in biased turbulence", Phys. Rev. E, 66, 038301.
[26] Jaekel, U., and H. Vereecken (1997), Renormalization group analysis of macrodispersion in a directed random flow, Water Resour. Res., 33, 2287-2299.
[27] Kitanidis, P. K. (1988), Prediction by the method of moments of transport in a heterogeneous formation, J. Hydrol., 102, 453-473.
[28] Kloeden, P. E. and E, Platen (1995), Numerical Solutions of Stochastic Differential Equations, Springer, Berlin.
[29] Komorowski, T., and G. Papanicolaou (1997), Motion in a Gaussian incompressible flow, Ann. Appl. Probab., 7(1), 229-264.
[30] Kramer, P. R., and A. J. Majda (2007), Lectures on Turbulent Diffusion, Springer (in preparation).
[31] LaBolle, E. M., J. Quastel, G. E. Fogg, and J. Gravner (2000), Diffusion processes in composite media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coefficients, Water Resour. Res., 36(3), 651-662.
[32] Le Doussal, P., and J. Machta (1989), Annealed versus quenched diffusion coefficient in random media, Phys. Rev. B, 40(12), 9427-9430.
[33] Lichtner, P. C., S. Kelkar, and B. Robinson (2002), New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking, Water Resour. Res., 38(8), 10.1029/2000WR000100.
[34] Lumley, J. L. (1962a), An approach to the Eulerian-Lagrangian problem, J. Math. Phys., 3(2), 309-312.
[35] Lumley, J. L. (1962b), The mathematical nature of the problem of relating Lagrangian and Eulerian statistical functions in turbulence, pp , 17-26 in Mécanique de la Turbulence (Coll. Intern. du CNRS à Marseille), Ed. CNRS, Paris.
[36] Lundgren, T. S. (1981), Turbulent pair dispersion and scalar diffusion, J. Fluid Mech. 111, 27-57.
[37] Matheron, G., and G. de Marsily (1980), Is transport in porous media always diffusive?, Water Resour. Res., 16, 901-917.
[38] Monin, A. S., and A. M. Yaglom (1975) Statistical Fluid Mechanics: Mechanics of Turbulence, MIT Press, Cambridge, M A.
[39] Morales-Casique, E., S.P. Neuman, and A. Gaudagnini (2006), Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: Theoretical framework, Adv. Water Resour. 29(8), 1238-1255, doi:10.1016/j.advwatres.2005.10.002.
[40] Morgado, R., F. A. Oliveira, G. G. Batrouni, and A. Hansen (2002), Relation between anomalous and normal difusion in systems with memory, Phys. Rev. Lett., 89(10), 100601.
[41] Mokshin, A. V., Yulmetyev, R. M., and P. Hnggi (2005), Simple measure of memory for dynamical processes described by a generalized Langevin equation, Phys. Rev. Lett., 95(20), 200601.
[42] Naff, R. L., D. F. Haley, and E. A. Sudicky (1998), High-resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media 2. Transport results, Water Resour. Res., 34(4), 679-697.
[43] Port, S. C., and C. J. Stone (1976), Random measures and their application to motion in an incompressible fluid. J. Appl. Prob., 13, 498506.
[44] Sposito, G. (2001), Topological groundwater hydrodynamics, Adv. Water Resour., 24,, 793801.
[45] Sposito, G. and G. Dagan (1994), Predicting solute plume evolution in heterogeneous porous formations, Water Resour. Res., 30(2), 585-589.
[46] Suciu N., C. Vamoş, J. Vanderborght, H. Hardelauf, and H. Vereecken (2006a), Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi:10.1029/2005WR004546.
[47] Suciu N., C. Vamoş, and J. Eberhard (2006b), Evaluation of the first-order approximations for transport in heterogeneous media, Water Resour. Res. 42, W11504, doi: 10.1029/2005WR004714.
[48] Suciu N., and C. Vamoş (2007), Comment on "Nonstationary flow and nonergodic transport in random porous media" by G. Darvini and P. Salandin, Water Resour. Res., 43, W12601, doi:10.1029/2007WR005946.
[49] Suciu N., C. Vamoş, and H. Vereecken (2007a), Multiple meanings of ergodicity in real life problems, pp.196-211 in Proceedings of the 5th Workshop on Mathematical Modelling of Environmental and Life Sciences Problems, Ed. G. Marinoschi, S. Ion and C. Popa, Rom. Acad. Publishing House, Bucharest.
[50] Suciu N., C. Vamoş, and K. Sabelfeld (2007b), Ergodic simulations of diffusion in random velocity fields, pp. 659-668 in Monte Carlo and Quasi-Monte Carlo Methods 2006, Ed. A. Keller, S. Heinrich, and H. Niederreiter, Springer Verlag, Heidelberg.
[51] Suciu N., C. Vamoş, K. Sabelfeld, and. C. Andronache (2007c), Memory effects and ergodicity for diffusion in spatially correlated velocity fields, Proc. Appl. Math. Mech., 7, 20100152010016, doi:10.1002/pamm. 20070057.
[52] Suciu N., F. A. Radu, and C. Vamoş (2008) Lagrangian-Eulerian statistics and numerical modeling of transport in homogeneous random media, in Proceedings of the 6th Workshop on Mathematical Modelling of Environmental and Life Sciences Problems, Ed. G. Marinoschi, S. Ion, and C. Popa, Rom. Acad. Publishing House, Bucharest (to appear).
[53] Trefry, M. G., F. P. Ruan, and D. McLaughlin (2003), Numerical simulations of preasymptotic transport in heterogeneous porous media: Departures from the Gaussian limit, Water Resour. Res., 39 (3), 1063, doi:10.1029/2001WRR001101.
[54] Vamoş, C., N. Suciu, and H. Vereecken (2003), Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186(2), 527-544, doi:10.1016/S0021-9991(03)00073-1.
[55] Vamoş, C., N. Suciu, H. Vereecken, J. Vanderborht, and O. Nitzsche (2001), Path decomposition of discrete effective diffusion coefficient, Internal Report ICG-IV.00501, Research Center Jülich.
[56] van Kampen, N. G. (1981), Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam.
[57] Wentzell, A. D. (1981), A Course in the Theory of Stochastic Processes, McGraw-Hill, New York.
[58] Yaglom, A. M. (1987) Correlation Theory of Stationary and Related Random Functions, Volume I: Basic Results, Springer-Verlag, New York.
[59] Zirbel, C., L. (2001), Lagrangian observations of homogeneous random environments, A d v A d v AdvA d vAdv. Appl. Prob., 33, 810-835.
[60] Zhang, Y.-K., and B.-m. Seo (2004), Stochastic analyses and Monte-Carlo simulations of nonergodic solute transport in three-dimensional heterogeneous statistically anisotropic aquifers, Water Resour. Res., 40, W05103, doi:10.1029/2003WR002871.
[61] Zwanzig, R. (1961), Memory effects and irreversible thermodynamics, Phys. Rev. 124(4), 983992.

Related Posts

No results found.