Let $\left\{X_{i,t}=X_{i}(t),i=1,2,3,t\geq 0\right\}$ be the advection-dispersion process with a space variable drift $\mathbf{V}(\mathbf{x})$, sample of a statistically homogeneous random

velocity field, and a local diffusion coefficient $D$, assumed constant, described by the Itô equation

where $x_{i,0}$ is a deterministic initial position and $W_{i}$ are the components of a Wiener process of mean zero and variance $\left\langle W_{i,t}^{2}\right\rangle=2Dt$.

The solution $g$ of the Fokker-Planck equation

which corresponds to Green’s function in Eulerian approaches to dispersion in random fields, is the density of the transition probability $P\left(d\mathbf{x},t,\mathbf{x}_{0}\right)=g\left(\mathbf{x},t\mid\mathbf{x}_{0}\right)d\mathbf{x}$ of the Itô process (1). Equation (2) is the starting point in Eulerian perturbation approaches $[1,7,8]$.

The following three centered processes are useful to describe the effective and ensemble dispersion and the fluctuations of the center of mass: $X_{i,t}^{eff}=X_{i,t}-\left\langle X_{i,t}\right\rangle$, the process centered about the mean $\left\langle X_{i,t}\right\rangle$ of $X_{i,t}$ in single realizations of the random velocity field, $\left\langle X_{i,t}^{eff}\right\rangle=0$, the process $X_{i,t}^{ens}=X_{i,t}-\overline{\left\langle X_{i,t}\right\rangle}$ centered about the average $\overline{\left\langle X_{i,t}\right\rangle}$ of $\left\langle X_{i,t}\right\rangle$ over the ensemble of realizations, $\overline{\left\langle X_{i,t}^{ens}\right\rangle}=0$, and the fluctuation of the center of mass $X_{i,t}^{cm}=\left\langle X_{i,t}\right\rangle-\overline{\left\langle X_{i,t}\right\rangle}$, with $\overline{X_{i,t}^{cm}}=0$. The half derivative of the variance of these tree processes defines, respectively, the effective and ensemble dispersion coefficients [1], and the center of mass coefficient [24],

related by $D_{ii}^{cm}(t)=D^{ens}(t)-D^{eff}(t)$.
According to (1), the process $X_{i,t}^{\text{ens }}$ starting from $x_{i,0}=0$ verifies

where $u_{i}=V_{i}-\overline{\left\langle V_{i}\right\rangle}$. The variance of the process (4) is obtained as follows,

where we used the statistical homogeneity of the velocity field, $\overline{\left\langle u_{i}\right\rangle}=0$, and the independence of the Wiener process on velocity statistics, which, together, cancel the cross-correlation term. The process $X_{i,t}^{cm}$ verifies (4) with $u_{i}=\left\langle V_{i}\right\rangle-\overline{\left\langle V_{i}\right\rangle}$ and the second term set to zero and has the variance

Further, the definitions (3) yield

The first iteration of (4) about the un-perturbed problem consisting of diffusion with constant coefficient $D$ in the mean flow field $U=\overline{\left\langle V_{i}\right\rangle}$, with trajectory $X_{i,t}^{(0)}=\delta_{j,1}Ut+W_{i,t}$, leads to replacing $\mathbf{X}_{t}$ by $\mathbf{X}_{t}^{(o)}$ in (5) and (6). Then, (5) gives the following first-order approximation of the ensemble dispersion coefficient

where $p\left(\mathbf{x},t;\mathbf{x}^{\prime},t^{\prime}\right)=\left\langle\delta\left(\mathbf{x}-\mathbf{X}_{t}^{(0)}\right)\delta\left(\mathbf{x}^{\prime}-\mathbf{X}_{t^{\prime}}^{(0)}\right)\right\rangle$ is the Gaussian joint probability density of the process $X_{i,t}^{(0)}$ expressed with the aid of the Dirac $\delta$-function [23]. Similarly, with $p(\mathbf{x},t)=\left\langle\delta\left(\mathbf{x}-\mathbf{X}_{t}^{(0)}\right)\right\rangle$, (6) gives the first-order approximation of the center of mass coefficient

