2. Flow and transport in the spinal canal

The SAS surrounding the spinal cord can be described in the first approximation as a thin annular channel whose characteristic width $h_c \sim 0.1\unicode{x2013}0.4\ {\rm cm}$ that is much smaller than the characteristic spinal cord perimeter $\ell _c \sim 2\unicode{x2013}3\ {\rm cm}$, which in turn is much smaller than the spine length $L \sim 60$ cm, so that the canal dimensions satisfy the inequalities $L \gg \ell _c \gg h_c$. The CSF moves along the canal with an oscillatory velocity that is synchronized with the cardiac and respiratory cycles. The CSF oscillatory flow is more pronounced near the canal entrance, where the characteristic velocities $u_c$ are of the order of a few ${\rm cm}\ {\rm s}^{-1}$, but become progressively smaller on approaching the closed end of the canal, as revealed by in vivo magnetic resonance measurements (Haughton & Mardal Reference Haughton and Mardal2014; Aktas et al. Reference Aktas, Kollmeier, Joseph, Merboldt, Ludwig, Gärtner, Frahm and Dreha-Kulaczewski2019; Coenen et al. Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019; Sincomb et al. Reference Sincomb, Coenen, Gutiérrez-Montes, Martínez Bazán, Haughton and Sánchez2022). The following analysis focuses specifically on the flow induced by the cardiac cycle, corresponding to angular frequencies $\omega \simeq 2 {\rm \pi}\ {\rm s}^{-1}$ and characteristic stroke lengths $L_s = u_c/\omega \sim 1$ cm that are much smaller than the canal length $L$.

The motion in the spinal canal is viscous, in that the characteristic viscous time across the canal $h_c^2/\nu$ based on the CSF kinematic viscosity $\nu \simeq 0.7 \times 10^{-3}\ {\rm cm}^2\ {\rm s}^{-1}$ is comparable to – although somewhat larger than – the characteristic flow oscillation time $\omega ^{-1}$, resulting in order-unity values $3 \lesssim \alpha \lesssim 12$ of the Womersley number $\alpha =(h_c^2 \omega /\nu )^{1/2}$. By way of contrast, effects of inertia associated with convective acceleration are very limited, as measured by the relevant Strouhal number $\omega L/u_c=L/L_s \gg 1$, the inverse of which defines an asymptotically small parameter $\varepsilon \sim L_s/L\simeq 0.02\unicode{x2013}0.04$. Thus, in the first approximation, the motion in the slender spinal canal is given by a balance between the pressure gradient, the local acceleration and the viscous forces. The resulting linear unsteady lubrication problem can be solved to give closed-form expressions for the leading-order oscillatory velocity (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), whose time-averaged value is identically zero. Corrections to this solution can be obtained by extending the asymptotic analysis to higher orders in $\varepsilon \ll 1$ (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019). The first-order velocity corrections, of order $\varepsilon u_c$, exhibit non-zero time-averaged values. This steady-streaming velocity, first identified in the seminal computational work of Kuttler et al. (Reference Kuttler, Dimke, Kern, Helmlinger, Stanski and Finelli2010), is due partly to the effect of convective acceleration and partly to the canal compliance (see e.g. Bhosale, Parthasarathy & Gazzola (Reference Bhosale, Parthasarathy and Gazzola2022a), Bhosale et al. (Reference Bhosale, Vishwanathan, Upadhyay, Parthasarathy, Juarez and Gazzola2022b) and Cui, Bhosale & Gazzola (Reference Cui, Bhosale and Gazzola2024) for recent analyses of steady-streaming flows stemming from boundary compliance). The associated residence times for the bulk flow in the canal $L/(\varepsilon u_c)= \varepsilon ^{-2} \omega ^{-1} \sim 30$ min are of the order of those observed in in vivo experiments employing radioactive tracers to mark the displacement of the CSF particles (Di Chiro Reference Di Chiro1964; Greitz & Hannerz Reference Greitz and Hannerz1996).

As shown by Lawrence et al. (Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), the disparity between the short time $\omega ^{-1}$ characterizing the oscillatory velocity fluctuations and the residence time $\varepsilon ^{-2} \omega ^{-1}$ associated with the bulk motion can be used in deriving a simplified transport equation for the drug. The analysis revealed that shear-enhanced diffusion (Watson Reference Watson1983), which is potentially important for solutes with order-unity values of the Schmidt number $S=\nu /\kappa$, is entirely negligible for the large Schmidt numbers $S \gg 1$ corresponding to the small molecular diffusivities $\kappa$ of typical ITDD drugs (e.g. for methotrexate, $\kappa = 5.26 \times 10^{-10}\ {\rm m}^2\ {\rm s}^{-1}$, yielding $S\simeq 1330$ for $\nu =0.7 \times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$). The evolution of the drug concentration in the long time scale $\varepsilon ^{-2} \omega ^{-1}$ was found to be governed by a transport equation involving molecular diffusion across the width of the canal and convective transport driven by the time-averaged Lagrangian motion resulting from the combined effects of steady streaming and Stokes drift. The use of this simplified equation effectively circumvents the need to describe the small concentration fluctuations occurring in the short time scale $\omega ^{-1}$, thereby drastically reducing computational times. The accuracy and limitations of this time-averaged description have been tested recently by means of comparisons with results of direct numerical simulations spanning hundreds of oscillation cycles (Gutiérrez-Montes et al. Reference Gutiérrez-Montes, Coenen, Lawrence, Martínez-Bazán, Sánchez and Lasheras2021), as needed to generate significant dispersion of the solute. The comparisons demonstrate clearly the accuracy of the time-averaged description, which is seen to provide excellent fidelity at a fraction of the computational cost involved in the direct numerical simulations. The present investigation extends our previous analyses of flow and transport in the spinal canal by accounting for the effects of the small density differences between the drug and the CSF. The mathematical development parallels that employed recently in our analysis of buoyant Lagrangian drift in a vertical wavy-walled channel (Alaminos-Quesada et al. Reference Alaminos-Quesada, Coenen, Gutiérrez-Montes and Sánchez2022).

3. Problem description

3.1. The Richardson number

As can be seen in table 1, the drug density $\rho _d$ of common intrathecal drug solutions is very close to that of the CSF ($\rho =1.00059\ {\rm g}\ {\rm cm}^{-3}$ at $37\,^{\circ }{\rm C}$) (Nicol & Holdcroft Reference Nicol and Holdcroft1992; Lui, Polis & Cicutti Reference Lui, Polis and Cicutti1998; McLeod Reference McLeod2004; Hejtmanek, Harvey & Bernards Reference Hejtmanek, Harvey and Bernards2011; Lynch Reference Lynch2014). The drug density can be modified by adding different diluents such as saline, glucose and dextrose. Even though the resulting relative differences are very small (i.e. $10^{-4}\lesssim |\rho -\rho _d|/\rho \lesssim 10^{-2}$), the associated buoyancy forces affect in a fundamental way the dispersion of the drug. Thus it has been seen that for hyperbaric (i.e. dense) drugs, the transport of the drug is restricted when the patient is seated for some time before moving to a supine position (Mitchell et al. Reference Mitchell, Bowler, Scott and Edström1988; Povey, Jacobsen & Westergaard-Nielsen Reference Povey, Jacobsen and Westergaard-Nielsen1989; Veering et al. Reference Veering, Immink-Speet, Burm, Stienstra and van Kleef2001; Loubert et al. Reference Loubert, Hallworth, Fernando, Columb, Patel, Sarang and Sodhi2011). Conversely, when a hypobaric (light) drug is injected, faster cephalic dispersion occurs in a seated injection position than in a lateral injection position (Richardson et al. Reference Richardson, Thakur, Abramowicz and Wissler1996). As expected, the density of the drug is inconsequential when injection occurs in the lateral position (Hallworth, Fernando & Columb Reference Hallworth, Fernando, Columb and Stocks2005) or when the solution density matches that of CSF (Wildsmith et al. Reference Wildsmith, McClure, Brown and Scott1981).

Table 1. A few common intrathecal drugs, with their densities (Nicol & Holdcroft Reference Nicol and Holdcroft1992; Lui et al. Reference Lui, Polis and Cicutti1998; Hejtmanek et al. Reference Hejtmanek, Harvey and Bernards2011) and associated Richardson numbers ${\textit {Ri}} = [g (\rho -\rho _d)]/(\rho \varepsilon ^2 \omega ^2 L)$, the latter evaluated with $g=9.81\ {\rm m}\ {\rm s}^{-2}$, $L=0.6$ m and $\rho =1.00059\ {\rm g}\ {\rm cm}^{-3}$ for two different values of the reduced stroke length $\varepsilon$.

To anticipate how the presence of buoyancy forces modifies drug dispersion for patients in a sitting or upright position, it is useful to compare the characteristic value of the buoyancy-induced acceleration $g (\rho -\rho _d)/\rho$ with the characteristic value of the convective acceleration along the canal $u_c^2/L$, their ratio defining the relevant Richardson number:

(3.1)\begin{equation} {\textit{Ri}} = \frac{g (\rho-\rho_d)/\rho}{u_c^2/L}=\frac{g (\rho-\rho_d)/\rho}{\varepsilon^2 \omega^2 L}. \end{equation}

Typical values of this number are evaluated in table 1 for a few common intrathecal drugs and two different values of the reduced stroke length $\varepsilon$. As can be seen, values of ${\textit {Ri}}$ of order unity characterize most situations of practical interest, so that in ITDD processes, buoyancy acceleration can be anticipated to be comparable to convective acceleration. As discussed previously, the motion of CSF at leading order is given by an unsteady lubrication balance involving the local acceleration and the viscous and pressure forces, with convective acceleration introducing small corrections of order $\varepsilon$, responsible for the steady-streaming motion. This leading-order balance is not altered in the relevant limit ${\textit {Ri}} \sim 1$ that applies to intrathecal drugs, in which the associated buoyancy-induced velocities are comparable to the steady-streaming velocities (and therefore a factor $\varepsilon$ smaller than the pulsating velocities).

3.2. The model problem

The problem is formulated in dimensionless form using the scales and notation employed in the previous buoyancy-free analysis of Lawrence et al. (Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), which can be consulted for details of the derivation. Attention is focused on the motion driven by the periodic intracranial pressure fluctuations associated with the arterial blood flow, to be described for simplicity with the simple sinusoidal function $(\Delta p)_c \cos (\omega t')$, where $(\Delta p)_c$ is the fluctuation amplitude and $\omega \simeq 2 {\rm \pi}\ {\rm s}^{-1}$ is the angular frequency of the cardiac cycle, with $t'$ representing the time. The spinal SAS is modelled as an annular canal bounded internally by the pia mater, surrounding the spinal cord, and externally by the dura membrane. The canal is compliant because of the presence of fatty tissue and venous blood. The displacement of the dura membrane at a given location is assumed to be equal to the product of the local pressure fluctuation and a compliance factor $\gamma '$ that may vary along the canal. Its mean value $\gamma '_c$ can be used to estimate the characteristic value of the dura displacement $\gamma '_c (\Delta p)_c$, which is much smaller than the canal width, with the ratio

(3.2)\begin{equation} \varepsilon=\frac{\gamma'_c (\Delta p)_c}{h_c} \sim \frac{L_s}{L} \end{equation}

defining the small asymptotic parameter representing the dimensionless stroke length.

As indicated in figure 1(c), the problem is described in terms of curvilinear coordinates, including the longitudinal distance to the canal entrance $x$ (scaled with $L$), the transverse distance from the spinal cord $y$ (scaled with the characteristic canal width $h_c$), and the azimuthal distance $s$ (scaled with the local spinal cord perimeter, so that $0\leqslant s \leqslant 1$). The corresponding streamwise, transverse and azimuthal velocity components $(u,v,w)$ are scaled with their characteristic values $u_c=\varepsilon \omega L$, $v_c=\varepsilon \omega h_c$ and $w_c=\varepsilon \omega \ell _c$, the last two of which follow from continuity. The geometry of the canal is defined by the dimensionless unperturbed canal width $\bar {h}(x,s)$ (scaled with $h_c$) and spinal cord perimeter $\ell (x)$ (scaled with $\ell _c$). The linear elastic equation for the canal takes the form

(3.3)\begin{equation} h'=\gamma (\cos t +k^2 p'), \end{equation}

where $h'$ is the dura–membrane displacement (scaled with $\varepsilon h_c$), $t=\omega t'$ is the dimensionless time, $p'(x,t)$ is the streamwise pressure variation (scaled with $\rho u_c \omega L$), $k=L\omega /[(h_c/\gamma _c')/\rho ]^{1/2}$ is a dimensionless elastic wavenumber, and $\gamma (x)=\gamma '/\gamma _c'$ is a dimensionless function describing the streamwise variation of the canal compliance.

3.3. Dimensionless formulation

In the thin-film approximation that applies in the limit $L \gg \ell _c \gg h_c$, the continuity, momentum and solute conservation equations take the simplified form

(3.4)$$\begin{gather} \frac{1}{\ell}\,\frac{\partial}{\partial x} (\ell u) +\frac{\partial v}{\partial y} + \frac{1}{\ell}\,\frac{\partial w}{\partial s}=0, \end{gather}$$

(3.5)$$\begin{gather}\frac{\partial u}{\partial t}+\varepsilon \left[u\,\frac{\partial u}{\partial x}+v\,\frac{\partial u}{\partial y}+\frac{w}{\ell}\, \frac{\partial u}{\partial s} \right]={-}\frac{\partial p'}{\partial x}+\frac{1}{\alpha^2}\,\frac{\partial^2 u}{\partial y^2} - \varepsilon \,{\textit{Ri}} \, c, \end{gather}$$

(3.6)$$\begin{gather}\frac{\partial w}{\partial t}+\varepsilon \left[\frac{u}{\ell}\,\frac{\partial}{\partial x}(\ell w)+v\,\frac{\partial w}{\partial y}+ \frac{w}{\ell}\,\frac{\partial w}{\partial s} \right]={-}\frac{1}{\ell}\,\frac{\partial \hat{p}}{\partial s}+\frac{1}{\alpha^2}\, \frac{\partial^2 w}{\partial y^2}, \end{gather}$$

(3.7)$$\begin{gather}\frac{\partial c}{\partial t}+\varepsilon \left(u\,\frac{\partial c}{\partial x}+v\,\frac{\partial c}{\partial y}+ \frac{w}{\ell}\,\frac{\partial c}{\partial s}\right)=\frac{\varepsilon^2}{\alpha^2 \sigma}\, \frac{\partial^2 c}{\partial y^2}, \end{gather}$$

where $c$ is the drug concentration and $\hat {p}$ is an auxiliary function describing the azimuthal pressure variations. The problem has been formulated using the Boussinesq approximation, as is appropriate for $|\rho -\rho _d| \ll \rho$. Since the spinal curvature is relatively small, for the case of a sitting patient considered here the streamwise coordinate $x$ is practically aligned with the vertical direction, so that the component of the buoyancy force acting in the azimuthal direction is small, and has been correspondingly neglected in writing (3.6). With the definition (3.1), the Richardson number ${\textit {Ri}}$ measuring the buoyancy force in (3.5) is positive/negative when the drug is lighter/heavier than the CSF, buoyancy driving the drug upwards/downwards, in the negative/positive $x$ direction. Following Lawrence et al. (Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), the diffusion term in (3.7) has been written in terms of the reduced Schmidt number $\sigma =\varepsilon ^2 S$, assumed to be of order unity, as is consistent with the values $S \sim 2000$ and $\varepsilon \sim 0.02\unicode{x2013}0.04$ that characterize drug dispersion in the spinal canal.

The velocity satisfies the no-slip condition $u=v=w=0$ at $y=0$, and $u=v-\partial h'/\partial t=w=0$ at $y=h$. Although drug uptake by the spinal nerve as well as through the dura membrane could be incorporated in the model by accounting for non-zero diffusive fluxes at the boundary, for simplicity the following analysis is restricted to non-permeable bounding surfaces, for which the boundary condition for the concentration reduces to $\partial c/\partial \eta =0$ at $y=0,h$. The pressure drop is negligible at the entrance of the canal, resulting in the condition $p'=0$ at $x=0$. The requirement that the axial volume flux $\int _0^1 ( \int _0^h u \,{\rm d}{y}) \,{\rm d} {s}$ must vanish at the closed end $x=1$ completes the set of boundary conditions needed to determine the flow in the canal.

Besides the Richardson number ${\textit {Ri}}$ defined in (3.1) and the compliance parameter $\varepsilon \ll 1$ defined in (3.2), the set of governing parameters includes the Womersley number $\alpha =h_c/(\nu /\omega )^{1/2}$, the dimensionless elastic wavenumber $k=L\omega /[(h_c/\gamma _c')/\rho ]^{1/2}$, and the rescaled Schmidt number $\sigma =S \varepsilon ^2$. The problem is to be solved in the limit $\varepsilon \ll 1$, with $\alpha \sim 1$ and $k \sim 1$, as is appropriate for describing CSF flow in the spinal canal, for solutes with $\sigma =S \varepsilon ^2 \sim 1$ and ${\textit {Ri}} \sim 1$, the distinguished limit of interest in intrathecal drug dispersion.

4. Solute transport in the presence of buoyancy

Following our previous analyses (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019; Alaminos-Quesada et al. Reference Alaminos-Quesada, Coenen, Gutiérrez-Montes and Sánchez2022), the problem defined above is solved by expressing the different variables as expansions in powers of $\varepsilon$ (e.g. $u=u_0+\varepsilon u_1+ \cdots$) and solving sequentially the equations that arise when collecting terms at different orders in $\varepsilon$. In the development, it is convenient to replace the transverse coordinate $y$ by its normalized counterpart $\eta =y/h$, with $0 \leqslant \eta \leqslant 1$. The velocity field depends on the solute concentration through the buoyancy term appearing in (3.5), although the dependence is weak, since $\varepsilon \ll 1$. The distribution of $c$ can be anticipated to vary over times of the order of the residence time associated with the bulk motion $\varepsilon ^{-2} \omega ^{-1}$, inducing slow changes in the velocity, to be described below by introducing the long time scale $\tau =\varepsilon ^2 t$ as an additional independent variable. In this two-time scale formalism, all variables are assumed to be $2 {\rm \pi}$ periodic in the short time scale $t$, slow changes in time being described by the additional time variable $\tau$, which is formally introduced in the equations by replacing the original time derivatives by $\partial /\partial t+ \varepsilon ^2\,\partial /\partial \tau$.

4.1. Leading-order solution

At leading order in the limit $\varepsilon \ll 1$, the problem reduces to the integration of

(4.1)$$\begin{gather} \frac{1}{\ell}\,\frac{\partial}{\partial x} (\ell u_0) - \frac{\eta}{\bar{h}}\,\frac{\partial \bar{h}}{\partial x}\,\frac{\partial u_0}{\partial \eta} + \frac{1}{\bar{h}}\,\frac{\partial v_0}{\partial \eta} + \frac{1}{\ell}\,\frac{\partial w_0}{\partial s}- \frac{\eta}{\bar{h}}\,\frac{1}{\ell}\,\frac{\partial \bar{h}}{\partial s}\,\frac{\partial w_0}{\partial \eta}=0, \end{gather}$$

(4.2)$$\begin{gather}\frac{\partial u_0}{\partial t}={-}\frac{\partial p'_0}{\partial x}+\frac{1}{\alpha^2 \bar{h}^2}\,\frac{\partial^2 u_0}{\partial \eta^2}, \end{gather}$$

(4.3)$$\begin{gather}\frac{\partial w_0}{\partial t}={-}\frac{1}{\ell}\,\frac{\partial \hat{p}_0}{\partial s}+\frac{1}{\alpha^2 \bar{h}^2}\,\frac{\partial^2 w_0}{\partial \eta^2}, \end{gather}$$

(4.4)$$\begin{gather}\frac{\partial c_0}{\partial t}=0 , \end{gather}$$

supplemented with $h'_0=\gamma (\cos t+k^2 p'_0)$, the leading-order form of (3.3), with boundary conditions $u_0=v_0=w_0=\partial c_0/\partial \eta =0$ at $\eta =0$, and $u_0=v_0-\partial h'_0/\partial t=w_0=\partial c_0/\partial \eta =0$ at $\eta =1$, $p'_0=0$ at $x=0$, and $\int _0^1 ( \bar {h} \int _0^1 u_0 \,{\rm d} \eta ) \,{\rm d} {s}=0$ at $x=1$. As indicated by (4.4), at leading order the solute concentration varies only in the long time scale $\tau$, variations with the short time scale $t$ affecting only higher-order corrections of relative order $\varepsilon$ and smaller. As shown previously (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018), the solution to the periodic linear lubrication problem (4.1)–(4.3) can be written as

(4.5)\begin{equation} \left.\begin{gathered}u_0={\rm Re}\left(\mathrm{i}\,\mathrm{e}^{\mathrm{i} t} U\right),\quad v_0={\rm Re}\left(\mathrm{i}\,\mathrm{e}^{\mathrm{i} t} V\right),\quad w_0={\rm Re} \left(\mathrm{i}\,\mathrm{e}^{\mathrm{i} t} W\right), \\ p'_0={\rm Re}\left(\mathrm{e}^{\mathrm{i} t} P' \right),\quad \hat{p}_0={\rm Re} \left(\mathrm{e}^{\mathrm{i} t} \hat{P} \right),\quad h'_0={\rm Re}\left(\mathrm{e}^{\mathrm{i} t} H' \right), \end{gathered}\right\} \end{equation}

where the complex functions $U(x,\eta,s)$, $V(x,\eta,s)$, $W(x,\eta,s)$, $P'(x)$, $\hat {P}(x,s)$ and $H'(x,s)$ are given in Appendix A for completeness. The leading-order solution (4.5), identical to that found in our earlier analyses (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), is buoyancy-free, and therefore independent of the long time scale $\tau$. Buoyancy will be seen to enter at the following order to modify the bulk motion.

4.2. Time-averaged Eulerian velocity

While the above harmonic functions (4.5) have zero mean values over an oscillation period, i.e. $\langle u_0 \rangle =0$, with $\langle\ \cdot\ \rangle =\int _t^{t+2{\rm \pi} } \cdot \,\textrm {d} t/(2{\rm \pi} )$, the velocity corrections $(u_1,v_1,w_1)$ contain non-zero cycle-averaged components $(\langle u_1 \rangle,\langle v_1 \rangle,\langle w_1 \rangle )$ that satisfy the quasi-steady conservation equations

(4.6)$$\begin{gather} \mathcal{F}=\frac{1}{\ell}\,\frac{\partial}{\partial x} (\ell \bar{h} \langle u_1 \rangle)+ \frac{1}{\ell}\, \frac{\partial }{\partial s}(\bar{h} \langle w_1 \rangle)-\frac{\partial}{\partial \eta} \left(\eta\,\frac{\partial \bar{h}}{\partial x} \langle u_1 \rangle + \frac{\eta}{\ell}\,\frac{\partial \bar{h}}{\partial s} \langle w_1 \rangle\right) + \frac{\partial \langle v_1 \rangle}{\partial \eta}, \end{gather}$$

(4.7)$$\begin{gather}\mathcal{F}_x ={-} \frac{\partial \langle p'_1 \rangle}{\partial x}+\frac{1}{\bar{h}^2 \alpha^2}\,\frac{\partial^2 \langle u_1 \rangle}{\partial \eta^2} - {\textit{Ri}} \, c_0, \end{gather}$$

(4.8)$$\begin{gather}\mathcal{F}_s ={-} \frac{1}{\ell}\,\frac{\partial \langle \hat{p}_1 \rangle}{\partial s}+\frac{1}{\bar{h}^2 \alpha^2}\,\frac{\partial^2 \langle w_1 \rangle}{\partial \eta^2}, \end{gather}$$

obtained by taking the time average of the equations that emerge when collecting terms of order $\varepsilon$ in (3.4)–(3.6). The functions $\mathcal {F}$, $\mathcal {F}_x$ and $\mathcal {F}_s$ appearing on the left-hand side of the above equations carry the combined effects of convective acceleration and canal deformation on the mean Eulerian motion. These functions involve time averages of products of the harmonic functions (4.5), with expressions given in Appendix A.

The velocity must satisfy the homogeneous boundary conditions $\langle u_1 \rangle =\langle v_1 \rangle =\langle w_1 \rangle =0$ at $\eta =(0,1)$ and $\int _0^1 \bigl (\bar {h} \int _0^1 \langle u_1 \rangle \,\textrm {d} {\eta }\bigr ) \,\textrm {d} {s}=0$ at $x=1$. Note that the condition $\langle v_1 \rangle =0$ at $\eta =1$ follows at this order from the general condition $v=\partial h'/\partial t$ written in the two-time scale formalism in the form $v=\partial h'/\partial t+ \varepsilon ^2\,\partial h'/\partial \tau$, so that $\langle v \rangle =\varepsilon ^2\,\partial \langle h' \rangle /\partial \tau$.

Observation of (4.6)–(4.8) reveals that the mean Eulerian motion has two different driving mechanisms, namely, the buoyancy force $- {\textit {Ri}} \, c_0$ appearing on the right-hand side of (4.7), which varies slowly in the long time scale $\tau$, and the steady functions $\mathcal {F}$, $\mathcal {F}_x$ and $\mathcal {F}_s$, associated with convective acceleration and canal deformation. Since the problem is linear, the two distinct driving mechanisms can be quantified separately by expressing the mean Eulerian velocity $(\langle u_1 \rangle,\langle v_1 \rangle,\langle w_1 \rangle )=(u_{ {SS}}+u_{ {B}},v_{ {SS}}+v_{ {B}},w_{ { SS}}+w_{ {B}})$ as the sum of the steady-streaming velocity $(u_{ {SS}},v_{ {SS}},w_{ {SS}})$ and the buoyancy-induced drift $(u_{ {B}},v_{ {B}},w_{ {B}})$. The former was obtained in our previous analyses (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019) by integration of the problem arising with ${\textit {Ri}}=0$, yielding the solution given in Appendix A, while the latter, the new contribution arising when the drug density differs from the CSF density (i.e. when ${\textit {Ri}} \ne 0$), can be obtained by integration of the reduced problem corresponding to $\mathcal {F}=\mathcal {F}_x=\mathcal {F}_s=0$. The resulting solution, involving integrals of the leading-order solute concentration $c_0$, can be cast in the form

(4.9)$$\begin{gather} \frac{u_{{B}}}{\alpha^2\,{\textit{Ri}} \,\bar{h}^2} = 3 \eta (1-\eta)\,\frac{\displaystyle\int\nolimits_0^1 \bar{h}^3 \mathcal{C} \,{\rm d}s} {\displaystyle\int\nolimits_0^1 \bar{h}^3 \,{\rm d} s} + \eta \int_0^\eta c_0 \,{\rm d} \tilde{\eta} - \int_0^\eta c_0 \tilde{\eta} \,{\rm d} \tilde{\eta} - \eta \int_0^1 c_0 (1-{\eta}) \, {\rm d} {\eta}, \end{gather}$$

(4.10)$$\begin{gather}\frac{w_{{B}}}{\alpha^2\,{\textit{Ri}} \,\bar{h}^2} = \frac{3 \eta (1-\eta)}{\bar{h}^3}\, \frac{\partial}{\partial x} \left[\ell \left(\int_0^s \bar{h}^3 \mathcal{C}\, {\rm d} \tilde{s} - \frac{\displaystyle\int\nolimits_0^1 \bar{h}^3 \mathcal{C} \,{\rm d} s}{\displaystyle\int\nolimits_0^1 \bar{h}^3 \,{\rm d} s} \displaystyle\int\nolimits_0^s \bar{h}^3 \,{\rm d} \tilde{s}\right) \right], \end{gather}$$

(4.11)$$\begin{gather}\frac{v_{{B}}}{\alpha^2\,{\textit{Ri}}} = \frac{\eta^2}{\ell}\left(\eta- \frac{3}{2}\right)\frac{\partial}{\partial x} \left( \ell \bar{h}^3 \mathcal{C} \right) -\frac{1}{\ell}\,\frac{\partial}{\partial x} \left( \ell \bar{h}^3 f_{{B}}\right) + \eta\,\frac{\partial \bar{h}}{\partial x}\,\frac{u_{{B}}}{\alpha^2\,{\textit{Ri}}} + \frac{\eta}{\ell}\,\frac{\partial \bar{h}}{\partial s}\,\frac{w_{{B}}}{\alpha^2\,{\textit{Ri}}}, \end{gather}$$

where

(4.12)\begin{equation} \mathcal{C}=\int_0^1 c_0 \eta (1-\eta) \,{\rm d} \eta \end{equation}

and

(4.13)\begin{equation} f_{{B}} = \frac{1}{2}\int_{0}^{\eta}c_0 \tilde{\eta}^2\,{\rm d}\tilde{\eta} + \left(\frac{\eta^2}{2}-\eta\right)\int_0^{\eta}c_0\tilde{\eta}\,{\rm d}\tilde{\eta} - \frac{\eta^2}{2}\int_{\eta}^{1}c_0(1-\tilde{\eta})\,{\rm d}\tilde{\eta}, \end{equation}

with tildes used to denote dummy integration variables.

4.3. The integro-differential transport equation

As shown by Lawrence et al. (Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), the transport equation that determines the slow spatiotemporal evolution of $c_0(x,\eta,s,\tau )$, given by

(4.14)\begin{equation} \frac{\partial c_0}{\partial \tau}+ u_{{L}}\,\frac{\partial c_0}{\partial x} + \left[\frac{v_{{L}}}{\bar{h}} - \frac{\eta}{\bar{h}}\left(u_{{L}}\,\frac{\partial \bar{h}}{\partial x} + \frac{w_{{L}}}{\ell}\,\frac{\partial \bar{h}}{\partial s}\right)\right]\frac{\partial c_0}{\partial \eta} + \frac{w_{{L}}}{\ell}\,\frac{\partial c_0}{\partial s} = \frac{1}{\alpha^2 \sigma \bar{h}^2}\,\frac{\partial^2 c_0}{\partial \eta^2}, \end{equation}

can be obtained by analysing terms of order $\varepsilon ^2$ in (3.7). The convective transport in the long time scale is found to be driven by the mean Lagrangian velocity

(4.15)\begin{equation} \left.\begin{array}{@{}l} u_{{L}} = u_{{SS}} + u_{{B}} + u_{{SD}}, \\ v_{{L}} = v_{{SS}} + v_{{B}} + v_{{SD}}, \\ w_{{L}} = w_{{SS}} + w_{{B}} + w_{{SD}}, \\ \end{array}\right\} \end{equation}

