Buoyancy-modulated Lagrangian drift in wavy-walled vertical channels as a model problem to understand drug dispersion in the spinal canal

Abstract

This paper investigates flow and transport in a slender wavy-walled vertical channel subject to a prescribed oscillatory pressure difference between its ends. When the ratio of the stroke length of the pulsatile flow to the channel wavelength is small, the resulting flow velocity is known to include a slow steady-streaming component resulting from the effect of the convective acceleration. Our study considers the additional effect of gravitational forces in configurations with a non-uniform density distribution. Specific attention is given to the slowly evolving buoyancy-modulated flow emerging after the deposition of a finite amount of solute whose density is different from that of the fluid contained in the channel, a relevant problem in connection with drug dispersion in intrathecal drug delivery (ITDD) processes, involving the injection of the drug into the cerebrospinal fluid that fills the spinal canal. It is shown that when the Richardson number is of order unity, the relevant limit in ITDD applications, the resulting buoyancy-induced velocities are comparable to those of steady streaming. As a consequence, the slow time-averaged Lagrangian motion of the fluid, involving the sum of the Stokes drift and the time-averaged Eulerian velocity, is intimately coupled with the transport of the solute, resulting in a slowly evolving problem that can be treated with two-time-scale methods. The asymptotic development leads to a time-averaged, nonlinear integro-differential transport equation that describes the slow dispersion of the solute, thereby circumventing the need to describe the small concentration fluctuations associated with the fast oscillatory motion. The ideas presented here can find application in developing reduced models for future quantitative analyses of drug dispersion in the spinal canal.

1. Introduction

The steady Lagrangian drift generated in oscillatory viscous flows in pipes and channels is known to play an important role in different heat and mass transport processes, including those occurring in extracorporeal membrane oxygenators (Bellhouse et al. 1973), pulmonary high-frequency ventilation devices (Grotberg 1994), compact heat exchangers (Mackley & Stonestreet 1995) and drug dispersion in the spinal canal (Lawrence et al. 2019). For configurations with slowly varying cross-section, the lubrication approximation can be used to derive insightful analytical results, with seminal analyses including those of Hall (1974), who considered flow in a pipe subject to a harmonic pressure difference, and Grotberg (1984), who investigated flow in a tapered channel subject to a prescribed oscillating stroke volume. More recent analytical studies pertaining to channels include those of Lo Jacono, Plouraboué & Bergeon (2005) and Guibert, Plouraboué & Bergeon (2010), involving three-dimensional wavy-walled configurations, and that of Larrieu, Hinch & Charru (2009), who considered an oscillating Couette flow over a wavy bottom. All of these analytical investigations of oscillating slender flows addressed configurations with weak inertia, corresponding to small values of the ratio $\varepsilon$ of the stroke length to the characteristic longitudinal length, with $\varepsilon ^{-1} \gg 1$ representing the relevant Strouhal number. The asymptotic analysis for $\varepsilon \ll 1$ reveals that the velocity, resulting from a balance between the local acceleration and the pressure and viscous forces, is harmonic at leading order, with the small corrections arising from the convective acceleration providing a small steady-streaming component of order $\varepsilon$ (Riley 2001). This steady streaming determines, together with the Stokes drift associated with the leading-order harmonic flow, the Lagrangian drift experienced by the fluid particles, with both contributions having in general comparable orders of magnitude (Larrieu et al. 2009).

As discussed by Guibert et al. (2010), the fundamental pulsatile-flow investigations mentioned previously are relevant in connection with the motion of cerebrospinal fluid (CSF) along the spinal subarachnoid space, a slender annular canal surrounding the spinal cord, depicted in figure 1. The flow features an oscillatory velocity driven by the pressure pulsations induced by the cardiac and respiratory cycles (Linninger et al. 2016). The dynamics of this flow and its associated Lagrangian transport are fundamental in understanding the role of CSF as a vehicle for metabolic-waste clearance (Klarica, Radoš & Orešković 2019) and also to quantify drug dispersion in intrathecal drug delivery (ITDD) (Onofrio 1981; Lynch 2014; Fowler et al. 2020), a medical procedure used for treatment of some cancers (Lee et al. 2017), infections (Remeš et al. 2013) and pain (Bottros & Christo 2014). The standard ITDD protocol involves the placement of a small catheter along the lumbar section of the spinal canal to continuously pump the drug or to release a finite dose at selected times. The transport of the drug depends fundamentally on its physical properties, including molecular diffusivities $\kappa$ that are much smaller than the kinematic viscosity $\nu$ of CSF.

Figure 1. A schematic view of the spinal canal, showing the location of (a) ITDD and (b,c) of the two-dimensional channel flow investigated here.

The flow in the spinal canal has been investigated computationally in numerous studies addressing different aspects of the problem, as summarised in a recent paper by Khani et al. (2018), although corresponding theoretical analyses are more scarce. Exact solutions for pulsatile viscous flow in a straight elliptic annulus have been proposed as a representation for the oscillatory flow in the spinal canal (Gupta, Poulikakos & Kurtcuoglu 2008). More recent studies, modelling the canal as a linearly elastic annular pipe of slowly varying section, have employed the lubrication approximation in the asymptotic limit $\varepsilon \ll 1$ to quantify steady streaming (Sánchez et al. 2018) and to investigate the dispersion of a solute (Lawrence et al. 2019). For the large values of the Schmidt number $\nu /\kappa \sim 1000$ corresponding to the molecular diffusivities of all ITDD drugs, the analysis of Lawrence et al. (2019) showed that the Lagrangian mean flow is the key mechanism responsible for the dispersion of the solute, whereas shear-enhanced Taylor dispersion has a negligibly small effect. An important outcome of the asymptotic analysis is a time-averaged transport equation that has been recently validated by means of comparisons with results of direct numerical simulations (DNS) spanning hundreds of oscillation cycles (Gutiérrez-Montes et al. 2021), as needed to generate significant dispersion of the solute. The comparisons demonstrate the accuracy of the reduced description, which is seen to provide excellent fidelity at a fraction of the computational cost involved in the DNS.

Our previous analysis of solute dispersion (Lawrence et al. 2019) assumed the density of the solute $\rho _s$ and the density of the carrier fluid $\rho$ to be identical, thereby neglecting the small density differences $\rho -\rho _s$ found in ITDD applications, for which the values of $\rho -\rho _s$ typically range from positive values of order $(\rho -\rho _s) \sim \rho /1000$ for drugs diluted with water to negative values of order $(\rho -\rho _s)/ \sim -\rho /100$ for drugs diluted with dextrose (Nicol & Holdcroft 1992; Lui, Polis & Cicutti 1998). These relative density differences $(\rho -\rho _s)/\rho \ll 1$ between the drug and the CSF, although extremely small, are known by clinicians to play an important role in the dispersion rate of ITDD drugs for patients in a sitting position, when buoyancy forces are nearly aligned with the canal. It has been seen that for a hyperbaric (dense) solution, injection while the patient is seated for some time before moving to a supine position leads to an initial restriction in the transport of the anaesthesia (Mitchell et al. 1988; Povey, Jacobsen & Westergaard-Nielsen 1989; Veering et al. 2001). (In spinal anaesthesia, baricity is the term used to refer to the density of the anaesthetic relative to that of the CSF. Thus, an anaesthetic is said to be hyperbaric/hypobaric when its density is higher/lower than that of the CSF, whereas the term isobaric describes anaesthetics whose density matches exactly that of the CSF.) Conversely, for a hypobaric (light) solution, the sitting injection position leads to more rapid cephalad spread of the anaesthesia as compared to a lateral injection position (Richardson et al. 1996). As could be expected, the density of the drug is inconsequential when injection occurs in the lateral position (Hallworth et al. 2005) and, similarly, positioning has no effect on the spread rate when the solution density matches that of CSF (Wildsmith et al. 1981). Given the abundance of clinical evidence on the importance of buoyancy forces on the drug dispersion rate, there is interest in developing a quantitative description; the present paper, focused on a simplified geometry, is a necessary first step in that direction.

In looking for a simplified geometrical model, we follow Guibert et al. (2010) in noting that the width $h_o \sim 1\unicode{x2013}2$ mm of the annular spinal canal is smaller than the spinal-cord diameter ${\sim}1$ cm, with the consequence that a two-dimensional channel can be used to describe many aspects of the flow. The channel is placed in a vertical position, as is appropriate in describing buoyancy effects for a patient in a sitting position. As indicated in figure 1, the quasi-periodic variation of the canal section, associated with the presence of the vertebrae, will be modelled by including a wavy boundary whose wavelength $\lambda$ mimics the inter-vertebral distance. The channel will be assumed to be slender in that $\lambda \gg h_o$, a good approximation in the spinal canal, where $\lambda \sim 2\unicode{x2013}4$ cm and $h_o/\lambda \simeq 0.05$. For simplicity, the total channel length is taken to be an integer multiple of the wavelength, so that the channel contains a finite number of identical cells. As in the seminal analysis of Hall (1974), an oscillating pressure difference with angular frequency $\omega$ will be imposed between the channel ends, resulting in a pulsating flow. We investigate the buoyancy-modulated dispersion of a bolus of solute released inside the channel when the buoyancy-induced acceleration is comparable to the convective acceleration of the pressure-driven flow, those being the conditions of interest in ITDD applications, as explained later below (2.6).

The rest of the paper is organised as follows. The problem is formulated in dimensionless form in § 2, which includes the identification of the relevant non-dimensional parameters and a discussion of the essential features of the subsequent asymptotic analysis, including the existence of a long time scale $\varepsilon ^{-2} \omega ^{-1}$ for solute dispersion, additional to the much smaller oscillation time $\omega ^{-1}$. The asymptotic description of the velocity field is presented next in § 3, with the time-averaged Eulerian velocity including the familiar steady-streaming contribution stemming from the convective acceleration along with an additional buoyancy-induced component that depends on the distribution of solute. This velocity field is used in § 4 to analyse solute dispersion with use of a two-time scale asymptotic analysis, resulting in a time-averaged transport equation that describes the evolution of the flow in the long-time scale $\varepsilon ^{-2} \omega ^{-1}$. The reduced description stemming from the asymptotic analysis is validated in § 5 through comparisons with DNS. In addition, the model is used to quantify effects of buoyancy-induced motion on the solute dispersion for different values of the controlling parameters. Finally, concluding remarks are given in § 6.

2. Problem formulation

2.1. Governing equations

Consider a vertical wavy channel of average gap size $h_o$ filled with a Newtonian fluid of density $\rho$ and kinematic viscosity $\nu$ (for CSF, $\rho \simeq 10^{3}$ kg m$^{-3}$ and $\nu \simeq 0.7 \times 10^{-6}$ m$^{2}$ s$^{-1}$). The channel, open at both ends, is bounded by a flat surface and a wavy wall of wave length $\lambda \gg h_o$, so that the resulting channel width is $h=h_o[1+\beta \cos (2 {\rm \pi}x^{*}/\lambda )]$, where $x^{*}$ is the longitudinal distance measured from the upper end and $\beta < 1$ is the relative amplitude of the wall undulation. The total channel length is $n \lambda$, with $n$ representing a general integer number, so that the channel comprises $n$ identical cells. The flow is described using cartesian coordinates $(x^{*},y^{*})$, with $y^{*}$ measured from the flat surface, and corresponding velocity components $\boldsymbol{v}^{*}=(u^{*},v^{*})$. The Navier–Stokes equations describing the planar unsteady flow are written in the Boussinesq approximation

(2.1)\begin{gather} \boldsymbol{\nabla}^{*} \boldsymbol{\cdot} \boldsymbol{v}^{*} =0, \end{gather}

(2.2)\begin{gather}\frac{\partial \boldsymbol{v}^{*}}{\partial t^{*}}+\boldsymbol{v}^{*} \boldsymbol{\cdot} \boldsymbol{\nabla}^{*} \boldsymbol{v}^{*} ={-}\frac{1}{\rho} \boldsymbol{\nabla}^{*} p^{*}+ \nu \boldsymbol{\nabla}^{*2} \boldsymbol{v}^{*} - \frac{\rho-\rho_s}{\rho} g c \boldsymbol{e}_x, \end{gather}

(2.3)\begin{gather}\frac{\partial c}{\partial t^{*}}+\boldsymbol{v}^{*} \boldsymbol{\cdot} \boldsymbol{\nabla}^{*} c=\kappa \boldsymbol{\nabla}^{*2} c, \end{gather}

where $p^{*}$ is the sum of the pressure difference from the upper end and the constant-density hydrostatic component $-\rho g x^{*}$, $c$ is the solute volume concentration, $\kappa$ is the solute diffusivity, $\boldsymbol {\nabla }^{*}=(\partial /\partial x^{*},\partial /\partial y^{*})$ and $\boldsymbol{e}_x$ is the unit vector aligned with the gravitational acceleration.

A pressure difference $n \Delta p \cos (\omega t^{*})$ oscillating harmonically in time is prescribed between the upper and lower ends of the canal, driving a periodic fluid motion with angular frequency $\omega$. The resulting slender flow is characterised by longitudinal velocities of order $u_c=\Delta p/(\rho \omega \lambda )$, as follows from a balance between the local acceleration $\partial u^{*}/\partial t^{*} \sim u_c \omega$ and the pressure gradient $\rho ^{-1} \partial p^{*}/\partial x^{*} \sim \Delta p/(\lambda \rho )$, and much smaller transverse velocities of order $v_c=(h_o/\lambda )u_c \ll u_c$, as follows from the continuity balance $\partial u^{*}/\partial x^{*} \sim \partial v^{*}/\partial y^{*}$.

2.2. Controlling parameters

The analysis assumes that the viscous time across the channel $h_o^{2}/\nu$ is comparable to the characteristic oscillation time $\omega ^{-1}$, resulting in Womersley numbers

(2.4)\begin{equation} \alpha=\left(\frac{\omega h_o^{2}}{\nu}\right)^{1/2}, \end{equation}

of order unity. The limit $\alpha \sim 1$ is instrumental in analysing cardiac-driven CSF flow ($\omega =2{\rm \pi}$ s$^{-1}$) in the spinal canal, for which typical values of $\alpha$ are in the range $3 \lesssim \alpha \lesssim 6$, as can be seen by evaluating the above expression with $h_o \simeq 1\unicode{x2013}2$ mm and $\nu =0.7 \times 10^{-6}$ m$^{2}$ s$^{-1}$. In the lumbar region, the typical drug-delivery site in ITDD procedures, the average CSF speeds are of order $u_c \sim 1$ cm s$^{-1}$, so that the associated stroke lengths $u_c/\omega$ are much smaller than the characteristic longitudinal distance $\lambda \simeq 2\unicode{x2013}4$ cm. Their ratio

(2.5)\begin{equation} \varepsilon = \frac{u_c/\omega}{\lambda}, \end{equation}

of order $\varepsilon \simeq 0.05$ for spinal CSF flow, defines the small parameter employed in the following asymptotic description. As shown earlier (Hall 1974; Grotberg 1984), the solution at leading order is determined by a balance between the pressure gradient, the local acceleration and the viscous forces, with the convective acceleration introducing small corrections of relative magnitude $\varepsilon$. Although the leading-order motion is harmonic, the velocity corrections include a non-zero steady-streaming component.

The familiar periodic channel flow described previously is altered by gravitational forces when a solute of density $\rho _s \ne \rho$ is introduced in the channel. The extent of the resulting buoyancy-induced motion can be measured by the associated Richardson number

(2.6)\begin{equation} {{\textit{Ri}}}=\frac{\rho-\rho_s}{\rho} \frac{g \lambda}{u_c^{2}}, \end{equation}

which compares the order of magnitude of the effective gravitational acceleration $g (\rho -\rho _s)/\rho$ with that of the convective acceleration $\boldsymbol {v}^{*} \boldsymbol {\cdot} \boldsymbol {\nabla }^{*} \boldsymbol {v}^{*} \sim u_c^{2}/\lambda$. Our analysis addresses the limit ${{\textit {Ri}}} \sim 1$, which is relevant for drug dispersion in ITDD procedures, as can be seen by evaluating (2.6) with $\lambda \simeq 2$ cm and $u_c \sim 1$ cm s$^{-1}$ for density differences in the range $10^{-3} \lesssim |\rho -\rho _s|/\rho \lesssim 10^{-2}$.

Also motivated by ITDD applications, we consider solutes with diffusivities $\kappa$ much smaller than the kinematic viscosity, that always being the case of diffusion in liquid phase. As $\nu /\kappa \sim 1000$ and $\varepsilon \simeq 0.05$ in ITDD applications, the following analysis of solute dispersion will specifically address the distinguished limit $\kappa /\nu \sim \varepsilon ^{2}$, with solute diffusion correspondingly characterised by the reduced Schmidt number

(2.7)\begin{equation} \sigma=\varepsilon^{2} \frac{\nu}{\kappa}, \end{equation}

assumed to be of order unity.

2.3. Non-dimensional formulation

We address the motion that follows from the deposition of the solute inside an intermediate cell along the channel. The problem is non-dimensionalised using the scales identified previously to give the dimensionless variables

(2.8a–f)\begin{equation} t=\omega t^{*}, \quad x=\frac{x^{*}}{\lambda}, \quad y=\frac{y^{*}}{h_o}, \quad u=\frac{u^{*}}{u_c}, \quad v=\frac{v^{*}}{v_c}, \quad p=\frac{p^{*}}{\Delta p} \end{equation}

and associated conservation equations

(2.9)\begin{gather} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, \end{gather}

(2.10)\begin{gather}\frac{\partial u}{\partial t}+\varepsilon \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} u={-}\frac{\partial p}{\partial x}+\frac{1}{\alpha^{2}} \frac{\partial^{2} u}{\partial y^{2}}- \varepsilon \,{{\textit{Ri}}} \,c, \end{gather}

(2.11)\begin{gather}\frac{\partial p}{\partial y}=0, \end{gather}

