2.1. Non-dimensionalisation

We denote the characteristic length scale of the channel cross-section as $L_{y}$ . We then proceed with the following change of variables for non-dimensionalisation:

(2.4) \begin{equation} \begin{aligned} &L_{y} \tilde {x}=x, \quad L_{y} \tilde {\boldsymbol{y}}=\boldsymbol{y},\quad U \tilde {u} =u, \quad U \tilde {\boldsymbol{v}}=\boldsymbol{v}, \quad \frac {\rho _{0}\nu U}{L_{y}} \tilde {p}=p, \quad \frac {L_{y}^{2}}{\kappa }\tilde {t}=t,\quad c_{c}\tilde {c}=c. \end{aligned} \end{equation}

Equation (2.2) becomes

(2.5) \begin{equation} \begin{aligned} &\partial _{\tilde {t}} \tilde {c}+\frac {U}{L_{y}}\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {c} = \frac {\kappa }{L_{y}^{2}} \Delta \tilde {c}, \quad \left . \frac {1}{L_{y}} \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {c} \right |_{w\textit{all}}=0 ,\end{aligned} \end{equation}

and (2.3) becomes

(2.6) \begin{equation} \begin{aligned} & \frac {U \kappa }{L_{y}^{2}} \partial _{\tilde {t}}\tilde {u}+ \frac {U^{2}}{L_{y}} \tilde {u} \partial _{\tilde {x}}\tilde {u} + \frac {U^{2}}{L_{y}} \tilde {\boldsymbol{v}} \boldsymbol{\cdot }\partial _{\tilde {\boldsymbol{y}}}\tilde {u}= -\frac {\nu U}{L_{y}^{2}} \partial _{\tilde {x}}\tilde {p} + \nu U \left ( \frac {1}{L_{y}^{2}} \partial _{\tilde {x}}^{2}\tilde {u} + \frac {1}{L_{y}^{2}} \Delta _{\tilde {\boldsymbol{y}}}^{2}\tilde {u} \right )\!,\\ & \frac {U \kappa }{L_{y}^{2}} \partial _{\tilde {t}}\tilde {v}_{i}+ \frac {U^{2}}{L_{y}} \tilde {u} \partial _{\tilde {x}}\tilde {v}_{i} + \frac {U^{2}}{L_{y}} \tilde {\boldsymbol{v}}\boldsymbol{\cdot }\partial _{\tilde {\boldsymbol{y}}}\tilde {v}_{i}= -\frac {\nu U}{L_{y}^{2}} \partial _{\tilde {y}_{i}}\tilde {p} + \nu U \left ( \frac {1}{L_{y}^{2}} \partial _{\tilde {x}}^{2}\tilde {v}_{i} + \frac {1}{L_{y}^{2}} \Delta _{\tilde {\boldsymbol{y}}}^{2}\tilde {v}_{i} \right )\!,\\ &\frac {U}{L_{y}}\partial _{\tilde {x}}\tilde {u}+\frac {U}{L_{y}}\boldsymbol{\nabla} _{\tilde {\boldsymbol{y}}} \boldsymbol{\cdot }\tilde {\boldsymbol{v}}=0. \end{aligned} \end{equation}

Now the velocity field $\tilde {\boldsymbol{u}}$ has the dimensionless fundamental period $\tilde {L}_{p} = L_{p} / L_{y}$ . After dropping the tilde, we obtain the equation for the scalar field

(2.7) \begin{equation} \begin{aligned} &\partial _{t} c+ {\textit{Pe}} \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }c = \Delta c ,\quad c (x,\boldsymbol{y},0) = c_{I} \left ( x,\boldsymbol{y}\right )\!,\quad \left . \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }c \right |_{w\textit{all}}=0, \end{aligned} \end{equation}

and the equation for the velocity field

(2.8) \begin{equation} \begin{aligned} &\frac {1}{{Sc}}\partial _{t}u+ {\textit{Re}}\left ( u \partial _{x}u + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}}u \right )= -\partial _{x}p + \big ( \partial _{x}^{2}u + \Delta _{\boldsymbol{y}}u \big ),\\&\frac {1}{{Sc}}\partial _{t}v_{i}+ {\textit{Re}} \left ( u \partial _{x}v_{i} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}}v_{i} \right )= - \partial _{y_{i}}p + \big ( \partial _{x}^{2}v_{i} + \Delta _{\boldsymbol{y}}v_{i} \big ),\\&\partial _{x}u+\boldsymbol{\nabla} _{\boldsymbol{y}} \boldsymbol{\cdot }\boldsymbol{v}=0,\quad \left . \boldsymbol{u} \right |_{w\textit{all}}=0, \end{aligned} \end{equation}

where we have defined the Pélect number ${\textit{Pe}}={L_{y}U}/{\kappa }$ , the Reynolds number ${\textit{Re}}={ L_{y}U}/{\nu }$ and the Schmidt number ${Sc}={\nu }/{ \kappa }$ . Unless otherwise noted, all quantities in the subsequent sections, including $\tilde {L}_{p}$ , are taken to be dimensionless, and for simplicity, we drop the tilde.

2.2. Variance and effective longitudinal diffusivity

Since the solution of the advection–diffusion equation in a channel domain converges to a Gaussian distribution at long times, the variance serves as a key quantity for characterising the solution. In this section we present its definition.

As the following discussion involves complex-valued functions, we define the inner product as

(2.9) \begin{equation} \left \langle f, g \right \rangle = \frac {1}{|\varOmega (L_p)|} \int \limits _{\varOmega } f(x, \boldsymbol{y}) g^{*}(x, \boldsymbol{y}) \,\mathrm{d}x \,\mathrm{d}\boldsymbol{y}, \end{equation}

where the asterisk denotes the complex conjugate. For simplicity, we use the notation $\left \langle u \right \rangle$ as a shorthand for $\left \langle u, 1 \right \rangle$ .

Since the solution approaches a Gaussian function at long times, its variance is given by

(2.10) \begin{equation} \begin{aligned} \text{Var}_{c}(t) &= \frac {\left \langle (x-\left \langle xc \right \rangle )^{2}c \right \rangle }{\left \langle c \right \rangle }=\frac {\left \langle x^{2}c \right \rangle - \left \langle xc \right \rangle ^{2}}{\left \langle c \right \rangle }. \end{aligned} \end{equation}

Then the effective longitudinal diffusivity is defined as

(2.11) \begin{equation} \kappa _{\textit{eff}} = \lim \limits _{t\rightarrow \infty } \frac {\text{Var}_{c}(t)}{2t} = \lim \limits _{t\rightarrow \infty } \frac {1}{2} \partial _{t} \text{Var}_{c}(t). \end{equation}

In some studies on time-dependent flows, the term $({1}/{2}) \partial _{t} \text{Var}(t)$ is directly defined as the effective diffusivity, as seen in Vedel & Bruus (Reference Vedel and Bruus2012), Vedel et al. (Reference Vedel, Hovad and Bruus2014), Ding & McLaughlin (Reference Ding and McLaughlin2023).

2.3. Eigenfunction expansion

Since the advection--diffusion equation (2.2) is linear, the eigenfunction expansion of the scalar field is a convenient tool for analysing scalar dynamics. The eigenvalue $\lambda$ and eigenfunction $\phi$ satisfy the following equation:

(2.12) \begin{equation} \begin{aligned} &\lambda \phi = \mathcal{L} \phi = {\textit{Pe}} \left ( u \partial _{x} \phi + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \phi \right ) - \partial _{x}^{2} \phi - \Delta _{\boldsymbol{y}} \phi , \quad \left . n_{x} \partial _{x} \phi + \partial _{\boldsymbol{n}_{y}} \phi \right |_{w\textit{all}} = 0. \end{aligned} \end{equation}

It forms a basis for square-integrable functions defined on the infinite channel domain that satisfy the no-flux boundary condition $\left . n_{x} \partial _{x} \phi + \partial _{\boldsymbol{n}_{y}} \phi \right |_{w\textit{all}} = 0$ .

There are three properties of eigenvalues and eigenfunctions that will be used later. First, given that the velocity field is real, $\lambda ^{*}$ is also an eigenvalue and $\phi ^{*}$ is the associated eigenfunction. Second, the real part of the eigenvalue is non-negative (see the proof in Appendix A.1). Third, $\lambda = 0$ is the smallest eigenvalue and $\phi = 1$ is the associated eigenfunction.

Since the advection operator is a skew-adjoint operator, $\mathcal{L}$ is not a self-adjoint operator when the velocity field is non-zero. Hence, we introduce $\lambda ^{*}$ and $\varphi (x, \boldsymbol{y})$ , which are the eigenvalues and eigenfunctions of the adjoint eigenvalue problem given by

(2.13) \begin{equation} \begin{aligned} &\lambda ^{*} \varphi = \mathcal{L}^{*} \varphi = -{\textit{Pe}} \left ( u \partial _{x} \varphi + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \varphi \right ) - \partial _{x}^{2} \varphi - \Delta _{\boldsymbol{y}} \varphi , \quad \left . n_{x} \partial _{x} \varphi + \partial _{\boldsymbol{n}_{y}} \varphi \right |_{w\textit{all}} = 0. \end{aligned} \end{equation}

In the adjoint eigenvalue problem, ${\textit{Pe}}$ is replaced with $-{\textit{Pe}}$ , effectively reversing the flow direction.

