Diffusion-driven flows in a nonlinear stratified fluid layer

Abstract

Diffusion-driven flow is a boundary layer flow arising from the interplay of gravity and diffusion in density-stratified fluids when a gravitational field is non-parallel to an impermeable solid boundary. This study investigates diffusion-driven flow within a nonlinearly density-stratified fluid confined between two tilted parallel walls. We introduce an asymptotic expansion inspired by the centre manifold theory, where quantities are expanded in terms of derivatives of the cross-sectional averaged stratified scalar (such as salinity or temperature). This technique provides accurate approximations for velocity, density and pressure fields. Furthermore, we derive an evolution equation describing the cross-sectional averaged stratified scalar. This equation takes the form of the traditional diffusion equation but replaces the constant diffusion coefficient with a positive-definite function dependent on the solution's derivative. Numerical simulations validate the accuracy of our approximations. Our investigation of the effective equation reveals that the density profile depends on a non-dimensional parameter denoted as $\gamma$ representing the flow strength. In the large $\gamma$ limit, the system is approximated by a diffusion process with an augmented diffusion coefficient of $1+\cot ^{2}\theta$, where $\theta$ signifies the inclination angle of the channel domain. This parameter regime is where diffusion-driven flow exhibits its strongest mixing ability. Conversely, in the small $\gamma$ regime, the density field behaves like pure diffusion with distorted isopycnals. Lastly, we show that the classical thin film equation aligns with the results obtained using the proposed expansion in the small $\gamma$ regime but fails to accurately describe the dynamics of the density field for large $\gamma$.

1. Introduction

The density of a fluid is influenced by several factors, including temperature, the concentration of solute and pressure profiles. Typically, these factors lead to a non-uniform distribution of density in the fluid. Density stratified fluids are commonly found in various natural environments, such as lakes, oceans, the Earth's atmosphere, making them a subject of great interest in numerous research endeavours (Linden 1979; Cenedese & Adduce 2008; Magnaudet & Mercier 2020; Camassa et al. 2022; More & Ardekani 2023). In the context of density-stratified fluids, achieving hydrostatic equilibrium is contingent upon aligning the density gradients with the direction of gravitational force. In scenarios involving impermeable boundaries and the modelling of diffusive stratified scalars, a no-flux boundary condition applies, demanding that the density gradient be orthogonal to the boundary's normal vector. If the boundary's normal vector is perpendicular to the gravitational direction, the no-flux condition prevents the density gradient from aligning with gravity, thus disrupting hydrostatic equilibrium. The disruption of hydrostatic equilibrium results in the emergence of a boundary layer flow, a phenomenon termed diffusion-driven flow in the field of physics (Phillips 1970; Wunsch 1970), or mountain and valley (katabatic or anabatic) winds (Prandtl, Oswatitsch & Wieghardt 1942; Oerlemans & Grisogono 2002) in meteorology.

The phenomenon of diffusion-driven flow has garnered significant attention across various fields. First, in the context of oceanography, the presence of salt in seawater leads to density stratification. The continental shelf has a gentle slope extending from the coast to the deeper waters of the ocean. The sloping boundaries induce a mean upwelling velocity along the boundaries to satisfy the no-flux condition. This upwelling flow plays a crucial role in facilitating the vertical exchange of oceanic properties (see, for example, Phillips 1970; Wunsch 1970; Dell & Pratt 2015; Holmes, de Lavergne & McDougall 2019; Drake, Ferrari & Callies 2020). Second, fluids confined within long, narrow fissures can also exhibit density stratification. This stratification can result from non-uniform concentration distributions of stratified scalars or vertical temperature gradients. For example, Earth's geothermal gradient induces convective flows that lead to significant solute dispersion on geological time scales (Woods & Linz 1992; Shaughnessy & Van Gilder 1995; Heitz, Peacock & Stocker 2005). Third, when a wedge-shaped object is immersed in a stratified fluid, the resulting diffusion-driven flow can propel the object forward (Allshouse, Barad & Peacock 2010; Mercier et al. 2014). The flows generated by such immersed wedge-shaped objects have been further investigated in various contexts (Chashechkin & Mitkin 2004; Zagumennyi & Dimitrieva 2016; Chashechkin 2018; Dimitrieva 2019; Levitsky, Dimitrieva & Chashechkin 2019). Fourth, diffusion-driven flow is one of the mechanisms that induce particle attraction and self-assembly in stratified fluids (Camassa et al. 2019; Thomas & Camassa 2023). Fifth, it has been demonstrated that this type of flow can be employed to measure the molecular diffusivity of stratified scalars (Allshouse 2010). Sixth, flows induced by the presence of insulating sloping boundaries play an important role in the layer formation in double-diffusive systems (Linden & Weber 1977).

Despite the significant findings in this field, there are two points that have not been adequately addressed in the existing literature. First, most existing theories primarily pertain to linearly stratified fluids, and there is a scarcity of theoretical investigations into diffusion-driven flow in nonlinearly stratified fluids. Nonlinearly stratified fluids are more prevalent in various natural environments, making theoretical analysis in this context highly applicable. Second, there can exist two different type of scalars in the fluid: the stratified scalar and the passive scalar. The stratified scalar causes non-uniform density distributions in the fluid and drives diffusion-driven flow. In contrast, the passive scalar represents the concentration of a different solute (in some cases, the temperature field) that does not contribute to density variation but is instead passively advected by the fluid flow. It is well-established that fluid flow enhances solute mixing in the fluid (Lin, Thiffeault & Doering 2011; Thiffeault 2012; Aref et al. 2017). Many studies focus on how diffusion-driven flow enhances the dispersion of a passive scalar (Woods & Linz 1992; Ding & McLaughlin 2023) and the corresponding analysis for the stratified scalar is rare. Experiments, such as the one documented in Ding (2022), illustrate that the dispersion of the stratified scalar in a tilting capillary pipe exhibits a higher dispersion rate compared with a vertically oriented pipe. This qualitative demonstration highlights the role of diffusion-driven enhancement at microscales. At diffusion time scales, the concentration of a diffusing passive scalar under the advection of shear flows can be described by an effective diffusion equation with an effective diffusion coefficient (Taylor 1953; Aris 1956). Consequently, it is of interest to establish the evolution equation for the concentration of the stratified scalar. This would enable us to quantitatively describe how diffusion-driven flow enhances the spreading rate of the stratified scalar.

To address these research gaps, this paper delves into the study of an incompressible viscous density stratified fluid layer confined between two infinite parallel walls. Investigating this fundamental domain geometry can enhance our comprehension of more intricate shapes found in real-world scenarios, such as rock fissures and capillary pipes. Furthermore, it provides a versatile framework for comprehending fluid dynamics within confined spaces, which holds substantial practical significance.

Our primary objectives are to derive approximations for the velocity field and density field and to determine the equation governing the dynamics of the fluid density field at diffusion time scales. One conventional approach, such as thin film approximation, involves introducing a small parameter and utilizing power series expansions for all relevant quantities with respect to this parameter. However, we demonstrate that results obtained through thin film approximation are only valid within specific parameter ranges and do not provide uniformly accurate approximations across all parameters. To obtain an approximation that accurately describes density dynamics across a wider parameter range, we adopt an alternative expansion method inspired by centre manifold theory (Mercer & Roberts 1990; Roberts 1996; Ding & McLaughlin 2022a). In this alternative approach, we assume that the density field varies slowly in the longitudinal direction of the channel. Then all quantities are represented in terms of derivatives of the averaged stratified scalar. This innovative approach enhances our ability to achieve accuracy in a broader range of flow scenarios.

This paper is organized as follows. In § 2, we formulate the governing equations for diffusion-driven flow and outline the procedure for non-dimensionalization. Section 3 introduces an asymptotic expansion technique for the velocity, density and pressure fields. Using this approach, we derive the leading-order approximation for the velocity field and establish the evolution equation for the cross-sectional averaged stratified scalar, referred to as the effective equation in subsequent sections. In § 3.3, we conduct numerical simulations of the full governing equations to demonstrate the validity of the asymptotic approximations. Section 4 delves deeper into the properties of the effective equation, providing a more comprehensive understanding of its behaviour in various scenarios. Section 5 documents the results obtained using the thin film approximation and compares them with the expansion proposed in § 3. Finally, in § 6, we summarize our findings and explore potential avenues for future research.

2. Governing equation and non-dimensionalization

This section summarizes the mathematical formulation of the problem and documents the non-dimensionalization procedure for the governing equation.

2.1. Governing equation

Figure 1 sketches two coordinate systems for a tilted parallel-plate channel domain with an inclination angle $\theta$ relative to the horizontal plane, which satisfies $0\leq \theta \leq {{\rm \pi} }/{2}$. In this set-up, the $x_3$-direction is parallel to the direction of gravity. Here $\varOmega =\{ y_{3}| y_3 \in [0,H] \}$ is the cross-section of the channel, where $H$ is the distance between plates. The longitudinal direction of the channel is denoted as the $y_1$-direction, with $y_{1}\in [-L,L]$, where $L$ may be either infinite or a finite number that is significantly larger than $H$. The relation between the laboratory frame coordinates $(x_{1},x_{2},x_{3})$ and the coordinates $(y_{1},y_{2},y_{3})$ is

(2.1)\begin{equation} \begin{bmatrix} y_{1}\\ y_{3}\\ \end{bmatrix} = \begin{bmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\\ \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{3}\\ \end{bmatrix}, \quad y_{2}=x_{2}. \end{equation}

In the $(y_{1},y_{2},y_{3})$ coordinate system, gravity acts along the direction $(-\sin \theta, 0, -\cos \theta )$. Experimental methods described in Allshouse (2010) and Heitz et al. (2005) provide feasible set-ups for this study. Another promising experimental configuration involves using temperature-stratified liquid gallium (Braunsfurth et al. 1997).

Figure 1. This schematic illustrates the set-up for the diffusion-driven flow problem. Tilted impermeable boundaries induce the diffusion-driven flow, while there is no flow in the vertically oriented channel ($\theta ={{\rm \pi} }/{2}$). Pure molecular diffusion governs the scalar distribution in the latter case. As depicted in (a), at $t=0$, the concentration profiles are approximately the same in both channels. However, as time progresses, the dispersion of the scalar becomes more pronounced in the tilted channel compared with the vertically oriented one. In addition, the flow strength decreases as the density gradient decreases.

We consider a scalar $c$ responsible for fluid density stratification, referred to as the stratified scalar in the subsequent context, such as temperature or solute concentration. The density is assumed to be a specified function of the stratified scalar. For example, Abaid et al. (2004) obtained the formula for the density of a sodium chloride solution as $\rho = 0.9971 + 0.00065S + 0.000322(25-T)\ {\rm g}\ {\rm cm}^{-3}$, where salinity $S$ is measured in parts per thousand and temperature $T$ is in degrees Celsius. In cases with a broader range of salinity, the relationship may exhibit nonlinearity. Hall (1924), through empirical data fitting, established $\rho =0.997071+0.00070109 S+ 1.3268\times 10^{-7}S^{2}+3.535\times 10^{-10} S^3 \ {\rm g}\ {\rm cm}^{-3}$ for sodium chloride solutions at $25\,^{\circ}$C. The theoretical method presented in this work can be easily extended to cases with variable fluid viscosity. However, for simplicity, we assume the fluid dynamic viscosity is constant throughout the paper.

The stratified scalar $c$ and the fluid flow $\boldsymbol {u}=(u_{1},u_{2},u_{3})$ in the $\boldsymbol {x}= (x_{1},x_{2},x_{3})$ direction satisfy the incompressible Navier–Stokes equation,

(2.2)\begin{align} \left.\begin{array}{c} \rho (c) ( \partial_{t} u_{i}+ \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} u_{i} ) = \mu \Delta u_{i} - \partial_{x_{i}} p- \rho (c) g \delta_{i3},\quad u_{i} |_{\boldsymbol{x}\in \partial \varOmega}=0,\quad i=1,2,3,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0,\quad \partial_{t}c+ \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} c= \kappa \Delta c,\quad \partial_{\boldsymbol{n}}c |_{\boldsymbol{x}\in\partial \varOmega }=0, \quad c|_{t=0}=c_{I}(\boldsymbol{x}), \end{array}\right\} \end{align}

where $\boldsymbol {n}$ is the outward normal vector of the boundary, $\delta _{ij}$ is the Kronecker delta, $g$ (cm s$^{-2}$) is the acceleration of gravity, $\rho$ (g cm$^{-3})$ is the density, $\mu$ (g cm$^{-1}$ s$^{-1})$ is the dynamic viscosity, $p$ (g cm$^{-1}$ s$^{-2})$ is the pressure and $\kappa$ (cm$^{2}$ s$^{-1})$ is the molecular diffusivity of the stratified scalar.

It is convenient to consider the problem in the $(y_{1},y_{2},y_{3})$ coordinate system. We denote $v_i$ as the velocity component along the $y_{i}$-direction. Since the initial condition and the boundary condition are independent of $y_2$, (2.2) reduces to a two-dimensional problem

(2.3)\begin{equation} \left.\begin{array}{c} \rho ( \partial_{t} v_{1}+ v_{1}\partial_{y_{1}} v_{1} + v_{3}\partial_{y_{3}} v_{1}) =\mu \Delta v_{1} - \partial_{y_{1}} p- \sin \theta g \rho ,\\ \rho ( \partial_{t} v_{3}+ v_{1}\partial_{y_{1}} v_{3} + v_{3}\partial_{y_{3}} v_{3}) =\mu \Delta v_{3} - \partial_{y_{3}} p- \cos \theta g\rho ,\\ \partial_{t}c+ \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} c= \kappa \Delta c, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}=0,\quad \boldsymbol{v} |_{\boldsymbol{y}\in \partial \varOmega }=\boldsymbol{0}, \quad \partial_{\boldsymbol{n}}c |_{\boldsymbol{y}\in\partial \varOmega}=0. \end{array} \right\} \end{equation}

In this study, we assume the initial density profile is a stable density stratification, where the density decreases monotonically as the height increases, namely, $\partial _{y_{1}}\rho (\bar {c}_{I}) \leq 0$. For convenience, we use the overline to represent the cross-sectional average of a quantity and use the tilde to represent the fluctuation, for example, $\overline {c_{I}}(y_{1})= \int _0^1 c_{I}(y_{1},y_{3}) \,\mathrm {d} y_{3}$.

2.2. Non-dimensionalization

As the flow is driven by molecular diffusion, we select the diffusion time scale as the characteristic time, ${H^{2}}/{\kappa }$. Specifically, it refers to the time it takes for solute molecules to diffuse across the cross-sectional area of the channel domain. Utilizing the following change of variables for non-dimensionalization:

(2.4)\begin{equation} \rho_{0}\rho'=\rho , \quad c_{0}c'=c, \quad\frac{H^2}{\kappa} t' =t, \quad H\boldsymbol{y}'=\boldsymbol{y}, \quad U\boldsymbol{v}'=\boldsymbol{v}, \quad \frac{\mu U}{H} p' =p, \end{equation}

then (2.3) becomes

(2.5) \begin{align} \left.\begin{array}{c} \displaystyle \rho' \left( \dfrac{U \kappa}{H^2} \partial_{t'} v_{1}'+ \dfrac{U^{2}}{H} v_{1}'\partial_{y_{1}'} v_{1}' + \dfrac{U^{2}}{H} v_{3}'\partial_{y_{3}'} v_{1}'\right) =\dfrac{ \mu U}{H^2\rho_{0}} \varDelta_{\boldsymbol{y}'} v_{1}' - \dfrac{\mu U}{H^2\rho_{0}}\partial_{y_{1}'} p'- \sin \theta g \rho' , \\ \displaystyle \rho' \left( \dfrac{U \kappa}{H^2} \partial_{t'} v_{3}'+ \dfrac{U^{2}}{H} v_{1}'\partial_{y_{1}'} v_{3}' +\dfrac{U^{2}}{H} v_{3}'\partial_{y_{3}'} v_{3}'\right) =\dfrac{ \mu U}{H^2\rho_{0}} \varDelta_{\boldsymbol{y}'} v_{3}' - \dfrac{\mu U}{H^2\rho_{0}}\partial_{y_{3}'} p'- \cos \theta g\rho' , \\ \displaystyle \dfrac{c_0\kappa}{H^2} \partial_{t'}c'+ \dfrac{U c_{0}}{H} \boldsymbol{v}' \boldsymbol{\cdot} \boldsymbol{\nabla} _{\boldsymbol{y}'}c'= \dfrac{\kappa c_{0}}{H^2} \varDelta_{\boldsymbol{y}'}c',\quad \dfrac{c_{0}}{H} \partial_{\boldsymbol{n}'}c' |_{\boldsymbol{y}'\in \partial \varOmega'}=0. \end{array} \right\} \end{align}

We can drop the primes without confusion and obtain the non-dimensionalized version

(2.6)\begin{align} \left.\begin{array}{c} \displaystyle \rho \left( \dfrac{1}{Sc} \partial_{t} v_{1}+ {Re} v_{1}\partial_{y_{1}} v_{1} + {Re} v_{3}\partial_{y_{3}} v_{1}\right) =\Delta v_{1} - \partial_{y_{1}} p- {Re}\, {Ri}\rho \sin \theta, \\ \displaystyle \rho \left( \dfrac{1}{Sc} \partial_{t} v_{3}+ {Re} v_{1}\partial_{y_{1}} v_{3} +{Re} v_{3}\partial_{y_{3}} v_{3}\right) = \Delta v_{3} - \partial_{y_{3}} p- {Re}\, {Ri}\rho \cos \theta, \\ \partial_{t}c+ {Pe} \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} c= \Delta c, \end{array} \right\} \end{align}

where the non-dimensional parameters are Péclet number ${Pe}= {UH}/{\kappa }$, Reynolds number ${Re}= {\rho _{0} H U}/{\mu }$, Richardson number ${Ri}= {g H}/{U^{2}}$ and Schmidt number ${Sc}= {\mu }/{\rho _{0} \kappa }={Pe}/{Re}$. If the scalar field is the temperature field, then $\kappa$ is the thermal diffusivity and ${Pr}={\mu }/{\rho _{0} \kappa }={Pe}/{Re}$ is the Prandtl number.

We proceed by examining a combination of experimental physical parameters, which can provide us with the order of magnitude of the non-dimensional parameters and assist in our perturbation analysis. The previous experiments primarily focused on diffusion in linearly stratified fluids. It is feasible to adapt these set-ups for experiments involving nonlinearly stratified fluids, making the parameters therein useful as reference points. For instance, in Camassa et al. (2019), a linear density stratification of the fluid was achieved using sodium chloride. The relevant physical parameters in this context are as follows: $g=980\ \textrm {cm}\ \textrm {s}^{-2}$, $\mu _{0}=0.008903\ \textrm {g}\ (\textrm {cm}~\textrm {s})^{-1}$, $\kappa = 1.5\times 10^{-5}\ \textrm {cm}^{2}\ \textrm {s}^{-1}$ (Vitagliano & Lyons 1956), $\rho _{0}= 0.9971\ \textrm {g}\ \textrm {cm}^{-3}$, the vertical density gradient $\partial _{x_{3}}\rho$ is a constant $0.007\ \textrm {g}\ \textrm {cm}^{-4}$. The scaling relation for the characteristic velocity and the physical parameters varies depending on the boundary geometries. For a linear stratified fluid, the characteristic velocity and characteristic boundary layer thickness of steady diffusion-driven flow in a parallel-plate channel can be calculated using the formula from (Phillips 1970; Heitz et al. 2005) as follows:

(2.7a,b)\begin{equation} U =2 \kappa \cot (\theta ) \left( \frac{ g \sin (\theta) \partial_{x_{3}}\rho}{4 \kappa \mu } \right)^{{1}/{4}},\quad H_{b}=\left( \frac{ g \sin (\theta )\partial_{x_{3}}\rho }{4 \kappa \mu } \right)^{-{1}/{4}}. \end{equation}

Using these formulae, we obtain the values $U =0.00164684\ \textrm {cm}\ \textrm {s}^{-1}$, $H_{b}= 0.0182167$ cm for $\theta ={{\rm \pi} }/{4}$. When the gap thickness between two walls is $H=0.1$ cm, this leads to the following non-dimensional parameters:

(2.8a––d)\begin{equation} {Re}=0.0184512,\quad {Pe}=10.9799, \quad {Sc}=595.076, \quad {Re} \, {Ri}=666613. \end{equation}

It is evident that the Reynolds number is small, indicating that viscous effects dominate the flow, and the gravity term is significant in the governing equation.

Last, with the above parameters and $c_{0}=1\ \textrm {mol}\ \textrm {kg}^{-1}$, the relation between $\rho$ and $c$ can be expressed as

(2.9)\begin{equation} \rho (c)=1+0.0410921 c+ 0.000454464 c^2+ 0.000070761 c^3. \end{equation}

3. Asymptotic analysis and numerical simulation

In this section, our goal is to derive approximations for the velocity and density fields and validate them through numerical simulations.

The hydrostatic body force terms in the governing equation contribute only to the hydrostatic pressure field. Since they do not influence the velocity and stratified scalar field, we can absorb them into a modified pressure. For convenience, we define $p=p_{0}+\tilde {p}$, where

(3.1)\begin{equation} \!\!\left.\begin{array}{c} \displaystyle p_{0}={-} {Re}\, {Ri} \left(\sin \theta \int_{0}^{y_{1}}\rho(\bar{c} (y_{1}))\, \mathrm{d} y_{1}+ \cos \theta\left(y_{3}-\dfrac{1}{2}\right)\rho(\bar{c}) \right),\\ \displaystyle \!\!\partial_{y_{1}}p_{0}={-} {Re}\, {Ri} \left(\sin \theta \rho( \bar{c}) + \cos \theta\left(y_{3}-\dfrac{1}{2}\right)\partial_{y_{1}}\rho (\bar{c}) \right) , \quad \partial_{y_{3}}p_{0}={-} {Re}\, {Ri} \cos \theta \rho(\bar{c}). \end{array}\! \right\} \end{equation}

Then the governing equation (2.6) becomes

(3.2)\begin{equation} \left.\begin{array}{c} \begin{aligned} & \rho \left( \dfrac{1}{Sc} \partial_{t} v_{1}+ {Re} v_{1}\partial_{y_{1}} v_{1} + {Re} v_{3}\partial_{y_{3}} v_{1}\right) \\ & \quad =\Delta v_{1} - \partial_{y_{1}} \tilde{p}- {Re}\, {Ri} \left( \sin \theta (\rho - \rho(\bar{c})) -\cos \theta\left(y_{3}-\dfrac{1}{2}\right)\partial_{y_{1}}\rho (\bar{c})\right), \end{aligned}\\ \displaystyle \rho \left( \dfrac{1}{Sc} \partial_{t} v_{3}+ {Re} v_{1}\partial_{y_{1}} v_{3} +{Re} v_{3}\partial_{y_{3}} v_{3}\right) = \Delta v_{3} - \partial_{y_{3}} \tilde{p}- {Re}\, {Ri} \cos \theta (\rho - \rho(\bar{c})), \\ \partial_{t}c+ {Pe} ( v_{1}\partial_{y_{1}}c+v_{3}\partial_{y_{3}}c ) = \Delta c,\quad \partial_{y_{1}}v_{1}+\partial_{y_{3}}v_{3}=0. \end{array} \right\} \end{equation}

3.1. Expansions for slow varying density profile

To address the scenario where the longitudinal length scale of the channel domain significantly exceeds the transverse length scale, our objective is to develop a simplified model that relies solely on the longitudinal variable of the channel domain. To achieve this, we adopt the following ansatz for the velocity and density:

(3.3) \begin{align} \left.\begin{array}{c}\!\! v_{1} (y_{1},y_{3},t) =v_{1,0}(y_{3},t, \bar{c})+v_{1,1} (y_{3},t, \partial_{y_{1}}\bar{c} )\partial_{y_{1}}\bar{c}+v_{1,2} (y_{3},t, \partial_{y_{1}}\bar{c},\partial_{y_{1}}^{2}\bar{c} ) \partial_{y_{1}}^{2}\bar{c}+\cdots, \\ \!\! v_{3}(y_{1},y_{3},t)=v_{3,0}(y_{3},t,\bar{c})+v_{3,1} (y_{3},t, \partial_{y_{1}}\bar{c} )\partial_{y_{1}}\bar{c}+v_{3,2} (\kern0.08em y_{3},t, \partial_{y_{1}}\bar{c},\partial_{y_{1}}^{2}\bar{c} ) \partial_{y_{1}}^{2}\bar{c}+\cdots, \\ \!\! c(y_{1},y_{3},t)=\bar{c}(y_{1},t)+ c_{1} (y_{3},t, \partial_{y_{1}}\bar{c} )\partial_{y_{1}}\bar{c}+c_{2} (\kern0.08em y_{3},t, \partial_{y_{1}}\bar{c},\partial_{y_{1}}^{2}\bar{c} ) \partial_{y_{1}}^{2}\bar{c}+\cdots,\\ \!\! p(y_{1},y_{3},t)=p_{0}(y_{1},y_{3},t)+ p_{1} (y_{3},t, \partial_{y_{1}}\bar{c} )\partial_{y_{1}}\bar{c}+p_{2} (\kern0.08em y_{3},t, \partial_{y_{1}}\bar{c},\partial_{y_{1}}^{2}\bar{c} ) \partial_{y_{1}}^{2}\bar{c}+\cdots. \end{array}\!\!\! \right\} \end{align}

In this expansion, $v_{i,1},c_{1},p_{1}$ are functions of $\partial _{y_{1}}\bar {c}$. Here $v_{j,2},c_{2},p_{2}$ are functions of $\partial _{y_{1}}\bar {c}$ and $\partial _{y_{1}}^{2}\bar {c}$. The $j$th term in the expansion is contingent on the derivative of $\bar {c}$ up to the $j$th order. By definition, $\partial _{y_{1}}\bar {c}$ is independent of $y_{3}$ and is solely a function of $y_{1}$ and $t$. For the specific problem we investigate here, both $v_{i,j}$ and $c_{i}$ are functions of $y_{3}$ and $t$, while their dependence on $y_{1}$ is captured through $\partial _{y_{1}}\bar {c}$. In more general cases, such as the domain with non-flat boundary, $v_{i,j}$ and $c_{i}$ may depend on $y_{1}$.

In our analysis, we consider the derivative of the averaged stratified scalar, $\partial _{y_{1}}^{n}\bar {c}$, as a small parameter within the framework of standard asymptotic calculations. Additionally, we assume that higher-order derivatives exhibit smaller magnitudes compared with lower-order derivatives. This assumption is reasonable for systems where diffusion dominates.

To provide an intuitive justification, we assume that flow effects are negligible, and diffusion is the dominant process in the system. Under these conditions, the scalar field exhibits a self-similarity solution represented as $\bar {c} = \frac {1}{2} \textrm {erfc}({y_{1}}/{2 \sqrt {t}})$, where $\mathrm {erfc}(z) = ({2}/{\sqrt {{\rm \pi} }})\int _z^\infty \textrm {e}^{-t^2} \, \textrm {d} t$. The derivatives of this solution are as follows:

(3.4)\begin{align} \left.\begin{array}{c} \displaystyle \partial_{y_{1}}\bar{c} ={-}\dfrac{{\rm e}^{-{y_{1}^2}/{4 t}}}{2 \sqrt{\rm \pi} \sqrt{t}},\quad \partial_{y_{1}}^{2}\bar{c} ={-} \dfrac{y_{1}\, {\rm e}^{-{y_{1}^2}/{4 t}}}{4 \sqrt{\rm \pi} t^{3/2}},\\ \displaystyle \partial_{y_{1}}^{n}\bar{c}={-}\dfrac{1}{\sqrt{{\rm \pi} }}\left(-\dfrac{1}{2 \sqrt{t}}\right)^n H_{n-1}\left(\dfrac{y_{1}}{2 \sqrt{t}}\right) {\rm e}^{-(\kern0.08em {y_{1}}/{2 \sqrt{t}})^2}, \end{array} \right\} \end{align}

where $H_{n}$ is the Hermite polynomial of degree $n$. In general, as $t$ approaches infinity, we have $\partial _{y_{1}}^{n}\bar {c} = O(t^{-n+{1}/{2}})$. Consequently, as $t\rightarrow \infty$, we observe the hierarchy $\bar {c} \gg \partial _{y_{1}}\bar {c} \gg \partial _{y_{1}}^{2}\bar {c} \gg \cdots$. Due to the orthogonality of the Hermite polynomial, $\partial _{y_{1}}^{n}\bar {c}$ forms a good basis for approximating the function. In addition, it is important to note that in the presence of a diffusion-driven flow, as we will demonstrate in the following sections, the scaling relationship can differ. For instance, in some parameter limit, the self-similarity variable in the solution can be $y_{1}t^{-{1}/{4}}$ for finite $t$. However, in this case, the derivatives of the averaged density still tend to be zero, and the higher-order derivatives become much smaller than the lower-order derivatives as $t\rightarrow \infty$.

The similar form of expansion (3.3) has found applications in various fields, including the modelling of chromatograph and reactors (Balakotaiah, Chang & Smith 1995), thin film fluid flows (Van Dyke 1987; Roberts 1996; Roberts & Li 2006) and shear dispersion of passive scalars (Gill 1967; Mercer & Roberts 1990; Young & Jones 1991; Ding & McLaughlin 2022a). A more rigorous foundation for this expansion can be established through centre manifold theory (Carr 1979; Roberts 1988, 2014, 2015; Aulbach & Wanner 1996, 1999). For in-depth discussions on constructing centre manifolds, we refer interested readers to the cited literature, and we will not delve into the details here.

Another possible asymptotic expansion approach is a power series expansion involving a small parameter, denoted as $\epsilon$, such as $c = \sum _{i=0}^{\infty }c_{i}\epsilon ^{i}$. This approach is commonly used in multiscale analysis (Kondic 2003; Pavliotis & Stuart 2008; Wu & Chen 2014; Ding 2023). In this context, we highlight two advantages of the ansatz provided by (3.3) over the classical power series expansion for the specific problem addressed in this study.

First, in the classical ansatz, the coefficients remain independent of the small parameter $\epsilon$. The quantity is expanded as a linear combination of $\epsilon$ powers. In contrast, in (3.3), the coefficients depend on small quantities, namely, the derivatives of the averaged scalar field. This nonlinear dependency enables us to achieve a more accurate approximation using fewer terms.

Second, when applying the standard multiscale analysis method to this problem, the results depend on the scaling relation between parameters. In the thin film limit, as we will discuss in § 5, the scalar field can be approximated as $c=c_{0} (y_{1},t) +\epsilon c_{2}(y_{1},y_{3},t)$, where $c_{0}$ is the solution to a diffusion equation with a diffusivity of unity. In this scaling relation, the diffusion process is dominant, and the diffusion-driven flow does not have a first-order contribution. For cases where the diffusion-driven flow significantly enhances the dispersion of the stratified agent, an asymptotic analysis using a different scaling relation for physical parameters become necessary. In some cases, selecting the most appropriate characteristic parameter can be challenging. Additionally, since certain scales in the problem are time-dependent, the choice of scale may lead to time-dependent small parameters, which can complicate the development of a model that is uniformly valid across a wide range of parameter regimes.

In § 5, we will provide detailed explanations of the thin film approximation and discuss the differences between these two approaches.

3.2. Leading-order term in asymptotic expansion

The expansion (3.3) suggests that the deviation of the stratified scalar from its cross-sectional average is small. Consequently, this motivates us to examine the evolution equation for the cross-sectional averaged stratified scalar $\bar {c}$. Taking the cross-sectional average on both sides of the advection–diffusion equation (3.2) and utilizing the incompressibility condition and non-flux boundary condition yields

(3.5)\begin{equation} \partial_{t} \bar{c}+ {Pe}\partial_{y_{1}} \overline{v_{1}c}= \partial_{y_{1}}^{2}\bar{c}. \end{equation}

Therefore, to obtain the leading-order approximation for the cross-sectional averaged stratified scalar, we need to compute $v_{1,0}$, $v_{1,1}$ and $c_{1}$ in the expansion (3.3).

There are some observations that can simplify the calculation of asymptotic expansion. Since the flow is induced via the density gradient, the fluid flow vanishes as $\partial _{y_{1}}\bar {c}$ vanishes. Therefore, $v_{1,0}=v_{3,0}=0$. Substituting the expansion (3.3) into the continuity equation and noticing that $\partial _{y_{1}} ( v_{1,1}\partial _{y_{1}}\bar {c} )= O( \partial _{y_{1}}^{2} \bar {c})$, we obtain

(3.6)\begin{equation} \partial_{y_{1}}v_{1,0}+ \partial_{y_{3}}v_{3,1}=0. \end{equation}

Then to satisfy the no-slip boundary condition, $v_{3,1}$ must be zero.

Collecting the terms that are comparable to $\partial _{y_{1}}\bar {c}$ yields the following equation:

