# Emergence of phytoplankton patchiness at small scales in mild turbulence

^{a}Department of Dynamics of Complex Fluids, Max Planck Institute for Dynamics and Self-Organization, 37077 Göttingen, Germany;^{b}Max Planck Institute for Dynamics and Self-Organization, 37077 Göttingen, Germany;^{c}Interdisciplinary Centre for Mathematical Modelling, Loughborough University, Loughborough, Leicestershire LE11 3TU, United Kingdom;^{d}Department of Mathematical Sciences, Loughborough University, Loughborough, Leicestershire LE11 3TU, United Kingdom

See allHide authors and affiliations

Edited by David R. Nelson, Harvard University, Cambridge, MA, and approved October 3, 2018 (received for review May 21, 2018)

## Significance

Recent observations found that swimming phytoplankton species have a patchy spatial distribution down to the millimeter scale in oceans and lakes. This is rather surprising because it goes against the intuitive expectation that the turbulent flow mixes the microorganisms and spreads them rather homogeneously. Here, we identify the relevant scales that rationalize the observation of patchiness. We find that a mild turbulent field can lead to the formation of small dense patches of motile cells when the typical scales of the turbulence roughly match the scales of the interaction among microorganisms, which is mediated by hydrodynamics.

## Abstract

Phytoplankton often encounter turbulence in their habitat. As most toxic phytoplankton species are motile, resolving the interplay of motility and turbulence has fundamental repercussions on our understanding of their own ecology and of the entire ecosystems they inhabit. The spatial distribution of motile phytoplankton cells exhibits patchiness at distances of decimeter to millimeter scales for numerous species with different motility strategies. The explanation of this general phenomenon remains challenging. Furthermore, hydrodynamic cell–cell interactions, which grow more relevant as the density in the patches increases, have been so far ignored. Here, we combine particle simulations and continuum theory to study the emergence of patchiness in motile microorganisms in three dimensions. By addressing the combined effects of motility, cell–cell interaction, and turbulent flow conditions, we uncover a general mechanism: The coupling of cell–cell interactions to the turbulent dynamics favors the formation of dense patches. Identification of the important length and time scales, independent from the motility mode, allows us to elucidate a general physical mechanism underpinning the emergence of patchiness. Our results shed light on the dynamical characteristics necessary for the formation of patchiness and complement current efforts to unravel planktonic ecological interactions.

The spatial distribution of plankton is key to understanding myriad aspects of marine ecology: from the flux of nutrients in the water column, or predation strategies (1, 2), to the global repercussions of harmful algal blooms (3). It is one of the oldest observations in biological oceanography that phytoplankton exhibit a heterogeneous spatial distribution (patchiness) at length scales of kilometers (4, 5). Recent technological advances have ushered in measurements at small scales (decimeters to millimeters) that revealed patchiness also at these smaller scales (6). These small distances form the natural arena where biological and physical processes such as nutrient uptake, growth, and community composition take place (2, 7⇓⇓–10). Due to various factors such as wind or tidal currents, the local, small-scale environment of phytoplankton is nearly invariably turbulent (11), and turbulence has a profound impact on the biology and ecology of plankton (12). For example, turbulent diffusion reduces the size of the concentration boundary layer of depleted nutrients and thus facilitates the nutrient uptake of the microorganism (13⇓–15).

Numerous phytoplankton species are motile, and the vast majority of harmful algal blooms are caused by motile phytoplankton (3, 11). Thus, a comprehension of the consequences of active motion in a turbulent environment can clarify important aspects of the phytoplankton ecology, such as division rates and nutrient retrieval. Recent work has established that gyrotactic plankton swimming in mildly turbulent environments form dense patches (16⇓⇓–19). However, patchiness is exhibited by a myriad of motile species with different motility modes, such as phototaxis, chemotaxis, pattern swimming, and autoregulated aggregation (3, 20⇓–22). Furthermore, flow visualization techniques have demonstrated that motile plankton create microflows which in turn modify the fluid environment (2) and will introduce effective hydrodynamic interactions among individual cells at the scale of about

Here, we combine molecular dynamics simulations and a continuous description of the incompressible Navier–Stokes flow to determine a general mechanism for patchiness in motile phytoplankton cells swimming in a turbulent environment. Unlike previous investigations of gyrotactic plankton, we include a minimal model of cell–cell interactions among the swimmers. Our work explores the fundamental spatial and temporal scales associated to the reorientation of the swim direction when cells encounter each other and the scales associated to the mild turbulent motion. We identify dimensionless numbers that compare the characteristic time scales associated to the small-scale motion of the fluid with the cell–cell interactions and that compare the characteristic length scales associated to the turbulent motion with the active motion of the cells.

Associating time and length scales to turbulent dynamics is a familiar operation from classical studies of turbulence, but the role of these scales in the context of phytoplankton with hydrodynamic coupling is not known. We consider two complementary methods to generate a turbulent flow: (*i*) a simplified but realistic representation of turbulent flow via random Fourier modes (23, 24), the so-called kinematic simulation method first proposed by Kraichnan (25), and (*ii*) direct numerical simulations (DNS) of the Navier–Stokes equations with a pseudospectral method. Our investigations are carried out in a parameter range close to the ecologically relevant Taylor-based Reynolds number (26)

We find that the correlation length scale associated to the vorticity of the flow constitutes a characteristic scale, and when it matches the length scale associated to the cell–cell interactions, strong patchiness emerges. Experimental evidence at such a level of precision is scarce. However, the current evidence points to the fact that minuscule changes in the cellular morphology can dramatically affect the hydrodynamic cell–cell interactions (27) and how different organisms respond to turbulent microflows (28), thus giving further motivation for our work.

## Results

### Modeling Phytoplankton Dynamics in Mildly Turbulent Environments.

Consider a motile microorganism, such as phytoplankton, swimming in a turbulent aqueous milieu. A fundamental length scale characterizing the scale where the energy associated with turbulent motion dissipates is the Kolmogorov length scale

Because of their sizes, the motion of microorganisms in aquatic environments is dominated by viscous forces, which are two to five orders of magnitude larger than the inertial forces (32). This means that the microorganisms’ motion and the small-scale perturbations to the fluid produced by their motility are described by the Stokes equation *Chlamydomonas* (34, 35). Beyond this length scale ϵ the stochastic nature of the flow prevails. Short- to medium-range cell–cell interactions with up–down symmetry are the key ingredient for our modeling.

We model the cells as particles with swimming speed *i*) *ii*) a torque representing the short-range cell–cell interactions with up–down symmetry. The interaction is restricted to particles within the range ϵ and takes place on a time scale *Materials and Methods* for more details.

To model the turbulent flow, we combine results from two complementary techniques to generate a turbulent flow. In the first approach, we use kinematic simulations; that is, the turbulent velocity field is modeled by a sum of unsteady random Fourier modes which obey a prescribed energy spectrum. We refer to this model as the “Kraichnan fluid” after its first proponent (25). In the second approach, we perform DNS of the incompressible Navier–Stokes equations with a pseudospectral method to model the turbulent flow directly (more details in *SI Appendix*). (We note that the Kolmogorov scale as a measure for the small-scale turbulent motions is meaningful only in the context of 3D Navier–Stokes flow. In that sense, comparisons to real-world data are more sensible for this setting than for simplified 2D flows.)

The orientational dynamics can be characterized by two dimensionless numbers. The first one is the rotational Péclet number

### Phase Diagram of Patchiness.

In our simulations, we observe phytoplankton patchiness as a consequence of the cell motility and the cell–cell interaction for a wide range of parameters. This can be reasoned as follows: The dynamics described by the equations of motion of the swimming cells (Eq. **2** in *Materials and Methods*) produce a compressing phase-space volume

For a more detailed picture of the mechanism we map out the nonequilibrium phase diagram of cell patchiness, identifying the relevant length and time scales. To quantify the degree of patchiness in the spatial distribution of swimming cells we calculate the patchiness factor Q, which is a normalized local density based on the Voronoi tessellation of the 3D domain (16) (more details in *Materials and Methods*). Fig. 1*A* shows the nonequilibrium phase diagram obtained from our molecular dynamics simulation of swimming cells immersed in a Kraichnan fluid. At fixed Péclet number, as the vortical Stokes number

The relevant length scales associated to phytoplankton patchiness warrant the investigation of the impact of the cells’ finite size—an aspect so far largely neglected by models of motile particles. The finite extension of the cells induces correlations in their spatial distribution, and these correlations will affect the characteristics of the patches. To study whether our results are robust upon inclusion of more realistic spatial correlations, we perform simulations of extended particles immersed in a Kraichnan fluid. Fig. 1*B* (circles) shows that patchiness emerges also for extended cells for

We can derive some limits on the location of the maximum for patchiness in the following way. First, for the cell–cell interaction to overcome the rotational diffusion, the Péclet number P must be larger than unity (this condition corresponds to the region above the dotted line in Fig. 1*A*). Second, to have an effect on the dynamics, the vorticity associated to the turbulent flow should be stronger than the rotational diffusion, which is quantified by the condition *A*). This characterization of patchiness in terms of time scales can be complemented by considering the relevant length scales. To this end, we define a dimensionless number *A*. In the current setting, we find that the condition

Kinematic simulations provide a synthetic flow field with realistic spatiotemporal correlations. However, they lack a number of hallmark features of real hydrodynamic turbulence such as a self-consistent energy spectrum with realistic large- and small-scale cutoffs and non-Gaussian velocity and vorticity fluctuations. Although the degree of non-Gaussianity is low for the mild turbulence considered here, departures could have a sensitive impact on the clustering behavior. We therefore turn to full DNS calculations (*SI Appendix*) to further corroborate our findings and compare them with the results of the Kraichnan fluid. Fig. 1*B* (squares) shows the patchiness factor Q for the same range of

This coincidence corroborates our argument above about the mechanism of patch formation: The cell–cell interactions lead to a patchy distribution of swimming cells provided that the scales associated to turbulent motion and cell–cell interactions match. The height and precise position of the maximum in Q depend on the details of the fluid flow. We find that both the Taylor-based Reynolds number and the ratio of cell–cell interaction length scale to Kolmogorov length scale influence the height of the maximum in Q (details in *SI Appendix*).

Fig. 2 shows a typical snapshot of the system. The visualization of vorticity fields of the Navier–Stokes flow exhibits the filamentary vortex structures, which are characteristic of turbulent flows. The position and orientation of the phytoplankton cells are indicated with arrows. It appears that the locations of the phytoplankton patches are subtly influenced by the joint effect of active motion and flow vorticity. While there are some events of considerable reorientation of the cells within the vortex cores, in a more typical situation cells are mildly deflected by the local vorticity.

Finally, to bolster our discussion of the influence of the correlation length of the turbulent field compared with the interaction range ϵ, we perform simulations of the Kraichnan flow at fixed

Fig. 3 shows the dependence of the maximum in patchiness factor Q on the ratio of

## Discussion

Noninertial particles in an ergodic incompressible flow can cluster only if they are motile because they can cross the streamlines of the flow. Passively advected particles which initially are homogeneously distributed will be advected following volume-conserving dynamics. We find that small-scale patchiness emerges when hydrodynamic cell–cell interactions couple with the vorticity in a turbulent flow field.

The following physical picture emerges from our results. When the reorientational strength of the cell–cell interaction is stronger than the rotational diffusion of the fluid, the cell–cell interaction leads to (some degree of) alignment of neighboring cells. The turbulent fluctuations may then statistically bring two motile cells closer together. From this point on the cells will experience a similar physical environment and will be very likely to remain at close distance.

When a balance is struck between the length scales where reorientational interactions and turbulence act, marked by

The results of Fig. 1, together with Fig. 3, also explain the maximum in clustering observed upon increasing

Let us now put our findings in the context of experimental observations. Recent technological advances have expanded the range of length scales we can access and consequently also our view of the problem of plankton patchiness. In situ fluorometry (36, 37), rapid freezing (38, 39), and pneumatically operated sampling (40, 41) allow one to reach submillimeter scales (6) in measurements of planktonic distributions.

Experiments using high-resolution gradient multisamplers can give access to spatial resolutions of the order of decimeters which correspond to the scale of patchiness predicted by our results. Measurements at the pycnocline in the Kattegat found an increase of a factor of 20 with respect to the surface layers and a variation of a factor of 3–5 within 30 cm (7) in *Gyrodinium aureolum*. *Dinobryon* species concentrations were observed to vary by up to a factor of 10 in a river estuary within a meter scale (42). Variations in concentration of factors up to 2 or 3 below the meter scale were observed for different species of dinoflagellates in the Aarhus Bay (8). These numbers are in agreement with the predictions of our theory (Fig. 1*B*).

A combination of traditional turbulence sensors [e.g., shear probes and fast response thermistors (43, 44)] and a backscatter fluorometer measured microstructure in the velocity and in the intensity of chlorophyll concentrations at the centimeter scale (45), so-called biophysical microstructures. As high-frequency acoustic scattering techniques grow in precision (46, 47), additional tools can be deployed to map marine and freshwater environments to increasing accuracy, which may facilitate future experimental confirmations of our findings.

Our results ignore the presence of stratification. While this is an important feature of real-world settings, we do not expect that our conclusions will change qualitatively, as Q will result in a weighted average of the individual strata.

## Conclusion

Aiming to identify a minimal model for motile phytoplankton, we investigated the collective behavior of a large number of both point-like and extended, motile cells immersed in a mildly turbulent flow field. We quantify the influence of cell–cell interactions via a dimensionless number, the vortical Stokes number

We explore the microhydrodynamics characterized by a Taylor-based Reynolds number of the order of

We confirmed that both models (Kraichnan flow and Navier–Stokes flow) show a near quantitative agreement. We find, however, an influence of the kinematic details, such as the ratio *SI Appendix*). This means that the occurrence of small-scale patches strongly depends on the specific turbulent field which acts on the cells.

This finding can complement discussions on the influence of physical forces on marine ecology. Although for smaller phytoplankton species it is conceivable that the local environment is a linear shear field, the microflows change rapidly on the scale of

By addressing the combined effect of turbulent flow and cell–cell interactions, our work paves the way to unravel the complex interplay of physical and biological interactions determining phytoplankton’s microenvironment.

## Materials and Methods

### Molecular Dynamics Simulations.

Our system is composed of N identical, finite-size particles which are self-propelled and have an up–down symmetric interaction with each other. These self-propelled particles inhabit a cubic box with side length L and periodic boundary conditions. We integrate the equations of motion for particle i at position

All simulations are performed in a cubic simulation domain of side length ranging *SI Appendix* for discussion on influence of

We perform two sets of simulations: In the first one we use

### Analysis of Patchiness.

To quantify the degree of clustering in a given configuration, we use the patchiness factor defined in ref. 16. To derive this factor, the Voronoi tessellation of a given configuration is calculated, which attaches a volume

## Acknowledgments

The authors thank Jens Elgeti, Stephan Herminghaus, and Anupam Sengupta for useful discussions. The authors gratefully acknowledge support of the Max Planck Society.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: marco.mazza{at}ds.mpg.de.

Author contributions: M.G.M. designed research; R.E.B., C.C.L., D.W., M.W., and M.G.M. performed research; R.E.B. and C.C.L. analyzed data; and R.E.B., C.C.L., M.W., and M.G.M. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1808711115/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵
- ↵
- Smayda TJ

- ↵
- ↵
- ↵
- Franklin RB,
- Mills AL

- Pinel-Alloul B,
- Ghadouani A

- ↵
- ↵
- ↵
- ↵
- Robinson AR,
- McCarthy JJ,
- Rothschild BJ

- Yamazaki H,
- Mackas DL,
- Denman KL

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Santamaria F,
- De Lillo F,
- Cencini M,
- Boffetta G

- ↵
- Gustavsson K,
- Berglund F,
- Jonsson PR,
- Mehlig B

- ↵
- ↵
- ↵
- Torney C,
- Neufeld Z

*Phys Rev Lett*, 101:078105. - ↵
- Mariani P,
- MacKenzie BR,
- Visser AW,
- Botte V

- ↵
- ↵
- ↵
- Tennekes H,
- Lumley JL

- ↵
- ↵
- ↵
- ↵
- Anis A,
- Moum JN

- ↵
- ↵
- ↵
- Baskaran A,
- Cristina Marchetti M

- ↵
- Drescher K,
- Dunkel J,
- Cisneros LH,
- Ganguly S,
- Goldstein RE

- ↵
- ↵
- ↵
- Ghadouani A,
- Smith REH

- ↵
- Krembs C,
- Juhll AR,
- Long RA,
- Azam F

- ↵
- Krembs C,
- Juhl AR,
- Rudi Strickler J

- ↵
- ↵
- Lunven M, et al.

- ↵
- Viličić D, et al.

- ↵
- Lueck RG,
- Wolk F,
- Yamazaki H

- ↵
- Nash JD,
- Moum JN

- ↵
- ↵
- DeGrâce DL,
- Ross T

- ↵
- ↵
- Maar M,
- Gissel Nielsen T,
- Stips A,
- Visser AW

- ↵
- Lebwohl PA,
- Lasher G

- ↵
- Breier RE,
- Selinger RLB,
- Ciccotti G,
- Herminghaus S,
- Mazza MG

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences