## Abstract

Shear flow significantly affects the transport of swimming algae in suspension. For example, viscous and gravitational torques bias bottom-heavy cells to swim towards regions of downwelling fluid (gyrotaxis). It is necessary to understand how such biases affect algal dispersion in natural and industrial flows, especially in view of growing interest in algal photobioreactors. Motivated by this, we here study the dispersion of gyrotactic algae in laminar and turbulent channel flows using direct numerical simulation (DNS) and a previously published analytical swimming dispersion theory. Time-resolved dispersion measures are evaluated as functions of the Péclet and Reynolds numbers in upwelling and downwelling flows. For laminar flows, DNS results are compared with theory using competing descriptions of biased swimming cells in shear flow. Excellent agreement is found for predictions that employ generalized Taylor dispersion. The results highlight peculiarities of gyrotactic swimmer dispersion relative to passive tracers. In laminar downwelling flow the cell distribution drifts in excess of the mean flow, increasing in magnitude with Péclet number. The cell effective axial diffusivity increases and decreases with Péclet number (for tracers it merely increases). In turbulent flows, gyrotactic effects are weaker, but discernable and manifested as non-zero drift. These results should have a significant impact on photobioreactor design.

## 1. Introduction

In natural bodies of water and in industrial bioreactors, microscopic algae experience laminar and turbulent flows that play a critical role in their dispersion, proliferation and productivity [1,2]. The dispersion of passive tracers in fluid flows is well understood, particularly in the industrially relevant cases of flows in pipes and channels [3]. Taylor first realized that the dispersion of passive tracers in a pipe, caused by the combination of fluid shear and molecular diffusion, could be described in terms of an effective axial diffusivity [4]. A similar result also holds for turbulent pipe [5] and channel [6] flows. Taylor's pioneering analyses of dispersion in a pipe inspired a series of studies that placed the understanding of the dispersion of neutrally buoyant tracers on a firm footing [7,8]. Particles whose density differs from the suspending medium exhibit more complex behaviour (e.g. they can accumulate at walls in turbulent channel flows; see [9,10]).

Swimming single-celled algae are known to respond non-trivially to flows. For example, the mean swimming direction of biflagellates, such as *Chlamydomonas* and *Dunaliella* spp., is biased by imposed flows [11,12]. This bias, known as gyrotaxis, results from the combination of viscous torques on the cell body, owing to flow gradients, and gravitational torques, arising from bottom-heaviness and sedimentation [13]. In the absence of flow gradients, the gravitational torque leads cells to swim upwards on average (gravitaxis). Recent simulations and laboratory experiments have shown how inertia and gyrotaxis can play significant roles in the formation of patchiness and/or thin layers of toxic algae in the ocean [14–16]. Gravitaxis and gyrotaxis can also couple in a complex fashion to other biases owing to chemical gradients (chemotaxis) and light (phototaxis [17]). Therefore, the fact that cells can actively swim across streamlines and accumulate in specific regions of a flow forewarns that the resultant dispersion of biased swimming algae in a complex flow field is non-trivial.

In a recent study, Bees & Croze [18] extend the classical Taylor–Aris analysis of dispersion in a laminar flow in a tube to the case of biased swimming algae. The Bees & Croze dispersion theory provides a general theoretical framework to describe the dispersion of biased micro-swimmers in confined geometries, such as pipes and channels. To make predictions for the dispersion of particular organisms, specific functional expressions for the swimming characteristics need to be obtained from microscopic models. For example, using expressions for the mean swimming orientation **q** and diffusivity **D** from the so-called Fokker–Planck (FP) model [12], Bees & Croze [18] found that gyrotaxis can significantly modify the axial dispersion of swimming algae (such as *Chlamydomonas* spp.) in a pipe flow relative to the case of passive tracers (or dissolved chemicals). More recently, Bearon *et al.* [19] obtained predictions for swimming dispersion in laminar flows in tubes of circular cross section by using an alternative microscopic model known as generalized Taylor dispersion (GTD), formulated for swimming algae [20,21] and bacteria [22]. The GTD model is considered superior, as it incorporates correlations in cell positions owing to cell locomotion within a local flow, as well as the gyrotactic orientational biases considered in the FP description. The qualitatively distinct predictions, however, have yet to be tested. Direct tests of the models might involve microscopic tracking of isolated cells swimming in prescribed flows. The Bees & Croze theory allows the predictions of the two models to be tested for entire populations in macroscopic experiments [23] and individual-based numerical simulations, as we shall demonstrate in this paper.

Microalgae currently play a prominent role in biotechnology. For example, they are cultured commercially as nutritious fish and crustacean feed in aquaculture, and for high-value products that they can synthesize, such as β-carotene [24]. Microalgae also hold significant potential for future exploitation; they can provide a sustainable biofuel feedstock that does not compete for arable land with food crops. Despite significant efforts, however, current production systems remain too inefficient for algal biofuels to be commercially viable [25]. Algal culture commonly occurs in open or closed bioreactors, with production systems employing three stages: cell culture, harvest and downstream processing. Open bioreactors consist of one or more artificial ponds, quiescent or stirred, where stirring allows one to obtain greater concentrations per reactor area [24]. A common design is the raceway pond, in which rotating paddles generate flow (figure 1*a*). Typically, raceway ponds are shallow (depths *L* = 10–30 cm), which aids cell exposure to light. Typically, the flow is turbulent (speeds of 30 cm s^{−1}), facilitating light exposure and gas exchange. Pond bioreactors are relatively cheap to operate, but are limited by their low areal productivity and susceptibility to contamination. This restricts the range of viable species in open pond culture to robust, high-value extremophiles, which grow in prohibitive conditions such as high salinity or low pH. A greater variety of fast-growing algae can be cultured in closed bioreactors, where single species are propagated within sealed transparent containers. Common geometries are flat, cylindrical (columnar or tubular) and annular (figure 1*b,c*). As for raceways, the constraint of light penetration limits the smallest dimension to 10 cm. In tubular bioreactors, the flow speed typically is 50 cm s^{−1}, while flat panel and columnar/annular set-ups operate at 1–10 cm s^{−1} [26]. Although closed bioreactors allow for rapid growth conditions with small footprints, current designs are relatively expensive to manufacture, operate and maintain, and the algae are susceptible to invasion by other less-useful species.

