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

TL;DR

This paper introduces stochastic structure-preserving numerical schemes for accurately computing effective diffusivity in particles moving through random flows, with rigorous convergence analysis and numerical validation.

Contribution

The paper develops novel stochastic structure-preserving schemes for SDEs in random flows, providing a rigorous convergence analysis and connecting discrete and continuous corrector problems.

Findings

Schemes accurately compute effective diffusivity in random flows.

Convergence analysis uses probabilistic approach and Markov process interpretation.

Numerical results demonstrate efficiency and reveal convection-enhanced diffusion phenomena.

Abstract

In this paper, we propose stochastic structure-preserving schemes to compute the effective diffusivity for particles moving in random flows. We first introduce the motion of particles using the Lagrangian formulation, which is modeled by stochastic differential equations (SDEs). We also discuss the definition of the corrector problem and effective diffusivity. Then we propose stochastic structure-preserving schemes to solve the SDEs and provide a sharp 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 using the central limit theorem for the solution process, we obtain the convergence analysis of our method in computing long time solutions. Most importantly our convergence analysis reveals the connection of discrete-type and continuous-type corrector problems, which is fundamental and interesting. 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.

Figures

Equations

$$d\mathbf{X}(t)=\mathbf{b}(t,\mathbf{X}(t),\omega )dt+\sigma d\mathbf{w}(t),\mathbf{X}(0)=0,$$

$$U^{\mathbf{x}}f(\chi )=f({\tau }_{\mathbf{x}}(\chi )),f\in L^2(\mathcal{X}),\mathbf{x}\in R^d.$$

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

$$-(Lf,f{)}_{L^2(\mathcal{X})}\ge c_1||f|{|}_{L^2(\mathcal{X})}^2,\text{where}{c}_1>0.$$

$$||P^{t}f|{|}_{L^2(\mathcal{X})}\le \exp (-2c_{1}t)||f|{|}_{L^2(\mathcal{X})},f\in L_0^2(\mathcal{X}).$$

$$\{\begin{array}{rl}& d{\mathbf{X}}_t^{\omega }=\mathbf{b}(t,{\mathbf{X}}_t^{\omega },\omega )dt+\sigma d{\mathbf{w}}_t,\\ & {\mathbf{X}}_0^{\omega }=0,\end{array}$$

$$\{\begin{array}{rl}& \eta (t)={\tau }_{{\mathbf{X}}_t^{\omega }}\omega (t),\\ & \eta (0)=\omega (0).\end{array}$$

$$S^{t}f(\chi )={\mathbf{ME}}_{\chi }f(\eta (t)),t\ge 0,f\in L^{\infty }(\mathcal{X}),$$

$$\int S^{t}fd{P}_0=\int fd{P}_0,t\ge 0,f\in L^2(\mathcal{X}).$$

$$\mathcal{L}f=Lf+\frac{{\sigma }^2}{2}\Delta f+\mathbf{b}\cdot \nabla f,$$

$$||S^{t}f|{|}_{L^2(\mathcal{X})}\le \exp (-2c_{1}t)||f|{|}_{L^2(\mathcal{X})},\text{where}{c}_1>0.$$

$$(-\mathcal{L}f,f{)}_{L_0^2(\mathcal{X})}\ge (-Lf,f{)}_{L_0^2(\mathcal{X})}\ge c_1||f|{|}_{L_0^2(\mathcal{X})}^2$$

$$\frac{d}{dt}||S^{t}f|{|}_{L^2(\mathcal{X})}^2=2(\mathcal{L}{S}^{t}f,S^{t}f{)}_{L^2(\mathcal{X})}\le -2c_1||S^{t}f|{|}_{L^2(\mathcal{X})}^2,$$

$$||S^{t}f|{|}_{L^2(\mathcal{X})}^2\le \exp (-2c_{1}t)||f|{|}_{L^2(\mathcal{X})}^2,\forall t\ge 0,$$

$$\mathbf{\psi }={\int }_0^{\infty }{S}^t\mathbf{b}dt$$

$$\mathcal{L}\mathbf{\psi }=-\mathbf{b}$$

$${\mathbf{v}}^T{D}^E\mathbf{v}=\frac{{\sigma }^2}{2}+(-\mathcal{L}{\psi }_{\mathbf{v}},{\psi }_{\mathbf{v}}{)}_{L^2(\mathcal{X})},$$

$$d{\mathbf{X}}_t^{\omega }$$

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

$$(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{)},$$

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

$${\mathbf{X}}_{(n+1)\Delta t}^{\omega }-{\mathbf{X}}_{n\Delta t}^{\omega }={\int }_{n\Delta t}^{(n+1)\Delta t}\mathbf{b}(t,{\mathbf{X}}_t^{\omega },\omega )dt.$$

$${\mathbf{X}}_{n+1}^{\omega }={\mathbf{X}}_n^{\omega }+{\mathbf{B}}_{\Delta t}^{\omega (n\Delta t)}({\mathbf{X}}_n^{\omega })+\sigma {\mathbf{\xi }}_n,$$

$${\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{)}.$$

$${\mathbf{X}}_{(n+1)\Delta t}^{\omega }={\mathbf{X}}_n^{\omega }+{\int }_{n\Delta t}^{(n+1)\Delta t}\mathbf{b}(t,{\mathbf{X}}_t^{\omega },\omega )dt.$$

$$\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.$$

$$\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{|}$$

$$\le$$

$$\{\begin{array}{rl}& {\eta }_n={\tau }_{{\mathbf{X}}_n^{\omega }}\omega (n\Delta t),\\ & {\eta }_0=\omega (0).\end{array}$$

$$S_{n}f(\chi )={\mathbf{ME}}_{\chi }f({\eta }_n),$$

Full text

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

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,

[TABLE]

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:

$\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

[TABLE]

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.

$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

[TABLE]

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,

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

In addition, environment process generates a semigroup of transformation

[TABLE]

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

Proposition 2.2**.**

[[9], Prop. 3]

$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,*

[TABLE]

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.,

[TABLE]

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

Proposition 2.3**.**

[[9],Prop 4]

$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})$

[TABLE]

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

[TABLE]

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,

[TABLE]

thus

[TABLE]

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

[TABLE]

which satisfies the following continuous-type corrector problem

[TABLE]

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

[TABLE]

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.

[TABLE]

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,

[TABLE]

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)

[TABLE]

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

[TABLE]

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,

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Moreover,

[TABLE]

Proof.

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

[TABLE]

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

[TABLE]

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*

[TABLE]

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

Proof.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Proof.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

We arrive that

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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.

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.,

[TABLE]

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

[TABLE]

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.

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.

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.

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.

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).