COMPUTING EFFECTIVE DIFFUSIVITIES IN 3D TIME-DEPENDENT CHAOTIC FLOWS WITH A CONVERGENT LAGRANGIAN NUMERICAL METHOD

. In this paper, we study the convergence analysis for a robust stochastic structure-preserving Lagrangian numerical scheme in computing effective diffusivity of time-dependent chaotic flows, which are modeled by stochastic differential equations (SDEs). Our numerical scheme is based on a splitting method to solve the corresponding SDEs in which the deterministic subproblem is discretized using a structure-preserving scheme while the random subproblem is discretized using the Euler-Maruyama scheme. We obtain a sharp and uniform-in-time convergence analysis for the proposed numerical scheme that allows us to accurately compute long-time solutions of the SDEs. As such, we can compute the effective diffusivity for time-dependent chaotic flows. Finally, we present numerical results to demonstrate the accuracy and efficiency of the proposed method in computing effective diffusivity for the time-dependent Arnold-Beltrami-Childress (ABC) flow and Kolmogorov flow in three-dimensional space


Introduction
In this paper, we study the convection-enhanced diffusion phenomenon for particles moving in time-dependent chaotic flows, which is defined by the following passive tracer model, i.e., a stochastic differential equation (SDE), dX() = v(, X)d + dw(), X ∈ R  , ( where X() = ( 1 (), . . .,   ())  ∈ R  is the position of the particle,  > 0 is the molecular diffusivity, and {w()} ≥0 is the standard -dimensional Brownian motion.The velocity field v(, x) is time-dependent and divergence free, i.e., ∇ x • v(, x) = 0, for all  ≥ 0. In order to guarantee the existence of the solution to the SDE (1.1), we also assume that v(, x) is Lipschitz in x.The passive tracer model (1.1) has many applications in physical and engineering sciences, including atmosphere science, ocean science, chemical engineering and combustion [23].
We will study the long-time large-scale behavior of the particle X() in the passive tracer model (1.1), i.e., whether the motion of the particle X() has a long-time diffusive limit.Let X  () ≡ X(/ 2 ) denote the rescaled process of (1.1).We aim to investigate whether X  () converge in law to a Brownian motion with a covariance matrix   ∈ R × as  → 0, where   is called the effective diffusivity matrix.The dependence of   on the velocity field of the passive tracer model is complicated and highly nontrivial.There are many theoretical works, where the homogenization theory was applied to study the effective diffusivity matrix   of the passive tracer model with spatial periodic velocity fields or random velocity fields with short-range correlations; see e.g., [3,12,16,28] and references therein.
For many complicated velocity fields of physical interests, one cannot apply the homogenization theory to compute the corresponding effective diffusivity matrix   , or even determine its existence.Therefore, many numerical methods were developed to compute   .These results include, among others, for time-independent Taylor-Green flows, the authors of [29] proposed a stochastic splitting method and calculated effective diffusivity in the limit of vanishing molecular diffusivity.For time-dependent chaotic flows, an efficient model reduction method based on the spectral method was developed to compute   using the Eulerian framework [22].The reader is referred to [23] for an extensive review of many existing mathematical theories and numerical simulations for the passive tracer model with different velocity fields.
Recently, we developed a robust structure-preserving Lagrangian scheme to compute the effective diffusivity for chaotic and stochastic flows in [32].We also obtained a rigorous error estimate for the numerical scheme in [32].Specifically, let  ,num denote the numerical effective diffusivity obtained by our method.We got the error estimate, || ,num −   || ≤  + ( )() 2 , where the computational time  should be greater than the diffusion time (also known as mixing time).This error estimate is not sharp in the sense that the pre-factor ( ) may grow fast with respect to  , since the error estimation is based on a Gronwall inequality technique.Later, we obtained a sharp convergence rate for our numerical scheme and got rid of the term ( ) in the error estimate in [33].However, this technique can only be used to study passive tracer models in time-independent (steady) flows, and cannot be applied to study passive tracer models in time-dependent (unsteady) flows.
In this paper, we aim to obtain a sharp convergence analysis for our numerical scheme in computing effective diffusivity of passive tracer models in spatial-temporal periodic velocity fields.These types of flow fields are wellknown for exhibiting chaotic streamlines and have many applications in turbulent diffusion [23].Since in this case the velocity field depends on the temporal variable, the generator associated with the stochastic process, i.e., the solution X() in equation (1.1) becomes non-autonomous.The generator is now a parabolic-type operator (see Eq. (2.2)), instead of an elliptic-type operator studied in [33] when the flows are time-independent.The cell problem is then defined in a space-time periodic domain; see equation (2.4).Hence the extension of the analysis developed in [33] to time-dependent flows is not straightforward.We will develop new techniques to overcome the difficulty arising from time dependence; see Theorem 4.3 and Lemma 4.5 in Section 4. We also emphasize that when the flows are time-independent, one can construct ballistic orbits of the ABC and Kolmogorov flows [17,24,34] and study their dynamic behaviors.When the flows are time-dependent however, their more complicated streamlines make it challenging to construct and study ballistic orbits even if they persist.
Though there are several prior works on structure-preserving schemes for ODEs and SDEs, see e.g., [1,14,15,21,31] and references therein, our work has several novel contributions.The first novelty is the convergence analysis, where we develop new techniques to deal with time-dependent flows.To handle the parabolic-type generator, we pile up snapshots of each time step within a single time period together.By viewing the numerical solutions as a Markov process and exploring the ergodicity of the solution process, we succeed in obtaining a sharp convergence analysis for our method in computing the effective diffusivity, where the error estimate does not depend on the computational time.Therefore, we can compute the long-time solutions of passive tracer models without losing accuracy; see Figures 1 and 2. If we choose the Gronwall inequality in the error estimate, we cannot get rid of the exponential growth pre-factor in the error term, which makes the convergence analysis not sharp.Most importantly, our convergence result reveals the equivalence of the definition of effective diffusivity using the Eulerian framework and the Lagrangian framework; see Theorem 4.10, which is fundamental and important.For 3D time-dependent flows, the Eulerian framework has good theoretical value yet the Lagrangian framework is mesh-free and computationally more accessible.
Another novelty is that the stochastic structure-preserving Lagrangian scheme is robust and quite cheap in computing the long-time solutions of the passive tracer model (1.1), especially for problems in three-dimensional space.If one adopts the Eulerian framework to compute the effective diffusivity of the passive tracer model (1.1), one needs to solve a convection-diffusion-type cell problem; see equation (2.4).When the molecular diffusivity  is small and/or the dimension of spatial variables is high, say  = 3, it is exorbitantly expensive to solve the cell problems.As indicated in equation (2.5), the effective diffusivities depend on the integration of the gradient of the solution to the cell problem.In many cases, e.g., time-dependent ABC flow, the effective diffusivities grow rapidly as  decreases; see Figure 3.In our Lagrangian approach, we can overcome the difficulties of long-time integration of the SDEs (raised as  decreases) by using robust structure-preserving schemes.However, for the Eulerian approach, one needs to solve a four-dimensional PDE (three variables in spatial dimension and one variable in the temporal dimension) and solutions have sharp gradients as the diffusivity decreases, which makes the Eulerian approach expensive for computing effective diffusivities.
Numerical results show that our Lagrangian scheme is insensitive to the molecular diffusivity  and computational cost linearly depends on the dimension  of spatial variables in the passive tracer models (1.1).Thus, we are able to investigate the convection-enhanced diffusion phenomenon for several typical time-dependent chaotic flows of physical interests, including the time-dependent ABC flow and the time-dependent Kolmogorov flow in three-dimensional space.We discover that the maximal enhancement is achieved in the former case, while a submaximal enhancement is observed in the latter case; see Figures 3 and 4b, respectively.In addition, we find that the level of chaos and the strength of diffusion enhancement seem to compete with each other in the time-dependent ABC flow; see Figure 5.To the best of our knowledge, our work is the first in the literature to develop a convergent Lagrangian method to study convection-enhanced diffusion phenomenon in 3D time-dependent chaotic flows.
The rest of the paper is organized as follows.In Section 2, we give the definition of the effective diffusivity matrix using the Eulerian framework and the Lagrangian framework.In Section 3, we propose the stochastic structure-preserving Lagrangian scheme in computing effective diffusivity for the passive tracer model (1.1).In Section 4, we provide a sharp convergence analysis for the proposed method based on a probabilistic approach.In addition, we shall show that our method can be used to solve high-dimensional flow problems and the error estimate can be obtained naturally.In Section 5, we present numerical results to demonstrate the accuracy and efficiency of our method.We also investigate the convection-enhanced diffusion phenomenon for time-dependent chaotic flows.Concluding remarks are made in Section 6.

Effective diffusivity of the passive tracer models
There are two frameworks to define the effective diffusivity of the passive tracer models.We first discuss the Eulerian framework.One natural way to study the expectation of the paths for the SDE given by equation (1.1) is to consider its associated backward Kolmogorov equation [27].Due to the time-dependence nature of the velocity field, we need to deal with a space-time ergodic random flow.Specifically, given a sufficiently smooth function and X() be the solution to equation (1.1).Then, (, , x) satisfies the following backward Kolmogorov equation (2.1) In equation (2.1), the generator ℒ is defined as where  0 =  2 /2 is the diffusion coefficient, v is the velocity field, and ∇ x and  x denote the gradient operator and Laplace operator, respectively.When v is incompressible (i.e., ∇ x • v(, x) = 0, ∀), deterministic and space-time periodic in (1) scale, where we assume the period of v is 1 in each dimension of the physical and temporal space, the formula for the effective diffusivity matrix is [3,28] where we have assumed that the fluid velocity v(, x) is smooth and the (vector) corrector field  satisfies the cell problem, and ⟨•⟩  denotes temporal and spatial average over T × T  .Since v is incompressible, the solution  to the cell problem (2.4) is unique up to an additive constant by the Fredholm alternative.By multiplying  to equation (2.4), integrating the corresponding result in T × T  and using the periodic conditions of  and v, we get an equivalent formula for the effective diffusivity as follows: (2.5) The correction to  0  in equation (2.5) is nonnegative definite.We can see that e    e ≥  0 for all unit column vectors e ∈ R  , which is called convection-enhanced diffusion.By using a variational principle for time-periodic velocity flows, one can find a upper bound for the effective diffusivity, i.e., there exists a nonzero unit column vector e ∈ R  , such that e    e ∼ 1  0 , as  0 → 0, (2.6) which is known as the maximal enhancement.More details of the derivation can be found in [4,9,25].We point out that many theoretical results were built upon the passive tracer models with time-independent flows.We are interested in studying the convection-enhanced diffusion phenomenon for time-dependent chaotic flows in this paper.Especially, whether the time-dependent chaotic flows still have the maximal enhancement.
In practice, the cell problem (2.4) can be solved using numerical methods, such as finite element methods and spectral methods.However, when  0 becomes small, the solutions of the cell problem (2.4) develop sharp gradients and demand a large number of finite element basis or Fourier basis functions to resolve, which makes the Eulerian framework expensive.In addition, when the dimension of spatial variables is high, say  = 3, the Eulerian framework becomes expensive too.
Alternatively, one can use the Lagrangian framework to compute the effective diffusivity matrix, which is defined as follows: where X() = ( 1 (), . . .,   ())  is the position of a particle tracer at time  and the average ⟨•⟩ is taken over an ensemble of particles.If the above limit exists, which means the transport of particles is a standard diffusion process, at least on a long-time scale.For example, when the velocity field is the Taylor-Green velocity field [9,29], the long-time and large-scale behavior of the passive tracer model is a diffusion process.However, there are cases showing that the spreading of particles does not grow linearly with time but has a power-law   , where  > 1 and  < 1 correspond to super-diffusive and sub-diffusive behaviors, respectively; see e.g., [2,4,23].
We shall adopt the Lagrangian framework in this paper.The Lagrangian framework has the advantages that: (1) it is easy to implement; (2) it does not directly suffer from a small molecular diffusion coefficient  during the computation; and (3) its computational cost only scales linearly with the dimension of spatial variables in the passive tracer models.However, the major difficulty in solving equation (1.1) is that the computational time should be long enough to approach the diffusion time scale.To address this challenge, we shall develop robust numerical schemes, which are structure-preserving and accurate for long-time integration.Moreover, we aim to develop the convergence analysis of the proposed numerical schemes.Finally, we shall investigate the relationship between parameters of the time-dependent chaotic flows and the corresponding effective diffusivity.
In [32], we proposed a stochastic structure-preserving scheme based on a Lie-Trotter splitting scheme to solve the SDE (3.1).Specifically, we split the problem (3.1) into a deterministic subproblem, which is solved by using a symplectic-preserving scheme (e.g., the symplectic Euler scheme for deterministic equations with frozen time), and a stochastic subproblem, which is solved by using the Euler-Maruyama scheme [27].When  is a constant in (3.6), the Euler-Maruyama scheme exactly solves equation (3.6)Now we discuss how to discretize equation (3.1).From time  =   to time  =  +1 , where  +1 =   + ,  0 = 0, and  is the time step, we assume the numerical solution X  = (  1 ,   2 )  is given, which approximates the exact solution X(  ) to the SDE (3.1) at time   = .Then, we apply the Lie-Trotter splitting method to solve the SDE (3.1) and obtain, where √  2 , and  1 ,  2 ∼  (0, 1) are i.i.d.normal random variables.In this paper, we view the solution sequence X  = (  1 ,   2 )  ,  = 1, 2, 3, . .., generated by the scheme (3.7) as a discrete Markov stochastic process, which enables us to use techniques from stochastic process to obtain a sharp convergence analysis for the numerical solutions; see Section 4.
In a 2D Hamiltonian system, when the system contains an additive temporal noise, for each path of the strong solution of SDE (3.1), the additive noise itself is considered to be a symplectic transform [26].Therefore, we state that the scheme (3.7) is stochastic symplectic-preserving since it preserves symplecticity.Specifically, the scheme (3.7) can be viewed as a composition of two symplectic transforms.In addition, we know that the numerical solution converges to the exact one as the time step  approaches zero.In high-dimensional systems, a structure-preserving scheme refers to a volume-preserving scheme; see Section 4.5.

The backward Kolmogorov equation and related results
We first define the backward Kolmogorov equation associated with equation (3.1) as where the generator ℒ associated with the Markov process in equation (3.1) is given by Recall that the solution (, , x) to equation (3.8) satisfies, (, , where X  is the solution to equation (3.1) and  is a smooth function in R 1 × R 2 .In other words, (, , x) is the flow generated by the original SDE (3.1).
Similarly, we can study the flow generated by the stochastic structure-preserving scheme (3.7).According to the splitting method used in the derivation of the scheme in Section 3.1, we respectively define (3.10) Then,  5 (, •, •) will be the flow at time  =  generated by our stochastic structure-preserving scheme (3.7) and it approximates the solution (, •, •) to equation (3.8) well when  is small.It is also worth mentioning that,  3 (, •, •) is the exact flow generated by the deterministic symplectic Euler scheme in solving equation (3.5).We repeat this process to compute the flow equations of our scheme at other time steps, which approximate the solution (, •, •),  = 2, 3, . . . to equation (3.8) at different time steps.
Remark 3.1.Given the operators ℒ  ,  = 1, 2, 3, 4, there are many possible choices in setting the coefficients for each operator ℒ  and designing the splitting method; see Section 2.5 of [14].Equation (3.10) is a simple choice that was used in this paper.
To analyze the error between the flow operator in equation (3.8) and the composition of operators in equation (3.10), we shall resort to the Baker-Campbell-Hausdorff (BCH) formula, which is widely used in non-commutative algebra [13].For example, in the matrix theory, exp() exp() = exp where  is a scalar,  and  are two square matrices of the same size, [, ] is the Lie-Bracket, and the remaining terms on the right hand side are all nested Lie-brackets.
In our analysis, we replace the matrices in equation (3.11) by differential operators and the BCH formula yields critical insights into the particular structure of the splitting error.Let   denote the composite flow operator associated with equation (3.10), i.e., After propagating for time  = , the exact solution to equation (3.8) started at any  can be represented as Therefore, we can apply the BCH formula to analyze the error between the original flow and the approximated flow.Moreover, we find that computing the th order modified equation associated with equation (3.1) in the backward error analysis (BEA) [7,30] is equivalent to computing the terms of BCH formula up to order ()  in equation (3.12).To show that the solution generated by equation (3.7) follows a perturbed Hamiltonian system (with divergence-free velocity and additive noise) at any order , we only need to consider the ( + 1)nested Lie bracket consisting of }︀ and we can easily see that they generate divergence-free fields.
We remark that given any explicit splitting scheme for deterministic systems, by adding additive noise we shall obtain a similar form of flow propagation.And we shall see in later proof that, the representation of flow operator in equation (3.12) is very effective in analyzing the order of convergence and volume-preserving property.

Convergence analysis
In this section, we prove the convergence rate of our stochastic structure-preserving schemes in computing effective diffusivity based on a probabilistic approach, which allows us to get rid of the exponential growth factor in the error estimate.We first limit our analysis to 2D separable Hamiltonian velocity fields.Then, in Section 4.5 we will show that all the derivations can be generalized to high-dimensional cases.

Convergence to an invariant measure
To compute the effective diffusivity of a passive tracer model using a Lagrangian numerical scheme is closely related to study the limit of a solution sequence (a stochastic process) generated by the numerical scheme.Therefore, we can apply the results from ergodic theory to study the convergence behaviors of the solution.
Let (, ) be a probability space, on which a family  (x, ), x ∈ ,  ∈ , of probability measure is defined.We assume x →  (x, ) is measurable, ∀ ∈ .This corresponds to a linear bounded operator on ℬ(), which is the space of bounded measurable functions on .This operator, denoted by  , is defined by, Clearly || || ≤ 1.One of the main objectives of ergodic theory is to study the limit of the operator sequence   as  → +∞.The result can be summarized into the following proposition, which plays a fundamental role in our convergence analysis.
Then, there exists one and only one invariant probability measure  on (, ) and one has, where  = log denote the transform of the density on Ỹ during [, 1 +  ] (time period is 1) using the numerical scheme (3.7).In addition, let  ,1+ denote the adjoint operator (i.e., the flow operator) of  * ,1+ in the space of ℬ( Ỹ ), which is the set of bounded measurable functions on Ỹ .Then, there exists one and only one invariant probability measure on ( Ỹ , ), denoted by   , satisfying, where  > 0,  > 0 are independent of (•).Moreover, the kernel space of (  −  ,1+ ) is the constant functions in Ỹ , where   is the identity operator.
Proof.We shall verify that the transition kernel associated with the numerical scheme (3.7) satisfies the assumptions required by Proposition 4.1.First we know that in the space R 2 , the integration process associated with the numerical scheme (3.7) can be expressed as a Markov process with the transition kernel, where 2 )  are the numerical solutions at time  =   and  =  +1 , respectively.
Then, using the periodicity of v, we directly extend equation (4.4) to the torus space Ỹ as Let K,+ denote the transition kernel obtained by our scheme with starting time  for  steps.Then, we have We choose  = 1  and obtain K,+1 .One can see that the kernel K,+1 is essentially bounded above zero since K+ in (4.6) are all positive.Moreover, if 0 <  ≪ 1, K,+1 is a continuous function on the domain Ỹ × Ỹ .Then by noticing that the domain Ỹ × Ỹ is compact, the kernel K,+1 is strictly positive.Namely, there exists   > 0 such that K,+1 (X 0 , X  ) >   , ∀(X 0 , X  ) ∈ Ỹ × Ỹ .If we apply Proposition 4.1 to  ,1+ (whose kernel is K,+1 ), we prove the statement in (4.3).
Finally, we know that the operator  ,1+ is compact since it is an integral operator with a continuous kernel.By using the Fredholm alternative, we know that dim ker(  −  ,1+ ) = dim ker(  −  * ,1+ ) = 1.Therefore, it is easy to verify that the constant functions are in the kernel of (  −  ,1+ ).
Equipped with the Lemma 4.2, we study the convergence rate of the space-time transition kernel associated with our numerical scheme (3.7).
,  is a positive integer.The following properties hold: (a) Given , there exists  > 0 and  > 0, such that, where  and  do not depend on  and  .To prove the property (a), we need to show that the lower bound of the kernel K,+1 , which is defined in the proof of Lemma 4.2, does not depend on  .For all  ∈ T, 2 ⌋, where ⌊⌋ denotes the largest integer not greater than .From equation (4.5), we can see that According to the definition of the kernel K,+1 ; see equation (4.6), we know the minimal value of K,+1 is above zero and is independent of  .Now, we apply this observation to Lemma 4.2 and conclude the proof of the property (a).The property (b) is a simple conclusion of the exponential decay property proved in (a).For the property (c), we consider the equation     = .Then, for a given time  , we have  ,1+ (, •) = (, •).The results in Lemma 4.2 imply that the invariant space of  ,1+ is constant in the spatial variable.Thus, we obtain  = ( ).
Before we close this subsection, we provide a convergence result for the inverse of operator sequences, which will be useful in our convergence analysis.  ,  = 1, 2, . . .uniquely exist, then we have a convergence estimate as follows: The proof is quite standard.It can also be viewed as a modification of Theorem 1.16 in Section IV of [18].

