Clustering of gyrotactic microorganisms in turbulent flows

Clustering of gyrotactic microorganisms in turbulent flows


We study the spatial distribution of gyrotactic microorganisms transported by a three-dimensional turbulent flow generated by direct numerical simulations. We find that gyrotaxis combines with turbulent fluctuations to produce small scales (multi-)fractal clustering. We explain this result by showing that gyrotactic swimming cells behave like tracers in a compressible flow. The effective compressibility is derived in the limits of fluid acceleration much larger and smaller than the gravity.

05.45.-a, 47.63.Gd, 92.20.jf

Microbial patchiness in oceans is important for ecological and evolutionary dynamics Levin and Whitfield (1994); Azam and Malfatti (2007) and for biogeochemical processes Falkowski et al. (2000). In motile aquatic microorganisms, self-propulsion provides a mechanism to escape fluid pathlines, potentially leading to small-scale patchiness Mitchell et al. (1990); Torney and Neufeld (2007); Durham et al. (2011). Remarkably, motility combined with fluid flows can also generate large-scale inhomogeneities. For instance, spectacular aggregation of phytoplankton cells (in layers centimeters to meters thin, horizontally extending from hundreds of meters to kilometers) can result from vertical shears and gyrotactic swimming Durham et al. (2009). Gyrotaxis characterizes several species of motile microalgae whose swimming direction is determined by the balance of viscous and gravitational torques, due to the displacement between the cell center of mass and buoyancy. As an effect of such balance, for example, gyrotactic algae aggregate in the center (wall) of descending (ascending) vertical pipe flows Kessler (1985); Pedley and Kessler (1987). Gyrotaxis is observed in algae, e.g., of the genus Chlamydomonas, which can be engineered to transport microloads Weibel et al. (2005), or Dunaliella, employed in biofuels Chisti (2007). So far most studies focused on the dynamics of gyrotactic microorganisms in simple stationary flows or kinematic models Kessler (1985); Pedley and Kessler (1987); Hill and Bees (2002); Durham et al. (2009); Thorn and Bearon (2010); Durham et al. (2011).

In this Letter, we investigate the interplay between gyrotactic motility and realistic turbulent flows, as occurring in the sea. We find that turbulence and gyrotaxis combine to generate inhomogeneous distributions with small-scale (multi-)fractal statistics (see Fig. 1). We study the limit of gravitational acceleration much smaller or larger than turbulent accelerations to identify the mechanisms responsible for gyrotactic clustering in terms of an effective compressible velocity field.

Figure 1: (color online) Spatial distribution of gyrotactic swimmers (dots) in a slab of a 3D turbulent flow. Color code: yellow/blue corresponds to high/low vorticity values (). (Left) Limit of orientation dominated by local fluid acceleration (, see text) with aggregation in high vorticity regions. (Right) Limit of gravity dominated orientation (). Parameters correspond to circled symbols in Fig. 2c and Fig. 3a, respectively.

We consider dilute suspensions of non interacting motile microorganisms, much smaller than the smallest scale of turbulence, the Kolmogorov length . We can thus model them as self-propelled particles with velocity,


given by the sum of the fluid velocity at the particle position and the swimming contribution , where the swimming speed is assumed constant Pedley and Kessler (1987); Torney and Neufeld (2007). Cells are assumed spherical and neutrally buoyant, with the center of mass displaced by with respect to the geometric one. The swimming direction , determined by the total torque acting on the cell, evolves as


where is the fluid vorticity and is the orientation speed for spherical cells subject to the acceleration Pedley and Kessler (1987). In a fluid at rest, besides viscous forces, only gravity (and buoyancy) is acting and thus , while acceleration due to swimming is neglected Pedley and Kessler (1987). In presence of a flow, we have where


is the fluid acceleration given by the Navier-Stokes equations ruling the velocity of an incompressible () fluid with viscosity , pressure and stirred by an external forcing . Previous studies on gyrotactic swimmers disregarded fluid acceleration, as mainly focused on simple, non-turbulent flows where . In turbulence, fluid acceleration can locally exceed La Porta et al. (2001) and therefore its contribution has to be taken into account.

