# Convergence analysis of stochastic structure-preserving schemes for
computing effective diffusivity in random flows

Junlong Lyu<EMAIL_ADDRESS>Zhongjian Wang<EMAIL_ADDRESS>Jack Xin<EMAIL_ADDRESS>Zhiwen Zhang<EMAIL_ADDRESS>Department of
Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China.
Department of Mathematics, University of California at Irvine, Irvine, CA
92697, USA.

###### Abstract

In this paper, we develop efficient stochastic structure-preserving schemes to
compute the effective diffusivity for particles moving in random flows. We
first introduce the motion of a passive tracer particle in random flows using
the Lagrangian formulation, which is modeled by stochastic differential
equations (SDEs). Then, we propose stochastic structure-preserving schemes to
solve the SDEs and provide rigorous convergence analysis for the numerical
schemes in computing effective diffusivity. The convergence analysis follows a
probabilistic approach, which interprets the solution process generated by our
numerical schemes as a Markov process. By exploring the ergodicity of the
solution process, we obtain a convergence analysis of our method in computing
long-time solutions of the SDEs. Most importantly, our analysis result reveals
the equivalence of the definition of the effective diffusivity by solving
discrete-type and continuous-type (i.e. Eulerian) corrector problems, which is
fundamental and interesting. Finally, we present numerical results to
demonstrate the accuracy and efficiency of the proposed method and investigate
the convection-enhanced diffusion phenomenon in two- and three-dimensional
incompressible random flows.

AMS subject classification: 37M25, 60J60, 60H35, 65P10, 65M75, 76M50.

###### keywords:

Convection-enhanced diffusion; random flows; structure-preserving schemes;
corrector problem; ergodic theory; Markov process.

##  1 Introduction

Diffusion enhancement in fluid advection has been studied for nearly a century
since the pioneering work of Sir G. Taylor [32]. It is a fundamental problem
to characterize and quantify the large-scale effective diffusion in fluid
flows containing complex and turbulent streamlines, which is of great
theoretical and practical importance; see e.g. [10, 12, 4, 20, 22, 25] and
references therein. Its applications can be found in many physical and
engineering sciences, including atmosphere science, ocean science, chemical
engineering, and combustion.

In this paper, we study the diffusion enhancement phenomenon for particles
moving in random flows, which is described by the following passive tracer
model, i.e., a stochastic differential equation (SDE) with a random drift,

$\displaystyle d{\bf X}(t)={\bf b}(t,{\bf X}(t),\omega)dt+\sigma d{\bf
w}(t),\quad{\bf X}(0)=0,$ (1)

where ${\bf X}(t)\in R^{d}$ is the position of the particle, $\sigma>0$ is the
molecular diffusivity, and $\\{{\bf w}(t)\\}_{t\geq 0}$ is the standard
$d$-dimensional Brownian motion. Here the velocity field ${\bf b}(t,{\bf
x},\omega)$, i.e., the random drift is modeled by a random field in order to
mimic the energy spectra of the turbulent flow [19, 22]. Specifically, we
assume ${\bf b}(t,{\bf x},\omega)$ is a zero mean, jointly stationary, ergodic
vector random field over a certain probability space, where $\omega$ is an
element of the probability space describing all possible environments. The
randomness in ${\bf b}(t,{\bf x},\omega)$ is independent of the randomness in
the Brownian motion ${\bf w}(t)$. In addition, we assume that the realizations
of ${\bf b}(t,{\bf x},\omega)$ are almost surely divergence-free, i.e.
$\nabla_{{\bf x}}\cdot{\bf b}(t,{\bf x},\omega)=0$. To guarantee the existence
of the solution to (1), ${\bf b}(t,{\bf x},\omega)$ should be at least almost
surely locally Lipschitz in x. To design numerical schemes and carry our
convergence analysis, we assume ${\bf b}(t,{\bf x},\omega)$ has certain
regularity in the physical space; see Assumption 3.1. We emphasize that since
any statement, such as the effective diffusivity, involving statistical
properties of the solution ${\bf X}(t)$, requires only convergence in law.
Thus, the regularity assumption on the velocity field is natural and will
facilitate our algorithm design and convergence analysis in this paper.

We are interested in studying the long-time large-scale behavior of the
particles ${\bf X}(t)$ in (1). Namely, whether the motion of the particles
${\bf X}(t)$ has a long-time diffusive limit? More specifically, let
$X_{\epsilon}(t)\equiv\epsilon X(t/\epsilon^{2})$ denote the rescaled process
of (1). We want to find conditions under which $X_{\epsilon}(t)$ converges in
law, as $\epsilon\rightarrow 0$, to a new Brownian motion with a certain
covariance matrix $D^{E}\in R^{d\times d}$, where $D^{E}$ is called the
effective diffusivity matrix. This problem is referred to as the
homogenization of time-dependent flow problem.

Computing the effective diffusivity matrix $D^{E}$ (i.e., homogenization of
time-dependent flows) has been widely studied under various conditions on the
flows. For spatial-temporal periodic velocity fields and random velocity
fields with short-range correlations, one can apply the homogenization theory
[1, 14, 16, 26] to compute the effective diffusivity matrix $D^{E}$, where
$D^{E}$ can be expressed in terms of particle ensemble average (Lagrangian
framework) or an average of solutions to corrector problems (Eulerian
framework).

The dependence of $D^{E}$ on the velocity field of the problem is highly
nontrivial. For time-independent Taylor-Green flows, the authors of [27]
proposed a stochastic splitting method and calculated the effective
diffusivity in the limit of vanishing molecular diffusion. For time-dependent
chaotic flows, we proposed a Lagrangian-type numerical integrator to compute
the effective diffusivity using structure-preserving schemes [34]. In the
subsequent work [33], we provided a sharp and uniform-in-time error estimate
for the numerical integrator in computing the effective diffusivity. However,
we point out that the method and the convergence analysis obtained in [34, 33]
were designated for flows generated from separable and deterministic
Hamiltonian only.

For random flows with long-range correlations, the long-time large-scale
behavior of the particle motion is complicated and difficult to study in
general, since various forms of anomalous diffusion, such as super-diffusion
and sub-diffusion may exist. The interested reader is referred to the review
paper [22], where anomalous diffusion was obtained in exactly solvable models.
See also [8] for progress in understanding of the intermittency (i.e., the
occurrence of large fluctuations in the velocity field on the small scales)
for the passive scalar transport in a turbulent velocity field.

There are several theoretical works on homogenization of time-dependent random
flows. Such results include, among others, [4] where the authors proved the
existence of the effective diffusivity for a two-dimensional time-dependent
incompressible Gaussian velocity field. In [20, 18], the authors proved the
homogenization of convection-diffusion in a time-dependent, ergodic,
incompressible random flow. In [11, 9], the authors proved some necessary
conditions under which the long-time behavior for convection-diffusion in a
turbulent flow is diffusive. There are some recent works on studying the
effective diffusivity in random flows; see e.g. [2, 5, 3, 21, 29, 28]. Those
results show that the dependence of the effective diffusivity upon the
molecular diffusion $\sigma$ and the velocity field ${\bf b}$ in the random
flow is complicated and how to describe this dependence is very difficult in
general. Additionally, it is difficult to study the existence of residual
diffusivity for the passive tracer model (1). The residual diffusivity refers
to the non-zero effective diffusivity in the limit of zero molecular diffusion
$\sigma$.

This motivates us to develop efficient numerical schemes so that we can
compute the effective diffusivity of random flows. Notice that these random
flows are generated from non-separable Hamiltonian, which are much more
difficult than the problems studied in [34]. In this work, we first propose an
implicit structure-preserving scheme to solve the SDE (1), in order to deal
with the non-separable Hamiltonian. Second, we provide a sharp error estimate
for the numerical scheme in computing effective diffusivity. Our analysis is
based on a probabilistic approach. We interpret the solution process generated
by our numerical scheme as a discrete Markov process, where the transition
kernel can be constructed according to the numerical scheme in solving (1). By
exploring the ergodicity of the solution process, we obtain a sharp
convergence analysis for our method. Most importantly, our convergence
analysis reveals the equivalence of the definition of the effective
diffusivity by solving discrete-type and continuous-type (i.e. Eulerian)
corrector problems; see Theorem 4.3, which is fundamental and interesting.
Finally, we present numerical results to demonstrate the accuracy of the
proposed method in computing effective diffusivity for several incompressible
random flows in both two- and three-dimensional space.

To the best of our knowledge, this paper appears to be the first one in the
literature to develop Lagrangian numerical methods to compute effective
diffusivity in random flows through the connection with the Eulerian corrector
problem. The probabilistic approach in the convergence analysis takes into
account the ergodic nature of the solution process and leads to a sharp error
estimate. Notice that if one chooses the Gronwall inequality in the error
estimate, one cannot get rid of the exponential growth pre-factor in the error
term, which makes the estimate not sharp. Moreover, the stochastic structure-
preserving Lagrangian scheme enables us to investigate the convection-enhanced
diffusion phenomenon in random flows. Especially, we can numerically study the
dependence of the effective diffusivity in the regime of small molecular
diffusion $\sigma$ and the setting of the velocity field ${\bf b}$ in the
random flow.

The rest of the paper is organized as follows. In Section 2, we briefly review
some existing results for diffusion in random flows and introduce the
definition of effective diffusivity by solving a continuous-type corrector
problem. In Section 3, we propose our stochastic structure-preserving schemes
in computing effective diffusivity for the passive tracer model (1). In
Section 4, we provide the convergence analysis for the proposed method based
on a probabilistic approach. In addition, we show the equivalence of the
definition of effective diffusivity through the discrete-type and continuous-
type corrector problems. In Section 5, we present numerical results to
demonstrate the accuracy and efficiency of our method. Concluding remarks are
made in Section 6.

##  2 Preliminaries

To make this paper self-contained, we give a brief review of existing results
on convection-enhanced diffusion in random flows and the effective
diffusivity. Since these are standard results, we adopt the notations that
were used in [11, 9].

###  2.1 Some formulations and results for diffusion in random flows

We first define a function space that satisfies stationary and ergodic
property in $R^{d}$. Let $(\mathcal{X},\mathcal{H},P_{0})$ be a probability
space. Let $\tau_{\bf{x}}$, ${\bf x}\in R^{d}$ be an almost surely continuous,
jointly measurable group of measure preserving transformation on $\mathcal{X}$
with the following properties:

  1. 1.

$\tau_{\bf{0}}=Id_{\mathcal{X}}$ and $\tau_{\bf
x+y}=\tau_{\bf{x}}\tau_{\bf{y}}$, $\forall~{}{\bf x},{\bf y}\in R^{d}$.

  2. 2.

The mapping $(\chi,\bf{x})\mapsto\tau_{\bf{x}}\chi$ is jointly measurable.

  3. 3.

$P_{0}(\tau_{\bf{x}}(A))=P_{0}(A)$, for ${\bf x}\in R^{d},A\in\mathcal{H}$.

  4. 4.

$\lim_{\bf{x}\to
0}P_{0}\big{(}\chi:|f\circ\tau_{\bf{x}}(\chi)-f(\chi)|\geq\eta\big{)}=0,$
$\forall f\in L^{2}(\mathcal{X})$ and $\forall\eta>0$.

  5. 5.

If $P_{0}\big{(}A\Delta\tau_{\bf{x}}(A)\big{)}=0$, $\forall~{}{\bf x}\in
R^{d}$, then $A$ is a trivial event, i.e., $P_{0}(A)$ is either 0 or 1.

One can verify that $\tau_{\bf x}$ induces a strongly continuous group of
unitary mapping $U^{\bf x}$ on $L^{2}(\mathcal{X})$, which satisfies

$\displaystyle U^{\bf x}f(\chi)=f(\tau_{\bf x}(\chi)),\ f\in
L^{2}(\mathcal{X}),~{}{\bf x}\in R^{d}.$ (2)

In addition, it is easily to find that the group $U^{\bf x}$ has $d$
independent, skew-adjoint generators $D_{k}:\mathcal{D}_{k}\to
L^{2}(\mathcal{X})$ corresponding to directions $\bf{e}_{k}$, $k=1,\cdots,d$.

We introduce some function spaces that are useful in the analysis. Let
$C_{b}^{m}(\mathcal{X})$ be the space of functions $f$ in the intersection of
the domains of $D^{\alpha}$ with
$||D^{\alpha}f||_{L^{\infty}(\mathcal{X})}<+\infty$, where
$\alpha=(\alpha_{1},...,\alpha_{d})$ is a multi-index, each component
$\alpha_{i}$ is a nonnegative integer, $\sum_{i=1}^{d}\alpha_{i}\leq m$, and
the partial derivative operator $D^{\alpha}=D_{1}^{\alpha_{1}}\circ
D_{2}^{\alpha_{2}}\circ\cdots\circ D_{d}^{\alpha_{d}}$. It is well known that
$C_{b}^{\infty}(\mathcal{X})=\cap_{m\geq 1}C_{b}^{m}(\mathcal{X})$ is dense in
$L^{p}(\mathcal{X}),1\leq p<+\infty$; see [6]. Let
$L^{2}_{0}(\mathcal{X})=\\{f\in L^{2}(\mathcal{X})|{\bf E}_{0}f=0\\}$, where
${\bf E}_{0}$ is the expectation associated with the probability measure
$P_{0}$.

Next, we incorporate the time variable and study the Markov property.
Following setting is standard for a general Markov process.

Let $\Omega$ be the space of $\mathcal{X}$-valued continuous function
$C([0,\infty);\mathcal{X})$ and let $\ell$ be its Borel $\sigma-$algebra. Let
$P^{t},\ t\geq 0$, be a strongly continuous Markov semigroup on
$L^{2}(\mathcal{X})$, which satisfies the following properties.

  1. 1.

$P^{t}\bf{1}=\bf{1}$ and $P^{t}f\geq 0$, if $f\geq 0$.

  2. 2.

$\int P^{t}fdP_{0}=\int fdP_{0}$, for all $f\in L^{2}(\mathcal{X}),\ t\geq 0$.

  3. 3.

${\bf{\bf E}}_{\chi}[f(\theta_{t+h}(\omega))|\ell_{\leq
t}]=P^{h}F(\omega(t))$, where $F(\chi):={\bf E}_{\chi}f$, for any $f\in
L^{1}(\Omega)$, $t$, $h\geq 0$, $\chi\in\mathcal{X}$.

In the property P3, ${\bf E}_{\chi}$ is the expectation associated with the
probability measures $P_{\chi}$, which can be considered as the conditional
probability for all events in $\ell$ condition on that initial point lies on
$\chi$. $\ell_{\leq t}$ are the $\sigma$-algebras generated by events
measurable up to time $t$, and $\theta_{t}(\omega)(\cdot):=\omega(\cdot+t),\
t\geq 0$ is the standard shift operator on the path space $(\Omega,\ell)$.

Moreover, we can define a measure $P$ on the path space $(\Omega,\ell)$
through

$\displaystyle P(B)=\int P_{\chi}(B)P_{0}(d\chi),\ B\in\ell$ (3)

and define ${\bf E}$ to be the corresponding expectation operator with respect
to the measure $P$. As a direct consequence of (T3) and (P2), we know that $P$
is stationary:

######  Proposition 2.1.

$P$ is invariant under the action of $\theta_{t}$ and $\tau_{\bf x}$ for any
$(t,{\bf x})\in R^{+}\times R^{d}$.

Let $L:\mathcal{D}(L)\to L^{2}(\mathcal{X})$ be the generator of the semigroup
$P^{t}$. To establish the central limit theorem for the Markov process
associated with $P^{t}$, we assume the generator $L$ satisfies the following
time relaxation property, also known as the spectral gap condition,

$\displaystyle-(Lf,f)_{L^{2}(\mathcal{X})}\geq
c_{1}||f||^{2}_{L^{2}(\mathcal{X})},~{}\text{where}~{}c_{1}>0.$ (4)

The time relaxation property (4) is equivalent to the exponential decay
property

$\displaystyle||P^{t}f||_{L^{2}(\mathcal{X})}\leq\exp(-2c_{1}t)||f||_{L^{2}(\mathcal{X})},~{}f\in
L^{2}_{0}(\mathcal{X}).$ (5)

In addition, time relaxation property (4) is equivalent to $\rho$-mixing of
the process $X(t)$, $t\geq 0$. Specifically, let
$\rho(h)=\sup\\{Cor(Y_{1},Y_{2}):Y_{1}~{}\text{is}~{}\ell_{\geq
t+h}~{}\text{measurable},Y_{2}~{}\text{is}~{}\ell_{\leq
t}~{}\text{measurable}\\}$, where $Cor(Y_{1},Y_{2})$ is the correlation
function. Then, (4) or (5) implies that $\lim_{h\rightarrow\infty}\rho(h)=0$;
see [30, 7]. The time relaxation property (4) (or the exponential decay
property (5)) plays an important role in proving the existing of the effective
diffusivity. We will numerically investigate this property in Section 5.

###  2.2 The continuous-type corrector problem and effective diffusivity

Equipped with the necessary properties and notations, we are ready to study
the effective diffusivity of the random flows associated with the passive
tracer model (1). First we assume that the random flow ${\bf
b}=(b_{1},...,b_{d})\in(L^{2}(\mathcal{X}))^{d}$ is jointly continuous in
$(t,\bf x)$, locally Lipschitzian in $\bf x$, with finite second moments, and
is divergence-free.

For each fixed realization $\omega$ of the environment, we consider the
stochastic process generated by the following SDE,