Next, we choose $\phi _{n}$ and $\varphi _{n}$ such that the sets $\left \{ \phi _{n} \right \}_{n=0}^{\infty }$ and $\left \{ \varphi _{n} \right \}_{n=0}^{\infty }$ form a biorthogonal system, i.e. $\left \langle \phi _{n}, \varphi _{m} \right \rangle = \delta _{n,m}$ , where $\delta _{n,m}$ is the Kronecker delta. The solution of the advection--diffusion equation $\partial _{t} c = \mathcal{L} c$ can be formally represented as

(2.14) \begin{align} c(x, \boldsymbol{y}, t) = \sum _{n} \left \langle c_{I}(x, \boldsymbol{y}), \varphi _{n}(x, \boldsymbol{y}) \right \rangle \phi _{n}(x, \boldsymbol{y}) \text{e}^{-\lambda _{n} t}. \end{align}

The exact solution to this eigenvalue problem is typically unavailable, necessitating the use of numerical methods. However, the infinite domain presents challenges for numerical computations. One significant simplification of the problem is to consider the eigenfunctions in the form

(2.15) \begin{equation} \begin{aligned} &\phi (x, \boldsymbol{y}, k) =\frac {1}{\sqrt {2\pi }} e^{\mathrm{i} k x} \hat {\phi }(x, \boldsymbol{y}, k), \quad \varphi (x, \boldsymbol{y}, k) = \frac {1}{\sqrt {2\pi }} e^{\mathrm{i} k x} \hat {\varphi }(x, \boldsymbol{y}, k), \end{aligned} \end{equation}

where $\hat {\phi }(x, \boldsymbol{y}, k)$ solves the equation

(2.16) \begin{align} &\lambda \hat {\phi } = \kern2pt \hat{\kern-2pt \mathcal{L}}(k) \hat {\phi } = {\textit{Pe}} \big ( \mathrm{i} k u \hat {\phi } + u \partial _{x} \hat {\phi } + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \hat {\phi } \big ) - \big ( \partial _{x}^{2} \hat {\phi } + 2\mathrm{i} k \partial _{x} \hat {\phi } + (\mathrm{i} k)^{2} \hat {\phi } + \Delta _{\boldsymbol{y}} \hat {\phi } \big ),\nonumber\\& \left . n_{x} \mathrm{i} k \hat {\phi } + n_{x} \partial _{x} \hat {\phi } + \partial _{\boldsymbol{n}_{y}} \hat {\phi } \right |_{w\textit{all}} = 0, \quad \hat {\phi }(x, \boldsymbol{y}, k)=\hat {\phi }(x+L_{p}, \boldsymbol{y}, k), \end{align}

and $\hat {\varphi }(x, \boldsymbol{y}, k)$ solves the adjoint equation

(2.17) \begin{align} &\lambda ^{*} \hat {\varphi } = \kern2pt \hat{\kern-2pt \mathcal{L}^{*}}(k) \hat {\varphi } = -{\textit{Pe}} \left ( \mathrm{i} k u \hat {\varphi } + u \partial _{x} \hat {\varphi } + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \hat {\varphi } \right ) - \big ( \partial _{x}^{2} \hat {\varphi } + 2\mathrm{i} k \partial _{x} \hat {\varphi } + (\mathrm{i} k)^{2} \hat {\varphi } + \Delta _{\boldsymbol{y}} \hat {\varphi } \big ),\nonumber \\& \left . n_{x} \mathrm{i} k \hat {\varphi } + n_{x} \partial _{x} \hat {\varphi } + \partial _{\boldsymbol{n}_{y}} \hat {\varphi } \right |_{w\textit{all}} = 0, \quad \hat {\varphi }(x, \boldsymbol{y}, k)=\hat {\varphi }(x+L_{p}, \boldsymbol{y}, k). \end{align}

A key observation is that $\hat {\phi }$ and $\hat {\varphi }$ have a fundamental period of $L_p$ , whereas the functions $\phi$ and $\varphi$ may have infinite fundamental periods. The periodic component of the solution is captured by $\hat {\phi }$ , while the non-periodic part is represented by the plane wave $e^{\mathrm{i}kx}$ . This decomposition reduces the eigenvalue problem from an infinite domain to a finite domain of length $L_p$ .

Two additional properties further simplify the analysis. First, it is straightforward to verify that $\lambda (-k)=\lambda (k)^{*}$ , $\phi (x, \boldsymbol{y}, -k)=\phi (x, \boldsymbol{y}, k)^{*}$ and $\varphi (x, \boldsymbol{y}, -k)=\varphi (x, \boldsymbol{y}, k)^{*}$ are eigenvalues and eigenfunctions. Therefore, it suffices to consider positive real values of $k$ . Second, $\hat {\phi }(x,\boldsymbol{y},k)$ is periodic in $k$ with period $2\pi /L_p$ . To see this, note that

(2.18) \begin{equation} \begin{aligned} \kern2pt \hat{\kern-2pt \mathcal{L}} (k) (e^{\mathrm{i} Gx}\hat {\phi }(x, \boldsymbol{y},k))=e^{\mathrm{i} Gx}\kern2pt \hat{\kern-2pt \mathcal{L}} (k+G)\hat {\phi }(x, \boldsymbol{y},k). \end{aligned} \end{equation}

Consequently, multiplying (2.16) by $e^{-\mathrm{i}x ({2\pi }/{L_{p}})}$ leads to

(2.19) \begin{equation} \begin{aligned} &\lambda (k) e^{-\mathrm{i}x \frac {2\pi }{L_{p}}} \hat {\phi } (x, \boldsymbol{y},k)=e^{-\mathrm{i}x \frac {2\pi }{L_{p}}} \kern2pt \hat{\kern-2pt \mathcal{L}} (k)\hat {\phi } (x, \boldsymbol{y},k)= \kern2pt \hat{\kern-2pt \mathcal{L}}\left(k+\frac{2\pi}{L_p}\right) \left ( e^{\mathrm{i}x \frac {2\pi }{L_{p}}}\hat {\phi } (x, \boldsymbol{y},k) \right )\!. \\ \end{aligned} \end{equation}

Since $e^{\mathrm{i}x ({2\pi }/{L_{p}})}\hat {\phi } (x, \boldsymbol{y},k)$ is periodic in $x$ with period $L_{p}$ , it follows that

(2.20) \begin{equation} \begin{aligned} \lambda \left ( k+\frac {2\pi }{L_{p}} \right )=\lambda (k), \quad \hat {\phi }\left ( x, \boldsymbol{y},k+\frac {2\pi }{L_{p}} \right )= e^{-\mathrm{i}x \frac {2\pi }{L_{p}}}\hat {\phi } (x, \boldsymbol{y},k). \end{aligned} \end{equation}

There are two main reasons for considering the eigenfunctions in the form (2.15). First, the special form of the eigenfunctions defined in (2.16) form a complete biorthogonal system for all square-integrable functions defined on the infinite channel domain. A rigorous justification for the completeness of this biorthogonal system is presented in Appendix A.2. The second reason is that an alternative approach yields the same form of the eigenfunction. It is intuitive to first consider the eigenvalue problem on a channel domain of finite length $L$ and then investigate the behaviour of the eigenfunctions as $L$ approaches infinity. The multiscale asymptotic analysis presented in Appendix A.5 demonstrates that the eigenfunctions take the form (2.15) as $L$ approaches infinity, thereby connecting to the multiscale analysis discussed in Rosencrans (Reference Rosencrans1997).

As a side remark, the eigenfunction representation (2.15) is mathematically equivalent to the Floquet–Bloch form used in solving the Schrödinger equation with periodic potentials (Bloch Reference Bloch1929). This representation was employed by Zwanzig (Reference Zwanzig1983) to study pure diffusion in periodically modulated channels. Allaire, Briane & Vanninathan (Reference Allaire, Briane and Vanninathan2016) compared the two-scale asymptotic expansion and Bloch-wave approaches for the periodic homogenisation of the Poisson equation, showing that the two methods yield identical second-order effective tensors but differ in their fourth-order corrections. To the best of our knowledge, the application of this form to the Taylor dispersion is new.

With the eigenfunctions expressed in the form (2.15), the solution to the advection--diffusion equation can be represented as follows:

(2.21) \begin{equation} \begin{aligned} c(x, \boldsymbol{y}, t) &= \int \limits _{-\frac {\pi }{L_{p}}}^{\frac {\pi }{L_{p}}} \sum \limits _{n=0}^{\infty } \left \langle c_{I}(x, \boldsymbol{y}), \varphi _{n}(x, \boldsymbol{y}, k) \right \rangle \phi _{n}(x, \boldsymbol{y}, k) e^{-\lambda _{n}(k) t} \mathrm{d}k \\&= \frac {1}{2\pi }\sum \limits _{n=0}^{\infty } \int \limits _{-\frac {\pi }{L_{p}}}^{\frac {\pi }{L_{p}}} \big \langle c_{I}(x, \boldsymbol{y}), e^{\mathrm{i} k x} \hat {\varphi }_{n}(x, \boldsymbol{y}, k) \big \rangle \hat {\phi }_{n}(x, \boldsymbol{y}, k) e^{\mathrm{i} k x - \lambda _{n}(k) t} \mathrm{d}k. \end{aligned} \end{equation}

The ordering of the eigenvalues is determined by their real parts at $k=0$ , specifically $0 = \lambda _{0}(0) \lt \textrm{Re} \lambda _{1}(0) \leqslant \ldots \leqslant \textrm{Re} \lambda _{n}(0) \ldots$ . If the real parts are equal when $k=0$ , the ordering is then based on their real parts in the vicinity of zero wavenumber. For non-zero wavenumbers, we order the eigenvalues to ensure that $\lambda _{n}(k)$ remains continuous with respect to $k$ . In this context, it is possible to have $\textrm{Re} \lambda _{n}(k) \lt \textrm{Re} \lambda _{m}(k)$ for some $k$ with $n \gt m$ . Figure 5 shows such an example, which will be discussed in § 3. For the flat-walled channel ( $L_p = 0$ ), (2.21) reduces to the Fourier-transform representation used in Stokes & Barton (Reference Stokes and Barton1990), Phillips & Kaye (Reference Phillips and Kaye1996), Ding & McLaughlin (Reference Ding and McLaughlin2022b ).

2.4. Long-time asymptotic expansion

In this section we use the representation (2.21) to derive the long-time asymptotic expansion of the scalar field. The key observation is that $0 = \min _{k} \textrm{Re} \lambda _{0}(k) \lt \min _{k} \textrm{Re} \lambda _{1}(k) \leqslant \ldots$ implying that all higher-order terms in (2.21) decay exponentially in time. Hence, the leading term associated with $\lambda _{0} (k)$ dominates the long-time behaviour:

(2.22) \begin{align} c(x, \boldsymbol{y}, t) &= \frac {1}{2\pi |\varOmega (L_p)|} \int \limits _{-\frac {\pi }{L_{p}}}^{\frac {\pi }{L_{p}}} \phi _{0}(x, \boldsymbol{y}, k) e^{-\lambda _{0}(k) t} \left (\, \int \limits _{\varOmega } c_{I}(x, \boldsymbol{y}) \varphi _{0}^{*}(x, \boldsymbol{y}, k) \mathrm{d}x \mathrm{d}\boldsymbol{y} \right ) \mathrm{d}k \nonumber\\&\quad + \mathcal{O}\left (e^{-\mathop{\textit{min}} \limits _{k} \textrm{Re} \lambda _{1}(k)t}\right )\!. \end{align}

We denote the integral above as $c_{s}$ and refer to it as the slow manifold of the system for several reasons. First, the integral itself is a solution to the advection--diffusion equation, although with a different initial condition. Second, it captures the algebraically decaying dynamics at long times, while other modes decay exponentially. This integral represents a reduced description of the system’s long-time behaviour, which aligns with the typical characteristics of a slow manifold. The idea behind the approximation (2.22) is utilised in Bronski & McLaughlin (Reference Bronski and McLaughlin1997), Camassa et al. (Reference Camassa, Ding, Kilic and McLaughlin2021), Ding & McLaughlin (Reference Ding and McLaughlin2022b ) to investigate the long-time asymptotic properties in the scalar intermittency problem. Within the centre manifold framework adopted in Mercer & Roberts (Reference Mercer and Roberts1990), Roberts (Reference Roberts2014, Reference Roberts2015), Ding & McLaughlin (Reference Ding and McLaughlin2022a ), the cross-sectional average of the integral in (2.22) represents the centre manifold of the system. Finally, we note that $\lambda _{0}(k)$ could exceed $\lambda _{1}(k)$ for values of $k$ away from zero. If the initial condition is localised around such modes, the approximation (2.22) may no longer be valid. In this work, however, we focus on slowly varying initial conditions, for which the $k=0$ mode is non-zero and the approximation (2.22) provides the dominant contribution to the long-time behaviour.

Further asymptotic analysis can be conducted by noting that the integral in (2.22) is a Laplace-type integral with respect to the wavenumber $k$ . The asymptotic expansion for large $t$ can be derived using the steepest descent method, as outlined in Bender, Orszag & Orszag (Reference Bender, Orszag and Orszag1999), Ablowitz, Fokas & Fokas (Reference Ablowitz, Fokas and Fokas2003). Intuitively, $\textrm{Re} \lambda _{0} = 0$ only when the wavenumber is zero, while it is positive for non-zero wavenumbers. For sufficiently large $t$ , $e^{-\lambda _{0}(k) t}$ becomes exponentially small for non-zero wavenumbers. Therefore, the integrand near the zero wavenumber contributes dominantly when $t$ is large. To compute the long-time asymptotic expansion, we expand the eigenvalue and eigenfunction around $k=0$ using the perturbation procedure detailed in Appendix A.3:

(2.23) \begin{equation} \begin{aligned} \lambda _{0} &= \textrm{i}{\textit{Pe}} \left \langle u \right \rangle k + \left \langle 1 - {\textit{Pe}} \left ( u - \left \langle u \right \rangle \right ) \theta _{0}^{(1,1)} + \partial _{x} \theta _{0}^{(1,1)} \right \rangle k^{2} \\ & \quad - \mathrm{i} \left \langle \partial _{x} \theta _{0}^{(2,2)} - {\textit{Pe}} \left ( u - \left \langle u \right \rangle \right ) \theta _{0}^{(2,2)} \right \rangle k^{3} + \mathcal{O}(k^{4}), \\ \phi _{0} &= e^{\mathrm{i} k x} \left ( 1 + \mathrm{i} k \theta _{0}^{(1,1)} - k^{2} \left ( \theta _{0}^{(2,2)} + \frac {1}{2}\left \langle \theta _{0}^{(1,1)}, q_{0}^{(1,1)} \right \rangle \right ) + \mathcal{O}(k^{3}) \right )\!, \\ \varphi _{0} &= e^{\mathrm{i} k x} \left ( 1 + \mathrm{i} k q_{0}^{(1,1)} - k^{2} \left ( q_{0}^{(2,2)} + \frac {1}{2}\left \langle \theta _{0}^{(1,1)}, q_{0}^{(1,1)} \right \rangle \right ) + \mathcal{O}(k^{3}) \right ). \end{aligned} \end{equation}

The functions $\theta _{0}^{(1,1)}, \theta _{0}^{(2,2)}, q_{0}^{(1,1)}$ are periodic in $x$ with a fundamental period $L_{p}$ and have zero mean over the unit cell domain, satisfying $\langle \theta _{0}^{(1,1)}, 1 \rangle = \langle \theta _{0}^{(2,2)}, 1 \rangle = \langle 1, q_{0}^{(1,1)} \rangle = 0$ . The function $\theta _{0}^{(1,1)}$ satisfies

(2.24) \begin{equation} \begin{aligned} \kern2pt \hat{\kern-2pt \mathcal{L}}_{0}(0) \theta _{0}^{(1,1)} &= - {\textit{Pe}} \left ( u - \left \langle u \right \rangle \right )\!, \quad \left . n_{x} + n_{x} \partial _{x} \theta _{0}^{(1,1)} + \boldsymbol{n}_{\boldsymbol{y}} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \theta _{0}^{(1,1)} \right |_{w\textit{all}} = 0. \end{aligned} \end{equation}

The function $q_{0}^{(1,1)}$ satisfies

(2.25) \begin{equation} \begin{aligned} &\kern-1pt \hat{\kern1pt \mathcal{L}^{*}_{0}}(0) q_{0}^{(1,1)} = {\textit{Pe}} \left ( u - \left \langle u \right \rangle \right )\!, \; \left . n_{x} + n_{x} \partial _{x} q_{0}^{(1,1)} + \boldsymbol{n}_{\boldsymbol{y}} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} q_{0}^{(1,1)} \right |_{w\textit{all}} = 0. \end{aligned} \end{equation}

The function $\theta _{0}^{(2,2)}$ satisfies

(2.26) \begin{equation} \begin{aligned} &\kern2pt \hat{\kern-2pt \mathcal{L}}_{0}(0) \theta _{0}^{(2,2)} = {\textit{Pe}} \Big \langle u \theta _{0}^{(1,1)} \Big \rangle - {\textit{Pe}} \left ( u - \left \langle u \right \rangle \right ) \theta _{0}^{(1,1)} - \Big \langle \partial _{x} \theta _{0}^{(1,1)} \Big \rangle + 2 \partial _{x} \theta _{0}^{(1,1)}, \\&\left . n_{x} \theta _{0}^{(1,1)} + n_{x} \partial _{x} \theta _{0}^{(2,2)} + \boldsymbol{n}_{\boldsymbol{y}} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \theta _{0}^{(2,2)} \right |_{w\textit{all}} = 0. \end{aligned} \end{equation}

In addition to the eigenfunction expansion, we need the expansion of the integral related to the initial condition for small wavenumbers:

(2.27) \begin{align} &\frac {1}{|\varOmega (L_p)|}\int \limits _{\varOmega } c_{I}(x, \boldsymbol{y}) \varphi _{0}^{*}(x, \boldsymbol{y}, k) \mathrm{d}x \mathrm{d}\boldsymbol{y} \nonumber \\& \ = \frac {1}{|\varOmega (L_p)|} \!\int \limits _{\varOmega }\! c_{I}(x, \boldsymbol{y}) e^{-\mathrm{i} k x} \!\left (\! 1 - \mathrm{i} k q_{0}^{(1,1)} - k^{2} \!\left (\! q_{0}^{(2,2)} + \frac {1}{2}\left \langle \theta _{0}^{(1,1)}, q_{0}^{(1,1)} \right \rangle \!\right )\! + \mathcal{O}(k^{3})\! \right )\! \mathrm{d}x \mathrm{d}\boldsymbol{y} \nonumber \\& \ = \frac {1}{|\varOmega (L_p)|}\int \limits _{\varOmega } c_{I}(x, \boldsymbol{y}) \left ( 1 - \mathrm{i} k x + \frac {(-\mathrm{i} k x)^{2}}{2} + \mathcal{O}(k^{3}) \right ) \nonumber \\& \ \quad \times \left ( 1 - \mathrm{i} k q_{0}^{(1,1)} - k^{2} \left ( q_{0}^{(2,2)} + \frac {1}{2}\left \langle \theta _{0}^{(1,1)}, q_{0}^{(1,1)} \right \rangle \right ) + \mathcal{O}(k^{3}) \right ) \mathrm{d}x \mathrm{d}\boldsymbol{y} \nonumber \\& \ = c_{I}^{(0)} - \mathrm{i} k c_{I}^{(1)} - k^{2} c_{I}^{(2)} + \mathcal{O}(k^{3}). \end{align}

Here

(2.28) \begin{align} c_{I}^{(0)}&=\frac {1}{|\varOmega (L_p)|}\int \limits _{\varOmega }c_{I}(x,\boldsymbol{y}) \mathrm{d}x \mathrm{d}\boldsymbol{y}, \quad c_{I}^{(1)}=\frac {1}{|\varOmega (L_p)|} \int \limits _{\varOmega }c_{I}(x,\boldsymbol{y}) \left ( x+q_{0}^{(1,1)} \right )\mathrm{d}x \mathrm{d}\boldsymbol{y},\nonumber \\ c_{I}^{(2)}&=\frac {1}{|\varOmega (L_p)|} \int \limits _{\varOmega }c_{I}(x,\boldsymbol{y}) \left ( \frac {x^{2}}{2}+ q_{0}^{(2,2)} + \frac {1}{2}\left \langle \theta _{0}^{(1,1)}, q_{0}^{(1,1)} \right \rangle +q_{0}^{(1,1)} x \right )\mathrm{d}x \mathrm{d}\boldsymbol{y}. \end{align}

The steepest descent method requires expanding $\lambda _{0}(k)$ around the saddle point in the complex plane. Note that the expansion (2.23) implies that the saddle point of $\lambda _{0}(k)$ is not located at $k=0$ if $\left \langle u \right \rangle \neq 0$ . In fact, this expansion around this saddle point provides us with an asymptotic expansion of (2.22) as $t \rightarrow \infty$ while keeping other variables fixed. Notably, this type of expansion is not uniform with respect to the spatial variable. To achieve a uniform asymptotic expansion, we should notice that the whole scalar is transported by the flow with the non-zero mean flow speed. Therefore, we analyse the problem in the moving frame of reference $x_{m} = x - {\textit{Pe}} \left \langle u \right \rangle t$ . In this frame, $k=0$ is the saddle point of $\lambda _{0}(k)$ and the steepest decent method yields a uniform asymptotic expansion.

Substituting the expansions (2.23) and (2.27) into (2.22) yields the asymptotic expansion of the scalar field at long times:

(2.29) \begin{equation} \begin{aligned} &c (x,\boldsymbol{y},t) =\frac { e^{-\frac {x_{m}^2}{2 t \lambda ''(0)}}}{|\varOmega (L_p)|\sqrt {2 \pi t \lambda ''(0)}} \int \limits _{\varOmega }c_{I}(x,\boldsymbol{y}) \mathrm{d}x \mathrm{d}\boldsymbol{y} +\frac { x_{m} e^{-\frac {x_{m}^2}{2 t \lambda ''(0)}}}{\sqrt {2 \pi }\left (t \lambda ''(0)\right )^{3/2}} \!\left ( c_{I}^{(1)}-\theta _{0}^{(1,1)} c_{I}^{(0)}\! \right )\\ &+\frac { e^{-\frac {x_{m}^2}{2 t \lambda ''(0)}} \!\left (x_{m}^2-t \lambda ''(0)\right ) \!\left ( \!\left ( \!\theta _{0}^{(2,2)} \!+\! \dfrac {1}{2}\left \langle \theta _{0}^{(1,1)}, q_{0}^{(1,1)} \right \rangle \right ) c_{I}^{(0)}- \theta _{0}^{(1,1)} c_{I}^{(1)}+c_{I}^{(2)} \!\right )}{\sqrt {2 \pi }\left (t \lambda ''(0)\!\right )^{5/2}}+\mathcal{O} \!\left ( t^{-\frac {3}{2}} \!\right )\!. \\ \end{aligned} \end{equation}

Note that not all terms of order $\mathcal{O}(t^{-3/2})$ are included above. Capturing the full expansion to this order requires extending the eigenvalue expansion to $\mathcal{O}(k^{5})$ and the eigenfunction expansion to $\mathcal{O}(k^{3})$ , resulting in a lengthy expression that is omitted here for brevity.

In Fourier space, multiplication by $k$ corresponds to differentiation in physical space. Since derivatives of the Gaussian function yield Hermite polynomials multiplied by the Gaussian, each higher-order term in (2.29) takes the form of a Hermite polynomial times a Gaussian, with coefficients determined by the eigenfunction expansion.

A more accurate approximation can be obtained by incorporating information about the initial condition. For instance, if the initial condition is a Gaussian distribution $\mathcal{N}(x; x_0,\sigma )$ centred at $x=x_{0}$ with standard deviation $\sigma$ ,

(2.30) \begin{equation} \begin{aligned} & c_{I} \left ( x,\boldsymbol{y}\right )= \mathcal{N}(x; x_0,\sigma )\equiv \frac {1}{\sqrt {2 \pi \sigma ^2}} \exp \left (-\frac {(x - x_0)^2}{2\sigma ^2}\right )\!.\end{aligned} \end{equation}

In this case, instead of approximating $e^{-\mathrm{i} k x}$ using its small $k$ expansion as in (2.27), we can directly evaluate the integral involving the initial condition. This approach yields an expansion that is asymptotically equivalent to (2.29) at long times while significantly improving accuracy at short times.

The leading-order term in the approximation (2.29) is a Gaussian function, which corresponds to the solution of a one-dimensional advection–diffusion equation with an enhanced diffusion coefficient:

(2.31) \begin{equation} \partial _{t}c + {\textit{Pe}} \left \langle u \right \rangle \partial _{x}c = \kappa _{\textit{eff}} \partial _{x}^{2}c, \quad \kappa _{\textit{eff}} = \frac {\lambda ''(0)}{2} = \left \langle 1 - {\textit{Pe}} \left ( u - \left \langle u \right \rangle \right ) \theta _{0}^{(1,1)} + \partial _{x} \theta _{0}^{(1,1)} \right \rangle\!. \end{equation}

This equation serves as the effective long-time approximation of (2.7). The coefficient $\kappa _{\textit{eff}}$ represents the effective diffusivity, which is consistent with the definition given in (2.11). To show that the effective diffusivity is real and non-negative, we rewrite it as

(2.32) \begin{equation} \kappa _{\textit{eff}} = \left \langle \boldsymbol{\nabla} _{\boldsymbol{y}} \theta _{0}^{(1,1)} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{y}} \theta _{0}^{(1,1)} + (1 + \partial _{x} \theta _{0}^{(1,1)})^{2} \right \rangle \geqslant 0, \end{equation}

where the inequality follows from the fact that $\theta _{0}^{(1,1)}$ is real. This expression for the effective diffusivity aligns with previously reported results in Brenner (Reference Brenner1980), Rosencrans (Reference Rosencrans1997). (The function $\theta _0^{(1,1)}$ satisfies the same cell problem as $\chi _i$ in Rosencrans (Reference Rosencrans1997, p. 1225), and thus, (2.32) corresponds directly to the expression for the enhanced diffusion coefficient given on p. 1226 of the same reference.) Furthermore, the validity of this formula has been confirmed through benchmarking against particle tracking methods, as demonstrated in Richmond et al. (Reference Richmond, Perkins, Scheibe, Lambert and Wood2013), Haugerud et al. (Reference Haugerud, Linga and Flekkøy2022).

Finally, there is an intriguing question that is difficult to address in other frameworks but becomes more tractable in the proposed approach: How does the effective diffusivity change when the flow direction is reversed, i.e. when ${\textit{Pe}}$ is replaced by $-{\textit{Pe}}$ ? This question arose from examining the formula for the effective diffusivity (2.31). At first glance, the forcing term suggests that $\partial _{x} \theta _{0}^{(1,1)}$ scales as ${\textit{Pe}}$ , while ${\textit{Pe}} (u - \langle u \rangle )\theta _{0}^{(1,1)}$ scales as ${\textit{Pe}}^{2}$ . This scaling would imply the possible presence of a ${\textit{Pe}}$ -order term in the effective diffusivity, potentially making it a non-even function of ${\textit{Pe}}$ .

Replacing ${\textit{Pe}}$ by $-{\textit{Pe}}$ essentially switches the roles of the eigenvalue problem (2.12) and its adjoint counterpart (2.13). Consequently, the effective diffusivity changes from $\kappa _{\textit{eff}} = {\lambda ''(0)}/{2}$ to $\kappa _{\textit{eff}} = {(\lambda ''(0))^{*}}/{2}$ . Since $\lambda ''(0)$ is real, the effective diffusivity remains unchanged. One alternative approach to demonstrate this symmetry is to compute a power series expansion of the effective diffusivity in terms of ${\textit{Pe}}$ and show that all odd-order terms vanish, implying that the effective diffusivity is an even function of ${\textit{Pe}}$ . However, this calculation is relatively lengthy, whereas the symmetry arises more directly within the framework proposed in this work.