The first term on the rhs of Eq. (2) causes the direction of swimming to align with on a time scale . When the contribution of fluid acceleration can be neglected, cells tend to orient vertically () on a time scale . The alignment is contrasted by the vorticity term and, depending on being smaller or larger than , cells may swim along a resulting local equilibrium direction or tumble randomly as the orientation becomes unstable due to vorticity Pedley and Kessler (1987); Thorn and Bearon (2010). In principle, the swimming direction may be modified also by rotational Brownian motion Berg (1993) and tumbling due to flagella desynchronization during swimming Polin et al. (2009), which are here neglected. The former effect is very small for typical algae (having size ); the latter can be neglected whenever the tumbling time is longer than the reorientation one.

We study gyrotactic swimming in homogeneous and isotropic turbulent velocity fields of moderate intensity () by means of direct numerical simulations of Navier-Stokes equations. In particular, Eq. (3) is solved by means of a standard pseudospectral algorithm with 2 order Runge-Kutta time-stepping, on a tri-periodic cubic grid of size (for and ). Statistical stationarity is guaranteed by means of a zero mean, Gaussian and white in time random forcing restricted to large scales. Viscosity is such that the Kolmogorov length is of the order of the grid spacing, ensuring well resolved small-scale velocity dynamics. For different values of , several populations of swimmers, characterized by different values of and are injected with random positions and orientations. At each time step, velocity and acceleration at the swimmers positions, needed to integrate Eqs. (1-2), are obtained by interpolation. The self-propelled particles are then evolved, and their distribution and orientation studied in statistically steady conditions. In the sequel, we mostly focus on the dependence on the orientation speed by fixing , being the typical fluid velocity fluctuation at the Kolmogorov scale.

Figure 2: Swimmer properties as a function of the orientation parameter () in the limit . (a) Average alignment with fluid acceleration . (b) Average square vorticity at swimmers position normalized to the volume average value. (c) correlation dimension . Circled symbol in (c) corresponds to the distribution shown in Fig. 1a.

Formally, Eqs. (1-2) define a dissipative dynamical system evolving in the -dimensional (actually because and ) phase space (,) with phase-space contraction rate


As orients in the direction , is expected to be negative on average, meaning that swimmers will evolve onto a dynamical attractor of dimension smaller than the whole phase space, which explains why clustering can be observed: if the fractal dimension of the attractor is smaller than , clustering in position space (as in Fig. 1) is possible (see Ref. Bec (2005) for a conceptually similar phenomenon occurring for inertial particles). We remark that clustering is a consequence of swimming: indeed for Eqs. (1) and (2) decouple, thus cells become tracers advected by an incompressible velocity and cannot cluster. Moreover, in the limit we have and therefore swimmers cannot cluster. Nonetheless, even in this limit, if they deviate from fluid trajectories and generate interesting dynamics Khurana et al. (2011).

We now discuss the physical mechanisms of clustering which, as anticipated, depend on whether the dominating effect comes from the gravitational () or fluid acceleration (which we quantify in terms of its rms value ).

We start considering the case and therefore we take in Eq. (2). Figure 2 summarizes the behavior of the main observables as a function of the dimensionless number (now ) measuring the ratio of the alignment timescale to rotation timescale induced by vorticity. When the alignment is very fast, the swimming direction becomes parallel to the local direction of the fluid acceleration , as confirmed by Fig. 2a showing that for (here and in the following denotes average over particle distribution). In this limit, swimming cells behave like tracers advected by an effective velocity . While is incompressible, the effective velocity field is not: being negative (positive) in high vorticity (strain) regions. Therefore, as it occurs for inertial particles lighter than fluid Balkovsky et al. (2001); Calzavarini et al. (2008), the swimmers cluster inside vortical structures (Fig. 1a and Fig. 2b). The divergence of is proportional to , clustering is thus expected to increase with the swimming speed. In the opposite limit of slow alignment, when , random tumbling due to fluid vorticity dominates, hence swimming orientation cannot align to the local acceleration (, see Fig. 2a): the compressible effect is lost and particles distribute uniformly in the volume. To quantify clustering we measured the correlation dimension , ruling the small-distance () behavior of the probability to find two swimmers at separation less than : Paladin and Vulpiani (1987). For uniformly distributed particles , while when clustering is present the probability to find close pairs increases and (see e.g. Bec et al. (2007) for a similar study in the case of inertial particles). In Fig. 2c we show as a function of : for , , indicating strong clustering in almost filamental structures; conversely, when , the correlation dimension approaches the uniform-distribution value .

Figure 3: (color online) Clustering properties as a function of , for (). (a) Correlation dimension of the swimmer positions. Circled symbol corresponds to the data shown in Fig. 1b. (b) Variances of swimming direction components ( in red, and ). The dashed blue curve is the parabola . The solid horizontal line represents the random orientation value .
Figure 4: (color online) (a) Average fluid acceleration (positive, red filled circles) and velocity (negative, blue open circles) along the vertical at swimmer positions, expressed in percentage of and , respectively. (b) Correlation vs . The maximum for is understood noticing that must decrease for . In the limit , (solid line) is implied by , see text. In the random tumbling limit (Fig. 3a), which implies (dashed line).

We now consider the limit when we can take and Eq. (2) reads


with . Similarly to the previous case, when the cells orient in the preferred direction , which is now fixed in space. The effective velocity thus becomes which, unlike the previous case, is incompressible (). Therefore, now we expect that not only for but also for swimmers distribute uniformly, as confirmed by Fig. 3a showing that in both limits. Remarkably, Fig. 3a shows that also in this case gyrotactic swimmers cluster on a fractal set (see Fig. 1b) for intermediate values, with a well defined minimum of the correlation dimension () for . We remark than an optimal orientation timescale for aggregation is also observed in steady kinematic vortical flows Durham et al. (2011) where, however, a vast class of trajectories is integrable.

We can understand the origin of the observed clustering by considering the limit . In such limit, cell orientation being very fast we can assume that the swimming direction is always at an equilibrium orientation with (see Fig. 3b). In particular, solving Eq. (5) for , at first order in , one finds and (which is confirmed by simulations). As a consequence, the effective swimmer velocity field with has a compressible component with divergence


which, unlike the previous case, is unrelated to fluid acceleration so that swimmers will cluster in regions different from those of high vorticity (compare Fig. 1a and b). We notice that (6) generalizes the well known mechanism of cell focusing in the center (walls) of downward (upward) vertical pipe flows Kessler (1985). Notice that in the above argument the vertical component of the vorticity plays no role, as it does not change .

Another consequence of the expansion is that (resp. ) and () have locally the same (opposite) sign. Numerical simulations show that this remains true also for larger values of , on average. Indeed, at stationarity, by averaging Eq. (5) and using isotropy on the plane (guaranteed by the isotropy of the fluid velocity field) we obtain . The correlation between the horizontal components of and implies that the swimmers will stay longer in regions of the flow characterized by positive vertical velocity and negative vertical acceleration (Fig. 4a). This can be easily seen in a case with, say, a vortex aligned with the x-axis, where the above argument with implies , , so that the trajectories spend more time in regions where , as there the swimming velocity opposes that of the fluid. The preferential concentration in these regions of the flow will be maximal (and correspondingly the correlation dimension minimal, i.e. clustering stronger) for where the correlation between swimming direction and vorticity is also maximal (Fig. 4b). For such value of fluid regions with maximal deviation of the swimming direction from the vertical will balance vorticity dominated ones where tumbles randomly. We observe that the swimmer vertical migration can be strongly inhibited by the bias towards downwelling regions: in Fig. 4a, e.g., can reach of the swimming speed .

Figure 5: (color online) (a) Average square vorticity at swimmers position normalized to the volume average value at varying the ratio ( corresponds to data of Fig. 3); (b) correlation dimension vs ; (c) vs swimming speed , the circled symbol corresponds to the circled one in (b); (d) generalized dimensions vs for circled data in (b), notice that the case (filled black circles) appears to be less multifractal than when also the fluid acceleration is contributing to clustering (empty blue circles).

In the general case, the relative importance of fluid and gravitational accelerations for clustering depends on the ratio . Figure 5a indeed shows that the bias towards regions of high vorticity decreases with and is absent when only the gravitational torque is acting (). The correlation dimension , shown in Fig. 5b, smoothly varies with , interpolating from the two limits shown in Fig. 2c and 3b. We observe that, as anticipated, clustering is more effective for large swimming speeds as displayed in Fig. 5c, showing that, at fixed value of , decreases with . Finally, as one can expect from general considerations on dynamical attractors Bec (2005), Fig. 5d demonstrates that the spatial distribution of the gyrotactic self-propelled particles is multifractal, as the generalized dimensions (controlling the probability to find particles at small separation) depends on the moment Paladin and Vulpiani (1987).

Summarizing, we have shown that gyrotactic motility and realistic turbulent flows can generate small-scale patchiness (down to the Kolmogorov scale) in the distribution of bottom-heavy swimming microorganisms. We identified two mechanisms driving microorganism clustering: the focusing in vortical regions due to local adjustment of the swimming orientation with fluid acceleration, and the correlation between vorticity and swimming direction on the plane perpendicular to gravity leading particles to preferentially explore downwelling, upward accelerating regions. In general, gravity is expected to dominate when turbulent intensity is not very high and it is likely the most important effect in the ocean. Crucial parameters for observing clustering are in this case the ratio between swimming speed and small-scale fluid velocity fluctuations () and the reorientation time scale with respect to vorticity intensity (). For typical microalgae and Jones et al. (1994); Hill and Häder (1997); Kessler (1985). In the ocean, the turbulence intensity, measured in terms of kinetic energy dissipation , varies from in the upper mixing layer down to a few meters deeper Yamazaki and Squires (1996); Yamazaki et al. (2002). We can thus estimate that and therefore the effects discussed in this Letter are relevant in realistic conditions and can definitely be tested in laboratory by tuning turbulence characteristics.

We conclude by remarking that for non-spherical cells such as, e.g., prolate spheroids the term should be added to Eq. (2) ( being the eccentricity, and and the symmetric rate of strain tensor and identity matrix, resp.) Pedley and Kessler (1987). Such term is also contributing to the phase-space contraction rate (4) providing an additional mechanism for clustering Torney and Neufeld (2007). It will thus be interesting to study if and how gyrotactic clustering in turbulence is modified at varying the cell shape.

We thanks S. Musacchio for useful discussions. GB and MC acknowledge KITPC institute for hospitality during the program New Directions in Turbulence and support by MIUR PRIN-2009PYYZM5 “Fluttuazioni: dai sistemi macroscopici alle nanoscale”.


  1. S. Levin and M. Whitfield, Philos. T. Roy. Soc. B 343, 99 (1994).
  2. F. Azam and F. Malfatti, Nat. Rev. Microbiol. 5, 782 (2007).
  3. P. Falkowski et al., Science 290, 291 (2000).
  4. J. Mitchell, A. Okubo,  and J. Fuhrman, Limn. Ocean. 35, 123 (1990).
  5. C. Torney and Z. Neufeld, Phys. Rev. Lett. 99, 078101 (2007).
  6. W. Durham, E. Climent,  and R. Stocker, Phys. Rev. Lett. 106, 238102 (2011).
  7. W. M. Durham, J. O. Kessler,  and R. Stocker, Science 323, 1067 (2009).
  8. J. O. Kessler, Nature 313, 218 (1985).
  9. T. J. Pedley and J. O. Kessler, Proc. Royal Soc. B 231, 47 (1987); Annu. Rev. Fluid Mech. 24, 313 (1992).
  10. D. Weibel et al., Proc. Natl. Acad. Sci. USA 102, 11963 (2005).
  11. Y. Chisti, Biotech. Advances 25, 294 (2007).
  12. N. A. Hill and M. A. Bees, Phys. Fluids 14, 2598 (2002).
  13. G. J. Thorn and R. N. Bearon, Phys. Fluids 22, 041902 (2010).
  14. A. La Porta et al., Nature 409, 1017 (2001).
  15. H. Berg, Random walks in biology (Princeton Univ Pr, 1993).
  16. M. Polin et al., Science 325, 487 (2009).
  17. J. Bec, J. Fluid Mech. 528, 255 (2005).
  18. N. Khurana, J. Blawzdziewicz,  and N. T. Ouellette, Phys. Rev. Lett. 106, 198104 (2011).
  19. E. Balkovsky, G. Falkovich,  and A. Fouxon, Phys. Rev. Lett. 86, 2790 (2001).
  20. E. Calzavarini et al., Phys. Rev. Lett. 101, 84504 (2008).
  21. G. Paladin and A. Vulpiani, Phys. Rep. 156, 147 (1987).
  22. J. Bec et al., Phys. Rev. Lett. 98, 84502 (2007).
  23. M. S. Jones, L. Le Baron,  and T. J. Pedley, J. Fluid Mech. 281, 137 (1994).
  24. N. Hill and D. Häder, J. Theor. Biol. 186, 503 (1997).
  25. H. Yamazaki and K. Squires, Mar. Ecol. Prog. Ser. 144, 299 (1996).
  26. H. Yamazaki, D. L. Mackas,  and K. L. Denman, in The sea: biological-physical interaction in the ocean, Vol. 12 (John Wiley & Sons, New York, 2002) p. 51.
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
Request answer
The feedback must be of minumum 40 characters
Add comment
Loading ...