(2.12)\begin{gather}\frac{\partial c}{\partial t}+\varepsilon \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} c = \frac{\varepsilon^{2}}{\alpha^{2} \sigma} \frac{\partial^{2} c}{\partial y^{2}}, \end{gather}

where $\boldsymbol {v}=(u,v)$ and $\boldsymbol {\nabla }=(\partial /\partial x,\partial /\partial y)$. In writing (2.9)–(2.12) from (2.1)–(2.3) we have used the slender-flow approximation resulting from the limit $h_o \ll \lambda$. Thus, the terms representing longitudinal diffusion of momentum and mass have been neglected in (2.10) and (2.12), because they are a factor $(h_o/\lambda )^{2}$ smaller than those associated with transverse diffusion. At the same level of approximation, the transverse component of the momentum equation takes the reduced form (2.11). The velocity and concentration must satisfy the boundary conditions

(2.13)\begin{equation} u=v=\frac{\partial c}{\partial y}=0 \quad {\rm at} \ \left\{ \begin{array}{@{}l} y=0 \\ y=H=1+\beta \cos(2 {\rm \pi}x) \end{array} \right., \end{equation}

corresponding to non-permeable no-slip surfaces, whereas the reduced pressure $p(x,t)$, independent of $y$, is identically zero at $x=0$ and takes the value $p=n \cos t$ at the lower end $x=n$.

In the absence of buoyancy (i.e. for ${{\textit {Ri}}}=0$), the solution for the velocity is periodic in time, including a steady component of order $\varepsilon$, and also periodic in space, so that the velocity distribution found in each cell is identical. On the other hand, for ${{\textit {Ri}}} \ne 0$ the motion is coupled to the solute transport, albeit weakly, with the result that the velocity necessarily evolves in time following the dispersion of the solute, which is driven partly by the steady streaming motion, with characteristic velocities $\varepsilon u_c$. It can be anticipated that the characteristic time for the slow evolution is that associated with the dispersion of the solute inside the deposition cell $\lambda /(\varepsilon u_c)=\varepsilon ^{-2} \omega ^{-1}$, much larger than the characteristic oscillation time $\omega ^{-1}$. These considerations suggest the introduction of a second time variable $\tau = \varepsilon ^{2} t$ for describing the slow evolution, additional to the fast time-scale variable $t$ describing the oscillatory motion. In the two-time-scale formalism, the time derivatives in (2.10) and (2.12) are replaced with $\partial /\partial t+\varepsilon ^{2} \partial /\partial \tau$ and the different variables are expressed in terms of the power expansions

(2.14)\begin{gather} u=u_0(x,y,t,\tau)+\varepsilon u_1(x,y,t,\tau)+\cdots, \end{gather}

(2.15)\begin{gather}v=v_0(x,y,t,\tau)+\varepsilon v_1(x,y,t,\tau)+\cdots, \end{gather}

(2.16)\begin{gather}p=p_0(x,t,\tau)+\varepsilon p_1(x,t,\tau)+\cdots, \end{gather}

(2.17)\begin{gather}c=c_0(x,y,t,\tau)+\varepsilon c_1(x,y,t,\tau)+\cdots, \end{gather}

with all functions assumed to be $2{\rm \pi}$ periodic in the fast time scale $t$. The asymptotic procedure leads to a hierarchy of problems that can be solved sequentially, as shown in the following.

3. Velocity description

We begin by describing the velocity field in the asymptotic limit $\varepsilon \ll 1$, following the procedure used in previous steady-streaming investigations of slender flows (Larrieu et al. 2009; Guibert et al. 2010; Sánchez et al. 2018). The solution at leading order and also the first-order corrections associated with convective acceleration are similar to those found earlier in three-dimensional wavy-walled channels (Guibert et al. 2010) and annular canals (Sánchez et al. 2018). These previous analyses did not address, however, effects of buoyancy forces, which are investigated here for order-unity values of the Richardson number ${{\textit {Ri}}}$, leading to a velocity correction that will be seen to be expressible in terms of integrals of the solute concentration.

3.1. Leading-order oscillatory flow

Convective acceleration and buoyancy are negligible at leading order, so that the velocity $\boldsymbol{v}_0=(u_0,v_0)$ satisfies a linear problem that can be solved in terms of the reduced variables

(3.1a–c)\begin{equation} u_0={\rm Re}(\mathrm{i} \mathrm{e}^{\mathrm{i} t} U),\quad v_0={\rm Re} (\mathrm{i} \mathrm{e}^{\mathrm{i} t} V),\quad p_0={\rm Re} (\mathrm{e}^{\mathrm{i} t} P), \end{equation}

where the complex functions $U(x,y)$, $V(x,y)$, and $P(x)$ satisfy

(3.2a,b)\begin{equation} \frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0 \quad {\rm and} \quad -U ={-}\frac{{\rm d} P}{{\rm d} x}+\frac{\mathrm{i}}{\alpha^{2}} \frac{\partial^{2} U}{\partial y^{2}}. \end{equation}

The second equation above can be integrated with boundary conditions $U=0$ at $y=0,H$ to give

(3.3)\begin{equation} U=\frac{{\rm d} P}{{\rm d} x} G(x,y), \end{equation}

where

(3.4a,b)\begin{equation} G=1-\frac{\cosh[\varLambda(2y/H-1)]}{\cosh \varLambda} \quad {\rm and} \quad \varLambda=\frac{\alpha}{2} \frac{1+\mathrm{i}}{\sqrt{2}} H(x). \end{equation}

The result can be used to integrate the first equation in (3.2a,b) subject to $V=0$ at $y=0$, yielding

(3.5)\begin{equation} V={-}\frac{\partial}{\partial x} \left(\int_0^{y} U \,{\rm d} \hat{y} \right)={-}\frac{\partial}{\partial x} \left(\frac{{\rm d} P}{{\rm d} x} \int_0^{y} G\, {\rm d}\hat{y} \right), \end{equation}

where

(3.6)\begin{equation} \int_0^{y} G \,{\rm d}\hat{y}=y-\frac{H}{2 \varLambda} \left\{\frac{\sinh[\varLambda(2y/H-1)]+\sinh \varLambda}{\cosh \varLambda} \right\}, \end{equation}

with $\hat {y}$ representing a dummy integration variable. Note that both velocity components $U$ and $V$ are spatially periodic in $x$ through the function $H=1+\beta \cos (2 {\rm \pi}x)$. The determination of the longitudinal pressure gradient ${\rm d}P/{\rm d} x$ that completes the solution begins by using the condition $V=0$ at $y=H$ in the first equation of (3.5) to give

(3.7)\begin{equation} \frac{{\rm d}}{{\rm d} x} \left(\int_0^{H} U \,{\rm d} y \right)=0, \end{equation}

indicating that the reduced flow rate

(3.8)\begin{equation} Q=\int_0^{H} U \,{\rm d} y=\frac{{\rm d} P}{{\rm d} x} \int_0^{H} G\, {\rm d} y \end{equation}

is constant. Further progress requires use of the conditions $P(0)=P(n)-n=0$, consistent with the boundary values $p(0,t)=0$ and $p(n,t)=n \cos t$ stated below (2.13). Using (3.6) to evaluate the integral $\int _0^{H} G \,{\rm d} y$ leads to the equation

(3.9)\begin{equation} Q=\frac{{\rm d} P}{{\rm d} x} H (1-\varLambda^{{-}1} \tanh \varLambda), \end{equation}

which can be integrated subject to $P(0)=0$ to yield the pressure distribution

(3.10)\begin{equation} P=Q \int_0^{x} \frac{{\rm d} \hat{x}}{H (1-\varLambda^{{-}1} \tanh \varLambda)}. \end{equation}

Using now the condition $P(n)=n$ provides

(3.11)\begin{equation} Q=n \left[\int_0^{n} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)} \right]^{{-}1}. \end{equation}

Owing to the spatial periodicity of $H$ and $\varLambda$, it follows that

(3.12)\begin{equation} \int_0^{n} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}=n \int_0^{1} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}, \end{equation}

thereby finally yielding

(3.13)\begin{equation} Q=\left[\int_0^{1} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}\right]^{{-}1} \end{equation}

and, from (3.9),

(3.14)\begin{equation} \frac{{\rm d} P}{{\rm d} x}=\left[H (1-\varLambda^{{-}1} \tanh \varLambda) \int_0^{1} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}\right]^{{-}1}, \end{equation}

independent of $n$. It is worth pointing out that, because at this order the velocity is harmonic, the associated time-averaged values $\langle u_0 \rangle$ and $\langle v_0 \rangle$ with $\langle {\cdot } \rangle = ({1}/{2 {\rm \pi}}) \int _t^{t+2 {\rm \pi}}\,{\rm d}t$ are identically zero, so that the steady bulk motion of the fluid occurs through the velocity corrections at the following order.

3.2. First-order corrections

Collecting terms of order $\varepsilon$ in (2.9) and (2.10) yields

(3.15)\begin{gather} \frac{\partial u_1}{\partial x}+\frac{\partial v_1}{\partial y}=0, \end{gather}