Existing bioreactors operate under both laminar and turbulent flow conditions. The transition between the two flow regimes depends on the ratio of inertial and viscous forces in a fluid. This is expressed by the dimensionless Reynolds number *Re* = *UL*/*ν*, where *U* and *L* are the characteristic speed and length scale of the flow, respectively, and *ν* is the suspension kinematic viscosity [27]. In general, flows with *Re* > 2000 will be turbulent, though the transition to turbulence depends on geometry and particular flow conditions. Approximating the kinematic viscosity by that of water and using the flow speed and length scales above, we see that flows in air-lift reactors can be both laminar or turbulent (), while raceway ponds and tubular reactors are always turbulent (*Re* > 30 000). The mixing properties of turbulent flows are thought to be beneficial for cell growth, although clear evidence for this is lacking. With air-lifts the mixing is caused by rising bubbles, which also provide aeration, and the flow does not need to be turbulent. Both cells and nutrients will disperse in the above flows, with cell growth and productivity depending critically on such dispersion.

Current photobioreactor designs assume that cells disperse as tracers, whether they swim or not [28]. However, as detailed above, recent studies show that the distribution and consequent dispersion of swimming cells in flows, including biotechnologically interesting species such as *Dunaliella*, should be very different from that of passive tracers [18,19]. For example, if cells disperse differently to nutrients in a reactor they will separate, which could have catastrophic consequences for growth. Inspired by the possibility of taking advantage of the peculiarities of swimmer response to flow to improve bioreactor operation, we study here the dispersion of gyrotactic swimming algae in laminar and turbulent channel flows. We carry out time-resolved direct numerical simulations (DNSs), comparing results to predictions from the Bees & Croze dispersion theory applied to channel geometry. The paper is structured as follows. In §2, we present the simulation methods and outline a derivation of the theory for the new geometry. In §3, we present simulation results for the dispersion of biased swimming algae in laminar and turbulent flows, comparing theory and simulation results for laminar flows. In §4, we interpret the results, and discuss their implications for dispersion in photobioreactors and the environment.

## 2. Methods

### 2.1. Direct numerical simulations: governing equations

We simulate the dynamics of a population of biased swimming micro-organisms, typically *N* = 2 × 10^{5} and 10^{6} individuals per simulation, placed in laminar and turbulent flows. Each micro-swimmer is modelled as a spheroidal particle whose position **x*** _{i}*,

*i*= 1 …

*N*, evolves according to 2.1where

*V*

_{s}is the mean cell swimming speed,

**u**the local fluid velocity and

**p**the local particle orientation.

The orientation **p** of each swimmer evolves in response to the biasing torques acting upon it. Making the realistic simplifying assumption that swimming cells have a spheroidal geometry, the reorientation rate of the organisms is defined by the inertia-free balance of gravitational and viscous torques [12,29]
2.2where **k** is a unit vector in the vertical direction; *B* = *μ**α*_{⊥}/(2 *h**ρ**g*) is the gyrotactic reorientation time (*h* denotes the centre-of-mass offset relative to the centre-of-buoyancy, *α*_{⊥} is the dimensionless resistance coefficient for rotation about an axis perpendicular to **p**, *ρ* and *μ* are the fluid density and viscosity, respectively); *α*_{0} = (*a*^{2} − *b*^{2})/(*a*^{2} + *b*^{2}) is the eccentricity of the spheroids with major axis *a* and minor axis *b*; and, finally, **E** is the rate of strain tensor and ** ω** the vorticity. The noise

*Γ*

_{r}is added to simulate the stochastic rotational diffusivity of a swimming cell. In the simulations it is implemented with random angular steps of magnitude so that cells diffuse with rotational diffusivity

*d*

_{r}[30,31].

The flow field **u** is obtained by solving the Navier–Stokes equations for an incompressible viscous fluid, such that
2.3where D*/*D*t* ≡ ∂/∂*t* + (**u** · ∇) denotes the material derivative. Here, for simplicity, we shall assume cells to be neutrally buoyant, neglecting feedback from the particles to the flow (these effects are explored elsewhere [18]).

### 2.2. Direct numerical simulations: geometry, scaling and statistical measures of dispersion

We consider a standard channel geometry: two flat plates parallel to each other, infinite in extent, and separated by a gap of size 2*H*. When focusing on swimming cells it would seem natural to rescale length by the plate half-width *H* and time by where *V*_{s} is the mean swimming speed and *d*_{r} is the cell rotational diffusivity, the characteristic time a cell swimming with *V*_{s} takes to diffuse across *H* with diffusivity However, the parameter values obtained from this ‘cell-based’ rescaling are too large and thus numerically inconvenient for simulations. Thus, we adopt a ‘flow-based’ rescaling in terms of the characteristic length *H* and the flow-based time scale *H*/*U*_{c}, the time taken for a flow with centreline speed *U*_{c} to advect a cell by *H*. In terms of this rescaling, the dimensionless equations of motion for a biased swimming cell are
2.4
2.5
2.6where starred quantities denote dimensionless variables. In particular, we define the dimensionless swimming speed *v*_{s} = *V*_{s}/*U*_{c}, the gyrotactic parameter *η* = 2*BU*_{c}/*H* and the Reynolds number
2.7where *U*_{c} is the centreline speed, *H* is the channel half-width and *ν* the fluid's kinematic viscosity. Non-dimensionalizing the noise *Γ*_{r} also defines the dimensionless rotational diffusivity

We adopt a Cartesian coordinate system where the mean flow in the *x* (streamwise) direction varies with the wall-normal coordinate *y*, and is independent of *z* (spanwise direction). We integrate equations (2.4)–(2.6) numerically (see the electronic supplementary material) to find Enumerating the number of cells *N*(**x**) at a given position **x** in a bin of fixed volume *Δ**V*, we obtain the cell concentration *c*(**x**) = *N*(**x**)/*Δ**V*. From the cell coordinates we can further define statistical measures of the cell dispersion in a flow: the adimensional drift with respect to the mean flow *U* = (2/3)*U*_{c}; the effective streamwise diffusivity ; and the skewness of the cell distribution, *γ*. These measures of dispersion are given by
2.8where (*p* = 1,2,3) are the distribution moments and is the variance of the cell distribution. The statistical measures can be transformed from the flow-based scaling with characteristic time scale *H*/*U*_{c} to the cell-based scaling with time scale by the transformations
2.9where *Pe* is the cell Péclet number defined with respect to the centreline speed (see equation (2.13)). Note that the skewness *γ* does not depend on the scaling used. We shall present results in terms of the cell-based scaling and compare their limit at long times with analytical predictions.