given by the sum of the cycle-averaged Eulerian velocity $(\langle u_1 \rangle,\langle v_1 \rangle,\langle w_1 \rangle )=(u_{ {SS}}+u_{ { B}},v_{ {SS}}+v_{ {B}},w_{ { SS}}+w_{ {B}})$ and the Stokes drift $(u_{ {SD}},v_{ {SD}},w_{ {SD}})$, the latter being a purely kinematic contribution resulting from the spatial non-uniformity of the pulsatile flow (Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019). The steady-streaming and Stokes-drift contributions to the time-averaged Lagrangian motion, constant and independent of the drug concentration, were identified in our previous analysis (Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), with corresponding expressions given in Appendix A. The slowly varying buoyancy-induced velocity $(u_{ {B}},v_{ {B}},w_{ {B}})$ is a new contribution coupling the bulk motion with the drug concentration. Since the expressions for $(u_{ {B}},v_{ {B}},w_{ {B}})$, given in (4.9)–(4.11), contain spatial integrals of the solute concentration $c_0$, the transport equation (4.14), which is a linear partial differential equation in the buoyancy-free case ${\textit {Ri}}=0$ analysed earlier (Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019), adopts for ${\textit {Ri}} \ne 0$ a nonlinear integro-differential character that complicates the description.

The transport equation (4.14), supplemented with (4.9)–(4.11) for the evaluation of the slowly varying buoyancy-induced velocity $(u_{ {B}},v_{ {B}},w_{ {B}})$ and with the expressions given in Appendix A for the time-independent velocity components $(u_{ {SS}},v_{ {SS}},w_{ {SS}})$ and $(u_{ {SD}},v_{ {SD}},w_{ {SD}})$, can be integrated with boundary conditions $\partial c_0/\partial \eta =0$ at $\eta =(0,1)$ to determine the evolution of the solute. An additional condition must be prescribed at points across the entrance section $x=0$ where there exists inflow (i.e. positive values of $u_{ {L}}$). In the following integrations, it is assumed that the drug concentration of the incoming fluid particles is identically zero, as is consistent with drug delivery in the lumbar region. Bolus injection can be described by using as initial condition the solute distribution $c_0=c_i(x,\eta,s)$ existing at the end of the short injection phase. The description of continuous drug infusion is somewhat more complicated, in that it requires consideration of a localized solute source at the delivery location, a case to be addressed separately in § 7.

Although the reduced Schmidt number $\sigma =S \varepsilon ^2$ can be expected to take order-unity values for the drugs typically used in applications (e.g. $\sigma =0.532\hbox{--}2.128$ when evaluated with $\varepsilon =0.02\hbox{--}0.04$ for methotrexate), it is instructive to investigate simplifications arising for extreme values of this parameter. For example, for $\sigma \gg 1$, the transverse-diffusion term in (4.14) becomes negligible, with the result that the solute particles are transported by the mean Lagrangian velocity while maintaining its initial concentration. Numerical methods specifically tailored to describe Lagrangian particle dispersion can be instrumental to speed up the associated computations (Guan et al. Reference Guan, Jiang, Wang, Zeng, Li and Chen2023). In the opposite limit, $\sigma \ll 1$, diffusion rapidly uniformizes the composition in the transverse direction, so that the concentration becomes independent of $\eta$. The simplified equation applying in this limit can be derived by integrating (4.14) in $\eta$ with boundary conditions $\partial c_0/\partial \eta =0$ at $\eta =(0,1)$, to yield

(4.16)\begin{equation} \frac{\partial c_0}{\partial \tau}+ \bar{u}_{{L}}\,\frac{\partial c_0}{\partial x} + \frac{\bar{w}_{{L}}}{\ell}\,\frac{\partial c_0}{\partial s} =0, \end{equation}

where $\bar {u}_{ {L}}=\int _0^1 u_{ {L}} \,\textrm {d}\eta$ and $\bar {w}_{ {L}}=\int _0^1 w_{ {L}}\, \textrm {d}\eta$ are the width-averaged values of the longitudinal and azimuthal components of the mean Lagrangian velocity. It will be of interest in future work to assess the predictive capability of the above simple equation.

It is worth noting that, unlike direct numerical simulations (DNS) of drug delivery, which need to account for the small cumulative concentration changes that occur over subsequent cardiac cycles, the reduced description (4.14) targets directly the solute evolution in the long time scale $\varepsilon ^{-2} \omega ^{-1}$ that characterizes drug dispersion along the canal. Since the number of cardiac cycles required to achieve significant drug dispersion scales with $\varepsilon ^{-2}$, DNS computations accounting for realistic values of $\varepsilon \sim 0.02\hbox{--}0.04$ must in general consider hundreds of cycles, resulting in computational times that are orders of magnitude larger than those involved in integrating (4.14).

5. Validation of the reduced model

For buoyancy-free systems (i.e. ${\textit {Ri}}=0$), the mean Lagrangian velocity reduces to $(u_{ {L}},v_{ {L}},w_{ {L}}) =(u_{ {SS}}+u_{ {SD}},v_{ {SS}}+v_{ {SD}},w_{ {SS}}+w_{ {SD}})$, independent of the solute concentration, with the result that the associated transport equation (4.14) becomes a linear partial differential equation with time-independent coefficients. The accuracy of the resulting simplified description was tested previously (Gutiérrez-Montes et al. Reference Gutiérrez-Montes, Coenen, Lawrence, Martínez-Bazán, Sánchez and Lasheras2021) by comparing the model predictions with results of DNS computations involving integrations of the complete Navier–Stokes equations. The previous comparisons are extended here to cases with ${\textit {Ri}} \ne 0$, for which (4.14) displays its complicated nonlinear integro-differential character. As in the previous paper, results are given below for two different geometrical configurations with constant perimeter $\ell =1$, namely, a constant-eccentricity annular canal bounded by parallel cylindrical surfaces, yielding a canal width $\bar {h}(s)=1-0.5\cos (2{\rm \pi} s)$, and a variable-eccentricity configuration with canal width $\bar {h}(x,s)=1-0.5\cos (2{\rm \pi} s)\cos (2{\rm \pi} x)$. The latter geometry is selected as a simplified model to mimic changes in the position of the spinal cord relative to the dura mater existing along the human spinal canal, which are depicted in figures 1(b) and 1(c). As one traverses the spine caudally, the spinal cord, which is closer to the posterior side of the canal in the cervical region, moves closer to the anterior side in the thoracic region, eventually returning to the posterior side in the lumbar region. These changes in the spinal canal eccentricity are known to produce changes in the direction of the longitudinal mean Lagrangian velocity (Coenen et al. Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019), leading to the recirculating pattern of bulk CSF flow shown in figure 1(d).

The validation addresses the temporal evolution of the solute following the release of a finite dose, with the initial solute concentration described 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_0$ and having a saturated core flanked by thin layers across which the concentration decays to zero. The values $\delta =0.2$ and $x_0=0.65$ are selected in the sample computations shown below.

The numerical scheme for the integration of (4.14) utilizes a second-order centred finite-difference approximation for the spatial discretization of the viscous terms, and an upwind scheme for the nonlinear terms. A second-order explicit Runge–Kutta scheme is used for time marching, with the integral expressions (4.9)–(4.11) evaluated with a simple trapezoidal rule. A detailed account of the numerical scheme employed in the accompanying DNS computations can be found in Gutiérrez-Montes et al. (Reference Gutiérrez-Montes, Coenen, Lawrence, Martínez-Bazán, Sánchez and Lasheras2021). The DNS computations were performed for a dimensionless stroke length $\varepsilon =0.02$, so that every unit in the long time scale $\tau$ corresponds to $(2 {\rm \pi}\varepsilon )^{-2} \simeq 400$ oscillatory cycles in the DNS computations. The resulting concentration, which includes short-time fluctuations associated with the oscillatory flow, is cycled-averaged to give $\langle c \rangle =\int _t^{t+2{\rm \pi} } c \, \textrm {d} t/(2{\rm \pi} )$, to be compared with the associated model prediction $c_0$.

Results are shown in figures 2 (constant eccentricity) and 3 (variable eccentricity) for a canal with $\alpha =3$, $k=0.5$, $\gamma =1$ and $\sigma =0.4$. To illustrate effects of buoyancy on drug dispersion, in addition to the buoyancy-neutral case ${\textit {Ri}}=0$, the computations consider both a heavy solute with $\rho _d>\rho$ (${\textit {Ri}}=-1$) and a light solute with $\rho _d<\rho$ (${\textit {Ri}}=1$). The figures display three-dimensional views of the entire canal showing isosurfaces of solute concentration $c_0$ for several values of $\tau$. The quantitative comparisons between the model and the DNS include distributions of width-averaged concentrations $\int _0^1 c_0\,\textrm {d}\eta$ and $\int _0^1 \langle c\rangle \,\textrm {d}\eta$ as well as corresponding axial distributions of concentration per unit length of canal, computed according to $C_0 =\int _0^1\bar {h}\int _0^1 c_0\,\textrm {d}\eta \,\textrm {d}s$ and $\langle C\rangle =\int _0^1\bar {h}\int _0^1 \langle c\rangle \,\textrm {d}\eta \,\textrm {d}s$, with the dotted curves representing the initial distribution $C_i =\int _0^1\bar {h}\int _0^1 c_i \,\textrm {d}\eta \,\textrm {d}s$. For reference, the left-hand contour panels showing $\int _0^1 c_0\,\textrm {d}\eta$ include the streamlines corresponding to the width-averaged Lagrangian drift velocity $\bigl (\int _0^1u_{ {L}}\,\textrm {d}\eta,\int _0^1w_{ {L}}\,\textrm {d}\eta \bigr )$, which evolve in time under the action of buoyancy when ${\textit {Ri}} \neq 0$. Figures 2(a) and 3(a) indicate the fraction of the drug bolus that remains in the canal at time $\tau$, as computed with the reduced transport model according to $\chi =\int _0^1 C_0\,\textrm {d}\kern 0.06em x/\int _0^1 C_i \,\textrm {d}\kern 0.06em x$.

Figure 2. The temporal evolution of the solute concentration in a constant-eccentricity canal with $\ell =1$, $\bar {h}(s)=1-0.5\cos (2{\rm \pi} s)$, $\alpha =3$, $k=0.5$, $\gamma =1$ and $\sigma =0.4$ as obtained from the reduced transport equation (4.14) and from DNS computations for three different values of the Richardson number, (b) ${\textit {Ri}}=-1$, (c) ${\textit {Ri}}=0$ and (d) ${\textit {Ri}}=1$, with (a) showing the temporal evolution of the total amount of solute contained in the canal (normalized with its initial value) predicted with the reduced model, as computed from $\chi =\int _0^1 C_0 \,\textrm {d}\kern 0.06em x/\int _0^1 C_i\, \textrm {d}\kern 0.06em x$. The plots include three-dimensional isosurfaces of solute concentration $c_0$, distributions of width-averaged concentrations $\int _0^1 c_0\,\textrm {d}\eta$ and $\int _0^1 \langle c\rangle \,\textrm {d}\eta$, and corresponding axial distributions of concentration per unit length of canal $C_0 =\int _0^1\bar {h}\int _0^1 c_0\, \textrm {d}\eta \,\textrm {d}s$ (solid curves) and $\langle C\rangle =\int _0^1\bar {h}\int _0^1 \langle c\rangle \,\textrm {d}\eta \,\textrm {d}s$ (dashed curves), with the dotted curves representing the initial distribution $C_i =\int _0^1\bar {h}\int _0^1 c_i \,\textrm {d}\eta \,\textrm {d}s$. The streamlines shown in the plots of $\int _0^1 c_0\,\textrm {d}\eta$, corresponding to the width-averaged Lagrangian drift velocity $\bigl (\int _0^1u_{ {L}}\,\textrm {d}\eta,\int _0^1w_{ {L}}\,\textrm {d}\eta \bigr )$, are plotted using constant spacing 0.01 for the associated width-averaged stream function.

Figure 3. Same as figure 2 but for a variable eccentricity canal with $\bar {h}(x,s)=1-0.5\cos (2{\rm \pi} s)\cos (2{\rm \pi} x)$.

Observation of the plots displaying streamlines reveals that the solute moves predominantly following the width-averaged flow, thereby highlighting the important role of the Lagrangian drift in the dispersion of the drug. For a non-buoyant solute in a constant-eccentricity canal, investigated in figure 2(c), the mean Lagrangian flow exhibits a simple circulating pattern, in which the fluid enters along the wide part of the canal ($s=0.5$) and leaves along the narrow part ($s=0$), the motion being slower near the closed end $x=1$. As seen in figures 2(b) and 2(d), the presence of buoyancy alters the flow, with associated streamlines evolving in time as the spatial distribution of the solute changes. Buoyancy promotes rapid ascension of the light solute along the narrow part of the canal, that being the behaviour displayed in figure 2(d). Conversely, heavy solutes tend to sink to the bottom, progression towards the canal entrance being limited to a thin solute filament stretching along the narrow section $s=0$, as seen in figure 2(b). While the overall agreement between the model and the DNS is generally satisfactory, a notable deviation arises at $x=1$ in the heavy-solute results. Here, the model predicts a zero concentration for all times, whereas the DNS yield a concentration that increases over time. These disparities stem from the effect of axial diffusion (not present in the model), which, though negligible elsewhere, becomes significant in this terminal region as the velocity diminishes to zero.

Buoyancy effects are clearly visible in the axial distributions of concentration per unit length of canal $C_0$ and $\langle C\rangle$, and also in the curves representing in figure 2(a) the fraction $\chi$ of the initial bolus that remains inside the canal at time $\tau$. The results indicate that at the longest time computed ($\tau =3$), most of the light solute (91 %) has abandoned the canal, while approximately 82 % of the heavy solute remains inside. This behaviour is consistent with previous clinical observations pertaining to hyperbaric and hypobaric drugs (Mitchell et al. Reference Mitchell, Bowler, Scott and Edström1988; Povey et al. Reference Povey, Jacobsen and Westergaard-Nielsen1989; Richardson et al. Reference Richardson, Thakur, Abramowicz and Wissler1996; Veering et al. Reference Veering, Immink-Speet, Burm, Stienstra and van Kleef2001; Loubert et al. Reference Loubert, Hallworth, Fernando, Columb, Patel, Sarang and Sodhi2011).

For the variable-eccentricity canal shown in figure 3, the streamline patterns of the mean Lagrangian motion feature multiple recirculating regions. The flow direction is reversed between contiguous recirculating cells, as can be inferred from the maps of solute concentration. The solute, carried by the fluid particles, encircles the recirculating regions, thereby hindering the solute progression towards the canal entrance. The plots at $\tau =1$ show most of the light solute accumulating at the interface separating near $x=0.25$ the two top recirculating regions (see figure 3d), while the heavy solute accumulates around $x=0.75$, above the nearly stagnant bottom recirculating region, as shown in figure 3(b). As indicated by comparison of figures 2(a) and 3(a), the rate at which the solute reaches the canal entrance is significantly lower for canals with variable eccentricity, in accordance with previous results (Coenen et al. Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019; Gutiérrez-Montes et al. Reference Gutiérrez-Montes, Coenen, Lawrence, Martínez-Bazán, Sánchez and Lasheras2021).

The agreement between the model and the DNS results is very satisfactory, quantitative departures remaining consistently small regardless of the value of ${\textit {Ri}}$. The degree of agreement is particularly remarkable in connection with the dashed and solid curves representing the longitudinal distribution of the solute at different instants of time. In view of the comparisons shown in figures 2 and 3, it can be concluded that the reduced model provides a sufficiently accurate description for most purposes while requiring computational times that are a fraction of those involved in the DNS computations. For instance, to generate the results corresponding to each value of ${\textit {Ri}}$ in figures 2 and 3, the computations using the reduced model were completed in approximately 10 minutes using a laptop computer, whereas the DNS computations took approximately a week on a 24-core cluster.

6. Dispersion of a drug bolus

The reduced transport equation (4.14) can be used to generate predictions of drug dispersion based on subject-specific canal boundaries and dimensions, with the model parameters determined using magnetic resonance imaging (MRI) measurements, as explained in Coenen et al. (Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019). The sample computations shown below use measurements corresponding to a 25-year old woman (subject 1 in Coenen et al. Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019), with relevant anatomical and Lagrangian-flow details shown in figures 1(b)–1(d). High-resolution images of the entire spine were segmented to extract the three-dimensional position of the pia and dura mater, with the cauda equina (the group of roots branching off at the end of the spinal cord in the lumbar region) represented as an extension of the spinal cord with cross-sectional area tapering down to the end of the spinal canal. The resulting canal anatomy is shown in figure 1(c), with the transverse dimension scaled by a factor 3 to facilitate visualization. A Gaussian filter was used to generate smooth distributions of the perimeter and width of the canal, their mean values $\ell _c=21.8$ mm and $h_c=3.6$ mm employed to scale the geometrical functions $\ell (x)$ and $\bar {h}(x,s)$ used in the model, with the longitudinal distance $x$ being scaled with the total canal length $L=59$ cm. As explained in Coenen et al. (Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019), the compliance of the canal was determined by comparing predictions of oscillatory flow rate with phase-contrast MRI measurements, yielding the function $\gamma '(x)=14.3[0.8+0.3\tanh (4x - 0.2)]$ m MPa$^{-1}$ with mean value $\gamma '_c = 14.107\ \textrm {m}\ \textrm {MPa}^{-1}$. For this subject, the associated values of the Womersley number and elastic wavenumber were found to be $\alpha =h_c/(\nu /\omega )^{1/2}=10.8$ and $k=L\omega /[(h_c/\gamma _c')/\rho ]^{1/2}=0.73$, respectively.

As discussed earlier in connection with figures 2 and 3, the solute moves predominantly following the Lagrangian drift. Before computing drug dispersion, it is therefore of interest to investigate the structure of the mean Lagrangian flow in the absence of buoyancy forces for the anatomically correct canal shown in figure 1(c). To that end, streamlines corresponding to the width-averaged velocity $\bigl (\int _0^1u_{ {L}}\,\textrm {d}\eta,\int _0^1w_{ {L}}\,\textrm {d}\eta \bigr )$ with $(u_{ {L}},w_{ {L}}) =(u_{ {SS}}+u_{ {SD}},w_{ {SS}}+w_{ {SD}})$ are plotted in figure 1(d). The resulting flow pattern comprises three main recirculating regions that occupy approximately the cervical, thoracic and lumbar regions, along with smaller recirculating regions distributed along the posterior midline ($s=0$). The streamlines plotted correspond to evenly spaced values of the associated stream function, so that the physical distance between contiguous streamlines is a measure of the local flow velocity. As is clear from the plot, the fluid is nearly stagnant in the lumbar region, where drug delivery usually takes place, suggesting that neutrally buoyant or heavy drugs will tend to remain near the injection site. The extent to which buoyancy promotes the dispersion of light drugs is to be evaluated in figure 4(c).

Figure 4. Drug dispersion following delivery of a finite dose via the L3/L4 intervertebral space as predicted for $\sigma =1$ and three different values of the Richardson number, (a) ${\textit {Ri}}=-1$, (b) ${\textit {Ri}}=0$ and (c) ${\textit {Ri}}=1$, by integration of the reduced transport equation (4.14) subject to the initial condition (6.1). The plots include distributions of width-averaged concentrations $\int _0^1 c_0\,\textrm {d}\eta$ at $\tau =0.01,0.04,1,3$ along with three-dimensional isosurfaces of solute concentration $c_0$ at intermediate times $\tau =0.02,0.1,2$.

To mimic an intrathecal injection via the L3/L4 posterior intervertebral space, the description of drug dispersion utilizes as initial condition the Gaussian solute distribution

(6.1)\begin{equation} c_i = \exp\left\{-\left[\left(\frac{x-x_0}{\delta_x}\right)^2+\left(\frac{\eta-\eta_0}{\delta_\eta}\right)^2+\left(\frac{s-s_0}{\delta_s}\right)^2\right]\right\}, \end{equation}

with $(x_0,\eta _0,s_0)=(0.8,0.5,0)$ and $(\delta _x,\delta _\eta,\delta _s)=(1/16,500,2/7)$. The reduced Schmidt number is selected to be $\sigma =\varepsilon ^2 S =1$, corresponding to a drug Schmidt number in the range $625 < S < 2500$ for $\varepsilon =0.02\hbox{--}0.04$. Buoyancy effects are investigated for ${\textit {Ri}}=1$ and $-1$, taken as representative of Midazolam and Morphine. Their temporal evolution is compared in figure 4 with results corresponding to a neutrally buoyant drug. To facilitate visualization, besides three-dimensional distributions of drug concentration $c_0$, the figure shows two-dimensional maps of width-averaged concentration $\int _0^1 c_0\,\textrm {d}\eta$ at selected times, with particular attention given to the short time evolution. For the three cases considered, corresponding supplementary movies are available at https://doi.org/10.1017/jfm.2024.297, showing the evolution of the drug up to $\tau =5$.

The plots in figure 4(b) reveal that since the mean Lagrangian motion exhibits low velocities in the lumbar region, in the absence of buoyancy the initial drug evolution is very slow, with changes in the solute concentration distribution remaining virtually inappreciable for $\tau \leqslant 0.1$. For longer times, the drug spreads following the lumbar recirculating vortices, with the result that the drug concentrates in an elongated region about the $s=0$ axis. For the longest time shown in the figure ($\tau =3$), only a small amount of drug has moved into the thoracic region.

Buoyancy fundamentally alters this dispersion pattern, as seen in figures 4(a) and 4(c). For the localized drug distribution considered in the computations, a fast buoyancy-driven vortex is formed upon injection, as revealed by the closely spaced streamlines shown in the two-dimensional plots for $\tau =0.01$ and $0.04$, rapidly spreading the drug around the spinal cord from the initial injection site. The associated recirculatory motion is directed upwards/downwards along the $s=0$ axis for a light/heavy drug, thereby promoting drug dispersion towards the cranial cavity/sacrum region. The progression rate, very rapid for short times, when the buoyancy-induced velocities are larger as a result of the existing high solute concentrations, slows down for longer times, with the heavy drug adopting a stratified distribution that slowly sinks towards the bottom end of the canal, while the light drug continues to evolve upwards, spreading through the thoracic and cervical regions, and eventually reaching the cranial cavity. The behaviour revealed in the figure is therefore consistent with clinical observations regarding intrathecal injections in a seated position (Wildsmith et al. Reference Wildsmith, McClure, Brown and Scott1981; Mitchell et al. Reference Mitchell, Bowler, Scott and Edström1988; Povey et al. Reference Povey, Jacobsen and Westergaard-Nielsen1989; Richardson et al. Reference Richardson, Thakur, Abramowicz and Wissler1996; Veering et al. Reference Veering, Immink-Speet, Burm, Stienstra and van Kleef2001).

7. The description of continuous drug infusion

Medication by ITDD is often released by continuous infusion with use of a percutaneous catheter connected to an external pump or a totally implanted system. The delivery rates are usually small, with maximum values $\dot {Q}\lesssim 1\ \textrm {ml}\ \textrm {h}^{-1}$ (De Andres et al. Reference De Andres, Hayek, Perruchoud, Lawrence, Reina, De Andres-Serrano, Rubio-Haro, Hunt and Yaksh2022). Since drug dispersion is driven by the mean Lagrangian motion, it can be anticipated that the total volume of drug released in times of order of the characteristic bulk flow residence time $\varepsilon ^{-2}\omega ^{-1}$, given by $\dot {Q} \varepsilon ^{-2}\omega ^{-1}$, will be spread over the entire volume of the canal, $L \ell _c h_c \sim 40\hbox{--}60\ \textrm {ml}$, resulting in characteristic drug concentrations of order

(7.1)\begin{equation} c_c=\frac{\dot{Q}\varepsilon^{{-}2}\omega^{{-}1}}{L \ell_c h_c}, \end{equation}

with $c_c \lesssim 0.01$. As a result, in describing continuous drug infusion, it is appropriate to use an order-unity rescaled concentration $\varphi =c/c_c$. Also, since the density differences associated with the presence of the drug can be expected to be of order $c_c(\rho -\rho _d)$, the Richardson number (3.1), which was defined assuming solute concentrations of order unity, must be replaced with

(7.2)\begin{equation} {\textit{Ri}}^* = \frac{g (\rho-\rho_d) c_c/\rho}{\varepsilon^2 \omega^2 L}, \end{equation}

so that the buoyancy acceleration term $-\varepsilon \,{\textit {Ri}} \, c$ in (3.5) becomes $-\varepsilon \,{\textit {Ri}}^*\,\varphi$.

Drug injection will be modelled using a localized volume source. To evaluate the contribution of the source to the mass and momentum balance, we must compare the characteristic value of the velocity induced by the source $\dot {Q}/(\ell _c h_c)$, obtained by dividing the volumetric injection rate $\dot {Q}$ by the characteristic canal cross-section $\ell _c h_c$, with the characteristic bulk flow velocity $\varepsilon ^{2}\omega L$, the ratio of both quantities reducing simply to $[\dot {Q}/(\ell _c h_c)]/(\varepsilon ^{2}\omega L)=c_c \ll 1$, as can be seen from (7.1). Since drug infusion induces negligibly small velocities, the presence of the localized source can be neglected in the first approximation when writing the continuity and momentum balance equations (3.4)–(3.6), but not in the solute conservation equation (3.7), which takes the form

(7.3)\begin{equation} \frac{\partial \varphi}{\partial t}+\varepsilon \left(u\,\frac{\partial \varphi}{\partial x}+v\,\frac{\partial \varphi}{\partial y}+\frac{w}{\ell}\,\frac{\partial \varphi}{\partial s}\right)=\frac{\varepsilon^2}{\alpha^2 \sigma}\,\frac{\partial^2 \varphi}{\partial y^2} + \varepsilon^2 q, \end{equation}

where the dimensionless function $q(x,\eta,s)$ represents the delivery rate per unit volume, scaled with $\dot {Q}/(L \ell _c h_c)$, so that $\int _0^1\ell \int _0^1\bar {h} \int _0^1 q\,\textrm {d}\eta \,\textrm {d}s\,\textrm {d}\kern 0.06em x = 1$. The asymptotic analysis, which parallels that leading to (4.14), provides in this case the reduced transport equation

(7.4)\begin{equation} \frac{\partial \varphi_0}{\partial \tau}+ u_{{L}}\,\frac{\partial \varphi_0}{\partial x} + \left[\frac{v_{{L}}}{\bar{h}} - \frac{\eta}{\bar{h}}\left(u_{{L}}\,\frac{\partial \bar{h}}{\partial x} + \frac{w_{{L}}}{\ell}\,\frac{\partial \bar{h}}{\partial s}\right)\right]\frac{\partial \varphi_0}{\partial \eta} + \frac{w_{{L}}}{\ell}\,\frac{\partial \varphi_0}{\partial s} = \frac{1}{\alpha^2 \sigma \bar{h}^2}\,\frac{\partial^2 \varphi_0}{\partial \eta^2} + q \end{equation}

for the leading-order representation $\varphi _0$ of the reduced solute concentration $\varphi =\varphi _0+\varepsilon \varphi _1 + \cdots$, with the buoyancy-driven component $(u_{ {B}},v_{ {B}},w_{ {B}})$ of the Lagrangian drift velocity $(u_{ {L}},v_{ { L}},w_{ {L}})$ evaluated from (4.9)–(4.11), with ${\textit {Ri}}$ and $c_0$ replaced by ${\textit {Ri}}^*$ and $\varphi _0$.

To represent injection in the posterior intrathecal region through the L3/L4 intervertebral space, the sample computations shown in figure 5 consider a localized source with a normalized Gaussian distribution $q(x,\eta,s)=q_o/(\int _0^1\ell \int _0^1\bar {h} \int _0^1 q_o\, \textrm {d}\eta \,\textrm {d}s\,\textrm {d}\kern 0.06em x)$ centred at $(x_0,\eta _0,s_0)=(0.8,0.5,0)$, where the function $q_o$ is the exponential distribution found on the right-hand side of (6.1) with $(\delta _x,\delta _\eta,\delta _s)=(1/18,1/5,1/13)$. For the three cases considered, corresponding supplementary movies are available. The integrations, initiated with a zero drug concentration everywhere in the canal, describe transient drug infusion for three different reduced Richardson numbers ${\textit {Ri}}^*=c_c\,{\textit {Ri}}$, with the values ${\textit {Ri}}^*=-0.1$ and $0.1$ being comparable to, although somewhat larger than, those expected in connection with the dispersion of Meperidine and Fentanyl (see table 1). As in figure 4, figure 5 shows three-dimensional distributions of drug concentration $\varphi _0$ along with two-dimensional maps of width-averaged concentration $\int _0^1 \varphi _0\,\textrm {d}\eta$. Note that for each plot, the scale of the colour contours has been adjusted to accommodate the increasing concentration, which is found to be significantly larger for non-buoyant drugs.

Figure 5. Drug dispersion corresponding to continuous drug infusion via the L3/L4 intervertebral space as predicted for $\sigma =1$ and three different values of the rescaled Richardson number, (a) ${\textit {Ri}}^*=-0.1$, (b) ${\textit {Ri}}^*=0$ and (c) ${\textit {Ri}}^*=0.1$, by integration of the reduced transport equation (7.4) with a localized solute source centred at $(x_0,\eta _0,s_0)=(0.8,0.5,0)$. The plots include distributions of width-averaged concentrations $\int _0^1 \varphi _0\,\textrm {d}\eta$ at $\tau =0.02,0.1,0.5,2$ along with three-dimensional isosurfaces of solute concentration $\varphi _0$ at intermediate times $\tau =0.05,0.2,1$.

As can be seen in the plots of figure 5(b), the neutrally buoyant drug accumulates near the injection location while spreading longitudinally along the posterior axis $s=0$ at a small rate determined by the existing mean Lagrangian velocity. In contrast, the heavy drug with ${\textit {Ri}}^*=-0.1$, shown in figure 5(a), immediately begins to sink upon injection, driving a recirculatory motion that promotes simultaneous azimuthal spreading. At $\tau =0.2$, the drug has already reached the sacral end of the canal, where it accumulates, forming a stratified distribution that is continuously stirred by the persistent buoyancy-driven recirculatory flow. Up to the longest time considered ($\tau =2$), the heavy drug is confined to the lumbar region, with the result that the mean Lagrangian motion remains virtually unperturbed in the thoracic and cervical regions. On the other hand, infusion of light drugs, considered in figure 5(c), leads to the development of a plume. The light fluid rises until it reaches the boundary separating the lumbar and thoracic recirculating regions, forming a front at $x\simeq 0.6$, corresponding approximately to the T11/T12 intervertebral space. At that level, the drug spreads azimuthally to reach the anterior side, where it continues to flow upwards into the thoracic region, thereby resuming its progression towards the cranial cavity.

In analysing the transient results of figure 5, one should bear in mind that while the present computation assumes impermeable surfaces, leading to continuous drug accumulation, in ITDD processes drug uptake by the spinal nerve as well as through the dura membrane would eventually balance the infusion rate, leading to a steady drug distribution along the spine. For heavy drugs, the results shown in figure 5(a) suggest that the combined effects of buoyancy forces and drug uptake may limit drug dispersion to the lumbar and sacral regions. On the other hand, the results in figure 5(c) indicate that the ability of light drugs to reach the cranial cavity will depend on the competition of buoyancy-enhanced drug dispersion and drug absorption, whose quantification necessitates an extended reduced model accounting for pharmacokinetic effects.