(3.7)\begin{equation} \left.\begin{array}{c} \displaystyle \dfrac{1}{Sc}\rho(\bar{c})\partial_{t}v_{1,1}= \partial_{y_{3}}^2v_{1,1}- {Re}\, {Ri} \partial_{c}\rho (\bar{c}) \left( c_{1} \sin \theta -\cos \theta\left(y_{3}-\dfrac{1}{2}\right)\right),\\ \partial_{t}c_{1}+ {Pe} v_{1,1} \partial_{y_{1}}\bar{c}=\partial_{y_{3}}^{2} c_{1}, \quad \partial_{y_3}c_{1}|_{y_{3}=0,1 }=0,\quad v_{1,1}|_{y_{3}=0,1 }=0. \end{array} \right\} \end{equation}

Previous studies (Kistovich & Chashechkin 1993; Ding & McLaughlin 2023) have uncovered fascinating dynamics in the time-dependent solution, particularly during short time scales when the initial velocity field is zero. However, these transient dynamics decay exponentially, and the solution converges to a steady-state either at the diffusion time scale or the viscosity time scale, depending on which of the two is larger. In many cases, the diffusion time scale significantly exceeds or is comparable to the viscosity time scale. For instance, in solute–liquid systems, the diffusivity of the solute molecule typically falls within the range of $10^{-5}\ \textrm {cm}^{2}\ \textrm {s}^{-1}$, while the kinematic viscosity of the liquid is around $10^{-2}\ \textrm {cm}^{2}\ \textrm {s}^{-1}$, resulting in a Schmidt number (${Sc}$) of $10^{3}$. In the case of temperature-stratified systems, the thermal diffusivity typically registers at approximately $10^{-3}\ \textrm {cm}^{2}\ \textrm {s}^{-1}$, with ${Pr}=10$. In both cases, the diffusion time scale is much larger than the viscosity time scale. Therefore, given our specific interest in approximations at the diffusion time scale and to streamline our analysis without compromising accuracy, our primary focus will be on the steady-state solution of the aforementioned equation. By neglecting the time derivatives in (3.7), we arrive at the following non-homogeneous linear system for analysis:

(3.8)\begin{equation} \left.\begin{array}{c} 0= \partial_{y_{3}}^2v_{1,1}- {Re}\, {Ri}\partial_{c}\rho (\bar{c}) (\sin \theta c_{1} -\cos \theta (y_{3}-\tfrac{1}{2})), \quad v_{1,1} |_{y_{3}=0,1 }=0, \\ {Pe} v_{1,1} \partial_{y_{1}}\bar{c}=\partial_{y_{3}}^{2} c_{1}, \quad \partial_{y_3}c_{1} |_{y_{3}=0,1 }=0. \end{array} \right\} \end{equation}

Differencing the above equation twice with respect to $y_{3}$ yields

(3.9)\begin{equation} \left.\begin{array}{c} 0= \partial_{y_{3}}^4v_{1,1}- {Re}\, {Ri} \sin \theta \partial_{c}\rho \partial_{y_{3}}^{2} c_{1}, \quad {Pe} \partial_{y_{1}}\bar{c} \partial_{y_{3}}^{2}v_{1,1}=\partial_{y_{3}}^{4} c_{1}, \\ \partial_{y_{3}}^{3}v_{1,1} |_{y_{3}=0,1 }={-}{Re}\, {Ri} \partial_{c}\rho (\bar{c}) \cos \theta. \end{array}\right\} \end{equation}

Then we can decouple the scalar field and velocity as

(3.10)\begin{equation} \left.\begin{array}{c} 0=\partial_{y_{3}}^4v_{1,1} -{Re}\,{Pe}\,{Ri}\sin \theta \partial_{y_{1}}\rho(\bar{c}) v_{1,1}, \quad v_{1,1} |_{y_{3}=0,1 }=0,\\ \partial_{y_{3}}^{3}v_{1,1} |_{y_{3}=0,1 }={-}{Re}\,{Ri} \partial_{c}\rho (\bar{c}) \cos \theta, \\ \partial_{y_{3}}^{4} c_{1}= {Re}\,{Pe}\,{Ri} \partial_{y_{1}}\rho(\bar{c}) ( \sin \theta c_{1} -\cos \theta(y_{3}-\tfrac{1}{2}) ) , \\ \partial_{y_3}c_{1} |_{y_{3}=0,1 }=0, \quad \partial_{y_2}^{2}c_{1} |_{y_{3}=0,1 }=0, \end{array} \right\} \end{equation}

where we have used $\partial _{y_{1}}\rho (\bar {c})=\partial _{c}\rho (\bar {c}) \partial _{y_{1}} \bar {c}$ to simplify the equation. We can solve the equation and express $v_{1,1}$, $c_{1}$ in terms of $\partial _{y_{1}}\bar {c}$ as follows:

(3.11) \begin{align} \left.\begin{array}{c} \displaystyle v_{1,1}= \dfrac{2 \gamma \cot (\theta )}{{Pe}\partial_{y_{1}} \bar{c}} \dfrac{ \sin (\gamma y_{3}) \sinh (\gamma (1- y_{3}))-\sin (\gamma (1- y_{3})) \sinh (\gamma y_{3}) }{ \sin (\gamma )+\sinh (\gamma )},\\ \displaystyle c_{1}=\cot \theta \left( y_{3}- \dfrac{1}{2}-\dfrac{ \cos (\gamma(1- y_{3})) \cosh (\gamma y_{3})-\cos (\gamma y_{3}) \cosh (\gamma(1- y_{3}))}{\gamma (\sin (\gamma)+\sinh (\gamma))} \right), \end{array} \right\} \end{align}

where $\gamma = ({1}/{\sqrt {2}})( -{Re}\, {Pe} \,{Ri} \sin \theta \partial _{y_{1}}\rho (\bar {c}))^{{1}/{4}}$. Recall our assumption of the stable density stratification, $\partial _{y_{1}}\rho (\bar {c})\leq 0$, which implies that $\gamma$ is a real non-negative number. In the case of linear density stratification, the equation above aligns with the steady solution previously presented in Phillips (1970) and Heitz et al. (2005).

As depicted in figure 2, both the normalized velocity $v_{1} \gamma ^{-1}\approx v_{1,1}\partial _{y_{1}}\bar {c}\gamma ^{-1}$ and $c_{1}$ are tightly confined within a narrow region near the boundary, especially when $\gamma$ is large. The velocity is positive near $y_{3}=0$ indicating a upwelling flow near the upward facing boundary. Notably, the normalized velocity exhibits nearly uniform magnitude across different values of $\gamma$. Thus, $\gamma ^{-1}$ serves as an effective indicator of the boundary layer's thickness, and $\gamma$ can be considered the characteristic velocity of the system. Moreover, the functions illustrated in figure 2 display an inherent odd symmetry about $y_{3}=\frac {1}{2}$. This symmetry suggests that the leading-order approximations for both the velocity and density fields also possess odd symmetry with respect to $y_{3}=\frac {1}{2}$. As a result, their cross-sectional averages inherently amount to zero.

Figure 2. (a) Normalized leading-order approximation of the longitudinal velocity ${v_{1,1}\partial _{y_{1}}\rho (\bar {c})}/{\gamma }$ for various parameter $\gamma$. (b) Here $c_{1}$ for different $\gamma$.

As $\partial _{y_{1}}\bar {c}\rightarrow 0$, we have

(3.12)\begin{align} v_{1,1}&={-}\frac{{Re}\, {Ri} y_{3}(y_{3}-1) (2 y_{3}-1) \cos (\theta )\partial_{c}\rho}{12 }\nonumber\\ &\quad -\frac{ {Pe} \,{Re}^2 {Ri}^{2}y_{3} (2 (y_{3} (2 y_{3}-7)+7) y_{3}^4-7 y_{3}+3) \sin (2 \theta ) \partial_{y_{1}}\rho\partial_{c}\rho}{40\,320}+O( ( \partial_{y_{1}}\bar{c} )^{2}), \end{align}

and

(3.13)\begin{equation} c_{1}=\frac{{Pe}\, {Re}\, {Ri} (1-2 y_{3}^3 (3 y_{3} (2 y_{3}-5)+10)) \cos (\theta ) \partial_{y_{1}}\rho}{1440}+O( ( \partial_{y_{1}}\bar{c} )^{2}). \end{equation}

Therefore $v_{1,1}\sim O(1)$, but $c_{1}\sim O(\partial _{y_{1}}\bar {c})$. The leading-order approximation of the longitudinal velocity and the density field are given by $v_{1}\sim v_{1,1}\partial _{y_{1}}\bar {c} = O (\partial _{y_{1}}\bar {c})$ and $c \sim c_{1} \partial _{y_{1}}\bar {c}= O( ( \partial _{y_{1}}\bar {c} )^{2} )$, which confirms the ansatz of the asymptotic expansion (3.3).

Substituting the expansion of the velocity field and scalar field to (3.5) and noticing that $\bar {v}_{1,1}=0$, the leading-order approximation reads

(3.14)\begin{equation} \partial_{t}\bar{c}+ {Pe} \partial_{y_{1}} ( \overline{v_{1,1}c_{1}} ( \partial_{y_{1}}\bar{c} )^{2} ) = \partial_{y_{1}}^{2}\bar{c}. \end{equation}

Computing the average yields

(3.15)\begin{equation} \overline{v_{1,1}c_{1}}=\frac{-\cot ^2(\theta ) }{{Pe} \partial_{y_{1}}\bar{c}} \left( \frac{\sin (\gamma ) \sinh (\gamma )}{(\sin (\gamma )+\sinh (\gamma ))^2}+\frac{5 (\cos (\gamma )-\cosh (\gamma ))}{2 \gamma (\sin (\gamma )+\sinh (\gamma ))}+1 \right). \end{equation}

The effective equation (3.14) for the averaged stratified scalar becomes

(3.16)\begin{equation} \left.\begin{array}{c} \partial_{t}\bar{c}= \partial_{y_{1}} ( \kappa_{eff} \partial_{y_{1}}\bar{c}),\\ \displaystyle \kappa_{eff}=1 +\cot ^2(\theta ) \left( \dfrac{\sin (\gamma ) \sinh (\gamma )}{(\sin (\gamma )+\sinh (\gamma ))^2}+\dfrac{5 (\cos (\gamma )-\cosh (\gamma ))}{2 \gamma (\sin (\gamma )+\sinh (\gamma ))}+1 \right). \end{array} \right\} \end{equation}

This equation can be considered as an generalization of the diffusion equation with the constant diffusion coefficient replaced by a positive-definite function that depends on the derivative of the solution. This equation belongs to the family of equations that takes the form $\partial _{t}w=f (\partial _{x}w) \partial _{x}^{2}w$, which occurs in the nonlinear theory of flows in porous media and also governs the motion of a nonlinear viscoplastic medium (see p. 342 in Polyanin & Zaitsev (2012)).

To complete the leading-order approximation for the whole system, we next consider the approximation for the transverse velocity and pressure. Once the leading-order approximation the longitudinal velocity is known, we can calculate the transverse velocity via the continuity equation. Substituting the expansion (3.3) into the continuity equation and collecting the terms that are comparable to $\partial _{y_{1}}^{2}\bar {c}$ yields

(3.17)\begin{equation} \partial_{y_{1}}(v_{1,1} \partial_{y_{1}}\bar{c})+\partial_{y_{3}}v_{3,2} \partial_{y_{1}}^{2}\bar{c}=0. \end{equation}

Using no-slip boundary condition and the fact that $\bar {v}_{1,1}=0$ yields the expression

(3.18)\begin{align} v_{3,2}&={-}\frac{\displaystyle \int_0^{y_{3}} \partial_{y_{1}}( v_{1,1} \partial_{y_{1}}\bar{c} )\, \mathrm{d} y_{3} }{\partial_{y_{1}}^{2} \bar{c}}\nonumber\\ &=\frac{\cot (\theta )\gamma } {2{Pe} (\sin (\gamma )+\sinh (\gamma ))^2 \partial_{y_{1}}\rho } (\sin (\gamma(1 - y)) \sinh (\gamma y) (y \sin (\gamma )\nonumber\\ &\quad +(y-1) \sinh (\gamma ))-\sin (\gamma y) \sinh (\gamma (1- y)) ((y-1) \sin (\gamma )+y \sinh (\gamma ))). \end{align}

In the limit $\partial _{y_{1}}\bar {c}\rightarrow 0$, we have

(3.19)\begin{align} v_{3,2}&=\frac{1}{24} {Re} \,{Ri} (y_{3}-1)^2 y_{3}^2 \cos (\theta )\nonumber\\ &\quad +\frac{{Pe} \,{Re}^2 {Ri}^2 y_{3}^2 ( y_{3}-1 )^{2} (3 y_{3}^4-6 y_{3}^3-y_{3}^2+4 y_{3}+9) \sin (2 \theta ) \partial_{y_{1}}\bar{\rho}}{120\,960}\nonumber\\ &\quad + O (( \partial_{y_{1}}\bar{c} )^{2}). \end{align}

Therefore, $v_{3} \sim v_{3,2}\partial _{y_{1}}^{2}\bar {c} = O (\partial _{y_{1}}^{2} \bar {c})$, consistent with the assumption regarding the order of magnitude in the expansion given by (3.3).

The $p_{1}$ in the expansion of pressure provided in (3.3) satisfies

(3.20)\begin{equation} -\partial_{y_{3}}p_{1} - {Re}\, {Ri} \cos \theta c_{1}\partial_{c}\rho(\bar{c}) =0. \end{equation}

The solution is

(3.21)\begin{align} p_{1}& ={-}\cot (\theta ) \cos (\theta ) {Ri}{Re}\partial_{c}\rho \left(\frac{1}{12}+ \frac{1}{2} (y-1) y \right. \nonumber\\ &\quad +\left. \frac{\sin (\gamma -\gamma y) \cosh (\gamma y)-\cos (\gamma -\gamma y) \sinh (\gamma y)-\cos (\gamma y) \sinh (\gamma -\gamma y)+\sin (\gamma y) \cosh (\gamma -\gamma y)}{2 \gamma ^2 (\sin (\gamma )+\sinh (\gamma ))} \right). \end{align}

In the limit $\partial _{y_{1}}\bar {c}\rightarrow 0$, we have

(3.22)\begin{equation} p_{1}= \frac{ \cos^{2} \theta {Pe} \,{Re}^{2} {Ri}^{2} (28 y_{3}^6-84 y_{3}^5+70 y_{3}^4-14 y_{3}+3) \partial_{y_{1}}\rho \partial_{c}\rho}{20\,160} +O( ( \partial_{y_{1}}\bar{c} )^{2}). \end{equation}

So far, we have derived the leading-order approximation for the entire system. The numerical simulations in the next section demonstrate that the leading-order approximation has decent accuracy in many scenarios. However, in cases of highly nonlinear scalar distribution, or where the density depends highly nonlinearly on the stratified scalar, higher-order terms are required to achieve desirable accuracy. Similar cases can be found in the previous literature. For example, Grayer et al. (2021) studied a stably stratified square cavity subjected to horizontal oscillations, which despite not being purely diffusion-driven, showed that approximations become less regular as viscosity decreases and how higher-order terms are needed.

3.3. Numerical simulation

In this subsection, we perform simulations of the full governing equation (2.6) to demonstrate both the accuracy and validity of our asymptotic approximation and the effective equation (3.16). The details of the numerical methods are documented in the Appendix. The computational domain is defined as $\{ (y_{1},y_{3}) | y_{1} \in [-5,5], y_{3}\in [0,1]\}$. The no-slip boundary condition is imposed for the velocity field, the no-flux boundary condition is imposed for the stratified scalar. The density and stratified scalar relation (2.9) is used for all test cases in this section.

In the first numerical simulation, the initial condition is

(3.23a,b)\begin{equation} c_{I} = \frac{1}{2}\mathrm{erfc}(y_{1}) -\frac{{\rm e}^{{-}x^2}}{\sqrt{{\rm \pi} }} \cos \theta \left( y_{3}-\frac{1}{2} \right) ,\quad v_{1}=v_{3}=0, \end{equation}

and the dimensionless parameters are

(3.24a––e)\begin{equation} {Re}=1,\quad {Ri}\,{Re}=12\ 500,\quad {Pe}=40, \quad {Sc}=1,\quad \theta=\frac{\rm \pi}{4}. \end{equation}

The stratified scalar and velocity at $t=0.04$ are depicted in figures 3(a i) and 3(b i), respectively. Since density gradient is large around $y_{1}=0$, the flow is also localized in that regime. Figures 3(a iii) and 3(b iii) display the density field and velocity field obtained by the simulation at $t=2$. As time increases, the density field becomes smoother and a large fluid circulation is formed in this whole domain.

Figure 3. In (a), we display the density field at three time instances: $t=0.04$ (i), $t=0.1$ (ii) and $t=2$ (iii). Panel (b) shows a pseudocolour plot illustrating the magnitude of the velocity field at the same time points. The initial condition and dimensionless parameter are provided in (3.23a,b) and (3.24a–e), respectively.

To proceed, we simulate the effective equation (3.16) with the same parameters (3.24a–e) for the comparison with the simulation of the governing equation (2.6). Figure 4(a) offers a comprehensive comparison. Here, we superimpose the solution of the effective equation with the cross-sectional averaged stratified scalar obtained through simulation of the complete governing equation (2.6) at the time instance $t=2$. The remarkable alignment between these curves serves to underscore the accuracy of the effective equation as an adept approximation for the full system. Additionally, we include the solution corresponding to the pure diffusion equation (representing the case in the absence of fluid flow) within the same visual representation. Notably, the contrast between the pure diffusion solution and the effective equation's solution demonstrates a visible enhancement of fluid mixing through diffusion-driven flow.

Figure 4. In (a), we present the solution of the effective equation (3.16) (represented by the red solid curve), the cross-sectional averaged stratified scalar derived from simulating the complete governing equation (illustrated by the blue dashed curve) and the solution to the diffusion equation (depicted by the black dotted curve) at the time instant $t=2$ with parameters in (3.24a–e). In (b), we showcase the temporal evolution of the relative difference between the solution of the effective equation and the cross-sectional averaged stratified scalar obtained through the full simulation. This relative difference is defined as the maximum deviation between the functions across the domain divided by the maximum value of the averaged stratified scalar obtained from the full simulation.

Shifting our attention to figure 4(b), we delve into the temporal dynamics of the relative difference between the effective equation's solution and the full system's solution. At the inception ($t=0$), this disparity is nil, attributed to identical initial conditions. During the initial stages, the relative difference magnifies as the system has yet to transition into the regime well-captured by the effective equation. Nevertheless, even during this phase, the maximum relative difference remains at approximately 0.0056. As the system approaches an asymptotic state, which occurs after the diffusion time scale ($t=1$), the relative difference diminishes, stabilizing at around 0.001 until the simulation's culmination. It is worth noting that due to the assumption of the asymptotic analysis, the effective equation is valid for small density gradient $\partial _{y_{1}}\bar {c} \ll 1$. In this numerical test case, we have $\max _{y_{1},t} \partial _{y_{1}}\bar {c}= {1}/{\sqrt {{\rm \pi} }}\approx 0.56419$, indicating that the parameter regime for the effective equation to reach a good approximation is larger than previously thought.

When deriving the effective equation in the previous section, we neglected the time derivative in (3.7). To verify the validity of this assumption, we conducted simulations with the same initial condition as in (3.23a,b) and parameters in (3.24a–e), but with different values of ${Sc}=10,100$. Figure 5 presents the relative differences in the scalar field, averaged scalar field, and velocity field between the solution for ${Sc}=10, 100$ and the solution for ${Sc}=1$. The relative differences for all three quantities are large at the early stage but quickly decreases as time increases. At the diffusion time scale, the relative difference in scalar field is around $10^{-4}$ and the relative difference in velocity field is around $10^{-2}$, which approach the maximum accuracy achievable with the numerical scheme and resolution used in the simulation. This observation suggest that the value of ${Sc}$ does not have a significant impact on the scalar field and velocity field at the diffusion time scale, justifying the neglect of the time derivative in (3.7) for ${Sc}\geq 1$.

Figure 5. Comparison of the simulations with the initial condition in (3.23a,b), dimensionless parameter in (3.24a–e) and different ${Sc}=1,10,100$: (a) difference in $\bar{c}$; (b) difference in scalar fields; (c) difference in velocity fields. The red curve depicts the relative difference between the solution for ${Sc}=1$ and the solution for ${Sc}=10$. The blue curve depicts the relative difference between the solution for ${Sc}=1$ and the solution for ${Sc}=100$.

To demonstrate the validity of the effective equation across a wide range of parameters, we conducted a simulation with the same initial condition (3.23a,b) and a different set of non-dimensional parameters,

(3.25a––e)\begin{equation} {Re}=0.1,\quad {Ri}\,{Re}=100\ 000,\quad {Pe}=100, \quad {Sc}=1,\quad \theta=\frac{\rm \pi}{4}. \end{equation}

The velocity field results are presented in figure 6(a) and figure 6(b i). Comparing figure 3 with figure 6, we make three observations: first, as ${Ri}$ increases, the velocity field exhibits more intricate structures during the initial stages. A distinctive ‘$S$’ shape curve forms in the pseudocolour plot of the velocity field around $t=0.08$ and persists in subsequent time instants. Second, the flow is more confined near the boundary, which is consistent with the observation we made from figure 2. Third, at $t=1$, figure 3(b) illustrates that the velocity field is antisymmetric with respect to $y_{3}=\frac {1}{2}$. However, in figure 6(b), at $t=1$, the velocity field is not antisymmetric with respect to $y_{3}=\frac {1}{2}$ near the two ends of the domain $y_{1}=\pm 5$.

Figure 6. In (a), a pseudocolour plot illustrates the magnitude of the velocity field with parameters in (3.25a–e) at three time instances: $t=0.03$ (i), $t=0.04$ (ii) and $t=0.08$ (iii). In (b), the subpanel (i) displays the velocity field at $t=1$; (ii) showcases the velocity field approximation at $t=1$, calculated using (3.11) and (3.19); (iii) depicts the relative approximation error of the speed.

Figure 6(b ii,iii) depict the approximation of the velocity field and the corresponding relative error at $t=1$. While the error remains relatively small within the interior of the domain, it increases notably near the corners and becomes significantly large in the narrow regions adjacent to the left- and right-hand boundaries, with the layer thickness approximately 0.05. This discrepancy stems from two main reasons. First, the leading-order approximation of velocity exhibits antisymmetry across the entire domain, whereas the actual velocity field does not possess this property. Second, the asymptotic approximation outlined in the previous section holds true for an infinite domain. However, in a confined domain, the no-slip boundary condition for the velocity field at the end of the domain ($y_{3}=\pm 5$) is required but was not enforced during the derivation of the asymptotic approximation. While the no-flux boundary condition for the stratified scalar ensures a longitudinal velocity component of zero, the transverse velocity could deviate from zero where it should be, resulting in a significant error near the boundary. Achieving a more uniform approximation over the domain necessitates incorporating boundary layer corrections near the boundary into the velocity approximation.

