2 CSF motion in the spinal canal

2.1 Geometrical considerations

In analyzing the motion of CSF, the spinal SAS will be modelled as a thin annular gap of length $L$ whose cross-section varies slowly between the base of the cranial cavity and the lumbar region. The orthogonal curvilinear coordinates $(x,y,s)$ shown in the lower inset of figure 1 will be employed in the analysis, with $x$ denoting the distance from the canal entrance measured along the spinal cord, $y$ being the normal distance to the pia mater (the inner surface), and $s$ being the distance measured from the symmetry plane around the surface of the pia mater normalized with the spinal–cord perimeter $\ell (x)$ . The values of the coordinates are in the ranges $0\leqslant x\leqslant L$ , $0\leqslant y\leqslant h(x,s,t)$ and $0\leqslant s\leqslant 1$ , where $h(x,s,t)$ is the canal width. Characteristic values of $h$ and $\ell$ are $h_{c}\sim 0.1~\text{cm}$ and $\ell _{c}\sim 2~\text{cm}$ , respectively, while $L\sim 60{-}80~\text{cm}$ , resulting in the inequalities

(2.1) $$\begin{eqnarray}L\gg \ell _{c}\gg h_{c},\end{eqnarray}$$

which are to be used in simplifying the description, as shown below.

The conservation equation in terms of the curvilinear coordinates $(x,y,s)$ can be written following the general expressions given in Batchelor (Reference Batchelor2000). In particular, when the condition (2.1) is accounted for, the continuity equation takes the simplified form

(2.2) $$\begin{eqnarray}\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\ell u)+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}s}=0,\end{eqnarray}$$

where $u$ , $v$ and $w$ are the velocity components in the $x$ , $y$ and $s$ directions, respectively. A straightforward order-of-magnitude balance of (2.2) provides

(2.3) $$\begin{eqnarray}\frac{u_{c}}{L}\sim \frac{v_{c}}{h_{c}}\sim \frac{w_{c}}{\ell _{c}}\end{eqnarray}$$

relating the characteristic values $u_{c}$ , $v_{c}$ and $w_{c}$ of the three velocity components. These scalings correspond to streamlines aligned in the $x$ direction (i.e. $u\sim (L/\ell _{c})w\gg w$ ) with an accompanying slow transverse motion occurring predominantly in the azimuthal direction with $w\sim (\ell _{c}/h_{c})v\gg v$ , with $v$ being the smallest velocity component.

In the resulting slender flow, the characteristic values of the spatial pressure differences needed to establish the motion in the $x$ , $y$ and $s$ directions are given by $\unicode[STIX]{x1D6FF}_{x}p=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714}u_{c}L$ , $\unicode[STIX]{x1D6FF}_{y}p=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714}v_{c}h_{c}$ and $\unicode[STIX]{x1D6FF}_{s}p=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714}w_{c}\ell _{c}$ , respectively, as follows from a balance between the pressure gradient and the local acceleration. These values are very different in magnitude, as can be seen by using (2.3), yielding

(2.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}_{x}p}{L^{2}}\sim \frac{\unicode[STIX]{x1D6FF}_{y}p}{h_{c}^{2}}\sim \frac{\unicode[STIX]{x1D6FF}_{s}p}{\ell _{c}^{2}}.\end{eqnarray}$$

The scalings (2.4) can be taken into account when describing the variations of the pressure $p$ from the entrance value $p_{c}$ . Neglecting the small contribution to the pressure fluctuation associated with breathing, the value of $p_{c}$ can be written as

(2.5) $$\begin{eqnarray}p_{c}(t)=p_{o}+(\unicode[STIX]{x0394}p)_{c}\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D714}t),\end{eqnarray}$$

where $p_{o}$ is the time-averaged intracranial pressure, $(\unicode[STIX]{x0394}p)_{c}$ is the amplitude of the cranial pressure pulse, and $\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D714}t)$ is a dimensionless periodic function with angular frequency $\unicode[STIX]{x1D714}\simeq 2\unicode[STIX]{x03C0}$ Hz, associated with the cardiac cycle. According to (2.4), with small errors of order $(h_{c}/\ell _{c})^{2}$ one may neglect the pressure variations in the $y$ direction when writing

(2.6) $$\begin{eqnarray}p-p_{c}(t)=p^{\prime }(x,t)+\tilde{p}(x,s,t)\end{eqnarray}$$

for the pressure variation along the canal, with $p^{\prime }(x,t)$ denoting the value of the pressure difference $p-p_{c}(t)$ along the curve $s=y=0$ . The term $\tilde{p}\sim (\ell _{c}/L)^{2}p^{\prime }\ll p^{\prime }$ , measuring the small pressure differences around the canal, must be included to describe the azimuthal motion.

2.2 Constitutive equation and elastic parameters

The dura membrane deforms in response to the local pressure variations to accommodate the inflow and outflow of CSF into the canal, so that the canal width $h(x,s,t)$ is a function of time that must be determined as part of the solution. The deformations $h^{\prime }=h-h_{o}$ are to be defined relative to the unperturbed width distribution $h_{o}(x,s)$ describing the geometry of the canal when the pressure everywhere is uniform and equal to $p=p_{o}$ , with $A_{o}(x)=\ell (x)\,\int _{0}^{1}h_{o}(x,s)\,\text{d}s$ representing the unperturbed value of the cross-sectional area. Since the tidal volume is small compared with the total volume of CSF, the associated changes in the shape of the canal must satisfy

(2.7) $$\begin{eqnarray}(h-h_{o})/h_{o}\sim \unicode[STIX]{x0394}V/V\ll 1.\end{eqnarray}$$

The description of the coupling between the fluid motion and the deformation of the dura membrane in general amounts to a very complicated fluid–structure interaction problem. A simplified solution is pursued below on the basis of the presumed linear relation $h^{\prime }\propto (p-p_{o})$ . The description accounts also for the fact that the pressure at a given canal section is nearly uniform, with small relative variations, represented in (2.6) by the term $\tilde{p}$ , that are of order $(\ell _{c}/L)^{2}$ . If these are neglected, then one may write $p-p_{o}\simeq (\unicode[STIX]{x0394}p)_{c}\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D714}t)+p^{\prime }(x,t)$ , which, together with the additional assumption that the elastic properties of the dura membrane are independent of $s$ , yields the simplified constitutive equation

(2.8) $$\begin{eqnarray}\frac{p-p_{o}}{E_{e}}=\frac{h^{\prime }(x,t)}{A_{o}/\ell }.\end{eqnarray}$$

In writing (2.8) we have chosen to scale the deformation with the average canal width $A_{o}/\ell \sim h_{c}$ and include for generality an effective elastic modulus $E_{e}(x)$ , the latter embodying the overall effect of different microanatomical features such as the distribution of veins and fatty epidural tissue in the dura membrane. Variations of the elastic properties of the dura membrane at a given section, not considered in the present development, could be incorporated in the model by allowing the effective elastic modulus $E_{e}$ to be a function of $s$ , resulting in the local deformation $h^{\prime }$ being also a function of $s$ , according to (2.8).

According to (2.7) and (2.8) the value of $E_{e}$ must be much larger than $p-p_{o}$ , so that the resulting deformations remain small according to $(p-p_{o})/E_{e}\sim \unicode[STIX]{x0394}V/V\ll 1$ . The ratio of the characteristic value of the local pressure variation $p-p_{o}\sim (\unicode[STIX]{x0394}p)_{c}$ to the characteristic value $E_{c}$ of $E_{e}$ defines the small parameter

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D700}=\frac{(\unicode[STIX]{x0394}p)_{c}}{E_{c}}\sim \frac{\unicode[STIX]{x0394}V}{V}\ll 1,\end{eqnarray}$$