$\left\\{\begin{aligned} &d{\bf X}_{t}^{\omega}={\bf b}(t,{\bf
X}^{\omega}_{t},\omega)dt+\sigma d{\bf w}_{t},\\\ &{\bf
X}_{0}^{\omega}=0,\end{aligned}\right.$ (6)

where ${\bf X}^{\omega}_{t}\in R^{d}$ is the position of the particle, the
superscript in ${\bf X}_{t}^{\omega}$ means that it depends on the realization
of the environment $\omega$, ${\bf w}_{t}$ is a standard Brownian motion
starting at the origin. Its corresponding probability space is denoted by
$(\Sigma,\mathcal{B},Q)$ and the associated expectation operator is denoted by
${\bf M}$. The SDE (6) is well-defined [9]. Moreover, the random flow in (6)
means ${\bf b}(t,{\bf x},\omega)={\bf b}(\tau_{{\bf x}}\omega(t))$. Viewed
from a particle at any instant of time $t$, we can define an environment
process $\eta:[0,\infty)\times\Omega\to\mathcal{X}$ as

$\left\\{\begin{aligned} &\eta(t)=\tau_{{\bf X}_{t}^{\omega}}\omega(t),\\\
&\eta(0)=\omega(0).\end{aligned}\right.$ (7)

In addition, environment process generates a semigroup of transformation

$S^{t}f(\chi)={\bf ME}_{\chi}f(\eta(t)),~{}t\geq 0,~{}f\in
L^{\infty}(\mathcal{X}),$ (8)

where $\eta(t)$ is defined by (7). And $S^{t}$ satisfies the following
properties,

######  Proposition 2.2.

[[9], Prop. 3]

  1. 1.

$S^{t},t\geq 0$ is a strongly continuous, Markov semigroup of contraction on
$L^{2}(\mathcal{X})$.

  2. 2.

$S^{t},t\geq 0$ is measure-preserving, that is,

$\int S^{t}fdP_{0}=\int fdP_{0},~{}t\geq 0,~{}f\in L^{2}(\mathcal{X}).$ (9)

Let $D_{1}=\mathcal{D}(L)\cap C_{b}^{2}(\mathcal{X})$ and $\mathcal{L}$ denote
the generator of the semigroup $S^{t}$, $t\geq 0$, i.e.,

$\mathcal{L}f=Lf+\frac{\sigma^{2}}{2}\Delta f+{\bf b}\cdot\nabla f,$ (10)

where $L$ is the generator of the semigroup $P^{t}$. One can easily verify the
following properties.

######  Proposition 2.3.

[[9],Prop 4]

  1. 1.

$D_{1}$ is dense in $L^{2}(\mathcal{X})$ and is invariant under the semigroup
$P^{t}$, $t\geq 0$, i.e., $P^{t}(D_{1})\subseteq D_{1}$ for all $t\geq 0$.

  2. 2.

Assume that the random flow ${\bf b}$ is bounded. Then, $D_{1}$ is invariant
under the semigroup $S^{t}$, $t\geq 0$, i.e., $S^{t}(D_{1})\subseteq D_{1}$
for all $t\geq 0$.

######  Lemma 2.4.

From the spectral gap condition (4), we obtain that for any $f\in
L_{0}^{2}(\mathcal{X})$

$\displaystyle||S^{t}f||_{L^{2}(\mathcal{X})}\leq\exp(-2c_{1}t)||f||_{L^{2}(\mathcal{X})},~{}\text{where}~{}c_{1}>0.$
(11)

###### Proof.

We first assume $\bf b$ is bounded and $f\in
D_{1}\subseteq\mathcal{D}(\mathcal{L})$. Using the spectral gap condition and
$\bf b$ is divergence-free, we have

$(-\mathcal{L}f,f)_{L^{2}_{0}(\mathcal{X})}\geq(-Lf,f)_{L^{2}_{0}(\mathcal{X})}\geq
c_{1}||f||^{2}_{L^{2}_{0}(\mathcal{X})}$ (12)

for all $f\in D_{1}\cap L_{0}^{2}(\mathcal{X})$. By Proposition 2.3,
$S^{t}f\in D_{1},t\geq 0$ for any $f\in D_{1}$. Consequently,

$\frac{d}{dt}||S^{t}f||_{L^{2}(\mathcal{X})}^{2}=2(\mathcal{L}S^{t}f,S^{t}f)_{L^{2}(\mathcal{X})}\leq-2c_{1}||S^{t}f||_{L^{2}(\mathcal{X})}^{2},$
(13)

thus

$||S^{t}f||_{L^{2}(\mathcal{X})}^{2}\leq\exp(-2c_{1}t)||f||_{L^{2}(\mathcal{X})}^{2},\quad\forall
t\geq 0,$ (14)

for all $f\in D_{1}\cap L_{0}^{2}(\mathcal{X})$. Then, the statement in (11)
is extended to $L_{0}^{2}(\mathcal{X})$ by using an approximation argument.
Finally, the boundedness of the random flow ${\bf b}$ is removed by using
another approximation argument. ∎

Given the semigroup of transformation $S^{t}$ in (8) and its associated
properties (see Proposition 2.2), we can define

${\bm{\psi}}=\int_{0}^{\infty}S^{t}{\bf b}dt$ (15)

which satisfies the following continuous-type corrector problem

$\displaystyle\mathcal{L}{\bm{\psi}}=-{\bf b}$ (16)

where $\mathcal{L}$ is the generator of $S^{t}$ defined in (10). By solving
the corrector problem (16), we are able to define the effective diffusivity.
This can be summarized into the following result.

######  Proposition 2.5.

Let ${\bf X}(t)$ be the solution to (1) and ${\bf
X}_{\epsilon}(t)\equiv\epsilon{\bf X}(t/\epsilon^{2})$. For any unit vector
${\bf v}\in R^{d}$, let $\psi_{\bf v}=\bf\bm{\psi}\cdot v$ denote the
projection of the vector solution ${\bm{\psi}}$ along the direction ${\bf v}$,
where ${\bm{\psi}}$ is the solution to corrector problem (16). Then, the law
of the process ${\bf X}_{\epsilon}(t)\cdot{\bf v}$ converges weakly in
$C[0,+\infty)$ to a Brownian motion with diffusion coefficient given by

$\displaystyle{\bf v}^{T}D^{E}{\bf
v}=\frac{\sigma^{2}}{2}+(-\mathcal{L}\psi_{\bf v},\psi_{\bf
v})_{L^{2}(\mathcal{X})},$ (17)

where $D^{E}$ is the effective diffusivity associated with the passive tracer
model (1).

The proof of Prop. 2.5 relies on an approximation of the additive functional
of an ergodic Markov process by a martingale and applying the central limit
theorem to continuous-time Markov process, which is very useful in studying
the long-time behavior of random dynamics; see Lemma 1 of [9] or Theorem of
[4]. We shall prove in Theorem 4.3 that the numerical solutions obtained by
our Lagrangian numerical scheme recover the definition of the effective
diffusivity in (17).

##  3 Stochastic structure-preserving schemes and related properties

###  3.1 Derivation of numerical schemes

In this section, we construct numerical schemes for the passive tracer model
(6), which is based on an operator splitting method [31]. For each fixed
realization $\omega$ of the environment, we first split the original problem
(6) into two sub-problems.

$\displaystyle d{\bf X}_{t}^{\omega}$ $\displaystyle={\bf b}(t,{\bf
X}^{\omega}_{t},\omega)dt,$ (18) $\displaystyle d{\bf X}_{t}^{\omega}$
$\displaystyle=\sigma d{\bf w}_{t},$ (19)

where we assume ${\bf w}_{t}$ in (19) is the same process as in (6). Let ${\bf
X}^{\omega}_{n}$ denote the numerical solution of ${\bf X}_{t}^{\omega}$ at
time $t=t_{n}$, $n=0,1,2,...$. From time $t=t_{n}$ to time $t=t_{n+1}$, where
$t_{n+1}=t_{n}+\Delta t$, $t_{0}=0$, assuming the solution ${\bf
X}^{\omega}_{n}$ is given, we now discuss how to discretize the above two sub-
problems (18)-(19), separately.

In the sub-problem (18), the velocity ${\bf b}(t,{\bf x},\omega)$ is almost
surely divergence-free and has certain regularity in the physical space. Thus,
we apply a volume-preserving scheme to discretize (18). Let ${\bf\Phi}_{\Delta
t}$ denote the numerical integrator associated with the volume-preserving
scheme during $\Delta t$ time and let ${\bf D_{{\bf x}}\Phi}_{\Delta t}$
denote the corresponding Jacobian matrix. The volume-preserving property
requests $\det({\bf D_{{\bf x}}\Phi}_{\Delta t})=1$. We obtain the numerical
integrator for the sub-problem (18) as follows,

$\displaystyle{\bf X}^{\omega}_{n+1}={\bf\Phi}^{\omega(n\Delta t)}_{\Delta
t}\big{(}{\bf X}^{\omega}_{n}\big{)},$ (20)

where the superscript in ${\bf\Phi}^{\omega(n\Delta t)}_{\Delta t}$ means that
the numerical integrator implicitly depends on the realization of ${\bf b}$ at
different computational times. Suppose $\bf b$ has bounded first derivatives
with respect to $\bf x$ for almost all $\omega$, it is easy to verify that the
volume-preserving integrator ${\bf\Phi}^{\omega(n\Delta t)}_{\Delta t}$ also
has bounded first derivatives for $\Delta t$ small enough. Thus,
${\bf\Phi}^{\omega(n\Delta t)}_{\Delta t}$ is well defined.

In addition, we assume that the numerical scheme only relies on the
information of ${\bf X}$ and ${\bf b}$ at the beginning of each computational
time, in order to make sure the solution process generated by our method is a
Markov process. For instance, to compute ${\bf X}^{\omega}_{n+1}$ the
numerical scheme only relies on the information of ${\bf X}$ and ${\bf b}$ at
$t=t_{n}$.

We illustrate this idea by constructing a volume-preserving scheme for a two-
dimensional problem. Let ${\bf
X}^{\omega}_{n}=(X^{\omega}_{n,1},X^{\omega}_{n,2})^{T}$ denote the numerical
solution at time $t=t_{n}$ and the velocity ${\bf b}(t,{\bf
x},\omega)=(b_{1}(t,{\bf x},\omega),b_{2}(t,{\bf x},\omega))^{T}$. Then, we
use the following numerical scheme to discretize (18)

$(X^{\omega}_{{n+1},1},X^{\omega}_{{n+1},2})^{T}=(X^{\omega}_{n,1},X^{\omega}_{n,2})^{T}+\Delta
t{\bf
b}\big{(}t_{n},(\frac{X^{\omega}_{n,1}+X^{\omega}_{{n+1},1}}{2},\frac{X^{\omega}_{n,2}+X^{\omega}_{{n+1},2}}{2})^{T},\omega\big{)},$
(21)

where we evaluate the velocity ${\bf b}(t,{\bf x},\omega)$ at $t=t_{n}$ to
ensure the Markov property. By solving Eq.(21) to get
$(X^{\omega}_{{n+1},1},X^{\omega}_{{n+1},2})^{T}$, we implicitly define a
numerical integrator ${\bf\Phi}^{\omega(n\Delta t)}_{\Delta t}$; see Eq.(20).
Since ${\bf b}(t,{\bf x},\omega)$ is almost surely divergence-free, we can
easily verify that the scheme (21) is volume-preserving, i.e, $\det({\bf
D_{{\bf x}}\Phi}^{\omega(n\Delta t)}_{\Delta t})=1$. As we will demonstrate in
the proof of Theorem 3.4 that using a volume-preserving numerical scheme to
discretize (18) is essential.

For a $d$-dimensional sub-problem (18), we split the velocity field ${\bf
b}(t,{\bf x},\omega)$ into a summation of $d-1$ velocity fields, where each of
them will generate a two-dimensional problem and thus we can design the
volume-preserving scheme accordingly. By applying a splitting method [23], we
can construct volume-preserving schemes for the original $d$-dimensional sub-
problem (18). More details can be found in [13, 15].

Given the numerical integrator $\bm{\Phi}^{\omega(n\Delta t)}_{\Delta t}$, we
define the mapping

$\displaystyle{\bf B}_{\Delta t}^{\omega(n\Delta t)}({\bf
x})={\bf\Phi}^{\omega(n\Delta t)}_{\Delta t}({\bf x})-{\bf x}.$ (22)

One can easily verity that ${\bf B}_{\Delta t}^{\omega(n\Delta t)}({\bf
X}^{\omega}_{n})$ is an approximation of the increment for the exact solution
of the sub-problem (18) as follows,

$\displaystyle{\bf X}^{\omega}_{(n+1)\Delta t}-{\bf X}^{\omega}_{n\Delta
t}=\int_{n\Delta t}^{(n+1)\Delta t}{\bf b}(t,{\bf X}_{t}^{\omega},\omega)dt.$
(23)

The sub-problem (19) can be exactly solved by many numerical schemes for SDEs,
including the Euler-Maruyama scheme [17].

Finally, we apply the Lie-Trotter splitting method and get the stochastic
structure-preserving scheme as follows,

$\displaystyle{\bf X}^{\omega}_{n+1}={\bf X}^{\omega}_{n}+{\bf B}_{\Delta
t}^{\omega(n\Delta t)}({\bf X}^{\omega}_{n})+\sigma\bm{\xi}_{n},$ (24)

where $\bm{\xi}_{n}=(\xi_{1},...,\xi_{d})^{T}$ is a $d$-dimensional i.i.d.
mean-free Gaussian random vector with ${\bf
E}{\bm{\xi}}_{n}\otimes{\bm{\xi}}_{n}=\Delta t{\bf I}_{d}$. Here ${\bf I}_{d}$
is an identity matrix.

The volume-preserving schemes for the sub-problem (18) are implicit in
general. Compared with explicit schemes, however, they allow us to choose a
relatively larger time step to compute. In practice, we find that a few steps
of Newton iterations are good enough to maintain accurate results. Therefore,
the computational cost is controllable. To design adaptive time-stepping
method for the passive tracer model (6) is an interesting issue, which will be
studied in our future work.

In general, the second-order Strang splitting [31] is more frequently used in
developing numerical schemes. In fact, the only difference between the Strang
splitting method and the Lie-Trotter splitting method is that the first and
last steps are half of the time step $\Delta t$. For the SDEs, however, the
dominant source of error comes from the random subproblem (19). Thus, it is
not necessary to implement the Strang splitting scheme here.

###  3.2 Some properties of the numerical schemes

In this subsection, we shall prove some properties of the proposed stochastic
structure-preserving scheme. Especially, we shall show that some important
properties of the random flows are maintained after numerical discretization.
Before proceeding to the analysis, we first introduce some notations and
assumptions. To emphasize the properties in spatial-domain, for any $f\in
L^{1}(\mathcal{X})$, we use $f^{\chi}({\bf x})$ to represent $f(\tau_{\bf
x}\chi)$. Moreover, we denote ${\bf b}(t,{\bf x},\omega)={\bf b}(\tau_{\bf
x}\omega(t))$, where $\tau_{\bf x}\omega(t)\in\mathcal{X}$.

######  Assumption 3.1.

Suppose the velocity field has certain regularity with respect to spatial
variables, i.e., ${\bf b}\in(C_{b}^{m}(\mathcal{X}))^{d}$ for some $m\geq 1$,
and has first-order partial derivative bounded with respect to temporal
variable, i.e., $||{\bf D}_{t}{\bf b}||_{L^{\infty}(\mathcal{X})}\leq
c<\infty$.

######  Assumption 3.2.

${\bf B}_{\Delta t}^{\chi}({\bf x})$ defined in (22) is a stationary process
with respect to ${\bf x}$, i.e., we can write ${\bf B}_{\Delta t}^{\chi}({\bf
x})={\bf B}_{\Delta t}(\tau_{\bf x}\chi)$.

######  Assumption 3.3.

If $\Delta t$ is small enough, we have ${\bf B}_{\Delta
t}\in(C_{b}^{m}(\mathcal{X}))^{d}$ provided that ${\bf
b}\in(C_{b}^{m}(\mathcal{X}))^{d}$. In addition, $||{\bf B}_{\Delta
t}||_{C_{b}^{m}(\mathcal{X})}=K||{\bf b}||_{C_{b}^{m}(\mathcal{X})}\Delta t$,
where $K$ is a constant that does not depend on $\Delta t$.

Under Assumption 3.1, we compute the local truncation error of the numerical
scheme (21). Recall that the numerical solution ${\bf
X}^{\omega}_{n}=(X^{\omega}_{n,1},X^{\omega}_{n,2})^{T}$ at time $t=n\Delta
t$. We rewrite (21) into a compact form as follows,

${\bf X}^{\omega}_{n+1}={\bf X}^{\omega}_{n}+\Delta t{\bf
b}\big{(}t_{n},\frac{{\bf X}^{\omega}_{n}+{\bf
X}^{\omega}_{n+1}}{2},\omega\big{)}.$

Assume ${\bf X}^{\omega}_{n}$ is equal to the exact solution ${\bf
X}^{\omega}_{t}$ at time $t=n\Delta t$. Then, we can obtain the exact solution
of the sub-problem (18) at time $t=(n+1)\Delta t$ as

${\bf X}^{\omega}_{(n+1)\Delta t}={\bf X}^{\omega}_{n}+\int_{n\Delta
t}^{(n+1)\Delta t}{\bf b}(t,{\bf X}_{t}^{\omega},\omega)dt.$

Let ${\bf T}^{\omega}_{n+1}$ denote the local truncation error at time
$t=(n+1)\Delta t$. We have

$\displaystyle{\bf T}^{\omega}_{n+1}={\bf X}^{\omega}_{(n+1)\Delta t}-{\bf
X}^{\omega}_{n+1}=\int_{n\Delta t}^{(n+1)\Delta t}\big{(}{\bf b}(t,{\bf
X}_{t}^{\omega},\omega)-{\bf b}(t_{n},\frac{{\bf X}^{\omega}_{n}+{\bf
X}^{\omega}_{n+1}}{2},\omega)\big{)}dt.$ (25)

We know that ${\bf b}$ satisfies an inequality of the following form

$\displaystyle\big{|}\big{|}{\bf b}(t,{\bf X}_{t}^{\omega},\omega)-{\bf
b}(t_{n},\frac{{\bf X}^{\omega}_{n}+{\bf
X}^{\omega}_{n+1}}{2},\omega)\big{|}\big{|}$ $\displaystyle\leq$
$\displaystyle\big{|}\big{|}{\bf b}(t,{\bf X}_{t}^{\omega},\omega)-{\bf
b}(t,\frac{{\bf X}^{\omega}_{n}+{\bf
X}^{\omega}_{n+1}}{2},\omega)\big{|}\big{|}+\big{|}\big{|}{\bf b}(t,\frac{{\bf
X}^{\omega}_{n}+{\bf X}^{\omega}_{n+1}}{2},\omega)-{\bf b}(t_{n},\frac{{\bf
X}^{\omega}_{n}+{\bf X}^{\omega}_{n+1}}{2},\omega)\big{|}\big{|},$
$\displaystyle\leq$ $\displaystyle||{\bf D_{{\bf x}}}{\bf
b}||_{L^{\infty}(\mathcal{X})}\big{|}\big{|}{\bf X}_{t}^{\omega}-\frac{{\bf
X}^{\omega}_{n}+{\bf X}^{\omega}_{n+1}}{2}\big{|}\big{|}+||{\bf D}_{t}{\bf
b}||_{L^{\infty}(\mathcal{X})}\big{|}t-t_{n}\big{|},$ (26)

where $t_{n}\leq t\leq t_{n+1}$ and $||\cdot||$ denotes the Euclidean norm of
a vector. If Assumption 3.1 holds true, we can easily obtain that the local
truncation error ${\bf T}^{\omega}_{n+1}=O(\Delta t)^{2}$, where the constant
in $O(\Delta t)^{2}$ depends on $||{\bf D_{{\bf x}}}{\bf
b}||_{L^{\infty}(\mathcal{X})}$ and $||{\bf D}_{t}{\bf
b}||_{L^{\infty}(\mathcal{X})}$.

We restrict ourselves to the convergence analysis based on Assumption 3.1 in
this paper. In fact, when ${\bf b}$ satisfies a Hölder-$\gamma$ continuous
assumption in time domain with $0<\gamma<1$, the local truncation error of
(21) becomes $O(\Delta t)^{1+\gamma}$. We can still prove the convergence
analysis of our method for computing effective diffusivity in such kind of
flows; see Remark 4.4.

As an analogy to the continuous-time case (7), we define the environment
process as viewed from the numerical solution ${\bf X}_{n}^{\omega}$ at
different time steps

$\left\\{\begin{aligned} &\eta_{n}=\tau_{{\bf X}_{n}^{\omega}}\omega(n\Delta
t),\\\ &\eta_{0}=\omega(0).\end{aligned}\right.$ (27)

The above environment process is defined on the space of trajectories
$(\tilde{\Omega},\ell)$, where $\tilde{\Omega}=C([0,\infty)\cap\Delta
t\mathbb{Z};\mathcal{X})$ is a subspace of $\Omega$ with time parameter lies
only on $\Delta t\mathbb{Z}$. The corresponding expectation operator is still
denoted by ${\bf E}_{\chi}$, which is the same as the one defined in property
P3 in Section 2.1. Under this process, we can write ${\bf B}_{\Delta
t}(\eta_{n})={\bf B}_{\Delta t}^{\omega(n\Delta t)}(X^{\omega}_{n})$. In
addition, we define

$\displaystyle\ {S_{n}f(\chi)={\bf ME}_{\chi}f(\eta_{n})},$ (28)

where ${\bf M}$ denotes the expectation respect to ${\bf w}_{t}$. We shall
prove that $S_{n}$ is a discrete-time Markov semi-group of contraction on
$L^{2}(\mathcal{X})$ and is measure-preserving with respect to $P_{0}$ defined
in Section 2.1. For clarity, we denote $\mathbb{E}$ the total expectation with
respect to all randomness, i.e., $\mathbb{E}={\bf ME}$, in the remaining part
of this paper.

######  Theorem 3.4.

$P_{0}$ is an invariant probability measure of $\eta_{n}$, i.e., $P_{0}$ is an
invariant measure of the Markov semigroup $\\{S_{n}\\}$.

###### Proof.

Let $p^{1}_{\chi}({\bf x,y})$ denote the transition probability density of the
solution process, which is defined by applying the numerical scheme (24) for
one time step. For simplicity of notation, let ${\bf x}$ be the current
solution and ${\bf y}$ be the solution obtained by applying the scheme (24)
with time step $\Delta t$. Notice that $\bm{\xi}_{n}$ in (24) is a mean-free
Gaussian random vector. We have

$\displaystyle p^{1}_{\chi}({\bf x,y})=\frac{1}{(2\pi\sigma^{2}\Delta
t)^{d/2}}\exp\left(-\frac{||{\bf y-x-B}^{\chi}_{\Delta t}({\bf
x})||^{2}}{2\sigma^{2}\Delta t}\right)=\frac{1}{(2\pi\sigma^{2}\Delta
t)^{d/2}}\exp\left(-\frac{||{\bf y-\Phi}^{\chi}_{\Delta t}({\bf
x})||^{2}}{2\sigma^{2}\Delta t}\right).$ (29)

Let us define $p_{0}({\bf x,y})=\frac{1}{(2\pi\sigma^{2}\Delta
t)^{d/2}}\exp\left(-\frac{||{\bf y-x}||^{2}}{2\sigma^{2}\Delta t}\right)$.
Then, we can verify that

$\displaystyle\int p^{1}_{\chi}({\bf x,y})d{\bf x}$ $\displaystyle=\int
p_{0}({\bf x+B}^{\chi}_{\Delta t}({\bf x),y})d{\bf x},$ $\displaystyle=\int
p_{0}({\bf z,y})\det({\bf D\Phi}^{\chi}_{\Delta t})^{-1}d{\bf z}=\int
p_{0}({\bf z,y})d{\bf z}=1,~{}a.e.~{}\chi,$ (30)

where we have used the fact that the numerical scheme (20) for sub-problem
(18) is volume-preserving, i.e., $\det({\bf D\Phi}^{\chi}_{\Delta t})=1$.
Thus, for all $f\in L^{2}(\mathcal{X})$, we have

$\displaystyle\int_{\mathcal{X}}S_{1}f(\chi)P_{0}(d\chi)$
$\displaystyle=\int_{\mathcal{X}}{\bf
E}_{\chi}f(\eta_{1})P_{0}(d\chi)=\int_{\mathcal{X}}P_{0}(d\chi)\int_{R^{d}}p^{1}_{\chi}({\bf
0,y}){\bf E}_{\chi}f(\tau_{\bf y}\omega(\Delta t))d{\bf y},$
$\displaystyle=\int_{\mathcal{X}}{\bf E}_{\chi}f(\omega(\Delta
t))P_{0}(d\chi)\int_{R^{d}}p^{1}_{\tau_{\bf-y}\chi}(\bf 0,y)dy,$
$\displaystyle=\int_{\mathcal{X}}{\bf E}_{\chi}f(\omega(\Delta
t))P_{0}(d\chi)\int_{R^{d}}p^{1}_{\chi}({\bf-y,0)dy},$
$\displaystyle=\int_{\mathcal{X}}{\bf E}_{\chi}f(\omega(\Delta
t))P_{0}(d\chi),$ (31)

where we have used the facts that $p^{1}_{\tau_{\bf x}\chi}({\bf
y,z})=p^{1}_{\chi}({\bf x+y,x+z})$ and
$\int_{R^{d}}p^{1}_{\chi}({\bf-y,0)dy}=1$. The first equality is easy to
verify, since

$\displaystyle p^{1}_{\tau_{\bf x}\chi}({\bf y},{\bf z})$
$\displaystyle=p_{0}({\bf y}+{\bf B}_{\Delta t}^{\tau_{\bf x}\chi}({\bf
y}),{\bf z})=p_{0}({\bf y}+{\bf B}_{\Delta t}^{\chi}({\bf x}+{\bf y}),{\bf
z}),$ $\displaystyle=p_{0}({\bf x}+{\bf y}+{\bf B}_{\Delta t}^{\chi}({\bf
x}+{\bf y}),{\bf x}+{\bf z})=p^{1}_{\chi}({\bf x}+{\bf y},{\bf x}+{\bf z}).$

Thus, we obtain from (3.2) that ${\bf E}S_{1}f={\bf E}P^{\Delta t}f={\bf E}f$,
where $P^{\Delta t}$ is measure-preserving by property (P2) in Section 2.1.
Similar argument shows that ${\bf E}S_{n}f={\bf E}S_{n-1}f$ for all $n$. We
prove that $S_{n}$ is measure-preserving. ∎

######  Remark 3.1.

Theorem 3.4 plays an important role in the remaining part of our convergence
analysis. Throughout the proof, one can see that using a volume-preserving
numerical scheme for solving sub-problem (18) is essential.

######  Remark 3.2.

In the proof of Theorem 3.4, the probability measures $p^{1}_{\chi}({\bf
x,y})d{\bf y}$ and $p_{0}({\bf x,y})d{\bf y}$ are associated with the Brownian
motion in the passive tracer model. While $P_{0}(d\chi)$ is the probability
measure associated with the randomness in the velocity field and initial data.
In the remaining part of this paper, we shall keep the same notations.

The following lemma will be very useful in our analysis.

######  Lemma 3.5.

For any ${\bf y}\in R^{d}$ and $f\in L^{2}(\mathcal{X})$, we have that

$\displaystyle\mathbb{E}f(\tau_{\bf
y}\eta_{n})=\mathbb{E}f(\eta_{n-1})=\mathbb{E}f.$ (32)

Moreover,

$\displaystyle\mathbb{E}f(\eta_{n+1})=\mathbb{E}f\Big{(}\tau_{{\bf
X}^{\omega}_{n}+{\bf B}_{\Delta t}(\eta_{n})}\omega\big{(}(n+1)\Delta
t\big{)}\Big{)}=\mathbb{E}f.$ (33)

###### Proof.

We prove the above equations through direct calculations. For the equation
(32), we have

$\displaystyle\mathbb{E}f(\tau_{\bf y}\eta_{n})={\bf
EME}_{\eta_{n-1}}f(\tau_{\bf
y}\tilde{\eta}_{1})=\int_{\mathcal{X}}P_{0}(d\chi){\bf M}{\bf
E}_{\chi}\big{[}\int_{R^{d}}p^{1}_{\eta_{n-1}}({\bf 0,z}){\bf
E}_{\eta_{n-1}}f\big{(}\tau_{\bf y+z}\omega(\Delta t)\big{)}d{\bf z}\big{]}$
$\displaystyle=$ $\displaystyle\int_{\mathcal{X}}{\bf ME}_{\chi}\big{[}{\bf
E}_{\eta_{n-1}}f\big{(}\omega(\Delta
t)\big{)}P_{0}(d\chi)\int_{R^{d}}p^{1}_{\tau_{\bf-y-z}\eta_{n-1}}({\bf
0,z})d{\bf z}\big{]}$ $\displaystyle=$ $\displaystyle\int_{\mathcal{X}}{\bf
ME}_{\chi}\big{[}{\bf E}_{\eta_{n-1}}f\big{(}\omega(\Delta
t)\big{)}P_{0}(d\chi)\int_{R^{d}}p^{1}_{\eta_{n-1}}({\bf-y-z,-y})d{\bf
z}\big{]}$ $\displaystyle=$ $\displaystyle\int_{\mathcal{X}}{\bf
ME}_{\chi}\big{[}{\bf E}_{\eta_{n-1}}f\big{(}\omega(\Delta
t)\big{)}\big{]}P_{0}(d\chi)=\int_{\mathcal{X}}{\bf
ME}_{\chi}\big{[}f(\eta_{n-1})\big{]}P_{0}(d\chi),$ (34)

where $\tilde{\eta}_{1}$ is defined according to (27) but with initial
condition $\tilde{\eta}_{0}=\eta_{n-1}$. Thus, the first equation in (32) is
proved. The second equation in (32) is obvious according to the definition
(28) and $S_{n}$ is measure-preserving.

To prove the equation (33), let ${\bf Y}_{n}^{\omega}={\bf
X}^{\omega}_{n}+{\bf B}_{\Delta t}(\eta_{n})={\bf
X}^{\omega}_{n+1}-\sigma{\bm{\xi}}_{n}$. Then, we have

$\displaystyle\mathbb{E}f(\eta_{n+1})$ $\displaystyle={\bf
EME}_{\eta_{n}}f\big{(}\tau_{{\bf
Y}_{n}^{\omega}+\sigma\bm{\xi}_{n}}\omega(\Delta
t)\big{)}=\int_{\mathcal{X}}P_{0}(d\chi)\int_{R^{d}}p_{0}({\bf 0,z}){\bf
ME}_{\chi}{\bf E}_{\eta_{n}}f\big{(}\tau_{\bf z}\tau_{{\bf
Y}_{n}^{\omega}}\omega(\Delta t)\big{)}d{\bf z},$
$\displaystyle=\int_{\mathcal{X}}{\bf ME}_{\chi}{\bf
E}_{\eta_{n}}f\big{(}\tau_{{\bf Y}_{n}^{\omega}}\omega(\Delta
t)\big{)}P_{0}(d\chi)\int_{R^{d}}p_{0}({\bf 0,z})d{\bf z},$
$\displaystyle=\mathbb{E}f\Big{(}\tau_{{\bf X}^{\omega}_{n}+{\bf B}_{\Delta
t}(\eta_{n})}\omega\big{(}(n+1)\Delta t\big{)}\Big{)}.$ (35)

Notice that in the proof we use the property that $\tau$ is a measure-
preserving transformation. ∎

Equipped with these preparations, we can state the main results. The first
result is that the operator $S_{n}$ defined in (28) is a contractive map on
$L^{2}(\mathcal{X})$.

######  Theorem 3.6.

$S_{n}$ has the property that

$\displaystyle||S_{n}f||_{L^{2}(\mathcal{X})}\leq\exp(-2c_{1}n\Delta
t)||f||_{L^{2}(\mathcal{X})},$ (36)

for all $f\in L_{0}^{2}(\mathcal{X}).$

###### Proof.

We first consider the case when $n=1$. The key observation is that

$\displaystyle\int_{\mathcal{X}}S_{1}f(\chi)\cdot
S_{1}f(\chi)P_{0}(d\chi)=\int_{\mathcal{X}}{\bf E}_{\chi}f(\eta_{1})\cdot{\bf
E}_{\chi}f(\eta_{1})P_{0}(d\chi),$ $\displaystyle=$
$\displaystyle\int_{\mathcal{X}}P_{0}(d\chi)\int_{R^{d}}p^{1}_{\chi}({\bf
0,y}){\bf ME}_{\chi}f(\tau_{\bf y}\omega(\Delta t))d{\bf
y}\cdot\int_{R^{d}}p^{1}_{\chi}({\bf 0,y}){\bf ME}_{\chi}f(\tau_{\bf
y}\omega(\Delta t))d{\bf y},$ $\displaystyle\leq$
$\displaystyle\int_{\mathcal{X}}P_{0}(d\chi)\int_{R^{d}}p^{1}_{\chi}({\bf
0,y}){\bf E}_{\chi}f(\tau_{\bf y}\omega(\Delta t))\cdot{\bf
E}_{\chi}f(\tau_{\bf y}\omega(\Delta t))d{\bf y},$ $\displaystyle=$
$\displaystyle\int_{\mathcal{X}}{\bf E}_{\chi}f(\omega(\Delta t))\cdot{\bf
E}_{\chi}f(\omega(\Delta
t))P_{0}(d\chi)\int_{R^{d}}p^{1}_{\chi}({\bf-y,0})d{\bf y},$ $\displaystyle=$
$\displaystyle\int_{\mathcal{X}}P^{\Delta t}f(\chi)\cdot P^{\Delta
t}f(\chi)P_{0}(d\chi),$ (37)

where $P^{\Delta t}$ is a strongly continuous Markov semigroup on
$L^{2}(\mathcal{X})$. In the third line of (37), we use the fact that
$p^{1}_{\chi}({\bf 0,y})$ is a probability density function so we can easily
get the result by using the Cauchy-Schwarz inequality. Therefore, we obtain

$\displaystyle||S_{1}f||_{L^{2}(\mathcal{X})}\leq||P^{\Delta
t}f||_{L^{2}(\mathcal{X})}\leq\exp(-2c_{1}\Delta
t)||f||_{L^{2}(\mathcal{X})},$ (38)

where the exponential decay property (5) is used. The assertion in (36) can be
obtained if we repeat to use the above property $n$ times. ∎

Next, we define $\bar{\bf B}_{\Delta t}=\mathbb{E}{\bf B}_{\Delta t}$ and
$\tilde{\bf B}_{\Delta t}={\bf B}_{\Delta t}-\bar{\bf B}_{\Delta t}$. We aim
to get some estimates for the mean values $\bar{\bf B}_{\Delta t}$ and
$\mathbb{E}{\bf X}^{\omega}_{n}$, which are important in our convergence
analysis for the effective diffusivity later.

######  Theorem 3.7.

Under the Assumptions 3.1, 3.2 and 3.3, if we choose a volume-preserving
numerical scheme (20) to compute the sub-problem (18), where the local
truncation error is $O(\Delta t)^{2}$, then $\bar{\bf B}_{\Delta t}$ is of
order $O(\Delta t)^{2}$. In addition, $\mathbb{E}{\bf
X}^{\omega}_{n}-n\bar{\bf B}_{\Delta t}$ is bounded.

###### Proof.

By using a volume-preserving numerical scheme (with a local truncation error
$O(\Delta t)^{2}$) to compute (18), we have

$\displaystyle\mathbb{E}{\bf B}_{\Delta t}=\mathbb{E}\int_{0}^{\Delta t}{\bf
b}(t,X_{t}^{\omega},\omega)dt+O(\Delta t)^{2}=\mathbb{E}\int_{0}^{\Delta
t}{\bf b}(\eta^{0}_{t})dt+O(\Delta t)^{2},$ (39)

where $\eta_{t}^{0}$ is the environment process defined in (7) with
$\sigma=0$. Based on the regularity Assumption 3.1 for ${\bf b}$, although the
constant in the local truncation error $O(\Delta t)^{2}$ of the numerical
scheme (20) is random, it has a uniform upper bound. Thus, the error in
Eq.(39) is still of order $O(\Delta t)^{2}$ after taking the expectation.
Notice that when we define ${\bf B}_{\Delta t}$, we only consider the sub-
problem (18). Recall the fact that $S^{t}$ is measure-preserving, so we get

$\displaystyle\mathbb{E}\int_{0}^{\Delta t}{\bf
b}(\eta^{0}_{t})dt=\int_{0}^{\Delta t}\int_{\mathcal{X}}{\bf E}_{\chi}{\bf
b}(\eta^{0}_{t})dP_{0}(\chi)dt=\int_{0}^{\Delta t}\mathbb{E}S^{t}{\bf
b}dt=\int_{0}^{\Delta t}\mathbb{E}{\bf b}dt=0,$ (40)

where we have used the definition of $S^{t}$ in (8) and ${\bf b}$ is mean-
zero. Therefore, $\mathbb{E}{\bf B}_{\Delta t}$ is of the order $(\Delta
t)^{2}$. Moreover, from the numerical scheme (24) we have

$\displaystyle\mathbb{E}{\bf X}^{\omega}_{n}=\mathbb{E}{\bf
X}^{\omega}_{n-1}+\ {\mathbb{E}{\bf B}_{\Delta t}^{\omega((n-1)\Delta t)}({\bf
X}^{\omega}_{n-1})}=\mathbb{E}{\bf
X}^{\omega}_{0}+\sum_{i=0}^{n-1}\mathbb{E}S_{i}{\bf B}_{\Delta
t}=\mathbb{E}{\bf X}^{\omega}_{0}+\sum_{i=0}^{n-1}\mathbb{E}S_{i}\tilde{\bf
B}_{\Delta t}+n\bar{\bf B}_{\Delta t}.$ (41)

Under the Assumptions 3.1, 3.2 and 3.3, we know that $\bar{\bf B}_{\Delta t}$
and $\tilde{\bf B}_{\Delta t}$ are bounded. According to (36) in Theorem 3.6,
$||S_{i}\tilde{\bf B}_{\Delta t}||_{L^{2}(\mathcal{X})}$ decays exponentially
with respect to $i$, so we can easily verify that
$\sum_{i=0}^{n-1}S_{i}\tilde{\bf B}_{\Delta t}$ is bounded in
$L^{2}(\mathcal{X})$, which implies
$\big{|}\sum_{i=0}^{n-1}\mathbb{E}S_{i}\tilde{\bf B}_{\Delta
t}\big{|}<\infty$. Thus, we prove that $\mathbb{E}{\bf
X}^{\omega}_{n}-n\bar{\bf B}_{\Delta t}$ is bounded. ∎

###  3.3 A discrete-type corrector problem

The corrector problem (16) plays an important role in defining the effective
diffusivity for the random flow. To study the property of the numerical
solutions, we will define a discrete-type corrector problem and study the
property of its solution.

######  Theorem 3.8.

Let us define ${\bm{\psi}}_{\Delta t}=\sum_{i=0}^{\infty}S_{i}\tilde{\bf
B}_{\Delta t}$. Then, $\bm{\psi}_{\Delta t}$ is the unique solution of the
discrete-type corrector problem in $(L_{0}^{2}(\mathcal{X}))^{d}$ defined as
follows

$\displaystyle(S_{1}-I)\bm{\psi}_{\Delta t}=-\tilde{\bf B}_{\Delta t}.$ (42)

###### Proof.

The formulation of ${\bm{\psi}}_{\Delta t}$ solves the discrete-type corrector
problem (42) can be easily verified through simple calculations, i.e.,

$\displaystyle(S_{1}-I)\bm{\psi}_{\Delta t}=\sum_{i=1}^{\infty}S_{i}\tilde{\bf
B}_{\Delta t}-\sum_{i=0}^{\infty}S_{i}\tilde{\bf B}_{\Delta t}=-\tilde{\bf
B}_{\Delta t}.$ (43)

The property ${\bf E}{\bm{\psi}}_{\Delta t}=0$ is a straightforward result
from the formulation of ${\bm{\psi}}_{\Delta t}$. The uniqueness of the
solution comes from Theorem 3.6. Suppose the equation (42) has two different
solutions $\bm{\psi}_{1},\bm{\psi}_{2}\in L_{0}^{2}(\mathcal{X})$, we have
that $(S_{1}-I)(\bm{\psi}_{1}-\bm{\psi}_{2})=0$, then

$\displaystyle||\bm{\psi}_{1}-\bm{\psi}_{2}||_{L^{2}(\mathcal{X})}=||S_{1}(\bm{\psi}_{1}-\bm{\psi}_{2})||_{L^{2}(\mathcal{X})}\leq\exp(-2c_{1}\Delta
t)||\bm{\psi}_{1}-\bm{\psi}_{2}||_{L^{2}(\mathcal{X})},$

which implies that $\bm{\psi}_{1}-\bm{\psi}_{2}=0$. Thus, the uniqueness of
solution for Eq.(42) is proved. ∎

######  Remark 3.3.

The formulation of the discrete-type corrector problem (42) is equivalent to
the equation

$\displaystyle\mathbb{E}\big{[}\bm{\psi}^{\omega(i\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i})|\eta_{i-1}]-\bm{\psi}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})=-\tilde{\bf B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1}).$ (44)

This can be seen by replacing $\chi$ with $\eta_{n-1}$ in the definition of
$S_{1}$; see Eq.(28).

Finally, we study the regularity of the solution of the discrete-type
corrector problem (42). The following result is based on the regularity
assumption on the velocity field ${\bf b}$. Since we are interested in
statistical properties of the solution ${\bf X}(t)$, which only requires
convergence in law, we can choose smooth realizations of the velocity field
${\bf b}$.

######  Theorem 3.9.

Suppose ${\bf b}\in(C_{b}^{m}(\mathcal{X}))^{d}$, then $\bm{\psi}_{\Delta t}$
is in $(H^{m}(\mathcal{X}))^{d}$.

###### Proof.

First we prove that, under the assumption ${\bf
b}\in(C_{b}^{m}(\mathcal{X}))^{d}$ for $m\geq 1$, we have that for any $f\in
L^{2}(\mathcal{X})$, $S_{1}f\in H^{1}(\mathcal{X})$. Since

$\displaystyle S_{1}f(\tau_{\bf x}\chi)$
$\displaystyle=\int_{R^{d}}p_{\tau_{\bf x}\chi}^{1}({\bf 0,y})P^{\Delta
t}f(\tau_{\bf x+y}\chi)d{\bf y}=\int_{R^{d}}p_{\chi}^{1}({\bf x,x+y})P^{\Delta
t}f(\tau_{\bf x+y}\chi)d{\bf y},$ $\displaystyle=\int_{R^{d}}p_{\chi}^{1}({\bf
x,y})P^{\Delta t}f(\tau_{\bf y}\chi)d\bf y,$ (45)

where $p_{\chi}^{1}({\bf x,y})$ is the transition probability density defined
in (29). Notice that

$\displaystyle{\bf D}_{\bf x}p^{1}_{\chi}({\bf x,y})=\big{(}I+{\bf D}{\bf
B}_{\Delta t}^{\chi}({\bf x})\big{)}\big{(}{\bf y-x-B}_{\Delta t}^{\chi}({\bf
x})\big{)}p^{1}_{\chi}({\bf x,y})/\ {(\sigma^{2}\Delta t)},$ (46)

and ${\bf B}_{\Delta t}\in(C_{b}^{m}(\mathcal{X}))^{d}$, we can obtain that
$\int_{R^{d}}({\bf y-x-B}^{\chi}_{\Delta t}({\bf x}))_{i}^{2}p^{1}_{\chi}({\bf
x,y})d{\bf x}$ is uniformly bounded for almost all $\chi$. Here the indicator
$i$ represents the $i-th$ component. This concludes that

$\displaystyle\int_{R^{d}}{\bf D}_{\bf x}p_{\chi}^{1}({\bf x,y})P^{\Delta
t}f(\tau_{\bf y}\chi)d{\bf y}\in(L^{2}(\mathcal{X}))^{d}.$ (47)

The statement (47) implies that ${\bf D}S_{1}f\in(L^{2}(\mathcal{X}))^{d}$ by
the dominant convergence theorem. Thus $S_{1}f\in H^{1}(\mathcal{X})$.
According to the definition of the discrete-type corrector problem (42),
${\bm{\psi}}_{\Delta t}$ satisfies

$\displaystyle{\bm{\psi}}_{\Delta t}=S_{1}{\bm{\psi}}_{\Delta t}+\tilde{\bf
B}_{\Delta t}.$ (48)

Therefore, we obtain that $\bm{\psi}_{\Delta t}\in(H^{1}(\mathcal{X}))^{d}$.
Moreover, noticing that

$\displaystyle{\bf D}S_{1}f(\chi)$ $\displaystyle=\int_{R^{d}}{\bf D}_{\bf
x}p_{\chi}^{1}({\bf 0,y})P^{\Delta t}f(\tau_{\bf y}\chi)d\bf y,$
$\displaystyle=\int_{R^{d}}\big{(}I+{\bf D}{\bf B}_{\Delta
t}^{\chi}(0)\big{)}\big{(}{\bf y-0-B}_{\Delta
t}^{\chi}(0)\big{)}p^{1}_{\chi}({\bf x,y})P^{\Delta t}f(\tau_{y}\chi)d{\bf
y}/(\sigma^{2}\Delta t)$ $\displaystyle=\big{(}I+{\bf D}{\bf B}_{\Delta
t}^{\chi}(0)\big{)}\int_{R^{d}}-{\bf D}_{\bf y}p^{1}_{\chi}({\bf
0,y})P^{\Delta t}f(\tau_{\bf y}\chi)d{\bf y}/(\sigma^{2}\Delta t)$
$\displaystyle=\big{(}I+{\bf D}{\bf B}_{\Delta
t}^{\chi}(0)\big{)}\int_{R^{d}}p^{1}_{\chi}({\bf 0,y}){\bf D}_{\bf y}P^{\Delta
t}f(\tau_{\bf y}\chi)d{\bf y}/(\sigma^{2}\Delta t)$
$\displaystyle=(\sigma^{2}\Delta t)^{-1}\big{(}I+{\bf D}{\bf B}_{\Delta
t}^{\chi}(0)\big{)}S_{1}{\bf D}f(\chi).$ (49)

We arrive that

$\displaystyle{\bf D}\bm{\psi}_{\Delta t}=\ {(\sigma^{2}\Delta t)^{-1}}(I+{\bf
D}{\bf B}_{\Delta t})S_{1}{\bf D}\bm{\psi}_{\Delta t}+{\bf D}\tilde{\bf
B}_{\Delta t}.$ (50)

Similar argument shows that ${\bf D}\bm{\psi}_{\Delta
t}\in(H^{1}(\mathcal{X}))^{d\times d}$. Doing this argument recursively, we
prove that $\bm{\psi}_{\Delta t}$ is in $(H^{m}(\mathcal{X}))^{d}$. ∎

##  4 Convergence analysis

In this section, we shall prove the convergence rate of our stochastic
structure-preserving scheme in computing effective diffusivity. The
convergence analysis is based on a probabilistic approach, which allows us to
get rid of the exponential growth factor in the error estimate.

###  4.1 Convergence of the discrete-type corrector problem to the continuous
one

We first show that, if $\Delta t$ is small enough, $S^{\Delta t}$ will
converge to $S_{1}$. Moreover, the following statement holds.

######  Lemma 4.1.

If $f$ is a globally Lipschitz function with respect to ${\bf x}$, then we
have

$\displaystyle||S_{n}f-S^{n\Delta t}f||_{L^{2}(\mathcal{X})}\leq c_{2}L\Delta
t,$ (51)

where $L$ is the Lipschitz constant for $f$ and $c_{2}$ depends only on the
computational time $n\Delta t$.

###### Proof.

According to the definitions of the semigroups in (8) and (28), we have that
$(S_{n}-S^{n\Delta t})f(\chi)={\bf E}_{\chi}\big{(}f(\eta_{n})-f(\eta(n\Delta
t))\big{)}$, which implies

$\displaystyle(S_{n}-S^{n\Delta t})f(\chi)\leq L{\bf E}_{\chi}\big{|}{\bf
X}_{n}^{\omega}-{\bf X}_{n\Delta t}^{\omega}\big{|}.$ (52)

The error estimate for the Euler-Maruyama method has been intensively studied
in the literature (see e.g. [17, 24]). According to Assumption 3.1, the
regularity assumption for $\bf b$ is satisfied. Thus, the strong order of
accuracy of the Euler-Maruyama scheme for SDE driven by additive noise is 1,
i.e.,

$\displaystyle\sqrt{{\bf E}_{\chi}|{\bf X}_{n}^{\omega}-{\bf X}_{n\Delta
t}^{\omega}|^{2}}\leq c_{2}\Delta t.$ (53)

The proof is a simple application of Theorem 1.1 in [24]. We apply the
Jensen’s inequality for expectation and obtain

$\displaystyle{\bf E}_{\chi}|{\bf X}_{n}^{\omega}-{\bf X}_{n\Delta
t}^{\omega}|\leq\sqrt{{\bf E}_{\chi}|{\bf X}_{n}^{\omega}-{\bf X}_{n\Delta
t}^{\omega}|^{2}}\leq c_{2}\Delta t.$ (54)

Combining the estimate results in (52) and (54), we prove the assertion of
Lemma 4.1.

∎

Then, we show that under certain conditions the discrete-type corrector
problem converges to the continuous one, which facilitates the convergence
analysis of our numerical method in computing the effective diffusivity for
random flows.

######  Theorem 4.2.

The solution $\bm{\psi}_{\Delta t}$ converges to the solution $\bm{\psi}$ of
the continuous-type corrector problem defined in (15) in $L^{2}(\mathcal{X})$,
as $\Delta t\rightarrow 0$.

###### Proof.

Using the exponential decay properties of $S^{t}$ and $S_{n}$, we first choose
a truncation time $T_{0}$ and obtain the following two inequalities

$\displaystyle\big{|}\big{|}\int_{T_{0}-\Delta t}^{\infty}S^{t}{\bf
b}dt\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq\frac{1}{2c_{1}}\exp(-2c_{1}T_{0}),\quad\text{and}\quad\big{|}\big{|}\sum_{n=[T_{0}/\Delta
t]-1}^{\infty}S_{n}\tilde{\bf B}_{\Delta
t}\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq\frac{1}{2c_{1}}\exp(-2c_{1}T_{0}),$
(55)

where $c_{1}>0$ is defined in (11). Then, for any $\epsilon>0$, we choose
$T_{0}$ big enough such that $\frac{1}{c_{1}}\exp(-c_{1}T_{0})<\epsilon$.
Next, we estimate the error between $\sum_{n=0}^{N-1}S_{n}\tilde{\bf
B}_{\Delta t}$ and $\int_{0}^{N\Delta t}S^{t}{\bf b}dt$ for $N\leq
T_{0}/\Delta t$. We know that

$\displaystyle\big{|}\big{|}\int_{0}^{N\Delta t}S^{t}{\bf
b}dt-\sum_{n=0}^{N-1}S^{n\Delta t}{\bf b}\Delta
t\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq C_{1}\Delta t$ (56)

due to the strongly continuity of $S^{t}$ (see Prop. 2.2) and

$\displaystyle\big{|}\big{|}\sum_{n=0}^{N-1}S_{n}\tilde{\bf B}_{\Delta
t}-\sum_{n=0}^{N-1}S^{n\Delta t}{\bf b}\Delta
t\big{|}\big{|}_{L^{2}(\mathcal{X})}$
$\displaystyle\leq\big{|}\big{|}\sum_{n=0}^{N-1}S_{n}\tilde{\bf B}_{\Delta
t}-\sum_{n=0}^{N-1}S_{n}{\bf b}\Delta t\big{|}\big{|}_{L^{2}(\mathcal{X})}$
$\displaystyle+\big{|}\big{|}\sum_{n=0}^{N-1}S_{n}{\bf b}\Delta
t-\sum_{n=0}^{N-1}S^{n\Delta t}{\bf b}\Delta
t\big{|}\big{|}_{L^{2}(\mathcal{X})}.$ (57)

We can estimate the two terms of the right hand side of the inequality (57)
separately. Since the local truncation error of the numerical scheme (20) is
at least second order, we have $\big{|}\big{|}\tilde{\bf B}_{\Delta t}-{\bf
b}\Delta t\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq O(\Delta t)^{2}$. From Lemma
4.1, we know $\big{|}\big{|}(S_{n}-S^{n\Delta t}){\bf b}\Delta
t\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq\ {O(\Delta t)^{2}}$ for all $n\leq
N$. This gives the estimate

$\displaystyle\big{|}\big{|}\sum_{n=0}^{N-1}S_{n}\tilde{\bf B}_{\Delta
t}-\sum_{n=0}^{N-1}S^{n\Delta t}{\bf b}\Delta
t\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq c_{2}N(\Delta t)^{2}\leq
c_{2}T_{0}\Delta t.$ (58)

Finally, we take $\Delta t\leq\epsilon/(c_{2}T_{0})$ and obtain

$\displaystyle\big{|}\big{|}\int_{0}^{\infty}S^{t}{\bf
b}dt-\sum_{n=0}^{\infty}S_{n}\tilde{\bf B}_{\Delta
t}\big{|}\big{|}_{L^{2}(\mathcal{X})}\leq 2\epsilon+O(\epsilon^{2}).$ (59)

We prove the assertion of the Theorem. ∎

######  Remark 4.1.

The constant $c_{2}$ in Lemma 4.1 actually exponentially depends on $T_{0}$,
i.e., $c_{2}=\exp(c_{3}T_{0})$ with $c_{3}>0$. To balance each value of
$\epsilon$, we have
$\frac{1}{2c_{1}}\exp(-2c_{1}T_{0})=\exp(c_{3}T_{0})T_{0}\Delta t$, which
requires $T_{0}\approx-1/(2c_{1}+c_{3})\log\Delta t$ and
$\epsilon\approx\frac{1}{c_{1}}\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}}$.

###  4.2 Convergence of the numerical method in computing effective
diffusivity

Now we are in a position to show the main results of our paper. We prove that
the effective diffusivity obtained by our numerical method converges to the
exact one defined in (17).

######  Theorem 4.3.

Let ${\bf X}^{\omega}_{n}$, $n=0,1,....$ be the numerical solution of the
stochastic structure-preserving scheme (24) and $\Delta t$ be the time step
that is fixed. Let $\bar{\bf X}^{\omega}_{n}={\bf X}^{\omega}_{n}-n\bar{\bf
B}_{\Delta t}$. We have the convergence estimate of the numerical method in
computing effective diffusivity as

$\displaystyle\frac{\mathbb{E}\bar{\bf X}^{\omega}_{n}\otimes\bar{\bf
X}^{\omega}_{n}}{n\Delta t}=\sigma^{2}{\bf I}_{d}+2{\bf
S}\int_{\mathcal{X}}\bm{\psi}\otimes{\bf b}dP_{0}+\rho(\Delta
t)+O(\frac{1}{\sqrt{n}\Delta t}),$ (60)

where $\rho(\Delta t)=O(\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}})$ is a function
satisfying $\lim_{\Delta t\to 0}\rho(\Delta t)=0$ and is independent of the
computational time $T$, and $\bf S$ represents the symmetrization operator on
a matrix, i.e., $\bf SA=\frac{A+A^{T}}{2}$.

###### Proof.

First of all, from direct computations we obtain that

$\displaystyle\mathbb{E}\bar{\bf X}^{\omega}_{n}\otimes\bar{\bf
X}^{\omega}_{n}=\mathbb{E}\big{(}\bar{\bf X}^{\omega}_{n-1}+\tilde{\bf
B}^{\omega((n-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n-1})+\sigma\bm{\xi}_{n-1}\big{)}\otimes\big{(}\bar{\bf
X}^{\omega}_{n-1}+\tilde{\bf B}^{\omega((n-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n-1})+\sigma\bm{\xi}_{n-1}\big{)},$ $\displaystyle=$
$\displaystyle\mathbb{E}\bar{\bf X}^{\omega}_{n-1}\otimes\bar{\bf
X}^{\omega}_{n-1}+\sigma^{2}{\bf I}_{d}\Delta t+2{\bf S}\mathbb{E}\bar{\bf
X}^{\omega}_{n-1}\otimes\tilde{\bf B}^{\omega((n-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n-1})+\mathbb{E}\tilde{\bf B}^{\omega((n-1)\Delta t)}_{\Delta
t}(\bar{\bf X}^{\omega}_{n-1})\otimes\tilde{\bf B}^{\omega((n-1)\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{n-1}),$ $\displaystyle=$
$\displaystyle\mathbb{E}\bar{\bf X}^{\omega}_{0}\otimes\bar{\bf
X}^{\omega}_{0}+\sigma^{2}{\bf I}_{d}n\Delta t+2\sum_{i=1}^{n}{\bf
S}\mathbb{E}\bar{\bf X}^{\omega}_{i-1}\otimes\tilde{\bf B}^{\omega((i-1)\Delta
t)}_{\Delta t}(\bar{\bf X}^{\omega}_{i-1})+\sum_{i=1}^{n}\mathbb{E}\tilde{\bf
B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf X}^{\omega}_{i-1})\otimes\tilde{\bf
B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf X}^{\omega}_{i-1}),$ (61)

where we have used the facts that $\bm{\xi}_{n-1}$ is independent with
$\bar{\bf X}^{\omega}_{n-1}$ and
$\mathbb{E}{\bm{\xi}}_{n-1}\otimes{\bm{\xi}}_{n-1}=\Delta t{\bf I}_{d}$.

The first two terms on the right hand side of Eq.(61) are easy to handle since
each entry in $\frac{\mathbb{E}\bar{\bf X}^{\omega}_{0}\otimes\bar{\bf
X}^{\omega}_{0}}{n\Delta t}$ is $O(\frac{1}{n\Delta t})$ and
$\frac{\sigma^{2}{\bf I}_{d}n\Delta t}{n\Delta t}=\sigma^{2}{\bf I}_{d}$. For
the forth term on the right hand side of Eq.(61), using the property that
$S_{i}$ is measure-preserving; see Theorem 3.4 and Assumption 3.3 , we can get

$\displaystyle\frac{1}{n\Delta t}\sum_{i=1}^{n}\mathbb{E}\tilde{\bf
B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf X}^{\omega}_{i-1})\otimes\tilde{\bf
B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf X}^{\omega}_{i-1})$
$\displaystyle=$ $\displaystyle\frac{1}{n\Delta
t}\sum_{i=1}^{n}\mathbb{E}S_{i-1}(\tilde{\bf B}_{\Delta t}\otimes\tilde{\bf
B}_{\Delta t})=\frac{1}{n\Delta t}n\mathbb{E}\tilde{\bf B}_{\Delta
t}\otimes\tilde{\bf B}_{\Delta t}=O(\Delta t).$ (62)