2.5. Connection to the time scale for dispersion, longitudinal normality and transverse uniformity

From the derivation in the previous section, we observe that the effective equation (2.31) is obtained through two successive approximations. The first approximation assumes that the scalar field can be described by its slow manifold with exponentially decaying correction terms. This approximation is valid on the time scale estimated by

(2.33) \begin{equation} t_s = \frac {1}{\min _k \textrm{Re} \lambda _1(k)}, \end{equation}

which depends solely on the domain geometry and velocity field. A larger $\min_k \textrm{Re} \lambda _1(k)$ leads to faster convergence, whereas a smaller $\min _k \textrm{Re} \lambda _1(k)$ results in slower convergence. The second approximation expresses the slow manifold as a Gaussian distribution plus algebraically decaying correction terms. Therefore, validity of the Gaussian distribution approximation depends on the decay of these correction terms, as given in (2.29), and the initial condition.

To better understand the relationship between the time scales for dispersion, longitudinal normality and transverse uniformity, we consider an explicit example: a flat channel domain

(2.34) \begin{equation} \varOmega = \left\{ (x, y) \mid x \in \mathbb{R}, -\frac {1}{2} \leqslant y \leqslant \frac {1}{2} \right\} \end{equation}

subject to shear flow $ \boldsymbol{u} = (({3}/{2})(1 - 4y^2),0)$ and an initial condition given by (2.30). To characterise the time scale for transverse uniformity, we analyse the long-time asymptotic expansion of the scalar field. For longitudinal normality, we evaluate the asymptotic behaviour of the cross-sectionally averaged scalar. To determine the dispersion time scale, we examine the long-time expansion of the scalar variance.

Using the method outlined in the previous section, we obtain the long-time asymptotic expansion of the scalar field:

(2.35) \begin{equation} \begin{aligned} c=&\frac {e^{-\frac {x^2}{2 v}}}{\sqrt {2 \pi v }}\left ( 1+\frac {a_{1}\left (240 y^4-120 y^2+7\right )}{480 } +\frac {a_{2} \left (240 y^2 \left (720 y^6-784 y^4+238 y^2-17\right )-413\right )}{3225600}\right .\\ & -\frac {a_{3}\left (1064960 t-520 \left (-15686 y^2+32 \left (10800 y^6-18568 y^4+10725 y^2-2079\right ) y^4+8447\right ) y^2+240261\right )}{73801728000 } \\ & \frac {a_{4}}{16862218813440000} \left ( 763467599-38993920 t \left (1560 \left (2 y^2-1\right ) y^2+283\right )+ 2720 y^2 \right . \\ &\hspace {2.3cm}\times \left (1127049 -26712502 y^2 +101146864 y^4-78215280 y^6 \right . \\ &\hspace {2.5cm} \left . \left . \left . -259440896 y^8 +779717120 y^{10} -897085440 y^{12} + 377395200 y^{14}\right ) \right ) \right )+\mathcal{O} \left(t^{-\frac {3}{2}}\right)\!. \end{aligned} \end{equation}

Here $ v = \sigma ^2 + 2 \kappa _{\textit{eff}} t$ and $ x_m = x - x_0 - {\textit{Pe}} t$ . The coefficients $ a_n$ are given by

(2.36) \begin{equation} a_n = \frac {{\textit{Pe}}^n}{(2 v)^{n/2}} H_n\left (\frac {x_m}{\sqrt {2 v}}\right )\!, \end{equation}

where $ H_n$ is the Hermite polynomial of degree $ n$ , $H_{1} (x)=2x, H_{2} (x)=4 x^2-2, H_{3} (x)=8 x^3-12 x, H_{4} (x)=16 x^4-48 x^2+12$ . Taking the cross-sectional average of (2.35) yields

(2.37) \begin{equation} \begin{aligned} \bar {c}=\int \limits _{-\frac {1}{2}}^{\frac {1}{2}} c\mathrm{d} y=\frac {e^{-\frac {x^2}{2 v}}}{\sqrt {2 \pi v }} \left ( 1-\frac {a_{2}}{8400}-\frac {a_{3}t}{69300 v^3} -\frac {(685440 t-54991)a_{4}}{1543782240000}\right )+\mathcal{O} \left(t^{-\frac {3}{2}}\right)\!. \end{aligned} \end{equation}

Next, we derive the long-time asymptotic expansion of the variance, which requires computing $ \langle x_m^2 c \rangle$ . The higher-order terms omitted from (2.35) are Hermite polynomials of degree greater than four. Due to the orthogonality of Hermite polynomials, these terms do not contribute to the variance. Thus, we obtain the exact variance of the slow manifold $ c_s$ and the asymptotic expansion of the scalar field variance with exponential correction terms, which is consistent with (3.22) in Camassa et al. (Reference Camassa, Lin and McLaughlin2010a ) (due to a different choice of characteristic velocity, comparison with the formula in Camassa et al. Reference Camassa, Lin and McLaughlin2010a requires replacing ${\textit{Pe}}$ in the present formula with $({2}/{3}){\textit{Pe}}$ ):

(2.38) \begin{equation} \begin{aligned} \text{Var}_{c_s}(t) &= \left \langle x_m^2 c_s \right \rangle = -\frac {{\textit{Pe}}^2}{4200} + \sigma ^2 + 2 \left ( \frac {{\textit{Pe}}^2}{210} + 1 \right )t, \\ \text{Var}_{c}(t) &= \left \langle x_m^2 c \right \rangle = -\frac {{\textit{Pe}}^2}{4200} + \sigma ^2 + 2 \left ( \frac {{\textit{Pe}}^2}{210} + 1 \right )t + \mathcal{O}\left (e^{-\mathop{\textit{min}} _k \textrm{Re} \lambda _1(k)t}\right )\!.\end{aligned} \end{equation}

Here $ \min _k \textrm{Re} (\lambda _1(k) )= \pi ^2$ . The variance of the slow manifold at $ t=0$ is given by $ \text{Var}_{c_{s}}(0)=-( {{\textit{Pe}}^2}/{4200})+ \sigma ^{2}$ , which differs from the variance of the initial condition of the scalar field, $ \text{Var}_{c}(0)=\sigma ^{2}$ . This suggests that the slow manifold corresponds to a solution with a different initial condition.

Now, we can analyse the three different time scales mentioned at the beginning. From (2.38), we observe that once the exponential terms have decayed, the variance increases linearly in time. Therefore, the characteristic time scale for dispersion is given by $t_{s} = 1/\min _{k} \textrm{Re} \lambda _{1}(k)$ .

Equation (2.37) becomes valid after this time scale is reached. The first term represents a Gaussian function, while the remaining terms describe deviations from Gaussianity. Since these correction terms have small coefficients, they become negligible unless $ {\textit{Pe}}$ is large or the initial variance is small. In such cases, after time $ t_s$ , the cross-sectionally averaged scalar field is well approximated by a Gaussian function.

Similarly, (2.35) is valid after time $ t_s$ . The first term represents the Gaussian function, while the second term describes transverse variations. Since the cross-sectional average of this term is zero, it does not appear in (2.37). As time increases, when the second term becomes negligible compared with the first, the transverse uniformity time scale is reached.

2.6. Validation by the numerical simulation

After analysing the example in a flat channel, we now use numerical methods to demonstrate that the proposed theory also holds in a curved channel. This is done by simulating the scalar transport (2.2) in the following sinusoidal channel:

(2.39) \begin{equation} \begin{aligned} \varOmega =\left \{ (x, y) | x \in \mathbb{R}, 0\leqslant y\leqslant h(x) \right \},\quad h(x) = 1 + A \cos \left (\frac {2\pi x}{L_{p}}\right )\!. \end{aligned} \end{equation}

Here $A\lt 1$ and $L_p$ is the fundamental period.

In the simulation we consider the sinusoidal channel (2.39) with $A = 0.2, L_{p}=1$ . The computational domain has a length of 42 with periodic boundary conditions at the left and right ends. To confirm that this domain provides an adequate approximation of an unbounded channel, we verify that the scalar concentration at both endpoints remains close to zero throughout the simulation. A steady pressure-driven flow with ${\textit{Re}} = 200$ and a pressure gradient of $f_x = -20.18$ is computed using the method described in Appendix B. The spatial mean of the resulting flow $\left \langle u \right \rangle =1$ . After obtaining the flow, we use it to advect the scalar by solving the advection--diffusion equation (2.2) with ${\textit{Pe}} = 6$ . The initial condition follows the Gaussian distribution given in (2.30), with $\sigma = 0.5$ and $x_{0} = 6$ . We find that $\max _k \textrm{Re} \lambda _1(k) = \lambda _1(0) = 9.0294$ , and the estimated time scale for the scalar field to converge to the slow manifold is $t_s = 1/\min _k \textrm{Re} \lambda _1(k) = 0.1107$ . The effective diffusivity is $\kappa _{\textit{eff}} = 1.6721$ .

We plot $({1}/{2}) \partial _{t} \text{Var}(t)$ as a function of time in figure 2, showing its convergence to the theoretical effective diffusivity. Once $({1}/{2}) \partial _{t} \text{Var}(t)$ stabilises to a constant value, it indicates that the variance grows linearly in time, marking the onset of Taylor dispersion. The relative differences between $({1}/{2}) \partial _{t} \text{Var}(t)$ and $\kappa _{\textit{eff}}$ are 0.1069, 0.0368 and 0.0132 at $t = t_s$ , $2t_s$ and $3t_s$ , respectively. The ratio of the adjacent relative differences is $2.7 \sim 2.8$ , indicating that $t_s$ is the e-folding time scale (the e-folding time is the time interval in which an exponentially growing or decaying quantity increases or decreases by a factor of $e$ ) for the system, thereby supporting the proposed theory.

Figure 2. The quantity $({1}/{2}) \partial _{t} \text{Var}(t)$ as a function of time is plotted as the red solid curve. The theoretical long-time limit predicted by (2.31) is indicated by the blue dashed curve. The vertical dotted lines indicate $t_{s}$ and $2t_{s}$ , where $t_{s}$ represents the estimated time scale for the scalar field to converge to the slow manifold.

Next, we investigate longitudinal normality. We measure the deviation of the cross-sectional average of the scalar field $\bar {c} (x) = ({1}/{h (x)})\int \nolimits _{0}^{h (x)}c (x,y)\mathrm{d}y$ to the Gaussian function measured from a Gaussian profile by the metric

(2.40) \begin{equation} E_{n,1} (t) =\frac {\lVert \bar {c} (x,t) - f (x,t)\rVert _{\infty } }{ \lVert \bar {c}\rVert _{\infty } }, \quad f(x,t)=\mathcal{N}\big ( x; \left \langle x c \right \rangle , \sqrt { \mathrm{Var} (t) } \big ), \end{equation}

where $\lVert \boldsymbol{\cdot }\rVert _{\infty }$ denotes the $L^{\infty }$ norm and $\mathcal{N}$ is the Gaussian function defined in (2.30). Here $\langle x c \rangle$ and $\mathrm{Var}(t)$ are the mean and variance of the scalar field, respectively. A smaller $E_{n,1}(t)$ indicates closer agreement between $\bar {c}$ and the corresponding Gaussian.

To quantify the approach to transverse uniformity, we define

(2.41) \begin{equation} E_{t,1} (t) =\frac {\lVert c (x,y,t) - f(x,t) \rVert _{\infty } }{ \lVert c\rVert _{\infty } }, \end{equation}

where the same Gaussian $f(x,t)$ is used as a reference. This measures the maximum pointwise deviation from a longitudinally varying but transversely uniform Gaussian state.

The Gaussian profile $f(x,t)$ corresponds to the first term in the asymptotic expansion of the scalar field (2.29). It is also useful to examine the approximation error when the first two terms in (2.29) are included:

(2.42) \begin{equation} \begin{aligned} &E_{n,2} (t) =\frac {\lVert \bar {c}- \left ( f+ \partial _{x}f \left ( c_{I}^{(1)}-\overline {\theta _{0}^{(1,1)}} c_{I}^{(0)} \right ) \right ) \rVert _{\infty } }{ \lVert \bar {c}\rVert _{\infty } }, \\ &E_{t,2} (t) =\frac {\lVert c- \left ( f+ \partial _{x}f \left ( c_{I}^{(1)}-\theta _{0}^{(1,1)} c_{I}^{(0)} \right ) \right ) \rVert _{\infty } }{ \lVert c\rVert _{\infty } }. \\ \end{aligned} \end{equation}

If $E_{n,1} \gg E_{n,2}$ and $E_{t,1} \gg E_{t,2}$ , the magnitude of the second term provides a practical estimate for the approximation error of the first term, and hence, for the time scales of longitudinal normality and transverse uniformity. In this case, $c_{I}^{(0)} = 1$ . Both $\theta _{0}^{(1,1)}$ and $c_{I}^{(1)}$ are computed numerically. The numerically obtained value of $c_{I}^{(1)}$ is small, of the order of $10^{-5}$ .

Figure 3. (a) Red solid curve: $\bar {c}$ ; blue dashed curve: Gaussian profile with the same mean and variance as $\bar {c}$ ; black dotted curve: first-order correction term. (b) Red solid curve: $E_{n,1}(t)$ ; blue dashed curve: $E_{n,2}(t)$ .

Figure 4. (a) Top: solution of the advection–diffusion equation (2.7) at $t = 0.2$ . Bottom: two-term approximation of the scalar field from (2.29). (b) Red solid curve: $E_{t,1}(t)$ ; blue dashed curve: $E_{t,2}(t)$ .

For a channel with flat walls, deviation of $\bar {c}$ from the Gaussian profile typically appears as non-zero skewness. In periodically modulated channels, in addition to skewness, $\bar {c}$ may exhibit periodic fluctuations due to the domain geometry. As shown in figure 3(a), $\bar {c}$ is nearly symmetric but still contains periodic fluctuations, which are captured by the second term in (2.29). This is confirmed by figure 3(b), which shows that $E_{n,1} \gt E_{n,2}$ .

Figure 4(a) demonstrates that including the second term in (2.29) significantly improves the capture of transverse variations of the entire scalar field. Figure 4(b) demonstrates that $E_{t,1} \gt E_{t,2}$ , supporting the use of the second term’s magnitude to estimate the time scales of both longitudinal normality and transverse uniformity.

3.1. Shear flow in the channel with flat boundaries

We first examine our theory using a fundamental case: the diffusion of a passive scalar transported by a shear flow given by $u = u(y), v = 0$ in a channel with flat boundaries, represented as $\varOmega = \{ (x, y) | x \in \mathbb{R}, 0 \leqslant y \leqslant 1 \}$ .

In this scenario, the fundamental period of the system is zero, and the unit cell reduces to $\varOmega (0) = \{ (0, y) | 0 \leqslant y \leqslant 1 \}$ . Therefore, $\hat {\phi }$ is independent of $x$ , and (2.16) simplifies to

(3.1) \begin{equation} \begin{aligned} &\lambda \hat {\phi } = \kern2pt \hat{\kern-2pt \mathcal{L}}(k) \hat {\phi } = {Pei} k u \hat {\phi } - \big(\!-k^{2} \hat {\phi } + \Delta _{\boldsymbol{y}} \hat {\phi } \big ). \end{aligned} \end{equation}

Figure 5. The first four eigenvalues for different shear flows with ${\textit{Pe}}=1$ .

Next, we numerically compute the eigenvalues for several different shear flows. An interesting family of shear flow is $u = \sqrt {2} \cos (n \pi y)$ , where $n \in \mathbb{N}^{+}$ . These functions serve as the eigenfunctions of $\kern2pt \hat{\kern-2pt \mathcal{L}}(0)$ . Many analyses regarding Taylor dispersion focus on slowly varying solutions, where the energy is concentrated in the vicinity of zero wavenumber. Therefore, the eigenvalues of $\kern2pt \hat{\kern-2pt \mathcal{L}}(0)$ can provide an estimate of the decay rate (Mercer & Roberts Reference Mercer and Roberts1990; Rosencrans Reference Rosencrans1997). However, it is important to note that $\textrm{Re} \lambda _{1}$ does not always reach its minimum value at $k = 0$ , which is a crucial consideration when the spectrum of the solution is not concentrated near zero wavenumber. For instance, $u = \sqrt {2} \cos (\pi y)$ exemplifies this behaviour. Figure 5(a,b) displays the real and imaginary parts of the first four eigenvalues for this shear flow with ${\textit{Pe}} = 1$ . The $k^{2}$ term arises from the contribution of diffusion in the longitudinal direction of the channel, and we subtract this term from the eigenvalue to highlight the contribution from the flow. Here $\textrm{Re} (\lambda _{1}) - k^{2}$ decreases for small wavenumbers, with its global minimum occurring at $k = 5.1252$ , yielding a value of 5.1537. The calculations in Appendix A.4 show the asymptotic expansion of $\lambda _{1}$ for small wavenumbers:

(3.2) \begin{equation} \lambda _{1} = \pi ^{2} + \left ( 1 - \frac {5{\textit{Pe}}^{2}}{6\pi ^{2}} \right ) k^{2} + \mathcal{O}(k^{3}). \end{equation}

This implies that if ${\textit{Pe}} \gt \sqrt {({6}/{5})} \pi \approx 3.4414$ then $k = 0$ may not be a local minimum of $\lambda _{1}(k)$ . In fact, when ${\textit{Pe}} \gt 2.36$ , the global minimum is not at $k = 0$ . Furthermore, as ${\textit{Pe}}$ increases, the location of the minimum approaches $k = 0$ , as shown in figure 6(a).

Figure 6. Plot of $\textrm{Re} (\lambda _{1} (k))$ as a function of the wavenumber, $k$ , for different ${\textit{Pe}}$ and different flows. The $y$ axis is in log scale. Results are shown for (a) $u=\sqrt {2}\cos \pi y$ and (b) $u=\sqrt {2}\cos 2\pi y$ .

Interestingly, when the shear flow is given by $u = \sqrt {2} \cos (n \pi y)$ with $n \geqslant 2$ , $k = 0$ is a local minimum of $\lambda _{1}(k)$ , as supported by the asymptotic analysis for small wavenumbers presented in Appendix A.4. Figure 5(c,d) displays the real and imaginary parts of the first four eigenvalues for the shear flow $u = \sqrt {2} \cos (2 \pi y)$ , further supporting this observation. Additionally, in conjunction with figure 6(b), we observe that $\lambda _{1}(0) = \pi ^{2}$ serves as the global minimum across different ${\textit{Pe}}$ values.

Therefore, for $u = \sqrt {2} \cos (n \pi y)$ with $n \geqslant 2$ , we have $ \min _k \textrm{Re} (\lambda _1(k) )= \pi ^2$ and $ t_s = 1/\min _k \textrm{Re} \lambda _1(k) = 1/\pi ^{2}$ , whereas for $u = \sqrt {2} \cos (\pi y)$ , it follows that $t_s \geqslant 1/\pi ^{2}$ . These results indicate that if the flow strength is sufficiently large, the scalar diffusing under the shear flow $u = \sqrt {2} \cos (\pi y)$ could converge more slowly to the slow manifold than when advected by $u = \sqrt {2} \cos (n \pi y)$ with $n \geqslant 2$ .

Next, we consider two practical flows: Couette flow given by $u = 2\sqrt {2} (y - ({1}/{2}) )$ and plane Poiseuille flow given by $u = 6\sqrt {5} (y (1-y) - ({1}/{6}) )$ . The coefficients in these expressions are chosen to ensure that $\left \langle u, u \right \rangle = 1$ , allowing for a fair comparison with the cosine shear flows. Couette flow describes the motion of a viscous fluid in the space between two surfaces, one of which moves tangentially relative to the other. Its velocity profile is similar to $u = \sqrt {2} \cos (\pi y)$ . In contrast, plane Poiseuille flow represents the steady, laminar flow of a viscous fluid between two stationary, parallel plates, driven by a pressure gradient, structurally resembling $u = \sqrt {2} \cos (2\pi y)$ .

Consequently, we can anticipate that the eigenvalues associated with these flows will be similar to those of the cosine flows. The results for ${\textit{Pe}} = 1$ are presented in figure 7. For Couette flow, $\textrm{Re} \lambda _{1} - k^{2}$ reaches a global minimum value of 5.2149 at $k = 5.2149$ . In contrast, for plane Poiseuille flow, $\textrm{Re} \lambda _{1} - k^{2}$ attains its global minimum value of $\pi ^{2}$ at $k = 0$ . These results imply that the diffusing scalar advected by the plane Poiseuille flow converges more slowly to the slow manifold compared with the scalar advected by the Couette flow.

Figure 7. The first four eigenvalues for different shear flows with ${\textit{Pe}}=1$ . Results are shown for (a) $u= 2\sqrt {2} (y- ({1}/{2}))$ , (b) $u=2\sqrt {2} (y- ({1}/{2}))$ , (c) $u=6\sqrt {5} (y (1-y)-({1}/{6}))$ , (d) $u=6\sqrt {5} (y (1-y)-({1}/{6}))$ .

Why may the global minimum of $\textrm{Re} (\lambda _1(k))$ occur at a non-zero wavenumber? While a complete explanation remains open, we offer two conjectures. From a mathematical point of view, the asymptotic expansion of the eigenvalue for a small wavenumber reveals that interactions among different eigenfunctions contribute to the derivative of $\lambda _1(k)$ and thereby determine whether $\textrm{Re} (\lambda _1(k))$ attains a local minimum at $k = 0$ . However, these interactions are generally complex and difficult to characterise analytically. For the flat parallel plate channel, the eigenfunctions at zero wavenumber take the form $\phi _{n}= \sqrt {2} \cos (n \pi y)$ . When the velocity field is expanded in this basis, calculations show that only the component corresponding to $\phi _1 = \sqrt {2} \cos (\pi y)$ causes $\textrm{Re} (\lambda _1(k))$ to fail to attain a local minimum at $k = 0$ for certain values of ${\textit{Pe}}$ . Since $\phi _1$ is also the eigenfunction associated with $\lambda _1$ , we conjecture that, in more general settings, the projection of the velocity field onto $\phi _1$ may similarly affect whether the global minimum of $\textrm{Re} (\lambda _1(k))$ occurs at a non-zero wavenumber.

From a physical perspective, flows such as $u = \sqrt {2} \cos (\pi y)$ or $u = ({1}/{2}) - y$ shear the scalar field in opposite directions across the channel: fluid near $y = 0$ is advected downstream, while fluid near $y = 1$ is advected upstream. This antisymmetric shearing creates a configuration where scalar filaments are stretched in opposing directions, potentially suppressing transverse dispersion and slowing convergence to the asymptotic state. In contrast, flows like $u = \sqrt {2} \cos (n\pi y)$ for $n \geqslant 2$ contain more than two extrema across the channel height. These more oscillatory profiles lead to multiple shearing layers, which may enhance stretching and increase transverse concentration gradients, thereby accelerating the decay to the slow manifold.

These interpretations remain conjectural, but they suggest interesting connections between the spectral structure of the advection--diffusion operator and the geometry of the underlying flow. A more detailed investigation of this phenomenon is left for future work.

3.2. Cellular flow in the channel with flat boundaries

Shear flow is a unidirectional flow and has a transverse velocity component of zero. To examine the effect of the transverse velocity component in a channel, we consider scalar transport in a channel domain with flat boundaries but subject to a cellular flow given by

(3.3) \begin{equation} \begin{aligned} &u = \frac {2 \sin (2 \pi x) \cos (\pi y)}{\sqrt {5}}, \quad v = -\frac {4 \cos (2 \pi x) \sin (\pi y)}{\sqrt {5}}. \end{aligned} \end{equation}

In this case, the fundamental period of the system is 1, and the unit cell is given by $\varOmega (1) = \left \{(x,y) \, | \, -({1}/{2}) \leqslant x \leqslant ({1}/{2}), \, 0 \leqslant y \leqslant 1 \right \}.$ In the case of shear flow, since the eigenfunction is independent of $x$ , $\kern2pt \hat{\kern-2pt \mathcal{L}}(k)$ in (2.12) reduces to the Laplace operator when $k = 0$ , yielding $\lambda _{1}(0) = \pi ^{2}$ for all ${\textit{Pe}}$ . In contrast, for cellular flow, the advection term persists at $k = 0$ , causing $\lambda _{1}(0)$ to depend on ${\textit{Pe}}$ . Figure 8 illustrates $\lambda _{1}(0)$ as a function of ${\textit{Pe}}$ . As ${\textit{Pe}}$ increases, $\lambda _{1}(0)$ initially grows rapidly, but the rate of increase slows down after ${\textit{Pe}} \approx 46$ . For example, at ${\textit{Pe}} = 200$ , $\lambda _{1}(0) = 47.6839$ , which is much larger than that in the shear flow case. This suggests that the scalar field under the advection of the cellular flow converges much faster to the slow manifold compared with the shear flow, a result attributed to the velocity component in the transverse direction.

Figure 8. Eigenvalues $\lambda _{n} (0), n=1,2,3$ as a function of ${\textit{Pe}}$ for the cellular flow in the channel domain with flat walls.

We note that for large ${\textit{Pe}}$ , the effective longitudinal diffusivity $\kappa _{\textit{eff}}$ induced by shear flow exceeds that induced by cellular flow. Specifically, $\kappa _{\textit{eff}} \sim {\textit{Pe}}^{2}$ for steady shear flow (Camassa et al. Reference Camassa, Lin and McLaughlin2010a ), whereas our numerics show $\kappa _{\textit{eff}} \sim {\textit{Pe}}^{1/2}$ for the cellular flow defined in (3.3). The same scaling was rigorously proved under periodic $y$ -boundary conditions in McLaughlin (Reference McLaughlin1994). It is a reasonable result since the effective longitudinal diffusivity measures the ability of the flow to spread the scalar in the longitudinal direction, but it does not fully capture the mixing ability of the flow. The eigenvalue, however, serves as a better measure of the flow’s mixing capability.

Figure 9. The first four eigenvalues as functions of the wavenumber. Here $\lambda _{n} (k)$ is periodic in $k$ with period $2\pi$ . Results are shown for (a) ${\textit{Pe}} = 1$ and (b) ${\textit{Pe}} = 50$ .

Similar to the shear flow case $u = \sqrt {2} \cos \pi y$ , the global minimum of the eigenvalue is not always located at $k = 0$ . Figure 9 shows $\lambda _{1}(k)$ as a function of the wavenumber for different values of ${\textit{Pe}}$ . For the ${\textit{Pe}}$ values considered in table 1, the global minimum of $\lambda _{1}(k)$ occurs at $k = \pm \pi$ . For large ${\textit{Pe}}$ , as ${\textit{Pe}}$ increases, $\lambda _{1}(0)$ grows slowly, while the ratio $\min _{k} \lambda _{1}(k) / \lambda _{1}(0)$ decreases.

Table 1. Eigenvalues for the cellular flow case.

3.3. Sinusoidal channel with pressure-driven flow

Next, we consider a more practical flow: the pressure-driven flow in the sinusoidal channel governed by (2.2). The flow is solved for different Reynolds number in the cell domain given by

(3.4) \begin{equation} \begin{aligned} \varOmega (L_{p}) =\left \{ (x, y) | - \frac {L_p}{2}\leqslant x \leqslant \frac {L_p}{2}, 0\leqslant y\leqslant h(x) \right \},\quad h(x) = 1 + A \sin \left (\frac {2\pi x}{L_{p}}\right )\!. \end{aligned} \end{equation}

It is important to note that, for a fixed pressure gradient, the fluid flux at the inlet decreases as the amplitude of the boundary variation increases. This indicates that the flow strength becomes weaker with larger boundary undulations. Figure 10 shows the relationship between the amplitude $A$ and the inlet flux for a fixed pressure gradient of $f_x = -12$ . When $A = 0$ , the flux is 1. In contrast, when $A = 0.8$ , the flux drops to approximately $0.025 \sim 0.026$ , which is about 40 times smaller. Therefore, to ensure that the characteristic velocity remains comparable across different values of $A$ , we adjust the pressure gradient so that the inlet flux is fixed to be 1 for all $A$ .

