# Determinism and invariant measures for diffusing passive scalars advected by
unsteady random shear flows

Lingyun Ding<EMAIL_ADDRESS>Richard M. McLaughlin<EMAIL_ADDRESS>Department of Mathematics, University of North Carolina, Chapel Hill, NC,
27599, United States

###### Abstract

Here we study the long time behavior of an advection-diffusion equation with a
general time varying (including random) shear flow imposing no-flux boundary
conditions on channel walls. We derive the asymptotic approximation of the
scalar field at long times by using center manifold theory. We carefully
compare it with existing time varying homogenization theory as well as other
existing center manifold based studies, and present conditions on the flows
under which our new approximations give a substantial improvement to these
existing theories. A recent study [30] has shown that Gaussian random shear
flows induce a deterministic effective diffusivity at long times, and
explicitly calculated the invariant measure. Here, with our established
asymptotic expansions, we not only concisely demonstrate those prior
conclusions for Gaussian random shear flows, but also generalize the
conclusions regarding determinism to a much broader class of random (non-
Gaussian) shear flows. Such results are important ergodicity-like results in
that they assure an experimentalist need only perform a single realization of
a random flow to observe the ensemble moment predictions at long time. Monte-
Carlo simulations are presented illustrating how the highly random behavior
converges to the deterministic limit at long time. Counterintuitively, we
present a case demonstrating that the random flow may not induce larger
dispersion than its deterministic counterpart, and in turn present rigorous
conditions under which a random renewing flow induces a stronger effective
diffusivity.

###### keywords:

Passive scalar , Scalar intermittency, Shear dispersion , Random shear flow ,
Turbulent transport , Ergodicity

###### MSC:

[2010]37A25, 37H10, 37N10, 82C70, 76R50

††journal: Physica D (Nonlinear Phenomena)

##  1 Introduction

An extremely important class of problems concerns how fluid motion transports
a diffusing scalar. Since G. I. Taylor [57] first introduced the calculation
showing that a steady pressure driven flow in a pipe leads to a greatly
enhanced effective diffusivity, the literature on this topic has exploded in
many directions spanning many disciplines. Shortly following G. I. Taylor,
Aris [4] presented an alternative approach for shear layers yielding a
hierarchy for the spatial moments of the scalar field. More recent moment
analysis shows how the boundary geometry of the pipe can be used to control
the distribution of solute which is advected by the pressure driven flow [1,
2, 3].

Unsteady flows typically generate different properties than their steady
counterparts. Practical examples of unsteady flow include pulsatile blood
flows [54] and tidal estuaries [26]. The first investigation of the Taylor
dispersion in time-dependent flow dates back to Aris [5], who presented the
study of a solute advected by pulsating flow in a circular tube. After that,
based on the Aris’ moment method, a number of studies reported on the enhanced
diffusivity induced by the single-frequency pulsating flow[15, 25, 63, 47,
38], the single frequency Couette-Poiseuille [10, 49, 12, 11] and the multi-
frequency flow [61, 62, 29]. Alternative approaches, using center manifold
theory, [45, 46, 42, 43] or homogenization methods [27, 28] not only predict
the effective diffusivity but also gives the direct expression for the full
concentration field at long times.

We notice that three points haven’t been addressed well in the literature
regarding to the shear dispersion in time-varying flows. First, most of those
theoretical studies focused on the cross-sectional averaged concentration,
while fewer studies have explored asymptotic corrections which capture cross-
channel variations. Here, with the center manifold theory, we present a
systematic procedure to construct an approximation to capture the traverse
variation of the scalar field. Second, several interesting articles [45, 46,
42] implemented center manifold theory for such unsteady problems employing
certain slowly varying assumptions to simplify the calculation. Such
assumptions restrict the applicability of the effective dynamics. Here we
relax this assumption by carefully incorporating the temporal fluctuation of
the flows into the analysis. Hence, our results can handle rapidly fluctuating
flows or even random flows. Third, recent results have explicitly calculated
using statistical moment closure the invariant measure for a diffusing passive
scalar advected by a class of random shear flows [30, 19] employing no-flux
boundary conditions on channel domains. These results generalize prior
turbulent intermittency in free space of Majda [39] and Kraichnan [41].
Interestingly, we establish here how center manifold theory can be used to
greatly extend these theories to a much broader class of random shear flows,
particularly regarding their temporal statistics. In doing so, we can extend
results which show how all the effective diffusion coefficients converge to a
deterministic value for this broader class of flows, in sharp contrast to the
free-space analog considered by Majda and others [44, 18] in which the
effective diffusivity is random at all times. Such results are important
ergodicity-like results in that they assure an experimentalist need only
perform a single realization of a random flow to observe the field moment
predictions at long time.

The paper is organized as follows. In section 2, we formulate the governing
equation of the shear dispersion problem and review the Aris moment method. In
section 3, we discuss the procedure of applying center manifold theory to the
Taylor dispersion problem with time-varying shear flow. By utilizing the
first-order approximation of the cross-sectional averaged concentration, we
present a nonnegative asymptotic expansion of the scalar field at long times
which captures the transverse variations. We document situations in which a
time varying cell problem produces a more accurate approximation than the
parametric (adiabatic) approach employed recently [45, 46, 42]. In section 4,
we demonstrate that a class of flows with finite correlation time will induce
a deterministic effective diffusivity at long times. Moreover, we establish
conditions which guarantee that the periodic in time problem always yields a
weaker effective diffusivity than the random counterpart. With the derived
effective equation, we computed the explicit formula of invariant measure of
the random passive field.

##  2 Setup and background of the problem

###  2.1 Governing Equation and Nondimensionalization

####  2.1.1  Advection-diffusion Equation

We consider the problem in a channel domain
$(x,\mathbf{y})\in\mathbb{R}\times\Omega$, where the $x$-direction is the
longitudinal direction of the channel and $\Omega\subset\mathbb{R}^{d}$ stands
for the cross section of the channel. Some practical examples of the boundary
geometry includes the parallel-plate channel
$\Omega=\left\\{y|y\in[0,L]\right\\}$, the circular pipe
$\Omega=\left\\{\mathbf{y}|\mathbf{y}^{2}\leq L\right\\}$, the rectangular
duct $\Omega=\left\\{\mathbf{y}|\mathbf{y}\in[0,L]^{2}\right\\}$, bowed
rectangular channels [40]. The passive scalar is governed by the advection-
diffusion equation with a general time-varying shear flow $u(\mathbf{y},t)$
and no-flux boundary condition which takes the form

$\partial_{t}T+v(\mathbf{y},t)\partial_{x}T=\kappa\Delta T,\quad
T(x,\mathbf{y},0)=T_{I}(x,\mathbf{y}),\quad\left.\frac{\partial
T}{\partial\mathbf{n}}\right|_{\mathbb{R}\times\partial\Omega}=0,$ (1)

where $\kappa$ is the diffusivity, $T_{I}(x,\mathbf{y})$ is the initial data,
$\mathbf{n}$ is the outward normal vector of the boundary
$\mathbb{R}\times\partial\Omega$ and $\partial\Omega$ is the boundary of
$\Omega$.

####  2.1.2 Nondimensionalization

With the change of variables

$\displaystyle Lx^{\prime}=x,\quad
L\mathbf{y}^{\prime}=\mathbf{y},\quad\frac{L^{2}}{\kappa}t^{\prime}=t,\quad
Uv^{\prime}(\mathbf{y}^{\prime},t^{\prime})=v(\mathbf{y},t),\quad
L\Omega^{\prime}=\Omega,$ (2) $\displaystyle
T_{I}^{\prime}(x^{\prime},\mathbf{y}^{\prime})L^{-d-1}\int\limits_{\mathbb{R}\times\Omega}T_{I}(x,\mathbf{y})\mathrm{d}x\mathrm{d}\mathbf{y}=T_{I}(x,\mathbf{y}),$
$\displaystyle
T^{\prime}(x^{\prime},\mathbf{y}^{\prime},t^{\prime})L^{-d-1}\int\limits_{\mathbb{R}\times\Omega}T_{I}(x,\mathbf{y})\mathrm{d}x\mathrm{d}\mathbf{y}=T(x,\mathbf{y},t),\;$

after dropping the primes, we obtain the nondimensionalized advection-
diffusion equation

$\partial_{t}T+\mathrm{Pe}v(\mathbf{y},t)\partial_{x}T=\Delta
T,\;T(x,\mathbf{y},0)=T_{I}(x,\mathbf{y}),\;\left.\frac{\partial
T}{\partial\mathbf{n}}\right|_{\mathbb{R}\times\partial\Omega}=0,$ (3)

where $\mathrm{Pe}={UL}/{\kappa}$ is the Péclet number.

###  2.2 Aris Moment Hierarchy

Aris showed in [4] that one could write down a recursive system of partial
differential equations (3) for the spatial moments of the tracer $T$. The
$n$th Aris moment is defined by
$T_{n}(\mathbf{y},t)=\int\limits_{-\infty}^{\infty}x^{n}T(x,\mathbf{y},t)\mathrm{d}x$.
With the assumption $T(\pm\infty,\mathbf{y},t)=0$, the Aris moments satisfy
the recursive relationship called the Aris equations,

$\displaystyle(\partial_{t}-\Delta)T_{n}=n(n-1)T_{n-2}+n\mathrm{Pe}v(\mathbf{y},t)T_{n-1},\;T_{n}(\mathbf{y},0)=\int\limits_{-\infty}^{\infty}x^{n}T_{I}(x,\mathbf{y})\mathrm{d}x,\;\left.\frac{\partial
T}{\partial\mathbf{n}}\right|_{\partial\Omega}=0,$ (4)

where $T_{n}=0$ if $n\leq-1$. The full moments of $T$ are then obtained though
the cross-sectional average of the moments
$\bar{T}_{n}=\frac{1}{\left|\Omega\right|}\int\limits_{\Omega}T_{n}\mathrm{d}\mathbf{y}$,
where $\left|\Omega\right|$ is the area of $\Omega$. In this following
context, we use the overline to denote the cross sectional average. Applying
the divergence theorem and boundary conditions yield

$\displaystyle\frac{\mathrm{d}\bar{T}_{n}}{\mathrm{d}t}=n(n-1)\bar{T}_{n-2}+n\mathrm{Pe}\overline{v(\mathbf{y},t)T_{n-1}},\;\bar{T}_{n}(0)=\frac{1}{\left|\Omega\right|}\int\limits_{\Omega}\int\limits_{-\infty}^{\infty}x^{n}T_{I}(x,\mathbf{y})\mathrm{d}x\mathrm{d}\mathbf{y}.$
(5)