A discrete cell problem
In the Eulerian framework, the periodic solution of the cell problem (2.4) and the corresponding formula for the effective diffusivity (2.3) play a key role in studying the behaviors of chaotic and stochastic flows.In the Lagrangian framework, we shall define a discrete analogue of cell problem that enables us to compute the effective diffusivity.Let X 0 = ( 0 1 ,  0 2 )  be the initial data and X  = (  1 ,   2 )  denote the numerical solution at   =  that is generated by the scheme (3.7).
First of all, we show that the solutions   1 and   2 obtained by the scheme (3.7) have bounded expectations if the initial values are bounded.Taking expectation of the first equation of equation (3.7) on both sides, we obtain As a symplectic scheme in 2D, the numerical scheme (3.7) admits the uniform measure as its invariant measure.
Applying the results (a) and (b) of Theorem 4.3 and using the fact that v is a periodic function with zero mean, we know that, sup Here  1 ( + 1 2 , X  ) is equivalent to  1 ( + 1 2 ,   2 ), since  1 is independent of   1 .By applying triangle inequalities in equation (4.11) and using the result in equation (4.12), we arrive at, where  1 does not depend on .Using the same approach, we know that expectation of the second component E  2 is also bounded.Now, we are in the position to define the discrete cell problem.Starting at time  with time step  = 1  , we denote the starting time index to be   .Then, we define where the summation is well defined due to the fact stated in equation (4.12).We will show that v1, (, x) satisfies the following properties.Namely, v1, (, x) is the solution of the discrete cell problem defined in equation (4.15).
Proof.Throughout the proof, we shall use the fact that if ,  are random processes and  is measurable under a filtration ℱ, then with appropriate integrability assumption, we have Some simple calculations will give that v1, (, x) which indicates that v1, has the same regularity as  1 does.We know that the kernel K,+ (x, y) has a fast decay property, which guarantees the order of the differentiation and summation is interchangeable.
Remark 4.6. 1 and v1, only depend on the second component of the numerical solution X  = (  1 ,   2 )  .However, we will write  1 and v1, as functions of X  when we view X  as a Markov process in the convergence analysis.
Remark 4.7.When the flow is time-independent, we obtain Therefore, the discrete cell problem defined in (4.15) is a generalization of the discrete cell problem for timeindependent flow problems studied in our previous work [33], although technically it is more involved.

