1 Introduction
The behaviour of microbial biofluids is important to a broad range of diverse fields, from medicine and biotechnology to biogeochemistry, aquatic ecology and climate science. When suspended cells self-propel, the contrast with passive particle hydrodynamics is particularly dramatic. For example, dilute suspensions of swimming bacteria, algae and ciliates display self-concentration and hydrodynamic instabilities (Pedley & Kessler Reference Pedley and Kessler1992). The dilute and concentrated suspension rheology of swimmers is also of much interest: for example, swimming algae and bacteria can enhance or reduce effective suspension viscosity, respectively (Marchetti et al. Reference Marchetti, Joanny, Liverpool, Prost, Rao and Simha2013). Many motile microorganisms respond to their environment and migrate towards nutrients (chemotaxis) or light (phototaxis). Intriguingly they also respond to flow gradients. In this work we focus on gyrotaxis, the passive orientational bias on swimming microalgae due to a combination of viscous and gravitational torques, but the results may be generalised to other taxes or combinations thereof. In downwelling pipe flow, gyrotactic swimmers self-focus in structures known as plumes (Kessler Reference Kessler1985). Though the discovery of gyrotaxis in algae is now over 30 years old, many open questions remain. Here, we shall address the dispersion of gyrotactic swimmers, and how well experiments can be described by existing theoretical models of gyrotaxis.
The classic papers by Taylor (Reference Taylor1953) and Aris (Reference Aris1956) considered the dispersion of passive solutes in laminar pipe flow, and showed that the cross-sectional averaged axial solute dispersion can be mapped to an effective diffusion, with diffusivity
$D_{e}=D_{m}(1+Pe^{2}/48)$
. Taylor tested this prediction experimentally himself (Taylor Reference Taylor1953), as well as its extension to turbulent pipe flow (Taylor Reference Taylor1954). Recently, we revisited the Taylor–Aris analytical theory for laminar flow to describe the dispersion of swimmers in confined flow (Bees & Croze Reference Bees and Croze2010) and used it to predict the dispersion of gyrotactic swimmers in laminar and turbulent flows (Bearon, Bees & Croze Reference Bearon, Bees and Croze2012; Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013). The Bees & Croze (Reference Bees and Croze2010) theory is general and applicable to any microswimmer. It requires as input the mean response to flow of swimmers from stochastic microscopic models. We have applied it to gyrotactic microalgae, where stochastic models of self-propelled spheroids in flows have been developed (Pedley & Kessler Reference Pedley and Kessler1990, Reference Pedley and Kessler1992; Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003; Bearon et al.
Reference Bearon, Bees and Croze2012), comparing predictions from analytical theory with individual-based numerical simulations (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013) and numerical solutions of the model advection–diffusion equation (Bearon et al.
Reference Bearon, Bees and Croze2012).
The present work on dispersion was inspired by the notion that biotechnologically useful microalgal species, such as the
$\unicode[STIX]{x1D6FD}$
-carotene-producing biflagellated alga Dunaliella salina, are cultured in pipe flow within industrial tubular photobioreactors (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013; Bees & Croze Reference Bees and Croze2014). As we shall discuss below, gyrotaxis at the individual cell level dramatically modifies the population-scale behaviour of suspensions of biflagellates such as D. salina, and this has important consequences for the engineering design of photobioreactors, for both laminar and turbulent flow regimes (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013). Evidence is also mounting that gyrotactic behaviour is relevant to environmental flows. Durham, Kessler & Stocker (Reference Durham, Kessler and Stocker2009) proposed gyrotactic trapping to be one of the mechanisms leading to the formation of oceanic thin layers (Durham & Stocker Reference Durham and Stocker2012). Theory and simulation also predict similar effects in turbulent flows, where gyrotactic swimmers ‘unmix’ in downwelling regions of turbulent flows (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013; Durham et al.
Reference Durham, Climent, Barry, Lillo, Boffetta, Cencini and Stocker2013), or even due to acceleration in strong turbulence (De Lillo et al.
Reference De Lillo, Cencini, Durham, Barry, Stocker, Climent and Boffetta2014). Furthermore, the study of bioconvection (Pedley & Kessler Reference Pedley and Kessler1992; Hill & Pedley Reference Hill and Pedley2005), an instability driven by upswimming and gyrotaxis, has recently seen some interesting developments. Weak shear flow has been shown to distort but not destroy bioconvection patterns (Croze, Ashraf & Bees Reference Croze, Ashraf and Bees2010; Hwang & Pedley Reference Hwang and Pedley2014a
). Light from below can frustrate the patterns (Willams & Bees Reference Willams and Bees2011a
; Williams & Bees Reference Williams and Bees2011b
), while in microchannels light from the side gives rise to horizontal accumulation (Garcia, Rafaï& Peyla Reference Garcia, Rafaï and Peyla2013). The cell excess density (the difference between cell and surrounding medium density) can give rise to instabilities in gyrotactic plumes known as blips (Kessler Reference Kessler1986; Pedley & Kessler Reference Pedley and Kessler1992; Dennisenko & Lukaschuk Reference Dennisenko and Lukaschuk2007). The existence of blips as a function of imposed flow and swimmer concentration can be predicted theoretically using linear stability analysis (Hwang & Pedley Reference Hwang and Pedley2014b
). Even in the absence of a gravitational torque, algae can display peculiar behaviour, such as limit cycle oscillations in pipe flow (Zöttl & Stark Reference Zöttl and Stark2012). Finally, while it is well known that algae swim with a helical trajectory, only recently have the mechanics been addressed (Bearon Reference Bearon2013), qualitatively explaining resonant alignment in oscillatory shear flow (Hope et al.
Reference Hope, Croze, Poon, Bees and Haw2016).
Quantitative comparison between mathematical theory and experiment is essential in fluid mechanics. Taylor excelled and delighted in this comparison (Croze & Peaudecerf Reference Croze and Peaudecerf2016). Without it, theory has no constraint and experiments no mechanism. Furthermore, a quantitatively tested theory is of greater engineering use. In this paper, we set out to test the theory of gyrotactic swimmer dispersion experimentally for vertical pipe flow. In § 2 we summarise the theory of axial cell dispersion, which links the statistical individual-level response to flow, predicted from two distinct microscopic stochastic models, with the population-level transport measures: the effective axial drift and diffusivity. In the process of experiment–theory iteration, we identified that cell concentration affected dispersion. Therefore we explore how algal cell negative buoyancy perturbs the flow field from simple Poiseuille flow, which further modifies the radial distribution of cells, altering axial dispersion. Dispersion with negative buoyancy had previously only been considered for a linearised response to flow (Bees & Croze Reference Bees and Croze2010). Here, we employ expressions for the nonlinear response to flow predicted by microscopic stochastic descriptions. In § 3 we describe experimental methods for cell culture, fluorescent staining and imaging, as well as data analysis and comparison with theoretical predictions. We also explain the challenges due to having to maximise fluorescent image contrast while avoiding blip instabilities. In spite of these issues, we establish a robust experimental protocol to study gyrotactic dispersion and obtain preliminary results allowing one to compare predicted and measured axial drift for various average cell concentrations. Finally, in § 4 we propose further experiments to quantify dispersion and discuss the broader relevance of our results for growing algae in photobioreactors.
2 Theory
In this paper, we test theoretical predictions of biased microswimmer dispersion with experimental data. To enable a direct comparison, we derive new results from the general axial dispersion framework (Bees & Croze Reference Bees and Croze2010; Bearon et al. Reference Bearon, Bees and Croze2012; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013) to include both the effect of swimmer negative buoyancy and the complete nonlinear response of swimmers to shear. We summarise the theory below, referring the reader to the above papers for further mathematical details.
Consider a suspension of gyrotactic swimming algae (such as D. salina) subject to imposed flow in a pipe. Gyrotactic swimmers actively cross streamlines and focus in downwelling flows, which is at the root of their peculiar non-Taylor dispersion. Following Pedley & Kessler (Reference Pedley and Kessler1992) and the assumptions discussed therein, an incompressible, dilute suspension of gyrotactic swimming cells in an axial flow in a tube can be described by the equations
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=\frac{1}{\unicode[STIX]{x1D70C}}[-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+nv_{c}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\boldsymbol{g}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[n(\boldsymbol{u}+\boldsymbol{q}V_{s})-\frac{V_{s}^{2}}{d_{r}}\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D735}n\right], & \displaystyle\end{eqnarray}$$
where (2.1) is the Navier–Stokes equation for the pipe flow
$\boldsymbol{u}$
, incorporating a driving pressure gradient term (with pressure
$p$
), a Newtonian viscous stress (medium viscosity
$\unicode[STIX]{x1D707}$
, approximately equal to that of water) and a negative buoyancy term (cell volume
$v_{c}$
, density excess
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$
, gravitational acceleration
$\boldsymbol{g}$
, cell concentration
$n$
). Negative buoyancy causes cell accumulations to sink, which can result in instabilities such as blips (Hwang & Pedley Reference Hwang and Pedley2014b
); see §§ 3.2 and 3 below. Equation (2.2) expresses conservation of cells, with advection contributions from flow and swimming (velocity
$\boldsymbol{q}V_{s}$
, mean orientation
$\boldsymbol{q}$
, and speed
$V_{s}$
) and swimming diffusion (with anisotropic diffusivity tensor
$\unicode[STIX]{x1D63F}$
, and magnitude
$V_{s}^{2}/d_{r}$
, where
$d_{r}$
is the rotational diffusivity). No-slip and no-flux boundary conditions are applied at the tube wall.
Equations (2.1) and (2.2) describe the dynamics of a suspension of gyrotactic swimming algae. As mentioned in § 1, such suspensions self-concentrate into plumes; cells focus in a plume at the centre of a vertically aligned tube and drive axial flow. However, we are free to label (or dye) locally some of these cells (which we call a ‘slug’) without affecting the cell concentration profile or the flow. The theory presented here describes the dispersion of a labelled slug of algae in the flow associated with an existing plume. The dynamics of the moments of a distribution of labelled cells are evaluated, providing information on axial drift and diffusion. We summarise the derivation of the main results as follows. Steady axisymmetric solutions will be obtained for (2.1) and (2.2) (steady plume solutions). Then, moment equations from (2.2) will be used to evaluate dispersion measures, using the mean cell response to flow from two alternative stochastic microscopic models. Finally, the predictions will be analysed for direct comparison with experiments.
2.1 The steady plume background
First, consider the background flow within which the slug of cells will disperse, provided by the steady solutions to (2.1) and (2.2) (a steady plume). As in Bees & Croze (Reference Bees and Croze2010), we assume a plume with no blips or varicose instabilities has formed in a long pipe, so that the dynamics are translationally invariant along the axial direction
$z$
. We further assume axisymmetric solutions for the cell concentration
$n(r)$
and flow deviation about the mean
$\unicode[STIX]{x1D712}(r)$
. In cylindrical polars, the flow field is
$\boldsymbol{u}(r)=u(r)\boldsymbol{e}_{z}=U[1+\unicode[STIX]{x1D712}(r)]\boldsymbol{e}_{z}$
, with
$U$
the mean flow. We non-dimensionalise lengths by the pipe radius,
$\hat{x}=x/a$
, times by
$\unicode[STIX]{x1D70F}_{d}=a^{2}d_{r}/V_{s}^{2}$
, the time to diffuse across it, such that
$\hat{t}=t/\unicode[STIX]{x1D70F}_{d}$
, and concentrations by the background mean cell concentration
$n_{av}$
, such that
$R_{0}^{0}=n/n_{av}$
(this notation is chosen for consistency with Bees & Croze (Reference Bees and Croze2010); see next section). Equations (2.1) and (2.2) then become (omitting hats for clarity)
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{r}\frac{\text{d}}{\text{d}r}\left(r\frac{\text{d}\unicode[STIX]{x1D712}}{\text{d}r}\right)=P_{z}-Ri\,R_{0}^{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}R_{0}^{0}}{\text{d}r}=\unicode[STIX]{x1D6FD}\frac{q^{r}}{\unicode[STIX]{x1D60B}^{rr}}R_{0}^{0}, & \displaystyle\end{eqnarray}$$
subject to no-slip and no-flux boundary conditions at the tube surface, and the integral constraints
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D712}}=2\int _{0}^{1}r\unicode[STIX]{x1D712}(r)\,\text{d}r=0,\quad \overline{R_{0}^{0}}=2\int _{0}^{1}rR_{0}^{0}(r)\,\text{d}r=1,\end{eqnarray}$$
where overbars denote cross-sectional averages,
$\overline{f}(z,t)=(1/\unicode[STIX]{x03C0})\!\int _{0}^{2\unicode[STIX]{x03C0}}\!\!\int _{0}^{1}f(r,\unicode[STIX]{x1D703},z,t)r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}$
. If
$f$
is independent of
$\unicode[STIX]{x1D703}$
, then
$\overline{f}(z,t)=2\int _{0}^{1}f(r,z,t)\,\text{d}r$
. Equations (2.5) follow from the definition of
$\unicode[STIX]{x1D712}$
and normalisation of the cross-sectional average concentration. In (2.3) and (2.4), we have defined the dimensionless parameters
$$\begin{eqnarray}P_{z}=\frac{1}{\unicode[STIX]{x1D707}}\frac{\text{d}p}{\text{d}z}\frac{a}{U},\quad Ri=\frac{v_{c}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g}{\unicode[STIX]{x1D707}}\frac{a^{2}}{U}n_{av},\quad \unicode[STIX]{x1D6FD}=\frac{ad_{r}}{V_{s}}.\end{eqnarray}$$
Here
$P_{z}$
represents the dimensionless pressure gradient in the
$z$
direction;
$Ri$
is a Richardson number quantifying the relative importance of buoyancy and viscous flow (
$\unicode[STIX]{x1D707}$
is the dynamic viscosity);
$\unicode[STIX]{x1D6FD}$
is the swimming Péclet number gauging the relative importance of advection by swimming to diffusion by swimming (see also (2.7)). Alternatively,
$\unicode[STIX]{x1D6FD}=1/Kn$
, where
$Kn=V_{s}d_{r}^{-1}/a$
is a swimming Knudsen number: the ratio of the cell mean free path to the tube radius. To close (2.3) and (2.4), we require the radial components of the swimming direction,
$q^{r}=q^{r}(\unicode[STIX]{x1D712}^{\prime })$
, and diffusivity tensor,
$\unicode[STIX]{x1D60B}^{rr}(\unicode[STIX]{x1D712}^{\prime })$
, where
$\unicode[STIX]{x1D712}^{\prime }\equiv \text{d}\unicode[STIX]{x1D712}/\text{d}r$
. These are known functions of shear obtained from stochastic models of gyrotactic response to flow and are discussed in § 2.3.
2.2 Slug dispersion
Next, we consider the dynamics of a slug of dyed cells dispersing within the plume flow background provided by the solutions of (2.3) and (2.4). A suspension of dyed swimmers replaces cells in the pipe within an existing plume, forming a dyed slug. The dyed swimmers are identical to the background cells in the plume and the slug dynamics are also governed by (2.2). However, since the cells are initially localised and subsequently disperse, the solutions are no longer independent of the axial coordinate or time. The slug concentration, denoted by
$n_{s}$
(non-dimensionalised by
$n_{av}$
as above, but not yet cross-sectionally averaged) is governed by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}n_{s}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D735}n_{s})-Pe(1+\unicode[STIX]{x1D712})\frac{\unicode[STIX]{x2202}n_{s}}{\unicode[STIX]{x2202}z}-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(n_{s}\boldsymbol{q}),\end{eqnarray}$$
where
$$\begin{eqnarray}Pe=\frac{Uad_{r}}{V_{s}^{2}}\end{eqnarray}$$
is the flow Péclet number quantifying the relative importance of advection by the flow to swimming diffusion.
Except for a few special cases, equation (2.7) is not amenable to analytical solution. It can be solved numerically using a spatially adaptive finite-element method, as described in Bearon et al. (Reference Bearon, Bees and Croze2012). This method allows one to probe the transient dispersion of a slug and the approach to steady state. To characterise dispersion, it is convenient to transform to a frame of reference moving with the mean flow,
$z\rightarrow z-Pe\,t$
, and consider axial moments of concentration. The
$p$
th radially varying moment is defined as
$c_{p}(r,t)=\int _{-\infty }^{+\infty }z^{p}n_{s}(z,r,t)\,\text{d}z$
, and its cross-sectional average as
$m_{p}(t)=\overline{c_{p}}=2\int _{0}^{1}rc_{p}(r,t)\,\text{d}r$
. We define the transient drift and diffusivity of the dispersing slug as
$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{0}(t)\equiv \frac{\text{d}}{\text{d}t}m_{1}(t),\quad D_{e}(t)\equiv \frac{1}{2}\frac{\text{d}}{\text{d}t}[m_{2}(t)-m_{1}^{2}(t)].\end{eqnarray}$$
An alternative to solving (2.7) is to solve for the axial moments
$c_{p}(r,t)$
and cross-sectionally averaged moments
$m_{p}(t)$
(Bees & Croze Reference Bees and Croze2010). The theory then allows one to predict that at long times the drift,
$\unicode[STIX]{x1D6EC}_{0}$
, and effective diffusivity,
$D_{e}$
, of a dyed slug of swimmers are given by (with the notation
$f^{\prime }=\text{d}f/\text{d}r$
)
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}_{0} & \equiv & \displaystyle \lim _{t\rightarrow \infty }\unicode[STIX]{x1D6EC}_{0}(t)=-\overline{\unicode[STIX]{x1D60B}^{rz}R_{0}^{0\prime }}+\overline{(Pe\unicode[STIX]{x1D712}+\unicode[STIX]{x1D6FD}q^{z})R_{0}^{0}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle D_{e} & \equiv & \displaystyle \lim _{t\rightarrow \infty }D_{e}(t)=-\overline{\unicode[STIX]{x1D60B}^{rz}g^{\prime }}+\overline{(Pe\unicode[STIX]{x1D712}+\unicode[STIX]{x1D6FD}q^{z}-\unicode[STIX]{x1D6EC}_{0})g}+\overline{\unicode[STIX]{x1D60B}^{zz}R_{0}^{0}},\end{eqnarray}$$
where
$R_{0}^{0}(r)$
is the background normalised concentration of cells from (2.3) and (2.4), and
$$\begin{eqnarray}g(r)=R_{0}^{0}\int _{0}^{r}\left(\frac{\unicode[STIX]{x1D60B}^{rz}(s)}{\unicode[STIX]{x1D60B}^{rr}(s)}-\frac{1}{2}\frac{\unicode[STIX]{x1D6EC}_{0}^{\ast }(s),-\unicode[STIX]{x1D6EC}_{0}m_{0}^{\ast }(s)}{s\unicode[STIX]{x1D60B}^{rr}(s)R_{0}^{0}(s)}\right)\text{d}s,\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{0}^{\ast }(r)=2\int _{0}^{r}s(-\unicode[STIX]{x1D60B}^{rz}{R_{0}^{0}}^{\prime }+(Pe\unicode[STIX]{x1D712}+\unicode[STIX]{x1D6FD}q^{x})R_{0}^{0})\,\text{d}s,\quad m_{0}^{\ast }(r)=2\int _{0}^{r}sR_{0}^{0}(s)\,\text{d}s.\end{eqnarray}$$
The function
$g(r)$
(first axial moment) gauges how the centre-of-mass (mean axial position) of a dyed slug varies radially. This controls the effective diffusivity
$D_{e}$
of the slug together with the cell concentration distribution
$R_{0}^{0}(r)$
(zeroth moment), which prescribes the slug drift
$\unicode[STIX]{x1D6EC}_{0}$
. The function
$\unicode[STIX]{x1D6EC}_{0}^{\ast }(r)$
in (2.13) is the radially averaged slug cell flux (due to anisotropic diffusion, flow and swimming) up to radius
$r$
. Evaluated at the wall (
$r=1$
), this provides the full non-dimensional drift:
$\unicode[STIX]{x1D6EC}_{0}^{\ast }(1)=\unicode[STIX]{x1D6EC}_{0}$
. Similarly,
$m_{0}^{\ast }(r)$
represents the radially averaged slug cell concentration up to
$r$
. At the wall, this gives
$m_{0}^{\ast }(1)=1$
(equivalent to (2.5)), which is the total non-dimensional number of cells reflecting cell conservation.
To obtain experimental predictions from (2.10) and (2.11), we require functional forms for the flow about the mean,
$\unicode[STIX]{x1D712}$
, the average swimming direction,
$\boldsymbol{q}$
, and diffusivity tensor,
$\unicode[STIX]{x1D63F}$
. Expressions for
$\boldsymbol{q}$
and
$\unicode[STIX]{x1D63F}$
for gyrotactic algae can be obtained from microscopic stochastic models, as detailed in the next section. If we neglect reciprocal coupling of the flow to the cells via cell negative buoyancy, we have Poiseuille flow, and expressions for
$\boldsymbol{q}$
and
$\unicode[STIX]{x1D63F}$
are given in Bearon et al. (Reference Bearon, Bees and Croze2012). Accounting for reciprocal coupling requires solving (2.3) and (2.4) for
$\unicode[STIX]{x1D712}$
and
$R_{0}^{0}$
, subject to the constraints (2.5). From this we can make predictions for dispersion using the integrals (2.10) and (2.11). Numerically it is easier to solve the equivalent ordinary differential equation (ODE) system (see next section and appendix E for details).
Dimensional and non-dimensional parameter values and ranges used in the dispersion model predictions. The D. salina gyrotactic motility and buoyancy parameters were estimated from tracking video microscopy measurements (see text and appendix A).
2.3 Shear dependence of transport, model parameters and solution
The components of the mean cell orientation,
$\boldsymbol{q}$
, and diffusivity,
$\unicode[STIX]{x1D63F}$
, in (2.7) and (2.10)–(2.12) gauge average gyrotactic cell response to shear flow. They can be obtained from the orientational moments of the probability density function (p.d.f.) for cell orientation and possibly position in a flow from stochastic continuum descriptions. As in previous work, we consider two such models, denoted F and G (described as Fokker–Planck and generalised Taylor dispersion models in the literature). In model F (Pedley & Kessler Reference Pedley and Kessler1990), the cell p.d.f. depends only on orientation, while in model G (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003; Bearon et al.
Reference Bearon, Bees and Croze2012) it depends on position and orientation. Models F and G were originally derived for unbounded linear shear flows (Bees, Hill & Pedley Reference Bees, Hill and Pedley1998; Hill & Bees Reference Hill and Bees2002) (see appendix C for a model equation summary). However, they also provide a good approximation for the swimmers pipe (radius
${\sim}a$
) flow, if the flow shear is linear on the scale of a cell (radius
$r_{c}$
):
$r_{c}\ll a$
(Bearon et al.
Reference Bearon, Bees and Croze2012). A separate requirement is that the cell diffusion is not in the Knudsen regime,
$Kn\,(=\unicode[STIX]{x1D6FD}^{-1})\sim 1$
, in which the cell random walk is strongly affected by walls. Using parameters from table 1 and
$r_{c}\sim R$
, where
$R\sim 10^{-4}$
cm is the cell hydrodynamic radius (see appendix C), we see that
$r_{c}/a\sim 10^{-4}$
and
$Kn\sim 0.1$
, so that both conditions are well-satisfied in our experiments. In both models F and G, a key dimensionless parameter is the stochasticity parameter
$\unicode[STIX]{x1D706}=1/(2d_{r}B)$
, expressing the relative importance of random reorientation from rotational diffusion
$d_{r}$
to reorientation by gravity at a rate
$1/(2B)$
. The p.d.f. can be solved using a Galerkin method for a given value of
$\unicode[STIX]{x1D706}$
. This has been achieved previously for Chlamydomonas augustae (Bees et al.
Reference Bees, Hill and Pedley1998; Bearon et al.
Reference Bearon, Bees and Croze2012). Here, we have obtained solutions for D. salina (see below and appendix A), allowing the evaluation of the functions
$q^{i}(\unicode[STIX]{x1D70E})$
and
$\unicode[STIX]{x1D60B}_{m}^{ij}(\unicode[STIX]{x1D70E})$
, where
$m=F$
or
$G$
and
$\unicode[STIX]{x1D70E}=-\unicode[STIX]{x1D712}^{\prime }Pe/(2\unicode[STIX]{x1D6FD}^{2})$
, the dimensionless shear rate. The model subscript in the average orientation components has been omitted intentionally: the predictions are the same for both models. The difference is in the diffusivity components, reflecting the different physics in the models. In model F the cell orientation distribution alone determines spatial dispersion, but in model G both orientation and position (leading to differential advection by swimming and flow) contribute to dispersion via spatial moments. Indeed, predictions for diffusion from model F are strictly valid at low shear rates, and break down for large shear rates (see Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013). Model G is valid for all shear rates, provided the shear flow can be locally approximated as linear. The description is appropriate for the Poiseuille flows considered here, but not straining flows (Bearon, Hazel & Thorn Reference Bearon, Hazel and Thorn2011). Both models provide approximations to the full Smoluchowski advection–diffusion equation for the cell p.d.f. (Doi & Edwards Reference Doi and Edwards1986), which could be solved numerically directly for the cell p.d.f. (Saintillan & Shelley Reference Saintillan and Shelley2006). However, this route does not allow one to take advantage of the simplification and analytical results that have been obtained for the F and G models (Pedley & Kessler Reference Pedley and Kessler1992; Bearon et al.
Reference Bearon, Bees and Croze2012).
The dispersion measurements described in § 2.5 provide the first experimental test of the validity of the microscopic stochastic models described above. The dispersion model was parametrised from independent tracking video microscopy measurements, carried out on suspensions of D. salina in vertically oriented capillaries. The model parameters thus obtained and used in predictions are summarised in table 1. The stochasticity parameter
$\unicode[STIX]{x1D706}$
from tracking was used to determine the specific forms of
$q^{i}(\unicode[STIX]{x1D70E})$
and
$\unicode[STIX]{x1D60B}_{m}^{ij}(\unicode[STIX]{x1D70E})$
for D. salina (see appendix A). Tracking microscopy was also used to quantify the settling of heat-immobilised cells, which allows one to infer the buoyant mass
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}v_{c}$
of D. salina cells. More details of the tracking measurements and error estimation can be found in appendix C. The model was solved numerically using MATLAB (Mathworks, Natick, MA, USA). As detailed in appendix D, it is easier to recast the dispersion prediction (integral equations (2.10) and (2.11) and dependent equations) as a boundary value ODE problem in which the coupled flow, cell concentration and dispersion are evaluated together. It should be noted that the dimensionless parameters of the dispersion and microscopic models are not all independent. In particular, the Richardson number depends on
$Pe$
and
$\unicode[STIX]{x1D6FD}$
:
$Ri=\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D6FD}^{3}/Pe)$
. The coupling of cells and flow thus depends on the imposed flow and the coupling strength
$\unicode[STIX]{x1D6FE}\equiv (v_{c}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g/\unicode[STIX]{x1D707})(V_{s}/d_{r}^{2})n_{av}$
, which we can control by changing the background mean cell concentration
$n_{av}$
.
2.4 Dispersion predictions
Before turning to a comparison between dispersion theory and experiment, it is useful to consider predictions for the distribution of gyrotactic algae in pipe flow and the ensuing dispersion. Similar predictions have been obtained using the swimming parameters of C. augustae (Bearon et al.
Reference Bearon, Bees and Croze2012), and also for channel geometries (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013), comparing predictions for the F and G microscopic models. These studies assumed that the coupling via buoyancy between cells and flow dynamics is negligible. Reciprocal coupling was considered by Bees & Croze (Reference Bees and Croze2010), but only with asymptotic results from the F model. Here, we present the first solutions of the dispersion equations with reciprocal coupling by buoyancy and nonlinear expressions from the stochastic microscopic descriptions F and G, and use them to investigate the dispersion of D. salina, the fully parametrised candidate species for experimental tests of the theory. We consider long-time dispersion, i.e. when axial moments are stationary in time. This is true when
$t\gg 1/\unicode[STIX]{x1D701}_{1}^{2}$
, where
$\unicode[STIX]{x1D701}_{1}$
is the smallest positive eigenvalue arising in the moment equations for cell conservation;
$\unicode[STIX]{x1D701}_{1}(Pe)$
is a monotonically increasing function of
$Pe$
, which is laborious to evaluate (see section 4 of Bees & Croze (Reference Bees and Croze2010) for details). Simulations of gyrotactic dispersion in channels confirm that steady dispersion is reached faster for larger
$Pe$
and, further, that the drift (measured in this work) becomes steady faster than the diffusivity (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013). For
$Pe=670$
, close to the largest
$Pe=560$
value studied here, these simulations show that the drift is steady at
$t\approx 0.5$
. That is to say, restoring dimensions, at a time
$t_{c}\approx \unicode[STIX]{x1D70F}_{d}/2$
, where
$\unicode[STIX]{x1D70F}_{d}=a^{2}d_{r}/V_{s}^{2}$
is the swimming diffusivity time scale. Multiplying by the mean flow speed
$U$
, and recalling the definition of
$Pe$
(2.8), this implies that dispersion will be steady beyond a distance
$Z_{c}\approx Pe\,a/2$
from the inoculation point (for general
$Pe$
,
$Z_{c}\approx Pe\,a/\unicode[STIX]{x1D701}_{1}^{2}$
). Thus, with
$a=0.35$
cm (see table 1) we estimate that for the largest
$Pe=560$
used in our experiments,
$Z_{c}\approx 100$
cm, the distance downstream from dyed cell injection at which our imaging was carried out (see § 3.1). We can thus assume that our experimental drift data are steady and it is reasonable to compare them with long-time dispersion predictions (see figure 4).
(a) Dimensionless normalised cell concentration profiles for different values of Péclet number,
$Pe$
(from (2.4)) using the G model for neutrally buoyant cells (
$Ri=0$
, no cell–flow coupling). The distribution is increasingly focused with increasing
$Pe$
, unlike the F model predictions (inset). (b,c) Reciprocally coupled cell concentration and flow as a function of mean background cell concentration
$n_{av}$
for the G model and
$Pe=1000$
. The coupling generates sharper focusing and distorts the flow from the parabolic Poiseuille shape. Insets: for the F model the coupling hardly affects the cell and flow distributions. Note that, since
$V_{s}$
and
$d_{r}$
are fixed, changing
$Pe$
corresponds to changing the mean speed
$U$
(or flow rate), as in our dispersion experiments.
Theoretical dispersion measures: drift and diffusivity in the absence (a,c) and presence (b,d) of reciprocal coupling via buoyancy effects. Models F (dashed) and G (solid) agree at low values of
$Pe$
for neutrally buoyant cells (a,c insets, where overlapping F and G predictions are shown), but predict different drift behaviour for large
$Pe$
. Coupling via negative cell buoyancy quantitatively changes the predictions of model G: negative buoyancy increases cell drift and enhances diffusivity (greater differential advection due to modified vorticity profile, see text) with respect to the neutrally buoyant case. On the other hand, coupling affects model F predictions very little (b,d insets).
In figure 1 we consider predicted cell distributions. For the case with no reciprocal coupling (neutral buoyancy; Richardson number,
$Ri=0$
), shown in figure 1(a), both the F and G models predict gyrotactic focusing. However, the distributions are different. For the G model, as the Péclet number
$Pe$
is increased, the distribution goes from almost uniform to sharply focused (the width diminishes monotonically). On the other hand, as
$Pe$
increases, the F model (figure 1
a, inset) displays a broader distribution (with a non-monotonic width progression). The differences between the two models, which are noticeable for
$Pe=1000$
, arise from the different average microscopic transport these models predict (Bearon et al.
Reference Bearon, Bees and Croze2012). In model G, local cell dispersion depends on reorientation by the flow and differential advection due to the spatial p.d.f.; in the F model the local spatial distribution is neglected (a good approximation only for weak flows). As a result, the models predict the same mean cell orientation
$q^{i}$
, but different diffusivities
$\unicode[STIX]{x1D60B}^{ij}$
as a function of dimensionless shear
$\unicode[STIX]{x1D70E}$
(Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003; Bearon et al.
Reference Bearon, Bees and Croze2012). The cell distribution is determined by the ratio of the radial components of these transport parameters,
$q^{r}/\unicode[STIX]{x1D60B}^{rr}$
, as is clear from (2.4). In model F, because
$\unicode[STIX]{x1D60B}^{rr}$
asymptotes to a constant value at large
$\unicode[STIX]{x1D70E}$
(Bearon et al.
Reference Bearon, Bees and Croze2012), gyrotactic advection in high shear (near the pipe walls) is weak compared to diffusion, so the cell distribution has significant tails at high
$Pe$
(see figure 1
a, inset). In model G, on the other hand,
$\unicode[STIX]{x1D60B}^{rr}\rightarrow 0$
at large
$\unicode[STIX]{x1D70E}$
(Bearon et al.
Reference Bearon, Bees and Croze2012): gyrotactic advection is strong compared to diffusion near the walls, and the distribution is sharp, with low cell concentration close to the walls.
Next, we consider the distributions of cells and flow with the additional coupling via the cells’ negative buoyancy (figure 1
b). As discussed above, the strength of coupling
$\unicode[STIX]{x1D6FE}$
for a given
$Pe$
depends on the cell concentration
$n_{av}$
. In model G, for a given
$Pe$
, increasing
$n_{av}$
results in a sharper distribution near the centre of the tube with respect to the uncoupled case. Concomitantly, cells that have accumulated close to the tube centre due to gyrotaxis drive flow faster there and cause the vorticity profile to deviate from zero around the origin more strongly than for the solutions with
$Ri=0$
(figure 1
c). The larger the value of
$n_{av}$
, the sharper the accumulation and the stronger the flow deviation,
$\unicode[STIX]{x1D712}$
. This does not occur in model F because cells on average spend more time by the walls (in the radial cell distribution tails). This prevents accumulations that would lead to significant reciprocal coupling via negative buoyancy. The model F distributions for concentrated cell suspensions (with non-zero reciprocal coupling) differ little from the
$Ri=0$
case (figure 1
b,c, inset).
The distribution of flow and cells determines dispersion, measured by the steady-state drift above the mean flow and effective diffusivity shown in figure 2; the skewness vanishes at long times (Bees & Croze Reference Bees and Croze2010). Figure 2(a) presents a comparison of the dimensionless drift above the mean flow,
$\unicode[STIX]{x1D6EC}_{0}$
, predicted for models F and G in the absence of reciprocal cell–flow coupling. Both models predict
$\unicode[STIX]{x1D6EC}_{0}>0$
for sufficiently large
$Pe$
, the distinguishing characteristic of gyrotactic dispersion (Bees & Croze Reference Bees and Croze2010; Bearon et al.
Reference Bearon, Bees and Croze2012; Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013). As expected, at low
$Pe$
(figure 2
a, inset) the F and G predictions coincide and, for very low
$Pe$
,
$\unicode[STIX]{x1D6EC}_{0}$
is negative because cells are oriented upwards (Bees & Croze Reference Bees and Croze2010). At higher
$Pe$
, the predictions of the two models diverge: in G,
$\unicode[STIX]{x1D6EC}_{0}\sim Pe$
; in F, it tends to a constant. The reason for this is evident from the distributions of figure 1(a). The increased concentration around the centreline with
$Pe$
in model G causes the larger drift (more cells in fast flow). In model F, on the other hand, the number of cells in the fast central zone of the pipe saturates, and with it the drift. Reciprocal cell–flow coupling enhances the focusing of the distribution for model G, as well as the magnitude of the flow, so a greater value of
$\unicode[STIX]{x1D6EC}_{0}$
is possible for the same
$Pe$
. In model F, coupling makes little difference to the drift (figure 2
b, inset).
The predicted effective diffusivity of the dispersing population is shown in figure 2(c) in the absence of reciprocal coupling. It increases monotonically in this
$Pe$
range for both the F and G models, which, as previously, agree for very low
$Pe$
. However, we see crucial differences. In model G,
$D_{e}$
appears to saturate, while in model F it continues increasing with
$Pe$
. To understand this, note that the axial dispersion measured by
$D_{e}$
is caused principally by diffusion across streamlines and differential advection of cells distributed by gyrotaxis. In model G, as
$Pe$
is increased, the cell distribution sharpens and fewer cells are exposed to the high shear close to the pipe wall; the increase of
$D_{e}$
thus drops at high
$Pe$
. In model F, on the other hand, the distribution tails ensure that cells will reside in regions of high shear, so
$D_{e}$
increases with
$Pe$
. The dispersion in the presence of coupling is shown in figure 2(d). We see that in model G, larger cell concentrations provide a larger
$D_{e}$
for the same
$Pe$
. The greater vorticity in the flow close to the tube centre (providing greater differential advection) contributes to a greater dispersion in the reciprocal coupling case. The effect of reciprocal coupling on diffusion is not very significant in model F, whose predictions for
$D_{e}$
are not very different from the uncoupled model (figure 2
d, inset), in contrast to model G.
2.5 Connection with experiment
Experimentally we measure the absolute drift of a slug of fluorescently labelled cells within a plume, given by
$V_{d}=\text{d}Z/\text{d}t$
, where
$Z$
is the vertical position of the peak intensity of recorded plume (see experimental methods, § 3.1). A useful measure of the drift is the fractional drift above the mean flow,
$\unicode[STIX]{x1D6FF}=(V_{d}-U)/U$
, where
$U$
is the mean flow speed. This fractional drift can be compared with the prediction from the dispersion model. Recall that the theoretical drift is defined in (2.9) as
$\unicode[STIX]{x1D6EC}_{0}=\text{d}\hat{m}_{1}/\text{d}\hat{t}$
, with
$\hat{t}=t(a^{2}D_{r}/V_{s}^{2})^{-1}$
,
$\hat{m}_{1}=(\langle Z\rangle -Ut)/a$
and
$Pe$
defined in (2.8). It then follows that the theoretical fractional drift is simply given by
$\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6EC}_{0}/Pe$
.
3 Experiments
3.1 Methods
3.1.1 Algae and culture
Batch cultures of D. salina (CCAP 19/18) were grown on a 12/12 light/dark cycle on a modification of the medium of Pick, Karni & Avron (Reference Pick, Karni and Avron1986), which we will denote herein as DSM (D. salina medium). Healthy motile cells were gravitactically concentrated overnight on rafts of cotton wool at the surface of a culture, in a similar fashion to Croze et al. (Reference Croze, Ashraf and Bees2010). These cultures were kept under the same light cycle as batch cultures to preserve circadian rhythms. Light intensity, however, was reduced to minimise the disruption of upswimming by phototaxis during concentration. Cells were harvested using a Pasteur pipette. Cell concentrations were obtained from readings of
$A$
, the suspension absorbance at
$590~\text{nm}$
, measured using a spectrophotometer (WPA CO7500), calibrated to haemocytometer counts; concentrations are given by
$c=5\times 10^{6}\,A~\text{cells}~\text{cm}^{-3}$
, with fractional error
$\unicode[STIX]{x0394}c/c=0.07$
.
3.1.2 Fluorescent staining
Cells were stained with fluorescein diacetate (FDA, Sigma), a vital dye. We mixed 8
$\unicode[STIX]{x03BC}$
l of dye from a 5
$\unicode[STIX]{x03BC}$
M stock (acetone) into 20 ml of cell suspension with absorbance typically in the range
$A=0.7$
–
$1.5$
. After letting the dye act for approximately
$5$
min, dyed cells were ‘washed’ by letting them swim up through a raft of cotton overlaid by a fresh layer of DSM medium. This step ensures that the dye is only inside the cells. However, since fluorescein cleaved from FDA inside the cell slowly leaks out, dispersion experiments were performed within
$20$
min of washing to allow macroscopic imaging with good contrast between dyed cells and the background medium.
3.1.3 Experimental set-up and imaging
A schematic of the apparatus is shown in figure 3(a). A Graseby 3500 syringe pump (Graseby Medical Ltd, Watford, UK) attached at the top of the set-up drives flow in a vertically oriented Perspex tube that is
$2$
m in length and
$7$
mm in diameter. Experiments were performed with flow rates in the range 100–400 ml h
$^{-1}$
at a stable ambient temperature of
$T=24\pm 1\,^{\circ }$
C. Subpopulations of D. salina, fluorescently labelled as above, were introduced into the flow using a
$1$
ml syringe at the top of the tube,
$2$
cm downstream from the pump; dyed cells were injected into a steady flow of unlabelled cells. The ensuing dispersion does not depend on the rate of injection because this rate,
${\sim}10^{3}$
ml h
$^{-1}$
, was much larger than the largest advection flow rate used. The absorbance of labelled and unlabelled populations was matched. We imaged the ensuing dispersion using a GE680C charge-coupled device (CCD) camera (Prosilica) and excitation/emission filters (MF475-35/MF525-39, Thorlabs) suitable for the fluorescence of fluorescein from FDA whose absorption/emission peaks are 495 nm/520 nm (blue/green). The imaging was centred on a point
$L=60$
or
$100$
cm below the inoculation point. Sequences were captured at
$2.1$
Hz.
(a) Experimental set-up: fluorescently labelled slugs of D. salina dispersing in pipe flow, driven by a syringe pump. Labelled cells are injected near the top of the tube, recruited to an existing plume, and imaged downstream. Image (b) shows a blippy slug for a flow rate
$\unicode[STIX]{x1D6F7}=100$
ml h
$^{-1}$
,
$n_{av}=3.65\times 10^{6}$
cells cm
$^{-3}$
at
$L=60$
cm from injection,
$[\unicode[STIX]{x1D6F7},n_{av},L]=[100,3.6,60]$
. The next images show blip-free plumes: (c)
$[200,4.5,60]$
, (d)
$[200,3.65,100]$
, (e)
$[300,2.95,100]$
, and (f)
$[400,2.95,100]$
(faint). Images are
$1.1$
cm
$\times$
11.0 cm.
3.1.4 Image analysis
Dispersion measures were obtained from captured image sequences of descending slugs. The total pixel intensity in four regions of interest (ROI) of an image was calculated for each image in a sequence of a given plume. The time
$t$
at which the intensity peaked corresponded to the time at which the plume passed the vertical position
$Z$
(on which the ROI was centred). Thus, a linear fit to a plot of
$Z$
versus
$t$
provides a value for the plume drift speed
$V_{d}=\text{d}Z/\text{d}t$
. As discussed in the next section, it was very challenging to obtain non-blippy, high-contrast plumes. For low-contrast plumes, it was possible to estimate the drift from the displacement of the plume endpoints. The error in the drift was estimated as the standard deviation over measurements with the same sample (same mean concentration
$n_{av}$
).
3.2 Plumes and blips
In order for the experiments to test meaningfully the theory of gyrotactic dispersion, the fluorescently stained cells have to be recruited to the plumes formed by non-stained cells. It is clear from the images shown in figure 3(b–f) that fluorescent D. salina swimmers are indeed recruited to the plumes. A blip instability has occurred in figure 3(b): fluorescence is localised in the blip nodules. The theory of dispersion described in § 2 assumes axial invariance of the background cell concentration and flow, which is broken when blip instabilities arise. The instability occupies certain regions of the flow–concentration parameter space, as predicted for channel flows for C. augustae parameters by Hwang & Pedley (Reference Hwang and Pedley2014b ). A similar stability analysis has not been carried out for D. salina in pipe flow. The dispersion experiments are challenging because of competing constraints. In particular, moderately concentrated suspensions are necessary for good fluorescent image contrast, but high concentrations are more likely to produce blips (by increasing the suspension Richardson number). However, with some care, it is possible to find blip-free plumes over the length and time scales of dispersion experiments.
3.3 Results: drift above the mean flow
The dispersion data were used to evaluate the drift speed, and hence the fractional drift above the mean flow,
$\unicode[STIX]{x1D6FF}=(V_{d}-U)/U$
, where
$V_{d}$
is the measured drift and
$U$
is the mean flow speed, as described in the methods. The predictions from the reciprocally coupled dispersion theory for the F and G models are shown in figure 4 for different values of the background mean cell concentration,
$n_{av}$
. In the absence of reciprocal cell–flow coupling (
$Ri=0$
), the model G fractional drift
$\unicode[STIX]{x1D6FF}$
increases monotonically with
$Pe$
, increasing to
$1$
for large
$Pe$
. The model F prediction on the other hand at first increases before reaching a peak and decreasing to zero (see figure 6 in appendix F). We can interpret these results in terms of the distribution of cells, in a similar manner to how we interpreted the drift in figure 2(a). The cells in model G can be sharply focused to travel close to the maximum flow speed
$V_{d}=2U$
(
$\unicode[STIX]{x1D6FF}=1$
). In model F, on the other hand,
$\unicode[STIX]{x1D6FF}$
tends to zero at large
$Pe$
(see figure 6 in appendix F), reflecting the uniformity of the cell distribution. Coupling for model G increases the value of
$\unicode[STIX]{x1D6FF}$
significantly above
$1$
at intermediate
$Pe$
, but, since
$Ri\sim 1/Pe$
, at large
$Pe$
the fractional drift tends to
$1$
. For model F, reciprocal coupling has little effect due to the uniformity of the cell distribution.
Comparison of theoretical predictions and preliminary experimental observations of fractional drift in dispersing slugs of D. salina, for different background cell concentrations,
$n_{av}$
. Same colour lines (theory) and symbols (experiment) denote equal concentrations. Quantitatively, predictions using microscopic model G (generalised Taylor dispersion) and accounting for cell–flow coupling agree much better with experiment than using model F (Fokker–Planck model), even for the lowest
$Pe$
data (inset, symbols/line types as in main panel). Dashed black lines indicate predictions neglecting reciprocal coupling due to negative buoyancy (
$Ri=0$
). For clarity, the model G prediction for
$n_{av}=10^{6}$
cells cm
$^{-3}$
(pink line) is shown only in the inset. Vertical error bars (from same-run repeats) are smaller than the data point symbol size. Horizontal error bars stem from uncertainties in
$d_{r}$
(see appendix C).
The experimental results provide good evidence for
$\unicode[STIX]{x1D6FF}>0$
: the dispersion of gyrotactic swimmers is not the same as that of passive tracers. Furthermore, the experimental drift values are consistent with the predictions of model G in the presence of reciprocal cell–flow coupling. Model F (with or without coupling) agrees poorly with the experimental observations, even for the lowest experimental values of
$Pe$
studied (figure 4, inset). The data suggest that it is necessary to include reciprocal cell–flow coupling in theoretical descriptions of fast-flowing suspensions of gyrotactic swimmers.
4 Discussion and outlook
We have carried out the first experimental test of the theory of gyrotactic dispersion in pipe flow and associated stochastic microscopic models. We injected D. salina algae, fluorescently labelled by fluorescein diacetate, into an unlabelled population in downwelling pipe flow, and imaged the dispersing ‘slugs’ to measure their fractional drift above the mean flow. Experimental drift values were compared with predictions from an axial dispersion model (Bees & Croze Reference Bees and Croze2010; Bearon et al. Reference Bearon, Bees and Croze2012) with swimming and buoyancy parameters measured using tracking microscopy. The dispersion model requires as inputs transport functions (local mean orientation and cell diffusivity tensor) dependent on the local shear rate, which can be evaluated from stochastic microscopic models. Predictions using two such models, referred to as models F (Pedley & Kessler Reference Pedley and Kessler1990, Reference Pedley and Kessler1992) and G (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003; Bearon et al. Reference Bearon, Bees and Croze2012), were employed with parameters for D. salina. We derived new solutions of the axial dispersion model combining the effects of both negative buoyancy and nonlinear solutions of the two alternative stochastic microscopic models. Bees & Croze (Reference Bees and Croze2010) had previously considered only linearised transport functions for the F model with negative buoyancy.
The experimental results provide clear evidence for the non-zero drift above the mean flow expected in the dispersion of gyrotactic swimmers (Bees & Croze Reference Bees and Croze2010; Bearon et al.
Reference Bearon, Bees and Croze2012; Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013) (see figure 4). This is in stark contrast to the zero drift of passive, neutrally buoyant tracers (Taylor Reference Taylor1953; Aris Reference Aris1956). Furthermore, in the
$Pe$
range studied, we find quantitative agreement between experiment and theoretical predictions obtained using model G and accounting for negative cell buoyancy (
$Ri>0$
). As shown in table 3 of appendix F, three out of the five experimental points agree well with theory within the stated errors. For the other two points (corresponding to the lowest cell concentrations studied), theory underestimates the drift. Agreement with model F, on the other hand, is poor with or without buoyancy, even for the lowest values of
$Pe$
(figure 4, inset). This should be expected for the relatively large
$Pe$
used (recall that models F and G agree only at low
$Pe$
), as local spatial stochastic effects in strong shear flow, neglected in model F, become important. The superior performance of model G over model F expected at large
$Pe$
has been shown previously by Croze et al. (Reference Croze, Sardina, Ahmed, Bees and Brandt2013), who compared analytical dispersion predictions for models F and G with Lagrangian simulations. This work, however, ignored the effect of negative cell buoyancy. Without the latter, agreement between the experimental results and model G is not very good. This highlights the important role of reciprocal flow–cell coupling in the dispersion dynamics.
Given the reciprocal coupling between cells and flow, it is very difficult to design experimental apparatus to test microscopic stochastic models of macroscopic suspensions of biased swimming cells. In this paper we have presented one method that has shown great promise, demonstrating the shortcomings of the F model whilst providing evidence in support of the G model with the inclusion of negative buoyancy effects. However, there were some significant experimental challenges. The experiments were limited by the constraints of good image contrast, steady dispersion and the emergence of blip instabilities. The latter proved antithetic to the former two, as high concentrations provided better contrast but also a greater propensity for blip formation for the range of
$Pe$
investigated. In addition, the use of long tubes to guarantee steady dispersion made blips more likely (giving slow-growing unstable modes time to develop). Blips are known to depend sensitively on concentration and flow rate (Kessler Reference Kessler1986; Dennisenko & Lukaschuk Reference Dennisenko and Lukaschuk2007), but we lack a quantitative guide to avoid the instability. The work of Hwang & Pedley (Reference Hwang and Pedley2014b
) provides only a qualitative guide, as this study focused on C. augustae, which has different gyrotactic parameters (Hill & Häder Reference Hill and Häder1997).
Future experiments should refine the drift measurements and test other dispersion predictions, such as the dependence of the effective diffusivity
$D_{e}$
of the swimming algae on key flow and swimming parameters. Fluorescent imaging with at least two downstream cameras will allow measurements of the axial variance
$\text{Var}(Z)$
of a dispersing slug. Theory predicts
$\text{Var}(Z)=2D_{e}V_{s}^{2}d_{r}^{-1}t$
from which
$D_{e}$
can be evaluated. The practical issues associated with fluorescent dye diffusion and loss of contrast may be circumvented using mutants (e.g. Chlamydomonas reinhardtii, expressing fluorescent protein). Our experience in experimenting with such approaches provides two important caveats: mutants from ‘normally swimming’ wild-type background strains should be employed; and motility should not be compromised by the insertion of fluorescent gene constructs.
The dispersion models we have tested experimentally in this work are predicated on swimmer suspension theories developed in the early 1990s (see Pedley & Kessler (Reference Pedley and Kessler1992) for a review). While several improvements have been made to these theories, they have not been subject to extensive experimental tests since their formulation. In particular, basic assumptions of the theories, such as the dilute nature of the swimmer suspensions, have not been challenged experimentally. In this work, we follow the traditional assumption that swimmer suspensions are dilute, but, considering the formation of plumes, one can question if the dilute assumption is appropriate for the description of self-focusing suspensions of gyrotactic swimmers and axial dispersion. The cell concentration at the centre of the pipe depends on the applied flow rate. From figure 1(b), we find that for
$Pe=1000$
(and the parameters in table 1), the concentration at the pipe centre is approximately 100 times that of the average, i.e.
$10^{8}$
cells cm
$^{-3}$
. This corresponds to a volume fraction of 0.02 and a mean distance between cells of 20
$\unicode[STIX]{x03BC}$
m (two cell diameters). Therefore, we might expect cell–cell interactions between cells to play a role (Drescher et al.
Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011). However, we note that the cell concentration drops very rapidly from the centre, so high-concentration effects should match the dilute theory except for a narrow core at the centre of the pipe, which is advected at close to the maximum flow rate. Furthermore, axial cell diffusivity is dominated by cell motion in the high-shear regions and not within the narrow core. Consequently, the axial dispersion predictions of the dilute theory are likely to hold even for sharply focused plumes. The rationale of this work was precisely to ascertain how well the predictions from dilute theory perform by testing them experimentally. The experimental results suggest that the dilute assumption is not a major issue. However, it would be worthwhile to extend the theory to explore a semidilute regime.
Another assumption made in the original theories (again, refer to Pedley & Kessler (Reference Pedley and Kessler1992)) is that the swimmer population heterogeneity is unimportant. It is assumed that averages of swimming parameters over a population (e.g. mean swimming speed) are adequate to describe the suspension behaviour and fluid mechanics. This will clearly only be a good approximation if parameter distributions are sufficiently narrow. Bees et al. (Reference Bees, Hill and Pedley1998) studied theoretically the effect of deviations from narrow distributions on gyrotactic transport in the F model. They evaluated how the ratio of speed variance to mean,
$\langle V^{2}\rangle /\langle V\rangle ^{2}$
(
$V$
denotes speed, angled brackets an average over the population), changes the components of the diffusivity tensor
$\unicode[STIX]{x1D63F}$
. We note that, as for the dilute suspension assumption, the agreement we find between theory and experiments seems to indicate that the neglect of swimming population heterogeneity is not a serious one in the dispersion of D. salina. This cannot, however, be generalised to other swimming species, which may have broad swimming parameter distributions. Furthermore, even within a particular species, biological changes, such as metabolism, can drive significant heterogeneity. It would be of interest to consider how swimming parameter distributions and their changes affect suspension dynamics and dispersion.
Accurately predicting the unusual dispersion of swimming algae in pipe flow is important for the design of more efficient and ‘considerate’ photobioreactors for their cultivation (Bees & Croze Reference Bees and Croze2014). Closed photobioreactor designs comprise channels or pipes in which growing suspensions of algae are flowed. In tubular external loop air-lifts, cells are bubbled with air in a riser tube and then recirculated in a downcomer connected to the riser by two horizontal tubes. For gyrotactic biotechnological species of interest, such as the
$\unicode[STIX]{x1D6FD}$
-carotene producer D. salina considered here, the dispersion observed in this work is very relevant. For example, the non-zero drift above the mean flow we have observed will cause nutrients or waste (passive tracers) to become separated from the cells in the downcomer. Experiments are in progress quantifying the behaviour of gyrotactic swimmers in laboratory and pilot-scale photobioreactors. More complex dispersion may ensue in photobioreactors operated in turbulent pipe flow. Simulations predict gyrotactic dispersion due to cells focusing transiently in downwelling regions of the turbulent channel flows (Croze et al.
Reference Croze, Sardina, Ahmed, Bees and Brandt2013), but experiments on swimmers in bounded turbulent flows have not yet been carried out. It will be of considerable ecological interest to explore the consequences of gyrotactic dispersion in river, lake and oceanic flows.
Acknowledgements
We thank P. Dominy and J. Christie for advice on fluorescent vital dyes; G. Gibson for help with fluorescent imaging; and J. Arlt and D. Dell’Arciprete for help with video microscopy of D. salina. We also acknowledge discussions with P. Cicuta, S. Furlan, Y. Hwang, D. Jin, V. Martinez, T. Pedley and E. Silvester. O.A.C. and M.A.B. gratefully acknowledge support from the Carnegie Trust for the Universities of Scotland and EPSRC (EP/D073398/1). O.A.C. acknowledges support from a Royal Society Research Grant and the Winton Programme for the Physics of Sustainability.
Appendix A. Stochastic models and shear-dependent transport
We summarise very briefly the mathematical structure of models F (the orientation-only Fokker–Planck model) (Pedley & Kessler Reference Pedley and Kessler1990, Reference Pedley and Kessler1992) and G (generalised Taylor dispersion) (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003). The reader is referred to the original literature for more details.
Consider
$P(\boldsymbol{p},\boldsymbol{x},t)$
, the probability of finding a cell with orientation
$\boldsymbol{p}$
at position
$\boldsymbol{x}$
at time
$t$
. This evolves according to (Hill & Bees Reference Hill and Bees2002)
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }[(\boldsymbol{u}+V_{s}\boldsymbol{p})P]+d_{r}\unicode[STIX]{x1D735}_{\boldsymbol{p}}\boldsymbol{\cdot }[\dot{\boldsymbol{p}}P-\unicode[STIX]{x1D735}_{\boldsymbol{p}}P]=0, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}$
is the fluid velocity,
$V_{s}$
is the constant cell swimming speed,
$d_{r}$
is the rotational diffusivity due to the intrinsic randomness in cell swimming, and
$\unicode[STIX]{x1D735}_{\boldsymbol{x}}$
and
$\unicode[STIX]{x1D735}_{\boldsymbol{p}}$
are physical and orientational gradients, respectively. For gyrotactic cells,
$$\begin{eqnarray}\dot{\boldsymbol{p}}=\unicode[STIX]{x1D706}(\boldsymbol{k}-(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{p})\boldsymbol{p})-\unicode[STIX]{x1D70E}\boldsymbol{j}\wedge \boldsymbol{p},\end{eqnarray}$$
where, as in the main text,
$\unicode[STIX]{x1D706}=1/(2d_{r}B)$
is the stochasticity parameter and
$\unicode[STIX]{x1D70E}=-\unicode[STIX]{x1D712}^{\prime }Pe/(2\unicode[STIX]{x1D6FD}^{2})$
is the dimensionless shear strength. The unit vectors
$\boldsymbol{k}$
and
$\boldsymbol{j}$
are in the upwards (antiparallel to gravity) and positive vorticity direction, respectively.
A.1 Model G
By calculating the moments of the distribution function
$P$
, Frankel & Brenner (Reference Frankel and Brenner1991, Reference Frankel and Brenner1993) showed that, on time scales long compared to
$1/d_{r}$
, the concentration of cells,
$n(\boldsymbol{x},t)$
, in a homogeneous shear flow satisfies an advection–diffusion equation (Hill & Bees Reference Hill and Bees2002)
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\left[\left(\boldsymbol{u}+V_{s}\boldsymbol{q}\right)n-D_{c}\,\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}n\right]=0, & & \displaystyle\end{eqnarray}$$
where
$D_{c}=V_{s}^{2}/d_{r}$
is the magnitude of the swimming diffusivity. The mean swimming direction,
$\boldsymbol{q}$
, and non-dimensional diffusivity tensor,
$\unicode[STIX]{x1D63F}$
, can be written as integrals over cell orientation,
$\boldsymbol{p}$
(Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003), such that
$$\begin{eqnarray}\displaystyle \boldsymbol{q} & = & \displaystyle \int _{\boldsymbol{p}}\boldsymbol{p}\,f(\boldsymbol{p})\,\text{d}\boldsymbol{p},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63F} & = & \displaystyle \int _{\boldsymbol{p}}\left[\boldsymbol{b}\boldsymbol{p}+\frac{2\unicode[STIX]{x1D70E}}{f(\boldsymbol{p})}\boldsymbol{b}\boldsymbol{b}\boldsymbol{\cdot }\hat{\boldsymbol{G}}\right]^{sym}\,\text{d}\boldsymbol{p}.\end{eqnarray}$$
Here
$[~]^{sym}$
denotes the symmetric part of the tensor, and
$\hat{\boldsymbol{G}}$
denotes the non-dimensional fluid velocity gradient. The equilibrium orientation,
$f(\boldsymbol{p})$
, and vector
$\boldsymbol{b}(\boldsymbol{p})$
are given by (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003)
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{p}}\boldsymbol{\cdot }[\dot{\boldsymbol{p}}\,f-\unicode[STIX]{x1D735}_{\boldsymbol{p}}\,f] & = & \displaystyle 0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{p}}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\boldsymbol{b}-\unicode[STIX]{x1D735}_{\boldsymbol{p}}\boldsymbol{b})-2\unicode[STIX]{x1D70E}\boldsymbol{b}\boldsymbol{\cdot }\hat{\boldsymbol{G}} & = & \displaystyle f(\boldsymbol{p})(\boldsymbol{p}-\boldsymbol{q}),\end{eqnarray}$$
subject to the integral constraints
$\int _{\boldsymbol{p}}f\,\text{d}\boldsymbol{p}=1$
and
$\int _{\boldsymbol{p}}\boldsymbol{b}\,\text{d}\boldsymbol{p}=0$
. Equations (A 6) and (A 7) follow from evaluating the moments of the distribution function
$P$
(Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003).
The orientation distribution
$f(\boldsymbol{p})$
represents the steady-state distribution of swimmer orientations. The meaning of the vector
$\boldsymbol{b}(\boldsymbol{p})$
is less intuitive. It can be written as
$\boldsymbol{b}=f\boldsymbol{B}d_{r}/V_{s}$
, where
$\boldsymbol{B}(\boldsymbol{p})$
quantifies the difference between the average position of a swimmer, given its instantaneous orientation is
$\boldsymbol{p}$
, and its average position averaged over all values of
$\boldsymbol{p}$
(Hill & Bees Reference Hill and Bees2002).
A.2 Model F
In model F, the same advection–diffusion equation (A 3) is assumed, but the transport functions are obtained from a p.d.f. in terms of orientation only. The distribution
$f$
for the orientation is governed by (A 6) with (A 2). The solution for the mean orientation (which depends only on
$f$
) is the same as in model G. However, in model F the diffusivity tensor is given by the approximation
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63F}_{F}=D_{c}\int _{\boldsymbol{p}}(\boldsymbol{p}-\boldsymbol{q})^{2}f(\boldsymbol{p})\,\text{d}\boldsymbol{p}, & & \displaystyle\end{eqnarray}$$
where now
$D_{c}=V_{s}^{2}\unicode[STIX]{x1D70F}$
, and
$\unicode[STIX]{x1D70F}$
is a direction correlation time estimated from experiment. This is not as straightforward as measuring the rotational diffusivity
$d_{r}$
required for the G model. However, following Bearon et al. (Reference Bearon, Bees and Croze2012), we can use the fact that models F and G must coincide asymptotically at low shear, so that
$\unicode[STIX]{x1D70F}=J_{1}/(K_{1}\unicode[STIX]{x1D706}d_{r})$
, where values of
$\unicode[STIX]{x1D706}$
and
$d_{r}$
can be found from experiment (see appendix C), and
$J_{1}(\unicode[STIX]{x1D706})$
and
$K_{1}(\unicode[STIX]{x1D706})$
are rationally derived functions that can be calculated (for more details, please refer to appendix D of Bearon et al. (Reference Bearon, Bees and Croze2012) and Pedley & Kessler (Reference Pedley and Kessler1990)).
A.3 Model solutions
The dispersion model requires knowledge of the functional dependence on the dimensionless shear
$\unicode[STIX]{x1D70E}$
of the components of the transport tensors
$\unicode[STIX]{x1D666}$
and
$\unicode[STIX]{x1D63F}$
. We follow the approach of Bearon et al. (Reference Bearon, Bees and Croze2012) where these components are evaluated from solutions of models F and G using a Galerkin method. The solutions depend on the bias parameter
$\unicode[STIX]{x1D706}=1/(2d_{r}B)$
evaluated from tracking experiments. Previously, solutions were obtained for models parametrised for C. augustae algae. Here we present new results for D. salina, for which
$\unicode[STIX]{x1D706}=0.21$
(see table 1 in the main text). For analytical convenience, these solutions are then fitted with rational functions. The relevant components for dispersion are
$q^{r}(\unicode[STIX]{x1D70E})=-\unicode[STIX]{x1D70E}P(\unicode[STIX]{x1D70E};\boldsymbol{a}^{r},\boldsymbol{b}^{r})$
,
$q^{z}(\unicode[STIX]{x1D70E})=-P(\unicode[STIX]{x1D70E};\boldsymbol{a}^{z},\boldsymbol{b}^{z})$
,
$\unicode[STIX]{x1D60B}_{m}^{rr}(\unicode[STIX]{x1D70E})=P(\unicode[STIX]{x1D70E};\boldsymbol{a}^{rr},\boldsymbol{b}^{rr})$
,
$\unicode[STIX]{x1D60B}_{m}^{zz}(\unicode[STIX]{x1D70E})=P(\unicode[STIX]{x1D70E};\boldsymbol{a}^{zz},\boldsymbol{b}^{zz})$
,
$\unicode[STIX]{x1D60B}_{m}^{rz}(\unicode[STIX]{x1D70E})=-\unicode[STIX]{x1D70E}P(\unicode[STIX]{x1D70E};\boldsymbol{a}^{rz},\boldsymbol{b}^{rz})$
, where the rational function
$P(\unicode[STIX]{x1D70E};\boldsymbol{a},\boldsymbol{b})$
is given by
$$\begin{eqnarray}\displaystyle P(\unicode[STIX]{x1D70E};\boldsymbol{a},\boldsymbol{b})=\frac{a_{0}+a_{2}\unicode[STIX]{x1D70E}^{2}+a_{4}\unicode[STIX]{x1D70E}^{4}}{1+b_{2}\unicode[STIX]{x1D70E}^{2}+b_{4}\unicode[STIX]{x1D70E}^{4}}. & & \displaystyle\end{eqnarray}$$
The parameters
$a_{0}$
,
$a_{2}$
,
$a_{4}$
,
$b_{2}$
and
$b_{4}$
(table 2) are chosen to ensure the polynomials match asymptotic solutions to the stochastic equations of models F and G (Bearon et al.
Reference Bearon, Bees and Croze2012).
Parameters used in the rational function fits of solutions for models F and G used in our dispersion theory predictions.
Appendix B. Modified Pick medium
The D. salina cells were grown in a modification of Pick’s medium (Pick et al.
Reference Pick, Karni and Avron1986), with final salt concentration of
$0.5$
M. The medium consists of (J. Polle, private communication): a
$100\times$
concentrated nutrient mix (2 M Tris-HCl, pH 7.5; 1 M
$\text{KNO}_{3}$
; 1 M
$\text{MgSO}_{4}$
;
$60$
mM
$\text{CaCl}_{2}$
;
$20$
mM
$\text{KH}_{2}\text{PO}_{4}$
; 0.4 mM
$\text{FeCl}_{3}$
in 4 mM EDTA, pH 7.5); a
$1000\times$
concentrated trace micronutrients (150 mM
$\text{H}_{3}\text{BO}_{3}$
; 10 mM
$\text{MnCl}_{2}$
; 0.8 mM
$\text{ZnCl}_{2}$
; 0.3 mM
$\text{CuCl}_{2}$
; 2 mM
$\text{Na}_{2}\text{MoO}_{4}$
; 2 mM
$\text{NaVO}_{3}$
; 0.2 mM
$\text{CoCl}_{2}$
); a
$4\times$
concentrated salt solution (
$2$
M
$\text{NaCl}$
); and a
$20\times$
concentrated bicarbonate solution (0.5 M
$\text{NaHCO}_{3}$
, filter sterilised and added after autoclaving the nutrients and salt mixed together).
Appendix C. Model parameters from microscopy
The swimming speed
$V_{s}$
and bias parameter
$\unicode[STIX]{x1D706}$
were obtained by tracking video microscopy of a suspension of D. salina swimming in a vertically oriented capillary. The mean speed
$V_{s}$
is obtained straightforwardly from the distribution of cell displacements. The distribution of angles to the vertical was inferred from the tracking data and fitted with a Von Mises distribution, as in Hill & Häder (Reference Hill and Häder1997):
$$\begin{eqnarray}f(\unicode[STIX]{x1D703})=\frac{\exp [\unicode[STIX]{x1D706}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{0})]}{2\unicode[STIX]{x03C0}I_{0}(\unicode[STIX]{x1D706})}.\end{eqnarray}$$
The fit allowed us to estimate
$\unicode[STIX]{x1D706}$
. Alternatively,
$\unicode[STIX]{x1D706}$
can be found from the ratio of the mean upswimming speed
$V_{u}$
to the mean swimming speed
$V_{s}$
, since
$V_{u}/V_{s}=K_{1}(\unicode[STIX]{x1D706})=\coth (\unicode[STIX]{x1D706})-1/\unicode[STIX]{x1D706}$
(Pedley & Kessler Reference Pedley and Kessler1992). The video data was captured at 100 frames per second (f.p.s) in videos lasting 40 s. These captured D. salina cells swimming in a 400
$\unicode[STIX]{x03BC}$
m deep flat glass chamber (CM Scientific, 50 mm
$\times$
8 mm
$\times$
0.4 mm). Particle tracking using the established algorithms of Crocker & Grier (Reference Crocker and Grier1996) was implemented via the Kilfoil laboratory MATLAB routines (http://people.umass.edu/kilfoil/downloads.html). This typically provided
${\approx}2500$
tracks per video that were analysed to obtain the values of
$V_{s}$
and
$\unicode[STIX]{x1D706}$
shown in table 1. Values of
$\unicode[STIX]{x1D706}$
obtained from the alternative methods discussed above from the same data were averaged; values from independent datasets were also averaged. The final error in the parameter was estimated by evaluating the standard error in the mean of these independent data.
Tracking also provides an estimate from the gravitational reorientation time
$B$
. It was estimated from the reorientation of cells immobilised by heating the suspension to
$40\,^{\circ }$
C. Tracking the reorientation of the sedimenting cells provided
$B$
by comparing with the analytical solution of (A 2) (with
$\unicode[STIX]{x1D70E}=0$
) for the reorientation of
$\unicode[STIX]{x1D703}$
:
$\unicode[STIX]{x1D703}(t)=2\arctan [\exp (-(1/2B)(t-t_{0}))]$
. We find
$B=10.5\pm 0.8$
s. With
$B$
and
$\unicode[STIX]{x1D706}$
known, the value for the rotational diffusivity follows from the definition of
$\unicode[STIX]{x1D706}=1/(2Bd_{r})$
. The errors in the tracking parameters determine the error in derived quantities. In particular, by standard error propagation it can be shown that the error in
$Pe$
is dominated by the error in the rotational diffusivity. This uncertainty
$\unicode[STIX]{x0394}Pe/Pe\approx \unicode[STIX]{x0394}d_{r}/d_{r}=0.18$
is sizeable, and is represented as horizontal error bars in figure 4.
Tracking microscopy of sedimenting cells (the same data that were used to estimate
$B$
) also allows one to quantify their negative buoyancy. The drag force on a sedimenting cell can be estimated from
$F=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}RU_{s}$
(generalised Stokes law), where
$R$
is the effective hydrodynamic radius of a sphere with the same viscous resistance to flow as the cell and
$U_{s}$
is the average sedimentation speed of the cells. These were heat-immobilised as above and their settling was captured in 100 s long videos at 2.52 fps and tracked (as above) to provide
$U_{s}=(1.44\pm 0.37)\times 10^{-4}~\text{cm}~\text{s}^{-1}$
. For sedimenting negatively buoyant cells, the buoyant force
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}v_{c}g$
, where
$v_{c}$
is cell volume and
$g$
is the gravitational acceleration, balances the viscous drag. Thus we can estimate
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}v_{c}=6\unicode[STIX]{x03C0}\frac{\unicode[STIX]{x1D707}}{g}RU_{s},\end{eqnarray}$$
if
$R$
is provided. This radius can be estimated from the angularly averaged friction factor of a prolate ellipsoid, our chosen approximation for the D. salina cells. Using an electrostatic analogy it can be shown that for ellipsoids
$R=2/\unicode[STIX]{x1D709}_{0}$
(Hubbard & Douglas Reference Hubbard and Douglas1993), where
$$\begin{eqnarray}\unicode[STIX]{x1D709}_{0}=\int _{0}^{\infty }\frac{1}{\sqrt{(a^{2}+\unicode[STIX]{x1D701})(b^{2}+\unicode[STIX]{x1D701})(c^{2}+\unicode[STIX]{x1D701})}}\,\text{d}\unicode[STIX]{x1D701}\underset{b=c<a}{=}\frac{2\cosh ^{-1}(a/b)}{\sqrt{a^{2}-b^{2}}}.\end{eqnarray}$$
By image analysis of micrographs, we estimate
$a=(4.8\pm 2.7)\times 10^{-4}$
cm and
$b~(=c)=(3.2\pm 1.9)\times 10^{-4}~\text{cm}$
. With these values
$\unicode[STIX]{x1D709}_{0}=(5.4\pm 1.6)\times 10^{-4}$
cm
$^{-1}$
and
$R=(3.7\pm 1.1)\times 10^{-4}~\text{cm}$
, so that, with
$\unicode[STIX]{x1D707}$
and
$g$
from table 1, we can use (C 2) to find
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}v_{c}=(0.92\pm 0.36)\times 10^{-11}$
g. This product is all that is required to parametrise the model: knowledge of the independent values of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$
and
$v_{c}$
is not necessary. Calculating the product
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}v_{c}$
rather than its factors has the added benefit of reducing the experimental uncertainty in our estimate with respect to the one evaluated from the factors. Apart from general values for biflagellates quoted by Kessler (Reference Kessler1986), we could not find measurements of the excess density of D. salina in the literature. However, a fluorometric study of phytoplankton sinking rates included measurements of Dunaliella tertiolecta (a close relative of D. salina), providing its excess density of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.1$
and estimated cell volume
$2.7\times 10^{-10}~\text{cm}^{3}$
(Eppley, Holmes & Strickland Reference Eppley, Holmes and Strickland1967). This gives
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}v_{c}=2.7\times 10^{-11}$
g, which is of the same order as our estimate for D. salina.
Appendix D. Model numerical solution
To obtain predictions from the steady-state dispersion model, we need to solve (2.10) and (2.11) with the coupled flow and cell concentrations given by (2.3) and (2.4), subject to the constraints (2.5). As mentioned in the main text, this problem is more simply tackled by writing down the equivalent differential equation system (a boundary value problem). Recalling that primes denote derivatives with respect to
$r$
, the dispersion ODE system is given by
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D712}^{\prime }}{\text{d}r}=-\frac{1}{r}\unicode[STIX]{x1D712}^{\prime }+P_{z}-Ri\,R_{0}^{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}R_{0}^{0}}{\text{d}r}=\unicode[STIX]{x1D6FD}\frac{q_{r}}{\unicode[STIX]{x1D60B}^{rr}}R_{0}^{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}m_{0}^{\ast }}{\text{d}r}=2rR_{0}^{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}X}{\text{d}r}=2r\unicode[STIX]{x1D712}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}P_{z}}{\text{d}r}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D6EC}_{0}^{\ast }}{\text{d}r}=2r[\unicode[STIX]{x1D60B}^{rz}R_{0}^{0\prime }+(Pe\unicode[STIX]{x1D712}+\unicode[STIX]{x1D6FD}q^{z})R_{0}^{0}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}g}{\text{d}r}=\unicode[STIX]{x1D6FD}\frac{q^{r}}{\unicode[STIX]{x1D60B}^{rr}}g+R_{0}^{0}\frac{\unicode[STIX]{x1D60B}^{rz}}{\unicode[STIX]{x1D60B}^{rr}}-\frac{\unicode[STIX]{x1D6EC}_{0}^{\ast }-m_{0}^{\ast }\unicode[STIX]{x1D6EC}_{0}}{2r\unicode[STIX]{x1D60B}^{rr}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}D_{e}}{\text{d}r}=2[g^{\prime }+(Pe\unicode[STIX]{x1D712}+\unicode[STIX]{x1D6FD}q^{z}-\unicode[STIX]{x1D6EC}_{0})g+\unicode[STIX]{x1D60B}^{zz}R_{0}^{0}], & \displaystyle\end{eqnarray}$$
subject to the boundary conditions
$\unicode[STIX]{x1D712}^{\prime }(0)=0$
,
$\unicode[STIX]{x1D712}(1)=-1$
,
$m^{\ast }(0)=0$
,
$m^{\ast }(1)=1$
,
$X(0)=0$
,
$X(1)=0$
,
$\unicode[STIX]{x1D6EC}_{0}^{\ast }(0)=0$
,
$g(0)=0$
and
$D_{e}(0)=0$
.
Equations (D 1) and (D 2) correspond to (2.3) and (2.4), respectively, for the coupled flow deviation above the mean
$\unicode[STIX]{x1D712}$
and (normalised) cell concentration
$R_{0}^{0}$
. Equation (D 3) corresponds to (2.13) for
$m_{0}^{\ast }$
, the cumulative cell concentration within a radius
$r$
. Equation (D 4) implements the constraint
$\overline{\unicode[STIX]{x1D712}}=0$
(from the definition of
$\unicode[STIX]{x1D712}$
) and (D 5) states that no pressure gradient is imposed. Equation (D 6) corresponds to (2.13) defining the partial drift. The full drift is then simply evaluated as
$\unicode[STIX]{x1D6EC}_{0}^{\ast }(1)=\unicode[STIX]{x1D6EC}_{0}$
. Equation (D 7) corresponds to (2.12) for the diffusivity weight function
$g$
and, finally, (D 8) corresponds to (2.11) for the effective axial diffusivity. Because this equation requires the value of the full drift, the system (D 1)–(D 6) needs to be solved once, before solving for the diffusivity. Numerical solutions of the ODEs were obtained using the MATLAB bvp4c routine. So that this routine could handle the singular boundary value problem in (D 7), we carried out the change of variables
$\tilde{g}=g/r$
and
$\tilde{D_{e}}=D_{e}/r$
. Solutions to the model equations were found using both model F and G solutions for
$q^{r}(\unicode[STIX]{x1D70E})$
and
$\unicode[STIX]{x1D60B}^{rr}(\unicode[STIX]{x1D70E})$
(see appendix A).
Appendix E. Blip instability
As indicated in the main text, several competing constraints made the observation of high-contrast non-blipping plumes very challenging. Figure 5 shows a diagram charting the parameter values where blips were observed. Note that sometimes blips were simply observed and image sequences not recorded for analysis (either for practical reasons or deliberately as we were not initially interested in blip data). The plot should not be considered a stability diagram, which would include a careful observation of the instability growth rates.
Chart recording the concentration–flow values giving blips or clean plumes at the time of observation of the dyed slug. This should not be considered a stability diagram: the syringe pump flow rate did not allow exploration of slow-growing instabilities.
Appendix F. Drift data and full model F predictions
The experimental values of the fractional drift above the mean flow
$\unicode[STIX]{x1D6FF}$
(figure 4) are provided in table 3 together with theoretical values predicted using model G, with which they compare favourably (a similar quantitative comparison with model F is not worthwhile, as is clear from the insets of figures 4 and 6). Bounds on the theoretical values were estimated by calculating the variations in the drift for upper and lower estimates of
$Pe$
and
$n_{av}$
, i.e.
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6FF}_{Pe}=\unicode[STIX]{x1D6FF}(Pe+\unicode[STIX]{x0394}Pe,n_{av})-\unicode[STIX]{x1D6FF}(Pe,n_{av})$
and similarly for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6FF}_{n_{av}}$
. The resulting uncertainties were then combined by standard error propagation. The other experimental uncertainties have been previously discussed (those in the drift in § 3.1, while those in
$Pe$
in appendix C). The uncertainties in
$n_{av}$
derive from the spectrophotometric calibration to concentration, discussed in § 3.1. As discussed in the text, the drift values predicted by model F are far from the experimentally observed values, whether reciprocal coupling due to buoyancy is included or not. This is clearly shown in figure 6.
Predictions of model F for the fractional drift above the mean flow
$\unicode[STIX]{x1D6FF}$
as a function of
$Pe$
. The theoretical curves are far from experimental drift values, as expected for a model that breaks down at large
$Pe$
(see discussion in § 4).
Experimental and theoretical values obtained using model G of the fractional drift compared for the values of
$Pe$
and
$n_{av}$
at which experiments were carried out. The evaluation of the uncertainties of the values shown is discussed in the text. No uncertainty is shown for experimental drift on the fourth line, as it was only possible to acquire one measurement of this drift (the single measurement error is negligible compared to the theoretical uncertainty estimates).