Figure 10. The averaged fluid flux at the inlet, $({1}/{|\varOmega _{{inlet}}|} )\int _{{inlet}} u (x,\boldsymbol{y}) \mathrm{d}\boldsymbol{y}$ , as a function of $A$ for different Reynolds numbers, where $|\varOmega _{{inlet}}|$ is the area of the inlet.

The constant-flux and constant-pressure-gradient cases are related by a nonlinear rescaling of ${\textit{Re}}$ and ${\textit{Pe}}$ via the pressure–flux relation shown in figure 10: at larger amplitudes, a fixed pressure drop yields a much smaller mean flux, which would correspond to smaller effective ${\textit{Re}}$ and ${\textit{Pe}}$ in the fixed-flux formulation. In other words, results for a constant pressure gradient can be interpreted as lying along the fixed-flux curves at appropriately reduced ${\textit{Re}}$ and ${\textit{Pe}}$ .

From the previous section, we see that $\lambda _1(k)$ may not reach the global minimum at $k = 0$ . However, $\lambda _1(0)$ does not differ significantly from $\min _k \lambda (k)$ and can serve as an estimate for the time scale $t_s$ . Therefore, to reduce computational cost, we focus primarily on $\lambda _1(0)$ in this section for simplicity.

The real part of $\lambda _1(0)$ for different Reynolds numbers, $A$ , and ${\textit{Pe}}$ is shown in figure 11. For $A = 0$ , where the channel boundary is flat, the shear flow is the solution to the Navier–Stokes equation for any ${\textit{Re}}$ . In this case, $\lambda _1(0) = \pi ^2$ for all ${\textit{Pe}}$ . The black curve in figure 11 indicates the contour for $\lambda _1(0) = \pi ^2$ , which separates the parameter regimes where the scalar field converges to the slow manifold faster than in the flat channel and where it converges more slowly.

Figure 11. The real part of $\lambda _{1} (0)$ as a function of ${\textit{Pe}}$ and $A$ . The axis for ${\textit{Pe}}$ is in log scale. The black line indicates the contour line for $\lambda _{1} (0)=\pi ^{2}$ . Results are shown for (a) ${\textit{Re}}=0$ , (b) ${\textit{Re}}=10$ , (c) ${\textit{Re}}=100$ , (d) ${\textit{Re}}=200$ .

For small ${\textit{Pe}}$ , increasing the amplitude $A$ causes $\textrm{Re} \lambda _{1}(0)$ to decrease below $\pi ^2$ . This suggests that, even though the non-flat boundary induces a non-zero transverse velocity component, the scalar may not converge to the slow manifold faster than in the case of pure shear flow. In contrast, for all Reynolds numbers considered here, if ${\textit{Pe}}$ is sufficiently large (above $10^{2}$ ), increasing $A$ leads to an increase in $\textrm{Re} \lambda _{1}(0)$ , making it larger than $\pi ^2$ .

This behaviour is reasonable, as two competing effects are at play. On one hand, the non-flat boundary restricts diffusion compared with the flat case, which reduces the effective diffusivity and tends to slow down convergence. On the other hand, the induced transverse velocity enhances convergence. The relative importance of these effects depends on ${\textit{Pe}}$ . When ${\textit{Pe}}$ is small, diffusion dominates, and the scalar converges more slowly. When ${\textit{Pe}}$ is large, advection becomes dominant, accelerating the convergence to the slow manifold.

If we view the contour as a curve ${\textit{Pe}} = f(A)$ in the $A-{\textit{Pe}}$ plane, its shape depends on the Reynolds number. When ${\textit{Re}} = 0$ , the curve is convex, but it becomes non-convex for larger ${\textit{Re}}$ . For non-zero Reynolds numbers ( ${\textit{Re}} \gt 0$ ), there exists a range of ${\textit{Pe}}$ where the non-flat boundary enhances convergence to the slow manifold for either small or sufficiently large values of $A$ .

We next investigate the asymptotic behaviour in the limit of small $A$ . Owing to symmetry, the expansions for both the effective diffusivity and the eigenvalue contain only even powers of $A$ :

(3.5) \begin{equation} \begin{aligned} &\kappa _{\textit{eff}} = 1 + \frac {{\textit{Pe}}^{2}}{210} + \kappa _{2}A^{2} + \kappa _{4}A^{4} + \mathcal{O}(A^{6}), \lambda _{1}(0) = \pi ^{2} + \lambda _{1,2}A^{2} + \lambda _{1,4}A^{4} + \mathcal{O}(A^{6}). \end{aligned} \end{equation}

We fit the data presented in figure 11 to determine the coefficients in (3.5). Figure 12 shows $\lambda _{1,2}$ and $\kappa _{2}$ as functions of ${\textit{Pe}}$ .

Figure 12. Coefficients in the asymptotic expansion of the effective diffusivity (a) and the eigenvalue (b) in (3.5). In (b) the inset shows a log–log plot of $\lambda _{1,2}$ over a wider range of ${\textit{Pe}}$ . The black dashed line represents ${\textit{Pe}}^{2}$ . The largest relative fitting error across all investigated parameters is $3\,\%$ . Results are shown for (a) $\lambda _{1,2}$ and (b) $\kappa _{2}$ .

Several trends can be observed for $\lambda _{1,2}$ . First, $\lambda _{1,2}$ approaches a limiting value as ${\textit{Pe}}$ becomes large or small, with distinct limits in each regime. Second, for all ${\textit{Re}}$ considered, $\lambda _{1,2}$ is negative at small ${\textit{Pe}}$ and positive at large ${\textit{Pe}}$ . This implies that, for small amplitudes, $\lambda _{1}(0) \lt \pi ^{2}$ when ${\textit{Pe}}$ is small and $\lambda _{1}(0) \gt \pi ^{2}$ when ${\textit{Pe}}$ is large, consistent with figure 11. The ${\textit{Pe}}$ values at which $\lambda _{1,2}$ changes sign are $18.5$ , $26.1$ , $108$ and $183$ for ${\textit{Re}} =$ 0, 10, 100, 200, respectively. For the ${\textit{Re}}$ values studied here, $\lambda _{1,2}$ decreases as ${\textit{Re}}$ increases, suggesting that higher ${\textit{Re}}$ slows convergence toward the slower manifold when the amplitude is small.

A similar trend holds for $\kappa _{2}$ : it is negative at small ${\textit{Pe}}$ and positive at large ${\textit{Pe}}$ . The sign-change ${\textit{Pe}}$ values are $3.99$ , $3.81$ , $3.15$ and $2.93$ for ${\textit{Re}} =$ 0, 10, 100, 200, respectively. For the range of ${\textit{Re}}$ considered, $\kappa _{2}$ increases with ${\textit{Re}}$ . For large ${\textit{Pe}}$ , $\lambda _{1,2}$ scales as ${\textit{Pe}}^{2}$ , which is consistent with the fact that $\kappa _{\textit{eff}}$ scales as ${\textit{Pe}}^{2}$ when $A=0$ . The coefficient in the small-amplitude expansion is therefore expected to exhibit the same scaling behaviour.

The large-amplitude limit is challenging to simulate because maintaining a constant fluid flux requires a very large pressure drop, which in turn necessitates a highly refined mesh for accuracy. Studying this regime would require specialised analytical techniques, such as those in Alexandre, Guérin & Dean (Reference Alexandre, Guérin and Dean2025).

To investigate the asymptotic behaviour of the effective diffusivity $\kappa _{\textit{eff}}$ and the eigenvalue $\lambda _{1}(0)$ in the long-wavelength limit $L_{p} \to \infty$ , we compute their values numerically and present the results in table 2. Several observations can be made. First, neither $\kappa _{\textit{eff}}$ nor $\lambda _{1}(0)$ converges to the flat channel values as $L_p$ increases, which is expected since boundary variations persist even in the long-wavelength limit. Notably, the approach to the limiting values can be non-monotonic in some cases. Second, for small $L_p$ , inertia can increase $\kappa _{\textit{eff}}$ substantially. However, the limiting values of $\kappa _{\textit{eff}}$ depend only on the geometry and are independent of the Reynolds number. This follows from the fact that $\partial _{x} u$ and $v$ scale as $L_{p}^{-1}$ , so inertial effects become negligible as $L_p \to \infty$ . Third, large boundary amplitudes slow the convergence to the limiting values dramatically.

Table 2. Effective diffusivity $\kappa _{\textit{eff}}$ and eigenvalue $\lambda _{1}(0)$ for the domain defined in (3.4) are shown for ${\textit{Pe}} = 50$ with $A = 0.3$ (first two rows) and $A = 0.8$ (last two rows). Columns correspond to $L_p = 1, 2, 5, 10, 20$ . The left block lists $\kappa _{\textit{eff}}$ and the right block lists $\lambda _{1}(0)$ . For reference, a flat-walled channel yields $\kappa _{\textit{eff}} = 1 + {\textit{Pe}}^2/210 \approx 12.905$ and $\lambda _{1}(0) = \pi ^2 \approx 9.870$ .

The dependence of $\kappa _{\textit{eff}}$ on flow and the sinusoidal geometry in this type of domain has been extensively studied in previous work Bouquain et al. (Reference Bouquain, Méheust, Bolster and Davy2012), Richmond et al. (Reference Richmond, Perkins, Scheibe, Lambert and Wood2013). For conciseness, we do not reproduce these results here, but they can be computed from the present framework in a straightforward manner.