Finally, $D_{ii}^{eff}(t)$ is obtained from (7). Together, (7-9) give the (stochastic) Lagrangian formulation of the first-order approximation equivalent to that obtained in the Eulerian perturbation approach [1,7,8].

Remark 1. As follows from (8-9), for a superposition $u_{i}(\mathbf{x})=u_{i}^{1}(\mathbf{x})+\cdots+u_{i}^{n}(\mathbf{x})$ of statistically independent velocity fields the dispersion coefficients estimated to the first-order are given by the sum of terms determined by the corresponding correlations $\overline{u_{i}^{1}(\mathbf{x})u_{i}^{1}\left(\mathbf{x}^{\prime}\right)},\cdots,\overline{u_{i}^{n}(\mathbf{x})u_{i}^{n}\left(\mathbf{x}^{\prime}\right)}$.

In the following, by $X_{t}$ we denote one of the processes $X_{i,t},X_{i,t}^{eff},X_{i,t}^{ens}$, or $X_{i,t}^{cm}$. Further, we consider the uniform time partition $0<\Delta<\cdots<k\Delta<\cdots<(K-1)\Delta<K\Delta$ and for $k>0$ we define the increments $X_{k,k-1}=X_{t_{k}}-X_{t_{k-1}}$ and $s_{k,k-1}=\left(X_{k,k-1}\right)^{2}$. For a sum of two increments $X_{k,0}=X_{k-1,0}+X_{k,k-1}$ we have

where $m_{k,k-1,0}=2X_{k-1,0}X_{k,k-1}$. By iterating the binomial formula (10) one obtains:

Summing up these relations we get

which, using the relation

can be further expressed as

For $X_{t}=X_{t}^{\text{ens }}$, the expectation $\Sigma_{k,0}=\overline{\left\langle s_{k,0}\right\rangle}$ of (10) verifies

where $M_{k,k-1,0}=\overline{\left\langle m_{k,k-1,0}\right\rangle}$. For continuous time $t=k\Delta$, (13) becomes

According to the definition of the ensemble coefficient in (3), the half time derive of (14) yields

The last term of (14), $M_{t,t-\Delta,0}$ is a memory term quantifying the departure from normal-diffusion behavior of the process $X_{t}^{ens}$ [25,22] and its half derivative $MD_{t,t-\Delta,0}^{ens}=(1/2)dM_{t,t-\Delta,0}/dt$ is a memory dispersion coefficient. Relations similar to (14-15) also hold for $S,D^{eff}$, and for $R,D^{cm}$. The expectation of the second term of Eq. (11),

is a cumulative memory term which, according to Eq. (12), contains a hierarchy of correlations between increments along paths of the transport process. Estimations of such cumulative memory terms by simulations of transport in single realizations of the random velocity field are presented in Ref. [29].

According to (15) and (5), the memory dispersion coefficient corresponding to the increment of the process $X_{t}$ over the time interval $\tau$ has the following Lagrangian expression

With (8), the first order approximation of (17) is given by

Memory coefficients (18) as well as their time integrals (i.e. memory terms) were computed in [25] for short range correlations of the velocity field as well as for diffusion in perfectly stratified fields (model of Matheron and de Marsily). In [22], the time behavior of the memory terms for fractional Gaussian noise has been also analyzed. Jeon and Metzler [16] used correlations of increments to quantify the memory of two special anomalous diffusion processes (fractional Brownian motion and processes governed by fractional Langevin equation). Their correlations are in fact memory terms related to the variance of the process by Eq. (13). For instance, in case of fractional Brownian motion with variance $\Sigma_{t_{2},t_{1}}=\left(t_{2}-t_{1}\right)^{2H}$ and Hurst coefficient $0<H<1$, for two successive time intervals $\Delta$, from (14) one obtains

which is just Eq. (4.4) of Jeon and Metzler [16].

Approximations of the flow and transport equations for small variance of the logarithm of the hydraulic conductivity $K$ of the aquifer system lead to explicit functional dependencies of the dispersion coefficients on the covariance $C_{Y}$ of the random field $Y=\ln K$ [14].

With the common convention for the Fourier transform

the first-order approximation in $\sigma_{Y}^{2}$ yields (e.g., [14, 21, 7,8 ])

where $p_{j}(\mathbf{k})=\left(\delta_{j,1}-k_{1}k_{j}/k^{2}\right)$ are projectors which ensure the incompressibility of the flow.

Remark 2. The Dirac function $\delta\left(\mathbf{k}+\mathbf{k}^{\prime}\right)$ which ensures the statistical homogeneity in (22) has to be changed into $\delta\left(\mathbf{k}-\mathbf{k}^{\prime}\right)$ to obtain a similar first-order relation for $\overline{\hat{u}}_{j}(\mathbf{k})\hat{u}_{j}^{*}\left(\mathbf{k}^{\prime}\right)$, where $\hat{u}_{j}^{*}$ denotes the complex conjugate [5, 6, 21]. Even though it leads to the same results, this relation render the computation more complicated.

Using the inverse Fourier transform (21) and (22), the Eulerian velocity correlation from (8) and (18) can be computed as follows

With (23), the integrand of the time integral in (8) and (18) becomes

The passage to the last line in (24) was achieved by using the expression

of the characteristic function of the Gaussian increment $\left(\mathbf{X}_{t}^{(0)}-\mathbf{X}_{t^{\prime}}^{(0)}\right)$ of the un-perturbed process $X_{j,t}^{(0)}=\delta_{j,1}Ut+W_{j,t},j=1\cdots d$ (see e.g. [18]).

Remark 3. Changing the sign convention of the Fourier transform,

changes the sign of the first term of the characteristic function without altering the sign of the second term, $\exp\left[-ik_{1}U\left(t-t^{\prime}\right)-k^{2}D\left(t-t^{\prime}\right)\right]$ (see detailed computation of characteristic function for Gaussian random variables of Pa poulis and Pillai [18], Example 5-28, p. 153]. However, changing the sign of the imaginary exponent is equivalent to changing $U$ into $-U$ (this also changes the sign of $p_{j}$ (see Gelhar and Axness [14], Eq. (28a)), but (24) depends only on $p_{j}^{2}$ ). Since second moments and dispersion coefficients do not depend on the sign of the mean flow, the expression (26) gives the same result irrespective of the sign of the imaginary exponent.

Inserting (24) into (8) one obtains the explicit expression of the advective contribution to the ensemble dispersion coefficient,

With the change of variable $t^{\prime\prime}=t-t^{\prime}$, (25) becomes

Remark 4.1. The time derivative of (26) is identical with the result of Dagan ([5], Eq. (3.14); [6], Eq. (17)), and Schwarze et al. ([21], Eq. (13)). Dagan [6] derived the result by using the changed sign in the homogeneity condition (Remark 2) and opposite Fourier transform convention (Remark 3).

Remark 4.2. The opposite Fourier transform convention and the same homogeneity condition as in (22) change $ik_{1}Ut^{\prime\prime}$ into $-ik_{1}Ut^{\prime\prime}$ and let $-k^{2}Dt^{\prime\prime}$ unchanged (see Remark 3). This is the case in the derivation of the variance
$\Sigma_{jj}$ by Fiori and Dagan [12] (Eq. (14)), which is exactly the integral of (25) from above, with $ik_{1}Ut^{\prime\prime}$ replaced by $-ik_{1}Ut^{\prime\prime}$. Dentz et al. [7] used the same conventions as Fiori and Dagan [12], (i.e. Eqs. (16) and (25) in [7]) and obtained the same result. The result given by their Eq. (40) with $\tilde{g}_{0}$ replaced by the characteristic function of the un-perturbed problem ([7], Eq. (25)), is just the relation (26) from above, with exponent $ik_{1}Ut^{\prime\prime}$ replaced by $-ik_{1}Ut^{\prime\prime}$.

The corresponding memory coefficient is similarly obtained by inserting (24) into (18),

Making the change of variable $t^{\prime\prime}=t-t^{\prime}$ in (27), we see that the ensemble memory coefficients can be computed, after a simple change of the integration limits in the time integral, with the explicit expression (26) of the advective contributions $\delta\left\{D_{ii}^{ens}\right\}$ to the ensemble coefficients:

Remark 5. The memory coefficient (27) can also be obtained within the Eulerian approach of Dentz et al. [7] if in Eq. (18) from [8] the initial concentration $\hat{\rho}$ is replaced by the perturbation solution for point source at time $\tau$ given in Eq. (26) from [7], i.e. $\hat{\rho}(k,\tau)=g(k,\tau)$. Then Eq. (18) from [8] gives $\hat{g}(k,t-\tau)$. Further, inserting $\hat{g}(k,t-\tau)$ in Eq. (14) from [8] and "expanding the logarithm consistently up to second order in the fluctuations" (as described at the beginning of Sect. 3.2 of [8]), should reproduce the expression (27) of the memory coefficient $MD^{\text{ens }}(t,\tau,0)$.

To derive the explicit form of the effective coefficients (7) we first have to compute the center of mass coefficients (9). Similarly to (24), we obtain

The passage to the last line in (29) is now achieved through the expressions

of the characteristic functions of the Gaussian un-perturbed process $X_{j,t}^{(0)}=\delta_{j,1}Ut+W_{j,t},j=1\cdots d[18]$, where we used the property pointed out in Remark 3. Inserting (29) into (9), we obtain the center of mass coefficient

With the change of variable $t^{\prime\prime}=t-t^{\prime}$ (which implies $t+t^{\prime}=2t-t^{\prime\prime}$ ) the expression (30) of the center of mass coefficients becomes

The center of mass coefficients (30) and the advective contribution (26) to the ensemble dispersion coefficients determine the advective contributions to the effective coefficients according to (7):

Remark 6. The ensemble contribution (26) and the center of mass coefficient (30) particularize to isotropic local diffusion coefficient the general relations $M^{-}$and $M^{+}$, respectively, derived by Dentz et al. ([7] Eq. (47)). We also note that (30) is the time derivative of the relation (15) of Fiori and Dagan [12], which, together with Remarks 4.1 and 4.2, proves the equivalence of Eulerian expressions for the ensemble and effective dispersion coefficients of Dentz et al. [7, 8] and the Lagragian expressions derived by Fiori and Dagan [12] and by Vanderborght [31] for the corresponding variances. We note here that the Eulerian approach has the important advantage of providing first-order approximations of the concentration fields (e.g., Eq. (26) in [7], on which the derivation of the dispersion coefficients is actually based).

The memory coefficient of the center of mass is readily obtained, similarly to (28), by a change of the integration limits in the expression (31) of the center of mass dispersion coefficients:

Finally, the memory coefficient of the effective dispersion is obtained from the relations (32), (31), and (28),

Remark 7. The memory coefficients (28), (33), and (34) can be computed with the explicit expressions of the advective contribution to the ensemble dispersion coefficients (25) and of the center of mass coefficients (31) by changing the lower limit of the time integrals from 0 to $t-\tau$.

The first-order results for log-K fields with finite correlation lengths can be further used to develop approaches for a superposition of scales [9, 20]. To proceed, let us consider the particular case of an isotropic power-law covariance of $Y=\ln K$,

where $L$ is some length scale associated with the spatial extension of the system. Power-law covariances can be constructed by integrating covariances with continuously increasing finite integral scales $\lambda$ ([15], p. 370, expression 3.478 (1.)) by the relation

where $\mu>0$ and $\beta>0$. For $\mu=z^{2}$ and $p=2$, (36) yields a superposition of Gaussian covariances: