2 Governing equations and numerical method

2.1 Fluid phase

The incompressible turbulent velocity field, $\boldsymbol{u}$ , obeys the Navier–Stokes and continuity equations:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$

Here $t$ denotes the time, $p$ the pressure, $Re$ the Reynolds number and $\boldsymbol{f}$ the random forcing field necessary to maintain the turbulent velocity in a statistically steady state. The equations are discretized in a cubic domain with periodic conditions at the boundaries. We employ direct numerical simulation (DNS) to solve all the relevant flow scales without any artificial model at the smallest scales. Since we evolve the equations in a three-periodic domain, it is natural to do so in Fourier space with a pseudo-spectral method. In this numerical scheme the nonlinear terms are evaluated in physical space using the classic $2/3$ rule to minimize aliasing error. Time integration is performed with a third-order, low-storage Runge–Kutta method where the diffusive terms are analytically calculated (note that in Fourier space the Laplacian operator is proportional to the square of the modulus of the wave vector $k^{2}$ ) while an Adam–Bashforth-like approximation is employed for the nonlinear terms (Rogallo Reference Rogallo1981). The stochastic forcing is evaluated in Fourier space and acts isotropically on the first shell of wave vectors. The forcing amplitude is constant in time and the field is delta-correlated in time and uniformly distributed in phase and directions (Vincent & Meneguzzi Reference Vincent and Meneguzzi1994). We employed a resolution of $128^{3}$ grid points with a Taylor Reynolds number $Re_{\unicode[STIX]{x1D706}}=100$ and $64^{3}$ grid points with a Taylor Reynolds number $Re_{\unicode[STIX]{x1D706}}=60$ . The turbulence is characterized by the Taylor Reynolds number $Re_{\unicode[STIX]{x1D706}}\equiv u^{\prime }\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ , where $u^{\prime }$ is the root-mean-square of velocity fluctuations, $\unicode[STIX]{x1D708}$ is the kinematic viscosity and $\unicode[STIX]{x1D706}$ is the Taylor microscale, $\unicode[STIX]{x1D706}\equiv \sqrt{\unicode[STIX]{x1D716}/15\unicode[STIX]{x1D708}u^{\prime 2}}$ , with $\unicode[STIX]{x1D716}$ the kinetic energy dissipation rate. The resolution was chosen so that the velocity gradients are well resolved since they drive the the equation of motion of the particles. Particle statistics are computed from independent samples separated by $1.6$ turbulent Kolmogorov time scales, $\unicode[STIX]{x1D70F}_{f}\equiv \sqrt{\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716}}$ , corresponding to about $300$ configurations per particle. The numerical code has been already employed to study gyrotactic micro-organisms in homogeneous turbulence (Zhan et al. Reference Zhan, Sardina, Lushi and Brandt2014) and rain droplet evaporation in clouds (Sardina et al. Reference Sardina, Picano, Brandt and Caballero2015).

2.2 Solution of the particle motion

The prolate spheroids studied here are small compared to the Kolmogorov length scale of the turbulence and have densities close to that of the carrier fluid. The translational relaxation time, $\unicode[STIX]{x1D70F}_{p}$ , for prolate spheroids, assuming an isotropic particle orientation distribution is given by

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{p}=\frac{2\unicode[STIX]{x1D70C}_{p}a^{2}}{9\unicode[STIX]{x1D707}_{f}}\frac{{\mathcal{A}}{\mathcal{R}}\ln ({\mathcal{A}}{\mathcal{R}}+\sqrt{{\mathcal{A}}{\mathcal{R}}^{2}-1})}{\sqrt{{\mathcal{A}}{\mathcal{R}}^{2}-1}},\end{eqnarray}$$

with $a$ the semi-minor axis of the prolate spheroid and ${\mathcal{A}}{\mathcal{R}}$ the aspect ratio (the ratio of polar over equatorial radii of the spheroid) (Shapiro & Goldenberg Reference Shapiro and Goldenberg1993). The time scale for the particles studied here is around three to four orders of magnitude smaller than the relaxation time of the turbulence flow around the particles. Furthermore, rotational relaxation times, concerning spinning and tumbling motions, are reported in the literature to be typically shorter than the translational relaxation time (Zhao et al. Reference Zhao, Challabotla, Andersson and Variano2015; Marchioli, Zhao & Andersson Reference Marchioli, Zhao and Andersson2016). With these operating parameters a tracer model for the particle motion is a suitable assumption (Shin & Koch Reference Shin and Koch2005; Voth & Soldati Reference Voth and Soldati2017) and we, therefore, assume that the cells behave like passive tracers with a correction given by their Stokes settling velocity.

The particle translational motion obeys

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{x}}{\text{d}t}=\boldsymbol{u}|_{x}+\boldsymbol{v}_{s}(\boldsymbol{p}), & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}_{s}\,(\boldsymbol{p})=v_{s}^{min}\hat{\boldsymbol{e}}_{g}+\left(v_{s}^{max}-v_{s}^{min}\right)\left(\hat{\boldsymbol{e}}_{g}\boldsymbol{\cdot }\boldsymbol{p}\right)\boldsymbol{p}, & \displaystyle\end{eqnarray}$$

with $\boldsymbol{x}$ denoting spheroid position, $\boldsymbol{u}|_{x}$ the fluid velocity at spheroid position, $\boldsymbol{p}$ the spheroid orientation vector, $\boldsymbol{v}_{s}$ the Stokes settling velocity and $\hat{\boldsymbol{e}}_{g}$ a unit vector in the direction of gravity. Here $v_{s}^{min}$ and $v_{s}^{max}$ are the minimum and maximum settling speeds in quiescent flow, corresponding to the particles falling with their major axis perpendicular or parallel to the gravity direction. These velocities can be written as

(2.6a,b ) $$\begin{eqnarray}v_{s}^{max}=\frac{\left(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f}\right)gl^{2}}{\unicode[STIX]{x1D707}_{f}}\unicode[STIX]{x1D6FE}_{0}\left({\mathcal{A}}{\mathcal{R}}\right);\quad v_{s}^{min}=\frac{\left(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f}\right)gl^{2}}{\unicode[STIX]{x1D707}_{f}}\unicode[STIX]{x1D6FE}_{1}\left({\mathcal{A}}{\mathcal{R}}\right),\end{eqnarray}$$

with $g$ the gravitational acceleration, $l$ the length of the spheroid’s major axis, $\unicode[STIX]{x1D707}_{f}$ the dynamic viscosity of the fluid, $\unicode[STIX]{x1D70C}_{p}$ and $\unicode[STIX]{x1D70C}_{f}$ the densities of the particle and the fluid. Here $\unicode[STIX]{x1D6FE}_{0}$ and $\unicode[STIX]{x1D6FE}_{1}$ are functions of the aspect ratio ${\mathcal{A}}{\mathcal{R}}$ , whose full expression can be found in Dahlkild (Reference Dahlkild2011) among others.

The spheroid rotation follows the inertialess Jeffery equation

(2.7) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{p}}{\text{d}t}=\frac{1}{2}\unicode[STIX]{x1D74E}|_{x}\times \boldsymbol{p}+\unicode[STIX]{x1D6FC}\left[\unicode[STIX]{x1D644}-\boldsymbol{p}\boldsymbol{p}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D640}|_{x}\boldsymbol{\cdot }\boldsymbol{p},\end{eqnarray}$$

with $\unicode[STIX]{x1D74E}|_{x}$ the flow vorticity at the particle’s position, $\unicode[STIX]{x1D644}$ the second-order identity tensor, $\unicode[STIX]{x1D640}|_{x}$ the symmetric part of the velocity deformation tensor and $\unicode[STIX]{x1D6FC}$ a function of the spheroid aspect ratio ${\mathcal{A}}{\mathcal{R}}$ , defined as

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\frac{{\mathcal{A}}{\mathcal{R}}^{2}-1}{{\mathcal{A}}{\mathcal{R}}^{2}+1}.\end{eqnarray}$$

To solve the equations above, the flow velocity, vorticity and velocity gradients at the particle positions are calculated with a second-order interpolation scheme and a high enough resolution to capture the flow velocity gradients. The same Runge–Kutta scheme used for the carrier phase is employed to perform the time integration of (2.4) and (2.7) with a small enough time step ( ${\approx}0.001\unicode[STIX]{x1D70F}_{f}$ ) to capture the time evolution of particle orientations via a formulation based on quaternions. We evolve a total of $200\,000$ spheroids for each simulation.

Table 1.

Maximum, minimum and average of settling velocities in quiescent flow, scaled with the Kolmogorov velocity for the investigated cases versus different aspect ratios. Here $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1.05$ and $D_{eq}=\unicode[STIX]{x1D702}/6$ .

4 Final remarks

We report results from simulations of elongated, rigid particles (prolate spheroids) with different aspect ratios sinking in homogenous, isotropic turbulence. The parameters used in this study represent marine or freshwater diatoms that are characterized by very small relaxation times compared to that of the fluid but with settling velocities on the order of Kolmogorov velocities characteristic of the mixed layer in the ocean. Given the small size of diatoms relative to the Kolmogorov length-scale and their short relaxation time, cells are assumed to behave like passive tracers with a correction given by their Stokes settling velocity. We show that turbulence does increase the settling velocities of inertia-less elongated particles, compared to their Stokes settling velocities, and that the effect increases as the aspect ratio of the particle increases, but saturates at ${\mathcal{A}}{\mathcal{R}}\geqslant 10$ . The magnitude of this effect, however, is small ( ${\sim}5\,\%$ , figure 2 $c$ ). It should be noted that despite the enhanced mean settling velocity for diatoms reported in this study, there could be dispersion of particles in turbulent flow that increases the residence time of some individuals.

Further, we study homogenous isotropic turbulence at the dissipation scale, and thus neglect large-scale turbulence that might play a role in suspending diatoms. Indeed, the question ‘How local enhancement is related to residence time in the mixed layer of the ocean?’ is still unanswered and might be very well related to the large scale turbulence that is not captured in the current model. In our model, turbulence dissipation rate, $\unicode[STIX]{x1D716}$ , was set to $10^{-6}~\text{m}^{2}~\text{s}^{-3}$ . This choice is within the range observed in the ocean’s surface boundary layer in a coastal environment (Gerbi et al. Reference Gerbi, Trowbridge, Terray, Plueddemann and Kukulka2009). For this turbulence level, settling is slightly enhanced. Ruiz et al. (Reference Ruiz, Macías and Peters2004) used a Couette device to generate turbulence and showed that sinking speeds of different species of phytoplankton increased significantly (by an order of magnitude for some species) as $\unicode[STIX]{x1D716}$ increased from $\unicode[STIX]{x1D716}<10^{-6}~\text{m}^{2}~\text{s}^{-3}$ to $\unicode[STIX]{x1D716}>10^{-4}~\text{m}^{2}~\text{s}^{-3}$ . They also showed that downward velocities increased as the ratio between the size of the cell (described as equivalent spherical diameter in their study) and the Kolmogorov length-scale increased (Ruiz et al. Reference Ruiz, Macías and Peters2004).

Here we show that non-random distributions can arise within the parameter space relevant to diatoms, with weak clustering occurring only for non-spherical cells. Preferential sampling of downwelling flow is due to the interactions between shear, particle rotation, drag, and lift. Elongated spheroids tend to align with the local flow strain (more frequently with the second eigen-direction) and with the local flow vorticity, in agreement with previous studies on neutrally buoyant spheroids in homogenous isotropic turbulence (e.g. Ni et al. (Reference Ni, Kramel, Ouellette and Voth2015), Voth & Soldati (Reference Voth and Soldati2017) and references therein). By considering the simple case of one particle sedimenting in linear shear flow we demonstrate that the time-averaged horizontal position of a particle in the flow (over a complete Jeffery orbit) is biased towards downwelling regions. Spheroids with higher average Stokes settling speeds ( $v_{s}^{iso}$ ) show greater tendency to cluster (see figure 4 $a$ ). However, when $v_{s}^{iso}$ is greater than the Kolmogorov velocity, the increase in clustering does not correspond to an increase in relative average settling speed. We also show that clustering increases significantly with drag anisotropy (i.e. the ratio $v_{s}^{max}/v_{s}^{min}$ ) in a way that does not appear to saturate. Drag anisotropy increases with increasing ${\mathcal{A}}{\mathcal{R}}$ and forces particles to sample regions of downwelling flow to a further extent, even when the isolated effect of aspect ratio on the rotational motion is saturated.

Mechanisms that affect sinking and dispersion of phytoplankton hold significant implications to carbon and energy fluxes in aquatic systems. Underlying distributions of phytoplankton at the micro-scale are not well understood as we still lack the tools to measure them in the field. Few available observations suggest that at scales of centimeters, the seascape is more heterogeneous than previously thought (Malkiel, Alquaddoomi & Katz Reference Malkiel, Alquaddoomi and Katz1999; Seuront, Schmitt & Lagadeuc Reference Seuront, Schmitt and Lagadeuc2001; Mitchell et al. Reference Mitchell, Yamazaki, Seuront, Wolk and Li2008) but mechanisms remain poorly understood. This study provides insights into one mechanism by which weak clustering of chains and elongated cells of diatoms can occur at the millimetre to centimetre scale. A typical fractal dimension of $2.8{-}2.9$ shows a weak clustering, and it might be that other mechanisms exist behind the clustering that our model cannot capture.

Cell shape is an important morphological trait affecting many aspects in the ecology of phytoplankton (Smayda Reference Smayda1970; Lewis Reference Lewis1976; Sournia Reference Sournia1982; Reynolds Reference Reynolds1989; Reynolds & Irish Reference Reynolds and Irish1997; Karp-Boss & Boss Reference Karp-Boss and Boss2016) and results from our numerical model further highlight the complex interactions of cell shape with ambient flows. Interestingly, the observed increase in settling speeds does not increase indefinitely, and saturates for ${\mathcal{A}}{\mathcal{R}}>10$ . Increasing the aspect ratio, while keeping the particle volume and density ratio constant, varies other characteristics of the particle simultaneously. For example, the mean settling speed in quiescent flow decreases as the ratio between maximum (when particle is oriented parallel to the gravity direction) and minimum (when particle is oriented normal to the gravity direction) settling speed increases. We have therefore performed a series of numerical experiments to isolate the effect of the aspect ratio ${\mathcal{A}}{\mathcal{R}}$ , the magnitude of settling velocity $v_{s}^{iso}$ and the ratio between the maximum and minimum settling speed $v_{s}^{max}/v_{s}^{min}$ in quiescent flow. When $v_{s}^{iso}$ and $v_{s}^{max}/v_{s}^{min}$ are held constant, the increase in clustering and settling speed begins to saturates at ${\mathcal{A}}{\mathcal{R}}=5$ . This indicates that the relative enhancement that we observe for ellipsoids with aspect ratios larger than $5$ is due to an increase in the ratio between the maximum and minimum settling speed ( $v_{s}^{max}/v_{s}^{min}$ ) that is associated with a change in particle shape ( ${\mathcal{A}}{\mathcal{R}}$ at constant volume). Interestingly, an analysis of typical aspect ratios among phytoplankton shows a peak in distribution at ${\mathcal{A}}{\mathcal{R}}\sim 5$ (Clavano, Boss & Karp-Boss Reference Clavano, Boss and Karp-Boss2007; Karp-Boss & Boss Reference Karp-Boss and Boss2016), but this may be a coincidence because in nature multiple processes (grazing, diffusion, mechanical breakage by shear) act simultaneously as selective pressures on shape and size in phytoplankton (Karp-Boss & Boss Reference Karp-Boss and Boss2016).

The model presented here is an idealized case of turbulence and phytoplankton in the ocean. Direct numerical simulations are computationally expensive and cannot be used to capture the full range of scales characteristic of oceanic turbulence, and are limited to low Reynolds numbers. We consider homogeneous isotropic turbulence, but in the ocean most turbulent mixing occurs near boundaries and is subjected to stratification and vertical shear. It was found, however, that in stratified regions of the ocean, isotropic motions do exist over some range of scale bounded at the lower end by the Kolmogorov length scale, providing that the non-dimensional isotropy parameter $I=\unicode[STIX]{x1D716}/\unicode[STIX]{x1D708}N^{2}$ exceeds ${\sim}200$ , where $\unicode[STIX]{x1D716}$ is the local mean kinetic energy dissipation, $\unicode[STIX]{x1D708}$ the kinematic viscosity and $N$ the local mean value of the buoyancy frequency (Thorpe Reference Thorpe2007).

Modelling phytoplankton as rigid prolate spheroids with an even distribution of mass is another simplification of the model. Diatoms exhibit complex morphologies (see figure 1) and often bear spines and other cell projections that can alter dynamics in the flow. Another attribute of chain-forming diatoms that is dynamically important is flexural stiffness, as it spans four orders of magnitude between different species (Young Reference Young2011). Variations in rotation rates between rigid and flexible particles were observed in a laboratory setting (Karp-Boss & Jumars Reference Karp-Boss and Jumars1998), and it is likely that rigid and flexible chains will also exhibit different dynamics in turbulence, a subject that requires more investigations in the future.