(3.16)\begin{gather}\frac{\partial u_1}{\partial t}+ \frac{\partial}{\partial x}(u_0^{2})+ \frac{\partial}{\partial y}(u_0 v_0)={-}\frac{\partial p_1}{\partial x}+\frac{1}{\alpha^{2}} \frac{\partial^{2} u_1}{\partial y^{2}}- {{\textit{Ri}}} \, c_0, \end{gather}

to be integrated with boundary conditions $u_1=v_1=0$ at $y=0,H$ and $p_1=0$ at $x=0,n$. There is interest in computing the corresponding time-averaged velocity correction $\langle \boldsymbol{v}_1 \rangle =(\langle u_1 \rangle,\langle v_1 \rangle )$. Taking the time average of (3.15) and (3.16) provides

(3.17a,b)\begin{equation} \frac{\partial \langle u_1 \rangle}{\partial x}+\frac{\partial \langle v_1 \rangle}{\partial y}=0 \quad {\rm and} \quad F(x,y)={-}\frac{\partial \langle p_1 \rangle}{\partial x}+\frac{1}{\alpha^{2}} \frac{\partial^{2} \langle u_1 \rangle}{\partial y^{2}}- {{\textit{Ri}}} \, c_0, \end{equation}

where the known function $F=\partial \langle u_0^{2} \rangle /\partial x+\partial \langle u_0 v_0 \rangle /\partial y$ can be expressed in terms of the complex velocities $U$ and $V$ defined above in the form

(3.18)\begin{equation} F=\frac{1}{2} {\rm Re}\left[\frac{\partial}{\partial x}(U \bar{U})+\frac{\partial}{\partial y}(V \bar{U}) \right], \end{equation}

a result following from the identity $\langle {\rm Re}(\mathrm {i} \mathrm {e}^{\mathrm {i} t} f_1) \, {\rm Re}(\mathrm {i} \mathrm {e}^{\mathrm {i} t} f_2) \rangle = {\rm Re}(f_1 \bar {f}_2)/2$, which applies to any generic time-independent complex functions $f_1$ and $f_2$, with the bar denoting complex conjugates. In writing (3.17a,b), we have anticipated that, at leading order, the solute concentration is independent of the fast time scale $t$, as follows from (2.12) when $\varepsilon \ll 1$, so that its time-averaged value $\langle c_0 \rangle$ reduces simply to $\langle c_0 \rangle =c_0$. Also of interest is that, because of the symmetry of $H(x)$, the periodic function $F$ defined in (3.18) is antisymmetric with respect to $x=1/2$, so that $F(x,y)=-F(1-x,y)$.

As can be concluded from (3.16), the velocity corrections arise partly owing to the convective acceleration and partly owing to the buoyancy force. In computing the corresponding time average, it is convenient to compute both contributions separately by introducing $\langle \boldsymbol{v}_1 \rangle =\boldsymbol{v}_{SS}+\boldsymbol{v}_{B}$ and $\langle p_1 \rangle =p_{SS}+p_{B}$, where $\boldsymbol{v}_{SS}(x,y)=(u_{SS},v_{SS})$ and $p_{SS}(x)$ describe the familiar steady-streaming associated with the nonlinear convective terms, which is independent of time and periodic in $x$, and $\boldsymbol{v}_{B}(x,y,\tau )=(u_{B},v_{B})$ and $p_{B}(x,\tau )$ describe the buoyancy-induced corrections, which evolve in the long time scale $\tau$ as the solute spreads in the channel.

3.3. Steady streaming

The solution procedure needed to compute the velocity corrections parallels that followed earlier at leading order. Thus, the longitudinal component of the steady-streaming velocity

(3.19)\begin{equation} \frac{u_{SS}}{\alpha^{2}}={-}\frac{{\rm d} p_{SS}}{{\rm d} x} \frac{1}{2} (H-y)y+y\int_0^{y} F \, {\rm d} \hat{y} -\int_0^{y} F \hat{y} \, {\rm d} \hat{y}-y \int_0^{H} F\left(1-\frac{y}{H} \right) {\rm d} {y} \end{equation}

is determined by integrating the momentum equation (3.16) written for ${{\textit {Ri}}}=0$ subject to the boundary conditions $u_{SS}=0$ at $y=0,H$. The result can be substituted into (3.15) to give

(3.20)\begin{align} \frac{v_{SS}}{\alpha^{2}}&=\frac{\partial}{\partial x} \left[\frac{{\rm d} p_{SS}}{{\rm d} x} \frac{y^{2}}{2} \left(\frac{H}{2}-\frac{y}{3}\right)+\frac{y^{2}}{2}\int_y^{H} F \left(1-\frac{\hat{y}}{H} \right) {\rm d} \hat{y} \right.\nonumber\\ &\quad \left.+\,y\left(1-\frac{y}{2H}\right)\int_0^{y} F \hat{y} \, {\rm d} \hat{y}-\frac{1}{2} \int_0^{y} F \hat{y}^{2} \, {\rm d} \hat{y}\right] \end{align}

upon integrating with $v_{SS}=0$ at $y=0$. To determine the unknown pressure gradient ${\rm d} p_{SS}/{\rm d} x$ we begin by using $v_{SS}=0$ at $y=H$ in the above equation to give

(3.21)\begin{equation} \frac{\rm d}{{\rm d} x} \left(\int_0^{H} \frac{u_{SS}}{\alpha^{2}} \,{\rm d} y\right)=\frac{\rm d}{{\rm d} x} \left(\frac{{\rm d} p_{SS}}\,{{\rm d} x} \frac{H^{3}}{12}+\frac{1}{2} \int_0^{H} F y (H-y) \, {\rm d} y \right)=0, \end{equation}

which indicates that the flow rate

(3.22)\begin{equation} Q_{SS}=\int_0^{H} u_{SS} \,{\rm d} y=\alpha^{2} \left[ \frac{{\rm d} p_{SS}}{{\rm d} x} \frac{H^{3}}{12}+\frac{1}{2} \int_0^{H} F y (H-y) \, {\rm d} y\right] \end{equation}

is constant. Its value can be determined by integrating a second time (3.22) subject to $p_{SS}(n)=p_{SS}(0)=0$ to yield

(3.23)\begin{equation} Q_{SS}= \frac{\displaystyle\alpha^{2} \int_0^{1} H^{{-}3} \left[\int_0^{H} F y (H-y) \,{\rm d} y \right] {\rm d} x}{2 \displaystyle\int_0^{1} H^{{-}3} \,{\rm d} x}, \end{equation}

when account is taken of the spatial periodicity of $H$ and $F$. Since $H(x)$ is symmetric about $x=1/2$ while $F(x,y)$ is antisymmetric, the double integral in the numerator of the above equation is identically zero, so that

(3.24)\begin{equation} Q_{SS}=\int_0^{H} u_{SS} \,{\rm d} y= 0, \end{equation}

in agreement with previous findings regarding steady streaming in tubes (Hall 1974) and channels (Lo Jacono et al. 2005; Guibert et al. 2010). Using the condition $Q_{SS}=0$ in (3.22) finally yields

(3.25)\begin{equation} \frac{{\rm d} p_{SS}}{{\rm d} x}={-}6 H^{{-}3} \int_0^{H} F y (H-y) \, {\rm d} y, \end{equation}

for the pressure gradient, thereby completing the solution.

3.4. Buoyancy-induced velocity

The corresponding solution for the buoyancy-induced velocity can be obtained by simply replacing $F(x,y)$ with ${{\textit {Ri}}} \, c_0(x,y,\tau )$ in (3.19) and (3.20), yielding

(3.26)\begin{equation} \frac{u_{B}}{\alpha^{2} {{\textit{Ri}}}}={-}\frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x} \frac{1}{2} (H-y)y+y\int_0^{y} c_0 \, {\rm d} \hat{y} -\int_0^{y} c_0 \hat{y} \, {\rm d} \hat{y}-y \int_0^{H} c_0\left(1-\frac{y}{H} \right) {\rm d} {y} \end{equation}

and

(3.27)\begin{align} \frac{v_{B}}{\alpha^{2} {{\textit{Ri}}}}&=\frac{\partial}{\partial x} \left[\frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x} \frac{y^{2}}{2} \left(\frac{H}{2}-\frac{y}{3}\right)+\frac{y^{2}}{2}\int_y^{H} c_0 \left(1-\frac{\hat{y}}{H} \right) {\rm d} \hat{y} \right.\nonumber\\ &\quad \left.+\,y\left(1-\frac{y}{2H}\right)\int_0^{y} c_0 \hat{y} \, {\rm d} \hat{y}-\frac{1}{2} \int_0^{y} c_0 \hat{y}^{2} \, {\rm d} \hat{y}\right]. \end{align}

Using the condition $v_{B}=0$ at $y=H$ in the above equation and integrating once gives

(3.28)\begin{equation} \frac{Q_{B}}{\alpha^{2} {{\textit{Ri}}}}=\frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x} \frac{H^{3}}{12}+\frac{1}{2} \int_0^{H} c_0 y (H-y) \, {\rm d} y \end{equation}

for the buoyancy-induced flow rate $Q_{B}=\int _0^{H} u_{B} \,{\rm d} y$. Integrating a second time with $p_{B}=0$ at $x=0,n$ to give

(3.29)\begin{equation} Q_{B}(\tau) =\frac{\displaystyle\alpha^{2} {{\textit{Ri}}} \int_0^{n} H^{{-}3} \left[\int_0^{H} c_0 y (H-y) \,{\rm d} y \right] {\rm d} x}{\displaystyle 2 n \int_0^{1} H^{{-}3}\, {\rm d} x}, \end{equation}

finally determines the pressure gradient

(3.30)\begin{equation} \frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x}=\frac{6}{H^{3}} \left\{\frac{\displaystyle\int_0^{n} H^{{-}3} \left[\int_0^{H} c_0 y (H-y) \,{\rm d} y \right] {\rm d} x}{\displaystyle n \int_0^{1} H^{{-}3} \,{\rm d} x}-\int_0^{H} c_0 y (H-y) \, {\rm d} y \right\}, \end{equation}

which can be used in (3.26) and (3.27) to complete the determination of the buoyancy-induced velocity. Note that, because $c_0$ is not spatially periodic, the solution carries a dependence on the channel length $n$ through the pressure gradient $\partial p_{B}/\partial x$.

4. Solute dispersion

The flow velocity is coupled with the solute concentration $c$ through the dependence on $c_0$ present in $\boldsymbol{v}_{B}=(u_{B},v_{B})$. The computation of $c_0$ involves substitution of the expansion $c=c_0+\varepsilon c_1 + \varepsilon ^{2} c_2 + \cdots$ into (2.12) with the time derivative replaced by the two-time-scale expression $\partial c/\partial t+\varepsilon ^{2} \partial c/\partial \tau$. At leading order we find the result $\partial c_0/\partial t=0$, anticipated earlier when writing (3.17a,b). Collecting terms of order $\varepsilon$ yields

(4.1)\begin{equation} \frac{\partial c_1}{\partial t}+\boldsymbol{v}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} c_0=0, \end{equation}

which can be integrated to provide the concentration correction

(4.2)\begin{equation} c_1=\langle {c}_1\rangle (x,y,\tau)-\int \boldsymbol{v}_0 \,{\rm d}t \boldsymbol{\cdot} \boldsymbol{\nabla} c_0, \end{equation}

where $\int \boldsymbol{v}_0 \,{\rm d}t={\rm Re}[\mathrm {e}^{\mathrm {i} t}(U,V)]$, as follows from (3.1a–c). It is worth noting that, because the solute diffusivity takes small values of order $\kappa /\nu \sim \varepsilon ^{2}$, effects of diffusion are absent at the first two orders in the asymptotic analysis. These effects are present in the equation that arises at the following order,

(4.3)\begin{equation} \frac{\partial c_2}{\partial t}+\frac{\partial c_0}{\partial \tau}+\boldsymbol{v}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} c_1+\boldsymbol{v}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} c_0=\frac{1}{\alpha^{2} \sigma} \frac{\partial^{2} c_0}{\partial y^{2}}, \end{equation}

which can be time-averaged to give

(4.4)\begin{equation} \frac{\partial c_0}{\partial \tau}+\left(\left\langle \int \boldsymbol{v}_0 \,{\rm d}t \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_0 \right\rangle + \langle \boldsymbol{v}_1 \rangle \right) \boldsymbol{\cdot} \boldsymbol{\nabla} c_0=\frac{1}{\alpha^{2} \sigma} \frac{\partial^{2} c_0}{\partial y^{2}}. \end{equation}

In deriving the second term in (4.4) from the third term in (4.3) use has been made of (4.2). Since $\langle {c}_1\rangle$ is independent of $t$ and $\langle \boldsymbol{v}_0 \rangle =0$, the contribution of the former to the resulting time average $\langle (\boldsymbol{v}_0 \boldsymbol {\cdot} \boldsymbol {\nabla } \langle {c}_1\rangle ) \rangle = \langle \boldsymbol {v}_0 \rangle \boldsymbol {\cdot} \boldsymbol {\nabla } \langle {c}_1\rangle$ is identically zero. The leading-order solute concentration $c_0$ is also independent of $t$, so that the contribution arising from the second term in (4.2) can be written in the form

(4.5)\begin{align} &-\left\langle \boldsymbol{v}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} \left( \int \boldsymbol{v}_0 \,{\rm d}t \boldsymbol{\cdot} \boldsymbol{\nabla} c_0\right) \right\rangle={-}\left\langle \boldsymbol{v}_0 \boldsymbol{\cdot} \int \boldsymbol{\nabla} u_0 \,{\rm d}t \right\rangle \frac{\partial c}{\partial x}-\left\langle \boldsymbol{v}_0 \boldsymbol{\cdot} \int \boldsymbol{\nabla} v_0 \,{\rm d}t \right\rangle \frac{\partial c}{\partial y} \nonumber\\ &\quad -\left\langle u_0 \int u_0 \,{\rm d}t \right\rangle \frac{\partial^{2} c}{\partial x^{2}} -\left\langle v_0 \int u_0 \,{\rm d}t +u_0 \int v_0 \,{\rm d}t \right\rangle \frac{\partial^{2} c}{\partial x \partial y} -\left\langle v_0 \int v_0 \,{\rm d}t \right\rangle \frac{\partial^{2} c}{\partial y^{2}}. \end{align}

With the time averages of any two harmonic functions $A$ and $B$ satisfying $\langle A (\int B \,{\rm d}t)\rangle =-\langle (\int A \,{\rm d}t) B \rangle$ and $\langle A (\int A \,{\rm d}t)\rangle =\langle B (\int B \,{\rm d}t) \rangle =0$, it follows that the terms in the second line of the above equation are identically zero, whereas the remaining two terms on the right-hand side can be cast in the compact form shown in (4.4).

As seen in (4.4), convective transport in the long time scale relies on the time-averaged Lagrangian velocity, given by the sum of the time-averaged Eulerian velocity $\langle \boldsymbol{v}_1 \rangle =\boldsymbol{v}_{SS}+\boldsymbol{v}_{B}$ and the Stokes drift $\boldsymbol{v}_{SD}=(u_{SD},v_{SD})=\langle \int \boldsymbol{v}_0 \,{\rm d}t \boldsymbol {\cdot} \boldsymbol {\nabla } \boldsymbol{v}_0 \rangle$ (see, e.g., Larrieu et al. 2009 for a discussion on Lagrangian transport in a similar wall-bounded flow). The latter contribution can be evaluated in terms of the complex functions $U$ and $V$ with use of the expressions

(4.6a,b)\begin{equation} u_{SD}=\frac{1}{2}{\rm Im}\left[\frac{\partial}{\partial x}(U \bar{U})+ \frac{\partial}{\partial y}(V \bar{U})\right] \quad {\rm and} \quad v_{SD}=\frac{1}{2}{\rm Im}\left[\frac{\partial}{\partial x}(U \bar{V})+ \frac{\partial}{\partial y} (V \bar{V})\right], \end{equation}

which follow from the identity $\langle {\rm Re}(\mathrm {e}^{\mathrm {i} t} f_1) \, {\rm Re}(\mathrm {i} \mathrm {e}^{\mathrm {i} t} f_2) \rangle = {\rm Im}(f_1 \bar {f}_2)/2$. The function $u_{SD}$, which is related to the function $F$ defined earlier in (3.18), is identically zero at $x=0,1/2,1,3/2,\dots$, so that the associated constant volumetric flow rate is simply

(4.7)\begin{equation} Q_{SD}=\int_0^{H} u_{SD} \,{\rm d} y= 0. \end{equation}

As our asymptotic description is limited to the leading-order term in the asymptotic expansion (2.17) for the solute concentration, to summarise the results of the asymptotic analysis one may replace $c_0$ with $c$ when rewriting the final transport (4.4) in the form

(4.8)\begin{equation} \frac{\partial c}{\partial \tau}+(u_{SD}+u_{SS}+u_{B}) \frac{\partial c}{\partial x} +(v_{SD}+v_{SS}+v_{B})\frac{\partial c}{\partial y}=\frac{1}{\alpha^{2} \sigma} \frac{\partial^{2} c}{\partial y^{2}}. \end{equation}

The description of the solute dispersion following its deposition in the channel reduces to the integration of the above equation with initial condition $c=c_i(x,y)$ at $\tau =0$ and boundary conditions $\partial c/\partial y=0$ at $y=0,H$. In the integration, the time-independent Stokes-drift and steady-streaming velocities are computed with use of (4.6a,b) and of (3.19), (3.20) and (3.25), respectively, whereas the time-varying buoyancy-induced velocity is evaluated in terms of the solute concentration through the integral expressions (3.26), (3.27) and (3.30) with $c_0=c$. Observation of (4.8) reveals that gravitational forces modify the character of the solution in a non-trivial way. Owing to the dependence of $u_{B}$ and $v_{B}$ on the concentration distribution $c$, the time-averaged transport equation that governs the dispersion of the solute, which for ${{\textit {Ri}}}=0$ reduces to a linear partial-differential equation with time-independent coefficients, turns into a complicated nonlinear integro-differential equation in the presence of buoyancy.

It is worth noting that, whereas the volumetric flow rates $Q_{SS}=\int _0^{H} u_{SS} \,{\rm d} y$ and $Q_{SD}=\int _0^{H} u_{SD} \,{\rm d} y$ associated with steady streaming and Stokes drift are identically zero, as discussed previously, that induced by buoyancy is in general non-zero, its value $Q_{B}=\int _0^{H} u_{B} \,{\rm d} y$ evolving in time according to (3.29). Note that writing (4.8) in conservative form and integrating across the channel with use of $\partial c/\partial y=0$ and $v_{SD}=v_{SS}=v_{B}=0$ at $y=0,H$ yields the relation

(4.9)\begin{equation} \frac{\partial C}{\partial \tau}+\frac{\partial \phi}{\partial x}=0 \end{equation}

between the amount of solute per unit channel length $C(x,\tau )=\int _0^{H} c \, {\rm d} y$ and the solute flux $\phi (x,\tau )=\int _0^{H} (u_{SD}+u_{SS}+u_{B})c \, {\rm d} y$, whereas a second integral between $x=0$ and $x=1$ gives

(4.10)\begin{equation} \frac{\partial}{\partial \tau}\left(\int_0^{1} C \,{\rm d} x \right)+\phi(1,\tau)-\phi(0,\tau)=0, \end{equation}

which naturally reduces to the expected conservation law

(4.11)\begin{equation} \int_0^{1} \left(\int_0^{H} c \, {\rm d} y\right) {\rm d} x=\int_0^{1} \left(\int_0^{H} c_i \, {\rm d} y\right) {\rm d} x, \end{equation}

when the solute flux at the two ends is zero.

An important aspect of the reduced description developed above is that the nonlinear integro-differential equation (4.8) targets directly the evolution of the flow in the long time scale $\tau \sim 1$, relevant for solute dispersion over distances of order $\lambda$ (i.e. dimensionless distances $x$ of order unity), thereby circumventing the need to describe the small concentration fluctuations occurring in the short time scale $t=\varepsilon ^{-2} \tau$. As a consequence, the model predictions involve computational times that can be expected to be a factor $\varepsilon ^{2}$ smaller than those required in DNS, because to describe solute dispersion the DNS must track the flow over a large number of cycles ${\sim }\varepsilon ^{-2}$, larger for smaller $\varepsilon$.

5. Selected numerical results

The reduced flow description is to be utilised to investigate the influence of buoyancy on solute dispersion. To facilitate the computation, the conservation (4.8) was written in terms of the normalised coordinate $\eta =y/H(x)$, so that the integration domain becomes $0< x< n$ and $0<\eta < 1$. The numerical scheme utilises a fourth-order compact centred finite-difference approximation for the spatial discretisations of the viscous terms and a second-order upwind scheme for the nonlinear terms. A third-order TVD Runge–Kutta scheme is used for the time marching, whereas the integral expressions (3.26), (3.27) and (3.30) are evaluated with a simple trapezoidal rule.

The accuracy of the model predictions, derived in the asymptotic limit $\varepsilon \ll 1$ for a slender channel with $h_o/\lambda \ll 1$, is tested through comparisons with two-dimensional, unsteady simulations of the fluid motion and solute dispersion for small but finite values of $h_o/\lambda$ and $\varepsilon$. The DNS, involving the complete (2.1)–(2.3) written in dimensionless form, span thousands of oscillation cycles, as needed to generate significant dispersion of the solute. The numerical integration was performed with the finite-volume solver Ansys Fluent (release 20.2), assuring second-order accuracy in time and in space. Computations employing upwind and central-differencing schemes for the convective terms were found to yield virtually indistinguishable results, with the former discretisation used in generating the figures shown in the following. A coupled algorithm was used for the pressure–velocity coupling. In addition to the boundary conditions used in integrating the slender flow (2.9)–(2.12), additional conditions of developed flow $\partial u/\partial x=\partial c/\partial x=0$ at the upper and lower open boundaries are incorporated in integrating (2.1)–(2.3). The computations presented in the following correspond to a canal of total length $n=3$ and aspect ratio $h_o/\lambda =1/20$ with the imposed pressure difference yielding a dimensionless stroke length $\varepsilon =0.02$. The time-periodic DNS results were averaged in time to determine the mean Eulerian velocity $\langle \boldsymbol{v} \rangle =({1}/{2 {\rm \pi}}) \int _t^{t+2 {\rm \pi}} \boldsymbol{v} \, {\rm d}t$, of order $\varepsilon$, to be compared with the steady-streaming velocity $\boldsymbol{v}_{SS}$, as explained later. In addition, tracer particles are used to compute the Lagrangian velocity $\boldsymbol{v}_{L}$ by following their displacement over a cycle, i.e. if the particle located at $(x,y)$ at time $t$ moves to occupy the new location $(x+ \delta x,y+ \delta y)$ at time $t + 2{\rm \pi}$, then the Lagrangian velocity at $(x,y)$ and time $t$ is defined as $\boldsymbol{v}_{L}=(\delta x,\delta y)/(2 {\rm \pi})$.

5.1. Buoyancy-free flow

As mentioned previously, in the absence of buoyancy, the flow induced by the imposed pressure gradient is periodic in time and space. The steady Lagrangian motion for $\varepsilon \ll 1$ is given in this case by the sum of the steady-streaming velocity $\boldsymbol{v}_{SS}$ and the Stokes-drift velocity $\boldsymbol{v}_{SD}$. These two contributions as well as their sum are shown in figure 2 for $\beta =0.4$ and two different values of $\alpha$. As the flow in each cell is identical, it suffices to show the solution for $0\le x \le 1$, symmetric with respect to the centre line $x=0.5$. For each value of $\alpha$, streamlines are plotted using a fixed increment $\delta \psi$ of the streamfunction $\psi$, defined in the usual way (e.g. $\partial \psi /\partial y=u_{SS}$ and $\partial \psi /\partial x=-v_{SS}$ for steady streaming) with $\psi =0$ along the wall, so that the interline spacing provides a measure of the local velocity. To further quantify the motion, colour contours are used to represent the associated vorticity $\varOmega =(h_o/\lambda )^{2} \partial v/\partial x-\partial u/\partial y$, which reduces to $\varOmega =-\partial u/\partial y$ in the slender flow approximation.

Figure 2. Streamlines and colour contours of vorticity corresponding to the steady-streaming velocity $\boldsymbol{v}_{SS}$, Stokes-drift velocity $\boldsymbol{v}_{SD}$ and steady mean Lagrangian velocity $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ in a canal with $\beta =0.4$ for (a) $\alpha =4$ and (b) $\alpha =16$. Results of DNS of a non-buoyant flow (i.e. ${{\textit {Ri}}}=0$) with $h_o/\lambda =1/20$ and $\varepsilon =0.02$ are also shown, including the rescaled time-averaged Eulerian velocity field $\langle \boldsymbol{v} \rangle /\varepsilon$ (first column) and the rescaled Lagrangian velocity $\boldsymbol{v}_{L}/\varepsilon$ (fifth column), the latter determined by following tracer particles, as explained in the text. To facilitate the comparisons, fixed constant streamline spacings $\delta \psi =0.003$ and $\delta \psi =0.03$ are used for the upper and lower plots, respectively.

The spatially periodic, time-independent, steady-streaming velocity computed with (3.19) and (3.20) supplemented with (3.25) is shown in the second column of figure 2. The results are qualitatively similar to those presented in Guibert et al. (2010). For $\alpha =4$, the flow structure of each half cell exhibits two counter-rotating vortices, whereas for $\alpha =16$ the flow develops an additional, much weaker vortex, located near the section with largest width. As expected, the vorticity, having peak values of order unity for $\alpha =4$, increases with increasing flow frequency as a result of augmented wall production to reach peak values exceeding $\varOmega =40$ for $\alpha =16$.

The steady-streaming results are compared with time-averaged velocity fields obtained in DNS with $\varepsilon =0.02$. In the comparison, the time-averaged DNS velocity is expressed in the rescaled form $\langle \boldsymbol{v} \rangle /\varepsilon \sim 1$, consistent with the scaling employed in defining $\boldsymbol{v}_{SS}$. The two functions $\boldsymbol{v}_{SS}$ and $\langle \boldsymbol{v} \rangle /\varepsilon$ are seen to be almost identical, thereby giving additional confidence in the mathematical development. For instance, the peak values of the stream function and vorticity corresponding to the time-averaged DNS velocity $\langle \boldsymbol{v} \rangle /\varepsilon$ are $\psi _{PEAK}=\pm (0.0115, 0.1680)$ and $\varOmega _{PEAK}=\pm (1.4465, 40.786)$ for $\alpha =(4,16)$, whereas the corresponding values for the steady-streaming motion are $\psi _{PEAK}=\pm (0.0115,0.1699)$ and $\varOmega _{PEAK}=\pm (1.4474,40.787)$. The small relative differences remain below about 1 %, as is consistent with the order of the asymptotic development.

The third column in figure 2 displays the Stokes-drift velocity field evaluated with (4.6a,b). As it is clear from a quantitative comparison with the corresponding steady-streaming results, both bulk-flow velocities have comparable magnitude for $\alpha =4$, whereas for $\alpha =16$ the Stokes drift provides a much smaller relative contribution to the Lagrangian drift. The dominant role of steady streaming in flows at high Womersley numbers is found also away from the wall in oscillating flow over circular cylinders (Holtsmark et al. 1954; Raney, Corelli & Westervelt 1954), for example. The mean Lagrangian velocity $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ corresponding to the asymptotic limit $\varepsilon \ll 1$ compares favourably with the velocity $\boldsymbol{v}_{L}/\varepsilon$ obtained numerically by following tracer particles in the DNS computation for $\varepsilon =0.02$, shown in the last column of figure 2, although the relative errors are somewhat larger than those of the Eulerian velocity. For instance, the peak values of the stream function and vorticity corresponding to $\boldsymbol{v}_{L}/\varepsilon$ are $\psi _{PEAK}=\pm (0.0235, 0.1674)$ and $\varOmega _{PEAK}=\pm (2.0454, 45.9497)$ for $\alpha =(4,16)$, while the corresponding values for $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ are $\psi _{PEAK}=\pm (0.0235, 0.1491)$ and $\varOmega _{PEAK}=\pm (1.9906, 40.787)$.

5.2. Buoyancy-free solute dispersion

The reduced transport (4.8) resulting from the two-time-scale asymptotic analysis indicates that the solute relies on the Lagrangian drift for longitudinal dispersion. As a consequence, the existence of the closed recirculating vortices displayed in figure 2 implies that in oscillatory buoyancy-free channel flow a solute released in a given cell would be unable to reach their neighbouring cells, thereby precluding its progression along the canal. To illustrate this important feature of the flow, we show in figure 3 the temporal evolution of a bolus of solute with reduced Schmidt number $\sigma =1$ released at the initial instant of time in the central cell of a three-cell canal. The initial concentration is given by the truncated Gaussian distribution

(5.1)\begin{equation} c_i = \min\left\{1,\frac{3}{2}\exp\left[{-}16\left(\frac{x-x_0}{\delta}\right)^{2}\right]\right\}, \end{equation}

which represents a band of solute with characteristic width $\delta$ centred at $x_o$ having a saturated core flanked by thin layers across which the concentration decays to zero. Results obtained by integration of (4.8) for $x_0=1.75$ and $\delta =0.2$ are compared in figure 3 at different instants of time in the interval $0\le \tau \le 8$ with DNS computations. Note that, with $\tau =\varepsilon ^{2} t$, for the value $\varepsilon =0.02$ used in the DNS, this interval of time corresponds to $0\le t \le 20\,000$ (i.e. about $20\,000/2{\rm \pi} \simeq 3200$ oscillatory cycles).

Figure 3. Snapshots of solute concentration for ${{\textit {Ri}}}=0$, $\beta =0.2$, $\alpha =4$ and $\sigma =1$ as obtained at five different instants of time from the reduced model and from DNS for $h_o/\lambda =1/20$ and $\varepsilon =0.02$. In addition to colour contours of local concentration, the figure shows distributions of solute per unit channel length $\int _0^{H} c \,{\rm d} y$ for the model (solid curve) and for the DNS (dot-dashed curves), with the initial distribution $\int _0^{H} c_i \,{\rm d} y$ shown as a dotted curve. For reference, the figure shows streamlines with constant spacing $\delta \psi =0.003$ for the Lagrangian mean drift, which is characterised using the asymptotic prediction $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ for the model results and the value of $\boldsymbol{v}_{L}/\varepsilon$ determined numerically for the DNS results.

As seen in figure 3 the model accurately reproduces the DNS results. To facilitate the quantitative comparisons, in addition to colour contours showing the solute concentration, the figure includes side plots for the amount of solute per unit channel length $C=\int _0^{H} c \, {\rm d} y$ at different instants of time, with the initial distribution $C_i=H(x) c_i(x)$ included for reference as a dotted curve. The model predictions lie very close to the DNS results, in that the normalised value of the integrated departure $(\int _0^{n} |C_{DN}-C_{MODEL}|\,{\rm d} x)/(\int _0^{n} C_i \,{\rm d} x)$, which provides a metric for the accuracy of the model, remains below $0.003$ over the entire range of times considered in the figure.

For the buoyancy-free conditions considered in figure 3, the steady Lagrangian motion is seen to stir the solute about the deposition location, uniformising its concentration within the recirculating cell. The effect of longitudinal diffusion, present in the DNS results, is found to be rather limited, in that, even at the latest instant of time considered, the presence of the solute in the adjacent cells is negligibly small. This tendency of the solute to remain trapped inside Lagrangian vortices has potential implications concerning the drug-dispersion rate in ITDD procedures. Although the Lagrangian flow in the spinal canal does not exhibit the spatial periodicity of the canonical configuration investigated here, closed recirculating vortices, associated with the changes in the eccentricity of the spinal cord along the canal, have been found to characterise the CSF bulk motion (Coenen et al. 2019). Typically, there are three main vortices, extending along the cervical, thoracic and lumbar regions. As ITDD injection occurs in the lumbar region, the buoyancy-free results in figure 3 seem to indicate that, when the density of the drug matches exactly the CSF density, the drug is bound to linger in the lumbar vortex near the injection site without reaching the thoracic region. This could be advantageous in applications involving pain medication, which is meant to be delivered to the spinal cord, but not in applications involving anticancer drugs targeting brain tumors, for example. As shown in the following, buoyancy-induced motion has the potential to drastically change the associated transport rate, in accordance with clinical observations (Mitchell et al. 1988; Povey et al. 1989; Richardson et al. 1996; Veering et al. 2001).

5.3. Slowly varying buoyancy-induced motion

As reasoned previously, buoyancy forces, acting on solutes with density $\rho _s \ne \rho$, alter the steady Lagrangian drift by adding an additional component that varies slowly in the long time scale $\tau$ following the solute dispersion, so that the flow and the solute transport are intimately coupled, as described by (4.8) supplemented with (3.26), (3.27) and (3.30). The corresponding behaviour is characterised in figure 4 for a light solute with ${{\textit {Ri}}}=1$ spreading upwards. Note that, because of the problem symmetry, results corresponding to a heavy solute with ${{\textit {Ri}}}=-1$ can be generated by simply reversing the direction of the gravity vector, i.e. by rotating the figure $180^{\circ }$.

Figure 4. Snapshots of solute concentration for ${{\textit {Ri}}}=1$, $\beta =0.2$, $\alpha =4$ and $\sigma =1$ as obtained at five different instants of time from the reduced model and from DNS for $h_o/\lambda =1/20$ and $\varepsilon =0.02$. In addition to colour contours of local concentration, the figure shows distributions of solute per unit channel length $\int _0^{H} c \,{\rm d} y$ for the model (solid curve) and for the DNS (dot-dashed curves), with the initial distribution $\int _0^{H} c_i \,{\rm d} y$ shown as a dotted curve. The plots include streamlines with constant spacing $\delta \psi =0.003$ for the varying Lagrangian mean drift, which is evaluated with use of $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}+\boldsymbol {v}_{B}$ (model results) and from the displacement of the tracer particles (DNS results).

As in the buoyancy-free flow depicted in figure 3, the solution in figure 4 includes Lagrangian streamlines, colour contours of solute concentration, and streamwise distributions of integrated solute concentration $C=\int _0^{H} c \, {\rm d} y$ along the canal. Buoyancy has a dramatic effect on the dispersion of the solute, as is apparent by comparing the results in both figures. Gravitational forces acting on the light solute induce a longitudinal pressure gradient that modifies drastically the resulting Lagrangian drift, as can be seen by comparing the streamlines in figure 3 with those in figure 4. The pattern of symmetric recirculating cells with unconnected streamlines existing for ${{\textit {Ri}}}=0$ is replaced for ${{\textit {Ri}}}=1$ by a more complicated streamline pattern featuring a net upward flow rate $Q_{B}(\tau )$ (see the solid curves in figure 5, to be discussed later). As can be seen, although the flow rate and the associated streamline pattern vary slowly in time, the observed changes are not very pronounced. Also of interest is that the quantitative agreement between the streamlines predicted by the model and the DNS results is again remarkable, thereby giving additional confidence in our development.

Figure 5. Influence of (a) the contraction ratio $\beta$, (b) reduced Schmidt number $\sigma$, (c) Richardson number ${{\textit {Ri}}}$ and (d) Womersley number $\alpha$ on the temporal evolution of the buoyancy-induced flow rate $Q_{B}$. The values of the parameters in each case are: (a) $\alpha =4$, $\sigma =Ri=1$; (b) $\alpha =4$, $Ri=1$ and $\beta =0.2$; (c) $\alpha =4$, $\sigma =1$ and $\beta =0.2$; (d) $\sigma =Ri=1$ and $\beta =0.2$.

The changes in the Lagrangian motion have a dramatic reflection in the solute dispersion. As shown in figure 4, the solute is transported upwards following the Lagrangian streamlines connecting the cells, enabling its upward progression. The bolus of solute is distorted by the recirculating flow as it travels upwards, driven by the buoyancy-induced draft. The variation with time of the concentration distribution predicted with the model is in excellent agreement with the DNS results. The model is shown to predict not only the mean location of the bolus but also its shape and elongation. The relative departures, measured by $(\int _0^{n} |C_{DNS}-C_{MODEL}|\, {\rm d} x)/(\int _0^{n} C_i \,{\rm d} x)$, remain below $0.009$ over the entire range of times shown in the figure. In assessing the potential benefits of the time-averaged formulation, it is important to emphasise that, although the integration of the integro-differential equation (4.8) over times $\tau \sim 1$ can be completed in a few hours using a laptop computer, generating the DNS results shown in figure 4, spanning over 3000 cardiac cycles, required 10 days in a computational cluster using a total of 72 cores.

5.4. Parametric dependence of the results

The case shown in figure 4, corresponding to $\beta =0.2$, $\alpha =4$, $\sigma =Ri=1$ is used in figures 5 and 6 as a basis to investigate the influence of the different parameters on the dispersion of a buoyant solute. To that end, results are generated with use of (4.8) by modifying one of the four controlling parameters at a time, while keeping the other three fixed at the values selected earlier. Figure 5 shows the variation with time of the buoyancy-induced flow rate $Q_{B}$, whereas figure 6 shows instantaneous solute-concentration distribution and associated Lagrangian streamlines at a fixed time, namely, $\tau =(6,6,6,2)$ for figures 6(a)–6(d), with corresponding results for the base case at these times shown in two of the subpanels of figure 4.

Figure 6. Influence of (a) the contraction ratio $\beta$, (b) reduced Schmidt number $\sigma$, (c) Richardson number ${{\textit {Ri}}}$ and (d) Womersley number $\alpha$ on the solute-concentration distribution and associated Lagrangian streamlines. The snapshots are taken at $\tau =6$ for (a–c) and at $\tau =2$ for (d). The values of the parameters in each case are: (a) $\alpha =4$, $\sigma =Ri=1$; (b) $\alpha =4$, $Ri=1$ and $\beta =0.2$; (c) $\alpha =4$, $\sigma =1$ and $\beta =0.2$; (d) $\sigma =Ri=1$ and $\beta =0.2$. The streamline spacing for the mean Lagrangian velocity is $\delta \psi =0.003$ for (a–c) and $\delta \psi =0.005$ for (d).

We begin by discussing the effect of the channel geometry. As shown in figure 6(a), increasing the undulation of the channel from $\beta =0.1$ to $\beta =0.4$ tends to increase the magnitude of the buoyancy-induced longitudinal velocity in the region of minimum cross-sectional area, where streamlines are closely spaced together for larger $\beta$, but these changes result in only moderately small variations of the flow rate $Q_{B}$, as shown in figure 5(a). As a result, the bolus of solute becomes more elongated as $\beta$ increases, but advances upward at approximately the same rate, so that the maximum concentration occupies approximately the same location at $\tau =6$, as shown in figure 6(a).

The effect of the Schmidt number $\sigma$, entering in the formulation only through the factor affecting the transverse diffusion rate on the right-hand side of (4.8), is investigated in figures 5(b) and 6(b). The changes observed in streamline pattern and flow rate when changing the Schmidt number from $\sigma =0.25$ to $\sigma =8$ are not very significant, so that the differences in solute evolution in figure 6(b) are attributable to the direct effect of transverse diffusion (or lack thereof). The snapshot corresponding to $\sigma =8$ displays the behaviour expected at large Schmidt numbers, for which fluid particles maintain a nearly constant concentration in their slow Lagrangian evolution, as described by the limiting form of (4.8) for $\sigma \gg 1$. In the opposite limit $\sigma \ll 1$, solute diffusion leads to a rapid uniformisation of the concentration, as can be seen by integrating $\partial ^{2} c/\partial y^{2}=0$ (i.e. the reduced form of (4.8) when $\sigma \ll 1$) subject to the non-permeability condition $\partial c/\partial y=0$ at $y=0,H$ to give $c=c(x,\tau )$. As a result, the bolus remains relatively compact as it moves along the channel with a velocity determined by the flow rate. The reduced transport equation governing the transport of solutes with $\sigma \ll 1$ can be derived from (4.9) by noting that in this limit the integrated solute concentration $C=\int _0^{H} c \, {\rm d} y$ becomes $C=H(x) c(x,\tau )$ while the solute flux $\phi (x,\tau )=\int _0^{H} (u_{SD}+u_{SS}+u_{B})c \, {\rm d} y$ reduces to $\phi =Q_{B} c$. Substituting these simplified expressions into (4.9) and using (3.29) to evaluate $Q_{B}$ finally yields

(5.2)\begin{equation} \frac{\partial c}{\partial \tau}+\frac{\alpha^{2} {{\textit{Ri}}} \displaystyle\int_0^{n} c \,{\rm d} x}{12 n H \displaystyle\int_0^{1} H^{{-}3} \,{\rm d} x} \frac{\partial c}{\partial x}=0 \end{equation}

as the limiting form of (4.8) for $\sigma \ll 1$.

As is to be expected from the velocity expressions derived in § 3.4, the Richardson number and the Womersley number have a pronounced effect on the mean Lagrangian motion. As shown in figures 5(c) and 5(d), the flow rate exhibits dependences on $\textit{Ri}$ and $\alpha$ that are approximately linear and approximately quadratic, respectively, consistent with (3.29). These dependences have a reflection on the evolution of the solute bolus shown in figures 6(c) and 6(d). With limited updraft for ${{\textit {Ri}}}=0.25$, the bolus is seen to spread about the injection location, without significant upward progression at the instant of time $\tau =6$ considered in the figure. An increase in ${{\textit {Ri}}}$ promotes the displacement of the bolus, but its longitudinal extent remains approximately equal in all three cases. By way of contrast, an increase in $\alpha$ increases the upward displacement and also enhances bolus distortion. The reason for the latter is that larger values of $\alpha$ hinder transverse diffusion, as can be inferred from (4.8), with the result that fluid particles travel following the Lagrangian recirculating paths with a nearly constant concentration, rapidly deforming the compact concentration distribution of the initial bolus.

6. Conclusions

Solute dispersion in a wavy-walled vertical channel subject to an oscillating pressure gradient has been used as a canonical model to investigate the effect of buoyancy on the transport of ITDD drugs, characterised by large values of the Schmidt number and order-unity values of the Richardson number. The mean Lagrangian velocity determined in the asymptotic limit of small stroke lengths, responsible for the convective transport of the solute, displays a buoyancy component whose local value depends on the solute concentration through integral expressions, resulting in a nonlinear integro-differential transport equation. The predictive capabilities of the reduced description are tested through comparisons with DNS computations involving thousands of oscillating cycles. The validation exercise reveals that the model provides accurate predictions of solute dispersion at a fraction of the computational cost involved in the DNS. In contrast to the motion observed in the buoyancy-free case investigated in figures 2 and 3, characterised by the existence of a series of closed Lagrangian vortices distributed periodically along the channel, the buoyancy-modulated mean Lagrangian flow shown in figure 4 includes a streamwise draft connecting neighbouring cells that promotes the longitudinal dispersion of the solute. The buoyancy-enhanced transport rate revealed in our channel computations is consistent with previous clinical observations pertaining to dispersion of light drugs for patients in a sitting injection position (Richardson et al. 1996).

The simple canonical flow considered here has served to unveil some of the key aspects of the solute-dispersion problem, including the enhanced transport associated with the buoyancy-modulated mean Lagrangian velocity. Future work should consider application of the two-time-scale asymptotic analysis delineated previously to the description of ITDD processes, with account taken of the three-dimensional morphology of the spinal canal, possibly including the effect of microanatomical features such as arachnoid trabeculae, which are thin strands of connective tissue that form a web-like structure stretching across the spinal canal. The presence of these fine anatomical structures, which has been shown to have an important effect on pressure loss (Tangen et al. 2015), can be accounted for by treating the spinal subarachnoid space as a Brinkman porous medium, as done in previous investigations (Gupta et al. 2009; Kurtcuoglu, Jain & Martin 2019; Salerno, Cardillo & Camporeale 2020; Sincomb et al. 2022).

The future developments envisioned here can potentially provide a reduced transport equation, possibly similar to (4.8), to be used in combination with magnetic resonance imaging characterisations of the canal anatomy (Coenen et al. 2019) to describe the transport of the drug in the relevant dispersion time scale. The ultimate goal of such efforts is the development of computationally effective subject-specific predictive tools for drug delivery to a target site from injection by a lumbar puncture with account taken of the specific anatomy and physiological conditions of the individual patient as well as for the molecular characteristics and injection rate of the drug, as needed in guiding clinical treatments.