The homogenization method in [29, 20] suggests that, assuming a scale
separation in the initial data, the solution of equation (1) can be
approximated by a diffusion equation with an effective diffusion coefficient.
The effective longitudinal effective diffusivity could be computed through the
Aris moments

$\kappa_{\mathrm{eff}}=\lim\limits_{t\rightarrow\infty}\frac{\mathrm{Var}(\bar{T})}{2t},$
(6)

where $\mathrm{Var}(\bar{T})=\bar{T}_{2}-\bar{T}_{1}^{2}$ is the variance of
the cross-sectional average $\bar{T}$. In this paper, we use
$\kappa_{\mathrm{eff}}$ to denote the dimensional effective diffusivity
computed by the dimensional Aris moment and use the
$\tilde{\kappa}_{\mathrm{eff}}=\kappa_{\mathrm{eff}}/\kappa$ to denote the
non-dimensional effective diffusivity.

The effective diffusivity characterizes the symmetric property of the
longitudinal distribution. We are also interested in the asymmetry properties
of $\bar{T}$. Skewness is the lowest order integral measure of the asymmetry
of a real-valued probability distribution, which is defined as

$\mathrm{S}(\bar{T})=\frac{\bar{T}_{3}-3\bar{T}_{2}\bar{T}_{1}+2\bar{T}_{1}^{3}}{\left(\bar{T}_{2}-\bar{T}_{1}^{2}\right)^{\frac{3}{2}}}.$
(7)

The information of shape provided by the skewness could improve the design of
microfluidic flow injection analysis [3, 58] and chromatographic separation
[14].

##  3  Center manifold description of the shear dispersion problem

###  3.1 Center manifold and reduction principle

In pioneering work, Mercer and Roberts [45] interpreted the long time
asymptotic of the shear dispersion problem as the center manifold of a
dynamical system, which provides a systematic and near rigorous approach to
derive the approximation. Besides the shear dispersion problem, the pratical
applications of center manifold theory include chromatographic model and
reactors[9], elastic beam deformations[50], and thin fluid flows dynamics [51,
53]. To explain the center manifold method, let’s consider an autonomous
differential system of the form

$\displaystyle\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=A\mathbf{x}+f(\mathbf{x},\mathbf{y}),\quad\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}=B\mathbf{y}+g(\mathbf{x},\mathbf{y}),$
(8)

where $\mathbf{x}\in\mathbb{R}^{m},\mathbf{y}\in\mathbb{R}^{n}$. $A,B$ are
matrices whose eigenvalues have vanishing and negative real parts,
respectively. $f(\mathbf{x},\mathbf{y})$, $g(\mathbf{x},\mathbf{y})$ and their
first order partial derivatives are zero at
$\mathbf{x}=\mathbf{0},\mathbf{y}=\mathbf{0}$. These conditions grantee the
existence of a center manifold $\mathbf{y}=h(\mathbf{x})$ which has two
important features. First, the stability properties of the dynamical system
(8) at the origin are shared by the following lower dimensional equation

$\displaystyle\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=A\mathbf{x}+f(\mathbf{x},h(\mathbf{x})).$
(9)

Second, in case of a stable equilibrium
$(\mathbf{x},\mathbf{y})=(\mathbf{0},\mathbf{0})$ each solution of system (8)
which starts close to the origin exponentially decays to a particular solution
on the center manifold [21, 22]. With these two features of the center
manifold, one can reduce the original $m+n$-dimensional system (8) to a
$m$-dimensional system (9) with only the price of exponential corrections.

This classical center manifold theory and reduction principle could be
generalized in many directions. First, the dynamical system (8) could be an
infinite dimensional system where the matrices $A$,$B$ become linear operators
[23]. Second, similar results hold for a more general dynamical system
$\frac{\mathrm{d}\mathbf{x}_{i}}{\mathrm{d}t}=A_{i}\mathbf{x}_{i}+f(\mathbf{x}_{1},...,\mathbf{x}_{N},t)$,
$1\leq i\leq N$ and the restriction of eigenvalues could be weakened [59, 7].
This generalization leads to a so-called two-mode invariant manifold model for
the shear dispersion problem [65, 64, 56]. Third, more related to our topic,
the system could be non-autonomous, where the center manifold becomes time-
dependent $\mathbf{y}=h(\mathbf{x},t)$ [6, 8, 52]. For further details
regarding center manifold theory, we refer to [21, 7] and references therein.

Notice that the advection-diffusion equation (3) is linear, while the center
manifold theory applies to a system with nonlinear terms. To fit the center
manifold theory, we first apply the Fourier transform on equation (3) and
obtain

$\displaystyle\frac{\partial\hat{T}}{\partial
t}-\mathrm{i}k\mathrm{Pe}v(\mathbf{y},t)\hat{T}=-k^{2}\hat{T}+\Delta_{\mathbf{y}}\hat{T},\;\left.\frac{\partial\hat{T}}{\partial\mathbf{n}}\right|_{\mathbb{R}\times\partial\Omega}=0,\;\hat{T}(k,\mathbf{y},0)=\hat{T}_{I}(k,\mathbf{y}),$
(10)

Second, we conceptually non-linearize equation (10) by treating the wavenumber
as a dependent variable of the dynamical system. Notice that
$\Delta_{\mathbf{y}}$ has a null space which consists of all function
independent on $y$. To fit the form of equation (8), we rewrite equation (10)
as

$\displaystyle\partial_{t}\begin{bmatrix}k\\\
\hat{\overline{T}}\end{bmatrix}=\begin{bmatrix}0&0\\\
0&\Delta_{\mathbf{y}}\end{bmatrix}\begin{bmatrix}k\\\
\hat{\overline{T}}\end{bmatrix}+\begin{bmatrix}0\\\
\mathrm{i}k\mathrm{Pe}\overline{v(\mathbf{y},t)\hat{T}}-k^{2}\hat{\overline{T}}\end{bmatrix},$
(11)
$\displaystyle\hat{T^{\prime}}=\Delta_{\mathbf{y}}\hat{T^{\prime}}+\mathrm{i}k\mathrm{Pe}v(\mathbf{y},t)\hat{T^{\prime}}-k^{2}\hat{T}-\overline{\mathrm{i}k\mathrm{Pe}v(\mathbf{y},t)\hat{T}},$
$\displaystyle\left.\frac{\partial\hat{T}}{\partial\mathbf{n}}\right|_{\mathbb{R}\times\partial\Omega}=0,\quad\hat{T}(k,\mathbf{y},0)=\hat{T}_{I}(k,\mathbf{y}),$

where $T^{\prime}(x,\mathbf{y},t)$ and $\overline{T}(x,t)$ are the fluctuation
and average of $T(x,\mathbf{y},t)$ with respect to $\mathbf{y}$. This system
admits a center manifold $\hat{T}^{\prime}=h(\hat{\overline{T}},k,t)$. Based
on center manifold theory, $\hat{T}$ converges to
$h(\hat{\overline{T}},k,t)+\hat{\overline{T}}$ exponentially as
$t\rightarrow\infty$. Due to the diffusion effect, $T$ is a decaying scalar
field. The energy concentrates near the neighborhood of $k=0$ at long times.
Hence, we can seek the expansion of $h(\hat{\overline{T}},k,t)$ for small $k$
and $\hat{\overline{T}}$,
$h=\sum\limits_{n=1}^{\infty}h_{n}(\mathbf{y},t)k^{n}\hat{\overline{T}}+\mathcal{O}(\hat{\overline{T}}^{2})$.
That is equivalent to approximating the scalar field $T$ by the derivatives of
its cross-sectional average $\bar{T}$ with respect to $x$. This idea dated
back to Gill [32, 33] and also has been discussed in [68].

For simplicity, we rewrite all equations in term of physical variables. The
governing equations are

$\displaystyle\partial_{t}\bar{T}$
$\displaystyle=\partial_{x}^{2}\bar{T}-\mathrm{Pe}\overline{v(\mathbf{y},t)\partial_{x}T},$
(12) $\displaystyle\partial_{t}T$
$\displaystyle=\Delta_{\mathbf{y}}T+\partial_{x}^{2}T-\mathrm{Pe}v(\mathbf{y},t)\partial_{x}T.$

The expansion becomes

$\displaystyle T$
$\displaystyle=T^{\prime}+\bar{T}=\bar{T}+h(\bar{T})=\sum\limits_{n=0}^{\infty}\theta_{n}(\mathbf{y},t)\partial_{x}^{n}\bar{T}.$
(13)

The fluctuation is mean zero, $\int_{\Omega}T^{\prime}\mathrm{d}\mathbf{y}=0$,
which implies $\bar{\theta}_{0}=1$ and $\bar{\theta}_{n}=0$ if $n\geq 1$ at
long times. We have
$\left.\frac{\partial}{\partial\mathbf{n}}\theta_{n}\right|_{\mathbf{y}\in\partial\Omega}=0$
from the no-flux boundary conditions of $T$. Substituting expansion (13) into
equation (12), we have

$\displaystyle\partial_{t}\bar{T}=\partial_{x}^{2}\bar{T}-\mathrm{Pe}\overline{v\partial_{x}T},$
(14a)
$\displaystyle\sum\limits_{n=0}^{\infty}\partial_{t}\theta_{n}\partial_{x}^{n}\bar{T}+\sum\limits_{n=0}^{\infty}\theta_{n}\partial_{x}^{n}\partial_{t}\bar{T}=\Delta_{\mathbf{y}}T+\partial_{x}^{2}T-\mathrm{Pe}v\partial_{x}T.$
(14b)

Grouping all terms of the same order, namely $\partial_{x}^{n}\bar{T}$, we
find that we have to solve the sequence of equations

$\displaystyle\left(\partial_{t}-\Delta_{\mathbf{y}}\right)\theta_{0}=0,$ (15)
$\displaystyle\left(\partial_{t}-\Delta_{\mathbf{y}}\right)\theta_{1}=-\mathrm{Pe}\theta_{0}\left(v-\overline{\theta_{0}v}\right),$
$\displaystyle\left(\partial_{t}-\Delta_{\mathbf{y}}\right)\theta_{n}=-\mathrm{Pe}v\theta_{n-1}+\mathrm{Pe}\sum\limits_{m=0}^{n-1}\theta_{n-m-1}\overline{v\theta_{m}},$

where $\theta_{n}=0$ if $n<0$. After we solve $\theta_{n}$ successively, we
obtain the closed evolution equation of $\bar{T}$ by substituting
$T=\sum\limits_{n=0}^{\infty}\theta_{n}(\mathbf{y},t)\partial_{x}^{n}\bar{T}$
into equation (12),

$\displaystyle\partial_{t}\bar{T}$
$\displaystyle=\partial_{x}^{2}\bar{T}-\mathrm{Pe}\sum\limits_{n=0}^{\infty}\overline{v\theta_{n}}\partial_{x}^{n+1}T.$
(16)

Finlay, once we solve equation (16) for $\bar{T}$, we obtain the approximation
of the scalar field $T$ via expansion (13).

###  3.2 The first and second order effective equation

In this subsection, we will compute equation (15) and (16) for the flow
$v(\mathbf{y},t)=\xi(t)u(\mathbf{y})$. For more general non-separable flow
$v(\mathbf{y},t)$, one could reduce it to a separable form by utilizing the
Fourier transform in time. To simplify the calculation, we assume
$T_{I}(x,\mathbf{y})=\delta(x)$. Otherwise, the general initial condition only
creates extra exponential decaying terms and yields the same asymptotic
expansion at long times.

With the constraints of the average and boundary conditions of $\theta_{n}$,
we have $\theta_{0}=1$. Therefore, the equation of $\theta_{1}$ becomes

$\displaystyle\left(\partial_{t}-\Delta_{\mathbf{y}}\right)\theta_{1}=-\mathrm{Pe}(v-\bar{v}),\quad\left.\frac{\partial}{\partial_{\mathbf{n}}}\theta_{1}\right|_{\partial\Omega}=0,$
(17)

which is identical to equation (4) in the Aris moments calculation. Since the
theory concerns the long time dynamics of the scalar field and the long time
limit of $\theta_{1}$ doesn’t depend on the initial condition, in principle,
one can solve equation (17) with arbitrary initial condition. To obtain a
better approximation at earlier stage, one can choose suitable initial
condition of $\theta_{n}$ to match both sides of the expansion (13) at $t=0$.
Then when $v(\mathbf{y},t)=\xi(t)u(\mathbf{y})$, the solution of equation (17)
is

$\displaystyle\theta_{1}(\mathbf{y},t)=-\mathrm{Pe}\sum\limits_{n=1}^{\infty}\phi_{n}\left\langle
u,\phi_{n}\right\rangle\int_{0}^{t}e^{\lambda_{n}(s-t)}\xi(s)ds.$ (18)

where $\left\langle
f,g\right\rangle=\frac{1}{|\Omega|}\int\limits_{\Omega}fg\mathrm{d}\mathbf{y}$.
$\phi_{n},\lambda_{n}$ are the eigenfunctions and eigenvalues of the Laplace
operator in the cross section of the channel $\Omega$ with no-flux boundary
condition, i.e.,

$\displaystyle-\Delta\phi_{n}=\lambda_{n}\phi_{n},\quad\left.\frac{\partial}{\partial_{\mathbf{n}}}\phi_{n}\right|_{\partial\Omega}=0,\quad\left\langle\phi_{n},\phi_{n}\right\rangle=1.$
(19)

Substituting $T=\bar{T}+\theta_{1}\partial_{x}\bar{T}$ into the evolution
equation of $\bar{T}$, we obtain the first order effective equation

$\displaystyle\partial_{t}\bar{T}+\mathrm{Pe}\bar{v}\partial_{x}\bar{T}=a_{2}\partial_{x}^{2}\bar{T},\quad
a_{2}=\left(1-\mathrm{Pe}\overline{v\theta_{1}}\right).$ (20)

The classical homogenization approach relies on the Fredholm alternative which
involves a space-time average. As a result, the effective equation is a
constant coefficient equation even for the time-varying flow case [29, 27].
Here, with the center manifold approach, we obtain the effective equation (20)
with time-dependent coefficients which could approximate the scalar field
better in an earlier stage. Comparing with the definition of Aris moments and
variance of the cross-sectional average, we have

$\displaystyle\mathrm{Var}(\bar{T})=\mathrm{Var}(\bar{T_{I}})+2\int\limits_{0}^{t}a_{2}(s)\mathrm{d}s$
(21)

For a periodic time-varying flow [61] and a class of random flows [19, 30], we
have $\mathrm{Var}(\bar{T})=2\kappa_{\mathrm{eff}}t+\mathcal{O}(1)$, where
$\kappa_{\mathrm{eff}}$ is the effective diffusivity

$\displaystyle\quad\kappa_{\mathrm{eff}}=\lim\limits_{t\rightarrow\infty}\frac{\mathrm{Var}(\bar{T})}{2t}=\lim\limits_{t\rightarrow\infty}\frac{1}{t}\int\limits_{0}^{t}a_{2}(s)\mathrm{d}s.$
(22)

In other words, $a_{2}$ can be approximated by its time average at long times.

With the expression (18), the effective diffusivity induced by the flow
$v(\mathbf{y},t)=\xi(t)u(\mathbf{y})$ is

$\displaystyle
a_{2}=\left(1-\mathrm{Pe}\overline{v\theta_{1}}\right)=1+\mathrm{Pe}^{2}\sum\limits_{n=1}^{\infty}\left\langle
u,\phi_{n}\right\rangle^{2}\xi(t)\int\limits_{0}^{t}e^{\lambda_{n}(s-t)}\xi(s)\mathrm{d}s,$
(23)
$\displaystyle\kappa_{\mathrm{eff}}=1+\mathrm{Pe}^{2}\lim\limits_{t\rightarrow\infty}\frac{1}{t}\int\limits_{0}^{t}\sum\limits_{n=1}^{\infty}\left\langle
u,\phi_{n}\right\rangle^{2}\xi(s_{2})\int\limits_{0}^{s_{2}}e^{\lambda_{n}(s_{1}-s_{2})}\xi(s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}.$

With the initial condition $T(x,\mathbf{y},0)=\delta(x)$, the solution of
equation (20) gives an approximation of $\bar{T}$ as $t\rightarrow\infty$,

$\displaystyle\bar{T}(x,t)=\frac{1}{\sqrt{4\pi
b_{2}}}\exp\left(\frac{-\tilde{x}^{2}}{4b_{2}}\right)+\mathcal{O}(t^{-\frac{3}{2}}),$
(24) $\displaystyle
b_{2}=\int\limits_{0}^{t}a_{2}(s)\mathrm{d}s,\quad\tilde{x}=x-\mathrm{Pe}\int\limits_{0}^{t}\bar{v}(s)\mathrm{d}s.$

For steady flow, we have $b_{2}=\kappa_{\mathrm{eff}}t$. Then equation (24)
reduces to the classical Gaussian approximation [24]. Since the scalar field
will be homogenized across the channel at long times, $\bar{T}$ itself could
be an approximation of $T$. In fact, we could obtain a more accurate
approximation of $T$,

$\displaystyle
T\approx\bar{T}+\theta_{1}\partial_{x}\bar{T}=\left(1-\frac{\theta_{1}(\mathbf{y},t)\tilde{x}}{2b_{2}}\right)\frac{1}{\sqrt{4\pi
b_{2}}}\exp\left(\frac{-\tilde{x}^{2}}{4b_{2}}\right)+\mathcal{O}(t^{-\frac{3}{2}}).$
(25)

Since $\partial_{x}\bar{T}$ is an odd function with respect to $x$, the error
of approximation (25) is still $\mathcal{O}(t^{-\frac{3}{2}})$. However,
equation (25) practically performs better than $\bar{T}$ because it depicts
the across channel variation of the scalar. Figure 1 presents the relative
errors of different approximations for the solution studied in right panel of
figure 3, where the metric of error is ${\lVert
T-T_{approx}\rVert_{\infty}}/{\lVert T\rVert_{\infty}}$. As shown in figure 1,
the relative error of approximation (24) (red curve) is around $0.1$ at $t=1$,
while, the relative error of approximation (25) (blue curve) is around
$10^{-3}$. Since two approximations are of the same asymptotic order at long
times, presumably the differences between the two approximations will reduce
as time is further increased.

Figure 1: The relative errors of different approximations for the solution
studied in right panel of figure 3. The red solid curve, blue dashed curve,
black dashed dot curve and purple curve are the relative error of
approximation (24), (25), (27) and (26), respectively.

In many applications, the scalar field usually stands for the concentration
which must be nonnegative for all times. However, this approximation (25)
could be negative for some $x$ and $t$, which may not be desirable in those
applications. [67] proposed the following nonnegative approximation to study
the transverse distribution of concentration distribution for laminar tube
flow,

$\displaystyle T\approx$ $\displaystyle\frac{1}{\sqrt{4\pi
b_{2}}}\exp\left(\frac{-(\tilde{x}-\theta_{1}(\mathbf{y},t))^{2}}{4b_{2}}\right).$
(26)

The approximation (26) converges asymptotically to approximation (25) as
$\theta\rightarrow 0$. However, the relative difference between them doesn’t
vanish as $t\rightarrow\infty$. As shown in figure 1, there is a visible
difference between the approximation (26) (purple curve) and (25) (blue
curve).

Here, we propose a nonnegative asymptotic expansion

$\displaystyle T\approx$
$\displaystyle\left(1-\frac{\theta_{1}(\mathbf{y},t)\tilde{x}}{4b_{2}}\right)^{2}\frac{1}{\sqrt{4\pi
b_{2}}}\exp\left(\frac{-\tilde{x}^{2}}{4b_{2}}\right),\quad
t\rightarrow\infty.$ (27)

Since the difference between equation (25) and (27) is
$\mathcal{O}(t^{-\frac{5}{2}})$, the relative difference between them vanishes
as $t\rightarrow\infty$. From figure 1, we can see that the relative error of
approximation (25) and (27) is almost indistinguishable after $t=0.1$.

Next, we study the second order approximation of the scalar field. We have to
solve the equation for $\theta_{2}$,

$\displaystyle\left(\partial_{t}-\Delta_{\mathbf{y}}\right)\theta_{2}=-\mathrm{Pe}(v\theta_{1}-\theta_{1}\bar{v}-\overline{v\theta_{1}}),\quad\left.\frac{\partial}{\partial_{\mathbf{n}}}\theta_{2}\right|_{\partial\Omega}=0,$
(28)

We have the expansion of
$v\theta_{1}-\theta_{1}\bar{v}-\overline{v\theta_{1}}$,

$\displaystyle v\theta_{1}-\theta_{1}\bar{v}-\overline{v\theta_{1}}=$
$\displaystyle\sum\limits_{n_{2},n_{1}=1}^{\infty}\left\langle\theta_{1},\phi_{n_{1}}\right\rangle\left\langle\phi_{n_{1}}(v-\bar{v}),\phi_{n_{2}}\right\rangle\phi_{n_{2}}.$
(29)

That leads to the solution

$\displaystyle\theta_{2}=$
$\displaystyle\mathrm{Pe}^{2}\sum\limits_{n_{2},n_{1}=1}^{\infty}\left\langle
u,\phi_{n_{1}}\right\rangle\left\langle\phi_{n_{1}}(u-\bar{u}),\phi_{n_{2}}\right\rangle\phi_{n_{2}}\int\limits_{0}^{t}\left(e^{\lambda_{n_{2}}(s_{2}-t)}\xi(s_{2})\int_{0}^{s_{2}}e^{\lambda_{n_{1}}(s_{1}-s_{2})}\xi(s_{1})ds_{1}\right)\mathrm{d}s_{2}$
(30)

Substituting
$T=\bar{T}+\theta_{1}\partial_{x}\bar{T}+\theta_{2}\partial^{2}_{x}\bar{T}$
into the evolution equation of $\bar{T}$, the approximated evolution equation
for $\bar{T}$ becomes a linearized Burgers-Korteweg-de Vries equation

$\displaystyle\partial_{t}\bar{T}+\bar{v}\partial_{x}\bar{T}=a_{2}\partial_{x}^{2}\bar{T}-a_{3}\partial_{x}^{3}\bar{T},\quad
a_{3}=\mathrm{Pe}\overline{v\theta_{2}}.$ (31)

Then we can consider two cases based on
$a_{3}=\mathrm{Pe}\overline{v\theta_{2}}$. First, we consider the case in
involving $a_{3}=\mathrm{Pe}\overline{v\theta_{2}}=0$, which implies the
skewness of $\bar{T}$ is zero. One such example is the linear shear flow
created by moving one boundary of parallel-plate channel [29]. In this case,
the evolution equation for $\bar{T}$ reduces to the diffusion equation, where
the Gaussian approximation (24) is still valid. Then we obtain the
approximation of the whole scalar field

$\displaystyle T$
$\displaystyle=\bar{T}+\theta_{1}\partial_{x}\bar{T}+\theta_{2}\partial^{2}_{x}\bar{T}$
(32)
$\displaystyle=\left(1-\frac{\theta_{1}\tilde{x}}{2b_{2}}+\frac{\theta_{2}\left(\tilde{x}^{2}-2b_{2}\right)}{4b_{2}^{2}}\right)\frac{1}{\sqrt{4\pi
b_{2}}}\exp\left(\frac{-\tilde{x}^{2}}{4b_{2}}\right)+\mathcal{O}(t^{-2}).$

Since $\partial_{x}^{2}\bar{T}$ is an even function with respect to $x$, the
error of approximation (32) is $\mathcal{O}(t^{-2})$ which is more accurate
than the approximation (25).

To verify the validity of the approximation (32), we compare it with the
numerical solution of equation (3) with the flow $v(y,t)=\cos\pi y$. Solving
equation (17) and (28), we have

$\displaystyle\theta_{1}=-\mathrm{Pe}\frac{\cos\pi
y}{\pi^{2}},\quad\theta_{2}=\frac{\mathrm{Pe}^{2}\cos(2\pi y)}{8\pi^{4}}.$
(33)

To fit the initial condition $T_{I}$, we can impose the initial condition
$\theta_{1}(y,0)=\theta_{2}(y,0)=0$ and obtain the time-dependent solutions,

$\displaystyle\theta_{1}=-\mathrm{Pe}\frac{\cos\pi
y}{\pi^{2}}\left(1-e^{-\pi^{2}t}\right),\quad\theta_{2}=\frac{\mathrm{Pe}^{2}\cos(2\pi
y)}{8\pi^{4}}\left(1-e^{-4\pi^{2}t}\right).$ (34)

Figure 2 shows the relative error of various different approximations. The
numerical solution is obtained via the method described in detail in appendix
7.1. We have three observations. First, the formula (25) and (32) retaining
cross-sectional variation provide more accurate approximation than (24).
Second, we can see that the second order approximation (32) has smaller error
than the first order approximation (25) at larger time. We expect this
difference will be more pronounced at longer times. Third, if we impose the
initial condition on $\theta_{1}$ and $\theta_{2}$, then we obtain a more
accurate approximation at earlier stage.

Figure 2: The relative error of various different approximations of numerical
solution of equation (3) with the flow $v(y,t)=\cos\pi y$, $\mathrm{Pe}=1$ and
initial condition
$T_{I}=\left(\sqrt{2\pi}\sigma\right)^{-1}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right),\sigma=1/20$.
The red solid, blue dashed and black dot dashed curve represent the relative
error of approximation (24), (25) and (32). Panel (a) $\theta_{1}$ and
$\theta_{2}$ are provided in (33). Panel (b) $\theta_{1}$ and $\theta_{2}$ are
provided in (34)

Next, we consider the case $a_{3}\neq 0$. When the initial condition is
$T_{I}(x,\mathbf{y})=\delta(x)$, the integral representation of the solution
is

$\displaystyle\bar{T}(x,t)=$
$\displaystyle\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\exp(-b_{2}k^{2}-\mathrm{i}b_{3}k^{3}+\mathrm{i}xk)\mathrm{d}k$
(35) $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int\limits_{0}^{\infty}\exp(-b_{2}k^{2})\cos(-b_{3}k^{3}+xk)\mathrm{d}k,$

where $b_{3}=\int\limits_{0}^{t}a_{3}(s)\mathrm{d}s$. We are interested in the
asymptotic expansion of solution (35) at long times. It is a hard task for a
very general time varying flow. Therefore, we restrict our attention to the
case where $a_{2},a_{3}$ are constant. For some time varying flows, we can
approximate $a_{2},a_{3}$ with their time average at long times, for example,
periodic time-varying flow. Hence, the asymptotic expansion we derived in the
section also applies to these cases.

If $x\ll t$ and $t\rightarrow\infty$, the integrand in equation (35) is
localized around $k=0$. Hence, we have the approximation

$\displaystyle\bar{T}(x,t)=$
$\displaystyle\int\limits_{-\infty}^{\infty}\left(1-\mathrm{i}a_{3}k^{3}t+\frac{(-\mathrm{i}a_{3}k^{3}t)^{2}}{2}\right)\exp(-a_{2}k^{2}t+\mathrm{i}xk)+\mathcal{O}(t^{-2})$
(36) $\displaystyle=$
$\displaystyle\left(1-\frac{a_{3}}{2^{3}a_{2}^{\frac{3}{2}}t^{\frac{1}{2}}}H_{3}\left(\frac{x}{2\sqrt{a_{2}t}}\right)+\frac{a_{3}^{2}}{2^{7}a_{2}^{3}t}H_{6}\left(\frac{x}{2\sqrt{a_{2}t}}\right)\right)\frac{\exp\left(\frac{-x^{2}}{4a_{2}t}\right)}{\sqrt{4\pi
a_{2}t}}+\mathcal{O}(t^{-2}),$

where $H_{n}$ is the degree $n$ Hermite polynomial associated with the weight
function $e^{-x^{2}}$. The approximation (36) is identical to the Hermite
polynomial representation proposed in equation (5.7) in [55].

###  3.3 Improvements compared with previous studies

We remark that there are two subtle differences compared with the previous
studies [45, 46, 42]. First, the previous studies made not only the ansatz of
the expansion of $T$, but also the expansion of $\bar{T}$. Therefore, the
recursive equations involve not only $\theta_{n}$, but also the coefficients
in the expansion of $\bar{T}$. Here, we avoid making the expansion ansatz for
$\bar{T}$ by utilizing equation (14a), which simplifies the calculation of
$\theta_{n}$.

Second, in the previous studies, the center manifold are assumed to be time
independent. Hence, the equation for the auxiliary function $\theta_{1}$
derived in [45, 46, 42] takes the form

$\displaystyle-\Delta_{\mathbf{y}}\theta_{1}=-\mathrm{Pe}(u-\bar{u}),$ (37)

in which the time derivative term doesn’t appear. We think the justification
is that the flow $u(\mathbf{y},t)$ varies slowly in time so that the time
derivative term is negligible. However, we think this assumption limits many
naturally arising applications. Let’s consider a simple example,
$\Omega=[0,1]$, $u=-e^{\mathrm{i}\omega t}\cos\pi y$. The solution of equation
(17) is $\frac{e^{\mathrm{i}\omega t}\cos\pi y}{\pi^{2}+\mathrm{i}\omega}$,
while the solution of equation (37) is $\frac{e^{\mathrm{i}\omega t}\cos\pi
y}{\pi^{2}}$. The only difference between them is the wavenumber $\omega$ in
the denominator, which yields a $\mathcal{O}(\omega)$ difference. Hence, for
small wavenumber $\omega$, the two solutions are close. However, for any fixed
$\omega$, the corresponding approximations of the solute distribution $T$
diverge at long times, due to the variances having different growth rates.
Recall that the
variance,$\mathrm{Var}(\bar{T})=2(1-\mathrm{Pe}\overline{v\theta_{1}})t+\mathcal{O}(1)$,
grows linearly at long times. The difference between the two variances arising
from the two different cell problems accumulates and becomes an
$\mathcal{O}(1)$ difference at the frequency time scale
$\mathcal{O}(\frac{1}{\omega})$. Since the solute distribution is
characterized by the variance, the $\mathcal{O}(1)$ difference between
variances implies an $\mathcal{O}(1)$ difference in the distributions at that
time. Moreover, this difference in distributions will keep increasing as time
increases. Hence, we conclude that equation (37) should only be used with a
slow varying flow and before the frequency time scale. In addition, this can
be considered as an example of non-commutating limits.

We know the center manifold becomes a good approximation if the exponential
correction is small, i.e., after the diffusion time scale $L^{2}/\kappa$. If
the frequency time scale is less than the diffusion time scale, then equation
(37) is invalid for all time. That certainly limits the application of the
result based on equation (37). [42, 43] adopted equation (37) to study
dispersion induced by pulsating flows. One of their applications is to blood
flows. Consider the following practical example. The typical frequency time
scale in the human blood vessel is $1$s (60 heartbeats per min). The sodium
chloride ($\kappa\approx 1.6*10^{-5}cm^{2}/s$ in water [35]) diffuses cross
the blood vessel with diameter $0.2$ mm takes around $25$ s. In this case, the
result based on (37) is unlikely valid.

To demonstrate the validity of our analysis, we solve equation (3) numerically
and present the results in figure 3. For the time varying shear flow
$u(y)=\xi(t)y(1-y)/{2}$, [45] derived the effective equation

$\displaystyle\partial_{t}\bar{T}+\frac{\mathrm{Pe}\xi(t)}{12}\partial_{x}\bar{T}=$
$\displaystyle\left(1+\frac{\mathrm{Pe}^{2}\xi(t)^{2}}{30240}\right)\partial_{x}^{2}\bar{T}.$
(38)

If $\xi(t)=\cos\omega t$, the solution of equation (17) is

$\displaystyle\theta=$
$\displaystyle\mathrm{Pe}\sum\limits_{n=1}^{\infty}\frac{(-1)^{n}+1}{\pi^{2}n^{2}}\cos
n\pi
y\left(\frac{\omega\sin(t\omega)}{\pi^{4}n^{4}+\omega^{2}}+\frac{\pi^{2}n^{2}\cos(t\omega)}{\pi^{4}n^{4}+\omega^{2}}-\frac{\pi^{2}n^{2}\text{Pe}e^{-\pi^{2}n^{2}t}}{\pi^{4}n^{4}+\omega^{2}}\right)$
(39)

Hence the effective equation (20) derived by time-dependent center manifold
theory is

$\displaystyle\partial_{t}\bar{T}+\frac{\mathrm{Pe}\cos\omega
t}{12}\partial_{x}\bar{T}=$
$\displaystyle\left(1+\mathrm{Pe}^{2}\sum\limits_{n\in\text{even}^{+}}^{\infty}\frac{2\cos^{2}(t\omega)}{\pi^{2}n^{2}\left(\pi^{4}n^{4}+\omega^{2}\right)}+\frac{\omega\sin(2t\omega)}{\pi^{4}n^{4}\left(\pi^{4}n^{4}+\omega^{2}\right)}\right)\partial_{x}^{2}\bar{T},$
(40)

where we neglect the exponential term in the solution of equation (17). When
$t\gg\frac{1}{\omega}$, we could approximate the series in the effective
equation by its time average

$\displaystyle\partial_{t}\bar{T}+\frac{\mathrm{Pe}\cos\omega
t}{12}\partial_{x}\bar{T}=$
$\displaystyle\left(1+\mathrm{Pe}^{2}\left(\frac{1}{24\omega^{2}}-\frac{\sin\left(\frac{\sqrt{\omega}}{\sqrt{2}}\right)-\sinh\left(\frac{\sqrt{\omega}}{\sqrt{2}}\right)}{4\sqrt{2}\omega^{5/2}\left(\cos\left(\frac{\sqrt{\omega}}{\sqrt{2}}\right)-\cosh\left(\frac{\sqrt{\omega}}{\sqrt{2}}\right)\right)}\right)\right)\partial_{x}^{2}\bar{T}.$
(41)

which is identical to the result of standard homogenization theory [29].
Equation (41) is simpler and performs as well as equation (40) at sufficiently
large time scales. Of course, at intermediate times scales or in the case with
irregular fluctuating flows, equation (40) performs better.

Figure 3 shows the comparison of the numerical solution and various different
approximations at diffusion time scale $t=1$. The left column shows the result
for a small frequency, $\omega=\pi/5$. The cross-sectional average of the
numerical solution, the solution of effective equations (38) and (40) are
almost indistinguishable. Recall that the standard homogenization result (41)
requires $t\gg\mathcal{O}(\frac{1}{\omega})$. As we expected, the standard
homogenization result on this timescale is substantially worse than both
center manifold results. Alternatively, at higher frequency, with
$\omega=20\pi$, (38) performs visibly worse than both standard homogenization
(41) as well as the time-dependent center manifold results (40). These
observations from the numerical simulation are consistent with our previous
theoretical analysis.

Figure 3: First row: Numerical solution of equation (3) at $t=1$ with the
shear flow $u(y,t)=\left(\cos\omega t\right){y(1-y)}/{2}$, $\mathrm{Pe}=200$
and initial condition
$T_{I}=\left(\sqrt{2\pi}\sigma\right)^{-1}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right),\sigma=1/40$,
where $\omega=\pi/5$ in left panel, and $\omega=20\pi$ in right column. Second
row: red curve is the cross sectional average of the numerical solution. The
blue dash curve is the solution of equation (38). The black dot dash curve is
the solution of equation (40). The purple dot curve is the solution of
equation (41).

##  4  Time varying random flows

Most studies of Taylor dispersion focused on periodic time varying flows,
fewer studies have addressed irregularly fluctuating flows and even random
flows. In this section, we will show that the theory we developed in the
previous section can be applied to time varying random flows. Moreover, for
random flows involving a white noise process or renewing processes, we show
that the effective diffusivity is deterministic at long times.

This is also inspired by our work [30] which studied the advection-diffusion
equation with the shear flow $(v(y,\xi(t)),0)$ where $\xi(t)$ is a stationary
Ornstein-Uhlenbeck (OU) process in parallel-plate channels enforcing the no-
flux boundary conditions. We derived the effective equation at long times via
analyzing the $N$-point correlation function of the random scalar field. The
analysis shows an interesting result that, in this random system, the
effective diffusivity is deterministic at long times.

First, we consider the case that $\xi(t)$ is a Gaussian white noise process
which is a zero-mean, Gaussian random process whose correlation function is
given by $\left\langle\xi(t)\xi(s)\right\rangle=\delta(t-s)$. The center
manifold approach is clearly valid for a smoothly varying velocity field. As
for the Gaussian white noise which is non-differentiable, we can consider a
sequence of function which converges to the white noise process. The Wong-
Zakai theorem states [66, 31, 36] that the convergence of a process to white
noise process yields the convergence of the systems driven by them. That
justifies the application of the center manifold approach in the non-
differentiable case involving white noise.

By utilizing the ergodicity of the white noise process and equation (23), we
obtain the effective diffusivity

$\displaystyle\kappa_{\mathrm{eff}}$
$\displaystyle=1+\mathrm{Pe}^{2}\left(\frac{1}{\left|\Omega\right|}\int\limits_{\Omega}u^{2}(\mathbf{y})\mathrm{d}\mathbf{y}-\left(\frac{1}{|\Omega|}\int\limits_{\Omega}u(\mathbf{y})\mathrm{d}\mathbf{y}\right)^{2}\right).$
(42)

Equation (42) is identical to equation (18) in [30] which is derived via the
analysis of $N$-point correlation function and Hausdorff moment problem. For
the system with the random flows, in general, one has to repeat the experiment
with different realizations of the flows to obtain the properties of the
passive scalar via ensemble average. However, the deterministic diffusivity
presented in equation (42) implies that one need only observe a single
realization of the passive scalar to access some measurable quantities.

Second, we switch our attention to a class of stochastic flows with finite
correlation time. Consider a shear flow takes the form
$(A(t)\xi(t)u(\mathbf{y}),0)$, where $\xi(t)$ is periodic function with a base
frequency $\omega$, or equivalent, a period $L_{t}=\frac{2\pi}{\omega}$.
$A(t)$ is a piecewise-constant zero-mean random function of time,

$\displaystyle A(t)=A_{n},\quad nL_{t}\leq t<(n+1)L_{t},\;n\in\mathbb{Z},$
(43)

where $A_{n}$ is an independent and identically distributed random variable
with zero mean and finite variance. This type of flow is in the class of
renewing (renovating, innovation) flows, that is, flows that decorrelate
completely in a finite time, taken here to be the period $L_{t}$. Therefore,
it is a good approximation to a stationary process with a finite correlation
time. It has wide applications in the study of the dynamo [13, 69] as well as
in study of the intermittency of passive-scalar decay[60, 37]. For this type
of flow, the closed evolution equation for the statistical moment is unknown.
Hence, the Hausdorff moment problem approach proposed in [30] for rigorously
studying the white noise flow case doesn’t apply to this case. However, we
could apply the center manifold approach to near rigorously derive the
effective equation at long times.

In this case, the time averaged diffusion coefficients is

$\displaystyle\kappa_{\mathrm{eff}}$
$\displaystyle=1+\lim\limits_{t\rightarrow\infty}\frac{\mathrm{Pe}^{2}}{t}\sum\limits_{n=1}^{\infty}\left\langle
u,\phi_{n}\right\rangle^{2}\int\limits_{0}^{t}e^{-\lambda_{n}s}\xi(s)A(s)\int\limits_{0}^{s}e^{\lambda_{n}\tau}\xi(\tau)A(\tau)\mathrm{d}\tau\mathrm{d}s.$
(44)

We can further simplify this formula by utilizing the property of the renewing
process.

To take advantage of the periodicity, we tessellate the integral region by
squares. The double integral in equation (44) becomes

$\displaystyle\sum\limits_{m_{1}=0}^{\lfloor\frac{t}{L_{t}}\rfloor-1}A_{m_{1}}^{2}\int\limits_{m_{1}L_{t}}^{(m_{1}+1)L_{t}}\int\limits_{m_{1}L_{t}}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s+A_{\lfloor\frac{t}{L_{t}}\rfloor}\int\limits_{L_{t}\lfloor\frac{t}{L_{t}}\rfloor}^{t}\int\limits_{0}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)A(\tau)\xi(\tau)\mathrm{d}\tau\mathrm{d}s$
(45)
$\displaystyle+\sum\limits_{m_{1}=1}^{\lfloor\frac{t}{L_{t}}\rfloor-1}\sum\limits_{m_{2}=0}^{m_{1}}A_{m_{1}}A_{m_{2}}\int\limits_{m_{1}L_{t}}^{(m_{1}+1)L_{t}}\int\limits_{m_{2}L_{t}}^{(m_{2}+1)L_{t}}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s$

In fact, only the first term in equation (45) grows linearly on time. Thus it
has the dominant contribution at long times. To demonstrate this point, we
will show that the second and third term are bounded in time. The second term
is an integral over a bounded interval
$s\in[L_{t}\lfloor\frac{t}{L_{t}}\rfloor,t]$. It is enough to show the
integrand is a bounded function of $s$ on this interval. We have

$\displaystyle\int\limits_{0}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)A(\tau)\xi(\tau)\mathrm{d}\tau$
(46) $\displaystyle=$ $\displaystyle
A_{\lfloor\frac{t\omega}{2\pi}\rfloor}\int\limits_{\frac{2\pi}{\omega}\lfloor\frac{t\omega}{2\pi}\rfloor}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau+\sum\limits_{m_{2}=0}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}A_{m_{2}}\int\limits_{\frac{2\pi
m_{2}}{\omega}}^{\frac{2\pi(m_{2}+1)}{\omega}}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau$
$\displaystyle=$ $\displaystyle
A_{\lfloor\frac{t\omega}{2\pi}\rfloor}\int\limits_{\frac{2\pi}{\omega}\lfloor\frac{t\omega}{2\pi}\rfloor}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau+\sum\limits_{m_{2}=0}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}A_{m_{2}}e^{-\lambda_{n}\left(s-\frac{2\pi
m_{2}}{\omega}\right)}\int\limits_{0}^{\frac{2\pi}{\omega}}e^{\lambda_{n}\tau}\xi(s)\xi(\tau)\mathrm{d}\tau,$

where both terms in the last step are bounded functions of $s$. Next, we
consider the third term in equation (45). With rearranging the order of the
double summation, we have

$\displaystyle\sum\limits_{m_{1}=1}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}\sum\limits_{m_{2}=0}^{m_{1}}A_{m_{1}}A_{m_{2}}e^{-\lambda_{n}\frac{2\pi(m_{1}-m_{2})}{\omega}}\int\limits_{0}^{\frac{2\pi}{\omega}}\int\limits_{0}^{\frac{2\pi}{\omega}}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s$
(47) $\displaystyle=$
$\displaystyle\sum\limits_{q=1}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}\left(e^{-\lambda_{n}\frac{2\pi
q}{\omega}}\int\limits_{0}^{\frac{2\pi}{\omega}}\int\limits_{0}^{\frac{2\pi}{\omega}}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s\sum\limits_{m_{1}=q}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}A_{m_{1}}A_{m_{1}-q}\right)$

For a fixed $q$,
$\sum\limits_{m_{1}=q}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}A_{m_{1}}A_{m_{1}-q}\rightarrow
E(A_{q}A_{0})+\mathcal{O}(1)=\mathcal{O}(1)$ almost surely because of the law
of large numbers. In addition, the summand decays rapidly as $q$ increases.
Therefore, this summation is also bounded in time almost surely.

Now, we have the leading order approximation of equation (44) at long times,

$\displaystyle\kappa_{\mathrm{eff}}$
$\displaystyle=1+\lim\limits_{t\rightarrow\infty}\frac{\mathrm{Pe}^{2}}{t}\sum\limits_{n=1}^{\infty}\left\langle
u,\phi_{n}\right\rangle^{2}\int\limits_{0}^{\frac{2\pi}{\omega}}\int\limits_{0}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s\sum\limits_{m_{1}=0}^{\lfloor\frac{t\omega}{2\pi}\rfloor-1}A_{m_{1}}^{2}+\mathcal{O}(t^{-1})$
(48)
$\displaystyle=1+\frac{\mathrm{Pe}^{2}\omega}{2\pi}\mathrm{Var}(A_{0})\sum\limits_{n=1}^{\infty}\left\langle
u,\phi_{n}\right\rangle^{2}\int\limits_{0}^{\frac{2\pi}{\omega}}\int\limits_{0}^{s}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s.$

where the second step follows the law of large numbers.

It is natural to compare the renewing flow
$(A(t)\xi(t)u(\mathbf{y}),\mathbf{0})$ with its deterministic counterpart
$(\mathrm{var}(A_{0})\xi(t)u(\mathbf{y}),\mathbf{0})$, and ask the question
which one induces a larger effective diffusivity. One may expect the random
motion creates a larger dispersion. However, it is not always true. A counter
example is the $\xi(t)=\cos t$, $\mathrm{Pe}=1$ and $\mathrm{Var}(A)=1$, where
the effective diffusivity induced by the renewing flow is
$\kappa_{\mathrm{eff,r}}\approx 1.3993$, while the effective diffusivity
induced by its deterministic counterpart is $\kappa_{\mathrm{eff,d}}\approx
1.4124$.

Interestingly, if we only consider the continuous renewing flow, then we have
$\kappa_{\mathrm{eff,r}}\geq\kappa_{\mathrm{eff,d}}$. The continuity of
$A(t)\xi(t)$ implies $\xi(0)=\xi(L_{t})=0$. Hence, $\xi(t)$ admits a sine
expansion $\xi(t)=\sum\limits_{k=1}^{\infty}c_{k}\sin k\omega t$. Since the
only difference in the effective diffusivity formula between the periodic in
time case and the renewing process is the third integral in equation (45), it
is enough to establish that this third term is non-positive. We have

$\displaystyle\int\limits_{m_{1}L_{t}}^{(m_{1}+1)L_{t}}\int\limits_{m_{2}L_{t}}^{(m_{2}+1)L_{t}}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s$
(49) $\displaystyle=$ $\displaystyle
e^{-\lambda_{n}L_{t}(m_{1}-m_{2})}\int\limits_{0}^{\frac{2\pi}{\omega}}\int\limits_{0}^{\frac{2\pi}{\omega}}e^{-\lambda_{n}(s-\tau)}\xi(s)\xi(\tau)\mathrm{d}\tau\mathrm{d}s$
$\displaystyle=$ $\displaystyle
e^{-\lambda_{n}L_{t}(m_{1}-m_{2})}\sum\limits_{k_{1},k_{2=1}}^{\infty}c_{k_{1}}c_{k_{2}}C_{k_{1},k_{2}}.$

where
$C_{k_{1},k_{2}}=\int\limits_{0}^{\frac{2\pi}{\omega}}\int\limits_{0}^{\frac{2\pi}{\omega}}e^{-\lambda_{n}(s-\tau)}\sin(k_{2}\omega
s)\sin(k_{1}\omega\tau)\mathrm{d}\tau\mathrm{d}s$. It is enough to show the
(possibly infinite) matrix $C$ is semi-negative definite. In fact, we have

$\displaystyle C_{k_{1},k_{2}}=$
$\displaystyle-4\omega^{2}\sinh^{2}\left(\frac{\pi^{3}n^{2}}{\omega}\right)\frac{k_{1}}{k_{1}^{2}\omega^{2}+\lambda_{n}^{2}}\frac{k_{2}}{k_{2}^{2}\omega^{2}+\lambda_{n}^{2}}.$
(50)

Hence, for any $n$, $C$ is a rank one matrix with only one negative
eigenvalue, which implies $C$ is semi-negative definite. Now, we finished the
proof of $\kappa_{\mathrm{eff,r}}\geq\kappa_{\mathrm{eff,d}}$ for the
continuous renewing flow.

To verify our theoretical results regarding the deterministic effective
diffusivity, we solve equation (3) with different shear flows by using the
forward Monte-Carlo method described in [29]. The computational domain is
$(x,y)\in\mathbf{R}\times[0,1]$. The time step size is $10^{-3}$. The total
number of the random walkers is $2\times 10^{6}$. We divide a simulation into
$400$ parallel jobs on UNC’s Longleaf computing cluster. The shear flow takes
the form $v(y,t)=A(t)\sin t(y-1/2)$. In panel (a), $A(t)$ is a white noise
process. In panel (b), $A(t)=1$. In panel (c, d, e), $A(t)$ is a renewing
process with a coin-toss random variable taking values plus or minus one with
equal probability, a uniform distributed random variable on
$[-\sqrt{3},\sqrt{3}]$ and a standard Gaussian distributed random variable
respectively. We plot $\frac{\mathrm{Var}(\bar{T})}{2t}$ as a function of time
for $5$ independent flow realizations and different shear flows in figure 4.
The curves with the same color are generated with the same seed from the same
random number generator.

From figure 4, we have four1 observations. First, in panel (a,c,d,e), all
curves fluctuate randomly at the earlier stage but converge at later times to
a deterministic effective diffusivity $\kappa_{\mathrm{eff}}$ given by
equation (48). Second, since all distributions in panel (c, d, e) have the
same unit variance, all renewing flows induce the same effective diffusivity
at long times. Third, comparing panel (b) and panel (c, d, e), we can see that
renewing random flows induce a larger effective diffusivity than their
deterministic counterpart, as just proven above. Fourth, from the right column
of figure 4, we can see that if the distribution of $A(t)$ has a heavier tail,
then $\frac{\mathrm{Var}(\bar{T})}{2t}$ takes a longer time to converge to the
theoretical limit.

Figure 4: $\frac{\mathrm{Var}(\bar{T})}{2t}$ as a function of time for $5$
independent flow realizations and different random flows. Note the convergence
of this quantity to the deterministic effective diffusivity given in equation
(6). We compute $\kappa_{\mathrm{eff}}$ by equation (42) for panel (a), by
equation (23) for panel (b), by equation (48) for panel (c, d, e) and report
three significant digits of the effective diffusivity to the right of each
panel. Pictures in the right column are simply zoom-in of pictures in the left
column at a larger time scale.

###  4.1 Invariant measure

Equation (24) is an approximation of the scalar field at long times, which is
a powerful tool to compute the invariance measure of the random field. When
$v(y,t)=\xi(t)u(\mathbf{y})$ and $\xi(t)$ is the Gaussian white noise,
equation (24) becomes

$\displaystyle\bar{T}(x,t)=\frac{1}{\sqrt{4\pi\kappa_{\mathrm{eff}}}}\exp\left(\frac{-\tilde{x}^{2}}{4\kappa_{\mathrm{eff}}t}\right)+\mathcal{O}(t^{-\frac{3}{2}}),\quad\tilde{x}=x-\mathrm{Pe}\bar{u}B(t),$
(51)
$\displaystyle\kappa_{\mathrm{eff}}=1+\mathrm{Pe}^{2}\left(\frac{1}{\left|\Omega\right|}\int\limits_{\Omega}u^{2}(\mathbf{y})\mathrm{d}\mathbf{y}-\left(\frac{1}{|\Omega|}\int\limits_{\Omega}u(\mathbf{y})\mathrm{d}\mathbf{y}\right)^{2}\right).$

where $B(t)$ is the standard Brownian motion. Then we apply the inverse
transform method (we refer reader to [17] for details) to obtain the invariant
measure of $\bar{T}$, i.e, the probability density function at long times,
from the probability density function of $B(s)$. We consider the rescaling of
$T$, $\tilde{T}(x,y,t)=\sqrt{4\pi\kappa_{\mathrm{eff}}t}T$. Without loss of
generality, we focus on the scalar at point $x=0,y=0$, i.e.,
$\tilde{T}(0,0,t)$. Thus, the invariant measure is