Convergence estimate of the discrete cell problem
In this section, we shall prove that the solution v1, (, x) of the discrete cell problem (i.e., Eq. (4.15)) converges to the solution of a continuous cell problem in certain subspace.Here, we choose the space C 2,6, 0 (T 1 × Ỹ ) to carry out our analysis.However, there is no requirement that we have to choose this one.In fact, any space that has certain regularity (belongs to the domain of the operator ℒ) will work.The continuous cell problem (2.4) is defined for a vector function, whose first component satisfies For the two-dimensional problem, the operator ℒ is defined in equation (3.9).Given the fact that  1 is a smooth function defined on T 1 × Ỹ , which satisfies ∫︀ Ỹ  1 (, x)dx = 0, ∀ ∈ T 1 .Then, equation (4.21) admits a unique solution  1 in C 2,6, 0 (T 1 × Ỹ ).This is a standard result of parabolic PDEs in Hölder space (see, e.g., the Thm.8.7.3 in [19]).The following theorem states that under certain conditions the solution of the discrete cell problem converges to the solution of the continuous one.
where v1 =  1 + ().Combining equations (4.18) and (4.23), we obtain Equation (4.24) shows the connection between  1 and v1 .After some simple calculations, we get that For the operator L2 , noticing that ℒ = ℒ 1 + ℒ 2 + ℒ 3 + ℒ 4 and operator   is defined in (3.12), we can use the BCH formula and obtain as  approaches zero.Then, applying the Proposition 4.4, we get, In addition, combining the results of equations (4.23), (4.27)-(4.29)for the right hand side of equation (4.25), we know that when  is small enough, the assertion in (4.22) is proved.The constant in the () of (4.22) does not depend on the total computational time  , but may depend on the regularities of  1 ,  2 and the constant .