One might assume that the substantial error in velocity approximation near the left- and right-hand ends of the domain could significantly compromise the accuracy of other approximations. However, as depicted in figure 7(a), the solution of the effective equation perfectly overlaps with the averaged scalar field obtained from the full simulation. This alignment is further confirmed in figure 7(b), illustrating the temporal dynamics of the relative difference, which remains small throughout the entire time interval. We believe this is due to two possible reasons: first, the error in the velocity field approximation is localized near the ends of the domain. Since the domain is long and narrow, the boundary layer constitutes a relatively small portion of the overall domain, consequently having a limited impact on the entire dynamic process. For a domain with comparable length scales in both directions, the error in the velocity field approximation may lead to a substantial error in the scalar field approximation. Second, only the longitudinal component of the velocity contributes to the scalar transport in the leading-order approximation. Therefore, the error in the transverse component of velocity does not undermine the accuracy of the effective equation.

Figure 7. In (a), we present the solution of the effective equation (3.16) (represented by the red solid curve) with parameters in (3.25a–e), the cross-sectional averaged stratified scalar derived from simulating the complete governing equation (illustrated by the blue dashed curve) and the solution to the diffusion equation (depicted by the black dotted curve) at the time instant $t=1$. Panel (b) showcases the temporal evolution of the relative difference between the solution of the effective equation and the cross-sectional averaged stratified scalar obtained through the full simulation.

Last, we consider a staircase-like profile as the initial condition, a feature commonly encountered in certain oceanographic scenarios

(3.26) \begin{align} \left.\begin{array}{@{}c@{}} \displaystyle c_{I} =f (y_{1})+ \partial_{y_{1}}f(y_{1}) \cos \theta \left( y_{3}-\frac{1}{2} \right), \quad f(y_{1})= \frac{ \mathrm{erfc}(y_{1}+\tfrac{5}{2})+\mathrm{erfc}(y_{1}-\tfrac{5}{2}) }{4},\\ v_{1} = v_{3}=0. \end{array}\right\} \end{align}

The non-dimensional parameters are provided in (2.8a–d), and $\theta ={{\rm \pi} }/{4}$.

The velocity field results are depicted in figure 8. At the early phase ($t=0.01$), multiple convection cells emerge. The adjacent cells have differing orientations. These cells swiftly amalgamate into two larger convection cells with the same anticlockwise orientation, located where the initial profile exhibits significant gradients. As the density distribution becomes smoother, around $t=0.2$, the merging of the two convection cells initiates. By $t=1$, only a single convection cell remains. The velocity distribution closely resembles that of simulations with differing initial conditions. Figure 8(b iii) illustrates the disparity between the solutions of the effective equation and the cross-sectional averaged stratified scalar obtained through the full simulation. The small relative difference again demonstrates the validity of the effective equation.

Figure 8. A pseudocolour plot illustrates the magnitude of the velocity field with the initial condition (3.26) and parameters in (2.8a–d) at different time instances. Subpanel (b iii) showcases the temporal evolution of the relative difference between the solution of the effective equation and the cross-sectional averaged stratified scalar obtained through the full simulation.

4. Analysis of the effective equation

After confirming the accuracy of the effective equation (3.16) for the stratified scalar through numerical simulations, this section delves deeper into its properties.

The behaviour of (3.16) is primarily governed by $\kappa _{eff}$. Its asymptotic expansions for small and large values of $\gamma$ are provided as follows:

(4.1)\begin{align} \left.\begin{array}{c} \displaystyle \kappa_{eff}=1 \,{+}\,\cot ^2(\theta ) \left(\! \left(\! 1\,{-}\,\dfrac{5}{2 \gamma }\right)\,{+}\,\dfrac{(2 \gamma \,{+}\,5) \sin (\gamma )\,{+}\,5 \cos (\gamma )}{\gamma \,{\rm e}^{\gamma }}\,{+}\, O ({\rm e}^{{-}2\gamma})\right)\!,\,\,\, \gamma\rightarrow \infty, \\ \displaystyle \kappa_{eff}=1 +\cot ^2(\theta ) \left( \dfrac{\gamma ^8}{22680}-\dfrac{2879 \gamma ^{12}}{4086482400}+ O ( \gamma^{16} )\! \right)\! ,\,\,\, \gamma\rightarrow 0. \end{array}\!\! \right\} \end{align}

Figure 9 shows the graph of $({\kappa _{eff}-1})/{\cot ^{2} (\theta )}$ and its approximations as a function of $\gamma$. By examining both the graph and the asymptotic approximations, we can conclude that for $\gamma \geq 0$:

(4.2)\begin{equation} 1\leq \kappa_{eff} (\gamma) < 1+ \cot^{2} (\theta). \end{equation}

Figure 9. The red solid curve represents $({\kappa _{eff}-1})/{\cot ^{2} (\theta )}$ as a function of the non-dimensional parameter $\gamma = ({1}/{\sqrt {2}})( -{Re}\, {Pe}\, {Ri}\sin \theta \partial _{y_{1}}\bar {\rho } )^{{1}/{4}}$. For small values of $\gamma$, the corresponding approximation is ${\gamma ^8}/{22\,680}$, which is depicted by the blue dashed curve. Conversely, for large values of $\gamma$, the corresponding approximation is $1- {5}/{2\gamma }$, illustrated by the black dotted curve.

For large values of $\gamma$, (3.16) becomes a diffusion equation with an enhanced diffusion coefficient:

(4.3)\begin{equation} \partial_{t}\bar{c}=( 1+\cot^{2}\theta ) \partial_{y_{1}}^{2}\bar{c} . \end{equation}

This resembles the scenario of a passive scalar governed by an advection–diffusion equation with a prescribed velocity field. In the context of the channel domain and at diffusion time scales, the advection–diffusion equation can be effectively simplified to a diffusion equation with an enhanced effective diffusion coefficient (see Chatwin (1970), Smith (1982), Young & Jones (1991), Ding et al. (2021) and Ding & McLaughlin (2022b) for related work).

However, it is important to note that $\gamma$ is not a constant – it depends on the density gradient. Therefore, the diffusion equation (4.3) cannot provide a uniform approximation of the effective equation for all time instances. As time progresses, the stratified scalar becomes more homogeneous across the domain. As the density gradient decreases, $\gamma$ decreases, and the approximation given by (4.3) becomes less accurate.

To illustrate this, we conducted simulations with the parameters provided in (3.25a–e) as an example. Figure 10(a) shows the maximum value of $\gamma$ across the domain as a function of time. As we expect, it is a decreasing function. Figure 10(b) demonstrates that the solution to the diffusion equation (4.3) reasonably approximates the solution to the effective equation (3.16) initially. However, as time progresses and the density gradient decreases, $\max _{y_{1}}\gamma$ also decreases. At $t=5$, the difference between the solution of (4.3) and that of (3.16) becomes more pronounced. This example also demonstrates that a single diffusion equation is not enough to accurately describes the dynamics of the averaged stratified scalar at all different time scale.

Figure 10. Panel (a) shows $\max _{y_{1}}\gamma$ as a function of time in the simulation with the parameters provided in (3.25a–e). Panel (b) presents the solution of the effective equation (3.16) (represented by the red solid curve), the solution of diffusion equation with diffusivity $1+\cot ^{2} (\theta )$ (illustrated by the blue dashed curve) and the solution to the diffusion equation with diffusivity 1 (depicted by the black dotted curve) at the time instant $t=1$. Panel (c) shows the corresponding curves at the time instant $t=5$.

Equation (4.3) suggests that as the inclination angle approaches zero, the dispersion rate may increase, possibly reaching infinity. It is crucial to emphasize that this conclusion holds true only when other parameters are kept constant. It is worth noting that as the inclination angle decreases, maintaining the same value of $\gamma$ necessitates keeping the same density variation in the $y_{1}$ direction, which in turn requires a larger density variation in the $x_{3}$ direction. Achieving such a condition in practical applications can pose significant challenges.

Next, we consider the case with small $\gamma$. In this limit, the effective equation for the stratified scalar (3.16) reduces to

(4.4)\begin{equation} \partial_{t}\bar{c}=\partial_{y_{1}} \left( \left( 1+\left( \frac{ {Ri}\,{Pe} \,{Re} \cos (\theta )}{72 \sqrt{70} } \partial_{y_{1}} \rho(\bar{c}) \right)^{2} \right) \partial_{y_{1}}\bar{c} \right). \end{equation}

The closed form expression of the solution of this equation is not available. However, we can obtain the solution in some limits for unbounded domain. When the nonlinear term is much larger than the diffusion term, namely, $({{Ri}\,{Pe}\,{Re} \cos (\theta )}/{72 \sqrt {70} }) \partial _{y_{1}}\rho (\bar {c}) \gg 1$, and $\partial _{c}\rho$ is a constant, (4.4) can be approximated by

(4.5)\begin{equation} \partial_{t}\bar{c}=\partial_{y_{1}} ( \alpha^{2} ( \partial_{y_{1}}\bar{c} )^{3} ), \quad \alpha=\frac{{Ri}\,{Pe}\,{Re} \cos (\theta ) \partial_{c}\rho}{72 \sqrt{70}}. \end{equation}

We can look for the self-similarity solution, $c (y_{1},t)=f ( \xi )$, where $\xi =t y_{1}^{-{1}/{4}}$. Equation (4.5) becomes

(4.6)\begin{equation} f'(\xi ) (12 \alpha^{2} f'(\xi ) f''(\xi )+\xi )=0. \end{equation}

The solution is

(4.7)\begin{align} f (\xi)= \begin{cases} C_{2}+ {\rm \pi}C_{1} \alpha \sqrt{3} & \xi \geq 2 \alpha \sqrt{6 C_{1}},\\ C_{2}- \dfrac{\dfrac{1}{2} \xi \sqrt{24 \alpha^{2} c_1-\xi ^2}+12 \alpha^{2} c_1 \tan ^{{-}1}\left(\dfrac{\xi }{\sqrt{24 \alpha^{2} c_1-\xi ^2}}\right)}{2 \sqrt{3} \sqrt{\alpha }} & - 2 \alpha \sqrt{6 C_{1}} \leq \xi \leq 2 \alpha \sqrt{6 C_{1}} ,\\ C_{2}- {\rm \pi}c_1 \alpha \sqrt{3} & \xi \leq{-}2 \alpha \sqrt{6 C_{1}}, \end{cases} \end{align}

where $C_{1}$, $C_{2}$ are the constants that can be defined by the boundary condition at infinity. For the case $\bar {c} (\infty )=0$ and $\bar {c} (-\infty )=1$, we have $C_{1}={1}/{2 \sqrt {3} {\rm \pi}\alpha }$, $C_{2}=\frac {1}{2}$. Notice that the similarity variable for the pure diffusion equation is $\xi ={y_{1}}/{\sqrt {t}}$. Therefore, $t\rightarrow \infty$, the neglected diffusion term in the above calculation becomes dominant eventually.

Last, we can represent $\gamma$ in terms of the physical parameters as follows: $\gamma = ({1}/{\sqrt {2}})(({g H^3 \rho _{0}}/{\kappa \mu })\sin \theta \partial _{y_{1}}\rho (\bar {c}))^{{1}/{4}}$. Therefore, smaller diffusivity, lower viscosity, greater gravitational constant, larger gap thickness and higher density result in a larger value of $\gamma$. We can calculate $\gamma$ using the practical parameters provided in (2.8a–d). If we set $\theta ={{\rm \pi} }/{4}$ and estimate the density gradient as $\partial _{y_{1}}\rho (\bar {c}) =\sin \theta \partial _{x_{3}} \rho (\bar {c})$, we obtain $\gamma \approx 8.94506$ and $\kappa _{eff} (\gamma ) \approx 1.720605$. With the same parameters and a smaller inclination angle $\theta ={{\rm \pi} }/{10}$, we have $\gamma \approx 5.91333$ and a larger effective diffusivity $\kappa _{eff} (\gamma ) \approx 6.46129$. In these cases, we cannot employ the approximate equations (4.3) and (4.4). Instead, we must use the complete effective equation (3.16) to solve for the stratified scalar.

4.1. Long time behaviour of the solution

As demonstrated in our previous examples, diffusion-driven flow enhances the mixing of a stratified scalar at some time scales. One might expect the difference between the density profile with and without the fluid to increase as time progresses or at least persist at long times. However, intriguingly, this difference actually vanishes over extended periods.

To understand this, it is essential to recall an interesting observation that emerges when considering different solutions to the same diffusion equation. Under some conditions, these solutions converge to a self-similarity solution asymptotically at long times, as discussed in Newman (1984). To illustrate this point, consider two solutions: $f_{1}=\frac {1}{2} \textrm {erfc}({x}/{2 \sqrt {t}})$ and $f_{2}=\frac {1}{2} \textrm {erfc}({x}/{2 \sqrt {\sigma ^2+t}})$, both satisfying the diffusion equation $\partial _{t}f=\partial _{x}^{2}f$. Although they begin with different initial conditions, the relative difference between them diminishes as $t\rightarrow \infty$:

(4.8)\begin{equation} \frac{f_{1}-f_{2}}{f_{1}}= \frac{\sigma ^2 x t^{-{3}/{2}}}{2 \sqrt{{\rm \pi} }}+ O ( t^{{-}2} ). \end{equation}

The time scale for this convergence depends on the difference in the initial conditions and can be multiple times the diffusion time scale.

This observation implies that even if diffusion-driven flow initially amplifies the dispersion of the stratified agent, when the density gradient is sufficiently weak, the governing equation for the stratified scalar approximates a diffusion equation with a diffusivity of unity. While some disparity remains between the solution of the effective equation and the solution of the pure diffusion equation, this difference diminishes over longer time scales.

To confirm this conclusion, we solve the effective equation (3.16) with parameters provided in (3.24a–e). In figure 11(a), we plot $\max _{y_{1}}\gamma$ and $\kappa _{eff} (\max _{y_{1}}\gamma )-1$ as functions of time, providing an estimation of the enhancement of dispersion induced by diffusion-driven flow. We observe that the enhanced diffusivity decreases below 0.1 after $t>20$, which implies the effective equation is close to the pure diffusion equation after that time. In figure 11(b), we display the relative difference between the density profiles with and without flow. This difference initially increases, reaching its maximum around $t=10$, but subsequently decreases to 0.003 at $t=50$.

Figure 11. (a) The blue dashed curve represents $\max _{y_{1}}\gamma$ and is associated with the left-hand $y$ axis. The red solid curve represents $\kappa _{eff} (\max _{y_{1}}\gamma )-1$ and is associated with the right-hand $y$ axis. (b) The relative difference between the density profile with and without the fluid as a function of time. The parameters used in the simulation are provided in (3.24a–e).

This result demonstrates that in a confined domain with no-flux boundary condition, the density profile in a system with the diffusion driven flow asymptotically converges to the one without the flow at extremely large time scale.

5. Comparison with the thin film approximation

As previously mentioned, we can employ the classical asymptotic expansion to analyse the governing equation (2.6). In this section, we will utilize the thin film approximation and subsequently compare the results with those obtained previously.

In the thin film approximation, we select distinct length scales in different directions for the purpose of non-dimensionalization:

(5.1)\begin{align} \left.\begin{array}{c} \displaystyle L y_{1}'=y_{1}, \quad H y_{3}'=y_{3}, \quad V_{i} v_{i}'=v_{i}, \quad \dfrac{\rho_{0}L\nu U}{H^{2}}p'=p, \quad \dfrac{L^{2}}{\kappa} t'=t, \quad \rho_{0}\rho'=\rho, \\ \displaystyle c_{0}c'=c,\quad \epsilon=\dfrac{V_{3}}{V_{1}}=\dfrac{H}{L} \ll 1, \quad {Re}=\dfrac{ V_{1} H}{\nu},\quad {Pe}= \dfrac{V_{1}H}{\kappa}, \quad {Sc}= \dfrac{\mu}{\rho_{0} \kappa},\quad {Ri}= \dfrac{ H g}{V_{1}^{2}}. \end{array} \right\} \end{align}

The non-dimensionalized governing equation becomes

(5.2)\begin{align} \left.\begin{array}{c} \displaystyle \rho' \left( \dfrac{\epsilon^{2}}{Sc}\partial_{t'} v_{1}'+\epsilon {Re} ( v_{1}'\partial_{y_{1}'} v_{1}' + v_{3}'\partial_{y_{3}'} v_{1}') \right) = ( \epsilon^{2}\partial_{y_{1}'}^{2} v_{1}' + \partial_{y_{3}'}^{2} v_{1}' )- \partial_{y_{1}'} p'- {Re}\, {Ri}\rho' \sin \theta, \\ \displaystyle \rho'\left( \dfrac{\epsilon^{3}}{Sc}\partial_{t'} v_{3}'+ \epsilon^{2} {Re}(v_{1}'\partial_{y_{1}'} v_{3}' + v_{3}'\partial_{y_{3}'} v_{3}') \right)= ( \epsilon^{3} \partial_{y_{1}'}^{2} v_{3}' + \epsilon \partial_{y_{3}'}^{2} v_{3}' ) - \dfrac{1}{\epsilon}\partial_{y_{3}'} p'- {Re}\, {Ri} \rho' \cos \theta, \\ \epsilon^{2}\partial_{t'}c'+ \epsilon {Pe} v_{1}'\partial_{y_{1}'}c'+ \epsilon{Pe} v_{3}'\partial_{y_{3}'}c'= \epsilon^{2} \partial_{y_{1}'}^{2}c'+ \partial_{y_{3}'}^{2}c',\\ \partial_{y_{1}'}v_{1}'+\partial_{y_{3}'}v_{3}'=0. \end{array} \right\} \end{align}

Dropping primes and rearranging the equation results in

(5.3)\begin{align} \left.\begin{array}{c} \displaystyle \rho\left( \dfrac{\epsilon^{2}}{Sc}\partial_{t} v_{1}+\epsilon {Re} ( v_{1}\partial_{y_{1}} v_{1} + v_{3}\partial_{y_{3}} v_{1}) \right)= ( \epsilon^{2}\partial_{y_{1}}^{2} v_{1} + \partial_{y_{3}}^{2} v_{1} )- \partial_{y_{1}} p- {Re}\, {Ri} \rho \sin \theta,\\ \displaystyle \rho \left( \dfrac{\epsilon^{4}}{Sc}\partial_{t} v_{3}+ \epsilon^{3} {Re}( v_{1}\partial_{y_{1}} v_{3} + v_{3}\partial_{y_{3}} v_{3}) \right)= \epsilon^{2}( \epsilon^{2} \partial_{y_{1}}^{2} v_{3} + \partial_{y_{3}}^{2} v_{3} ) - \partial_{y_{3}} p- \epsilon {Re}\, {Ri} \rho \cos \theta, \\ \epsilon^{2}\partial_{t}c+\epsilon {Pe} ( v_{1}\partial_{y_{1}}c+ v_{3}\partial_{y_{3}}c )= \epsilon^{2} \partial_{y_{1}}^{2}c+ \partial_{y_{3}}^{2}c,\\ \partial_{y_{1}}v_{1}+\partial_{y_{3}}v_{3}=0. \end{array} \right\} \end{align}

We consider the formal power series expansions for the velocity, stratified scalar and pressure:

(5.4a––d)\begin{equation} v_{i}=\sum_{k=0}^{\infty}v_{i,k}\epsilon^{k}, \quad c=\sum_{k=0}^{\infty}c_{k}\epsilon^{k},\quad p=\sum_{k=0}^{\infty}p_{k}\epsilon^{k}, \quad \rho (c)= \sum_{k=0}^{\infty} \rho_{k}\epsilon^{k}. \end{equation}

We have $v_{i,k} |_{\boldsymbol {y}\in \partial \varOmega }=0$ from the no-slip boundary condition of the velocity field, and $\partial _{\boldsymbol {n}} c_{k}|_{\boldsymbol {y} \in \partial \varOmega }=0$ from the no-flux boundary condition of the stratified scalar. The coefficient in the expansion of $\rho$ can be calculated from $c_k$ and the Taylor expansion of $\rho (c)$.

Substituting the expansion into the governing equation and collecting $O(1)$ terms yield the following equation:

(5.5)\begin{equation} 0= \partial_{y_{3}}^{2} v_{1,0} - \partial_{y_{1}} p_{0}- {Re}\,{Ri} \rho_{0} \sin \theta, \quad 0=\partial_{y_{3}}p_{0},\quad 0 = \partial_{y_{3}}^{2}c_{0},\quad \partial_{y_{1}}v_{1,0}+\partial_{y_{3}}v_{3,0}=0. \end{equation}

There are two possible group of solutions,

(5.6a,b)\begin{equation} v_{1,0}=v_{3,0}=0, \quad p_{0}= {Re}\,{Ri} \sin \theta \int_0^{y_{1}}\rho_{0} \,\mathrm{d} y_{1}, \end{equation}

and

(5.7a––c)\begin{equation} v_{1,0}=C y_{3} (y_{3}-1), \quad v_{3,0}=0, \quad p_{0}={-}2Cy_{1}+ {Re}\, {Ri} \sin \theta \int_0^{y_{1}}\rho_{0} \, \mathrm{d} y_{1}, \end{equation}

where $C$ is a constant number and $c_{0}, \rho _0$ are independent of $y_{3}$. The second solution describes a flow generated by a constant pressure gradient, which does not represent the correct physical problem. Therefore, we regard the first solution as the correct leading-order approximation.

Collecting the terms that is comparable to $O(\epsilon )$ yields the following equation:

(5.8)\begin{equation} \left.\begin{array}{c} {Re}\rho_{0} ( v_{1,0}\partial_{y_{1}} v_{1,0}+ v_{3,0}\partial_{y_{3}} v_{1,0}) = \partial_{y_{3}}^{2} v_{1,1} - \partial_{y_{1}} p_{1}- {Re}\,{Ri} \rho_{1} \sin \theta, \\ 0={-}\partial_{y_{3}} p_{1}-{Re}\,{Ri} \rho_{0} \cos \theta, \quad {Pe} ( v_{1,0}\partial_{y_{1}}c_{0} + v_{3,0}\partial_{y_{3}}c_{0} )= \partial_{y_{3}}^{2}c_{1},\\ \partial_{y_{1}}v_{1,1}+\partial_{y_{3}}v_{3,1}=0. \end{array} \right\} \end{equation}

Using the conclusions (5.6a,b) reduces the above equation to

(5.9)\begin{equation} \left.\begin{array}{c} 0 = \partial_{y_{3}}^{2} v_{1,1} - \partial_{y_{1}} p_{1}- {Re}\,{Ri} \rho_{1} \sin \theta, \quad 0={-}\partial_{y_{3}} p_{1}-{Re}\,{Ri} \rho_{0} \cos \theta, \\ 0= \partial_{y_{3}}^{2}c_{1},\quad \partial_{y_{1}}v_{1,1}+\partial_{y_{3}}v_{3,1}=0. \end{array} \right\} \end{equation}

Similar to the equation for $O(1)$ terms, the solution is not unique. To eliminate the solution representing the pressure driven flow, we impose the condition $\bar {p}_{1}=0$. From the second equation, we have

(5.10)\begin{equation} p_{1}={-}(\kern0.08em y_{3}-\tfrac{1}{2}){Re}\, {Ri} \rho_{0} \cos \theta. \end{equation}

Since $c_{1}$ is independent of $y_{3}$, it follows that $\rho _{1}=c_{1}\partial _{c}\rho (c_{0})$ is also independent of $y_{3}$. Averaging the first equation yields $\bar {v}_{1,1}=0$. Differentiating the first equation in (5.9) with respect to $y_{3}$ twice give us

(5.11a––c)\begin{equation} 0 = \partial_{y_{3}}^{4} v_{1,1}, \quad v_{1,1} |_{y=0,1 }=0,\quad \partial_{y_{3}}^{3}v_{1,1} |_{y=0,1 }={-} {Re}\, {Ri} \cos \theta \partial_{y_{1}} \rho_{0}. \end{equation}

The solution is

(5.12)\begin{equation} v_{1,1}={-} {Re}\,{Ri} \cos \theta \partial_{y_{1}} \rho_{0} \frac{y_{3} (y_{3}-1) (2 y_{3}-1)}{12}, \end{equation}

which is consistent with the first term in (3.12). It is important to note that the definition of $v_{1,1}$ in (3.12) differs from the one presented here. To make a comparison, $v_{1,1}$ in (3.12) need to be multiplied with $\partial _{y_{1}}\bar {c}$.

From the third equation in (5.9), we conclude that $c_{1}$ is independent of $y_{3}$. Moreover, after substituting the expression of $v_{1,1}$ back to the first equation in (5.9), we obtain $\rho _{1}=0$. The continuity equation gives us the expression of $v_{3,1}$, as follows:

(5.13)\begin{equation} v_{3,1}={-}\int_0^{y_{3}} \partial_{y_{1}}v_{1,1}\,\mathrm{d} y_{3}= {Re}\, {Ri} \cos \theta \partial_{y_{1}}^{2} \rho_{0} \frac{y_{3} ^{2} (y_{3}-1)^{2}}{24}, \end{equation}

which is consistent with the first term in (3.19).

Collecting the terms that are comparable to $O(\epsilon ^{2})$ yields the following equation:

(5.14)\begin{align} \left.\begin{array}{c} \begin{aligned} & \dfrac{\rho_{0}}{Sc}\partial_{t} v_{1,0}+{Re} \rho_{0}( v_{1,0}\partial_{y_{1}} v_{1,1}+ v_{1,1}\partial_{y_{1}} v_{1,0} + v_{3,0}\partial_{y_{3}} v_{1,1}+ v_{3,1}\partial_{y_{3}} v_{1,0}) \\ & \quad + {Re} \rho_{1}( v_{1,0}\partial_{y_{1}} v_{1,0} + v_{3,0}\partial_{y_{3}} v_{1,0} )= (\partial_{y_{1}}^{2} v_{1,0} + \partial_{y_{3}}^{2} v_{1,2} )- \partial_{y_{1}} p_{2}- {Re}\, {Ri} \rho_{2} \sin \theta, \end{aligned}\\ 0 = \partial_{y_{3}}^{2} v_{3,0} - \partial_{y_{3}} p_{2}- {Re}\, {Ri} \rho_{1} \cos \theta, \\ \partial_{t}c_{0}+ {Pe} ( v_{1,0}\partial_{y_{1}}c_{1}+ v_{1,1}\partial_{y_{1}}c_{0}+ v_{3,0}\partial_{y_{3}}c_{1} + v_{3,1}\partial_{y_{3}}c_{0} )= \partial_{y_{1}}^{2}c_{0}+ \partial_{y_{3}}^{2}c_{2},\\ \partial_{y_{1}}v_{1,2}+\partial_{y_{3}}v_{3,2}=0. \end{array} \right\} \end{align}

Using the current conclusions for $O (\epsilon )$ terms, the above equation reduces to

(5.15)\begin{equation} \left.\begin{array}{c} 0 = \partial_{y_{3}}^{2} v_{1,2} -\partial_{y_{1}} p_{2}- {Re}\, {Ri} \rho_{2} \sin \theta, \quad 0 ={-} \partial_{y_{3}} p_{2}- {Re}\, {Ri} \rho_{1} \cos \theta, \\ \partial_{t}c_{0}+ {Pe} v_{1,1}\partial_{y_{1}}c_{0}= \partial_{y_{1}}^{2}c_{0}+ \partial_{y_{3}}^{2}c_{2}. \end{array} \right\} \end{equation}

Notice that $\int _0^{1}v_{1,1} \,\mathrm {d} y_{3}=0$ and $c_0$ is independent of $y_{3}$. Integrating on both sides of the last equation in (5.15) with respect to $y_{3}$ and using the no-flux boundary conditions yields the evolution equation for $c_{0}$: $\partial _{t}c_{0}= \partial _{y_{1}}^{2}c_{0}$.

Now the last equation in (5.15) becomes ${Pe} v_{1,1}\partial _{y_{1}}c_{0}= \partial _{y_{3}}^{2}c_{2}$ and implies

(5.16)\begin{equation} c_{2}=\frac{{Pe}\,{Re}\,{Ri} (1-2 y_{3}^3(10+3 y_{3} (2 y_{3}-5))) \cos (\theta ) \partial_{y_{1}}c_{0} \partial_{y_{1}} \rho_{0} }{1440}, \end{equation}

which is consistent with the second term in (3.13). It is important to note that the definition of $c_{1}$ in § 3 differs from the one presented here.

The stratified scalar approximation is expressed as $c=c_{0}+\epsilon ^{2}c_{2}+O (\epsilon ^{3})$, which serves as the solution to a pure diffusion equation augmented by a higher-order correction term. The fact that $\bar {c}_{2}=0$ implies that as $\epsilon$ tends to zero, the diffusion-driven flow only distorts the contour line of the stratified scalar without amplifying the dispersion of the stratified scalar in the longitudinal direction of the channel.

So far, we have observed that the results obtained from the thin film approximation are consistent with the results obtained in § 3 when the density gradient vanishes, i.e. $\partial _{y_{1}}\bar {\rho }\rightarrow 0$, or equivalently when $\gamma \rightarrow 0$. Therefore, the thin film approximation presented here may not accurately model systems where the diffusion-driven flow significantly enhances dispersion, as is the case with the parameters provided in (3.25a–b), where the corresponding $\gamma =9$–15 and the flow make visible enhancement of the scalar dispersion as shown in figure 10. For those parameter combinations, an asymptotic analysis using a scaling relation that differs from the one presented in this section become necessary.

6. Conclusion and discussion

This paper explores the diffusion-driven flow in a tilted parallel-plate channel domain with a nonlinear density stratification. By employing a novel asymptotic expansion provided in (3.3), we derive leading-order approximations for the velocity field (3.11) and (3.19), stratified scalar (3.11) and pressure field (3.21). Furthermore, we formulate an effective equation (3.16) to describe the cross-sectional averaged stratified scalar, and its accuracy is confirmed through numerical simulations of the governing equation (3.2).

The effective equation reveals that the dynamics of the stratified scalar depend on the dimensionless parameter $\gamma$, defined as $\gamma = ({1}/{\sqrt {2}})(-{Re}\, {Pe}\, {Ri} \sin \theta \partial _{y_{1}}\rho (\bar {c}))^{{1}/{4}}= ({1}/{\sqrt {2}}) (({g L^3 \rho _{0}}/{\kappa \mu }) \sin \theta \partial _{y_{1}} \rho (\bar {c}))^{{1}/{4}}$. When $\gamma$ is large, the system behaves akin to a diffusion equation with an enhanced diffusion coefficient of $1+\cot ^{2}\theta$, where $\theta$ is the inclination angle relative to the horizontal plane. This result reveals an upper bound for the mixing capability of the diffusion-driven flow. Conversely, in the small $\gamma$ limit, the behaviour of the stratified scalar approximates a pure diffusion process, with the diffusion-driven flow failing to amplify the dispersion of the stratified scalar significantly. Additionally, we demonstrate that in a confined domain with no-flux boundary conditions, the density profile in a system featuring diffusion-driven flow asymptotically converges to the one without flow over long times, although this convergence occurs on a time scale much larger than the diffusion time scale.

Moreover, we establish that the thin film approximation aligns with the results obtained using the novel expansion when the diffusion-driven flow is weak. In such scenarios, the stratified scalar can be modelled by a diffusion equation featuring molecular diffusivity, and the diffusion-driven flow primarily distorts the scalar distribution without significantly increasing longitudinal dispersion. Consequently, the thin film approximation falls short in describing systems with relatively large density gradients, where diffusion-driven flow markedly enhances dispersion. Importantly, we both numerically and theoretically demonstrate that the proposed expansion effectively addresses these situations.

Future research directions encompass several avenues. First, the passive scalar transport in the channel with non-flat boundaries has various practical applications (Chang & Santiago 2023; Roggeveen, Stone & Kurzthaler 2023). Therefore, our aim is to develop a reduced model of diffusion-driven flow in channels with rough boundaries, such as rock fissures or microfluidic devices. When the amplitude and wavelength of the boundary variation are small, the rough boundary with a no-slip boundary condition can be approximated by a flat wall with an effective slip boundary condition obtained through the multiscale method (Achdou, Pironneau & Valentin 1998; Carney & Engquist 2022). In this scenario, the method presented in our work can be readily applied. Alternatively, when the wavelength of the boundary variation is large, the asymptotic expansion presented by Mercer & Roberts (1990) could offer a solution. Second, while this work primarily focuses on analysing parallel plate domains, it is worth noting that the method proposed herein is applicable to channels with arbitrary cross-sections. Concentrating on the parallel plates domain stems from the availability of an exact solution for the auxiliary problem. Investigating this fundamental geometry can augment our understanding of domains with more complex shapes, such as tilted cylindrical cavities embedded in rocks (Sánchez, Higuera & Medina 2005), or tilted square containers (Page 2011a,b; French 2017; Grayer et al. 2020).