Itô equation model for dispersion of solutes in heterogeneous media

7 years ago

Abstract

Transport processes in heterogeneous media such as ionized plasmas, natural porous media, and turbulent atmosphere are often modeled as diffusion processes in random velocity fields.

Using the Itô formalism, we decompose the second spatial moments of the concentration and the equivalent effective dispersion coefficients in terms corresponding to various physical factors which influence the transport.

We explicitly define “ergodic” dispersion coefficients, independent of the initial conditions and completely determined by local dispersion coefficients and velocity correlations. Ergodic coefficients govern an upscaled process which describes the transport at large tine-space scales. The non-ergodic behavior at finite times shown by numerical experiments for large initial plumes is explained by “memory terms” accounting for correlations between initial positions and velocity fluctuations on the trajectories of the solute molecules.

Authors

N. Suciu
Friedrich-Alexander University of Erlangen-Nuremberg, Institute of Applied Mathematics

C. Vamos
Tiberiu Popoviciu Institute of Numerical Analysis (Romanian Academy)

H. Vereecken
Research Center Julich, ICG-IV: Agrosphere Institute

K. Sabelefld
Weierstrass Institute for Applied Analysis and Stochastics, Berlin
Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch of Russian Academy of Sciences

P. Knabner
Friedrich-Alexander University of Erlangen-Nuremberg, Institute of Applied Mathematics, Germany

Keywords

Ito equation; random fields; memory effects; ergodicity.

References

See the expanding block below.

Paper coordinates

N. Suciu, C. Vamoş, H. Vereecken, K. Sabelfeld P. Knabner, Itô equation model for dispersion of solutes in heterogeneous media, Rev. Anal. Numér. Théor. Approx., 37 (2008) no. 2, 221-238

PDF

https://ictp.acad.ro/jnaat/journal/article/view/2008-vol37-no2-art15/Suciu-Vamos-Vereecken-Sabelfeld-Knabner-pdf

About this paper

Journal

Rev. Anal. Numér. Théor. Approx.

Publisher Name

Editions de l’Academie Roumaine

Paper on journal website

https://ictp.acad.ro/jnaat/journal/article/view/2008-vol37-no2-art15

Print ISSN

1222-9024

Online ISSN

2457-8126

MR

ZBL

Google Scholar

[1] Balescu, R., Wang, H.-D. and Misguich, J. H., Langevin equation versus kinetic equation: Subdiffusive behavior of charged particles in a stochastic magnetic field, Phys. Plasmas 1(12), pp. 3826-3842, 1994.

[2] Bhattacharya, R. N. and Gupta, V. K., A Theoretical Explanation of Solute Dispersion in Saturated Porous Media at the Darcy Scale, Water Resour. Res., 19, pp. 934-944, 1983.

[3] Bouchaud, J.-P. and Georges, A., Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep., 195, pp. 127-293, 1990.

[4] Brouwers, J. J. H., On diffusion theory in turbulence, J. Eng. Mat., 44, pp. 277–295, 2002.

[5] Compte, A. and Cáceres, M. O., Fractional dynamics in random velocity fields, Phys. Rev. Lett., 81(15), pp. 3140-3143, 1998.

[6] Conlon, J. G. and Naddaf, A., Green’s functions for elliptic and parabolic equations with random coefficients, New York Journal of Mathematics, 6, pp. 153-225, 2000.

[7] Dagan, G., Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, pp. 183-215, 1987.

[8] Doob, J. L., Stochastic Processes, John Wiley & Sons, London, 1953.

[9] Eberhard J., Suciu, N. and C. Vamoş, On the self-averaging of dispersion for transport in quasi-periodic random media, J. Phys. A: Math. Theor., 40, pp. 597-610, doi: 10.1088/1751-8113/40/4/002, 2007.

[10] Fiori, A. and Dagan, G., Concentration fluctuations in aquifer transport: a rigorous first-order solution and applications, J. Contam. Hydrol., 45, pp. 139-163, 2000.

[11] Gardiner, C. W., Handbook of Stochastic Methods (for Physics, Chemistry and Natural Science), Springer, New York, 1985.

[12] Kloeden, P. E. and Platten, E., Numerical solutions of stochastic differential equations, Springer, Berlin, 1999.

[13] Jaekel, U. and Vereecken, H., Renormalization group analysis of macrodispersion in a directed random flow, Water Resour. Res., 33, pp. 2287-2299, 1997.

[14] LaBolle, E. M., Quastel, J., Fogg, G. E. and Gravner, J., Diffusion processes in composite media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coefficients, Water Resour. Res., 36(3), pp. 651-662, 2000.

[15] Le Doussal, P. and Machta, J., Annealed versus quenched diffusion coefficient in random media, Phys. Rev. B, 40(12), pp. 9427-9430, 1989.

[16] Lumley, J. L., An approach to the Eulerian-Lagrangian problem, J. Math. Phys., 3(2), pp. 309-312, 1962.

[17] Lundgren, T. S., Turbulent pair dispersion and scalar diffusion, J. Fluid Mech. 111, pp. 27-57, 1981.

[18] Monin, A. S. and Yaglom, A. M., Statistical Fluid Mechanics: Mechanics of Turbulence, MIT Press, Cambridge, M A, 1971.

[19] Sposito, G. and Dagan, G., Predicting solute plume evolution in heterogeneous porous formations, Water Resour. Res., 30(2), pp. 585-589, 1994.

[20] Suciu, N., Some Relations Between Microscopic and Macroscopic Modelling of Thermodynamic Processes (in Romanian), Ed. Univ. Piteşti, Appl. and Ind. Math. Series, No. 5, 2001.

[21] Suciu, N. and Vamoş, C., Effective diffusion in heterogeneous media, Internal Report ICG-IV.00303, Research Center Jülich, 2003.

[22] Suciu, N., Vamoş, C., Vanderborght, J., H. Hardelauf and Vereecken, H., Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi: 10.1029/2005WRR004546, 2006.

[23] Suciu, N., Vamoş, C. and J. Eberhard, Evaluation of the first-order approximations for transport in heterogeneous media, Water Resour. Res., 42, W11504, doi: 10.1029/2005WR004714, 2006.

[24] Suciu, N. and Vamoş, C., 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, 2007.

[25] Suciu N., Vamoş, C., Sabelfeld, K. and Andronache, C., Memory effects and ergodicity for diffusion in spatially correlated velocity fields, Proc. Appl. Math. Mech., 7, 2010015-2010016, doi: 10.1002/pamm.20070057, 2007.