Convergence analysis for the effective diffusivity
This section contains the main results of our convergence analysis.We first prove that the second-order moment of the solution obtained by using our numerical scheme has an (at most) linear growth rate.Secondly, we provide the convergence rate of our numerical method in computing the effective diffusivity.
Theorem 4.9.Let X  = (  1 ,   2 )  denote the solution of the two-dimensional passive tracer model (3.1) obtained by using our numerical scheme (3.7) with time step .If the Hamiltonian function (,  1 ,  2 ) is separable, periodic and smooth (in order to guarantee the existence and uniqueness of the solution to the SDE (3.1)), then we can prove that the second-order moment of the solution X  (which can be viewed as a discrete Markov process) grows at most linearly, i.e., Proof.We first estimate the second-order moment of the first component of X  = (  1 ,   2 )  , since the other one can be estimated in the same manner.Simple calculations show that The term )] corresponds to the strength of the convection-enhanced diffusion.Our goal here is to prove that it is bounded over , though it may depend on  1 ,  2 and .Noticing that we are calculating the expectation of (  1 ) 2 , which is not defined in the torus space, however in the following derivation we will show that it can be decomposed into sums of periodic functions acting on X  = (  1 ,   2 )  .Hence after the decomposition (see Eq. (4.34)) we can still apply the previous analysis on torus space.
We now directly compute the contribution of the term E[ −1 )] to the effective diffusivity with the help of equation (4.17), Let ℱ  denote the filtration generated by the solution process until X  , for example,   1 ∈ ℱ  .For equation (4.32), we have Hence, we obtain the following result By using the Cauchy-Schwarz inequality, we know the term ]︀ is also bounded.Thus, the assertion in equation (4.30) is proved.
In practice, we shall first choose a time step  and run our numerical scheme (3.7) to compute the effective diffusivity until the result converges to a constant, which may depend on .As such, we shall prove that the limit of the constant converges to the exact effective diffusivity of the original passive tracer model as  approaches zero.Namely, we shall prove that our numerical scheme is robust in computing effective diffusivity.More details on the numerical results will be given in Section 5.
Theorem 4.10.Let   1 ,  = 0, 1, . . .be the first component of the numerical solution obtained by using the scheme (3.7) and  be the time step.We have the convergence estimate of the effective diffusivity as where the constant in () may depend on the regularity of  1 ,  2 and the constant , but does not depend on the computational time  .
Proof.We will prove the statement by direct computation.We divide both sides of equation (4.34) by  and obtain First, we notice that for a fixed , the terms 1  ( 0 1 ) 2 and 2  v1 ( 0 , X 0 ) 0 1 converge to zero as  → ∞, where we have used the fact that v1 ( 0 , X 0 ) is bounded.Also we observe that the term 2 ) is (), because the term ( 1 ) 2 is bounded.Then, for a fixed , we have where the term E[ ] is bounded due to the Theorem 4.9 and ||v 1 || ∞ → || 1 || ∞ < ∞ due to the Theorem 4.8.Therefore, we only need to focus on the estimate of terms in the second line of equation (4.37), which corresponds to the convection-enhanced diffusion effect.Noticing that v1 ∈ C 2,6, , we compute the Ito-Taylor series approximation of v1 (  , X  ), where we have used the fact that  2 (︀ ) is uniformly bounded when  is small enough.Substituting the Taylor expansion of v1 (  , X  ) in equation ( 4.39) into the target term of our estimate (i.e., terms in the second line of Eq. (4.37)), we get Combining the terms with the same order of , we obtain where we have used the facts that: (1) X −1 is independent of  −1 ) 2 = .Finally, by using the Theorem 4.3 and noticing the invariant measure is the uniform measure, we obtain from equation (4.37) that lim Thus, our statement in equation (4.36) is proved using the facts that v1 converges to  1 (see Thm. 4.8) and ∫︀ v1,1 = 0.
Remark 4.11.If we divide two on both sides of equation ( 4.36), we can find that our result recovers the definition of the effective diffusivity   11 defined in equation (2.3).Recall that  0 =  2 /2.Theorem 4.10 reveals the connection of the definition of the effective diffusivity using the Eulerian framework and Lagrangian framework; see equations (2.3) and (2.7), which is fundamental in this context.For 3D time-dependent flow problems, the Eulerian framework has good theoretical values but the Lagrangian framework is computationally accessible.

