{'fraction_non_alphanumeric': 0.03841347355468825, 'fraction_numerical': 0.013650433545459788, 'mean_word_length': 4.56180496150556, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dispersion of swimming algae in laminar and turbulent channel flows: consequences for photobioreactors

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 the analytical swimming dispersion theory of Bees and Croze (2010). Time-resolved dispersion measures are evaluated as functions of the Peclet 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 Peclet number. The cell effective axial diffusivity increases and decreases with Peclet number (for tracers it merely increases). In turbulent flows, gyrotactic effects are weaker, but discernable and manifested as non-zero drift. These results should significantly impact photobioreactor design.


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 [2; 3]. The dispersion of passive tracers in fluid flows is well understood, particularly in the industrially relevant cases of flows in pipes and channels [4]. G. I. Taylor first realised 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 [5]. A similar result also holds for turbulent pipe [6] and channel [7] 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 [8; 9]. Particles whose density differs from the suspending medium exhibit more complex behaviour (for example, they can accumulate at walls in turbulent channel flows; see [10; 11] and references therein).
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 [12; 13]. This bias, known as gyrotaxis, results from the combination of viscous torques on the cell body, due to flow gradients, and gravitational torques, arising from bottom-heaviness and sedimentation [14]. 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 [15; 16; 17]. Gravitaxis and gyrotaxis can also couple in a complex fashion to other biases due to chemical gradients (chemotaxis) and light (phototaxis [18]). 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 nontrivial.
In a recent study, Bees & Croze extend the classical Taylor-Aris analysis of dispersion in a laminar flow in a tube to the case of biased swimming algae [1]. 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 model [13], Bees & Croze 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) [1]. 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 generalised 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 due to cell locomotion within local a flow, as well as the gyrotactic orientational biases considered in the Fokker-Planck 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 to obtain greater concentrations per reactor area [24]. A common design is the raceway pond, in which rotating paddles generate flow (see figure 1a). 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), 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, see figure 1b-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, whilst flat panel and columnar/annular set-ups operate at 1-10 cm/s [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 Article submitted to Royal Society 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 (Re 100), while raceway ponds and tubular reactors are always turbulent (Re> 30000). The mixing properties of turbulent flows are thought 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 like 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 biotechnological interesting species like Dunaliella, should be very different from that of passive tracers [1; 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, comparing results to predictions from the Bees & Croze dispersion theory applied to channel geometry. The paper is structured as follows. In Section 2, we present the simulation methods and outline a derivation of the theory for the new geometry. In Section 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 Section 4, we interpret the results, and discuss their implications for dispersion in photobioreactors and the environment. We simulate the dynamics of a population of biased swimming micro-organisms, typically N = 2 × 10 5 and 10 6 individuals per simulation (see table 1), placed in laminar and turbulent flows. Each microswimmer is modelled as a spheroidal particle whose position x i , i = 1..N , evolves according to

Methods
where V s is the mean cell swimming speed, u is the local fluid velocity and p is 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 [29; 13]

2)
Article submitted to Royal Society where k is a unit vector in the vertical direction; B = µα ⊥ /(2hρg) is the gyrotactic reorientation time (h denotes the centre-of-mass offset relative to the centre-ofbuoyancy, α ⊥ 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 √ 2d r , 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 where D Dt ≡ ∂ ∂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 [1]).

(b) DNS: 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 2H. When focusing on swimming cells it would seem natural to rescale length by the plate half-width H and time by H 2 d r /V 2 s , 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 V 2 s d −1 r . However, the parameter values obtained from this 'cellbased' 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 timescale H/U c , the time taken for a flow with centerline 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 where starred quantities denote dimensionless variables. In particular, we define the dimensionless swimming speed v s = V s /U c , the gyrotactic parameter η = 2BU c /H and the Reynolds number where U c is the centerline speed, H is the channel half-width and ν the fluid's kinematic viscosity. Non-dimensionalising the noise Γ r also defines the dimensionless rotational diffusivity d * r = d r H/U c . 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.5) and (2.6) numerically (see supplementary materials) to find Var(x * ) 3/2 , where Pe is the cell Péclet number defined with respect to the center-line 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.
(c) 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 [1] continuum dispersion theory for biased swimmers will be applied to the new channel geometry; it is valid in the long-time limit (t ≫ H 2 d r /V 2 s ). The derivation for channel flow is similar to the pipe flow example in Bees & Croze, so we provide an outline here. Those readers that 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.10) The 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 where 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-dimensionalise lengthscales by H and timescales by the time to diffuse across the channel H 2 d r /V 2 s . Thus,x = (x−U t)/H,ŷ = y/H, z = z/H andt = V 2 s t/(H 2 d r ), where hats denote dimensionless variables.
We assume unidirectional coupling of the cell dynamics to the flow (for bidirectional coupling due to non-neutrally buoyant cells see [1]): 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, ω = ∇ × u = −χ ′ e z .
Consider planar Poiseuille flow, χ = (1 − 3y 2 )/2 and U = (2/3)U c [32], where U c is the centerline (maximum) flow speed. Equation (2.10) becomes where 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 re-write β = HV 2 s d r /V 2 s and interpret it as a 'swimming' Péclet number, the ratio of the rates of transport by swimming and diffusion. Noflow and no-flux boundary conditions will be applied to (2.12), such that u = 0 and n · (D · ∇n − βqn) = 0, on Σ, (2.14) where 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 [8]. The pth moment with respect to the axial direction is defined as provided x p n(x, y, t) → 0 as x → ±∞. We denote cross-sectional averages by overbars. The cross-sectionally averaged axial moment is thus In Cartesian coordinates equation (2.12) becomes Multiplying by x p and integrating over the length of an infinite channel we obtain the moment evolution equation Averaging over the cross-section and applying the no-flux boundary conditions (2.14), yields from which we calculate measures of cell dispersion. The drift above the mean flow, Λ 0 ≡ lim t→∞ d dt m 1 (t), is given by is the zeroth axial moment (normalised concentration profile). Similarly the effective diffusivity, D e ≡ lim t→∞ The weighting function Y 0 0 (y) 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 Fokker-Planck (FP) model [33; 34] and generalized Taylor dispersion (GTD) [20; 21; 19]. For spherical cells, each model predicts how nondimensional swimming direction q and diffusivity D depend on two nondimensional quantities: , the ratio of reorientation time by rotational diffusion to that by shear (vorticity), and λ = 1/(2d r B), the ratio between the time of reorientation by rotational diffusion (1/2)d −1 r 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. 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) using q(σ) and D m (σ) from the two models, see supplementary materials.

(d ) Parameters, simulation time and flow profiles
In simulations and theoretical evaluations we used the following cell parameters based on C. 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) [1]. Further, 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 Since χ ′ (y) = −3y for plane Poiseuille flow, the maximum dimensionless shear is given by , the nondimensional bias parameter. Simulations are more readily interpreted and compared to 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 selfsimilar, so the nondimensional flow has the same parabolic profile for all Re. Timeaverged 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= 10000.

(a) Passive tracers: classical dispersion
The dispersion of passive tracers, such as molecular dyes or nonmotile 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 [32; 4]. As a benchmark test, we carried out direct numerical simulations for passive tracers (solving equation (2.4) with v s = 0 and translational noise to simulate molecular diffusivity, see methods). 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 = 76191, compares very well with simulation results. The skewness slowly tends to zero at long times as t −0.5 (see supplementary material), suggesting approach to a normal distribution [8]. The above results for passive    . The inset shows the expected zero drift, Λ0, above the mean flow and the skewness, γ, tending to zero at long times, suggesting an approach to the expected normally distributed axial concentration profile [8].
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 xx = 1/6 = D yy , q i = 0 and D xy = 0 (as Y 0 0 = 1 and 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 V 2 s /d r . Using the scaling (1/3)V 2 s /(2d r ) for random diffusers in 3D, we recover D ′ e = 1 + (8/945)Pe ′2 , the classical Taylor-Aris result for channels.
The dispersion of passive tracers in turbulent channel flows has also been elucidated [9]. 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 [7; 9]. 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 nondimensional 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 [6]. Notice 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. 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 [12; 13]. At the population scale, the long-time concentration distribution is a result of the balance between cross-stream cell diffusion and biased swimming. In section (2c), the cell concentration (normalised by its mean) is given by equation (2.21)): Y 0 0 = c(y)/c = exp β y 0 (q y /D yy m )ds . (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 β. Whilst 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    β/Pe. From asymptotic expressions for the ratio q y /D yy G [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 discussion), 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: t ≫ τ 1 . 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 supplementary materials), 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 [8].
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 due 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 upwards 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 section (3c) below. Prior to this, we shall consider time-dependent dispersion in turbulent flows. 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 (< 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, see 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, due to 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 7). In particular, notice that cells have a nonzero drift Λ 0 ; this is due to 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 whilst D e increases, indicative of increased mixing by turbulence.

(iii) 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 section (3i) 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 [12]. Interestingly, accumulation and dispersion depend critically on the flow direction! The inset of figure 9, bottom, shows the normalised cell concentration c/c in the wall normal direction y, demonstrating the Article submitted to Royal Society wall accumulation, which becomes more peaked with decreasing β/Pe (for more details, see supplementary materials). 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 section (3ii). 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 to the top figure 10 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 due 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 [12]. 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. (iv) 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 (28140, 55) (turbulent). We chose the relatively large value of α 0 = 0.8 to emphasise 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 (bottom) there is a much smaller broadening effect. The effect of cell elongation on other statistical measures for the parameters considered here is marginal, see supplementary materials. If realistic values of biflagellate eccentricity (α 0 = 0-0.3 [14]) are used, predictions for the dispersion of biased swimmers are not very different from those obtained here for spherical cells.

(c) Long-time dispersion of gyrotactic cells as a function of P e (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 Article submitted to Royal Society 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: Fokker-Planck (FP) and generalised Taylor dispersion (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 12a, 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  Simulation results compare well with GTD predictions for all β, but, for Pe 1, they are incompatible with FP predictions for small β. Also plotted in (b) are the classical results for passive tracers. In laminar flows, Λ0 = 0 and De ∼Pe 2 , in sharp contrast with our predictions for swimmers. In turbulent flow, we still have Λ0 > 0 for swimmers, though it appears to decay to zero with increasing Pe. The turbulent De for swimmers is close to the passive case; see text.
as a function of Pe in 12b. 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 12b is the diffusivity for passive tracers, which grows as D e ∼ Pe 2 [5], 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 summarised 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 /dt − 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.

Discussion
We have studied the dispersion of gyrotactic swimming algae in channel flows by direct numerical simulations (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 Peclet number, 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 ) [5], 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 [1; 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 Fokker-Planck (FP) [33; 34] and generalised Taylor dispersion (GTD) [20; 21; 19] approaches. Which model is more realistic has long been a matter of debate [20]. Here we find that the DNS are in excellent agreement with analytical predictions based on GTD, providing the first evidence that it is much better than the FP model at describing swimmers in flows. As well as validating analytical predictions for laminar flow, the DNS 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 direct numerical simulations of biased swimmers in turbulent channel flows; previous studies having concentrated on vortical flow [37] or synthetic approximations of homogeneous isotropic turbulence [35; 31].
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; see figure 1c). The Bees & Croze [1] analytical theory, confirmed by simulations, predicts that gyrotactic swimming algae will be focused more and more sharply at the center 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 Peclet number β = 33.5 (H= 5cm), we predict that the concentration at the walls is a vanishingly small fraction of the mean, given by c/c ≈ exp [−Peλ/(2β)] ∼ 10 −22 (using λ = 2.2 for Chlamydomonas augustae; see approximations in Appendix A). Non-swimming cells and molecular solutes, on the other hand, would be uniformly distributed across the tube, c/c = 1. 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 cells buildup 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 Chlamydomonas 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) 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% faster than passive tracers, such as nutrients. For the nondimensional effective diffusivity we predict D e = 0.8: the axial diffusivity of cells is smaller than the cell diffusivity scale (V 2 s /d r = 1.49 × 10 −3 cm 2 /s) 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), a vital nutrient for algae. In the flow under consideration, cells have Pe= 670, but CO 2 , with its smaller diffusivity, has Pe t =Pe V 2 s /(d r D t ) = 6.2 × 10 4 . Thus, the non-dimensional tracer diffusivity is D e,t ≈ (8/945)Pe 2 t ∼ 10 7 , from the Taylor-Aris result. Re-dimensionalising swimming and solute effective diffusivities by multiplying by diffusivity scales, we find that the axial diffusivity of CO 2 along a channel is ≈ 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 realised in air-lifts, the effects 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.4cm) and Pe= 28140 (Re= 4200, U c = 5cm/s), the DNS predict a fractional drift of δ = 0.015. This is very small, and may be neglected for short channels. As shown in figure 12b, 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 generalised Taylor dispersion theory, provides a first confirmation that the swimming dispersion theory of Bees & Croze [1] 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 Chlamydomonas 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% 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 [16]. However, the role of biased swimming in the dispersion and ecology of algae remains largely unexplored.  Table 2. Parameters used in the functional fits of solutions to the FP and GTD models used for our analytical theory predictions (from [19]).
Appendix A. Microscopic model solutions (a) Numerical solutions for the Fokker-Planck and generalised Taylor dispersion models q(σ) and D m (σ) predicted by the FP and GTD microscopic models have the following structure: where we recall 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 y ); D xx m (σ) = P (σ; a xx , b xx ); D yy m (σ) = P (σ; a yy , b yy );D xy m (σ) = −σP (σ; a xy , b xy ). Here the rational function P (σ; a, b) is given by and, for λ = 1/(2d r B) = 2.2, the coefficients a and b are provided in table (2) below.
(b) Approximate GTD 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 c(y)/c = exp β y 0 (q y /D yy m )ds , where m = G for GTD predictions. For σ ≪ 1 at leading order the GTD prediction asymptotes to q y (σ) = −σJ 1 /λ, where J 1 is a known constant for λ = 2.2. In the same limit, D yy G (σ) = J 1 /λ 2 , so that (q y /D yy G ) σ≪1 = −σλ. For σ ≫ 1 at leading order, q y (σ) = −(2/3)λ/σ and D yy m (σ) = d 1 /σ 2 , where d 1 = 0.68 for λ = 2.2. Thus q y /D yy G ≈ −σλ is a reasonable approximation for all σ. Recalling for channel flow σ = (Pe/β 2 )y, we see that βq y /D yy G ≈ λ(Pe/β)y. Inserting this expression in the equation above and integrating gives the Gaussian profile where is the width of the profile. This scaling is observed in the concentration profiles obtained from DNS; see main text.

S1. Numerical methods for the simulations
The numerical code is an efficient pseudo-spectral solver for the three-dimensional incompressible Navier-Stokes equations (equations 2.3 or 2.6 in the main text) with a particle tracking algorithm [11] adapted to swimming micro-organisms. The velocity components of the fluid phase are expanded in both x (streamwise) and z (spanwise) direction with Fourier modes, whereas Chebyshev polynomials are employed in the wall-normal y-direction. To advance the Navier-Stokes equations in time, we use a fourth order Runge-Kutta algorithm. Periodic boundary conditions are assumed in x and z, with no-slip at both walls y = ±1. More details about the code are given in [41]. For the DNS with laminar flow, the flow is not simulated, but is taken to be the analytic parabolic Poiseuille solution to the Navier-Stokes equations .
The micro-organisms are treated as point particles evolved by means of a Lagrangian solver. The fluid velocities and the components of the velocity gradient tensor are interpolated from the Eulerian grid onto the particle positions using a tri-linear interpolation. The time advancement of the particle uses the same Runge-Kutta algorithm as the time-advancement of the fluid except for the noise which is integrated by an Euler-Marayuma method [42]. Equation (2.5) of the main text is solved using quaternions rather than angles [43]. We perform numerical simulations of turbulent plane Poiseuille flow at constant mass flux. The Reynolds number, domain size and resolution used for the results presented here are reported in table 3. Note that resolution and box size are chosen following the classic results for turbulent channel flow in [44; 45].
In the table, we report both the large-scale Reynolds number defined by Re = U c H/ν where U c is the centerline streamwise velocity for the laminar flow of same mass flux, and the turbulent Reynolds number Re τ = u τ H/ν. The latter is defined by the friction velocity u τ = √ σ w and the half-channel width, where σ w ≡ ν[∂U/∂y] wall is the shear stress at the wall [27].
Note that the simulation domain is shorter than the distance travelled by the swimmers during dispersion; also, for large times, the area occupied by the cells becomes longer than our computational box. However, the computational domain can be considered long enough for the velocity correlation to be negligible at a distance of the order of the box length L x . Therefore, a very long domain for the swimmers can be assembled by means of copies of a single Eulerian computational box of length L x by using the method of repeated domains described in [10].

S2. Evalutation of the analytical dispersion measures
To evaluate the integrals in (2.20) and (2.22), it was simpler to numerically evaluate an equivalent system of steady differential equations (e.g. use the numerical solution to dY 0 0 /dy = β(q y /D yy )Y 0 0 , instead of (2.21)). The system of ODEs was solved using the Matlab bvp4c routine. Solutions to the system were found for using both the Fokker-Planck and GTD expressions given in Appendix A. For comparison with the analytical predictions, concentration profiles from DNS were normalised by dividing by the area under the profile (the number of cells), evaluated using the Matlab trapz routine.   [1] who calculated that the skewness for swimmers should approach zero as ∼ t −0.5 at long times, as it does for passive tracers [8].    Figure 15. Effect of cell elongation on dispersion measures for downwelling flows: (a) laminar (b) turbulent. Even for the unrealistically large value of eccentricity used (α0 = 0.8), the effect of elongation is small apart from the diffusivity in the laminar case, as discussed in the main text.