We will focus on the third term on the right hand side of Eq.(61), which
corresponds to the strengthen of the convection-enhanced diffusion and is the
most difficult term to deal with. Substituting the formulation of the
discrete-type corrector problem (44) into it, we obtain

$\displaystyle\sum_{i=1}^{n}\mathbb{E}\bar{\bf
X}^{\omega}_{i-1}\otimes\tilde{\bf B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})=-\sum_{i=1}^{n}\mathbb{E}\bar{\bf
X}^{\omega}_{i-1}\otimes\big{(}\mathbb{E}[\bm{\psi}^{\omega(i\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i})|\eta_{i-1}]-\bm{\psi}^{\omega((i-1)\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i-1})\big{)}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}\mathbb{E}\mathbb{E}\Big{[}\bar{\bf
X}^{\omega}_{i-1}\otimes\big{(}\bm{\psi}^{\omega(i\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i})-\bm{\psi}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})\big{)}\big{|}\eta_{i-1}\Big{]}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}\mathbb{E}(\bar{\bf X}^{\omega}_{i-1}-\bar{\bf
X}^{\omega}_{i})\otimes\bm{\psi}^{\omega(i\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i})+\mathbb{E}\bar{\bf
X}^{\omega}_{0}\otimes\bm{\psi}^{\omega(0)}_{\Delta t}({\bf
X}^{\omega}_{0})-\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\mathbb{E}\big{(}\tilde{\bf B}^{\omega((i-1)\Delta
t}_{\Delta t}({\bf
X}^{\omega}_{i-1})+\sigma\bm{\xi}_{i-1}\big{)}\otimes\bm{\psi}^{\omega(i\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i})+\mathbb{E}\bar{\bf
X}^{\omega}_{0}\otimes\bm{\psi}^{\omega(0)}_{\Delta t}({\bf
X}^{\omega}_{0})-\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n}).$ (63)

Here from the first row to the second row, we use the fact that $\bar{\bf
X}^{\omega}_{i-1}$ and $\bm{\psi}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})$ are measurable in the $\sigma$-algebra generated by
$\eta_{i-1}$. From the second row to the third row, we use the property of
conditional expectation and Abel’s summation formula.

Let us first estimate the summation term on the right hand side of Eq.(63).
For each index $i$, we have

$\displaystyle\mathbb{E}\big{(}\tilde{\bf B}^{\omega((i-1)\Delta t)}_{\Delta
t}({\bf
X}^{\omega}_{i-1})+\sigma{\bm{\xi}}_{i-1}\big{)}\otimes\bm{\psi}^{\omega(i\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i}),$ $\displaystyle=$
$\displaystyle\mathbb{E}\tilde{\bf B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})\otimes\bm{\psi}^{\omega(i\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i})+\mathbb{E}\sigma{\bm{\xi}}_{i-1}\otimes\bm{\psi}^{\omega(i\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i}).$ (64)

Through simple calculations, we can show that the second term of the right
hand side of Eq.(64) is zero. Specifically, we have

$\displaystyle\mathbb{E}\sigma{\bm{\xi}}_{i-1}\otimes\bm{\psi}^{\omega(i\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i})$ $\displaystyle=$
$\displaystyle\mathbb{E}\sigma{\bm{\xi}}_{i-1}\otimes\bm{\psi}^{\omega(i\Delta
t)}_{\Delta t}\big{(}{\bf X}^{\omega}_{i-1}+{\bf B}^{\omega((i-1)\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{i-1})+\sigma{\bm{\xi}}_{i-1}\big{)},$
$\displaystyle=$ $\displaystyle\int_{\mathcal{{X}}}\int_{R^{d}}p_{0}({\bf
0,y})\sigma{\bf y}\otimes{\bf ME}_{\chi}\bm{\psi}_{\Delta
t}\big{(}\tau_{\sigma\bf y}\tau_{{\bf X}^{\omega}_{i-1}+{\bf B}_{\Delta
t}(\eta_{i-1})}\omega(i\Delta t)\big{)}d{\bf y}P_{0}(d\chi),$ $\displaystyle=$
$\displaystyle\int_{R^{d}}p_{0}({\bf 0,y})\sigma{\bf
y}\otimes\int_{\mathcal{{X}}}{\bf ME}_{\chi}\bm{\psi}_{\Delta
t}\big{(}\tau_{\sigma\bf y}\tau_{{\bf X}^{\omega}_{i-1}+{\bf B}_{\Delta
t}(\eta_{i-1})}\omega(i\Delta t)\big{)}P_{0}(d\chi)d{\bf y},$ $\displaystyle=$
$\displaystyle\int_{R^{d}}p_{0}({\bf 0,y}){\sigma\bf y}\otimes{\bf
E}\bm{\psi}_{\Delta t}d{\bf y}=0.$ (65)

Here, the expectation is taken over all the randomness in the system. In the
third row of Eq.(65), ${\bf y}$ is a realization of ${\bm{\xi}}_{i-1}$ and
$p_{0}({\bf 0,y})d{\bf y}$ is the measure associated with the Brownian motion,
while $P_{0}(d\chi)$ is the measure associated with the randomness in the
velocity field and initial data. The Fubini’s theorem is used in the fourth
row of Eq.(65) to switch the order of integration. The fifth row of (65) is
derived from Lemma 3.5; see Eq.(33). And ${\bf E}\bm{\psi}_{\Delta t}=0$ since
the solution of the discrete-type corrector problem is mean-zero; see Theorem
3.8.

Then, we compute the first term on the right hand side of Eq.(64) as follows,

$\displaystyle\mathbb{E}\tilde{\bf B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})\otimes\bm{\psi}^{\omega(i\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i})=\mathbb{E}\tilde{\bf B}_{\Delta
t}(\eta_{i-1})\otimes\bm{\psi}_{\Delta t}(\eta_{i})$ $\displaystyle=$
$\displaystyle\mathbb{E}\mathbb{E}\big{(}\tilde{\bf B}_{\Delta
t}(\eta_{i-1})\otimes\bm{\psi}_{\Delta
t}(\eta_{i})\big{|}\eta_{i-1}\big{)}=\mathbb{E}\tilde{\bf B}_{\Delta
t}(\eta_{i-1})\otimes\mathbb{E}\big{(}\bm{\psi}_{\Delta
t}(\eta_{i})\big{|}\eta_{i-1}\big{)}$ $\displaystyle=$
$\displaystyle\mathbb{E}\tilde{\bf B}_{\Delta t}(\eta_{i-1})\otimes
S_{1}\bm{\psi}_{\Delta t}(\eta_{i-1})=\mathbb{E}\tilde{\bf B}_{\Delta
t}(\eta_{i-1})\otimes\big{(}\bm{\psi}_{\Delta t}(\eta_{i-1})-\tilde{\bf
B}_{\Delta t}(\eta_{i-1})\big{)}$ $\displaystyle=$
$\displaystyle\mathbb{E}\tilde{\bf B}_{\Delta
t}(\eta_{i-1})\otimes\bm{\psi}_{\Delta t}(\eta_{i-1})-\mathbb{E}\tilde{\bf
B}_{\Delta t}(\eta_{i-1})\otimes\tilde{\bf B}_{\Delta t}(\eta_{i-1})$
$\displaystyle=$ $\displaystyle\mathbb{E}S_{i-1}(\tilde{\bf B}_{\Delta
t}\otimes\bm{\psi}_{\Delta t})-\mathbb{E}S_{i-1}(\tilde{\bf B}_{\Delta
t}\otimes\tilde{\bf B}_{\Delta t}).$ (66)

Using the property that each $S_{i-1}$ is measure-preserving; see Theorem 3.4,
we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\tilde{\bf
B}^{\omega((i-1)\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i-1})\otimes\bm{\psi}^{\omega(i\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{i})=\mathbb{E}\tilde{\bf B}_{\Delta t}\otimes\bm{\psi}_{\Delta
t}-\mathbb{E}\tilde{\bf B}_{\Delta t}\otimes\tilde{\bf B}_{\Delta t}.$ (67)

The term $\mathbb{E}\tilde{\bf B}_{\Delta t}\otimes\bm{\psi}_{\Delta t}$ in
Eq.(67) is corresponding to the strengthen of the convection-enhanced
diffusion. The term $\mathbb{E}\tilde{\bf B}_{\Delta t}\otimes\tilde{\bf
B}_{\Delta t}$ in Eq.(67) is of the order $O(\Delta t)^{2}$ due to Assumption
3.3 . This completes the estimate of the first term on the right hand side of
Eq.(63).

Now, we estimate the second term and third term on the right hand side of
Eq.(63). The term $\mathbb{E}\bar{\bf
X}^{\omega}_{0}\otimes\bm{\psi}^{\omega(0)}_{\Delta t}({\bf X}^{\omega}_{0})$
does not depend on $n$ and is bounded. For the third term, we want to prove
that

$\frac{1}{n\Delta t}\big{|}\big{|}\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n})\big{|}\big{|}\leq O(\frac{1}{\sqrt{n}\Delta t}),$ (68)

where $||\cdot||$ is a matrix norm. By using the Holder’s inequality, we know
that each entry of the matrix $\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n})$ satisfies

$\big{|}\mathbb{E}(\bar{\bf
X}^{\omega}_{n})_{l}\big{(}\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n})\big{)}_{j}\big{|}\leq\big{(}\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2}\big{)}^{1/2}\big{(}\mathbb{E}[(\bm{\psi}^{\omega(n\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{n}))_{j}]^{2}\big{)}^{1/2},\quad 1\leq
l,j\leq d.$ (69)

Again, using the property that $S_{n}$ is measure-preserving; see Theorem 3.4,
we have

$\mathbb{E}\big{[}(\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n}))_{j}\big{]}^{2}=\mathbb{E}\big{(}\bm{\psi}_{\Delta
t,j}(\eta_{n})\big{)}^{2}=\mathbb{E}S_{n}(\bm{\psi}_{\Delta
t,j})^{2}=\mathbb{E}(\bm{\psi}_{\Delta t,j})^{2},$ (70)

which is bounded since $\bm{\psi}_{\Delta t}\in(L_{0}^{2}(\mathcal{X}))^{d}$
according to Theorem 3.8. Thus, if we can prove
$\frac{1}{n}\mathbb{E}[(\bar{\bf X}^{\omega}_{n})_{l}]^{2}$ is bounded, then

$\frac{1}{n\Delta t}\big{|}\mathbb{E}(\bar{\bf
X}^{\omega}_{n})_{l}(\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n}))_{j}\big{|}\leq\frac{1}{\sqrt{n}\Delta
t}(\frac{1}{n}\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2})^{1/2}(\mathbb{E}[(\bm{\psi}^{\omega(n\Delta
t)}_{\Delta t}({\bf
X}^{\omega}_{n}))_{j}]^{2})^{1/2}=O(\frac{1}{\sqrt{n}\Delta t}).$ (71)

In order to prove that $\frac{1}{n}\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2}$ is bounded, we apply the AM-GM inequality on the
diagonal entries of the matrix $\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n})$ and obtain,

$\displaystyle\big{|}\mathbb{E}(\bar{\bf
X}^{\omega}_{n})_{l}(\bm{\psi}^{\omega(n\Delta t)}_{\Delta t}({\bf
X}^{\omega}_{n}))_{l}\big{|}\leq\epsilon\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2}+(4\epsilon)^{-1}\mathbb{E}[(\bm{\psi}^{\omega(n\Delta
t)}_{\Delta t}({\bf X}^{\omega}_{n}))_{l}]^{2},\quad 1\leq l\leq d,$ (72)

where $0<\epsilon<1$. The choice of $\epsilon$ will be specified later.

According to (69), we only need to estimate the terms $\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2}$. We first substitute the estimated result (72) into
Eq.(63), and then substitute the estimated results of (63) (including
Eqns.(65)(67)) into Eq.(61). Combining all the estimate results for terms on
the right hand side of Eq.(61), we obtain an estimate for
$\mathbb{E}[(\bar{\bf X}^{\omega}_{n})_{l}]^{2}$ as follows,

$\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2}\leq(\textbf{R}_{n})_{l}+\epsilon\mathbb{E}[(\bar{\bf
X}^{\omega}_{n})_{l}]^{2},$ (73)

where $(\textbf{R}_{n})_{l}$ denotes all the remaining terms with
$(\textbf{R}_{n})_{l}=O(n)$. Notice that $(\textbf{R}_{n})_{l}$ also contains
the term $(4\epsilon)^{-1}\mathbb{E}[(\bm{\psi}^{\omega(n\Delta t)}_{\Delta
t}({\bf X}^{\omega}_{n}))_{l}]^{2}$, which is $O(1)$ due to Eq.(70) and the
choice of $\epsilon$. We choose $0<\epsilon<1$ (e.g. $\epsilon=1/3$), move the
term $\epsilon\mathbb{E}[(\bar{\bf X}^{\omega}_{n})_{l}]^{2}$ to the left hand
side of (73), and divide both sides of the inequality by $(1-\epsilon)n$. We
can obtain that $\frac{1}{n}\mathbb{E}||\bar{\bf X}_{n}^{\omega}||^{2}$ is
bounded. Therefore, we prove the claim in (68).

Finally, we combine the estimate results in Eqns.(61)(62)(63)(67)(71) and
obtain that

$\displaystyle\frac{\mathbb{E}\bar{\bf X}^{\omega}_{n}\otimes\bar{\bf
X}^{\omega}_{n}}{n\Delta t}=\sigma^{2}{\bf I}_{d}+2{\bf
S}\mathbb{E}\bm{\psi}_{\Delta t}\otimes\tilde{\bf B}_{\Delta t}/\Delta
t+O(\Delta t)+O(\frac{1}{\sqrt{n}\Delta t}).$ (74)

According to Theorem 4.2 and Remark 4.1, we have the estimate

$\displaystyle\big{|}\big{|}2{\bf S}\mathbb{E}\bm{\psi}_{\Delta
t}\otimes\tilde{\bf B}_{\Delta t}/\Delta t-2{\bf
S}\mathbb{E}\bm{\psi}\otimes{\bf
b}\big{|}\big{|}_{L^{2}(\mathcal{X})}=O(\Delta
t^{\frac{2c_{1}}{2c_{1}+c_{3}}}):=\rho(\Delta t),$ (75)

where $\lim_{\Delta t\to 0}\rho(\Delta t)=0$. Thus, the statement in (60) is
proved. ∎

######  Remark 4.2.

Theorem 4.3 shows that when the time step $\Delta t$ is given and fixed, we
have

$\displaystyle\lim_{n\to\infty}\frac{\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bar{\bf X}^{\omega}_{n}}{n\Delta t}=\sigma^{2}{\bf
I}_{d}+2{\bf S}\int_{\mathcal{X}}\bm{\psi}\otimes{\bf b}dP_{0}+\rho(\Delta
t),$ (76)

which reveals the connection of the definition of the effective diffusivity by
solving discrete-type and continuous-type corrector problems. Our result
appears to be the first one in the literature to build this connection.

Notice that in the Theorem 4.3, we assume $\bar{\bf X}^{\omega}_{n}={\bf
X}^{\omega}_{n}-n\bar{\bf B}_{\Delta t}$ are given, where we use Monte Carlo
method to compute $\bar{\bf B}_{\Delta t}$. In some cases, if we cannot
calculate the drift constant $\bar{\bf B}_{\Delta t}$ exactly, we can directly
estimate the term $\mathbb{E}{\bf X}^{\omega}_{n}\otimes{\bf X}^{\omega}_{n}$,
which is summarized in the following corollary.

######  Corollary 4.4.

Let ${\bf X}^{\omega}_{n}$, $n=0,1,....$ be the numerical solution of the
stochastic structure-preserving scheme (24) and $\Delta t$ be the time step
that is fixed. Suppose $n(\Delta t)^{3}$ and $\frac{1}{\sqrt{n}\Delta t}$ are
small enough, we have

$\displaystyle\frac{\mathbb{E}{\bf X}^{\omega}_{n}\otimes{\bf
X}^{\omega}_{n}}{n\Delta t}=\sigma^{2}{\bf I}_{d}+2{\bf
S}\int_{\mathcal{X}}\bm{\psi}\otimes{\bf b}dP_{0}+\rho(\Delta
t)+O(\frac{1}{\sqrt{n}\Delta t})+O\big{(}n(\Delta t)^{3}\big{)},$ (77)

where $\rho(\Delta t)=O(\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}})$ is a function
satisfying $\lim_{\Delta t\to 0}\rho(\Delta t)=0$ and is independent of the
computational time $T$, and ${\bf S}$ represents the symmetrization operator.

###### Proof.

Using the observation that

$\displaystyle\frac{\mathbb{E}{\bf X}^{\omega}_{n}\otimes{\bf
X}^{\omega}_{n}}{n\Delta t}=\frac{\mathbb{E}\bar{\bf
X}^{\omega}_{n}\otimes\bar{\bf X}^{\omega}_{n}}{n\Delta t}+\frac{2{\bf
S}\mathbb{E}\bar{\bf X}^{\omega}_{n}\otimes\bar{\bf B}_{\Delta t}}{\Delta
t}+\frac{n^{2}\bar{\bf B}_{\Delta t}\otimes\bar{\bf B}_{\Delta t}}{n\Delta t}$
(78)

and Theorem 3.7, we can straightforwardly get the proof. ∎

######  Remark 4.3.

In our convergence analysis, we interpret the solution process generated by
our numerical scheme as a Markov process. By exploring the ergodicity of the
solution process (i.e., Markov process), we give a sharp error estimate of the
proposed numerical scheme in computing effective diffusivity.

######  Remark 4.4.

If ${\bf b}$ satisfies a Hölder-$\gamma$ continuous condition in the time
domain with $0<\gamma<1$, we will obtain a weaker convergence rate in Lemma
4.1, i.e. $||S_{n}f-S^{n\Delta t}f||_{L^{2}(\mathcal{X})}\leq c_{2}L(\Delta
t)^{\gamma}$ with $0<\gamma<1$. Under such condition, we can still obtain
convergence analysis of the numerical methods for computing effective
diffusivity, e.g. Theorem 4.3 with a smaller convergence rate in $\rho(\Delta
t)=O(\Delta t^{\frac{2\gamma c_{1}}{2c_{1}+c_{3}}})$.

##  5 Numerical Results