Generalizations to high-dimensional cases
To show the essential idea of our probabilistic approach in proving the convergence rate of the numerical schemes, we have carried out our convergence analysis based on a two-dimensional model problem (3.1).In fact, the extension of our approach to higher-dimensional problems is straightforward.Now we consider a highdimensional problem as follow: where X = ( 1 ,  2 , . . .,   )  ∈ R  is the position of a particle, v = ( 1 ,  2 , . . .,   )  ∈ R  is the Eulerian velocity field at position X,  is a  ×  constant non-singular matrix, and w() is a -dimension Brownian motion vector.In particular, we assume the component   does not depend on   ,  = 1, . . ., .Thus, the incompressible condition for v(, X) (i.e.∇ X • v(, X) = 0) is easily guaranteed.
For a deterministic and divergence-free dynamical system, Feng et.al. proposed a volume-preserving method [10], which splits a -dimensional problem into  − 1 subproblems with each of them being a two-dimensional problem and thus being volume-preserving.We shall modify Feng's method (first-order case) by including the randomness as the last subproblem to take into account the additive noise, i.e., where X * = ( * 1 ,  * 2 , . . .,  *  )  , W  − W −1 is a -dimensional independent random vector with each component of the form √   ,   ∼  (0, 1), and X  = (  1 ,   2 , . . .,    )  is the numerical approximation to the exact solution X(  ) to the SDE (4.43) at time   = .
The techniques of the convergence analysis for the two-dimensional problem can be applied to highdimensional problems without much difficulty.For the high-dimensional problem (4.43), the smoothness and strict positivity of the transition kernel in the discrete process can be guaranteed if one assumes that the covariance matrix  is non-singular and the scheme (4.44) is explicit.According to our assumption for the velocity field, the scheme (4.44) is volume-preserving for each step.Thus, the solution to the first-order modified equation is divergence-free and the invariant measure on the torus (defined by R  /Z  , when the period is 1) remains uniform for all .Finally, the convergence of the cell problem can be studied by using the BCH formula (3.11) with  + 2 differential operators.Recall that in equation (3.12) we have four differential operators when we study the two-dimensional problem.
Therefore, our numerical methods are robust in computing effective diffusivity for high-dimensional problems, which will be demonstrated through time-dependent chaotic flow problems in three-dimensional space in Section 5.

