1 Introduction
Understanding the dynamical and statistical properties of self-propelled agents (such as motile microorganisms, microrobots or phoretic particles) in a flow is key to a range of fields encompassing aquatic ecology (Kiørboe Reference Kiørboe2008; Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012), biomedicine (Nelson, Kaliakatsos & Abbott Reference Nelson, Kaliakatsos and Abbott2010) and active matter (Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013). Recently, microswimmers have gathered much attention in virtue of their collective behaviour (Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013; Elgeti, Winkler & Gompper Reference Elgeti, Winkler and Gompper2015). However, in the presence of a background flow, as typical in natural and artificial aquatic environments, the dynamics of microswimmers can be highly non-trivial even at the level of the single swimmer or for dilute (non-interacting) suspensions. Indeed, complex behaviour can originate from the interplay between swimming, advection, the fluid velocity gradients that act on the particles by changing their direction of propulsion, and the possible presence of biases (such as chemotaxis, gravi/gyrotaxis, phototaxis, etc.) affecting the swimming direction (for a review see, for example, Guasto et al. Reference Guasto, Rusconi and Stocker2012). For instance, Torney & Neufeld (Reference Torney and Neufeld2007, Reference Torney and Neufeld2008) were among the first to put forward the possibility of fractal clustering and non-trivial mixing properties for elongated microswimmers, with or without phototaxis, in Taylor–Green vortical flows. Rusconi, Guasto & Stocker (Reference Rusconi, Guasto and Stocker2014) found that cell elongation is responsible for aggregation and transport suppression of bacteria in shear flows due to Jeffery orbits (Jeffery Reference Jeffery1922). Reduced transport and vortical trapping were observed for spherical microswimmers in a chaotic flow (Khurana, Blawzdziewicz & Ouellette Reference Khurana, Blawzdziewicz and Ouellette2011).
In this work we focus on gyrotactic swimmers in turbulent flows. Many motile phytoplankton species are gyrotactic, i.e. they swim in the direction resulting from the competition between stabilizing torque due to buoyancy/gravity and the shear-induced viscous torque (Kessler Reference Kessler1985; Pedley & Kessler Reference Pedley and Kessler1992). The stabilizing torque can originate either from bottom heaviness or from a fore–rear asymmetry (Ten Hagen et al. Reference Ten Hagen, Kümmel, Wittkowski, Takagi, Löwen and Bechinger2014; Sengupta, Carrara & Stocker Reference Sengupta, Carrara and Stocker2017), which tend to keep the swimming direction oriented upwards, favouring vertical migration towards well-lit waters near the surface. Like other forms of biased motility, gyrotaxis can impinge the transport properties and spatial distribution of swimming plankton. In both laminar and turbulent pipe flows, gyrotaxis leads to peculiar transport properties with respect to passive particles (Bearon, Bees & Croze Reference Bearon, Bees and Croze2012; Croze et al. Reference Croze, Sardina, Ahmed, Bees and Brandt2013). In laminar flows, it induces remarkable beam-like accumulations in downwelling pipe flows (Kessler Reference Kessler1985) and in vortical flows (Cencini et al. Reference Cencini, Franchino, Santamaria and Boffetta2016), and high-concentration layers in horizontal shear flows (Durham, Kessler & Stocker Reference Durham, Kessler and Stocker2009; Santamaria et al. Reference Santamaria, De Lillo, Cencini and Boffetta2014). In moderately turbulent flows, as in the oceans, gyrotaxis can generate intense microscale fractal patchiness (Durham et al. Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013; Zhan et al. Reference Zhan, Sardina, Lushi and Brandt2014; Fouxon & Leshansky Reference Fouxon and Leshansky2015; Gustavsson et al. Reference Gustavsson, Berglund, Jonsson and Mehlig2016) and, in strong turbulence, accumulation in vortical regions (De Lillo et al. Reference De Lillo, Cencini, Durham, Barry, Stocker, Climent and Boffetta2014). Moreover, gyrotactic cells in turbulence have been found to preferentially sample specific flow regions (Durham et al. Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013; Fouxon & Leshansky Reference Fouxon and Leshansky2015; Gustavsson et al. Reference Gustavsson, Berglund, Jonsson and Mehlig2016). Indeed, direct numerical simulations (DNS) in both turbulent (multiscale) and stochastic (single-scale) flows, corroborated by theoretical considerations, have shown that spherical gyrotactic cells preferentially visit downwelling regions of the flow (Durham et al. Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013; Fouxon & Leshansky Reference Fouxon and Leshansky2015; Gustavsson et al. Reference Gustavsson, Berglund, Jonsson and Mehlig2016), which can hinder the upward vertical cell migration. Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016) also found a dependence on cell morphology in stochastic flows: elongated swimming particles can preferentially sample either downwelling or upwelling flow regions, depending on their swimming speed and stability. This is rather intriguing since some gyrotactic cells can change their morphology and biasing direction when exposed to turbulence (Sengupta et al. Reference Sengupta, Carrara and Stocker2017) and in this way control their motility. Moreover, Zhan et al. (Reference Zhan, Sardina, Lushi and Brandt2014) and Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016) found that turbulent suspensions of elongated gyrotactic cells are generally less clustered than spherical ones, but for the regime in which gyrotaxis is very weak.
Here we aim at clarifying the differences in the dynamics of turbulent suspensions of spherical and elongated gyrotactic cells by means of DNS. The rest of the paper is organized as follows. In § 2, we present the model and the numerical procedure. In §§ 3.1 and 3.2 we present our results on preferential sampling of the vertical component of fluid velocity and on fractal clustering, respectively. We end in § 4 by summarizing and discussing our findings.
2 Model equations and previous results
Following Kessler (Reference Kessler1985) and Pedley & Kessler (Reference Pedley and Kessler1992), we model gyrotactic swimmers as axisymmetric ellipsoids swimming with constant speed
$v_{s}$
in the direction
$\boldsymbol{p}$
(
$|\boldsymbol{p}|=1$
), along the axis of symmetry of the ellipsoid. Owing to their size being typically smaller than the Kolmogorov scale and their density being very close to that of the carrier fluid, the microswimmers behave as tracers if not for the motility. Neglecting interactions between swimmers (as appropriate for dilute suspensions) and stochastic reorientations, the evolution of the position,
$\boldsymbol{x}$
, and orientation,
$\boldsymbol{p}$
, of each swimmers is ruled by the following equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{x}}=\boldsymbol{u}(\boldsymbol{x},t)+\unicode[STIX]{x1D6F7}\boldsymbol{p}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{p}}=\frac{1}{2\unicode[STIX]{x1D6F9}}[\hat{\boldsymbol{z}}-(\hat{\boldsymbol{z}}\boldsymbol{\cdot }\boldsymbol{p})\boldsymbol{p}]+\frac{1}{2}\unicode[STIX]{x1D74E}\times \boldsymbol{p}+\unicode[STIX]{x1D6FC}[\unicode[STIX]{x1D64E}\boldsymbol{p}-(\boldsymbol{p}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{p})\boldsymbol{p}], & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}$
,
$\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$
and
$\unicode[STIX]{x1D61A}_{ij}=(\unicode[STIX]{x2202}_{j}u_{i}+\unicode[STIX]{x2202}_{i}u_{j})/2$
are the fluid velocity, vorticity and the rate of strain tensor at the swimmer’s position. The three terms on the right-hand side of (2.2) are: the buoyancy/gravity torque tending to rotate the cell upwards (along the vertical unit vector
$\hat{\boldsymbol{z}}$
) with a characteristic time
$B$
; rotation due to vorticity and strain rate. The latter, so-called Jeffery term (Jeffery Reference Jeffery1922), is proportional to
$\unicode[STIX]{x1D6FC}=(l^{2}-d^{2})/(l^{2}+d^{2})$
(with
$l$
the length of the cell, along
$\boldsymbol{p}$
, and
$d$
its width) and is absent for spherical cells (
$\unicode[STIX]{x1D6FC}=0$
). In general
$-1<\unicode[STIX]{x1D6FC}<1$
, with the extremes corresponding to flat disks and rods, respectively. In what follows, we will consider only
$\unicode[STIX]{x1D6FC}\geqslant 0$
, thus restricting ourselves to the case of spherical or prolate cells, being the typical shapes of microorganisms. We remark that in (2.2) we have not included a rotational stochastic term (Hill & Häder Reference Hill and Häder1997) due to the biological reorientation of the swimming direction (Polin et al.
Reference Polin, Tuval, Drescher, Gollub and Goldstein2009) since this effect is expected to be negligible with respect the reorientation induced by the small-scale turbulent flow.
As for the fluid velocity field,
$\boldsymbol{u}$
, we will consider moderately turbulent flows described by the Navier–Stokes equation for incompressible fluids (
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$
)
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f},\end{eqnarray}$$
where
$\unicode[STIX]{x1D708}$
is the kinematic viscosity and
$p$
the pressure rescaled by the fluid density. The turbulent flow is maintained in a statistically stationary state by the volume force
$\boldsymbol{f}$
, which injects energy at rate
$\unicode[STIX]{x1D716}$
, statistically equal to the rate of dissipation. Based on the flow properties one can define the usual small-scale parameters, namely the Kolmogorov scale
$\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D716})^{1/4}$
, time
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})^{1/2}$
and velocity
$u_{\unicode[STIX]{x1D702}}=(\unicode[STIX]{x1D708}\unicode[STIX]{x1D716})^{1/4}$
. The large dynamics can be parameterized in terms of the r.m.s. value of the velocity
$u_{rms}=\sqrt{\langle |\boldsymbol{u}|^{2}\rangle /3}$
(with
$\langle \cdots \rangle$
denoting ensemble average). In the following, physical variables will be made dimensionless by rescaling spatial and temporal coordinates with the Kolmogorov scale
$\unicode[STIX]{x1D702}$
and time
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
, respectively, and velocities with the Kolmogorov velocity
$u_{\unicode[STIX]{x1D702}}$
. Correspondingly, swimmer motion will be governed by the dimensionless swimming number
$\unicode[STIX]{x1D6F7}=v_{s}/u_{\unicode[STIX]{x1D702}}$
and stability number
$\unicode[STIX]{x1D6F9}=B/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
, while the flow is parameterized in terms of the Taylor-scale Reynolds number
$Re_{\unicode[STIX]{x1D706}}=u_{rms}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$
, where
$\unicode[STIX]{x1D706}=\sqrt{15\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716}}u_{rms}$
.
Based on the model equations (2.1), (2.2), recent numerical and theoretical works have provided insights into the properties of clustering of gyrotactic swimmers in both turbulent and chaotic flows (Durham et al.
Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013; Zhan et al.
Reference Zhan, Sardina, Lushi and Brandt2014; Fouxon & Leshansky Reference Fouxon and Leshansky2015; Gustavsson et al.
Reference Gustavsson, Berglund, Jonsson and Mehlig2016). In particular, it was first shown in Durham et al. (Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013) that spherical gyrotactic swimmers in turbulence can form fractal clusters. A phenomenological argument for this goes as follows. In the limit of fast orientation
$\unicode[STIX]{x1D6F9}\rightarrow 0$
, cells tend to swim upwards (
$\boldsymbol{p}\rightarrow \hat{\boldsymbol{z}}$
) so that their trajectories would approximately follow the effective velocity field
$\boldsymbol{v}=\boldsymbol{u}+\unicode[STIX]{x1D6F7}\hat{\boldsymbol{z}}$
, which is incompressible,
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0$
, and thus they do not cluster. Also in the limit of strongly unstable cells,
$\unicode[STIX]{x1D6F9}\gg 1$
, there would be no accumulation, since the swimming direction would be essentially random. However, when
$\unicode[STIX]{x1D6F9}\ll 1$
(and
$\unicode[STIX]{x1D6F7}\ll 1$
) a perturbative solution of (2.2) suggests (in analogy to inertial particles in the limit of small Stokes number (Balkovsky, Falkovich & Fouxon Reference Balkovsky, Falkovich and Fouxon2001; Bec Reference Bec2003)) that the effective velocity field,
$\boldsymbol{v}=\boldsymbol{u}+\unicode[STIX]{x1D6F7}\boldsymbol{p}^{\ast }$
with
$\boldsymbol{p}^{\ast }=(\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D714}_{y},-\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D714}_{x},1)$
, is compressible with
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=-\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6FB}^{2}u_{z}$
. Hence, for
$\unicode[STIX]{x1D6F9}\ll 1$
, cells behave as tracers in a compressible velocity field, and evolve onto a fractal attractor. On the basis of this argument, one concludes that cells will preferentially sample downwelling portions of the flow, where
$u_{z}<0$
(Durham et al.
Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013) (see Fouxon & Leshansky Reference Fouxon and Leshansky2015, for a refined derivation). The same conclusion was obtained with a perturbative approach in a stochastic single-scale velocity field in Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016). In the latter paper, however, the authors consider the general case of elongated swimmers,
$\unicode[STIX]{x1D6FC}\geqslant 0$
. Analytical results and simulations of stochastic model indicate that while preferential sampling of downwards velocities is typical of the limit
$\max (\unicode[STIX]{x1D6F7},\unicode[STIX]{x1D6F9})\ll 1$
for small-aspect-ratio swimmers, this is reversed for
$\unicode[STIX]{x1D6FC}>\unicode[STIX]{x1D6FC}_{c}=3/5$
, i.e. sufficiently elongated swimmers concentrate in upwelling regions. Furthermore, the inversion is predicted to be present only for
$\unicode[STIX]{x1D6F7}>\unicode[STIX]{x1D6F7}_{c}(\unicode[STIX]{x1D6F9})$
, with a critical swimming parameter decreasing with
$\unicode[STIX]{x1D6F9}$
.
2.1 Numerical simulations
The Navier–Stokes equations (2.3) are solved in a three-dimensional, periodic cubic domain of size
$L=2\unicode[STIX]{x03C0}$
by means of a parallel, fully dealiased pseudospectral code at resolutions of
$64^{3}$
,
$256^{3}$
and
$1024^{3}$
grid points, corresponding to
$Re_{\unicode[STIX]{x1D706}}=21,68$
and 173. A constant energy input is provided by a deterministic, large-scale forcing
$\boldsymbol{f}$
acting on modes with wavenumbers in the spherical shell
$1<|\boldsymbol{k}|<3$
. Time integration is implemented by a second-order Runge–Kutta scheme with exact integration of the dissipative term. Accuracy at small scales is guaranteed by ensuring that
$k_{max}\unicode[STIX]{x1D702}\gtrsim 1.9$
for all simulations.
In each run, the trajectories of up to
$128\,000$
swimmers per parameter set are integrated according to (2.1), (2.2). The Eulerian velocity is interpolated at the particle positions by third-order polynomials. The derivatives needed for the calculation of
$\unicode[STIX]{x1D74E}$
and
$\unicode[STIX]{x1D61A}_{ij}$
are computed by using the derivatives of the interpolating polynomials, thus ensuring that they are second-order accurate. For each
$Re_{\unicode[STIX]{x1D706}}$
we have explored a wide range of cells’ parameters, with
$\unicode[STIX]{x1D6F9}\in [0.1:50]$
,
$\unicode[STIX]{x1D6F7}\in [0.1:30]$
and shape parameter
$\unicode[STIX]{x1D6FC}=0,0.5,1$
, thus considering up to 90 different types of swimmer per simulation. Swimmers are initialized with uniform random positions
$\boldsymbol{x}$
in the domain and orientations
$\boldsymbol{p}$
on the unit sphere. At stationarity, data are collected for several configurations (120–200) separated by about half a large-scale eddy turnover time to ensure statistical convergence.
3 Results
3.1 Preferential sampling of fluid velocities
Figure 1. Preferential sampling of the vertical component of fluid velocity by gyrotactic, rod-like swimmers (
$\unicode[STIX]{x1D6FC}=1$
). (a) Average vertical component of the fluid velocity measured along cell trajectories. Curves are computed at
$Re_{\unicode[STIX]{x1D706}}=173$
. (b,c) Distributions of swimmers for populations marked in (a) with
$\unicode[STIX]{x1D6F9}=1$
,
$\unicode[STIX]{x1D6F7}=5$
(b) and
$\unicode[STIX]{x1D6F9}=10$
,
$\unicode[STIX]{x1D6F7}=20$
(c) plotted on a colour map of the vertical velocity
$u_{z}$
. Slower (b) or faster (c) swimmers sample preferentially regions with negative (blue) or positive (red) vertical velocity.
The vertical displacement in the water column of a gyrotactic swimmer following (2.1) can be affected by fluid transport even in a flow with vanishing average velocity. Indeed, although the Eulerian average of the fluid velocity vanishes, the average swimmer velocity
$\langle \boldsymbol{v}\rangle =\langle \boldsymbol{u}\rangle _{S}+\unicode[STIX]{x1D6F7}\langle p_{z}\rangle \hat{\boldsymbol{z}}$
depends on the mean value of
$\boldsymbol{u}$
computed along particle trajectories (notice that in the above expression symmetry implies
$\langle p_{x,y}\rangle =0$
). While for the horizontal components one has
$\langle u_{x}\rangle _{S}=\langle u_{y}\rangle _{S}=0$
for symmetry reasons, gravitaxis breaks the up–down symmetry so that one can have
$\langle u_{z}\rangle _{S}\neq 0$
(Fouxon & Leshansky Reference Fouxon and Leshansky2015). Since the largest swimming speeds of motile algae are typically of the order of
$u_{\unicode[STIX]{x1D702}}$
or smaller and
$u_{rms}\gg u_{\unicode[STIX]{x1D702}}$
, preferential sampling of regions of upwards or downwards flow can in principle strongly enhance or hinder the swimmer’s drift towards the surface. Figure 1(a) shows the average value of the vertical component of the fluid velocity along swimmer trajectories
$\langle u_{z}\rangle _{S}$
. The different curves are computed for rod-like cells (
$\unicode[STIX]{x1D6FC}=1$
) and
$\langle u_{z}\rangle _{S}$
is plotted for several values of
$\unicode[STIX]{x1D6F9}$
as a function of the swimming parameter,
$\unicode[STIX]{x1D6F7}$
. It is clear that slow swimmers tend to preferentially sample downwards fluid velocities (figure 1
b), while faster ones (figure 1
c) spend more time where
$u_{z}>0$
.
Figure 2. Small-
$\unicode[STIX]{x1D6F7}$
behaviour of preferential sampling of fluid velocities and small-scale clustering, compared with the theoretical predictions in Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016) (dashed lines). (a) Vertical component of fluid velocity conditioned on swimmer trajectories for different elongations as labelled at
$\unicode[STIX]{x1D6F9}=0.1$
. Inset: same as main panel but for rods (
$\unicode[STIX]{x1D6FC}=1$
) at different
$\unicode[STIX]{x1D6F9}$
. (b) Accumulation index
$N$
rescaled with the stochastic prediction. Inset: the same curves without rescaling.
Figure 3. Average of the vertical fluid velocity,
$\langle u_{z}\rangle _{S}$
, along swimmer trajectories for rod-like (
$\unicode[STIX]{x1D6FC}=1$
) cells as a function of
$\unicode[STIX]{x1D6F7}=v_{s}/u_{\unicode[STIX]{x1D702}}$
(a–c), and of
$\unicode[STIX]{x1D6F7}_{L}=v_{s}/u_{rms}$
(d–f) for different
$Re_{\unicode[STIX]{x1D706}}$
, as labelled in (a). Whereas the small-
$\unicode[STIX]{x1D6F7}$
behaviour in (a–c) rescales consistently, the transition to preferential sampling of positive values of
$u_{z}$
is better rescaled by
$\unicode[STIX]{x1D6F7}_{L}$
.
For a more quantitative analysis, we focus at first on small values of
$\unicode[STIX]{x1D6F7}$
. In this limit the behaviour of the curves is consistent with what observed is in (Durham et al.
Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013), in that the average fluid velocity seen by the cells is increasingly negative as cell velocity grows. This behaviour has been discussed for a
$d$
-dimensional, random, single-scale flow in Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016), where it was predicted that in the small-
$\unicode[STIX]{x1D6F7}$
limit,
$\langle u_{z}\rangle _{S}\propto -\unicode[STIX]{x1D6F7}G_{d}(\unicode[STIX]{x1D6FC})K(\unicode[STIX]{x1D6F9})$
, where
$G_{d}(\unicode[STIX]{x1D6FC})=[d(1-\unicode[STIX]{x1D6FC})+2]/d$
and
$K(\unicode[STIX]{x1D6F9})=\unicode[STIX]{x1D6F9}/(2\unicode[STIX]{x1D6F9}+1)$
. That prediction is remarkably verified in our simulations, as shown in figure 2, where the average fluid velocity along the trajectories is plotted for different values of the shape parameter
$\unicode[STIX]{x1D6FC}$
, showing that the factor due to particle geometry
$G_{d}(\unicode[STIX]{x1D6FC})$
is indeed correct. When the shape parameter
$\unicode[STIX]{x1D6FC}$
is fixed, the dependence on
$\unicode[STIX]{x1D6F9}$
is also fairly well satisfied (see inset). A similar prediction is confirmed also on the small-
$\unicode[STIX]{x1D6F7}$
properties of clustering as shown in figure 2(b), which will be discussed in the next section. In general, we can also conclude that the small-
$\unicode[STIX]{x1D6F7}$
behaviour of this observable is independent of
$Re_{\unicode[STIX]{x1D706}}$
. Figure 3(a–c) show the vertical fluid speed curves for fixed
$\unicode[STIX]{x1D6F9}$
at various
$Re_{\unicode[STIX]{x1D706}}$
. Notice that in our non-dimensional units, based on the Kolmogorov velocity, the curves collapse in the small-
$\unicode[STIX]{x1D6F7}$
region for each value of
$\unicode[STIX]{x1D6F9}$
(see also Durham et al.
Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013, for the case of spherical cells).
We now consider larger values of the swimming velocity. As noted above, fast swimming cells tend to preferentially sample regions of upwards flow. As shown by figure 3, the critical speed
$\unicode[STIX]{x1D6F7}_{c}$
of transition from downwelling to upwelling regions (defined as the point at which
$\langle u_{z}\rangle _{S}=0$
) is smaller for larger values of
$\unicode[STIX]{x1D6F9}$
. However, numerical results show that
$\unicode[STIX]{x1D6F7}_{c}$
saturates to a finite value for large
$\unicode[STIX]{x1D6F9}$
. This has profound implications for the statistics of velocity and velocity gradients along particle trajectories. While the correlation time of gradients is of the order of
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
, the time it takes for a swimmer to cross the correlation length of gradients,
$\unicode[STIX]{x1D702}$
, is of the order of
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6F7}}=\unicode[STIX]{x1D702}/v_{s}=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6F7}$
, so that
$\unicode[STIX]{x1D6F7}>1$
implies that the gradients correlation time as seen by a swimmer is
$\unicode[STIX]{x1D70F}_{corr}=\min (\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6F7}},\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}})=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6F7}}$
(Fouxon & Leshansky Reference Fouxon and Leshansky2015). As a consequence, the preferential sampling of upwards velocities cannot be understood in terms of a quasiequilibrium solution, at variance with the
$\langle u_{z}\rangle _{S}<0$
case for spherical particles. However, an even more important difference is revealed if the actual observed values of
$\unicode[STIX]{x1D6F7}_{c}$
are considered. Indeed, when
$\langle u_{z}\rangle _{S}$
is plotted at fixed
$\unicode[STIX]{x1D6F9}$
upon changing the extension of the inertial range with simulations at different resolutions, one can see that
$\unicode[STIX]{x1D6F7}_{c}$
increases with
$Re_{\unicode[STIX]{x1D706}}$
, and the corresponding values of the swimming velocity are of the order of (or larger than) the large-scale velocity of the flow
$u_{rms}$
. Figure 3(d–f) shows the
$\langle u_{z}\rangle _{S}$
as a function of a large-scale swimming parameter
$\unicode[STIX]{x1D6F7}_{L}=v_{s}/u_{rms}$
. In this parameterization, the critical values
$\unicode[STIX]{x1D6F7}_{L,c}(\unicode[STIX]{x1D6F9})$
are independent of
$Re_{\unicode[STIX]{x1D706}}$
. For large
$\unicode[STIX]{x1D6F9}$
, the critical
$\unicode[STIX]{x1D6F7}_{L}$
appears to asymptotically converge from above to
$\unicode[STIX]{x1D6F7}_{L,\infty }\approx 0.86$
, which can be obtained as the
$\unicode[STIX]{x1D6F9}\gg 1$
limit of the analytical expression in Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016). This numerical result is summarized in figure 4 and is consistent with a picture in which the transition to upwards velocity sampling is essentially controlled by the integral scale of the flow.
Figure 4. Critical value of
$\unicode[STIX]{x1D6F7}_{L}$
(main panel) and
$\unicode[STIX]{x1D6F7}$
(inset) for various
$Re_{\unicode[STIX]{x1D706}}$
as a function of
$\unicode[STIX]{x1D6F9}$
. Dashed line: theoretical asymptotic value of
$\unicode[STIX]{x1D6F7}_{L,c}$
(see text).
The fact that the inversion in preferential sampling is controlled by the large-scale velocity has not been noted before, to our knowledge. Indeed it can be observed only if a sufficient separation between large and small scales is present in the velocity field. It is therefore important to note that this scale separation effect is not in contradiction to the stochastic results in Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016), and indeed the asymptotic value of the critical swimming parameter defined with respect to
$u_{rms}$
is in quantitative agreement with the prediction found in the same paper (dashed line in figure 4). The techniques adopted in that paper relied on a flow with a single scale, so that both gradients and velocities have the same correlation times. It is remarkable that those techniques are able to capture a qualitative behaviour that is observed in a multiscale flow. The observation that those features, in high-
$Re_{\unicode[STIX]{x1D706}}$
flows, are found for very large swimming velocities (of the order of
$u_{rms}$
instead of
$u_{\unicode[STIX]{x1D702}}$
) can perhaps suggest why the qualitative agreement between stochastic theory and simulations is so good. Indeed, for very fast swimmers, as observed above, fluid fields can be described as noise with correlation time inversely proportional to the swimming speed, losing the details of the multiscale structure of the flow.
3.2 Fractal clustering
We now focus on the small-scale, fractal clustering properties of gyrotactic suspensions by varying the shape parameter,
$\unicode[STIX]{x1D6FC}$
, with emphasis on the spherical (
$\unicode[STIX]{x1D6FC}=0$
) and rod-like (
$\unicode[STIX]{x1D6FC}=1$
) limits. We start by considering the limit of stable (
$\unicode[STIX]{x1D6F9}\ll 1$
) and slowly swimming (
$\unicode[STIX]{x1D6F7}\ll 1$
) cells.
In such limit, the perturbative approach originally developed in Durham et al. (Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013) for spherical cells (see § 2) can be extended to generic shapes by imposing stationarity in (2.2) (Gustavsson et al.
Reference Gustavsson, Berglund, Jonsson and Mehlig2016). At the first order in
$\unicode[STIX]{x1D6F9}$
one obtains that gyrotactic cells behave as tracers in an effective velocity field
$\boldsymbol{v}$
$$\begin{eqnarray}\dot{\boldsymbol{x}}=\boldsymbol{v}=\boldsymbol{u}+\unicode[STIX]{x1D719}\boldsymbol{p}^{\ast },\quad \text{with }\boldsymbol{p}^{\ast }=[\unicode[STIX]{x1D6F9}(\unicode[STIX]{x1D714}_{y}+2\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D61A}_{13}),\unicode[STIX]{x1D6F9}(-\unicode[STIX]{x1D714}_{x}+2\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D61A}_{23}),1].\end{eqnarray}$$
The resulting velocity field,
$\boldsymbol{v}$
, is compressible with divergence
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=-\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9}[(1+\unicode[STIX]{x1D6FC})\unicode[STIX]{x2202}_{z}^{2}u_{z}+(1-\unicode[STIX]{x1D6FC})(\unicode[STIX]{x2202}_{x}^{2}+\unicode[STIX]{x2202}_{y}^{2})u_{z}].\end{eqnarray}$$
For rods, this expression reduces to
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=-2\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9}\unicode[STIX]{x2202}_{z}^{2}u_{z}$
, while for spheres it gives
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=-\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6FB}^{2}u_{z}$
. Tracers advected by a weakly compressible flow field
$\boldsymbol{v}(\boldsymbol{x},t)$
form clusters of fractal codimension
$3-D\sim \langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}\rangle _{S}$
(Falkovich, Fouxon & Stepanov Reference Falkovich, Fouxon and Stepanov2002). Thus one can infer that swimmers should form clusters with codimension
$3-D\sim (\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9})^{2}$
. For a generic
$\unicode[STIX]{x1D6FC}$
, in a
$d$
-dimensional stochastic Gaussian flow
$\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}\rangle _{S}$
can be computed exactly at first order of a perturbative approach based on a Kubo expansion (Gustavsson et al.
Reference Gustavsson, Berglund, Jonsson and Mehlig2016), obtaining the prediction
$$\begin{eqnarray}3-D\sim (\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9})^{2}C_{d}(\unicode[STIX]{x1D6FC})\quad [C_{d}(\unicode[STIX]{x1D6FC})=[(d+2)(d+4)-2d(d+4)\unicode[STIX]{x1D6FC}+(4+2d+d^{2})\unicode[STIX]{x1D6FC}^{2}]/d].\end{eqnarray}$$
We measured the fractal dimension of the clusters in terms of the correlation dimension,
$D_{2}$
, ruling the small-scale scaling behaviour of the probability to find a pair of cells with distance less than
$r$
, i.e.
$p_{2}(r)\sim r^{D_{2}}$
when
$r\rightarrow 0$
. We estimated
$D_{2}$
from the local slopes of such probability. When
$D_{2}\approx 3$
, as expected in the limit
$\unicode[STIX]{x1D6F9},\unicode[STIX]{x1D6F7}\ll 1$
here considered, this kind of measure is affected by large errors, so we also measured the accumulation index
$N$
which is obtained as follows. The volume is divided in cubes of side
$\ell$
, and we count the average number of particles in such cubes,
$\langle n\rangle _{\ell }$
, and the standard deviation,
$\unicode[STIX]{x1D70E}^{2}=\langle n^{2}\rangle _{\ell }-\langle n\rangle _{\ell }^{2}$
. Then we define
$N=(\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D70E}_{P})/\langle n\rangle _{\ell }$
, where
$\unicode[STIX]{x1D70E}_{P}=\sqrt{\langle n\rangle _{\ell }}$
is the standard deviation of a Poissonian process having the same mean number of particles. As discussed in Durham et al. (Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013), one can show that if
$\ell$
is small enough (here we have chosen
$\ell \approx 2\unicode[STIX]{x1D702}$
) then
$N\propto 3-D_{2}$
.
In figure 2(b) we show
$N$
as a function of
$(\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6F7})^{2}C_{d}(\unicode[STIX]{x1D6FC})$
for different values of
$\unicode[STIX]{x1D6F9},\unicode[STIX]{x1D6F7}$
and
$\unicode[STIX]{x1D6FC}=0,1/2,1$
. We observe a fair collapse of the data (compare with the inset where the geometric constant
$C_{d}(\unicode[STIX]{x1D6FC})$
is removed) onto a curve that displays a linear behaviour for small values of the abscissa, as predicted by (3.3), providing again (see also the discussion of figure 2
a) strong evidence that the first-order Kubo expansion of Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016) reproduces the small-
$\unicode[STIX]{x1D6F9}$
and small-
$\unicode[STIX]{x1D6F7}$
limits of gyrotactic cells also in turbulent multiscale flows. We observe that the linear limiting behaviour seems to be more persistent for the spherical shapes. Similarly to inertial particles, which display the maximal clustering for Stokes number of order
$1$
(Bec Reference Bec2003), gyrotactic cells also have a more pronounced clustering for
$\unicode[STIX]{x1D6F9}\sim 1$
(Durham et al.
Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013). Zhan et al. (Reference Zhan, Sardina, Lushi and Brandt2014) showed that spherical cells are generally more clustered than elongated swimmers, except for very large
$\unicode[STIX]{x1D6F9}$
. In figure 5 we compare vis a vis the behaviour of
$D_{2}$
(estimated from the local slopes of
$p_{2}(r)$
) for spherical (
$\unicode[STIX]{x1D6FC}=0$
) and rod-like (
$\unicode[STIX]{x1D6FC}=1$
) as a function of the swimming number
$\unicode[STIX]{x1D6F7}$
for different values of
$\unicode[STIX]{x1D6F9}\geqslant 1$
. By comparing figure 5(a,c), it is clear that spherical cells cluster more than elongated ones when
$\unicode[STIX]{x1D6F9}$
is small enough, while the opposite is true for unstable cells
$\unicode[STIX]{x1D6F9}\gg 1$
. The latter limiting behaviour can be explained by noticing that at
$\unicode[STIX]{x1D6F9}\rightarrow \infty$
the dynamics of (2.1)–(2.2) in the (
$\boldsymbol{x},\boldsymbol{p}$
)-phase space has an average volume contraction rate equal to
$-\unicode[STIX]{x1D6FC}(d+2)\langle \,\boldsymbol{p}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\,\boldsymbol{p}\rangle _{S}$
, which is identically zero for spherical cells. In other terms, the dynamics of spherical cells with large
$\unicode[STIX]{x1D6F9}$
preserves volumes in phase space, and so does its projection in position space (i.e.
$D_{2}=3$
), while elongated swimmers obey a dissipative dynamics which leads to fractal clustering. For intermediate values of
$\unicode[STIX]{x1D6F9}$
(figure 5
b) a non-trivial role is played by the swimming speed: slow spherical swimmers cluster more than elongated ones, but the opposite is true for fast swimmers.
Figure 5. Correlation dimension of clusters for
$\unicode[STIX]{x1D6F9}=0.5$
(a), 1.0 (b) and 10 (c) as a function of the swimming parameter
$\unicode[STIX]{x1D6F7}$
. Spherical swimmers (
$\unicode[STIX]{x1D6FC}=0$
) concentrate more than elongated cells (in this case, ideal rods) for small values of
$\unicode[STIX]{x1D6F9}$
(a), whereas rod-like swimmers maintain observable clustering even for large
$\unicode[STIX]{x1D6F9}$
(c). For intermediate values of
$\unicode[STIX]{x1D6F9}$
, slow (fast) elongated swimmers cluster less (more) than spherical ones (b).
4 Conclusions and discussion
A systematic study of the statistical properties of gyrotactic swimmers, on varying their shape (from spherical to rod-like) and other cell parameters, in realistic turbulent flows has been presented. In particular, we focused on the characterization of their aggregation in clusters and preferential sampling of the vertical component of the fluid velocity field. The first property is relevant to many biological processes, such as reproduction and predation (Kiørboe Reference Kiørboe2008), while the second has an impact on the vertical migration of phytoplankton cells. The main conclusions are as follows: (i) Sufficiently slow swimmers (with swimming speeds
$v_{s}\lesssim u_{\unicode[STIX]{x1D702}}$
), independently of their shape, always sample downwelling regions, the effect being maximal for non-dimensional reorientation times,
$\unicode[STIX]{x1D6F9}$
, of order
$1$
and increasing with the swimming speed. (ii) Fast elongated swimmers can sample upwelling regions (thus acquiring an extra vertical migration speed), whereas this behaviour is absent for spherical cells (as previously shown in Durham et al. (Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013), Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016) and confirmed here by simulations at large
$\unicode[STIX]{x1D6F7}$
, not shown). Simulations at different Reynolds numbers show that the transition from downwelling to upwelling regions is controlled by the large-scale velocity, and thus requires
$v_{s}\gtrsim u_{rms}$
. (iii) Both spherical and elongated cells form fractal aggregates, which are more clustered for spherical (elongated) cells in the regime of fast (slow) reorientation. For intermediate reorientation, the relative strength of clustering depends on the swimming speed. (iv) We provided evidence that some analytical results derived by Gustavsson et al. (Reference Gustavsson, Berglund, Jonsson and Mehlig2016) for a stochastic flow are qualitatively in agreement with numerical data obtained from fully resolved turbulent flows. The agreement becomes quantitative in the regime of slow swimming (figure 2). This is remarkable since the analytical results were obtained by a first-order Kubo expansion for a single-scale flow.
Our results clarify several aspects of the dynamics and the statistics of gyrotactic suspensions and are of particular importance for potential application to swimming phytoplankton, also in virtue of recent evidence of possible morphological adaptations (Sengupta et al.
Reference Sengupta, Carrara and Stocker2017). Most gyrotactic microorganisms have reorientation times in the range 2–10 s, with swimming speeds around
$100{-}300~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$
and cell sizes around
$10~\unicode[STIX]{x03BC}\text{m}$
(Guasto, Johnson & Gollub Reference Guasto, Johnson and Gollub2010; Harvey, Menden-Deuer & Rynearson Reference Harvey, Menden-Deuer and Rynearson2015). Typical values for the energy dissipation rate in the ocean are in the range
$\unicode[STIX]{x1D716}\sim 10^{-4}{-}10^{-8}$
(Thorpe Reference Thorpe2007) and therefore phytoplankton are able to explore the parameter range
$0.1\lesssim \unicode[STIX]{x1D6F9}\lesssim 20{-}40$
and
$0.1\lesssim \unicode[STIX]{x1D6F7}\lesssim 3{-}5$
in which small-scale clustering is expected for both spherical and elongated cells. In contrast, given that the large-scale velocity is typically of order
$\text{cm}~\text{s}^{-1}$
or
$\text{m}~\text{s}^{-1}$
, it is unlikely that phytoplankton cells are able to reach sufficiently high
$\unicode[STIX]{x1D6F7}_{L}$
in order to exploit preferential concentration in upwelling regions. This regime could be in principle accessible to faster natural or artificial swimmers, as long as they are smaller than the Kolmogorov scale, as required for the validity of (2.1).
Acknowledgements
M.B., G.B. and F.D.L. acknowledge support by the Departments of Excellence grant (MIUR). HPC Center CINECA is gratefully acknowledged for computing resources, within the INFN-Cineca agreement INF18-fldturb.