The aim of this section is two-fold. First, we will verify the convergence
results obtained in Section 4.2. Second, we will use the proposed method to
compute effective diffusivity in random flows, where incompressible random
flows in two- and three-dimensional spaces will be studied. Without loss of
generality, we compute the quantity $\frac{\mathbb{E}[({\bf
X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$, which is used to approximate
$D^{E}_{11}$ in the effective diffusivity matrix $D^{E}$. Notice that ${\bf
X}^{\omega}_{n,1}$ is the first component of the solution vector ${\bf
X}^{\omega}_{n}$. One can obtain $D^{E}_{11}$ by choosing ${\bf v}=(1,0)^{T}$
in the equation (17) of the Prop. 2.5.

###  5.1 Numerical methods for generating random flows

To start with, we discuss how to generate random flows that will be used in
our numerical experiments. Assume the vector field ${\bf b}(t,{\bf x},\omega)$
has a spectral measure

$\displaystyle\exp(-r({\bf k})|t|)\Gamma({\bf k})(\bf I-\frac{k\otimes
k}{|k|^{2}}),$ (79)

where ${\bf k}=(k_{1},k_{2})^{T}$ or ${\bf k}=(k_{1},k_{2},k_{3})^{T}$,
$r({\bf k})>c_{0}$ for some positive constant $c_{0}$, and $\Gamma({\bf k})$
is integrable and decays fast for large $\bf k$. Under such settings, the
velocity field ${\bf b}(t,{\bf x},\omega)$ satisfies the $\rho$ mixing
condition and is stationary and divergence-free [30, 7]. In order to mimic the
energy spectrum of real flows, we assume $\Gamma({\bf k})\propto 1/|{\bf
k}|^{2\alpha+d-2}$ with ultraviolet cutoff $|{\bf k}|\leq K<\infty$ and
$r({\bf k})\propto|{\bf k}|^{2\beta}$. The spectral gap condition 4 requires
$\beta\leq 0$ and the integrability of $\Gamma({\bf k})$ requires $\alpha<1$.
Here for simplicity, we choose $\beta=0$.

Given the spectral measure (79), we use the randomization method [19, 22] to
generate realizations of the velocity field. Specifically, we approximate it
as

$\displaystyle{\bf b}(t,{\bf x})=\frac{1}{\sqrt{M}}\sum_{m=1}^{M}\big{[}{\bf
u}_{m}\cos({\bf k}_{m}\cdot{\bf x})+{\bf v}_{m}\sin({\bf k}_{m}\cdot{\bf
x})\big{]}.$ (80)

Notice that we have suppressed the dependence of the velocity on $\omega$ for
notation simplicity here. In fact, the parameters ${\bf k}_{m}$, ${\bf u}_{m}$
and ${\bf v}_{m}$ contain randomness. The spectrum points ${\bf k}_{m}$ were
chosen independently according to the spectral measure $\Gamma(\bf k)$. Due to
the isotropicity, we first generate a point uniformly distributed on the unit
sphere or unit circle, which represents the direction of the ${\bf k}_{m}$.
Then we generate the length $r$ of ${\bf k}_{m}$, which satisfies a density
function $\rho(r)\propto 1/r^{2\alpha-1},0<r\leq K$.

For the random flows in two-dimensional space, we have

$\displaystyle{\bf u}_{m}=\xi_{m}(t)\frac{{\bf k}^{\perp}_{m}}{|{\bf
k}^{\perp}_{m}|},\quad{\bf v}_{m}=\eta_{m}(t)\frac{{\bf k}^{\perp}_{m}}{|{\bf
k}^{\perp}_{m}|},\quad{\bf k}_{m}=(k_{m}^{1},k_{m}^{2}),\quad m=1,...,M,$ (81)

where ${\bf k}_{m}^{\perp}=(-k_{m}^{2},k_{m}^{1})$, $\xi_{m}(t)$ and
$\eta_{m}(t)$ are independent 1D Ornstein-Uhlenbeck (OU) processes with
covariance function

$\displaystyle
Cov(\xi_{m}(t_{1}),\xi_{m}(t_{2}))=Cov(\eta_{m}(t_{1}),\eta_{m}(t_{2}))=\exp(-\theta|t_{1}-t_{2}|).$

Here $\theta>0$ is a parameter to control the roughness of the OU process. To
obtain the OU path for $\xi_{m}(t)$, we generate a series of
$\\{\xi_{m}(n\Delta t)\\}$ satisfies

$\displaystyle\xi_{m}(n\Delta t)=e^{-\theta\Delta t}\xi_{m}((n-1)\Delta
t)+\sqrt{1-e^{-2\theta\Delta t}}\zeta_{m,n},\quad n=1,2,3,...$ (82)

where $\xi_{m}(0)$, $\zeta_{m,n}$, $m,n=1,2,3,...$ are i.i.d. $N(0,1)$
distributed random variables. One can easily verify that $Cov(\xi_{m}(i\Delta
t),\xi_{m}(j\Delta t))=\exp(-\theta|i-j|\Delta t)$. The OU path for
$\eta_{m}(t)$ can be generated by using the same approach.

For the random flows in three-dimensional space, we have

$\displaystyle{\bf u}_{m}={\bm{\xi}}_{m}(t)\times\frac{{\bf k}_{m}}{|{\bf
k}_{m}|},\quad{\bf v}_{m}={\bm{\eta}}_{m}(t)\times\frac{{\bf k}_{m}}{|{\bf
k}_{m}|},\quad{\bf k}_{m}=(k_{m}^{1},k_{m}^{2},k_{m}^{3}),$ (83)

where the samples ${\bm{\xi}}_{m}(t)$ and ${\bm{\eta}}_{m}(t)$ are independent
3D random vectors, whose components are independent stationary OU process
having the covariance function
$Cov(\bm{\xi}_{m}(t_{1}),\bm{\xi}_{m}(t_{2}))=Cov(\bm{\eta}_{m}(t_{1}),\bm{\eta}_{m}(t_{2}))=\exp(-\theta|t_{1}-t_{2}|){\bf
I}_{3}$. Each component of ${\bm{\xi}}_{m}(t)$ and ${\bm{\eta}}_{m}(t)$ can be
generated by using the method (82). One can easily verify that the random
velocity fields generated by Eq.(80) with the setting (81) in two-dimensional
space and (83) in three-dimensional space automatically satisfy the
divergence-free condition.

###  5.2 Verification of the convergence analysis

In this subsection, we study the convergence rate of our method in computing
incompressible random flow in 2D and 3D spaces.

For the random flow in 2D space, we solve the SDE (1), where the velocity
filed is chosen as (80) with the setting (81). The velocity field were
simulated with $M=1000$. The parameters in the spectral measure $\Gamma(\bf
k)$ are $K=10$ and $\alpha=0.75$. The time-mixing constant $\theta=10$ in the
covariance function. The molecular diffusivity $\sigma=0.1$. We use Monte
Carlo method to generate independent samples for the Brownian motion ${\bf
w}(t)$ and velocity field ${\bf b}(t,{\bf x})$. The sample number is denoted
by $N_{mc}$.

We choose time step $\Delta t_{ref}=0.001$ and $N_{mc}=100,000$ to solve the
SDE (1) and compute the reference solution, i.e., the “exact” effective
diffusivity, where the final computational time is $T=22$ so that the
calculated effective diffusivity converges to a constant. It takes about 24
hours to compute the reference solution on a 64-core server (Gridpoint System
at HKU). The reference solution for the effective diffusivity is
$D^{E}_{11}=0.1736$.

For the random flow in 3D space, we solve the SDE (1), where the velocity
field is chosen as (80) with the setting (83). The velocity field were
simulated with $M=100$. The parameters in the spectral measure $\Gamma(\bf k)$
are $K=10$ and $\alpha=0.75$. The time-mixing constant $\theta=10$ in the
covariance function. The molecular diffusivity $\sigma=0.1$. Again, we use
Monte Carlo method to generate dependent samples for the Brownian motion ${\bf
w}(t)$ and velocity field ${\bf b}(t,{\bf x})$.

We choose $\Delta t_{ref}=0.001$ and $N_{mc}=180,000$ to solve the SDE (1) and
compute the reference solution, i.e., the “exact” effective diffusivity, where
the final computational time is $T=25$ so that the calculated effective
diffusivity converges to a constant. It takes about 21 hours to compute the
reference solution on a 64-core server (Gridpoint System at HKU). The
reference solution for the effective diffusivity is $D^{E}_{11}=0.1137$. We
remark that in our numerical experiment, we choose $M=1000$ for 2D random flow
and $M=100$ for 3D random flow so that the velocity field numerically
satisfies the ergodicity assumption.

In Fig.1(a), we plot the convergence results of the effective diffusivity for
the 2D random flow using our method (i.e., $\frac{\mathbb{E}[({\bf
X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$) with respective to different time-step
$\Delta t$ at $T=22$, where the number of the Monte Carlo samples
$N_{mc}=50,000$. In addition, we show a fitted straight line with the slope
$1.17$, i.e., the convergence rate is about $O(\Delta t)^{1.17}$. Similarly,
we show the convergence results of $\frac{\mathbb{E}[({\bf
X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$ for the 3D random flow in Fig.1(b) with
respective to different time-step $\Delta t$ at $T=25$, where the number of
the Monte Carlo samples $N_{mc}=50,000$. We also show a fitted straight line
with the slope $0.98$, i.e., the convergence rate is about $O(\Delta
t)^{0.98}$. These numerical results agree with our error analysis.

(a) 2D random flow, fitted slope $\approx$ 1.17

(b) 3D random flow, fitted slope $\approx$ 0.98.

Figure 1: Error of $D^{E}_{11}$ for random flows with different time-steps.

###  5.3 Comparasion between the volume-preserving scheme and Euler scheme

To demonstrate the benefit of our method in computing effective diffusivity,
we compare the performance of the volume-preserving scheme and Euler-Maruyama
scheme (also called the Euler scheme).

For the random flow in 2D space, we solve the SDE (1), where the velocity
filed is chosen as (80) with the setting (81). The time-mixing constant
$\theta=1$ in the covariance function and other parameters are the same as
that were used in Section 5.2. We use the volume-preserving scheme with
$\Delta t_{ref}=0.003125$ and $N_{mc}=100,000$ to solve the SDE (1) and
compute the reference solution, i.e., the “exact” effective diffusivity, where
the final computational time is $T=54$ so that the calculated effective
diffusivity converges to a constant. It takes about 24 hours to compute the
reference solution on a 64-core server (Gridpoint System at HKU). The
reference solution for the effective diffusivity is $D^{E}_{11}=0.3610$.

In Fig.2(a), we plot the convergence results of the effective diffusivity for
the 2D random flow using the volume-preserving scheme and the Euler scheme
with respect to different time-step $\Delta t$ at $T=54$, where the number of
the Monte Carlo samples $N_{mc}=50,000$. The slopes of the fitted lines for
the volume-preserving scheme and the Euler scheme are $0.86$ and $0.62$,
respectively. In addition, we can see that the volume-preserving scheme
reduces the numerical error by more than one order of magnitude than that of
the Euler scheme by using the same time-step $\Delta t$.

For the random flow in 3D space, we solve the SDE (1), where the velocity
field is chosen as (80) with the setting (83). The time-mixing constant
$\theta=4$ in the covariance function and other parameters are the same as
that were used in Section 5.2. We use the volume-preserving scheme with
$\Delta t_{ref}=0.003125$ and $N_{mc}=100,000$ to solve the SDE (1) and
compute the reference solution, where the final computational time is $T=40$
so that the calculated effective diffusivity converges to a constant. It takes
about 32 hours to compute the reference solution on a 64-core server
(Gridpoint System at HKU). The reference solution for the effective
diffusivity is $D^{E}_{11}=0.2266$.

In Fig.2(b), we plot the convergence results of the effective diffusivity for
the 3D random flow using the volume-preserving scheme and the Euler scheme
with respect to different time-step $\Delta t$ at $T=40$, where the number of
the Monte Carlo samples $N_{mc}=50,000$. The slopes of the fitted lines for
the volume-preserving scheme and the Euler scheme are $0.91$ and $0.44$,
respectively. Again, we can see that the volume-preserving scheme
significantly reduces the numerical error by more than one order of magnitude
than that of the Euler scheme by using the same time-step $\Delta t$.

We remark that the volume-preserving scheme is an implicit scheme which needs
to use Newton’s iteration method to solve the corresponding nonlinear
equations. In our numerical experiments, we use the numerical solutions at
time $t=t_{n}$ as an initial guess for the solution at time $t=t_{n+1}$. We
find this approach is very efficient, i.e., three or four steps of iterations
will give convergent results. Thus, the computational cost for the volume-
preserving scheme is about three or four times of Euler scheme in the same
setting. However, the volume-preserving scheme is superior to the Euler scheme
due to its faster convergence rate and smaller magnitude in the numerical
error.

(a) 2D comparasion: fitted slope for Euler scheme is 0.62, fitted slope for
volume preserving scheme is 0.86.

(b) 3D comparasion: fitted slope for Euler scheme is 0.44, fitted slope for
volume preserving scheme is 0.91.

Figure 2: Comparasion between Euler scheme and volume preserving scheme

###  5.4 Verification of the exponential decay property.

The time relaxation property (4), which is equivalent to the exponential decay
property (5), plays an important role in the existence of the effective
diffusivity; see Prop. 2.5. In Theorem 3.6, we prove that the numerical
solutions inherit the exponential decay property. Based on this key fact, we
can define the discrete-type corrector problem and prove the convergence
analysis of our method. In this subsection, we will verify that the velocity
field propagated by the random flow (80) has the exponential decay property,
where both the 2D and 3D cases will be tested.

In the experiment for 3D random flow, we choose the time step size $\Delta
t=0.05$. The velocity field will be approximated by $M=100$ terms in (80) with
the setting (83). The parameters in the spectral measure $\Gamma(\bf k)$ are
$K=10$ and $\alpha=0.75$. The molecular diffusivity $\sigma=0.1$. We randomly
generate 200 samples $\\{{\bf
k}_{m}^{i},{\bm{\xi}}_{m}^{i}(0),{\bm{\eta}}_{m}^{i}(0),m=1,...,M\\}$,
$i=1,...,200$, which will be used to generate initial states for the velocity
field (80), i.e.,

$\displaystyle{\bf b}^{i}(0,{\bf
x})=\frac{1}{\sqrt{M}}\sum_{m=1}^{M}\big{[}{\bm{\xi}}_{m}^{i}(0)\times\frac{{\bf
k}^{i}_{m}}{|{\bf k}^{i}_{m}|}\cos({\bf k}_{m}^{i}\cdot{\bf
x})+{\bm{\eta}}^{i}_{m}(0)\times\frac{{\bf k}^{i}_{m}}{|{\bf
k}^{i}_{m}|}\sin({\bf k}_{m}^{i}\cdot{\bf x})\big{]},\quad i=1,...,200.$

Then, for each initial state ${\bf b}^{i}(0,{\bf x})$, we generate $5000$
different samples of the OU paths ${\bm{\xi}}^{i,p}_{m}(n\Delta t)$ and
${\bm{\eta}}^{i,p}_{m}(n\Delta t)$ and Brownian motion paths
$\bm{w}^{i,p}(n\Delta t)$, $1\leq p\leq 5000$. Given the sample data, we
calculate the corresponding solution paths $\\{{\bf X}_{n}^{i,p}\\}_{0\leq
n<\infty}$ and then calculate the value

$\displaystyle{\bf b}^{i,p}(n\Delta t,{\bf X}^{i,p}_{n})$
$\displaystyle=\frac{1}{\sqrt{M}}\sum_{m=1}^{M}[{\bm{\xi}}_{m}^{i,p}(n\Delta
t)\times\frac{{\bf k}^{i}_{m}}{|{\bf k}^{i}_{m}|}\cos({\bf k}_{m}^{i}\cdot{\bf
X}^{i,p}_{n})+{\bm{\eta}}^{i}_{m}(n\Delta t)\times\frac{{\bf k}^{i}_{m}}{|{\bf
k}^{i}_{m}|}\sin({\bf k}_{m}^{i}\cdot{\bf X}^{i,p}_{n})],$ $\displaystyle
i=1,...,200,\quad 1\leq p\leq 5000.$ (84)

Finally, we compute $\bar{\bf b}_{n}^{i}=\frac{1}{5000}\sum_{p=1}^{5000}{\bf
b}^{i,p}(n\Delta t,{\bf X}^{i,p}_{n})$ and the sample variance of $\bar{\bf
b}_{n}^{i}$ with respect to $i$. This is an approximation to the value
$||S_{n}{\bf b}||_{L^{2}(\mathcal{X})}$, which should satisfies exponential-
decay property according to our analysis. The experiment for 2D random flow is
almost the same except the setting of the velocity filed (83) is replaced by
(81) and we choose $M=1000$.

In Fig.3(a) and Fig.3(b), we plot the calculated sample variance of the first
component of $\bar{\bf b}_{n}^{i}$ for the 2D random flow and 3D random flow,
respectively. We observe exponential decay of the sample variance with respect
to time. Moreover, we find that larger $\theta$ leads to a faster decay in the
sample variance, since larger $\theta$ results in a faster decorrelation in
the random flow. Our numerical results show that the exponential decay
property (see Theorem 3.6) holds for the random flows we studied here.

(a) Calculated variance in the 2D flow along time.

(b) Calculated variance in the 3D flow along time.

Figure 3: Decay behaviors of the sample variance in 2D and 3D random flows.

###  5.5 Investigation of the convection-enhanced diffusion phenomenon

In the first experiment, we study the relation between the numerical effective
diffusivity $\frac{E[({\bf X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$ and the
parameter $\theta$, which controls the de-correlation rate in the temporal
dimension of the random flow. In this experiment, the setting of the velocity
field and the implementation of our method is the same as we used in Section
5.4. We only choose different parameter $\theta$ to compute the numerical
effective diffusivity.

In Fig.4(a), we plot the numerical effective diffusivity of 2D random flow
obtained at different computational times, where the flow is generated with
different $\theta$. The result for 3D random flow is shown in Fig.4(b). We
find that different $\theta$ affects the mixing time of the system. When we
increase the $\theta$, the system will quickly enter a mixing stage.

(a) The quantity $\frac{E[({\bf X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$ in the
2D flow along time.

(b) The quantity $\frac{E[({\bf X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$ in the
3D flow along time.

Figure 4: The relation between numerical effective diffusivity and $\theta$.

In the second experiment, we choose different molecular diffusivity $\sigma$
to compute the corresponding numerical effective diffusivity, which allows us
to study the existence of residual diffusivity for this random flow. The
residual diffusivity, a special yet remarkable convection-enhanced diffusion
phenomenon, refers to the non-zero and finite effective diffusivity in the
limit of zero molecular diffusivity as a result of a fully chaotic mixing of
the streamlines.

In the experiment for 2D random flow, we choose the time step $\Delta t=0.05$,
the velocity field were simulated with $M=1000$, the time-mixing constant
$\theta=0.1$ and the parameters in the spectral measure $\Gamma(\bf k)$ are
$K=10$ and $\alpha=0.75$. For the 3D random flow, we choose $M=100$ and keep
other parameters the same.

Let $\kappa=\sigma^{2}/2$. In Fig.5(a), we show the relation between numerical
effective diffusivity of 2D random flow obtained at different computational
times, where the result is generated with different $\sigma$. The result for
3D random flow is shown in Fig.5(b). We find that as $\kappa$ approaches zero,
the quantity $\frac{\mathbb{E}[(\bar{\bf X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$
converges to a non-zero (positive) constant, which indicates the existence of
residual diffusivity in the random flows here.

(a) The quantity $\frac{\mathbb{E}[({\bf X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$
in the 2D flow along time.

(b) The quantity $\frac{\mathbb{E}[({\bf X}^{\omega}_{n,1})^{2}]}{2n\Delta t}$
in the 3D flow along time.

Figure 5: The relation between numerical effective diffusivity and
$\kappa=\sigma^{2}/2$, where $\sigma$ is molecular diffusivity.

In Fig.6(a) and Fig.6(b), we plot the convergence behaviors of
$D^{E}_{11}(\kappa)$ approaching $D^{E}_{11}(0)$ for the 2D and 3D random
flows, respectively, when the systems enter a mixing stage. The convergence
behaviors when $\kappa$ approaches zero are slightly different though, both
figures show that residual diffusivity exists in the random flows we studied
here.

(a) Results for the 2D random flow, $D_{0}\approx 0.3584$.

(b) Results for the 3D random flow, $D_{0}\approx 0.5685$.

Figure 6: Convergence behaviors of $D^{E}_{11}(\kappa)$ approaching
$D^{E}_{11}(0)$.

##  6 Conclusion

In this paper, we studied the numerical homogenization of passive tracer
models in random flows. Based on a splitting method, we proposed stochastic
structure-preserving schemes to compute the effective diffusivity of the
random flows. In addition, we provided rigorous convergence analysis for the
numerical schemes. Our error analysis is new in the sense that it is based on
a probabilistic approach. Specifically, we interpreted the solution process
generated by our numerical schemes as a Markov process. By using the ergodic
theory for the solution process, we proved a sharp error estimate for our
numerical schemes in computing the effective diffusivity. Finally, we present
numerical results to verify the convergence rate of the proposed method for
incompressible random flows both in 2D and 3D spaces. In addition, we observed
the exponential decay property and investigated the residual diffusivity
phenomenon in the random flows we studied here.

There are two directions we plan to explore in our future work. First, we
shall extend the probabilistic approach to provide sharp convergence analysis
in computing effective diffusivity for quasi-periodic time-dependent flows.
This type of problem is more challenging since the corrector problem does not
exist in the $L^{2}$ space corresponding to the invariant measure. We shall
develop other techniques to address this problem. In addition, we shall
investigate the convection-enhanced diffusion phenomenon for general spatial-
temporal stochastic flows [20, 22] and develop convergence analysis for the
corresponding numerical methods.

## Acknowledgement

The research of J. Lyu and Z. Wang are partially supported by the Hong Kong
PhD Fellowship Scheme. The research of J. Xin is partially supported by NSF
grants DMS-1211179, DMS-1522383, and IIS-1632935. The research of Z. Zhang is
supported by Hong Kong RGC grants (Projects 27300616, 17300817, and 17300318),
National Natural Science Foundation of China (Project 11601457), Seed Funding
Programme for Basic Research (HKU), and Basic Research Programme
(JCYJ20180307151603959) of The Science, Technology and Innovation Commission
of Shenzhen Municipality. The computations were performed using the HKU ITS
research computing facilities that are supported in part by the Hong Kong UGC
Special Equipment Grant (SEG HKU09).

## References

  * [1] A. Bensoussan, J. L. Lions, and G. Papanicolaou, Asymptotic analysis for periodic structures, vol. 374, American Mathematical Soc., 2011.
  * [2] S. Boi, M. Afonso, and A. Mazzino, Anomalous diffusion of inertial particles in random parallel flows: theory and numerics face to face, Journal of Statistical Mechanics: Theory and Experiment, 2015 (2015), p. P10023.
  * [3] S. Boi, A. Mazzino, and P. Muratore-Ginanneschi, Eddy diffusivities of inertial particles in random gaussian flows, Physical Review Fluids, 2 (2017), p. 014602.
  * [4] R. Carmona and L. Xu, Homogenization for time-dependent two-dimensional incompressible Gaussian flows, The Annals of Applied Probability, 7 (1997), pp. 265–279.
  * [5] Z. Chai, C. Huang, B. Shi, and Z. Guo, A comparative study on the lattice Boltzmann models for predicting effective diffusivity of porous media, International Journal of Heat and Mass Transfer, 98 (2016), pp. 687–696.
  * [6] P. Dedik and M. Subin, Random pseudodifferential operators and the stabilization of solutions of parabolic equations with random coefficients., Sbornik: Mathematics, 41 (1982), pp. 33–52.
  * [7] P. Doukhan, Mixing: properties and examples, vol. 85, Springer Science & Business Media, 2012.
  * [8] W. E and E. V. Eijnden, Generalized flows, intrinsic stochasticity, and turbulent transport, Proceedings of the National Academy of Sciences, 97 (2000), pp. 8200–8205.
  * [9] A. Fannjiang and T. Komorowski, Turbulent diffusion in markovian flows, The Annals of Applied Probability, 9 (1999), pp. 591–610.
  * [10] A. Fannjiang and G. Papanicolaou, Convection-enhanced diffusion for periodic flows, SIAM J Appl. Math., 54 (1994), pp. 333–408.
  * [11] , Diffusion in turbulence, Probability Theory and Related Fields, 105 (1996), pp. 279–334.
  * [12] , Convection-enhanced diffusion for random flows, J. Stat. Phys., 88 (1997), pp. 1033–1076.
  * [13] K. Feng and Z. Shang, Volume-preserving algorithms for source-free dynamical systems, Numerische Mathematik, 71 (1995), pp. 451–463.
  * [14] J. Garnier, Homogenization in a periodic and time-dependent potential, SIAM Journal on Applied Mathematics, 57(1) (1997), pp. 95–111.
  * [15] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, Springer Science and Business Media, 2006.
  * [16] V. V. Jikov, S. Kozlov, and O. A. Oleinik, Homogenization of Differential Operators and Integral Functionals, Springer, Berlin, 1994\.
  * [17] P. E. Kloeden and E. Platen, Numerical solution of Stochastic Differential Equations, Springer-Verlag Berlin, Heidelberg, 1992.
  * [18] T. Komorowski and S. Olla, On homogenization of time-dependent random flows, Probability theory and related fields, 121 (2001), pp. 98–116.
  * [19] R. Kraichnan, Diffusion by a random velocity field, The physics of fluids, 13 (1970), pp. 22–31.
  * [20] C. Landim, S. Olla, and H. Yau, Convection–diffusion equation with space–time ergodic random flow, Probability theory and related fields, 112 (1998), pp. 203–220.
  * [21] A. Loisy, A. Naso, and P. Spelt, The effective diffusivity of ordered and freely evolving bubbly suspensions, Journal of Fluid Mechanics, 840 (2018), pp. 215–237.
  * [22] A. J. Majda and P. R. Kramer, Simplified models for turbulent diffusion: theory, numerical modelling, and physical phenomena, Phys. Rep., 314 (1999), pp. 237–574.
  * [23] R. I. McLachlan and G. R. Quispel, Splitting methods, Acta Numerica, 11 (2001), pp. 341–434.
  * [24] G. N. Milstein, Numerical integration of stochastic differential equations, vol. 313, Springer Science & Business Media, 1994.
  * [25] G. Pavliotis and A. Stuart, Homogenization for inertial particles in a random flow, Commun Math Sci., 5 (2007), pp. 507–531.
  * [26] , Multiscale methods: averaging and homogenization, Springer Science and Business Media, 2008.
  * [27] G. Pavliotis, A. Stuart, and K. Zygalakis, Calculating effective diffusivities in the limit of vanishing molecular diffusion, J. Comput. Phys., 228 (2009), pp. 1030–1055.
  * [28] A. Renaud and J. Vanneste, Dispersion and reaction in random flows: Single realization versus ensemble average, Physical Review Fluids, 4 (2019), p. 124502.
  * [29] R. Roque-Malherbe, Adsorption and diffusion in nanoporous materials, CRC press, 2018.
  * [30] M. Rosenblatt, Markov Processes: Structure and Asymptotic Behavior, vol. 184, Springer Science & Business Media, 2012.
  * [31] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal., 5 (1968), pp. 506–517.
  * [32] G. Taylor, Diffusion by continuous movements, Proceedings of The London Mathematical Society, (1922), pp. 196–212.
  * [33] Z. Wang, J. Xin, and Z. Zhang, Sharp uniform in time error estimate on a stochastic structure-preserving Lagrangian method and computation of effective diffusivity in 3D chaotic flows, arXiv:1808.06309, (2018).
  * [34] Z. J. Wang, J. Xin, and Z. W. Zhang, Computing effective diffusivity of chaotic and stochastic flows using structure-preserving schemes, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2322–2344.