Numerical results
In this section, we will present numerical examples to verify the convergence analysis of the proposed method in computing effective diffusivity for time-dependent chaotic flows.In addition, we will investigate the convectionenhanced diffusion phenomenon in 3D time-dependent flow, i.e., the time-dependent ABC flow and the timedependent Kolmogorov flow.Without loss of generality, we compute the quantity E(1( )) 2

2𝑇
, which is used to approximate    11 in the effective diffusivity matrix (2.3).

Verification of the convergence rate
We first consider a two-dimensional passive tracer model.Let ( 1 ,  2 )  ∈ R 2 denote the position of a particle.Its motion is described by the following SDE, where  = √ 2 × 0.1,  , ,  = 1, 2 are independent Brownian motions, and the initial data ( 0 1 ,  0 2 )  follows uniform distributions in [−0.5, 0.5] 2 .One can easily verify the velocity field in (5.1) is time-dependent and divergence-free.
In our numerical experiments, we use Monte Carlo samples to discretize the Brownian motions  1, and  2, .The sample number is denoted by   .We choose   = 1 2 12 and   = 3 200 000 to solve the SDE (5.1) to compute the reference solution, i.e., the "exact" effective diffusivity, where the final computational time is  = 3000 to guarantee the convergence of the calculated effective diffusivity along time.In fact, we find that the passive tracer model will enter a mixing stage if the computational time is bigger than  = 1000.It takes about 17 hours to compute the reference solution on a 80-core server (HPC2015 System at HKU).The reference solution for the effective diffusivity is   11 ≈ 0.219.In Figure 6a, we plot the convergence results of the effective diffusivity using our method (i.e., E(1( )) 2

2𝑇
) with respective to different time-step  at  = 3000.To minimize error involved in Monte Carlo simulation, the particle number is the same as in computation of reference solution.In addition, we show a fitted straight line with the slope 1.04, i.e., the convergence rate is about () 1.04 .This numerical result verifies the convergence analysis in Theorem 4.10.
To further study the accuracy and robustness of our method for long-time integration, we consider a 3D time-dependent Kolmogorov flow problem.Let ( 1 ,  2 ,  3 )  ∈ R 3 denote the position of a particle.The motion of a particle moving in the 3D time-dependent Kolmogorov flow is described by the following SDE, where  1, ,  2, and  3, are independent Brownian motions.When  = 0, the velocity field in (5.2) corresponds to the Kolmogorov flow [11].The Kolmogorov flow possesses very chaotic behaviors [6], which brings challenges for our method.
In our numerical experiment, we choose  = 10 −1 and  = √ 2 × 10 −3 in equation (5.2).We choose  ref = 1   2 11   and   = 480 000 to compute the reference solution for the SDE (5.2), i.e., the "exact" effective diffusivity.In our numerical tests, we find that the passive tracer model will enter a mixing stage if the computational time is bigger than  = 2000.To show the accuracy and robustness of our numerical scheme, we set  = 10 5 here.It takes about 59 hours to compute the reference solution on the server and the reference solution for the effective diffusivity is   11 ≈ 0.693.In Figure 6b, we plot the convergence results of the effective diffusivity using our method with respect to different time-step .To minimize error involved in Monte Carlo simulation, the particle number is the same as in computation of reference solution.In addition, we show a fitted straight line with the slope 1.22, i.e., the convergence rate is about () 1.22 .This numerical result again agrees with our error analysis.