### 2.3. Analytical theory: dispersion in laminar flows at long times

Here, we shall obtain analytical predictions for the dispersion of algae swimming in laminar channel flow that will be compared with the results from the simulations. The general Bees & Croze [18] continuum dispersion theory for biased swimmers will be applied to the new channel geometry; it is valid in the long-time limit (). The derivation for channel flow is similar to the pipe flow example in Bees & Croze, so we provide an outline here. Those readers who are less mathematically inclined may skip the derivation of the results in this section.

We shall begin with the continuity equation for the cell number density *n*:
2.10The flow velocity **u** is the solution of the Navier–Stokes equation (2.3). We adopt the same coordinate system described above for the DNS. For laminar flows such that the flow downstream is only a function of the wall-normal coordinate *y*, the flow velocity can be expressed as
2.11where *U* is the mean flow speed and *χ*(*y*) describes the variation in the streamwise direction about this mean. It is convenient to translate to a reference frame travelling with the mean flow, and non-dimensionalize length scales by *H* and time scales by the time to diffuse across the channel Thus, and where variables with hats denote dimensionless variables.

We assume unidirectional coupling of the cell dynamics to the flow (for bidirectional coupling owing to non-neutrally buoyant cells, see [18]): cells are biased by shear in the flow, so that the cell swimming velocity **V**_{c} = *V*_{s}**q**, where **q** is the mean swimming direction, and diffusivity tensor **D** are functions of local flow gradients. Here, we consider analytical predictions for spherical cells (*α*_{0} = 0), so cell orientation is only a function of vorticity,

Consider planar Poiseuille flow, *χ* = (1 − 3*y*^{2})/2 and *U* = (2/3)*U*_{c} [32], where *U*_{c} is the centreline (maximum) flow speed. Equation (2.10) becomes
2.12with
2.13where hats have been dropped for clarity. *Pe* is a Péclet number, the dimensionless ratio of the rates of transport by the flow and swimming diffusion, and *β* is the ratio of channel half-width to the length a cell swims before reorienting significantly. Alternatively, we can rewrite and interpret it as a ‘swimming’ Péclet number, the ratio of the rates of transport by swimming and diffusion. No-flow and no-flux boundary conditions will be applied to (2.12), such that
2.14where **n** is the unit vector normal to the channel boundary *Σ*.

As the flow is translationally invariant along *x*, the mean swimming direction and diffusivity tensor are independent of *x*: **D** = **D**(*y*), **q** = **q**(*y*). This permits a treatment of dispersion using moments in a similar vein to that of Aris [7]. The *p*th moment with respect to the axial direction is defined as
2.15provided *x ^{p}n*(

*x*,

*y*,

*t*) → 0 as

*x*→ ± ∞. We denote cross-sectional averages by overbars. The cross-sectionally averaged axial moment is thus 2.16

In Cartesian coordinates, equation (2.12) becomes 2.17

Multiplying by *x ^{p}* and integrating over the length of an infinite channel, we obtain the moment evolution equation
2.18Averaging over the cross section and applying the no-flux boundary conditions (2.14) yields
2.19from which we calculate measures of cell dispersion. The drift above the mean flow, is given by
2.20where
2.21is the zeroth axial moment (normalized concentration profile). Similarly, the effective diffusivity, is given by
2.22where with and The weighting function controls the drift and

*g*(

*y*) (related to the first axial moment) controls the value of the diffusivity.

To make predictions from (2.20) and (2.22), we require expressions for **D** and **q** from microscopic models of the statistical response of cells to flow. Two main models have been proposed for biased swimming algae: the FP model [33,34] and GTD [19–21]. For spherical cells, each model predicts how non-dimensional swimming direction **q** and diffusivity **D** depend on two non-dimensional quantities: σ(*y*) = −χ*‘*(*y*)*U*/(2*d*_{r}*H*) = −χ’(*y*)(1/3)(*Pe*/*β*^{2}), the ratio of reorientation time by rotational diffusion to that by shear (vorticity), and *λ* = 1/(2*d*_{r}*B*), the ratio between the time of reorientation by rotational diffusion and the gyrotactic reorientation time *B*. The functional forms for the transport parameters **q**(*σ*) and **D*** _{m}*(

*σ*), where subscripts

*m*= F and G denote the FP and GTD models, respectively, are given in appendix A [19]. In particular, the two models differ qualitatively in their predictions for the diffusivity as a function of the shear rate, and therefore they provide different predictions for cell distribution and dispersion. The dispersion predictions (2.20) and (2.22) were evaluated numerically with Matlab (Mathworks, Natick, MA, USA) using

**q**(

*σ*) and

**D**

*(*

_{m}*σ*) from the two models (see electronic supplementary material).

### 2.4. Parameters, simulation time and flow profiles

In simulations and theoretical evaluations we used the following cell parameters based on *Chlamydomonas augustae*: *V*_{s} = 0.01 cm s^{−}^{1} (swimming speed), *d*_{r} = 0.067 s^{−}^{1} (cell rotational diffusivity), and *B* = 3.4 s (gyrotactic reorientation time) [18]. Furthermore, we consider the flow of suspensions taken to have the same viscosity as water, *ν* = 0.01 cm^{2} s^{−}^{1}. The cell eccentricity *α*_{0} is also held fixed. The analytical theory presented above assumes that cells are spherical, *α*_{0} = 0, but laminar and turbulent simulations have also been performed for elongated cells with *α*_{0} = 0.8. With these parameters fixed, choosing the centreline speed *U*_{c} and channel width *H* gives the non-dimensional flow-based parameters used in the DNS: *v*_{s} = *V*_{s}/*U*_{c} (dimensionless swimming speed), *η* = 2*BU*_{c}/*H* (gyrotactic parameter), (dimensionless rotational diffusivity) and Reynolds number *Re* = *U*_{c}*H*/*ν*. For the test run with passive tracers an additional noise term was added to equation (2.4) to simulate the translational diffusivity *D _{t}*, non-dimensionalized as These flow-based parameters can be transformed into cell-based non-dimensional parameters for comparison with analytical predictions. From the definitions above and equation (2.13) it can be shown that (flow Péclet number), (swimming Péclet number) and (local dimensionless shear rate). Since

*χ′*(

*y*) = −3

*y*for plane Poiseuille flow, the maximum dimensionless shear is given by where we recall

*λ*= 1/(2

*Bd*

_{r}), the non-dimensional bias parameter. Simulations are more readily interpreted and compared with analytical theory in terms of these parameters (shown in table 1).

The dimensionless flow profiles corresponding to the Reynolds numbers used in this study are plotted in figure 2 for the benefit of the reader. Laminar flows are self-similar, so the non-dimensional flow has the same parabolic profile for all *Re*. Time-averged turbulent flows have distinct profiles that depend on the Reynolds number. Note that the number of degrees of freedom in the turbulent flow simulations scales as *Re*^{9} [27]. This makes large *Re* simulations computationally expensive. For this reason, we do not investigate dispersion for flows beyond *Re* = 10 000.

## 3. Results

### 3.1. Passive tracers: classical dispersion

The dispersion of passive tracers, such as molecular dyes or non-motile cells, is generally well understood. In laminar channel flow passive tracers are transported on average at the mean flow speed; there is no drift relative to mean flow: *Λ*_{0} = 0. The effective axial diffusivity *D*_{e} is given at long times by the Taylor–Aris result *D*_{e} = 1 + (8/945)*Pe*^{2} [3,32]. As a benchmark test, we carried out DNS for passive tracers (solving equation (2.4) with *v*_{s} = 0 and translational noise to simulate molecular diffusivity; see §2). A typical result is shown in figure 3 for *Pe* = 3000. As expected at long times *Λ*_{0} → 0 and the Taylor–Aris diffusivity prediction, *D*_{e} = 76 191, compares very well with simulation results. The skewness slowly tends to zero at long times as *t ^{−}*

^{0.5}(see the electronic supplementary material), suggesting approach to a normal distribution [7]. The above results for passive tracers will be compared with the dispersion of gyrotactic swimmers. They can also be used to check our analytical theory in the limiting case of unbiased swimmers with no gyrotactic response to flow. In this limit,

**q**= 0 and the diffusivity tensor is isotropic

**D**=

**I**/6, where

**I**is the identity matrix (see also [19]). Thus, in the calculation of equations (2.20) and (2.22),

*D*= 1/6 =

^{xx}*D*,

^{yy}*q*= 0 and

_{i}*D*= 0 (as and

^{xy}*J*(

*y*) = 0), so that

*Λ*

_{0}= 0, as expected for passive tracers, and

*D*

_{e}= 1/6 + 6(8/945)

*Pe*

^{2}. In the Bees & Croze theory,

*D*

_{e}and

*Pe*are scaled with respect to the diffusivity scale Using the scaling for random diffusers in three dimensions, we recover the classical Taylor–Aris result for channels.

The dispersion of passive tracers in turbulent channel flows has also been elucidated [8]. Elder derived the (dimensional) effective diffusivity *K* of passive tracers in turbulent open channel flow as *K* = 5.9 *u _{τ}H*, where

*u*is the friction velocity and

_{τ}*H*is the channel depth [6,8]. This open channel result applies equally to a closed channel with half-width

*H*. Using the approximation

*U*

_{c}/

*u*≈ 5 log

_{τ}_{10}

*Re*[27], this result can be written as

*K*≈ 2.72

*νRe*/ln(

*Re*). To compare results across the laminar–turbulent transition (see later in figure 12), we use equations (2.7) and (2.13) to obtain a non-dimensional turbulent diffusivity

*D*

_{e}=

*K*/

*D*

_{0}as a function of

*Pe*, where

*D*

_{0}is the characteristic diffusivity scale. However, we stress that the diffusivity depends only on

*Re*in the turbulent case: molecular diffusivity is negligible [5]. Note also from the result of Elder how the effective turbulent diffusivity grows less sharply with

*Pe*than in the laminar regime. Below, we will contrast the classical predictions for dispersion discussed in this section with our new results for biased swimmers.

### 3.2. Gyrotactic algae: a new, non-classical ‘swimming dispersion’

#### 3.2.1. Laminar downwelling flows: gyrotactic swimming strongly affects dispersion

It is illuminating to consider how gyrotactic cells distribute across a channel in downwelling laminar flows, as this determines their streamwise dispersion. Figure 4 displays the time evolution of the cross-sectional cell concentration profiles in the wall-normal direction *y* for a selection of values of (*Pe*,*β*). It is seen that an initially uniform concentration profile evolves to one focused at the channel centre. This is what one would expect for swimming gyrotactic cells with orientations biased by combined gravitational and viscous torques [11,12]. At the population scale, the long-time concentration distribution is a result of the balance between cross-stream cell diffusion and biased swimming. In §2.3, the cell concentration (normalized by its mean) is given by equation (2.21): (Recall subscripts *m* = G, F denote solutions of the GTD and FP microscopic models, respectively.) The two models predict a qualitatively different dependence of concentration distribution on *Pe* and *β*. While the FP approach allows for cells to diffuse more and focus less as the cells tumble in flows with large shear rates, the GTD model finds that with increasing local shear rate, *σ*, the diffusivity of tumbling cells decreases faster than the decrease in cross-stream cell focusing.

It has hitherto been unclear if GTD, mathematically more complicated than the FP approach, provides more accurate predictions. Thus, we compare the DNS results and theoretical predictions from the two models. Figure 5 displays the long-time cell distributions from DNS contrasted with theoretical predictions using GTD and FP. It is clear that DNS profiles get broader with increasing *β*/*Pe*, so profiles are shown for different values of this ratio (individual values of *Pe* and *β* are shown in the figure caption). GTD and DNS profiles agree very well for all the values of *β*/*Pe* simulated. The FP model only agrees well with the DNS for large values of *β*/*Pe*, where predictions coincide with those of GTD. The agreement is poor at low values (very poor for the sharply focused distribution *β*/*Pe* = 0.002). It is possible to quantify the scaling of the breadth of the profiles with the ratio *β*/*Pe*. From asymptotic expressions for the ratio [19], it follows that GTD profiles can be conveniently approximated by Gaussian distributions of width *y*_{0} = [(2/*λ*)(*β*/*Pe*)]^{0.5} (see appendix A). Clearly, these predictions need also to be tested experimentally (see §4), but the DNS results indicate that GTD is likely a superior model of gyrotactic response to flow than the FP approach.

As for the passive case of figure 3, we have quantified dispersion using the statistical measures: *Λ*_{0}(*t*), the drift above the mean flow; *D*_{e}(*t*), the effective diffusivity; and *γ*(*t*), the skewness. These are plotted in figure 6 for *Pe* = 27, 670 and 1675 at the fixed value of *β* = 33.5. In these simulations, statistically stationary values for *Λ*_{0}(*t*) and *D*_{e}(*t*) are achieved for dimensionless times *t* ∼ 1; steady values are reached earlier for larger *Pe*. In terms of the Bees & Croze dispersion theory, steady dispersion is achieved in the long-time limit when transient solutions to the moment equations have died down: . The analysis of transient solutions with DNS will be carried out in a future study, but it is reasonable that gyrotaxis makes the approach to steady dispersion faster than for passive tracers. The skewness is negative and approaches zero (with *γ* ∼ *t*^{−0.49} for *Pe* = 27; see the electronic supplementary material), suggesting a distribution tending to normality. Bees & Croze predicted the power-law decay *γ* = *γ*_{0}*t*^{−1/2} (the pre-factor *γ*_{0}(*β*, *Pe*) depends on gyrotactic swimming) with the same exponent as the passive case [7].

The steady gyrotactic swimmer dispersion displays some very surprising features, evident in the data presented in figure 6. For example, *Λ*_{0}, zero in the passive case, grows from a negative value to large positive values as *Pe* is increased. The effective diffusivity *D*_{e}, on the other hand, shows a non-monotonic behaviour for increasing *Pe* (recall *D*_{e} ∼ *Pe*^{2} for passive tracers). We can qualitatively account for this behaviour considering the concentration distributions analysed above. Cells are biased to swim to the centre of the channel. Here only the torque owing to gravity acts on cells, so cells swim upward at their mean swimming speed, which may be comparable with the mean flow speed for small *Pe* leading to *Λ*_{0} < 0. For large *Pe*, the upward swimming speed is negligible. As cell accumulation at the centre of the channel increases with *Pe* due to gyrotaxis they drift more relative to the mean flow, and so *Λ*_{0} is an increasing function of *Pe*. Cell accumulation at the centre of the channel removes them from regions of high shear rate, eventually leading to a decrease in their effective axial diffusivity with increasing *Pe*.

The Bees & Croze dispersion theory allows for the quantification of this intriguing dispersion behaviour; the results are summarized in §3.3. Prior to this, we shall consider time-dependent dispersion in turbulent flows.

#### 3.2.2. Turbulent downwelling flows: persistent but weaker swimming dispersion

Here, we describe the first DNS study of the dispersion of gyrotactic cells in turbulent channel flows. To compare turbulent and laminar results we first focus on downwelling flows. Simulations were performed for (*Pe*, *Re*) = (16 570, 2500), (28 140, 4200) and (67 000, 10 000), with corresponding values of *β* = 33.5, 55 and 67. The long-time wall-normal cell concentration profiles for these turbulent flows are shown in figure 7. It is clear that the mean cell concentration is uniform (barring small fluctuations about the mean) for almost the entire channel width, except for regions close to the wall that are gyrotactically depleted of cells (gyrotactically depleted regions occupy only a small fraction (less than 4%) of the channel width). Shear and advection experienced by a cell in turbulent flows are very different from the laminar case. The turbulent flow can be thought of as a time-averaged base profile on which are superposed turbulent fluctuations (the well-known Reynolds decomposition [27]). The shear rate of the base profile is close to zero in the middle of the channel, and large at the walls (figure 2) and deviates from the laminar case; this alone will lead to broader concentration profiles. On top of this, turbulent fluctuations perturb the flow, causing cell reorientation and advection. We can think of turbulence as endowing cells with an increased diffusivity [35] acting to make gyrotactic accumulations in this downwelling case less pronounced. Only close to the walls is the impact of mean shear sufficient to cause significant cell depletion; an effect that increases with *β*/*Pe* in a similar fashion to the laminar case. The small but measurable effects of gyrotaxis on the concentration profiles are reflected in the time-dependent dispersion measures for the turbulent case, plotted in figure 8 for the same values of (*Pe*, *β*) considered above. We leave a detailed analysis of transients to a future study but note that, because of the increased diffusivity by turbulence, the approach to the limiting behaviour is faster than in the laminar case. Less obviously, the long-time dispersion retains a rich behaviour (observe the order of the curves in figure 8). In particular, note that cells have a non-zero drift *Λ*_{0}; this is because of local focusing of cells in downwelling regions of the fluctuating flow. The dispersion of gyrotactic swimmers is thus qualitatively distinct from that of passive tracers even in a turbulent flow. As *Pe* is increased to the maximum value simulated, *Λ*_{0} decreases while *D*_{e} increases, indicative of increased mixing by turbulence.

#### 3.2.3. Dispersion in upwelling flows

The case of dispersion in flows directed vertically upwards (against the direction of gravity) is considered here briefly. The DNS results for the same values of *Pe* and *β* as the laminar downwelling case of §3.2.1 are shown in figure 9. The drift above the mean flow *Λ*_{0} is positive for small *Pe*, and grows more negative with increasing *Pe*. This behaviour is the result of the peculiar distribution of gyrotactic cells in the flow: cells in upwelling flow are biased to swim not to the channel centre, but to the walls [11]. Interestingly, accumulation and dispersion depend critically on the flow direction! The inset of figure 9*b* shows the cell concentration *c* in the wall-normal direction *y*, demonstrating the wall accumulation, which becomes more peaked with decreasing *β*/*Pe* (for more details, see the electronic supplementary material). Thus, for fixed *β* and increasing *Pe*, cells focus at the walls and experience slower flow and less of the shear profile, leading to a decrease in both the drift, *Λ*_{0}, and the diffusivity, *D*_{e}. The change of sign in *Λ*_{0} has a similar explanation as for the downwelling case.

The upwelling turbulent case, presented in figure 10, was also investigated for the same parameter values as the downwelling case of §3.2.2. We see that dispersion measures do not reach steady values: both *Λ*_{0} and *D*_{e} grow monotonically with time for the duration of the simulation run. The inset in figure 10*a* shows the cell concentration profiles displaying a thin band of gyrotactic accumulation at the walls. Cells that end up in this band are subject to strong dispersal by the large mean shear rate at the wall. Diffusion, whether owing to turbulence or swimming, is not strong enough to balance the shear-induced migration towards the wall, as is the case for laminar flow, leading to non-steady dispersion over the course of the simulations.

However, the upwelling dispersion results presented above are not realistic for swimming species that are negatively buoyant; accumulations of negatively buoyant cells at the walls can result in instability and the formation of downwelling flow and descending plumes near the walls [11]. These instabilities will differ in their extent in the laminar and turbulent upwelling cases, but we expect buoyancy-driven wall flows to ensue after accumulation in both regimes.

#### 3.2.4. The effects of cell elongation

So far we have assumed that the gyrotactic cells are spherical, setting the eccentricity parameter *α*_{0} to zero. Here, we present simulations obtained with a non-zero eccentricity, *α*_{0} = 0.8, for the cases (*Pe*, *β*) = (670, 6.7) (laminar) and (28 140, 55) (turbulent). We chose the relatively large value of *α*_{0} = 0.8 to emphasize the effect of eccentricity. Results for the effective diffusivity for downwelling flow are shown in figure 11 for the laminar and turbulent cases. The insets in the figure display the concentration profiles. We see that in the laminar case the effect of eccentricity is to broaden the profile, as observed in Bearon *et al.* [36]. This broader distribution results in an increased value of effective diffusivity (cells sample more of the shear profile). In the turbulent case of figure 11*b*, there is a much smaller broadening effect. The effect of cell elongation on other statistical measures for the parameters considered here is marginal (see the electronic supplementary material). If realistic values of biflagellate eccentricity (*α*_{0} = 0–0.3 [13]) are used, predictions for the dispersion of biased swimmers are not very different from those obtained here for spherical cells.

### 3.3. Long-time dispersion of gyrotactic cells as a function of *Pe* (*Re*) across the laminar–turbulent transition

Having analysed the full time dependence of gyrotactic cell dispersion, we concentrate here on its long-time behaviour. In this limit, laminar DNS results can be compared with predictions from analytical theory for drift *Λ*_{0} and diffusivity *D*_{e} as functions of *Pe* (for given *β*). The theory requires as inputs expressions for the mean swimming direction, **q**, and diffusivity tensor, **D**, obtained from microscopic stochastic models. We test two alternative microscopic models: FP and GTD. The models predict qualitatively different functional forms for the components of **D** as a function of the dimensionless shear *σ*. Correspondingly, predictions for *Λ*_{0}(*Pe*) and *D*_{e}(*Pe*) obtained with the two models also differ.

The predictions from the analytical theory are shown as curves in figure 12*a*,*b*, for fixed values of *β*; the corresponding DNS results for selected *Pe* and the same values of *β* are shown as points. It is clear that, for the GTD predictions, the agreement between theory and simulation for *Λ*_{0} is very good. The GTD prediction is that, for fixed *β*, *Λ*_{0} increases with *Pe*, tending to a linear behaviour for large enough *Pe*. This is in contrast to passive tracers, for which *Λ*_{0} = 0. For small *Pe* the FP prediction coincides with that of GTD, but, for larger *Pe*, *Λ*_{0} tends to a constant instead of growing with *Pe*. The smaller the value of *β*, the greater is the difference between the predictions of the two microscopic models (as for concentration profiles). A similar trend is seen for the effective diffusivity *D*_{e}, shown as a function of *Pe* in figure 12*b*. We see that for *β* = 33.5, both GTD and FP predict a diffusivity rising and then falling with *Pe*, consistent with DNS results. For smaller values of *β*, however, the FP and GTD predictions are different: GTD predicts a decrease of *D*_{e} at large *Pe*, while FP predicts a (second) rise. The DNS results are not compatible with this rise, but agree well with the GTD prediction. We remark that this dispersion behaviour is unique to swimming cells. It is due to gyrotactic bias and the ensuing accumulations that change the distribution of cells and thus their dispersion, as discussed above. Also shown in figure 12*b* is the diffusivity for passive tracers, which grows as *D*_{e} ∼ *Pe*^{2} [4], without the decrease at large *Pe* we predict for gyrotactic swimmers.

In earlier sections, we have seen how the effect of gyrotaxis on accumulation and dispersion is more subtle in turbulent flows. DNS results for long-time dispersion reflect these changes and are also summarized in figure 12. An analytical theory of turbulent swimming dispersion has yet to be formulated, so we compare only with theory for passive tracers. It is seen for all simulations with *β* = 33.5 that the value of *Λ*_{0}, which was growing with *Pe* in the laminar regime, drops a little just beyond the transition to turbulence. On the other hand, the diffusivity *D*_{e}, which was falling at large laminar *Pe*, suddenly acquires a sizeable value. The turbulent shear profile and fluctuations alter the distribution of cells. The effect of gyrotaxis is thus weakened but still measurable in the dispersion, which is similar but not identical to that of passive tracers. The fractional drift above the mean flow speed *U* = (2/3)*U*_{c} is given by *δ* ≡ lim_{t}_{→∞}(1/*U*)[*dm*_{1}/d*t* − *U*] = (3/2)*Λ*_{0}/*Pe*; this is rather small for large turbulent *Pe* values. Nevertheless, it is remarkable that drift should not be zero in a turbulent flow, as it is for passive tracers. For very large values of *Pe*, we expect gyrotaxis to have practically no effect on dispersion: the time-averaged concentration profile will be well approximated as uniform and *Λ*_{0} = 0. Indeed, our results suggest that *Λ*_{0} → 0 for increasing *Pe* in the turbulent regime.

## 4. Discussion

We have studied the dispersion of gyrotactic swimming algae in channel flows by DNS and analytical theory. This is the first study to evaluate cell distributions and statistical measures of gyrotactic cell dispersion (drift above mean flow *Λ*_{0}, effective diffusivity *D*_{e} and skewness *γ*) for flows on either side of the laminar–turbulence transition. We find unique cell accumulations and dispersion with non-zero drift and a non-monotonic diffusivity with *Pe*, with qualitative differences for upwelling and downwelling flows. The dispersion behaviour is remarkable and unique to gyrotactic swimming cells: passive tracers are transported at the mean flow rate (*Λ*_{0} = 0), the diffusivity increases with *Pe* (*D*_{e} ∼ *Pe*^{2}) [4], and dispersion is insensitive to channel orientation. In the laminar downwelling regime, simulation results were compared with the predictions derived from a recently formulated general analytical theory of swimming dispersion [18,19], applied here for the first time to channel geometry. The theory requires as inputs expressions for cell response to the local shear in a flow, determined from microscopic models.

We have evaluated predictions based on two alternative models: the FP [33,34] and GTD [19–21] approaches. Which model is more realistic has long been a matter of debate [20]. Here, we find that the DNSs are in excellent agreement with analytical predictions based on GTD, providing the first evidence that GTD is much better than the FP model at describing swimmers in flows. As well as validating analytical predictions for laminar flow, the DNSs allow us to explore the industrially relevant dispersion of gyrotactic algae in turbulent channel flows. In the turbulent regime, the effects of gyrtoaxis (accumulations resulting in non-zero drift) persist, but are much more subtle. Effective diffusivity in downwelling turbulent flows is similar to that of passive tracers. These are the first full DNSs of biased swimmers in turbulent channel flows; previous studies having concentrated on vortical flow [37] or synthetic approximations of homogeneous isotropic turbulence [31,35].

The fact that swimming algae in channel flows distribute very differently to passive tracers has important practical consequences for the culture of swimming species. The most dramatic implications of our findings are for photobioreactors that operate using laminar flows. For example, in draft tube air-lift bioreactors, bubbles up a central draft tube (riser) mix and aerate cells, which then circulate down the channel formed between the draft tube and the surface of the reactor (downcomer; figure 1*c*). The Bees & Croze [18] analytical theory, confirmed by simulations, predicts that gyrotactic swimming algae will be focused more and more sharply at the centre of the downcomer as *Pe* (the flow rate for fixed channel width, or *Re*) is increased. For example, considering a flow with *Pe* = 1675 (*Re* = 250) and the swimming Péclet number *β* = 33.5 (*H* = 5 cm), we predict that the concentration at the walls is a vanishingly small fraction of the mean, given by (using *λ* = 2.2 for *C. augustae*; see approximations in appendix A). Non-swimming cells and molecular solutes, on the other hand, would be uniformly distributed across the tube, As we can reasonably assume that the probability of cell adhesion to the walls is proportional to the concentration there, we predict that surface fouling by gyrotactic cells will be markedly reduced relative to non-swimming cells. Fouling can be a big problem in closed bioreactors because cell build-up can prevent light penetration and thus growth, and in extreme cases can clog reactor conduits [38].

The peculiar dispersion of biased swimmers will also affect growth in bioreactors. Our results indicate that gyrotactic swimmers in a downwelling flow drift faster and diffuse less than passive tracers, which travel at the mean flow. We can make experimentally measurable predictions for the dispersion of *C. augustae*, the alga whose swimming parameters (reorientation time *B*, rotational diffusivity *d*_{r}, swimming speed *V*_{s}) we have used in this study. For example, in the case with *β* = 6.7 and *Pe* = 670 (e.g. a realistic, small air-lift with *H* = 1 cm and *U*_{c} = 1 cm s^{−1}) the fractional drift above the mean flow, *δ* = (3/2)*Λ*_{0}/*Pe*, is estimated to be *δ* = 0.48. In other words, *C. augustae* cells drift about 50 per cent faster than passive tracers, such as nutrients. For the non-dimensional effective diffusivity we predict *D*_{e} = 0.8: the axial diffusivity of cells is smaller than the cell diffusivity scale () in the absence of flow. To compare these predictions with those for passive molecular solutes, we consider the dispersion of CO_{2} (molecular diffusivity *D _{t}* = 1.6 × 10

^{−}^{5}cm

^{2}s

^{−1}), a vital nutrient for algae. In the flow under consideration, cells have

*Pe*= 670, but CO

_{2}, with its smaller diffusivity, has Thus, the non-dimensional tracer diffusivity is , from the Taylor–Aris result. Re-dimensionalizing swimming and solute effective diffusivities by multiplying by diffusivity scales, we find that the axial diffusivity of CO

_{2}along a channel is approximately 10

^{5}times greater than that of

*C. augustae*. This dramatic differential dispersion of cells and nutrients could have important consequences for swimming cell growth in reactors.

In turbulent downwelling flows, which may also be realized in air-lifts, the effects of gyrotactic swimming are much less pronounced. Gyrotactic depletion is not as efficient as in laminar flows: the concentration at the walls is at best half of the mean concentration. We thus expect significant fouling in air-lifts under turbulent regimes. The influence of gyrotaxis on dispersion is also weaker than in the laminar case. For *C. augustae* in a turbulent channel flow with *β* = 55 (*H* = 8.4 cm) and *Pe* = 28 140 (*Re* = 4200, *U*_{c} = 5 cm s^{−1}), the DNSs predict a fractional drift of *δ* = 0.015. This is very small, and may be neglected for short channels. As shown in figure 12*b*, the effective diffusivity is of the same order as that of a passive tracer such as CO_{2}.

The excellent agreement between DNS and predictions for dispersion obtained using GTD theory provides a first confirmation that the swimming dispersion theory of Bees & Croze [18] is robust. If they are to be useful in bioreactor engineering design, however, the theoretical and numerical predictions for dispersion need to be tested against experiments. Work is in progress to test GTD predictions for pipe flow with the alga *Dunaliella* [23]. It would also be interesting to test experimentally the numerical and analytical predictions for channel flow presented herein, exploring in particular the laminar–turbulent transition region. Since tubes are often arranged horizontally in bioreactor designs, it will be interesting to study the effect of tube orientation on biased swimmer dispersion. Croze *et al.* [39] have carried out experiments on the dispersion of *C. augustae* in horizontal pipe flow for low *Pe*. They found a complex transport mediated by sheared bioconvection patterns and suggested that such cell-driven flows could alter the transition to turbulence. It would be interesting to study this transition experimentally and with the aid of simulations such as presented here, including the effects of cell buoyancy (especially for small *Pe* flows) [19].

More generally, our study could be extended to open channel flows, such as those present in pond bioreactors, channels, waterfalls and rivers. The swimming dispersion effects we have explored must also exist whenever biased swimmers are subject to shear in the ocean and lakes. Significantly, about 90 per cent of species implicated in the formation of harmful algal blooms swim using flagella [40]. A recent study has shown that *Heterosigma akashiwo* is gyrotactic and could be responsible for thin layer formation [15]. However, the role of biased swimming in the dispersion and ecology of algae remains largely unexplored.

## Acknowledgements

We thank Rachel Bearon for comments and discussions, and Miriam La Vecchia for help with implementing the formulation based on quaternions. Discussions with Victoria Adesanya, Stuart Scott and Alison Smith on photobioreactors are also acknowledged. Computer time provided by SNIC (Swedish National Infrastructure for Computing) is gratefully acknowledged. O.A.C. and M.A.B. are thankful for support from the Carnegie Trust and EPSRC (EP/D073398/1 and EP/J004847/1).

## Appendix A. Microscopic model solutions

**A.1. Numerical solutions for the Fokker–Planck and generalized Taylor dispersion models**

**q**(*σ*) and **D*** _{m}*(

*σ*) predicted by the FP and GTD microscopic models have the following structure: A1where we recall that

*m*= F represents solutions of the FP model and

*m*= G those of the GTD model. Such solutions can be obtained using approximations via a Galerkin method [19]. It is convenient to fit such solutions with the following expressions:

*q*(

^{x}*σ*) = −

*P*(

*σ*;

**a**

*,*

^{x}**b**

*);*

^{x}*q*(

^{y}*σ*) = −σ

*P*(σ;

**a**

*,*

^{y}**b**

*); Here, the rational function*

^{y}*P*(

*σ*;

**a**,

**b**) is given by A2and, for

*λ*= 1/(2

*d*

_{r}

*B*) = 2.2, the coefficients

**a**and

**b**are provided in table 2.

**A.2. Approximate generalized Taylor dispersion profiles and width scaling**

Using asymptotic results derived in [19], we derive an approximation for the concentration profiles predicted by the GTD model. These profiles are given by equation (2.21) as , where *m* = G for GTD predictions. For at leading order the GTD prediction asymptotes to , where *J*_{1} is a known constant for *λ* = 2.2. In the same limit, , so that For at leading order, *q ^{y}*(

*σ*) = −(2/3)

*λ*/

*σ*and , where

*d*

_{1}= 0.68 for

*λ*= 2.2. Thus, is a reasonable approximation for all

*σ*. Recalling for channel flow

*σ*= (

*Pe*/

*β*

^{2})

*y*, we see that Inserting this expression in the equation above and integrating gives the Gaussian profile A3where A4is the width of the profile. This scaling is observed in the concentration profiles obtained from DNS; see the main text.

- Received December 19, 2012.
- Accepted January 21, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.