[26] Suciu N., Vamos, C., Vereecken, H., Sabelfeld, K. and Knabner, P., Dependence on initial conditions, memory effects, and ergodicity of transport in heterogeneous media, Preprint No. 324, Institute of Applied Mathematics, Friedrich-Alexander University Erlangen-Nuremberg (available online at http://www.am.uni-erlangen.de/de/preprints2000.html), 2008.

[27] Suciu, N., Vamos, C., Vereecken, H., Sabelfeld, K. and Knabner, P., Memory effects induced by dependence on initial conditions and ergodicity of transport in heterogeneous media, Water Resour. Res., 44, W08501, doi: 10.1029/2007WR006740, 2008.

[28] Vamoş, C., Suciu, N., Vereecken, H., Vanderborht, J. and Nitzsche, O., Path decomposition of discrete effective diffusion coefficient, Internal Report ICG-IV.00501, Research Center Jülich, 2001.

[29] Vamoş, C., Suciu, N. and Vereecken, H., Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186(2), pp. 527–544, doi: 10.1016/S0021-9991(03)00073-1, 2003.

[30] van Kampen, N. G., Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1981.

[31] Yaglom, A. M., Correlation Theory of Stationary and Related Random Functions, Volume I: Basic Results, Springer-Verlag, New York, 1987.

[32] Zirbel, C. L., Lagrangian observations of homogeneous random environments, Adv. Appl. Prob., 33, pp. 810-835, 2001.

[33] Zwanzig, R., Memory effects and irreversible thermodynamics, Phys. Rev., 124(4), pp. 983-992, 1961.

jnaat,+Journal+manager,+2008-2-Suciu-Vamos-Vereecken-Sabelfeld-Knabner-15-10-29

ITÔ EQUATION MODEL FOR DISPERSION OF SOLUTES IN HETEROGENEOUS MEDIA

NICOLAE SUCIU $^{1}$ , CĂLIN VAMOŞ $^{2}$ , HARRY VEREECKEN $^{3}$ , KARL SABELFELD $^{4}$ and PETER KNABNER $^{5}$

Abstract

Transport processes in heterogeneous media such as ionized plasmas, natural porous media, and turbulent atmosphere are often modeled as diffusion processes in random velocity fields. Using the Itô formalism, we decompose the second spatial moments of the concentration and the equivalent effective dispersion coefficients in terms corresponding to various physical factors which influence the transport. We explicitly define "ergodic" dispersion coefficients, independent of the initial conditions and completely determined by local dispersion coefficients and velocity correlations. Ergodic coefficients govern an upscaled process which describes the transport at large tine-space scales. The non-ergodic behavior at finite times shown by numerical experiments for large initial plumes is explained by "memory terms" accounting for correlations between initial positions and velocity fluctuations on the trajectories of the solute molecules.

MSC 2000. Primary 82C70, 60K37; Secondary 60J70, 82C31.
Keywords. Itô equation, random fields, memory effects, ergodicity.

1. INTRODUCTION

Dispersion of solutes in spatially heterogeneous media, such as natural porous media, turbulent atmosphere, or ionized plasmas, is often described

as a two sale process, consisting of a local dispersion mechanism and random space functions describing the medium properties at larger scales. Let

X_{l}

be the

l

component of the trajectory of a solute molecule, governed by a dispersion process in a random velocity field

V

. To simplify matters, we consider only the diagonal components

l l

of the various second moments. For a fixed realization of the velocity, the second moment

s_{l l}

of the actual concentration is given by the dispersion, or the mean-square displacement of the molecules at a given time

\begin{matrix} (1) & s_{l l} = {⟨ {[X_{l} - {⟨ X_{l} ⟩}_{D X_{0}}]}^{2} ⟩}_{D X_{0}} \end{matrix}

The subscripts

D

and

X_{0}

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. A subscript

V

shall indicate averaging over the ensemble of realizations of the random velocity field and ensemble averages shall be noted by corresponding capital letters, e.g.

S_{l l} = {⟨ s_{l l} ⟩}_{V}

The dispersion (1) provides a measure for the spatial extension of the solute plume. Its dependence on time

s_{l l} \sim t^{α}

is commonly used to investigate whether the transport is diffusive,

α = 1

, or it has anomalous behavior,

α \neq 1

. Since it can be estimated, by either analytical approximations or numerical simulations, without solving the transport equations [23, 9] this quantity is particularly useful in investigations on pre-asymptotic transport regime, for which generally there are no close form solutions. The aim of this paper is to investigate the structure of various dispersion terms for transport processes described by the Itô equation. Numerical simulations for a typical case of contaminant transport in groundwater will be used to illustrate the time behavior of dispersion and its dependence of on initial conditions.

2. ITÔ FORMALISM

When the 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 [2]. For a given realization

V (x)

of a time independent random velocity field and a constant local dispersion coefficient

D

, the trajectories of the local advection-dispersion process starting at

t = 0

from

X_{0}

are solutions of the integral Itô equation

\begin{matrix} (2) & X_{l} (t) = X_{0 l} + \int_{0}^{t} V_{l} (X (t^{'})) d t^{'} + W_{l} (t) \end{matrix}

where

W_{l} (t) = \int_{0}^{t} d W_{l} (t^{'})

is a stochastic integral which defines a Wiener process with the properties

\begin{matrix} (3) & {⟨ W_{l} ⟩}_{D} = 0, {⟨ W_{l}^{2} ⟩}_{D} = 2 D t \end{matrix}

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}

of the Wiener process lead to the cancelation of the average of the Itô integral [11, 12]. For instance,

\begin{matrix} (4) & {⟨ \int_{0}^{t} X_{l} (t^{'}) d W_{l} (t^{'}) ⟩}_{D} = 0 . \end{matrix}

Taking into account that

{⟨ d W_{l} ⟩}_{D}^{2} = 2 D d t

, as follows from (3), and collecting terms of first-order in

d t

and

d W

, the differential of a function of time and of the

l

component of the Itô process (2),

f (t, x) = f (t, X_{l} (t))

, is computed as

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,

where the derivatives are evaluated at

x = X_{l} (t)

. This is a particular form of the stochastic chain rule known as the Itô formula (see e.g. [12], Section 3.3), which is mainly used in its integral form

\begin{aligned} (5) & f (t, X_{l} (t)) - f (0, X_{l} (0)) = \int_{0}^{t} \frac{\partial f (t^{'}, X_{l} (t^{'}))}{\partial t^{'}} d t^{'} \\ + \int_{0}^{t} [V_{l} (X_{l} (t^{'})) \frac{\partial f (t^{'}, X_{l} (t^{'}))}{d x} + D \frac{\partial^{2} f (t^{'}, X_{l} (t^{'}))}{\partial x^{2}}] d t^{'} \\ + \int_{0}^{t} \frac{\partial f (t^{'}, X_{l} (t^{'}))}{\partial x} d W (t^{'}) \end{aligned}

If the probability distribution of the Itô process (2) has a density

p (x, t)

with respect to the Lebesgue measure, the expectation

⟨ f (X (t)) ⟩_{D X_{0}} = \int f (x) p (x, t) d x

of some function with compact support

f (x)

can be written by using the Dirac

δ

distribution as

{⟨ \int f (x) δ [x - X (t)] d x ⟩}_{D X_{0}} = \int f (x) ⟨ δ [x - X (t)] ⟩_{D X_{0}} d x

. Hence, the probability density is given by

p (x, t) = ⟨ δ [x - X (t)] ⟩_{D X_{0}}

. Generalizing, one obtains a hierarchy of consistent joint

n

-times probability densities [20],

\begin{aligned} (6) & p (x_{1}, t_{1}; x_{2}, t_{2}; \dots; x_{n}, t_{n}) = \\ {⟨ δ [x_{1} - X (t_{1})] δ [x_{2} - X (t_{2})] \dots δ [x_{n} - X (t_{n})] ⟩}_{D X_{0}} \end{aligned}

The densities (6) are projections on cylindrical sets

{X (t_{1}) = x_{1}; X (t_{2}) = x_{2}; \dots X (t_{n}) = x_{n}}

of the probability measure on the space of Itô processes (2) starting from initial positions

X_{0}

[8, 33, 30, 11, 12]. In particular, the normalized concentration is defined by the one-time probability density,

c (x, t) = p (x, t)

Doob ([8], Chap. VI.3) has proved that with the Itô interpretation of the stochastic integral the solutions of the Itô equation, under suitable regularity properties of its coefficients, are trajectories of the diffusion process with densities of transition probabilities obeying a Fokker-Planck equation, and conversely, the trajectories of a diffusion process defined by the coefficients of the Fokker-Planck equation are solutions of the Itô equation. The normalized concentration

c

verifies the same Fokker-Planck equation. In the case of constant
coefficient

D

which we consider here, this is simply the advection-dispersion equation

\begin{matrix} (7) & \partial_{t} c + \nabla (V c) = D \nabla^{2} c \end{matrix}

The connection with the Stratonovich equation and Fick's law, for smoothly variable

D

, can be established by using the Itô formula ([11], Section 4.3.6). Itô formalism can also be used to account for variable porosity and discontinuous

D

[14]. The equivalence between (2) and (7) 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 (2) [12], or the GRW algorithm [29] used in simulations presented in this paper.

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 (6) for

n = 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 (7) ([11, p. 54). When the Lagrangian approach is extended to advection-dispersion processes, the velocity field is of the form

V + d W / d t

(e.g. [7, 10]). Since the derivative of the Wiener process does not exist ([11], p. 69), the transformation from

X_{0}

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 [17]. However, the smoothness and boundedness of the time derivative of the trajectory

X_{l} (t)

which permit the use of a Lagrangian description [19] can be ensured by modeling the local dispersion as a spatially correlated noise ([32], 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 4 . Even if the conceptual inconsistencies can be avoided in this way and equation (2) 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 [10]. 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 ([7], equation (3.11), [10], 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.

3. DISPERSION TERMS AND DISPERSION COEFFICIENTS

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
(2) and (3) one obtains the average of

X_{l}

, i.e. the

l

component of the plume center of mass,

\begin{matrix} (8) & {⟨ X_{l} ⟩}_{D X_{0}} (t) = {⟨ X_{0 l} ⟩}_{X_{0}} + \int_{0}^{t} {⟨ V_{l} (X (t^{'})) ⟩}_{D X_{0}} d t^{'} \end{matrix}

By applying the Itô formula (5) to the function

f (t, X_{l} (t) = [X_{l} (t) - {{⟨ X_{l} ⟩}_{D X_{0}} (t)]}^{2}

, with

X_{l} (0)) = X_{0 l}

, one obtains

\begin{aligned} {[X_{l} (t) - {⟨ X_{l} ⟩}_{D X_{0}} (t)]}^{2} = {[X_{0 l} - {⟨ X_{0 l} ⟩}_{X_{0}}]}^{2} \\ - 2 \int_{0}^{t} [X_{l} (t^{'}) - {⟨ X_{l} ⟩}_{D X_{0}} (t^{'})] {⟨ V_{l} (X (t^{'})) ⟩}_{D X_{0}} d t^{'} \\ + 2 \int_{0}^{t} [X_{l} (t^{'}) - {⟨ X_{l} ⟩}_{D X_{0}} (t^{'})] V_{l} (X (t^{'})) d t^{'} + 2 D t \\ (9) & + 2 \int_{0}^{t} [X_{l} (t^{'}) - {⟨ X_{l} ⟩}_{D X_{0}} (t^{'})] d W_{l} (t^{'}) \end{aligned}

After averaging with respect to

D

and

X_{0}

and using the property (4) of the Itô integral, the first term of (9) yields the deterministic initial second moment

S_{l l} (0) = {⟨ {[X_{0 l} - {⟨ X_{0 l} ⟩}_{X_{0}}]}^{2} ⟩}_{X_{0}}

, the second and fourth terms are averaged out, and one obtains the following explicit form of the dispersion (1):

\begin{aligned} s_{l l} (t) & = {⟨ {[X_{l} (t) - {⟨ X_{l} ⟩}_{D X_{0}} (t)]}^{2} ⟩}_{D X_{0}} \\ (10) & = S_{l l} (0) + 2 D t + 2 \int_{0}^{t} {⟨ [X_{l} (t^{'}) - {⟨ X_{l} ⟩}_{D X_{0}} (t^{'})] V_{l} (X (t^{'})) ⟩}_{D X_{0}} d t^{'} \end{aligned}

Useful information on dispersion of a non-singular initial concentration distribution (not concentrated at a point) is obtained if the relative displacement

{\tilde{X}}_{l} = X_{l} - X_{0 l}

is introduced into (10),

\begin{aligned} s_{l l} (t) & = S_{l l} (0) + 2 D t + 2 \int_{0}^{t} {⟨ [{\tilde{X}}_{l} (t^{'}) - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0}} (t^{'})] V_{l} (X (t^{'})) ⟩}_{D X_{0}} d t^{'} \\ (11) & + 2 \int_{0}^{t} {⟨ [X_{0 l} - {⟨ X_{0 l} ⟩}_{X_{0}}] V_{l} (X (t^{'})) ⟩}_{D X_{0}} d t^{'} \end{aligned}

As follows from (2),

{\tilde{X}}_{l}

is the solution of the Itô equation

\begin{matrix} (12) & {\tilde{X}}_{l} (t) = \int_{0}^{t} V_{l} (\tilde{X} (t^{'}) + X_{0}) d t^{'} + W_{l} (t) \end{matrix}

with the drift coefficient defined by

V_{l} (\tilde{X} (t) + X_{0}) = V_{l} (X (t), t)

and initial condition

{\tilde{X}}_{l} = 0

. From (12) and (3) we can see that the last term of (11) is twice the spatial correlation between initial positions measured from the center of mass of the initial plume

X_{0 l} - {⟨ X_{0 l} ⟩}_{X_{0}}

and the mean relative displacements

{⟨ {\tilde{X}}_{l} ⟩}_{D}

, which can also be expressed by twice the time integral of the spatial
correlation between initial positions and mean velocities

{⟨ V_{l} (X (t^{'})) ⟩}_{D}

sampled on the trajectories of the Itô process starting from initial positions

X_{0}

\begin{aligned} m_{l l} (t) & = 2 {⟨ [X_{0 l} - {⟨ X_{0 l} ⟩}_{X_{0}}] {⟨ {\tilde{X}}_{l} ⟩}_{D} (t) ⟩}_{X_{0}} \\ (13) & = 2 \int_{0}^{t} {⟨ [X_{0 l} - {⟨ X_{0 l} ⟩}_{X_{0}}] {⟨ V_{l} (X (t^{'})) ⟩}_{D} ⟩}_{X_{0}} d t^{'} \end{aligned}

Since the time behavior of

m_{l l}

defined by (13) tell us how much and how long the solute molecules "remember" their starting points, it has been called "memory term" 25, 24, 26, 27.

Considering the Itô process (12), applying the Itô formula (5) to the function

f (t, {\tilde{X}}_{l} (t) = {[{\tilde{X}}_{l} (t) - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0}} (t)]}^{2}

, and taking the expectation, one finds that the dispersion of the relative displacements

{⟨ {[{\tilde{X}}_{l} - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0}}]}^{2} ⟩}_{D X_{0}}

is just the sum of the second and third terms of (11). Hence, the total dispersion is

\begin{matrix} (14) & s_{l l} (t) = S_{l l} (0) + {⟨ {[{\tilde{X}}_{l} - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0}}]}^{2} ⟩}_{D X_{0}} (t) + m_{l l} (t) \end{matrix}

(A similar relation holds for the case

D = 0

as well [19, 24]). Further, from the identity

{⟨ {[{\tilde{X}}_{l} - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0}}]}^{2} ⟩}_{D X_{0}} = {⟨ {[{\tilde{X}}_{l} - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0} V}]}^{2} ⟩}_{D X_{0}} - {[{⟨ X_{l} ⟩}_{D X_{0}} - {⟨ X_{l} ⟩}_{D X_{0} V}]}^{2}

one obtains

\begin{matrix} (15) & s_{l l} (t) = S_{l l} (0) + {\tilde{x}}_{l l} (t) + m_{l l} (t) - r_{l l} (t) \end{matrix}

where the new terms

{\tilde{x}}_{l l}

and

r_{l l}

are defined by

\begin{aligned} {\tilde{x}}_{l l} (t) & = {⟨ {[{\tilde{X}}_{l} - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0} V}]}^{2} ⟩}_{D X_{0}} \\ (16) & = 2 D t + 2 \int_{0}^{t} {⟨ [{\tilde{X}}_{l} (t^{'}) - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0} V} (t^{'})] u_{l} (X (t^{'}), t^{'}) ⟩}_{D X_{0}} d t^{'} \\ r_{l l} (t) & = {[{⟨ X_{l} ⟩}_{D X_{0}} - {⟨ X_{l} ⟩}_{D X_{0} V}]}^{2} \\ (17) & = \int_{0}^{t} \int_{0}^{t} {⟨ u_{l} (X (t^{'}), t^{'}) ⟩}_{D X_{0}} {⟨ u_{l} (X (t^{''}), t^{''}) ⟩}_{D X_{0}} d t^{'} d t^{''} \end{aligned}

and

u_{l} (X (t), t) = V_{l} (X (t)) - {⟨ V_{l} (X (t)) ⟩}_{D X_{0} V}

is the velocity fluctuation with respect to the total average of the velocity sampled on trajectories

{⟨ V_{l} (X (t)) ⟩}_{D X_{0} V}

. Hereafter the subscript

V

denotes the average over the ensemble of velocity realizations. The explicit form of (17) follows directly from Itô equation (12), while (16) is derived from the Itô formula (5) applied to the function

f (t, {\tilde{X}}_{l} (t) = {[{\tilde{X}}_{l} (t) - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0} V} (t)]}^{2}

. Note that because the average velocity is independent of

X_{0}

, the memory term (13) can also be expressed by correlations initial positions-velocity fluctuations.

The relations (14) and (15) are simple identities and are strictly equivalent to the definition (1) of the dispersion. They are independent of the model for the trajectories of the transport process and apply for either Itô diffusions or for discrete/continuous time Markov chains [26, 27]. The Itô model discussed in this paper leads to the explicit forms (13), (16), and (17) of dispersion terms, which, as shown in the next Section can be expressed as averages with
respect to the probability densities (6) of the Itô process and lead to useful first order approximations.

Assuming all necessary joint measurability conditions which allow permutations of averages [32, the average over the ensemble of velocity realizations

S_{l l} = {⟨ s_{l l} ⟩}_{V}

can be expressed as [25, 24]

\begin{matrix} (18) & S_{l l} (t) = S_{l l} (0) + {⟨ X_{l l} ⟩}_{X_{0}} (t) + Q_{l l} (t) + M_{l l} (t) - R_{l l} (t) \end{matrix}

where

X_{l l} = {⟨ {[{\tilde{X}}_{l} - {⟨ {\tilde{X}}_{l} ⟩}_{D V}]}^{2} ⟩}_{D V}

is the one-particle dispersion and

Q_{l l} = {⟨ {[{⟨ {\tilde{X}}_{l} ⟩}_{D V} - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0} V}]}^{2} ⟩}_{X_{0}}

is the spatial variance of one-particle center of mass

{⟨ {\tilde{X}}_{l} ⟩}_{D V}

, which are related to the ensemble average

{\tilde{X}}_{l l} = {⟨ {\tilde{x}}_{l l} ⟩}_{V}

by the relation

\begin{matrix} (19) & {\tilde{X}}_{l l} (t) = {⟨ X_{l l} ⟩}_{X_{0}} (t) + Q_{l l} (t) \end{matrix}

The other two terms in (18) are respectively

M_{l l} = {⟨ m_{l l} ⟩}_{V}

, the mean memory term, and

R_{l l} = {⟨ r_{l l} ⟩}_{V}

, the variance of the center of mass.

It is readily to check that the sum of the first four terms of (18) gives the dispersion over the total ensemble consisting of realizations of the Wiener process, initial positions, and velocity realizations,

\begin{matrix} (20) & Σ_{l l} = S_{l l} (0) + {⟨ X_{l l} ⟩}_{X_{0}} + Q_{l l} (t) + M_{l l} = {⟨ {[X_{l} - {⟨ X_{l} ⟩}_{D X_{0} V}]}^{2} ⟩}_{D X_{0} V} \end{matrix}

The "ensemble dispersion" (20) is related to the ensemble average

S_{l l}

of the dispersion in single realizations of the random velocity field and to the variance of the center of mass

R_{l l}

by the identity

\begin{matrix} (21) & S_{l l} (t) = Σ_{l l} (t) - R_{l l} (t) \end{matrix}

We note that the identity (21) is not subject to restrictions permitting permutations of averages and has a single-realization correspondent,

\begin{matrix} (22) & s_{l l} (t) = σ_{l l} (t) - r_{l l} (t) \end{matrix}

The new dispersion term occurring in (22) is given, according to (15), by

σ_{l l} = S_{l l} (0) + {\tilde{x}}_{l l} + m_{l l}

Dividing both sides of (15) by

2 t

one obtains an equivalent relation which describes the structure of the effective dispersion coefficients

d_{l l}^{eff} = S_{l l} / (2 t)

for a given realization of the velocity field:

\begin{matrix} (23) & d_{l l}^{eff} (t) = \frac{S_{l l} (0)}{2 t} + D + d_{l l}^{adv} (t) + d_{l l}^{mem} (t) - d_{l l}^{cm} (t) \end{matrix}

The last three terms in (23) are respectively: the "advective coefficient"

d_{l l}^{adv} = {\tilde{x}}_{l l} / (2 t) - D

, associated to velocity fluctuations along trajectories, the "memory coefficient"

d_{l l}^{mem} = m_{l l} / (2 t)

, describing correlations between velocity fluctuations and initial positions, and the "center of mass coefficient"

d_{l l}^{cm} = r_{l l} / (2 t)

, a contribution of the fluctuating center of mass velocity. An identical relation holds for the corresponding ensemble averages. Among them, the mean advective coefficient

D_{l l}^{adv} = {⟨ d_{l l}^{adv} ⟩}_{V} = {\tilde{X}}_{l l} / (2 t) - D

plays a special role in modeling diffusion in random fields. If the random velocity field is statistically homogeneous and the solutions of the Itô equation (2) are pathwise
unique, then the one-particle center of mass

{⟨ {\tilde{X}}_{l} ⟩}_{D V}

and dispersion

X_{l l}

are independent of initial positions

X_{0}

, and the mean memory term

M_{l l}

vanishes [27]. Then, according to (19),

{\tilde{X}}_{l l} = X_{l l}

and one defines an "ergodic" coefficient

D_{l l}^{erg} = X_{l l} / (2 t) [22

. Under the above conditions, the ergodic coefficient is the sum between the local dispersion coefficient and the mean advective coefficient,

\begin{matrix} (24) & D_{l l}^{erg} (t) = D + D_{l l}^{adv} (t) \end{matrix}

Known mathematical results show that under suitable limiting conditions the ergodic coefficient (24) tends to a "macrodispersion coefficient" governing an up-scaled process for the ensemble averaged concentration (e.g. 22, equation (4)).

If the velocity field has finite correlation range then it is ergodic [31, a property which essentially determines the dispersion in case of large supports of the initial concentration distributions [15, 26, 27. Heuristically, under the ergodic assumption one can replace space averages by ensemble averages, the second term of (14) estimates the one-particle dispersion

X_{l l}

, and one obtains

\begin{matrix} (25) & s_{l l} (t) \approx S_{l l} (0) + X_{l l} (t) + m_{l l} (t) . \end{matrix}

Accordingly, the dispersion coefficient for single realizations of the velocity fields (23) can be approximated by

\begin{matrix} (26) & d_{l l}^{eff} (t) \approx \frac{S_{l l} (0)}{2 t} + D_{l l}^{erg} (t) + d_{l l}^{mem} (t) . \end{matrix}

Thus, if one proves that in the limit of large times velocity fluctuations decorrelate from initial positions of the diffusing molecules, then the memory terms become negligible and the single realizations of transport are governed by the memory-free dispersion

X_{l l}

, or equivalently by the ergodic dispersion coefficients

D_{l l}^{erg}

. Such a behavior has been indicated by recent numerical investigations [27.

Since, as shown by (24), the ergodic behavior is basically governed by the advective coefficients, we analyze them here in more detail. The structure of the advective coefficient

d_{l l}^{adv} = {\tilde{x}}_{l l} / (2 t) - D

is that of the second term of (16),

2 t d_{l l}^{adv} (t) = 2 \int_{0}^{t} {⟨ [{\tilde{X}}_{l} (t^{'}) - {⟨ {\tilde{X}}_{l} ⟩}_{D X_{0} V} (t^{'})] u_{l} (X (t^{'}), t^{'}) ⟩}_{D X_{0}} d t^{'}

which, using the Itô equation (12), can be rewritten as

\begin{aligned} (27) & 2 t d_{l l}^{adv} (t) & = \int_{0}^{t} \int_{0}^{t} {⟨ u_{l} (X (t^{'}), t^{'}) u_{l} (X (t^{''}), t^{''}) ⟩}_{D X_{0}} d t^{'} d t^{''} \\ + 2 \int_{0}^{t} {⟨ W_{l} (t^{'}) u_{l} (X (t^{'}), t^{'}) ⟩}_{D X_{0}} d t^{'} \end{aligned}

In the case of homogeneous random velocity fields, which have vanishing mean fluctuations, the ensemble average of the second term of (27) vanishes. Hence,
the mean advective coefficient

D_{l l}^{adv}

is solely determined by velocity correlations and does not contain contributions from the cross-correlation between velocity fluctuations and Wiener process [26]. We show in the following that the single-realization coefficient

d_{l l}^{adv}

has a similar structure. For that, we consider the solution of the Itô equation (2) constructed by successive approximations

\begin{matrix} (28) & X_{l}^{(k)} (t) = X_{0 l} + \int_{0}^{t} V_{l} (X^{(k - 1)} (t^{'})) d t^{'} + W_{l} (t) \end{matrix}

If the drift coefficient

V_{l} (x)

of the Itô equation is Lipschitz continuous, then

X_{l}^{(k)}

converges in mean-square for

k ⟶ \infty

to the exact solution

X_{l}

([12, p. 133). The proof of convergence involves expectations of the squares of the terms of (28) and does not require a prescribed Wiener process. So, we can construct weak solutions, for instance by Euler scheme, of the equations for successive iterations

k

of (28) with different Wiener processes independent of each other. In practice, this is always the case if the seed of the random number generator used in the weak Euler scheme is not resetted at the initial value after each iteration. This procedure yields a process that is not generally samplepath equivalent to

X_{l}

but has the same probability law ([12, p. 144). Under these conditions

u_{l} (X^{(k - 1)} (t^{'}), t^{'})

is independent of

W_{l} (t^{'})

at each iteration step, and since the successive approximations converge the second term of (27) vanishes. Thus, the advective coefficients in single-realizations of the velocity field are completely determined by correlations of velocity field sampled by the diffusion process,

\begin{matrix} (29) & 2 t d_{l l}^{adv} (t) = {\tilde{x}}_{l l} - 2 D t = \int_{0}^{t} \int_{0}^{t} {⟨ u_{l} (X (t^{'}), t^{'}) u_{l} (X (t^{''}), t^{''}) ⟩}_{D X_{0}} d t^{'} d t^{''} \end{matrix}

This result was also obtained within a more laborious approach based on Fokker-Planck equation and Kolmogorov's definitions of drift and diffusion coefficients, under the assumption that the first two moments of the concentration field are finite at finite times [21]. Numerical simulations also show that contributions arising from cross-correlations between velocity fluctuations and diffusion jumps, replacing the Wiener process in weak schemes, are of the order of the numerical precision [28].

4. FIRST-ORDER APPROXIMATIONS

To write the dispersion terms (13), (16), and (17) as moments of the joint probabilities (6) of the Itô process, we define the projection of the continuous Eulerian velocity

V_{l} (x)

on the trajectory

X (t)

by using the Dirac measure

δ

\begin{matrix} (30) & V_{l} (X (t)) = \int V_{l} (x) δ [x - X (t)] d x \end{matrix}

We note by

c (x_{0}) = c (x, 0)

the initial normalized concentration distribution given by (6) for

n = 1

and

t = 0

, express the two-times joint probability densities as transition probabilities multiplied by

c (x_{0}), p (x_{0}; x, t) = p (x, t ∣ x_{0}) c (x_{0})

, and the actual concentration as

c (x, t) = \int_{0}^{t} p (x_{0}; x^{'}, t^{'}) d x_{0}

[11. Considering the form given by (29) of the dispersion term (16) and using (30) one obtains

\begin{aligned} {\tilde{x}}_{l l} (t) = & 2 D t + 2 \int_{0}^{t} \int_{0}^{t} d t^{'} d t^{''} \int c (x_{0}) d x_{0} \\ (31) & (31) & \iint u_{l} (x^{'}, t^{'}) u_{l} (x^{''}, t^{''}) p (x^{'}, t^{'} ∣ x^{''}, t^{''}) p (x^{''}, t^{''} ∣ x_{0}) u_{l l} (x^{'}, x^{''}) d x^{'} d x^{''} \end{aligned}

Similarly, for the memory term (13) and center of mass fluctuations term (17), we have

\begin{aligned} (32) & m_{l l} (t) = 2 \int_{0}^{t} d t^{'} \int c (x_{0}) d x_{0} (x_{0 l} - \int x_{0 l}^{'} c (x_{0}^{'}) d x_{0}^{'}) \int u_{l} (x^{'}, t^{'}) p (x^{'}, t^{'} ∣ x_{0}) d x^{'} \\ r_{l l} (t) = \int_{0}^{t} \int_{0}^{t} d t^{'} d t^{''} \iint c (x_{01}) c (x_{02}) d x_{01} d x_{02} \\ (33) & \iint u_{l} (x^{'}, t^{'}) u_{l} (x^{''}, t^{''}) p (x^{'}, t^{'} ∣ x_{01}) p (x^{''}, t^{''} ∣ x_{02}) d x^{'} d x^{''} \end{aligned}

The relations (31-33) achieve the explicit form of the decomposition (15) of the dispersion. Such relations were also derived for the more general case of space-time dependent velocity and local dispersion coefficient, 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 [21].

With the exception of some solvable problems, as for instance in the case of transport in perfectly stratified flows [1, 5], the ensemble average of (31-33) cannot be computed exactly. The difficulty arises from the highly non-linear dependence of probability densities in these relations on the velocity

V_{l} (X (t))

sampled on trajectories which, in their turn, are determined by the velocity field [1]. 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 (see [18, p. 571). Even in the advective case (

D = 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 [16]. In principle, the problem of ensemble averaging can be handled by iterating indefinitely the transport equation around a particular solution independent of velocity realizations [3] 1, 13.

Iterations of the advection-dispersion equation (7) are usually done in La-place-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 [5]. Alternatively, the solution

X (t)

of the Itô equation (2) can be approximated
by successive approximations (28), starting with some trajectory

X^{(0)} (t)

that is independent of velocity statistics. The second iteration is already highly complex. In this approximation, the argument of

u_{l}

in (16), (17), and (29) has to be replaced by

X^{(1)}

. But, as shown by (28),

X^{(1)}

depends on the values of the velocity

V_{l}

X^{(0)} (t^{'})

for all times between 0 and

t

. Hence, an infinite-dimensional probability distribution is needed to compute ensemble averages.

The approximation by the first iteration of the Itô equation,

X_{l} (t) \approx X_{l}^{(1)} (t)

, which is equivalent to the most-used first-order approximations in velocity fluctuations [23], simplifies matters considerably. Now, as follows from (28), the argument of

V_{l}

is replaced by

X^{(0)}

, which is independent of the realization of the velocity field. The probability densities (6) computed from the trajectories

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. Then, the mean of the velocity fluctuations on trajectories

{⟨ u_{l} ⟩}_{V} = 0

and, accordingly to (32), the mean memory term

M_{l l} = {⟨ m_{l l} ⟩}_{V}

vanishes. Further, by taking the ensemble average of (31) one obtains the first-order approximation

\begin{aligned} {\tilde{X}}_{l l} (t) = 2 D t + & \int_{0}^{t} \int_{0}^{t} d t^{'} d t^{''} \int c (x_{0}) d x_{0} \\ (34) & \iint p (x^{'}, t^{'} ∣ x^{''}, t^{''}) p (x^{''}, t^{''} ∣ x_{0}) u_{l l} (x^{'}, x^{''}) d x^{'} d x^{''} \end{aligned}

where

u_{l l} (x^{'}, x^{''}) = {⟨ u_{l} (x^{'}) u_{l} (x^{''}) ⟩}_{V}

is the Eulerian velocity correlation. The approximated variance of the center of mass is similarly obtained by ensemble averaging (33),

\begin{aligned} R_{l l} (t) = \int_{0}^{t} \int_{0}^{t} d t^{'} d t^{''} & \iint c (x_{01}) c (x_{02}) d x_{01} d x_{02} \\ (35) & \iint p (x^{'}, t^{'} ∣ x_{01}) p (x^{''}, t^{''} ∣ x_{02}) u_{l l} (x^{'}, x^{''}) d x^{'} d x^{''} \end{aligned}

If the initial solution

X^{(0)}

is a diffusion in the ensemble mean velocity field

U = ⟨ V ⟩_{V}

, then the probability densities in (34-35) are those of a Gaussian process of mean

x_{0} + U t

and dispersion

2 D t

. Since the Gaussian transition probability and the Eulerian correlation of the homogeneous velocity field are invariant to space translations, (34) is independent of

x_{0}

. It is also easily to see that for homogeneous velocity fields the first order approximation of the term

Q_{l l}

from (19) vanishes. Hence, as it follows from (19),

{\tilde{X}}_{l l}

given by (34) is a first-order approximation of the one-particle dispersion

X_{l l}

. Due to the simplicity of this approximation which uses Gaussian probability densities, the averages (

34 - 35

), as well as the standard deviations of the single realization terms

{\tilde{x}}_{l l}

and

r_{l l}

, can be computed, at least by numerical integration, for given
two- and four-points Eulerian correlations. We note that such computations can be achieved in space-time domains, avoiding the tedious technique based Laplace and Fourier representations [3. The latter are however essential tools for higher order approximations 13 .

For the advection-dominated transport problem considered in simulations presented in the next Section, consistent asymptotic expansions of (2) are obtained when one starts with the solution

X^{(0)}

given by the deterministic trajectory

X_{0} + U t

of the ensemble mean velocity [23]. Then, as follows from (6), the conditional probabilities from (34-35) degenerate to Dirac delta measures,

p (x^{'}, t^{'} ∣ x_{0}) = δ [x^{'} - (x_{0} + U t^{'})]

and

p (x^{'}, t^{'}; x^{''}, t^{''} ∣ x_{0}) = δ [x^{'} - (x_{0} + U t^{'})] δ [x^{''} - (x_{0} + U t^{''})]

and one obtains

\begin{matrix} (36) & {\tilde{X}}_{l l} (t) = 2 D t + \int_{0}^{t} \int_{0}^{t} d t^{'} d t^{''} \int c (x_{0}) u_{l l} (x_{0} + U t^{'}; x_{0} + U t^{''}) d x_{0} \\ (37) & R_{l l} (t) = \int_{0}^{t} \int_{0}^{t} d t^{'} d t^{''} \iint c (x_{01}) c (x_{02}) u_{l l} (x_{01} + U t^{'}; x_{02} + U t^{''}) d x_{01} d x_{02} \end{matrix}

Obviously, (36) is also a first-order approximation of the one-particle dispersion

X_{l l}

. Since the velocity field is assumed to be statistically homogeneous, the one-particle velocity correlation from (36) depends only on time differences

t^{''} - t^{'}

. The two-particle correlation from (37) instead depends on both

t^{''} - t^{'}

and the separation lag

x_{02} - x_{01}

For either the inconsistent approximation (34-35) or the consistent one, (36-37), the ensemble average dispersion takes on the simpler form

\begin{matrix} (38) & S_{l l} (t) = S_{l l} (0) + X_{l l} (t) - R_{l l} (t) \end{matrix}

which was the starting point in many applications to transport in natural porous media presented in the hydrological literature (see references in [26]). For small velocity fluctuations, the approximations (34) and (36) of the memo-ry-free one-particle dispersion

X_{l l}

are very close to each other at all times ([23], Figure 2). The first-order approximations (34-35) also give reasonably good estimations for the ensemble average

S_{l l}

but they yield considerable overestimations, by

40 %

80 %

, of the fluctuations of the single-realization dispersion

s_{l l}

at early times ([23], Figure 4). Therefore, such approximations might not be accurate enough for investigations of memory effects caused by terms

m_{l l} (13)

, which manifest themselves in the pre-asymptotic regime.

5. NUMERICAL SIMULATIONS

A numerical experiment was conducted for a typical situation of contaminant transport in groundwater. We considered an isotropic two-dimensional aquifer system. The log-hydraulic conductivity was modeled as a statistically homogeneous random space function with exponentially decaying correlation, small variance

σ^{2} = 0.1

, and correlation length

λ = 1 m

. The realizations of the velocity field were generated numerically, as first order approximations in

Fig. 1. Longitudinal memory terms for different initial conditions and

P e

numbers.

Fig. 2. Transverse memory terms for different initial conditions and

P e

numbers.

σ

, by a superposition of 6400 periodic modes [9]. A constant local dispersion coefficient

D = 0.01 m^{2} /

day was chosen. For the fixed mean flow velocity

U = 1 m /

day the Péclet number got a typical value

P e = U λ / D = 100

. We obtained ensembles of transport realizations with the "global random walk algorithm" 29, namely by tracking simultaneously in every simulation

10^{10}

particles that were initially uniformly distributed in rectangular domains

L_{1} λ \times L_{2} λ

, (for details on the implementation of the numerical method see [22], Appendix A. To reduce the statistical oscillations (for both means and standard deviations of the memory terms) to less than

D t

at all times, we simulated the dispersive transport in

R = 1024

realizations of the velocity field for every initial condition. For narrow slab sources, statistical estimations for

P e = \infty

were also done by simulating advective transport using the same ensemble of velocities. To avoid artificial trapping phenomena ([22, Appendix

B 1)

, these simulations were limited to a period equal to

100 U t / λ

In the following we present some numerical results illustrating the memory effects on solute dispersion in the pre-asymptotic regime and the attainment of an ergodic behavior in the long time limit.

Comparisons for different dimensions of slab sources oriented across the

l

direction have shown that the difference between the second moment

Σ_{l l}

of the ensemble average concentration field and the initial second moment

S_{l l} (0)

is practically independent of the length of the initial slab source ([26], Figure 2). This indicates that, for the case of slab sources perpendicular to the

l

-direction and statistically homogeneous velocity fields, the relation (20) simplifies and the one-particle dispersion can be approximated by

X_{l l} \approx Σ_{l l} - S_{l l} (0)

. Using this approximation, we estimated the longitudinal and transverse one-particle dispersion,

X_{11}

and

X_{22}

, from

Σ_{11}

computed for the slab source (

λ, 100 λ

) and from

Σ_{22}

for the slab

(100 λ, λ)

. Further, since the spatial variance of the one-particle center of mass

Q_{l l}

was found to be negligible small ([24], Figure 2 ), we estimated the ensemble average memory terms for slabs oriented along

l

-axis, according to (20), by

M_{l l} \approx Σ_{l l} - S_{l l} (0) - X_{l l}

and the corresponding standard deviations by

S D (m_{l l}) \approx S D (σ_{l l})

. Actually, for the case of large slabs
considered here

Σ_{l l} \approx S_{l l}

and

S D (σ_{l l}) \approx S D (s_{l l})

([26, Figures 2, 3, 5) and both the mean and standard deviations of the memory terms can be computed with the aid of the ergodic estimate (25) [27]. The results for longitudinal and transverse memory coefficients are presented in Figures 1 and 2, respectively. The thin lines in Figures 1 and 2 correspond to

[M_{l l} \pm S D (m_{l l}) / R^{1 / 2}] / (2 D t)

, where

R = 1024

is the number of realizations used for statistical estimates.

The average memory coefficients

M_{11}

for longitudinal slabs (Figure 1) and

M_{22}

for transverse slabs (Figure 2) have significant non-vanishing values at finite times and approach to zero after hundreds of dimensionless times

U t / λ

. The decrease of the standard deviations indicates that the memory terms

m_{l l}

converge to zero for

t \to \infty

, in the mean square limit. The comparison with the case

P e = \infty

shows that there are no important differences with respect to the purely advective transport. For the square source

(100 λ, 100 λ)

the average memory terms are smaller than the local dispersion

2 D t

. In the pre-asymptotic regime the standard deviations are considerably large:

S D [m_{l l} (100, 100)] / (2 D t) \sim 8

for

l = 1, 2, S D [m_{22} (1, 100)] / (2 D t) \sim 20

, and

S D [m_{11} (100, 1)] / (2 D t) \sim 100

Apparently, the occurrence of non-vanishing mean memory terms disagrees with theoretical considerations and first order approximations which predict their cancelation for transport in statistically homogeneous velocity fields. We note that the fact that first order approximations presented in Section (4) cancel the mean memory terms is not surprising. The memory terms (32) depend only linearly on velocity fluctuations, unlike the other two dispersion terms, (31) and (33), which have a quadratic dependence. Hence, contributions of the order of velocity variance to the mean memory terms can still arise if a second iteration step is performed.

In fact, there are no sufficient theoretical arguments that mean memory terms should vanish for the velocity fields considered in the numerical experiment. The ensemble average of (13) vanish if and only if the one-particle statistics is independent of initial conditions. This can be the case if the velocity field is statistically homogeneous and if there are unique solutions of the transport equations [6, 32]. A prerequisite for the existence of pathwise unique strong solutions of Itô equation (2) is the Lipschitz continuity of the sample velocity field ([12], Theorem 4.5.3). Or, the sufficient condition for sample continuity, that is the existence of the second derivative of the correlation function at the origin [31], is not fulfilled for exponential correlations as those of the velocity fields generated numerically. Since this does not exclude the existence of continuous samples, we estimated numerically the Lipschitz constant

K = | V (x + Δ x) - V (x) | / Δ x

for velocity samples

V (x)

generated with 6400 periodic modes on a straight line (Figure 3) and for a sinusoidal velocity (Figure 4). The latter is continuous even for exponential correlations ([31], p. 67), and Figure 4 indicates indeed the existence of a Lipschitz constant. But for the noisy velocity sample from Figure 3 we found that

K

increases with

Fig. 3. Estimation of Lipschitz constant for noisy velocity sample with 6400 modes.

Fig. 4. Estimation of Lipschitz constant for smooth sinusoidal velocity sample.

decreasing

Δ x

. Thus, very likely at the limit of an infinite number of modes the resulting velocity field is not Lipschitz continuous.

6. DISCUSSION AND CONCLUSIONS

Since the Itô formalism leads to a consistent model of dispersion in heterogeneous media, we propose it as an alternative to questionable attempts, discussed at the end of Section 2, aiming at extending the Lagrangian description to diffusion processes.

Using the Itô model, we obtained an explicit decomposition of the second moments and dispersion coefficients in terms accounting for different physical factors that govern the process (Section 3) and we derived first order approximations of such dispersion quantities (Section 4). In particular, the influence of the location, shape, and dimensions of the support of the initial concentration in case of space variable coefficients of the transport equations is quantified by memory terms. The physical explanation of the memory terms consists of correlations between initial positions and the velocity experienced by solutes molecules when diffusion takes place in a space variable velocity field.

The results presented in Figures 1 and 2 supply a numerical evidence that velocities sampled on the trajectories of the Itô diffusion (2) decorrelate from initial positions and the memory terms (13) vanish after a sufficiently large time. Then, according to (25-26), the physically observable transport process, corresponding to single replicates of the random velocity field, is governed by the one-particle dispersion and by the ergodic dispersion coefficient.

At finite times, however, memory terms can be large and can produce nonergodic behavior with respect to the one-particle dispersion. The persistence of such memory effects, over hundreds of days in the conditions of the numerical experiment presented above, should be carefully considered when ergodic dispersion coefficients are used in practice. The definition (13) and the Cauchy-Schwartz inequality

m_{l l}^{2} (t) \leq 4 S_{l l} (0) {⟨ {⟨ {\tilde{X}}_{l} ⟩}_{D}^{2} ⟩}_{X_{0}}

show that memory terms increase with the increase of the dimension, and of the corresponding second moment

S_{l l} (0)

, of the initial concentration distribution. Numerical
results indicate that the increase is more pronounced for asymmetric initial concentration distributions ([27], Figure 2).

Non-vanishing mean memory terms shown in Figures 1 and 2 warn that the mean behavior of the pre-asymptotic dispersion is not always described by the well known relation (38). This happens for instance when the samples of the velocity fields, as for instance those with exponential correlations analyzed in Figures 3 and 4, do not fulfill the smoothness requirements which ensure the existence of strong unique solutions of the transport equations. In such circumstances, that are incompatible with analytical results derived from Itô or Fokker-Planck equations, the transport process can still be analyzed with general relations like (

14 - 15

) or (

21 - 23

). It is nevertheless noticeable that even in the worse case of velocity fields with exponential correlations, there is a large time behavior described by the ergodic coefficient (24) [22].

Acknowledgement. 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-1196, NATO Collaborative Linkage Grant ESP.NR.CLG 981426, and RFBR Grant 06-01-00498. The authors thank C. Andronache, J. G. Conlon, G. Grün, A. Millet, G. Keller, P. Kramer, O. Kurbanmuradov, F. Radu, U. Rüde, M. Trefry, W. Wagner, and C. L. Zirbel for stimulating discussions.

REFERENCES

[1] Balescu, R., Wang, H.-D. and Misguich, J. H., Langevin equation versus kinetic equation: Subdiffusive behavior of charged particles in a stochastic magnetic field, Phys. Plasmas 1(12), pp. 3826-3842, 1994.
[2] Bhattacharya, R. N. and Gupta, V. K., A Theoretical Explanation of Solute Dispersion in Saturated Porous Media at the Darcy Scale, Water Resour. Res., 19, pp. 934-944, 1983.
[3] Bouchaud, J.-P. and Georges, A., Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep., 195, pp. 127-293, 1990.
[4] Brouwers, J. J. H., On diffusion theory in turbulence, J. Eng. Mat., 44, pp. 277-295, 2002.
[5] Compte, A. and Cáceres, M. O., Fractional dynamics in random velocity fields, Phys. Rev. Lett., 81(15), pp. 3140-3143, 1998.
[6] Conlon, J. G. and Naddaf, A., Green's functions for elliptic and parabolic equations with random coefficients, New York Journal of Mathematics, 6 , pp. 153-225, 2000.
[7] Dagan, G., Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, pp. 183-215, 1987.
[8] Doob, J. L., Stochastic Processes, John Wiley & Sons, London, 1953.
[9] Eberhard J., Suciu, N. and C. Vamoş, On the self-averaging of dispersion for transport in quasi-periodic random media, J. Phys. A: Math. Theor., 40, pp. 597-610, doi: 10.1088/1751-8113/40/4/002, 2007.
[10] Fiori, A. and Dagan, G., Concentration fluctuations in aquifer transport: a rigorous first-order solution and applications, J. Contam. Hydrol., 45, pp. 139-163, 2000.
[11] Gardiner, C. W., Handbook of Stochastic Methods (for Physics, Chemistry and Natural Science), Springer, New York, 1985.
[12] Kloeden, P. E. and Platten, E, Numerical solutions of stochastic differential equations, Springer, Berlin, 1999.
[13] Jaekel, U. and Vereecken, H., Renormalization group analysis of macrodispersion in a directed random flow, Water Resour. Res., 33, pp. 2287-2299, 1997.
[14] LaBolle, E. M., Quastel, J., Fogg, G. E. and Gravner, J., Diffusion processes in composite media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coefficients, Water Resour. Res., 36(3), pp. 651-662, 2000.
[15] Le Doussal, P. and Machta, J., Annealed versus quenched diffusion coefficient in random media, Phys. Rev. B, 40(12), pp. 9427-9430, 1989.
[16] Lumley, J. L., An approach to the Eulerian-Lagrangian problem, J. Math. Phys., 3(2), pp. 309-312, 1962.
[17] Lundgren, T. S., Turbulent pair dispersion and scalar diffusion, J. Fluid Mech. 111, pp. 27-57, 1981.
[18] Monin, A. S. and Yaglom, A. M., Statistical Fluid Mechanics: Mechanics of Turbulence, MIT Press, Cambridge, M A, 1971.
[19] Sposito, G. and Dagan, G., Predicting solute plume evolution in heterogeneous porous formations, Water Resour. Res., 30(2), pp. 585-589, 1994.
[20] Suciu, N., Some Relations Between Microscopic and Macroscopic Modelling of Thermodynamic Processes (in Romanian), Ed. Univ. Piteşti, Appl. and Ind. Math. Series, No. 5, 2001.
[21] Suciu, N. and Vamos, C., Effective diffusion in heterogeneous media, Internal Report ICG-IV.00303, Research Center Jülich, 2003.
[22] Suciu, N., Vamoş, C., Vanderborght, J., H. Hardelauf and Vereecken, H., Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi: 10.1029/2005WRR004546, 2006.
[23] Suciu, N., Vamoş, C. and J. Eberhard, Evaluation of the first-order approximations for transport in heterogeneous media, Water Resour. Res., 42, W11504, doi: 10.1029/2005WR004714, 2006.
[24] Suciu, N. and Vamoş, C., 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, 2007.
[25] Suciu N., Vamoş, C., Sabelfeld, K. and Andronache, C., Memory effects and ergodicity for diffusion in spatially correlated velocity fields, Proc. Appl. Math. Mech., 7, 2010015-2010016, doi: 10.1002/pamm.20070057, 2007.
[26] Suciu N., Vamos, C., Vereecken, H., Sabelfeld, K. and Knabner, P., Dependence on initial conditions, memory effects, and ergodicity of transport in heterogeneous media, Preprint No. 324, Institute of Applied Mathematics, Friedrich-Alexander University Erlangen-Nuremberg (available online at http://www.am.uni-erlangen.de/de/preprints2000.html), 2008.
[27] Suciu, N., Vamos, C., Vereecken, H., Sabelfeld, K. and Knabner, P., Memory effects induced by dependence on initial conditions and ergodicity of transport in heterogeneous media, Water Resour. Res., 44, W08501, doi: 10.1029/2007WR006740, 2008.
[28] Vamoş, C., Suciu, N., Vereecken, H., Vanderborht, J. and Nitzsche, O., Path decomposition of discrete effective diffusion coefficient, Internal Report ICG-IV.00501, Research Center Jülich, 2001.
[29] Vamoş, C., Suciu, N. and Vereecken, H., Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comp. Phys., 186(2), pp. 527-544, doi: 10.1016/S0021-9991(03)00073-1, 2003.
[30] van Kampen, N. G., Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1981.
[31] Yaglom, A. M., Correlation Theory of Stationary and Related Random Functions, Volume I: Basic Results, Springer-Verlag, New York, 1987.
[32] Zirbel, C. L., Lagrangian observations of homogeneous random environments, Adv. Appl. Prob., 33, pp. 810-835, 2001.
[33] Zwanzig, R., Memory effects and irreversible thermodynamics, Phys. Rev., 124(4), pp. 983-992, 1961.

Received by the editors: May 8, 2008.

$^{1}$ Friedrich-Alexander University of Erlangen-Nuremberg, Institute of Applied Mathematics, Martensstrasse 3, 91058 Erlangen, Germany, e-mail: suciu@am.uni-erlangen.de, web: http://www.am.uni-erlangen.de.
$^{2}$ "T. Popoviciu" Institute of Numerical Analysis, Romanian Academy, 400320 Cluj-Napoca, P. O. Box 68-1, Romania, e-mail: cvamos@ictp.acad.ro, web: http://www.ictp.acad.ro.
$^{3}$ Research Center Jülich, ICG-IV: Agrosphere Institute, 52425 Jülich, Germany, e-mail: h.vereecken@fz-juelich.de, web: http://www.fz-juelich.de/icg/icg-4.
$^{4}$ Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstr., 39, 10117 Berlin; Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch of Russian Academy of Sciences, Prospect Lavrentjeva 6, 630090 Novosibirsk, Russia, e-mail: sabelfel@wias-berlin.de, karl@osmf.sscc.ru, web: http://www.wias-berlin.de, http://www.sscc.ru.
$^{5}$ Friedrich-Alexander University of Erlangen-Nuremberg, Institute of Applied Mathematics, Martensstrasse 3, 91058 Erlangen, Germany, email: knabner@am.uni-erlangen.de, web: http://www.am.uni-erlangen.de.

2008

Influence of deterministic trend on the estimated parameters of GARCH(1,1) model

August 10, 2018

Existence results for systems of nonlinear evolution equations

April 14, 2022

Nonresonance and existence for systems of nonlinear operator equations

April 14, 2022

Itô equation model for dispersion of solutes in heterogeneous media

Abstract

Authors

Keywords

References

Paper coordinates

PDF

About this paper

Journal

Publisher Name

Paper on journal website

Print ISSN

Online ISSN

MR

ZBL

Google Scholar

References

Paper in html format

ITÔ EQUATION MODEL FOR DISPERSION OF SOLUTES IN HETEROGENEOUS MEDIA

Abstract

1. INTRODUCTION

2. ITÔ FORMALISM

3. DISPERSION TERMS AND DISPERSION COEFFICIENTS

4. FIRST-ORDER APPROXIMATIONS

5. NUMERICAL SIMULATIONS

6. DISCUSSION AND CONCLUSIONS

REFERENCES

Related Posts