to be used in the following description as a measure of the compliance.

The effective elastic modulus $E_{e}$ determines the wave speed $(E_{e}/\unicode[STIX]{x1D70C})^{1/2}$ of the resulting pulsatile motion along the canal. Noninvasive measurements (Kalata et al. Reference Kalata, Martin, Oshinski, Jerosch-Herold, Royston and Loth2009) using MRI have shown this wave speed to be on the order of a few metres per second, giving associated wavelengths $(E_{c}/\unicode[STIX]{x1D70C})^{1/2}/\unicode[STIX]{x1D714}$ that are comparable to the canal length $L$ . Correspondingly, the ratio

(2.10) $$\begin{eqnarray}k=\frac{L}{(E_{c}/\unicode[STIX]{x1D70C})^{1/2}/\unicode[STIX]{x1D714}}\sim 1,\end{eqnarray}$$

defines a dimensionless elastic wavenumber of order unity, another important parameter in the description below.

2.3 Simplified conservation equations

We give below the equations that determine the velocity, whose $x$ , $y$ and $s$ components, $u(x,y,s,t)$ , $v(x,y,s,t)$ and $w(x,y,s,t)$ , must satisfy the non-slip boundary conditions $u=v=w=0$ at $y=0$ and $u=v-\unicode[STIX]{x2202}h^{\prime }/\unicode[STIX]{x2202}t=w=0$ at $y=h(x,s,t)$ . The problem also involves the displacement $h^{\prime }(x,t)=h-h_{o}(x,s)$ , and the pressure field, described in terms of the functions $p^{\prime }(x,t)$ and $\tilde{p}(x,s,t)$ .

The computation of $u$ and $w$ requires consideration of the $x$ and $s$ components of the momentum equation, which can be written for the curvilinear coordinates defined above by neglecting terms of order $(\ell _{c}/L)^{2}$ or $(h_{c}/\ell _{c})^{2}$ (or smaller) to give

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+v\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}+\frac{w}{\ell }\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}s}=-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}y^{2}}, & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}t}+\frac{u}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\ell w)+v\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}y}+\frac{w}{\ell }\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}s}=-\frac{1}{\unicode[STIX]{x1D70C}\ell }\frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}s}+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}w}{\unicode[STIX]{x2202}y^{2}}, & \displaystyle\end{eqnarray}$$

whereas the $y$ component of the momentum equation, which would be needed to compute the small spatial pressure differences $\unicode[STIX]{x1D6FF}_{y}p$ , can be omitted in the slender-flow approximation at the order considered here. At the same level of approximation, the viscous terms have been simplified in (2.11) and (2.12) by discarding the terms involving $\unicode[STIX]{x2202}^{2}/\unicode[STIX]{x2202}x^{2}$ and $\unicode[STIX]{x2202}^{2}/\unicode[STIX]{x2202}s^{2}$ , which are smaller than the terms that have been retained by a factor $(h_{c}/L)^{2}$ and $(h_{c}/\ell _{c})^{2}$ , respectively.

The transverse velocity component $v$ can be evaluated in terms of $u$ and $w$ from

(2.13) $$\begin{eqnarray}v+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\ell \int _{0}^{y}u\,\text{d}{\tilde{y}}\right)+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\int _{0}^{y}w\,\text{d}{\tilde{y}}\right)=0,\end{eqnarray}$$

involving the dummy integration variable ${\tilde{y}}$ . The above equation is obtained by integrating the continuity equation (2.2) in the $y$ direction with the condition $u=v=w=0$ at $y=0$ . Evaluating the above equation at $y=h$ gives

(2.14) $$\begin{eqnarray}\ell \frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\ell \int _{0}^{h}u\,\text{d}y\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\int _{0}^{h}w\,\text{d}y\right)=0,\end{eqnarray}$$

which is the corresponding Reynolds lubrication equation for the problem.

Further integrations of the continuity equation also become useful in the development, which assumes that the spinal canal is symmetric, such that $h_{o}(x,s)=h_{o}(x,1-s)$ , resulting in the conditions $\int _{0}^{h}w\,\text{d}y=0$ at $s=0$ and at $s=1/2$ . Multiplying (2.14) by $\text{d}s$ and integrating between $0$ and $s$ provides

(2.15) $$\begin{eqnarray}\ell \frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}s+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\ell \int _{0}^{s}\left(\int _{0}^{h}u\,\text{d}y\right)\text{d}s\right]+\int _{0}^{h}w\,\text{d}y=0\end{eqnarray}$$

after using the symmetry condition $\int _{0}^{h}w\,\text{d}y=0$ at $s=0$ , whereas the condition $\int _{0}^{h}w\,\text{d}y=0$ at $s=1/2$ (or $s=1$ ) yields

(2.16) $$\begin{eqnarray}\ell \frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\ell \int _{0}^{1}\left(\int _{0}^{h}u\,\text{d}y\right)\text{d}s\right].\end{eqnarray}$$

The displacement $h^{\prime }$ is related to the local pressure $p-p_{o}\simeq (\unicode[STIX]{x0394}p)_{c}\unicode[STIX]{x1D6F1}+p^{\prime }$ by the constitutive equation (2.8), which can be differentiated with respect to time to give

(2.17) $$\begin{eqnarray}\ell \frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}=\frac{A_{o}}{E_{e}}\left((\unicode[STIX]{x0394}p)_{c}\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}t}+\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}\right).\end{eqnarray}$$

The additional equation

(2.18) $$\begin{eqnarray}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\ell \int _{0}^{1}\left(\int _{0}^{h}u\,\text{d}y\right)\text{d}s\right]=\frac{A_{o}}{E_{e}}\left((\unicode[STIX]{x0394}p)_{c}\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}t}+\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}\right)\end{eqnarray}$$

obtained by combining (2.16) and (2.17) will be found below to be useful in computing $p^{\prime }$ as a function of the streamwise velocity $u$ .

2.4 Dimensionless formulation

An order-of-magnitude analysis provides the characteristic scales for the problem. The development begins by using (2.18) with $(\unicode[STIX]{x0394}p)_{c}/E_{e}\sim \unicode[STIX]{x1D700}$ to provide $u_{c}=\unicode[STIX]{x1D700}L\unicode[STIX]{x1D714}$ as an estimate for the characteristic value of $u$ , on the order of a few centimetres per second, as can be seen by using $\unicode[STIX]{x1D700}\sim \unicode[STIX]{x0394}V/V\sim 0.02$ , $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}$ and $L\sim 60{-}80~\text{cm}$ . The characteristic values $w_{c}=\unicode[STIX]{x1D700}\ell _{c}\unicode[STIX]{x1D714}$ and $v_{c}=\unicode[STIX]{x1D700}h_{c}\unicode[STIX]{x1D714}$ follow from (2.3). The spatial pressure changes needed to move the fluid can be estimated from the balance between the local acceleration and the pressure gradient in (2.11) and (2.12) to give $\unicode[STIX]{x1D6FF}_{x}p=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D700}(\unicode[STIX]{x1D714}L)^{2}$ and $\unicode[STIX]{x1D6FF}_{s}p=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D700}(\unicode[STIX]{x1D714}\ell _{c})^{2}$ . Additional order-of-magnitude balances in these two equations reveal that the convective acceleration is smaller than the local acceleration by a factor $\unicode[STIX]{x1D700}$ , whereas the comparison between the local acceleration and the viscous force yields

(2.19) $$\begin{eqnarray}\frac{O(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}t)}{O(\unicode[STIX]{x1D708}\unicode[STIX]{x2202}^{2}u/\unicode[STIX]{x2202}y^{2})}\sim \frac{O(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}t)}{O(\unicode[STIX]{x1D708}\unicode[STIX]{x2202}^{2}w/\unicode[STIX]{x2202}y^{2})}\sim \unicode[STIX]{x1D6FC}^{2}=\frac{\unicode[STIX]{x1D714}h_{c}^{2}}{\unicode[STIX]{x1D708}},\end{eqnarray}$$

where

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=h_{c}/(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714})^{1/2}\end{eqnarray}$$

is the Womersley number of the flow (the reciprocal of the square root of the relevant Stokes number $\unicode[STIX]{x1D708}/(h_{c}^{2}\unicode[STIX]{x1D714})$ ), typically of order unity, as can be seen by using $h_{c}\sim 10^{-3}~\text{m}$ , $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}~\text{ s}^{-1}$ and $\unicode[STIX]{x1D708}\sim 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ .

The relative scalings discussed above become clear when writing the above problem in dimensionless form by introduction of the dimensionless variables $u^{\ast }=u/u_{c}$ , $v^{\ast }=v/v_{c}$ , $w^{\ast }=w/w_{c}$ , $p^{\prime \ast }=p^{\prime }/\unicode[STIX]{x1D6FF}_{x}p$ , $\tilde{p}^{\ast }=\tilde{p}/\unicode[STIX]{x1D6FF}_{s}p$ , $t^{\ast }=\unicode[STIX]{x1D714}t$ , $x^{\ast }=x/L$ , $\ell ^{\ast }=\ell /\ell _{c}$ , $h_{o}^{\ast }=h_{o}/h_{c}$ , $h^{\ast }=h/h_{c}$ , $h^{\prime \ast }=(h-h_{o})/(\unicode[STIX]{x1D700}h_{c})$ , $A_{o}^{\ast }=A_{o}/(\ell _{c}h_{c})$ and $E_{e}^{\ast }=E_{e}/E_{c}$ . In addition, to facilitate the description of the temporal changes of the canal width it is convenient to scale the distance $y$ with the local width $h(x,s,t)$ through the normalized coordinate $\unicode[STIX]{x1D702}=y/h(x,s,t)$ . In what follows, the asterisks used to denote dimensionless variables will be dropped to simplify the notation.

Since we have a total of six unknowns (i.e. $u$ , $v$ , $w$ , $p^{\prime }$ , $\tilde{p}$ and $h^{\prime }$ , with $h=h_{o}+\unicode[STIX]{x1D700}h^{\prime }$ ) the description requires in principle six equations, which are given below in dimensionless form. We begin by writing (2.11) and (2.12) as

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D700}F_{x}=-\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}+\frac{1}{h^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}, & \displaystyle\end{eqnarray}$$

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D700}F_{s}=-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}s}+\frac{1}{h^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}w}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}, & \displaystyle\end{eqnarray}$$

where

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle F_{x}=\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\ell u^{2})+\frac{1}{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(uv)+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}(uw)-\frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x1D702}}{h}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x1D702}}{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(u^{2})-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x1D702}}{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(uw) & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and

(2.24) $$\begin{eqnarray}\displaystyle F_{s} & = & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(uw)+2\frac{uw}{\ell }\frac{\unicode[STIX]{x2202}\ell }{\unicode[STIX]{x2202}x}+\frac{1}{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(vw)+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}(w^{2})-\frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x1D702}}{h}\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x1D702}}{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(uw)-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x1D702}}{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(w^{2})\end{eqnarray}$$

represent the convective terms, including the apparent convection associated with $\unicode[STIX]{x2202}h^{\prime }/\unicode[STIX]{x2202}t$ . The streamwise pressure variation $p^{\prime }$ is related to the volume flux by (2.18), which can be written as

(2.25) $$\begin{eqnarray}\frac{A_{o}}{E_{e}}\left(\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}t}+k^{2}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}\right)=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\ell \int _{0}^{1}h\left(\int _{0}^{1}u\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s\right],\end{eqnarray}$$

whereas the deformation $h^{\prime }$ satisfies

(2.26) $$\begin{eqnarray}\ell \frac{\unicode[STIX]{x2202}h^{\prime }}{\unicode[STIX]{x2202}t}=\frac{A_{o}}{E_{e}}\left(\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}t}+k^{2}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}\right),\end{eqnarray}$$

as obtained from (2.17). The transverse velocity component $v$ can be evaluated from (2.13) in terms of $u$ and $w$ to yield

(2.27) $$\begin{eqnarray}v=-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\ell h\int _{0}^{\unicode[STIX]{x1D702}}u\,\text{d}\unicode[STIX]{x1D702}\right)-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(h\int _{0}^{\unicode[STIX]{x1D702}}w\,\text{d}\unicode[STIX]{x1D702}\right)+\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}\unicode[STIX]{x1D702}u+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}s}\unicode[STIX]{x1D702}w.\end{eqnarray}$$

Finally, the integral continuity balance

(2.28) $$\begin{eqnarray}h\int _{0}^{1}w\,\text{d}\unicode[STIX]{x1D702}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[s\,\ell \int _{0}^{1}h\left(\int _{0}^{1}u\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s-\ell \int _{0}^{s}h\left(\int _{0}^{1}u\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s\right],\end{eqnarray}$$

obtained by combining (2.15) and (2.16), will be useful in computing the azimuthal variation of the pressure $\tilde{p}$ .

The velocity components $u$ and $w$ must satisfy the non-slip boundary conditions $u=w=0$ at $\unicode[STIX]{x1D702}=0,1$ . At the canal entrance the pressure is $p=p_{c}(t)$ , resulting in the condition

(2.29) $$\begin{eqnarray}p^{\prime }=0\quad \text{at}~x=0.\end{eqnarray}$$

Also, the streamwise volume flux must vanish at the closed end of the canal, so that

(2.30) $$\begin{eqnarray}\int _{0}^{1}h\left(\int _{0}^{1}u\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s=0\quad \text{at}~x=1.\end{eqnarray}$$

Besides the Womersley number $\unicode[STIX]{x1D6FC}$ defined in (2.20) and the dimensionless elastic wavenumber $k$ defined in (2.10), both of order unity, the above problem depends on the elasticity parameter $\unicode[STIX]{x1D700}\ll 1$ , identified earlier in (2.9), with $\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D700}\ll 1$ representing the relevant Reynolds number for the flow. As can be seen in (2.21) and (2.22), at leading order in the limit $\unicode[STIX]{x1D700}\ll 1$ the flow in the thin canal corresponds to an oscillatory lubrication problem, where the motion is determined by a balance between the local acceleration, the pressure gradient, and the viscous forces associated with the transverse velocity derivatives.

3 Leading-order oscillatory flow

The solution can be obtained using perturbation analysis by introducing expansions for the different flow variables of the form $u=u_{0}+\unicode[STIX]{x1D700}u_{1}+\cdots \,,v=v_{0}+\unicode[STIX]{x1D700}v_{1}+\cdots \,,w=w_{0}+\unicode[STIX]{x1D700}w_{1}+\cdots \,,p^{\prime }=p_{0}^{\prime }+\unicode[STIX]{x1D700}p_{1}^{\prime }+\cdots \,$ and $\tilde{p}=\tilde{p}_{0}+\unicode[STIX]{x1D700}\tilde{p}_{1}+\cdots \,,$ together with the expansion $h^{\prime }(x,t)=(h-h_{o})/\unicode[STIX]{x1D700}=h_{0}^{\prime }+\unicode[STIX]{x1D700}h_{1}^{\prime }+\cdots \,$ for the deformation of the canal. It will be seen that the leading-order solution is purely oscillatory, whereas the corrections include a steady component corresponding to the steady bulk-flow motion, with characteristic velocities of order $\unicode[STIX]{x1D700}u_{c}=\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D714}L$ (i.e. a few centimetres per minute, in agreement with the observed transport rates in the spinal canal).

At leading order, (2.21), (2.22), and (2.25)–(2.28) reduce to the linear equations

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}p_{0}^{\prime }}{\unicode[STIX]{x2202}x}+\frac{1}{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}w_{0}}{\unicode[STIX]{x2202}t}=-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\tilde{p}_{0}}{\unicode[STIX]{x2202}s}+\frac{1}{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}w_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{A_{o}}{E_{e}}\left(\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}t}+k^{2}\frac{\unicode[STIX]{x2202}p_{0}^{\prime }}{\unicode[STIX]{x2202}t}\right)=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\ell \int _{0}^{1}h_{o}\left(\int _{0}^{1}u_{0}\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s\right], & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \ell \frac{\unicode[STIX]{x2202}h_{0}^{\prime }}{\unicode[STIX]{x2202}t}=\frac{A_{o}}{E_{e}}\left(\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}t}+k^{2}\frac{\unicode[STIX]{x2202}p_{0}^{\prime }}{\unicode[STIX]{x2202}t}\right), & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle v_{0}=-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\ell h_{o}\int _{0}^{\unicode[STIX]{x1D702}}u_{0}\,\text{d}\unicode[STIX]{x1D702}\right)-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(h_{o}\int _{0}^{\unicode[STIX]{x1D702}}w_{0}\,\text{d}\unicode[STIX]{x1D702}\right)+\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}x}\unicode[STIX]{x1D702}u_{0}+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}s}\unicode[STIX]{x1D702}w_{0}, & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle h_{o}\int _{0}^{1}w_{0}\,\text{d}\unicode[STIX]{x1D702}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[s\,\ell \int _{0}^{1}h_{o}\left(\int _{0}^{1}u_{0}\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s-\ell \int _{0}^{s}h_{o}\left(\int _{0}^{1}u_{0}\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s\right]. & \displaystyle\end{eqnarray}$$

The solution depends on the shape of the canal through the known functions $h_{o}(x,s)$ and $\ell (x)$ , on its elastic properties through the function $E_{e}(x)$ , and on the temporal variation of the intracranial pressure, defined by the periodic function $\unicode[STIX]{x1D6F1}(t)$ .

An analytic solution can be obtained for the case of a harmonic intracranial pressure $\unicode[STIX]{x1D6F1}(t)=\cos t$ by introducing the complex functions $U(x,\unicode[STIX]{x1D702},s)$ , $V(x,\unicode[STIX]{x1D702},s)$ , $W(x,\unicode[STIX]{x1D702},s)$ , $P^{\prime }(x)$ , $\tilde{P}(x,s)$ and $H^{\prime }(x)$ defined such that

(3.7) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\displaystyle u_{0}=\text{Re}(\text{i}\text{e}^{\text{i}t}U),\quad v_{0}=\text{Re}(\text{i}\text{e}^{\text{i}t}V),\quad w_{0}=\text{Re}(\text{i}\text{e}^{\text{i}t}W),\\ \displaystyle p_{0}^{\prime }=\text{Re}(\text{e}^{\text{i}t}P^{\prime }),\quad \tilde{p}_{0}=\text{Re}(\text{e}^{\text{i}t}\tilde{P}),\quad \text{and}\quad h_{0}^{\prime }=\text{Re}(\text{e}^{\text{i}t}H^{\prime }).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The functions

(3.8a,b ) $$\begin{eqnarray}U=\frac{\text{d}P^{\prime }}{\text{d}x}G\quad \text{and}\quad W=\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}s}G,\end{eqnarray}$$

with

(3.9) $$\begin{eqnarray}G(x,\unicode[STIX]{x1D702},s)=1-\frac{\cosh \left[{\displaystyle \frac{\unicode[STIX]{x1D6FC}h_{o}}{2}}{\displaystyle \frac{1+\text{i}}{\sqrt{2}}}\left(2\unicode[STIX]{x1D702}-1\right)\right]}{\cosh \left[{\displaystyle \frac{\unicode[STIX]{x1D6FC}h_{o}}{2}}{\displaystyle \frac{1+\text{i}}{\sqrt{2}}}\right]}\end{eqnarray}$$

are determined by integration of

(3.10) $$\begin{eqnarray}\frac{\text{i}}{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+U=\frac{\text{d}P^{\prime }}{\text{d}x};\quad U=0~\text{at}~\unicode[STIX]{x1D702}=(0,1)\end{eqnarray}$$

and

(3.11) $$\begin{eqnarray}\frac{\text{i}}{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}W}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+W=\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}s};\quad W=0~\text{at}~\unicode[STIX]{x1D702}=(0,1),\end{eqnarray}$$

derived from (3.1) and (3.2), respectively. On the other hand, the function $G$ can be integrated to give

(3.12) $$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x1D702}}G\,\text{d}\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}-\frac{1-\text{i}}{\sqrt{2}\unicode[STIX]{x1D6FC}h_{o}}\frac{\displaystyle \sinh \left[\frac{\unicode[STIX]{x1D6FC}h_{o}}{2}\frac{1+\text{i}}{\sqrt{2}}\left(2\unicode[STIX]{x1D702}-1\right)\right]+\sinh \left[\frac{\unicode[STIX]{x1D6FC}h_{o}}{2}\frac{1+\text{i}}{\sqrt{2}}\right]}{\displaystyle \cosh \left[\frac{\unicode[STIX]{x1D6FC}h_{o}}{2}\frac{1+\text{i}}{\sqrt{2}}\right]}\end{eqnarray}$$

which can be used in (3.5) to provide

(3.13) $$\begin{eqnarray}\displaystyle V & = & \displaystyle -\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\ell h_{o}\frac{\text{d}P^{\prime }}{\text{d}x}\int _{0}^{\unicode[STIX]{x1D702}}G\,\text{d}\unicode[STIX]{x1D702}\right)-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{h_{o}}{\ell }\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}s}\int _{0}^{\unicode[STIX]{x1D702}}G\,\text{d}\unicode[STIX]{x1D702}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}x}\frac{\text{d}P^{\prime }}{\text{d}x}\unicode[STIX]{x1D702}G+\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}s}\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}s}\unicode[STIX]{x1D702}G\end{eqnarray}$$

for the velocity component normal to the surface.

The expressions given in (3.8) and (3.13) for the velocity components involve the pressure gradients $\text{d}P^{\prime }/\text{d}x$ and $\unicode[STIX]{x2202}\tilde{P}/\unicode[STIX]{x2202}s$ , to be obtained for given values of $h_{o}(x,s)$ and $\ell (x)$ with use made of (3.3) and (3.6). In the computation it is convenient to introduce the functions

(3.14) $$\begin{eqnarray}q(x,s)=h_{o}\int _{0}^{1}G\,\text{d}\unicode[STIX]{x1D702}=h_{o}-\frac{\sqrt{2}(1-\text{i})}{\unicode[STIX]{x1D6FC}}\tanh \left(\frac{\unicode[STIX]{x1D6FC}h_{o}}{2}\frac{1+\text{i}}{\sqrt{2}}\right)\end{eqnarray}$$

and

(3.15) $$\begin{eqnarray}Q(x)=\ell \int _{0}^{1}q\,\text{d}s=\ell \int _{0}^{1}\left[h_{o}-\frac{\sqrt{2}(1-\text{i})}{\unicode[STIX]{x1D6FC}}\tanh \left(\frac{\unicode[STIX]{x1D6FC}h_{o}}{2}\frac{1+\text{i}}{\sqrt{2}}\right)\right]\text{d}s,\end{eqnarray}$$

which define the volume fluxes

(3.16a,b ) $$\begin{eqnarray}h_{o}\int _{0}^{1}u_{0}\,\text{d}\unicode[STIX]{x1D702}=\text{Re}\!\left(\text{i}\text{e}^{\text{i}t}\frac{\text{d}P^{\prime }}{\text{d}x}q\right),\quad h_{o}\int _{0}^{1}w_{0}\,\text{d}\unicode[STIX]{x1D702}=\text{Re}\!\left(\frac{\text{i}\text{e}^{\text{i}t}}{\ell }\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}s}q\right),\end{eqnarray}$$

and

(3.17) $$\begin{eqnarray}\ell \int _{0}^{1}h_{o}\left(\int _{0}^{1}u_{0}\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s=\text{Re}\!\left(\text{i}\text{e}^{\text{i}t}\frac{\text{d}P^{\prime }}{\text{d}x}Q\right)\end{eqnarray}$$

appearing in (3.3) and (3.6). In particular, equation (3.3) can be used, together with the boundary conditions (2.29) and (2.30), to write the boundary-value problem

(3.18) $$\begin{eqnarray}\frac{E_{e}(x)}{A_{o}(x)}\frac{\text{d}}{\text{d}x}\left[Q(x)\frac{\text{d}P^{\prime }}{\text{d}x}\right]+k^{2}P^{\prime }+1=0;\quad P^{\prime }=0\text{ at }x=0\text{ and }\frac{\text{d}P^{\prime }}{\text{d}x}=0\text{ at }x=1\end{eqnarray}$$

for the function $P^{\prime }(x)$ , which can be used in (3.4) to give

(3.19) $$\begin{eqnarray}H^{\prime }=\frac{A_{o}}{\ell E_{e}}(k^{2}P^{\prime }+1)\end{eqnarray}$$

for the deformation of the outer surface and in (3.6) to yield

(3.20) $$\begin{eqnarray}\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}s}=\frac{1}{q}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\frac{\text{d}P^{\prime }}{\text{d}x}\left(sQ-\ell \int _{0}^{s}q\,\text{d}s\right)\right]\end{eqnarray}$$

for the azimuthal pressure gradient.

In general, numerical integration is needed to solve the boundary-value problem (3.18). The closed-form solution

(3.21) $$\begin{eqnarray}P^{\prime }=\frac{1}{k^{2}}\left\{\frac{\cos [k(1-x)/\sqrt{Q}]}{\cos (k/\sqrt{Q})}-1\right\}\end{eqnarray}$$

is obtained for uniform elastic modulus $E_{e}=1$ when the cross-section has constant shape (i.e. $h_{o}=h_{o}(s)$ and $\ell =A_{o}=1$ ), so that $Q=\text{const}$ .

For a given canal geometry, defined by $h_{o}(x,s)$ , $\ell (x)$ and $A_{o}(x)=\ell (x)\int _{0}^{1}h_{o}\,\text{d}s$ , with given elastic properties, defined by $E_{e}(x)$ , the determination of the oscillatory motion begins by employing (3.15) to compute $Q(x)$ . The resulting function is to be used in (3.18) when computing $P^{\prime }(x)$ , which in turn provides through (3.19) the deformation $H^{\prime }$ . The associated gradient $\text{d}P^{\prime }/\text{d}x$ can be used in (3.8) to compute $U$ and in (3.20) to evaluate  $\unicode[STIX]{x2202}\tilde{P}/\unicode[STIX]{x2202}s$ , which allows us to determine $W$ and $V$ from (3.8) and (3.13). Finally, the complex functions $U$ , $V$ , $W$ , $P^{\prime }$ and $H^{\prime }$ can be used in (3.7) to provide $u_{0}$ , $v_{0}$ , $w_{0}$ , $p_{0}^{\prime }$ and $h_{0}^{\prime }$ .

4 Steady-streaming motion

Because of nonlinear interactions, associated with the convective terms and with the temporal and spatial variations of the canal width, the first-order corrections to the flow contain a steady component, additional to the oscillatory component. The computation of this steady-streaming flow begins by collecting terms of order $\unicode[STIX]{x1D700}$ in (2.21) and (2.22). Taking the time average $\langle \cdot \rangle =(1/2\unicode[STIX]{x03C0})\int _{0}^{2\unicode[STIX]{x03C0}}\cdot \,\,\text{d}t$ of the resulting equations yields

(4.1) $$\begin{eqnarray}\displaystyle {\mathcal{F}}_{x} & = & \displaystyle -\frac{\unicode[STIX]{x2202}\langle p_{1}^{\prime }\rangle }{\unicode[STIX]{x2202}x}+\frac{1}{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}\langle u_{1}\rangle }{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}},\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle {\mathcal{F}}_{s} & = & \displaystyle -\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\langle \tilde{p}_{1}\rangle }{\unicode[STIX]{x2202}s}+\frac{1}{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}\langle w_{1}\rangle }{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}},\end{eqnarray}$$

where

(4.3) $$\begin{eqnarray}\displaystyle {\mathcal{F}}_{x} & = & \displaystyle \frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\ell \langle u_{0}^{2}\rangle )+\frac{1}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\langle u_{0}v_{0}\rangle +\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\langle u_{0}w_{0}\rangle \nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D702}}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left\langle \frac{\unicode[STIX]{x2202}h_{0}^{\prime }}{\unicode[STIX]{x2202}t}u_{0}\right\rangle -\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x1D702}}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\langle u_{0}^{2}\rangle -\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x1D702}}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\langle u_{0}w_{0}\rangle +\frac{2}{h_{o}^{3}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}\langle h_{0}^{\prime }u_{0}\rangle \qquad\end{eqnarray}$$

and

(4.4) $$\begin{eqnarray}\displaystyle {\mathcal{F}}_{s} & = & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\langle u_{0}w_{0}\rangle +2\frac{\langle u_{0}w_{0}\rangle }{\ell }\frac{\unicode[STIX]{x2202}\ell }{\unicode[STIX]{x2202}x}+\frac{1}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\langle v_{0}w_{0}\rangle +\frac{1}{\ell }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\langle w_{0}^{2}\rangle -\frac{\unicode[STIX]{x1D702}}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left\langle \frac{\unicode[STIX]{x2202}h_{0}^{\prime }}{\unicode[STIX]{x2202}t}w_{0}\right\rangle \nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x1D702}}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\langle u_{0}w_{0}\rangle -\frac{1}{\ell }\frac{\unicode[STIX]{x2202}h_{o}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x1D702}}{h_{o}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\langle w_{0}^{2}\rangle +\frac{2}{h_{o}^{3}\unicode[STIX]{x1D6FC}^{2}}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}\langle h_{0}^{\prime }w_{0}\rangle\end{eqnarray}$$

can be evaluated in terms of the leading-order solution. The steady–streaming velocities must satisfy $\langle u_{1}\rangle =\langle w_{1}\rangle =0$ at $\unicode[STIX]{x1D702}=0,1$ .

The functions ${\mathcal{F}}_{x}$ and ${\mathcal{F}}_{s}$ drive the steady-streaming motion, in that, if they were zero, the solution to (4.1) and (4.2) would reduce to $\langle u_{1}\rangle =\langle w_{1}\rangle =0$ . Besides the well-known contributions arising from the time average of the nonlinear convective acceleration, present in any steady-streaming phenomenon, additional terms appear in (4.3) and (4.4) due to the deformation of the canal (i.e. the terms involving $h_{0}^{\prime }$ and $\unicode[STIX]{x2202}h_{0}^{\prime }/\unicode[STIX]{x2202}t$ ). Similar contributions have been identified earlier in studies of steady-streaming in elastic tubes (Wang & Tarbell Reference Wang and Tarbell1992). In principle, all terms in (4.3) and (4.4) can contribute significantly to the motion depending on the flow conditions and on the specific geometrical features of the canal. For the specific model problem considered below in § 5 the canal width $h_{o}$ and its perimeter $\ell$ are selected to be independent of $x$ , so that the terms involving $\unicode[STIX]{x2202}h_{o}/\unicode[STIX]{x2202}x$ and $\unicode[STIX]{x2202}\ell /\unicode[STIX]{x2202}x$ are identically zero. All other terms were found to be significant, with their relative importance varying along the canal.

Integrating (4.1) and (4.2) subject to $\langle u_{1}\rangle =\langle w_{1}\rangle =0$ at $\unicode[STIX]{x1D702}=0,1$ yields

(4.5) $$\begin{eqnarray}\frac{\langle u_{1}\rangle }{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}=-\frac{\text{d}\langle p_{1}^{\prime }\rangle }{\text{d}x}\frac{(1-\unicode[STIX]{x1D702})\unicode[STIX]{x1D702}}{2}+\unicode[STIX]{x1D702}\int _{0}^{\unicode[STIX]{x1D702}}{\mathcal{F}}_{x}\,\text{d}\bar{\unicode[STIX]{x1D702}}-\int _{0}^{\unicode[STIX]{x1D702}}{\mathcal{F}}_{x}\bar{\unicode[STIX]{x1D702}}\,\text{d}\bar{\unicode[STIX]{x1D702}}-\unicode[STIX]{x1D702}\int _{0}^{1}{\mathcal{F}}_{x}(1-\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}\end{eqnarray}$$

and

(4.6) $$\begin{eqnarray}\frac{\langle w_{1}\rangle }{h_{o}^{2}\unicode[STIX]{x1D6FC}^{2}}=-\frac{1}{\ell }\frac{\unicode[STIX]{x2202}\langle \tilde{p}_{1}\rangle }{\unicode[STIX]{x2202}s}\frac{(1-\unicode[STIX]{x1D702})\unicode[STIX]{x1D702}}{2}+\unicode[STIX]{x1D702}\int _{0}^{\unicode[STIX]{x1D702}}{\mathcal{F}}_{s}\,\text{d}\bar{\unicode[STIX]{x1D702}}-\int _{0}^{\unicode[STIX]{x1D702}}{\mathcal{F}}_{s}\bar{\unicode[STIX]{x1D702}}\,\text{d}\bar{\unicode[STIX]{x1D702}}-\unicode[STIX]{x1D702}\int _{0}^{1}{\mathcal{F}}_{s}(1-\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}.\end{eqnarray}$$

The average streamwise pressure gradient $\text{d}\langle p_{1}^{\prime }\rangle /\text{d}x$ , which completes the determination of $\langle u_{1}\rangle$ , is a function of $x$ that can be computed by imposing the condition

(4.7) $$\begin{eqnarray}\int _{0}^{1}h_{o}\left(\int _{0}^{1}\langle u_{1}\rangle \,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s+\int _{0}^{1}\int _{0}^{1}\langle h_{0}^{\prime }u_{0}\rangle \,\text{d}\unicode[STIX]{x1D702}\,\text{d}s=0\end{eqnarray}$$

obtained at this order from the time average of (2.25) with use made of (2.30). Similarly, the azimuthal pressure gradient $\unicode[STIX]{x2202}\langle \tilde{p}_{1}\rangle /\unicode[STIX]{x2202}s$ appearing in (4.6) can be determined from the condition

(4.8) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\ell \int _{0}^{s}\left(\int _{0}^{1}(h_{o}\langle u_{1}\rangle +\langle h_{0}^{\prime }u_{0}\rangle )\,\text{d}\unicode[STIX]{x1D702}\right)\text{d}s\right]+h_{o}\int _{0}^{1}\langle w_{1}\rangle \,\text{d}\unicode[STIX]{x1D702}+\int _{0}^{1}\langle h_{0}^{\prime }w_{0}\rangle \,\text{d}\unicode[STIX]{x1D702}=0\end{eqnarray}$$

obtained at this order from the time average of (2.28).

5 Selected results for a model problem

A simple geometrical model will be used to illustrate some of the salient features of the flow. Specifically, we consider the flow in the annular canal bounded between two cylindrical surfaces of radii $R$ and $R+h_{c}$ whose axes are separated by a distance $\unicode[STIX]{x1D6FD}h_{c}$ , with $h_{c}\ll R$ and $0\leqslant \unicode[STIX]{x1D6FD}<1$ (see the schematic view shown in figure 4(a), to be discussed later). Using $h_{c}$ and $\ell _{c}=2\unicode[STIX]{x03C0}R$ as characteristic scales leads to the dimensionless geometrical functions $h_{o}=1-\unicode[STIX]{x1D6FD}\cos (2\unicode[STIX]{x03C0}s)$ , $\ell =1$ and $A_{o}=1$ , resulting in a constant value of the reduced volume flux $Q$ , as can be computed from (3.15). Furthermore, a uniform elastic modulus $E_{e}=1$ is assumed, so that the complex function $P^{\prime }(x)$ for the axial pressure variation reduces to (3.21). Similar simple models involving coaxial flexible tubes have been considered in analyses of pressure-wave propagation in the spinal canal (Berkouk, Carpenter & Lucey Reference Berkouk, Carpenter and Lucey2003; Carpenter, Berkouk & Lucey Reference Carpenter, Berkouk and Lucey2003).

Besides the geometry, which is defined by the eccentricity parameter $\unicode[STIX]{x1D6FD}$ for the simplified problem, the solution depends on the values of $k$ and $\unicode[STIX]{x1D6FC}$ , both assumed to be of order unity, as corresponds to the scales characterizing the flow in the spinal canal. As can be seen from the definitions given in (2.10) and (2.20), the specific values $k=0.5$ and $\unicode[STIX]{x1D6FC}=3$ employed in the computations used for figures 3 through 5 correspond to a canal with width $h_{c}\sim 1$ mm and elastic wave speed $(E_{c}/\unicode[STIX]{x1D70C})^{1/2}\sim 10~\text{m}~\text{s}^{-1}$ .

Figure 2. The value of $|Q\,\text{d}P^{\prime }/\text{d}x|$ evaluated with use made of (3.21) and (5.1) for $h_{o}=1-\unicode[STIX]{x1D6FD}\cos (2\unicode[STIX]{x03C0}s)$ , $\ell =1$ and $\unicode[STIX]{x1D6FD}=0.5$ . Panel (a) shows the variation with $x$ , whereas (b) shows the parametric dependence of the entrance value $|Q\,\text{d}P^{\prime }/\text{d}x|_{x=0}$ .

The leading-order oscillatory flow is investigated in figures 2 and 3. In particular, figure 2 represents the variation of

(5.1) $$\begin{eqnarray}\left|Q\frac{\text{d}P^{\prime }}{\text{d}x}\right|=\left|\frac{\sqrt{Q}}{k}\frac{\sin [k(1-x)/\sqrt{Q}]}{\cos (k/\sqrt{Q})}\right|,\end{eqnarray}$$

which is the amplitude of the tidal volume flux across the canal, according to (3.17). The plot includes the variation of this quantity with $x$ as well as the dependences of the entrance value $|Q\,\text{d}P^{\prime }/\text{d}x|_{x=0}=|(\sqrt{Q}/k)\tan (k/\sqrt{Q})|$ on $\unicode[STIX]{x1D6FC}$ and $k$ for $\unicode[STIX]{x1D6FD}=0.5$ . Additional computations for other values of $\unicode[STIX]{x1D6FD}$ (not shown here) revealed that the eccentricity has a negligible effect on the tidal volume flux, so that the curves obtained with $\unicode[STIX]{x1D6FD}\neq 0.5$ are almost identical to those plotted in the figure.

Figure 3. The distribution of $|u_{0}|=|U|$ for $k=0.5$ and $\unicode[STIX]{x1D6FC}=3$ as obtained at different sections $x$ for varying eccentricities.

As can be seen in figure 2(a), the tidal volume flux decreases along the canal to vanish at the closed end. As expected, larger values of $\unicode[STIX]{x1D6FC}$ , associated with smaller viscous forces, result in increased flow rates. The dependence of $|Q\,\text{d}P^{\prime }/\text{d}x|_{x=0}$ on $k$ is shown in figure 2(b). For small $k$ , i.e. wavelengths that are much larger than the canal length, the pressure along the canal is almost uniform, and the dura membrane deforms uniformly in response to the pressure oscillations in the cranial cavity, following (2.26), resulting in an entrance flow rate that is independent of the flow conditions, so that the value of $|Q\,\text{d}P^{\prime }/\text{d}x|_{x=0}$ for $k\ll 1$ is identical for all $\unicode[STIX]{x1D6FC}$ . In the opposite limit $k\gg 1$ the wavelength is much shorter than the canal length, so that at a given instant of time we find regions of positive and negative deformation of the dura membrane along the canal, which tend to have a cancelling effect on the entrance flow rate, reflected in the decrease of $|Q\,\text{d}P^{\prime }/\text{d}x|_{x=0}$ for $k\gg 1$ . An interesting result from the model is that $|Q\,\text{d}P^{\prime }/\text{d}x|_{x=0}$ reaches a maximum at an intermediate value of the wavenumber $k$ (i.e. $k\simeq 1$ for $\unicode[STIX]{x1D6FC}=3$ ), associated with elastic wave speeds $(E_{c}/\unicode[STIX]{x1D70C})^{1/2}$ on the order of $\unicode[STIX]{x1D714}L$ .

Although variations of the eccentricity do not significantly alter the tidal volume flux, the velocity fields associated with different values of $\unicode[STIX]{x1D6FD}$ are very different, as seen in the comparisons in figure 3, which exhibit the distributions of streamwise velocity amplitude $|u_{0}|=|U|$ at different cross-sections and for varying eccentricities, including the concentric case, $\unicode[STIX]{x1D6FD}=0$ . Although the analytical solution given above applies only to configurations with $h_{c}\ll R$ , the width of the annular cross-section in the figure is arbitrarily enlarged to facilitate visualization. As expected, the decelerating effect of viscous forces becomes more (less) effective in regions of smaller (larger) canal width, leading to the nonuniform velocity distributions depicted in the figure.

Figure 4. Schematic view of the model geometry considered here (a), distribution of $\langle u_{1}\rangle$ corresponding to $k=0.5$ , $\unicode[STIX]{x1D6FC}=3$ and $\unicode[STIX]{x1D6FD}=0.5$ (b), and streamwise variation of the maximum (downward) and minimum (upward) velocities for different values of $\unicode[STIX]{x1D6FD}$ (c).

Steady streaming velocities $\langle u_{1}\rangle$ , evaluated for $k=0.5$ and $\unicode[STIX]{x1D6FC}=3$ from (4.5) with $\text{d}\langle p_{1}^{\prime }\rangle /\text{d}x$ computed using (4.7), are plotted in figures 4 and 5. The circular plot in figure 4(b) shows the distribution of $\langle u_{1}\rangle$ for $\unicode[STIX]{x1D6FD}=0.5$ and $x=0.25$ . Transverse profiles of $\langle u_{1}\rangle$ across the narrowest ( $s=0$ ) and widest ( $s=1/2$ ) sections are plotted below at different heights for this same eccentricity. As can be seen, the resulting velocities are negative (upwards) in the narrow part of the canal and predominantly positive (downwards) in the wider regions. The associated flow pattern is schematically represented by the arrows in figure 4(a).

The circular plots in figure 5 represent distributions of $\langle u_{1}\rangle$ across the canal section at different heights and eccentricities, as in figure 3. The computations for $\unicode[STIX]{x1D6FD}\neq 0$ reveal that the maximum velocity $\langle u_{1}\rangle _{max}$ always occurs in the symmetry plane at the centre point of the widest section. By way of contrast, negative (upward) velocities are found in the narrow part of the canal, in agreement with the transverse profiles shown in figure 4. The minimum velocity $\langle u_{1}\rangle _{min}$ occurs at the narrowest section $s=0$ for sufficiently small values of $\unicode[STIX]{x1D6FD}>0$ but moves away from the symmetry plane for large eccentricities, as can be clearly observed in the results for $\unicode[STIX]{x1D6FD}=0.75$ . Also of interest is the fact that, as can be inferred from the plots corresponding to $\unicode[STIX]{x1D6FD}=0$ , the apparent volume flux $\int _{0}^{1}h_{o}(\int _{0}^{1}\langle u_{1}\rangle \,\text{d}\unicode[STIX]{x1D702})\,\text{d}s$ associated with the steady-steaming flow is nonzero, that being a direct consequence of the deformation of the canal wall, which enters in the integral continuity balance (4.7).

The distributions of $\langle u_{1}\rangle _{max}$ and $\langle u_{1}\rangle _{min}$ along the canal are plotted on figure 4(b) for different values of $\unicode[STIX]{x1D6FD}$ . The bulk motion, characterized by these two bounding values, is found to be extremely weak for concentric cylinders, but becomes much more pronounced for eccentric canals, for which we find values of $\langle u_{1}\rangle$ of order unity, corresponding to dimensional steady-streaming velocities of order $\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D714}L$ . It is also of interest that the configuration with intermediate eccentricity ( $\unicode[STIX]{x1D6FD}=0.5$ ) exhibits velocities that are larger than those found with $\unicode[STIX]{x1D6FD}=0.25$ and $\unicode[STIX]{x1D6FD}=0.75$ . This non-monotonic dependence, associated with the nonlinear interactions of the steady-streaming flow, is an indication of the importance of the geometrical effects. The results suggest that accurate quantifications of CSF flow in the spinal canal should consider the specific geometric features of the SAS, described in our model by the functions $h_{o}(x,s)$ and $\ell (x)$ .

Figure 5. The distribution of $\langle u_{1}\rangle$ for $k=0.5$ and $\unicode[STIX]{x1D6FC}=3$ as obtained at different sections $x$ for varying eccentricities.

For completeness, the distributions of azimuthal velocity $\langle w_{1}\rangle$ across the canal section are shown in figure 6 for the case $k=0.5$ and $\unicode[STIX]{x1D6FC}=3$ corresponding to the plots in figure 5. As can be seen, the distributions exhibit the expected symmetry, with the azimuthal velocity vanishing at the symmetry plane, defined by the sections $s=0$ and $s=0.5$ . The circumferential motion, identically zero for concentric configurations, becomes more pronounced for increasing $\unicode[STIX]{x1D6FD}$ , with the configuration with $\unicode[STIX]{x1D6FD}=0.5$ exhibiting the largest values of $\langle w_{1}\rangle$ . It is also of interest that the azimuthal motion is faster near the canal end, as needed to accommodate the more rapid deceleration of the streamwise motion occurring there, which is apparent in figure 4(c). The location of the peak values of $\langle w_{1}\rangle$ , found near $s\approx 0.25$ at the entrance, progressively move towards the narrowest section ( $s=0$ ) for increasing values of $x$ , an effect that is more noticeable as $\unicode[STIX]{x1D6FD}$ increases.

Figure 6. The distribution of $\langle w_{1}\rangle$ for $k=0.5$ and $\unicode[STIX]{x1D6FC}=3$ as obtained at different sections $x$ for varying eccentricities.

6 Experiments of steady streaming in slender elastic annular tubes

We have conducted a series of in-vitro experiments in order to qualitatively validate the predictions of our asymptotic analysis. The experiments involve the simplified geometry depicted in figure 4(a), that is, a straight, slender, elastic annular canal of length $L=25~\text{cm}$ bounded between cylindrical tubes of radii $R=7~\text{mm}$ and $R+h_{c}=10.5~\text{mm}$ . Axisymmetric cases ( $\unicode[STIX]{x1D6FD}=0$ ) with coaxial tubes were investigated along with a case with eccentricity $\unicode[STIX]{x1D6FD}=0.1$ . The annular tube is connected at one end to a rigid container where the pressure $p_{c}$ is varied periodically through a peristaltic pump with a frequency of 1 Hz to produce an oscillatory flow at the entrance of the channel of $\pm 2~\text{cm}~\text{s}^{-1}$ . The values of $h_{c}$ and $\unicode[STIX]{x1D714}$ yield $\unicode[STIX]{x1D6FC}=h_{c}/(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714})^{1/2}\simeq 8.7$ for the associated Womersley number. The ratio of the displaced volume to that contained in the unperturbed annular tube was selected to be $\unicode[STIX]{x0394}V/V=1/50$ , consistent with the values in the clinical observations. The results shown below correspond to an elastic outer tube, representing the deformable dura membrane, with an inner rigid rod, representing the spinal cord. Additional experiments using an inner elastic tube surrounded by a rigid outer cylinder and two elastic tubes were found to give similar flow features.

The flow velocity was measured in a coaxial plane at various downstream locations along the tube by particle image velocimetry (PIV) using neutrally buoyant hollow, spherical, silver-coated micro particles with a diameter of $0.5~\unicode[STIX]{x03BC}\text{m}$ . Figure 7 shows the time variation of the oscillatory velocity, averaged over the cross-section of the annular gap at various downstream locations along the tube for the case of a concentric annular tube ( $\unicode[STIX]{x1D6FD}=0$ ). The flow corresponds to a pure harmonic motion of biphasic periodic tides of ebb-and-flow with frequency $\unicode[STIX]{x1D714}$ , with amplitudes monotonically decreasing as one moves from the open end to the closed end of the annular tube, consistent with the predictions of the zeroth-order solution of the asymptotic analysis shown in figures 2 and 3, and as also observed in many clinical MRI measurements.

Figure 7. Time variation of the oscillatory axial velocity averaged along the cross-section at various downstream locations along the tube for the case of an axisymmetric annular canal ( $\unicode[STIX]{x1D6FD}=0$ ) with flexible outer wall and rigid inner surface.

Figure 8. (a) Time evolution of the fluorescent marker initially injected in the annular gap at the midpoint between the entrance and closed end of the tube. The open entrance to the channel is located on the left-hand side. (b) Schematic representation of the evolution of the fluorescent marker injected at the midpoint along the tube. (c) Propagation velocity of the dye front towards the open entrance of the tube. Note the monotonic increase of the streaming velocity $\langle u_{1}\rangle$ approaching the open entrance. The measurements correspond to an eccentric case with $\unicode[STIX]{x1D6FD}=0.1$

The PIV measurements were sufficiently accurate to capture the decay in the amplitude of the leading-order oscillatory flow as shown in figure 7. Unfortunately, the polymer tubes that are necessary for simulation of the compliance of the canal are not optically clean and give rise to undesirable reflections and scattering effects that affected the accuracy of the PIV velocity measurements near the walls of the narrow annular gap. Because the resulting measurement error is much greater than the streaming velocities, integrating these measurements over 20–30 min (1200–1800 cycles) results in very noisy measurements. Thus, to quantify the velocity of the streaming motion we opted to use the far more reliable method of injecting a bolus of neutrally buoyant fluorescent dye in the annular gap of the canal at the midpoint between its entrance and closed end, and accurately measure over 20–30 min the velocity of the upward and downward propagation of the fronts as shown in the schematic panels of figure 8. We selected this method not only for its accuracy compared to the PIV measurements but also because it fully resembles the technique used in all the radiological measurements done in the clinical setting, where a bolus of marker is injected in the lumbar region, making our results easier to understand by the medical community. The tube was oriented horizontally with the inner tube offset slightly to introduce an eccentricity ( $\unicode[STIX]{x1D6FD}=0.1$ ), as shown in figure 8(b).

The annular tube was illuminated with an ultraviolet light and imaged using a CCD camera with a resolution of $2000\times 2000~\text{pixels}$ . Fluorescence images were recorded at the same oscillation phase every 30 s over 20 min to track the time evolution of the injected patch. Figure 8 shows the images at 150 s intervals for the first 15 min. These measurements clearly demonstrate that, in addition to the periodic harmonic oscillatory motion described in figure 7, there is a bulk motion that gradually transports the dye upwards (to the left in this experiment) into the tank and downwards (right) towards the closed end. Measurements of the speed at which these two fronts propagate allowed us to obtain order-of-magnitude estimates for the local velocity of the mean streaming motion (as an example, figure 8(c) represents the upward velocity of the front recorded in one of the experiments). Consistent with the analysis in § 5, we found that the speed of the dye front propagating towards the tank gradually increases as the open end is approached, while at the same time the speed of the opposite front moving to the right decreases on approaching the closed end of the tube, in agreement with the results in figures 4 and 5. It is important to emphasize that, as predicted by the analysis, these streaming velocities are a factor $\unicode[STIX]{x1D700}$ smaller (i.e. about two orders of magnitude) than the amplitude of the velocity fluctuations measured for the harmonic base flow described above. Our results are consistent with many radiological measurements; see, for example, the measurements of Greitz & Hannerz (Reference Greitz and Hannerz1996), who showed in their radionuclide cisternography (figure 4 in that paper) that tracers injected in the lumbar region reached the thoracic area much faster than they moved downwards into the sacral end.

We have further corroborated that the observed bulk transport is not due to any gravitational or thermal effects by repeating the above experiments without the imposed pressure pulsations in the tank. For this case, the dye diffused only several millimetres left and right over 10 h.