Investigation of the convection-enhanced diffusion phenomenon
As we have already demonstrated in Section 5.1, our method is very accurate and robust for long-time integration.Here, we will study the dependence of the effective diffusivity   11 on different parameters in the time-dependent flows.First of all, we solve equation (5.2) and carry out the test for the 3D time-dependent Kolmogorov flow.In Figure 1, we show the time evolution of E(1()) 2

2𝑡
for different  0 's (here  0 =  2 /2) and for four different 's, where the result in Figure 1d corresponds to the time-independent Kolmogorov flow (see Fig. 6 of [33]).The parameter  in equation (5.2) controls the strength of the time dependence.For each  0 and , we use   = 240 000 particles to solve the SDE (5.2).
In Figure 4a, we show the time evolution of [ 1 () 2 ]/(2) for different 's with  0 = 10 −5 .One can see that the effective diffusivity   11 converges as  approaches zero.Similar convergence behaviors were observed for other  0 's, which are not shown here.The convergence of the effective diffusivity with respect to  can be rigorously justified through analysis; see Appendix A.
In addition, in Figure 4b by fixing , we observe a certain amount of enhanced diffusion when  0 decreases.We find that for each given  0 as  decreases the corresponding effective diffusivity   11 converges to the effective diffusivity   11 associated with  = 0.This means the time dependency of  improves the chaotic property of Kolmorogov flow though, it does not change the pattern of convection-enhanced diffusion in the Kolmogorov flow.When  ≤ 1 the fitted slope within  0 ∈ [10 −5 , 10 −3 ] is −0.2, which indicates that   11 ∼ (1/ 0.2 0 ).However, the dependency of   11 on  0 is quite different from the pattern of the time-dependent ABC flow, which is known as the maximal enhancement and will be discussed later; see Figure 3.We call this behavior a sub-maximal enhancement, which may be explained by the fact that the Kolmogorov flow is more chaotic than the ABC flow [11].The chaotic trajectories in Kolmogorov flow enhance diffusion much less than channel-like structures such as the ballistic orbits of ABC flows [24,34].
Next, we use our stochastic structure-preserving scheme to solve time-dependent ABC flow problems.Let ( 1 ,  2 ,  3 )  ∈ R 3 denote the position of a particle in the 3D Cartesian coordinate system.The motion of a particle moving in the 3D time-dependent ABC flow is described by the following SDE, where  1, ,  2, and  3, are independent Brownian motions.For  = 0 and  = 0, the velocity field in (5.3) corresponds to the standard ABC flow [8].The ABC flow is a three-dimensional incompressible velocity field which is an exact solution to the Euler's equation.It is notable as a simple example of a fluid flow that can have chaotic trajectories.In our numerical experiments, we set  =  =  = 1.
In Figure 2, we show the time evolution of the E(1()) 2

2𝑡
for different  0 's (here  0 =  2 /2) and for four different 's, where the result in Figure 2d corresponds to the time-independent ABC flow (see Fig. 3 of [33]).Again the parameter  controls the strength of the time dependence.For each  0 and , we use   = 240 000 particles to solve the SDE (5.3).We find that for each given  0 , the time evolution of the E(1()) 2 2 converges when  converges to zero.However, we observe two different patterns compared with the results shown in Figure 1.First, when we decrease  0 , it takes a longer time for the system to enter a mixing stage.Second, we observe a large amount of enhanced diffusion when  0 decreases.
To further investigate the dependence of   11 on  0 and , in Figure 3, we show the dependence of effective diffusivity   11 on  0 and .We find that for each given  0 , as  decreases the corresponding effective diffusivity   11 converges to the effective diffusivity   11 associated with  = 0. Thus, the time-dependent ABC flow has a similar convection-enhanced diffusion behavior as the time-independent ABC flow.The fitted slope within  0 ∈ [10 −5 , 10 −1 ] is about −1.0, which indicates that   11 ∼ (1/ 1 0 ).This result indicates that the   11 of the time-dependent ABC flow achieves the upper-bound of equation (2.6), i.e. the maximal enhancement.This maximal enhancement phenomenon may be attributed to the ballistic orbits of the ABC flow, where the time-independent case was discussed in [24,34].Moreover, our result for  0 ∈ [10 −3 , 10 −1 ] and  = 0 recovers the same phenomenon as the Figure 2 in [4], which was obtained by using the Eulerian framework, i.e., solving a cell problem.In Figure 3, our method can be easily used to compute the effective diffusivity when  0 ∈ [10 −5 , 10 −4 ].It will be, however, extremely expensive for the Eulerian framework since one needs to solve a convection-dominated PDE (2.4)    where  =  =  = 1 and Ω is the frequency.Here we first choose  = 2 −7 ,   = 240 000 and  = 10 5 .Then, we choose different Ω and compute the corresponding effective diffusivity   11 .In Figure 5, we show the numerical results.We find that when Ω is near 0.1 the diffusion enhancement is weak.When Ω is away from 0.1, say Ω < 0.05 or Ω > 0.2, we observe the maximal enhancement phenomenon.A similar sensitive dependence on the frequency of time-dependent ABC flows was reported in [5], where the Lyapunov exponent of the deterministic time-dependent ABC flow problem (i.e.,  = 0 in Eq. (5.3)) was studied as the indicator of the extent of chaos; see Figures 2 and 3 of [5].
When Ω = 0, the flow of equation (5.4) is the same as that for  = 0 case in equation (5.3), which will give the maximal enhancement phenomenon.When Ω is positive, the flow becomes time-dependent and the regions of chaos expand until the extent of chaos (i.e. the Lyapunov exponent) appears to reach a maximum, which is corresponding to Ω = 0.1.It seems that the diffusion enhancement is significantly weakened in this range of Ω.When Ω continues to grow, the islands of the integrability regrow and the chaotic regions have shrunk significantly.We again observe the maximal enhancement phenomenon in this range of Ω.Our numerical results suggest that the level of chaos and the strength of diffusion enhancement seem to compete with each other.More intensive theoretic and numerical studies will be reported in our future work.

Conclusion
In this paper, we developed a stochastic structure-preserving Lagrangian scheme in computing effective diffusivity of passive tracer models in 3D time-dependent chaotic flows and provided a sharp convergence analysis on the proposed numerical scheme.Our convergence analysis is based on a probabilistic approach, which interprets the solution process generated by our numerical scheme as a Markov process.By exploring the ergodicity of the solution process, we gave a sharp and uniform-in-time error estimate for our numerical scheme, which allows us to compute the effective diffusivity over infinite time.Numerical results verify that the proposed method is robust and accurate in computing effective diffusivity of time-dependent chaotic flows.We observed the maximal enhancement phenomenon in time-dependent ABC flows and the sub-maximal enhancement phenomenon in time-dependent Kolmogorov flows, respectively.Moreover, we found that the time dependency in the velocity field improves the chaotic property of ABC flow and Kolmorogov flow though, it does not change the pattern of convection-enhanced diffusion in both flows.
There are two directions we plan to explore in our future work.First, we intend to study the convectionenhanced diffusion phenomenon and provide a sharp convergence analysis for general time-dependent chaotic flows, where the flows have a quasi-periodic property in the time domain.In addition, we shall investigate the convection-enhanced diffusion phenomenon for general spatial-temporal stochastic flows [20,23] and develop convergence analysis for the corresponding numerical methods.

Figure 3 .
Figure 3. Convection-enhanced diffusion with a maximal enhancement in the time-dependent ABC flow.

Figure 5 .
Figure 5. Dependence of   11 on the frequency of the time-dependent ABC flow.
.35) Noticing that if  1 and v1 are bounded in sup norm, right-hand-side of equation (4.35) is bounded for any .Other terms on the right-hand side of equation (4.34) are also bounded, which can be checked easily.Therefore, we can prove that 1  E in 3D space, whose Péclet number is proportion to 1 0 .Finally, we investigate the dependence of   11 on the frequency of the time-dependent ABC flow.Specifically, we solve the following SDE,