$\displaystyle
f_{\tilde{T}}(z)=\frac{z^{\frac{1}{\beta}-1}}{\sqrt{-\pi\beta\log(z)}},\quad
z\in[0,1],$ (52)

where
$\beta=\frac{\mathrm{Pe}^{2}\bar{u}^{2}v(t)}{2t\kappa_{\text{eff}}}=\frac{\mathrm{Pe}^{2}\bar{u}^{2}}{2\kappa_{\text{eff}}}+\mathcal{O}(t^{-1})$
and $v(t)$ is the variance of $\int\limits_{0}^{t}\xi(s)\mathrm{d}s$.
$f_{\tilde{T}}(z)$ always has the logarithmic singularity at $z=1$. It is
continuous at $z=0$ when $\beta\leq 1$, and singular when $\beta>1$ (see
figure 5). This property of the distribution implies that when the strength of
the input random signal exceeds some certain threshold, the value of scalar is
more likely to be zero. As a result of that, the distribution changes from
negatively-skewed to positively-skewed as $\beta$ increases.

Figure 5:  The invariant measure $f_{\tilde{T}}(z)$ in equation (52) for
different parameters $\beta$. The red solid curve, blue dashed curve and black
dot dash curve stands for the case $\beta=1/20,1,20$, respectively.
$f_{\tilde{T}}(z)$ changes from negatively-skewed to positively-skewed as
$\beta$ increases.

##  5 Conclusion and Discussion

We have studied the long time behavior of an advection-diffusion equation with
a general time-varying (including random) shear flow imposing no-flux boundary
conditions on channel walls using center manifold theory. Our study extends
and improves prior work of [45, 46, 42] to properly incorporate general time
variation into the effective theory. Comparisons with full simulations
document conditions when this improved approach gives a better approximation,
and also illustrates situations in which standard homogenization does not
perform on finite timescales. Convergence studies illustrate how the accuracy
of the different approximations. Armed with this improved time varying center
manifold theory, we derived new effective equations for random shear flows
involving both white in time statistics, as well as more correlated renewing
flows. For white in time, these predictions agree with our prior work [30],
which forecast a deterministic effective diffusivity on long times. For the
case of renewing flows, less is known, and our current work also a
deterministic effective diffusivity, with new explicit formulae. These
theories are demonstrated to be quantitatively accurate through Monte-Carlo
simulations. New conditions are derived which guarantee when the random
renewing flow generates a larger effective diffusivity than its deterministic
analog. Lastly, using inverse transform method and the effective equations, we
derived the invariant measure and study its Péclet number dependence.

In this study, we only considered constant diffusivity. Future immediate areas
of exploration include case with spatial variable dependent diffusivity or
even concentration dependent diffusivity. A practical example concerns the
shear-enhanced diffusion in colloidal suspensions explored in [34]. The
nonlinearity in those system imposes challenges to the traditional method. We
expect center manifold theory could overcome the difficulties. Further, center
manifold theory will apply nicely to study the mixing ability of time-varying
flow in a non-flat channel to generalize the conclusion in
[rosencrans1997taylor].

##  6 Acknowledgements

We acknowledge funding received from the following NSF Grant Nos.:DMS-1910824;
and ONR Grant No: ONR N00014-18-1-2490. Partial support for Lingyun Ding is
gratefully acknowledged from the National Science Foundation, award NSF-
DMS-1929298 from the Statistical and Applied Mathematical Sciences Institute.

##  7 Appendix

###  7.1 Numerical Method

In this section, we document details of the algorithm for the numerical
simulation of equation (1). The computational domain is $x\times
y\in[-H,H]\times[0,L]$. When $H$ is large enough, we can assume there is a
periodic boundary condition in the $x$-direction. Since there are non-
penetration conditions in the $y$-direction, we perform the even extension in
the $y$-direction to obtain the periodic condition on the extended domain.
Thus, we can use the standard Fourier spectral method to solve the advection-
diffusion equation with periodic boundary conditions on the rectangular domain
$[-H,H]\times[0,2L]$. In the dealiasing process at each time step, we apply
the all-or-nothing filter with the two-thirds rule to the spectrum, that is,
we set the upper one-third of the resolved spectrum to zero (see chapter 11 of
the book [16] for details).

The diffusion operator is stiff, which requires a very small time step size
for the explicit method to ensure numerical stability. In order to use a
larger time step size and improve the efficiency, we adopt the implicit-
explicit third-order Runge-Kutta method presented in table 6 in [48]. In our
application, we use the explicit Runge-Kutta method to integrate the advection
part and use the implicit diagonal Runge-Kutta method to integrate the
diffusion term. When the diffusivity is a constant, the diffusion operator is
a diagonal matrix in the Fourier space. Thus, the implicit equation can be
solved explicitly and efficiently. The implicit-explicit method is as
efficient as the explicit method at each iteration while allows a much larger
time step size. When the diffusivity is a function of spatial variables, the
implicit part requires solving a linear system, which is expansive. Therefore,
in this case or when the advection is dominant, we adopt the explicit 4th-
order Runge-Kutta method as the time-marching scheme.

We also present the Butcher tableau of the explicit-implicit Runge-Kutta
method in table 1 here for convenience. Unfortunately, [48] only reported
13-14 significant digits of parameters $(\alpha,\beta,\eta)$ which are the key
parameters defining the algorithm. That may potentially deteriorate the
accuracy of double-precision floating-point based or even higher precision
floating-point based algorithms. Hence, we documented the exact value for
those parameters,
$(\frac{9-\sqrt{57}}{6},\frac{9-\sqrt{57}}{24},\frac{-6+\sqrt{57}}{12})$. We
also find another two groups of parameters that achieve the same convergence
order and ensure the L-stable, $(1/2,1/8,0)$ and
$(\frac{9+\sqrt{57}}{6},\frac{9+\sqrt{57}}{24},\frac{-6-\sqrt{57}}{12})$.

0 | 0 | 0 | 0 | 0
---|---|---|---|---
0 | 0 | 0 | 0 | 0
1 | 0 | 1 | 0 | 0
1/2 | 0 | 1/4 | 1/4 | 0
| 0 | 1/6 | 1/6 | 2/3

$\alpha$ | $\alpha$ | 0 | 0 | 0
---|---|---|---|---
0 | -$\alpha$ | $\alpha$ | 0 | 0
1 | 0 | 1-$\alpha$ | $\alpha$ | 0
1/2 | $\beta$ | $\eta$ | $1/2-\beta-\eta-\alpha$ | $\alpha$
| 0 | 1/6 | 1/6 | 2/3

Table 1: Butcher tableau for the Explicit (left) Implicit (right) L-Stable
scheme, $(\alpha,\beta,\eta)$ could be $(1/2,1/8,0)$,
$(\frac{9-\sqrt{57}}{6},\frac{9-\sqrt{57}}{24},\frac{-6+\sqrt{57}}{12})$ or
$(\frac{9+\sqrt{57}}{6},\frac{9+\sqrt{57}}{24},\frac{-6-\sqrt{57}}{12})$.

###  7.2 Lists of abbreviations

See table 2.

Full Form | Abbreviation
---|---
Ornstein-Uhlenbeck | OU
Partial differential equation | PDE
Probability distribution function | PDF
Stochastic differential equation | SDE
Table 2: Lists of abbreviations.

## References

  * Aminian et al. [2016] Aminian, M., Bernardi, F., Camassa, R., Harris, D.M., McLaughlin, R.M., 2016.  How boundaries shape chemical delivery in microfluidics.  Science 354, 1252–1256.
  * Aminian et al. [2015] Aminian, M., Bernardi, F., Camassa, R., McLaughlin, R.M., 2015\.  Squaring the circle: Geometric skewness and symmetry breaking for passive scalar transport in ducts and pipes.  Physical review letters 115, 154503\.
  * Aminian et al. [2018] Aminian, M., Camassa, R., McLaughlin, R.M., 2018.  Mass distribution and skewness for passive scalar transport in pipes with polygonal and smooth cross sections.  Studies in Applied Mathematics 141, 399–417.
  * Aris [1956] Aris, R., 1956.  On the dispersion of a solute in a fluid flowing through a tube.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 235, 67–77.
  * Aris [1960] Aris, R., 1960.  On the dispersion of a solute in pulsating flow through a tube.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 259, 370–376.
  * Aulbach [1982] Aulbach, B., 1982.  A reduction principle for nonautonomous differential equations.  Archiv der Mathematik 39, 217–232.
  * Aulbach and Wanner [1996] Aulbach, B., Wanner, T., 1996\.  Integral manifolds for carathéodory type differential equations in banach spaces.  Six lectures on dynamical systems 2\.
  * Aulbach and Wanner [1999] Aulbach, B., Wanner, T., 1999\.  Invariant foliations for carathéodory type differential equations in banach spaces.  Advances of Stability Theory at the End of XX Century?, Gordon & Breach Publishers. http://citeseerx. ist. psu. edu/viewdoc/download .
  * Balakotaiah et al. [1995] Balakotaiah, V., Chang, H.C., Smith, F., 1995.  Dispersion of chemical solutes in chromatographs and reactors.  Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences 351, 39–75.
  * Bandyopadhyay and Mazumder [1999] Bandyopadhyay, S., Mazumder, B., 1999\.  On contaminant dispersion in unsteady generalised couette flow.  International journal of engineering science 37, 1407–1423.
  * Barik and Dalal [2017] Barik, S., Dalal, D., 2017.  On transport coefficients in an oscillatory couette flow with nonlinear chemical decay reactions.  Acta Mechanica 228, 2391–2412.
  * Barik and Dalal [2019] Barik, S., Dalal, D., 2019.  Multi-scale analysis for concentration distribution in an oscillatory couette flow.  Proceedings of the Royal Society A 475, 20180483.
  * Bhat and Subramanian [2015] Bhat, P., Subramanian, K., 2015\.  Fluctuation dynamos at finite correlation times using renewing flows.  Journal of Plasma Physics 81.
  * Blom et al. [2003] Blom, M.T., Chmela, E., Oosterbroek, R.E., Tijssen, R., Van Den Berg, A., 2003.  On-chip hydrodynamic chromatography separation and detection of nanoparticles and biomolecules.  Analytical chemistry 75, 6761–6768.
  * Bowden [1965] Bowden, K., 1965.  Horizontal mixing in the sea due to a shearing current.  Journal of Fluid Mechanics 21, 83–95.
  * Boyd [2001] Boyd, J.P., 2001.  Chebyshev and Fourier spectral methods.  Courier Corporation.
  * Bronski et al. [2007] Bronski, J.C., Camassa, R., Lin, Z., McLaughlin, R.M., Scotti, A., 2007.  An explicit family of probability measures for passive scalar diffusion in a random flow.  Journal of Statistical Physics 128, 927–968.
  * Bronski and McLaughlin [1997] Bronski, J.C., McLaughlin, R.M., 1997\.  Scalar intermittency and the ground state of periodic schrödinger equations.  Physics of fluids 9, 181–190.
  * Camassa et al. [2021] Camassa, R., Ding, L., Kilic, Z., McLaughlin, R.M., 2021\.  Persisting asymmetry in the probability distribution function for a random advection-diffusion equation in impermeable channels.  Physica D: Nonlinear Phenomena , 132930\.
  * Camassa et al. [2010] Camassa, R., Lin, Z., McLaughlin, R.M., 2010.  The exact evolution of the scalar variance in pipe and channel flow.  Communications in Mathematical Sciences 8, 601–626.
  * Carr [2012] Carr, J., 2012.  Applications of centre manifold theory. volume 35.  Springer Science & Business Media.
  * Carr and Muncaster [1983a] Carr, J., Muncaster, R.G., 1983a.  The application of centre manifolds to amplitude expansions. i. ordinary differential equations.  Journal of differential equations 50, 260–279.
  * Carr and Muncaster [1983b] Carr, J., Muncaster, R.G., 1983b.  The application of centre manifolds to amplitude expansions. ii. infinite dimensional problems.  Journal of differential equations 50, 280–288.
  * Chatwin [1970] Chatwin, P., 1970.  The approach to normality of the concentration distribution of a solute in a solvent flowing along a straight pipe.  Journal of Fluid Mechanics 43, 321–352.
  * Chatwin [1975] Chatwin, P., 1975.  On the longitudinal dispersion of passive contaminant in oscillatory flows in tubes.  Journal of Fluid Mechanics 71, 513–527.
  * Chatwin and Allen [1985] Chatwin, P., Allen, C., 1985\.  Mathematical models of dispersion in rivers and estuaries.  Annual Review of Fluid Mechanics 17, 119–149.
  * Chu et al. [2019] Chu, H.C., Garoff, S., Przybycien, T.M., Tilton, R.D., Khair, A.S., 2019.  Dispersion in steady and time-oscillatory two-dimensional flows through a parallel-plate channel.  Physics of Fluids 31, 022007\.
  * Chu et al. [2020] Chu, H.C., Garoff, S., Tilton, R.D., Khair, A.S., 2020\.  Dispersion in steady and time-oscillatory flows through an eccentric annulus.  AIChE Journal 66, e16831.
  * Ding et al. [2021] Ding, L., Hunt, R., McLaughlin, R.M., Woodie, H., 2021\.  Enhanced diffusivity and skewness of a diffusing tracer in the presence of an oscillating wall.  Research in the Mathematical Sciences 8, 1–29.
  * Ding and McLaughlin [2020] Ding, L., McLaughlin, R.M., 2020\.  Ergodicity and invariant measures for a diffusing passive scalar advected by a random channel shear flow and the connection between the kraichnan-majda model and taylor-aris dispersion.  arXiv preprint arXiv:2012.06610 .
  * Eugene and Moshe [1965] Eugene, W., Moshe, Z., 1965\.  On the relation between ordinary and stochastic differential equations.  International Journal of Engineering Science 3, 213–229.
  * Gill [1967] Gill, W., 1967.  A note on the solution of transient dispersion problems.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 298, 335–339.
  * Gill and Sankarasubramanian [1970] Gill, W., Sankarasubramanian, R., 1970\.  Exact analysis of unsteady convective diffusion.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 316, 341–350.
  * Griffiths and Stone [2012] Griffiths, I., Stone, H.A., 2012\.  Axial dispersion via shear-enhanced diffusion in colloidal suspensions.  EPL (Europhysics Letters) 97, 58005\.
  * Guggenheim [1954] Guggenheim, E., 1954.  The diffusion coefficient of sodium chloride.  Transactions of the Faraday society 50, 1048–1051.
  * Hairer and Pardoux [2015] Hairer, M., Pardoux, É., 2015\.  A wong-zakai theorem for stochastic pdes.  Journal of the Mathematical Society of Japan 67, 1551–1604.
  * Haynes and Vanneste [2005] Haynes, P.H., Vanneste, J., 2005\.  What controls the decay of passive scalars in smooth flows?  Physics of Fluids 17, 097103\.
  * Jimenez and Sullivan [1984] Jimenez, C., Sullivan, P., 1984\.  Contaminant dispersion in some time-dependent laminar flows.  Journal of Fluid Mechanics 142, 57–77.
  * Kraichnan [1968] Kraichnan, R.H., 1968.  Small-scale structure of a scalar field convected by turbulence.  The Physics of Fluids 11, 945–953.
  * Lee et al. [2021] Lee, G., Luner, A., Marzuola, J., Harris, D.M., 2021\.  Dispersion control in pressure-driven flow through bowed rectangular microchannels.  Microfluidics and Nanofluidics 25, 1–11.
  * Majda [1993] Majda, A.J., 1993.  The random uniform shear layer: an explicit example of turbulent diffusion with broad tail probability distributions.  Physics of Fluids A: Fluid Dynamics 5, 1963–1970.
  * Marbach and Alim [2019] Marbach, S., Alim, K., 2019\.  Active control of dispersion within a channel with flow and pulsating walls.  Physical Review Fluids 4, 114202\.
  * Masri et al. [2019] Masri, R., Puelz, C., Riviere, B., 2019.  A reduced model for solute transport in compliant blood vessels with arbitrary axial velocity profile.  arXiv preprint arXiv:1912.09587 .
  * McLaughlin and Majda [1996] McLaughlin, R.M., Majda, A.J., 1996\.  An explicit example with non-gaussian probability distribution for nontrivial scalar mean and fluctuation.  Physics of Fluids 8, 536–547.
  * Mercer and Roberts [1990] Mercer, G., Roberts, A., 1990\.  A centre manifold description of contaminant dispersion in channels with varying flow properties.  SIAM Journal on Applied Mathematics 50, 1547–1565.
  * Mercer and Roberts [1994] Mercer, G., Roberts, A., 1994\.  A complete model of shear dispersion in pipes.  Japan journal of industrial and applied mathematics 11, 499–521.
  * Mukherjee and Mazumder [1988] Mukherjee, A., Mazumder, B., 1988\.  Dispersion of contaminant in oscillatory flows.  Acta mechanica 74, 107–122.
  * Pareschi and Russo [2005] Pareschi, L., Russo, G., 2005\.  Implicit–explicit runge–kutta schemes and applications to hyperbolic systems with relaxation.  Journal of Scientific computing 25, 129–155.
  * Paul and Mazumder [2008] Paul, S., Mazumder, B., 2008\.  Dispersion in unsteady couette–poiseuille flows.  International journal of engineering science 46, 1203–1217.
  * Roberts [1993] Roberts, A., 1993.  The invariant manifold of beam deformations.  Journal of elasticity 30, 1–54.
  * Roberts [1996] Roberts, A., 1996.  Low-dimensional models of thin film fluid dynamics.  Physics Letters A 212, 63–71.
  * Roberts [2018] Roberts, A., 2018.  Backwards theory supports modelling via invariant manifolds for non-autonomous dynamical systems.  arXiv preprint arXiv:1804.06998 .
  * Roberts and Li [2006] Roberts, A., Li, Z., 2006.  An accurate and comprehensive model of thin fluid flows with inertia on curved substrates.  Journal of Fluid Mechanics 553, 33–73.
  * Salerno et al. [2020] Salerno, L., Cardillo, G., Camporeale, C., 2020.  Aris-taylor dispersion in the subarachnoid space.  Physical Review Fluids 5, 043102\.
  * SMITH [1982] SMITH, R., 1982.  Gaussian approximation for contaminant dispersion.  The Quarterly Journal of Mechanics and Applied Mathematics 35, 345–366.
  * Smith [1987] Smith, R., 1987.  Diffusion in shear flows made easy: the taylor limit.  Journal of Fluid Mechanics 175, 201–214.
  * Taylor [1953] Taylor, G.I., 1953.  Dispersion of soluble matter in solvent flowing slowly through a tube.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 219, 186–203.
  * Trojanowicz and Kołacińska [2016] Trojanowicz, M., Kołacińska, K., 2016\.  Recent advances in flow injection analysis.  Analyst 141, 2085–2139.
  * Van Minh [1993] Van Minh, N., 1993.  A reduction principle for topological classification of nonautonomous differential equations.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics 123, 621–632.
  * Vanneste [2006] Vanneste, J., 2006.  Intermittency of passive-scalar decay: Strange eigenmodes in random shear flows.  Physics of fluids 18, 087108\.
  * Vedel and Bruus [2012] Vedel, S., Bruus, H., 2012.  Transient taylor–aris dispersion for time-dependent flows in straight channels.  Journal of fluid mechanics 691, 95–122.
  * Vedel et al. [2014] Vedel, S., Hovad, E., Bruus, H., 2014.  Time-dependent taylor–aris dispersion of an initial point concentration.  Journal of fluid mechanics 752, 107–122.
  * Watson [1983] Watson, E., 1983.  Diffusion in oscillatory pipe flow.  Journal of Fluid Mechanics 133, 233–244.
  * Watt and Roberts [1996] Watt, S., Roberts, A., 1996\.  The construction of zonal models of dispersion in channels via matched centre manifolds.  The ANZIAM Journal 38, 101–125.
  * Watt and Roberts [1995] Watt, S.D., Roberts, A.J., 1995\.  The accurate dynamic modelling of contaminant dispersion in channels.  SIAM Journal on Applied Mathematics 55, 1016–1038.
  * Wong and Zakai [1965] Wong, E., Zakai, M., 1965.  On the convergence of ordinary integrals to stochastic integrals.  The Annals of Mathematical Statistics 36, 1560–1564.
  * Wu et al. [2016] Wu, Z., Fu, X., Wang, G., 2016\.  On spatial pattern of concentration distribution for taylor dispersion process.  Scientific reports 6, 1–13.
  * Young and Jones [1991] Young, W.a., Jones, S., 1991\.  Shear dispersion.  Physics of Fluids A: Fluid Dynamics 3, 1087–1101.
  * Zel’Dovich et al. [1984] Zel’Dovich, Y.B., Ruzmaikin, A., Molchanov, S., Sokoloff, D., 1984\.  Kinematic dynamo problem in a linear velocity field.  Journal of Fluid Mechanics 144, 1–11.