Table Of ContentA Lagrangian model of Copepod dynamics:
clustering by escape jumps in turbulence
H. Ardeshiri,1,2,∗ I. Benkeddad,2 F. G. Schmitt,2 S. Souissi,2 F. Toschi,3,4 and E. Calzavarini1
1Univ. Lille, CNRS, FRE 3723, LML, Laboratoire de M´ecanique de Lille, F 59000 Lille, France
2Univ. Lille, CNRS, Univ. Littoral Cote d’Opale, UMR 8187, LOG,
Laboratoire d’Oc´eanologie et de G´eoscience, F 62930 Wimereux, France
3Department of Applied Physics and Department of Mathematics and Computer Science,
Eindhoven University of Technology, 5600 MB, Eindhoven, The Netherlands
4Istituto per le Applicazioni del Calcolo CNR, Via dei Taurini 19, 00185 Rome, Italy
(Dated: April27, 2016)
Planktoniccopepodsaresmallcrustaceansthathavetheabilitytoswimbyquickpowerfuljumps.
Suchanaptnessisusedtoescapefrom high shearregions, which maybecaused eitherbyflowper-
6 turbations, produced by a large predator (i.e. fish larvae), or by the inherent highly turbulent
1 dynamicsof the ocean. Through a combined experimental and numerical study,we investigate the
0 impactofjumpingbehaviouronthesmall-scale patchinessofcopepodsinaturbulentenvironment.
2 Recorded velocity tracks of copepods displaying escape response jumps in still water are here used
to defineand tunea Lagrangian Copepod (LC) model. The model is further employed to simulate
r
p thebehaviourofthousandsofcopepodsinafullydevelopedhydrodynamicturbulentflowobtained
A bydirect numerical simulation of theNavier-Stokesequations. First, we show that theLCvelocity
statistics is in qualitative agreement with available experimental observations of copepods in tur-
6 bulence. Second,wequantify theclustering of LC,via thefractal dimension D2. Weshow that D2
2 can be as low as ∼ 2.3 and that it critically depends on the shear-rate sensitivity of the proposed
LC model, in particular it exhibits a minimum in a narrow range of shear-rate values. We further
] investigate the effect of jump intensity, jump orientation and geometrical aspect ratio of the cope-
n
podsonthesmall-scalespatialdistribution. Atlast,possibleecological implicationsoftheobserved
y
clustering on encounterrates and mating success are discussed.
d
-
u
l I. INTRODUCTION Copepods also consume the mosquito larvae, acting as
f
. control mechanism for malaria [26]. They are of great
s
c The study of swimming microorganisms and their inter- importance in fishery industry. A central issue in breed-
i action with fluid flows has attracted enormous attention ing fish species, is the external food supply. Most fishes
s
y in the last decade. A line of research has focused on prefer copepods to other zooplankton species (i.e. ro-
h characterizing individual swimming strategies by means tifers) and they grow bigger in shorter time when eating
p of experiments [1–3] as well as by theoretical and nu- copepods [27, 28].
[ merical modelling [4, 5]. A second direction of study Living in a fluid environment characterised by body-
2 devoted to the consequences of swimming on population scaleReynoldsnumberupto1000,they aresubjectedto
v dynamics, e.g., by focusing on encounter rates and other the physics of the flow field both in viscous and inertial
8 collectivebehaviours[6–10]. Athirddirectionfocusedon regime [29]. Copepods typically have a short, cylindri-
3 the mutualinteractionsofmicroorganismswiththe fluid cal body with antennas, few pairs of swimming legs and
4 flow environment, in particular bio-induced flow fluctu- tales. Using their antennas, copepods can sense the dis-
1
0 ations, sometimes dubbed as bacterial turbulence [11– turbance,whichiscausedeitherbythepresenceofpreda-
. 13], or, vice-versa, on active matter clustering induced tors or by high turbulent regions in the flow. Kiørboe
1 by non homogeneous flows or fluid turbulence [14–22]. et al. [30, 31] performed series of experiments, investi-
0
Thepresentstudywillfocusonthislatteraspect,inpar- gating the effect of non uniform flow motion on cope-
6
1 ticular on copepod’s dynamics in turbulent flow. pods. In order to find the component of the flow which
: Copepods are the most diversified crustaceans in the copepods react the most to, the copepods were put into
v aquatic environment whose length ranges from 0.1 mm a time dependent siphon flow (which ideally generates
i
X to few millimetres. They are important to global ecol- a pure longitudinal deformation rate), in an oscillating
r ogy and to the carbon cycle [23] (see also [24] & [25]). chamber where copepods experience only acceleration,
a Although copepods are not at the top of the food web, in a couette device producing shear deformation, and fi-
they have a major role in the marine ecosystem because nally in a rotating cylinder where acceleration and vor-
they are the secondary producers in the ecological food ticity areboth present. The conclusionofthis study was
weblinking phytoplanktoncells (the primaryproducers) that these small crustaceans react to the flow deforma-
tofishlarvaeandeventolargemammalssuchaswhales. tion rate. Kiørboe also reported [32], that there are two
threshold values of the deformation rate: the upper one,
around 10s−1, corresponds to either the presence of a
predators or to a region where turbulence intensity is
∗ [email protected]
2
high, and the lower one, 1s−1, corresponds to regions in works.
theflowwhereturbulenceintensityislowerorfoodabun-
danceisnotenoughforcopepods. Thesetinycrustaceans
find themselves at ease in regions in between these two II. METHODS
thresholds. To avoid uncomfortable regions, copepods
exhibit a rapid escape in the flow which is often dubbed
A. Experimental jump data analysis
a jump. Buskey et al. [33, 34] showed that copepod’s
velocity can reach the rate of 500 body length per sec-
We begin presenting an analysis of a new experimen-
ond (0.5 m/s) while jumping. The mechanical energy
tal trajectory data set of the estuarine copepod, Eu-
produced during their escape is reported to be very high
(8 10−5J/s) [35], which makes copepods, relative to rytemora affinis, recorded at LOG Laboratory between
× May and June 2015. Copepods originated from the
their size, among the fastest and the strongest animals
Seine river estuary (France) are maintained in the lab-
in the world.
oratoryunderoptimalconditionsforseveralgenerations.
Buskey [33, 34] also reported that males and females re-
The experimental set up is a shallow-depth aquarium,
spond differently to hydrodynamic stimulus in terms of
63 53 6mm3 inlength,heightanddepthrespectively,
response latency, jump speed, number of thrusts, dis- wit×h tw×o light sources on the lateral side (53 6mm2).
tance jumped andmany other parameters. According to × ◦
The water is kept still and at temperature of (18 1) C
their investigations, copepods jump in an unpredictable ±
and salinity of 15 psu. Copepods were introduced one
direction, but rarely in the backward direction of their
at a time in the aquarium and their dynamics filmed.
motion. Other studies have considered the mating be-
A total of 14 individuals were analysed (7 males and 7
haviourofcopepods[36]andtheeffectofsalinityoncope-
females). A copepod in the aquarium is lead to jump
pod’s dynamics and copepod’s encounter rate [37, 38].
preferentially along the horizontal direction by switch-
Copepods are also sensitive to light stimuli, being at-
ing on just one of the light sources. The copepod dy-
tracted by natural light sources [39].
namics in a vertical plane is recorded by a high speed
In the last two decades many studies have been con-
camera (1000 frames/second) and the single trajectory
ducted to quantify the dynamics of copepods. Most of
is extracted by means of a particle tracking velocimetry
them focused on their behaviour in still water [36–38],
software (TEMA Motion by Image Systems). In such a
whilelessstudieshavestudiedthedynamicsintheirnat-
way hundreds of trajectories are recorded, each with an
urallivingenvironmentbecauseofthedifficultiesofsuch
average time length of 19 s. A typical copepod velocity
experimental investigations. Few works have been de-
signal as a function of time is shown in figure 1(a). We
voted to the dynamics of copepods in turbulent flows
seeextremelyabruptspikes(jumps) alternatingtocalm,
[40–44]. However,thedensitiesofcopepodsusedinthese
nearly immobile, phases.
studies are often lowerthan the maximumdensities that
In order to see if the velocity signal of the jump events
can be encountered in the field.
share some common features, we zoom in on the signal
The numerical simulation can provide a tool that inte-
and superpose several jumps by a shift taking as refer-
grate our current knowledge on copepod dynamics and
ence their peak position. In figure 1(b) we can appre-
use high number of individuals. The objective of the
ciate that almost all of the jumps, after a steep rise,
presentstudyistosimulatecopepodsnumericallyintur-
display a similar decay. We associate such a decay to
bulence to characterise their dynamics induced by a be-
a purely hydrodynamical effect. It can be interpreted as
haviour model. To achieve this goal, our strategy is
a drag-induced decay of an instantaneous acceleration.
two-fold: on one hand, new experimental measurements
The inset of panel in figure 1, shows that the probabil-
andobservationsavailableintheliterature[45–51],along
ity density function (PDF) of the jump intensity, has a
with the aforementioned copepods properties, should be
maximum value around 0.07m/s. Note that spikes
consideredin detailsin orderto introduce arealistic and ∼
are identified based on a threshold on the time-averaged
physical model. On the other hand, fundamental knowl-
velocity of the copepods in each copepod’s trajectory.
edge on simulation of particles in turbulent flows, avail-
We then average the dataset of jumps in order to ob-
able in numerical and experimental studies on particles
tain an averaged shape of jump. This is shown in figure
inturbulence[52–56],isneededtocouplethephysicsand
2, from which we can deduce the average jump velocity
biology in the numerical model.
amplitude u =0.0939m/sand the mean decaying time
Thepaperisorganisedasfollows: sectionIIdescribesthe J
τ =8.87ms. We alsoseethatforlongtime the velocity
experimentalframeworkusedtostimulatecopepods. We J
reaches a very low plateau at 5 10−3m/s 1/20 u ,
thenanalysecopepod’s trajectoriestointroduceamodel × ∼ J
which we are tempted to associate to a weak random
equation describing copepods behaviour. Furthermore
wandering behaviour of the copepod.
similarityanalysisisperformedtotunetheLCmodeland
The distribution of jumps in time in the experimental
its numerical implementation is explained at the end of
dataset seems to deviate from an exponential distribu-
thissection. SectionIIIdetailsthesinglepointstatistics,
tion suggesting the existence of a memory effect. This
fractaldimension andorientationdynamics ofcopepods.
may however be dependent on the type of stimulus (the
The paper ends with conclusion and outlook on future
light source) which is continuous in time, very different
3
(a) 0.25 100
Standard Deviation
Average shape
m/s) 0.2 m/s) 10-1 uJ e(u-Jt/ τeJ()-,t /uτJJ) =+ 0n.o0i9s3e9, nmo/iss,e τ=J =0. 080.857 m m/ss
V| ( 0.15 V| (
y | y |
cit 0.1 cit
elo elo 10-2
V V
0.05
0 10-3
0 5000 10000 15000 20000 0 20 40 60 80 100 120
Time (ms) Time (ms)
(b) 102 FIG. 2. Average shape of copepod velocity (over about 730
Velocity |V| (m/s) 111000--021 PDF of Jump Intensity111000-011 0.04 Jum 0p.0 I8ntens 0it.y1 2(m/s) 0.16 0.2 jatbewauoxrneimeptdathpoph)τrn.steJee)hs:Nnueeistmnnoiatacttdeeheladbrefnttuiuthadntiviaoecnatnchtltauiayaooyesnnfiin(nbaerugeerJtnredtheonoi−lmeriisntoaee/iemdlτ)ovoJeinefantlwtgnttoeihdhcfidteiehcstratyehjeautehinomuroffodeJnp.rsaie.irzostdGoBf.tndhrljtueueeaeevmljniuladpimntliiirepnope:neecia:tSn(ikasotFshmenniamtedsdtieaeufitddyyet
10-3 the same way to external flow disturbances. Their jump
-20 0 20 40 60 80 100 120 reaction is embedded in their neural system. Further-
more,the stimulustriggeringthe jumpishighly stylised,
Time (ms)
we only take into account a mechanical signal with a
single-threshold, to be specified later on, and ignore any
FIG. 1. (a) The copepod velocity relative to temporal se-
other activity induced by light, food, or chemistry (e.g.
quencewithmultiplejumpsoccurredinresponsetostimulus.
(b) Several jumps superposed by a shift, taking as reference pheromones). On the mechanical side, we assume that
time that are associated with their peak position. Almost copepods are small enough that their centre of mass can
all of the jumps decay exponentially. (Inset) The probability be considered a perfect fluid tracer in a flow, except for
density function (PDF) of thejump intensity. the time when a jump event takes place. In hydrody-
namic terms this means that copepods are assumed to
be rigid, homogeneous, neutrally buoyant particles with
fromtheoneduetothepresenceofavariableflowshear- a size which is of the order of the dissipative scale of the
rate. This aspect will therefore not be taken into ac- flow. Gravity force has no role in producing accelera-
count in the model presented in the next section. We tion or torque. Only the drag force effect is taken into
plantoinvestigateinter-jumpstatisticsmorecarefullyin account during the jumps. Finally, copepods are cou-
the future, whenexperiments withmechanicallyinduced pledtothe fluidina one-wayfashion,they reactandare
stimulus may be available. carried by it, but they do not modify the surrounding
flow, copepods-copepods interactions are also neglected.
Adding all together the above hypothesis the LC equa-
tion of motion is as follows:
B. Model equation for copepods dynamics
x˙(t)=u(x(t),t)+J(t,t ,t ,γ˙,p) (1)
i e
In this section we introduce a simple model system of
copepod’sdynamics. Thisrepresentationisbasedonthe
ideathatthecopepod’strajectoriesinafluidcanbemim- where u(x(t),t) is the velocity of the carrying fluid at
icked by properly defined active particles. Similar mod- time t and position x(t) and where J is an added ve-
elshavebeensuccessfullyemployedforthedescriptionof locitytermthatdescribestheactivebehaviour(jump)of
thebehaviourofphytoplankton,suchaschlamydomonas thecopepod. J(t,t ,t ,γ˙,p)isafunctionoftimet,itde-
i e
[16, 57, 58] both in laminar and, more recently, in tur- pendsalsoonaninitialandafinaltimet andt ,onflow
i e
bulent flows [14, 15, 22]. Copepods, and zooplankton in shear rate value γ˙ and on orientation vector p. If cope-
general, display higher complexity compared e.g. to al- podsaretakento besphericalinshape,their orientation
gae species because of their higher motility. The model dynamics is given by:
relies both on biological and hydrodynamical assump-
tions. First,we assume thatcopepods respondalwaysin p˙(t)=Ω p(t) (2)
·
4
ocean turbulence: The dissipative length η = (ν3/ǫ)1/4,
where Ω is the fluid rotation rate antisymmetric tensor, time τ = (νǫ)1/2 and velocity u = (νǫ)1/4. The or-
η η
definedasΩ =1/2(∂ u ∂ u ). Amoregeneralformof der of magnitude estimate as from Ref. [61] for these
ij i j j i
−
the equation (2), valid for axisymmetric ellipsoidal par- quantities arereportedin table I. According to the same
ticles, is as follows: authorsthetypicalTaylor-scaleReynoldsnumberRe in
λ
the ocean can reach values up to (102).
O
Given that the typical size of copepods is of the order
p˙(t)= Ω+ α2−1 pT(t) p(t) p(t) (3) of millimetres, it is clear that the relevantflow scales for
(cid:16) α2+1(cid:0)S− ·S· (cid:1)(cid:17)· their dynamics are close to the Kolmogorov scale or be-
low in turbulence [29]. When the LC model is recast
whereα l/distheaspectratiooftheellipsoidsgivenby
≡ in a dimensionless form in terms of these scales we get
the ratio of length (l) to diameter (d), which is typically
three dimensionless groups of parameters: τ /τ , u /u
around 3 for Eurytemora affinis. The above equation J η J η
andτ γ˙ . Theseparameters,togetherwiththe flowRe
was first proposed by Jeffery, and its full derivation is η T λ
fully specify the working conditions (or tuning) of the
detailedin[59]. Itsphenomenologyinturbulentflowshas
copepods-in-turbulence model.
been investigatedmore recently in [54]. Notice that here
In this study we take as reference for the energy dissipa-
we designate by the fluid deformation rate symmetric
S tion rate the value ǫ = 10−6m2s−3, and by taking into
tensor as S = 1/2(∂ u +∂ u ) and the shear rate is
ij i j j i
account the dimensional values estimated for the cope-
thendefinedasγ˙ =√2 : . Wenote thatthe factthat
the jump term is assumSedSto depend on γ˙ represents a pods jump intensity uJ and jump decaying time τJ , the
ratiosu /u =93.9andτ /τ =0.00887canbededuced
generalizationtothe3DgeometryofKiørboe’sempirical J η J η
fromthesimilarityanalysis. Thistellsusthatinordinary
findings[32]. Forthejumptermweproposethefollowing
turbulenceconditionsthecopepodspossessanalmostin-
functional form:
stantaneous reaction, since their response time is about
ti−t one hundredth of the smallest scale of turbulence. On
J(t,ti,te,γ˙,p)=H[γ˙(ti)−γ˙T]H[te−t]uJe τJ p(ti) the opposite the velocity reached during a jump is of a
(4)
magnitude that is comparableif not higher to the one of
where H[x] denotes the Heaviside step function, γ˙ is
T turbulent velocity fluctuations. Finally, we note that we
a threshold value of the shear rate, u and τ are two
J J donothaveanyexperimentalguessforthe magnitudeof
characteristicparameterscharacterisingthe jump shape,
γ˙ , therefore the value τ γ˙ is a free parameter of our
T η T
its velocity amplitude (u ) and duration τ respectively.
J J model.
The first H step function models the fact that a jump
can begin only when the shear rate is above the given
Parameter Unit Range This study
threshold value, while the second step function accounts
ν m2s−1 10−6 10−6
for the factthat the jump time spanis finite. The initial
∼
and final time of a jump are defined as: ǫ m2s−3 10−8 10−4 10−6
ti =t if (γ˙(t)>γ˙T) (t>te) (5) η m 3 10−3 3 10−4 10−3
∩ × ×
te =ti+cτJ =ti+log(102)τJ (6) τη s 10 0.1 1
u ms−1 3 10−4 3 10−3 10−3
η
In other words we assume that a jump can not begin × ×
Re – (102) 80
if a previous jump has not finished (t > te) and that a λ O
jump terminates when its amplitude has decreased to a
TABLEI.Referencepropertiesoftheoceanturbulentflowas
negligible level, here taken as one percent of the initial
from [61]. ǫ is themean turbulent energy dissipation rate, η,
amplitude, i.e. |J(te)|=10−2|J(ti)|. τη anduη aretheturbulencespace,timeandvelocitydissipa-
tive scales. Re is the Taylor-scale based Reynolds number.
λ
Their approximaterange of variability is given together with
the reference values chosen for the similarity analysis in the
C. Model tuning for turbulent flows
present study.
We now take into account the presence of the oceanic
flow environment surrounding the copepods. The prop-
erties of oceanic turbulence relevant for our work have
been studied, among others by MacKenzie et al. [60] D. Numerical implementation of the LC model and
of the turbulent flow simulation
and Jimenez [61]. In these surveys it was observed that
themeanvalueoftheturbulentkineticenergydissipation
rate, ǫ=2ν : , varies from about 10−8m2s−3 in open The copepods-in-turbulence model system is conve-
oceanto 10−S4mS2s−3 in coastalzones (althoughit is also niently implemented via an Eulerian-Lagrangian ap-
sensitivetothewindspeedconditionsandonthedepth). proach,meaningthatthe trajectoryx(t)ofeachindivid-
The value of ǫ along with the kinematic viscosity of sea ual copepod is computed by means of Lagrangiantrack-
water, ν, allow to estimate the Kolmogorov scales of ingmethodappliedtoeq. (1)[62,63],whilethefluidflow
5
is obtained by solving the field equations of incompress- tionwasstartedandparticleswereletdisplaceforabout
ible fluid-dynamics, i.e. Navier-Stokes equations, in tur- 2eddyturnovertimes,afterthatduringthefollowing 2
∼
bulent conditions. All the particles are advancedin time eddyturnovertimesabout10instantaneousdistributions
using Adams-Bashforth method with a time step equal of LC particles were saved for analysis. Copepods are
to δt=1.4 10−3τ , the same time step as for the inte- modelledassolidsphereparticles,andorientationvector
η
×
gration of the Navier-Stokes equations. Such a choice of affected by fluid rotation rate (eq. 2), unless otherwise
timestepshallalsosatisfytheconstraintδt τ . noted. For comparison a set of passive fluid tracers are
J
≪
A Direct Numerical Simulation (DNS) approach was also included in all our simulations.
used to solve the Navier-Stokes equations for homoge-
neousisotropicturbulencebymeansofapseudo-spectral
method:
A. Single Point Statistics
∂ u+u u= p/ρ+ν∆u+f (7)
t
·∇ −∇ In order to see how the LC dynamics in turbulence dif-
fer from that of a fluid tracer, we first address the ve-
where u(x(t),t) is the incompressible ( u = 0) fluid
∇· locity single point statistics. The PDF of the absolute
velocity field, p is the pressure, ν is the kinematic vis-
value of single component velocity for the copepods, i.e.
cosity and ρ is the fluid density. The f is the forcing
x˙ , is shown in figure 3. Tracers, the particles which
which is applied on large scales to sustain the statisti- | i|
move along the streamlines, agree with a Gaussian dis-
cally stationary turbulence. The solution domain is a
tribution, while for copepods a slower decaying tail is
cube of length L = 2π with N3 = 1283 grid points,
found. This deviation becomes more pronounced at in-
subject to periodic boundary condition. Aliasing error
creasing the jump intensity for a given threshold value
is controlled by omitting the wavenumber larger than
of the shear rate, as shown in figure 3(a). It also ap-
k = 2/3 (2πN/L), to reach the Taylor Reynolds num-
× pearsthatlowjumpintensitiesu <10u arenotstrong
beroftheflowRe =√15u2 /(νǫ)1/2 80whereu J η
λ rms ≈ rms enough to make effective changes on the copepods PDF.
is the single component root mean square velocity fluc-
On the other hand, increasing the threshold value of the
tuation. k η > 1.4, in which k = N/3 and η is
max max shear rate leads to fewer jumps, therefore in this case
the Kolmogorov length scale, assures that small scales
copepods behave almost like tracers. Their deviation in
structures are well resolved.
velocitydistributionfromthe Gaussian,indeedincreases
by decreasing the shear rate threshold value as can be
seen in figure 3(b).
III. RESULTS AND DISCUSSION
The general trend of the observed deviation from Gaus-
sianitycanbepredictedbymeansofthefollowingproba-
As mentioned above, the LC model is characterised by
bilistic model. We suppose that the instantaneous single
three control parameters: the jump intensity u , the de-
J cartesian component velocity of LC particles can be ap-
caying time of the jump τ and the shear rate threshold
J proximated by the sum of three statistical independent
valueγ˙ ,whichareconvenientlypresentedindimension-
T randomvariables. Thefirstvariableaccountsforthetur-
less form in terms of turbulence dissipative scale units.
bulent velocity field contribution, therefore it is a Gaus-
Since the LC model is just one-way coupled to the fluid,
sian with zero mean and same standarddeviation as the
inthenumericswecanperformsimultaneoussimulations
onemeasuredintheDNS.Thesecondandthethirdvari-
of several families of copepods in the same turbulent
ablemimic respectivelythe jumpdirectionanditsinten-
flow, where each family is characterized by the triplet
sity: we assume that the orientation is random uniform
[u /u ,τ /τ ,γ˙ τ ].
J η J η T η inthesolidangleandthatthejumpshappenuniformlyin
In agreement with the experimental observation we al-
time. OnecanobtaintheresultingPDFfortheLCparti-
ways keep fixed the decaying time of the jump to the
cle velocityfromthe convolutionofthe three elementary
value τ /τ =10−2, while the other parameters are var-
J η PDFsassociatedtothethreedescribedrandomvariables.
ied independently in the ranges u /u [1,400] and
J η ∈ The resulting density distribution function when com-
τ γ˙ [0,4]. Note that if γ˙ = 0, according to the
η T ∈ T pared to the LC measurements at low threshold value
model, all the particles will jump in a synchronous way.
τ γ˙ =0.21 (i.e. when copepods jump very frequently),
η T
In order to avoidsuch an unphysical feature, the time t
e shows an overall qualitative agreement with a slight de-
for each particle is initialised by a random variable with
viation in the tails (see Fig. 3(b)).
homogeneous distribution in the interval [0,log(102)τ ].
J Such a discrepancy comes from the fact that in reality
Weperformaseriesofsimulationswithmultiplefamilies,
the jump directions develop some correlations with the
with about 2.56 105 particles per family1. The simula-
× underlying flow, via Eq. 2, while the probabilistic model
neglectsit. Onecanmakeuseofthe approximateproba-
bilistic modelto estimate the averagefractionofparticle
1 In physical dimension this corresponds to a number density of performing jumps as a function of the shear-rate thresh-
O(1) LC particles per cm3, a density comparable to the one old value. This is done by introducing an adjustable pa-
foundforrealcopepods estuarinewater. rameter accounting for the probability that a given par-
6
ticle is actually jumping, andby fitting the model to the 1
10
PDF curves. Figure 4 shows the fitted predictions ob- 102
tofiafgiunthreeeddpwisripothblaaysbusicltihhsteaicipnmrfeoorcrdeededlu)jr,uewm(hwpihlpeiectrhhcecenoitnnafisgreemtaostfhaetfhuveanlsciatdimiotnye |)η 1100-01 Inferredmp percentage111000-011
tonhfeentthjiueamlsdhpeeicnargreaprsaaetreatiscthlγ˙erTefrsisahcortalidiosnevdai.lsuFaeor.oruWthneedvo5ab0ls%ueer.vτeηγ˙aTn=ex0p.o5-, PDF(| ˙x/ui 1100--32 ju10-120 00%.5 o 1f j u1.m5τpη 2 -γ˙ Tτ 2η.55 γ˙0T T%3 =r a -30c .0.e52.r51s
We finally observe that the shape of the PDF displayed
25% - 0.7
-4
by the LC model, is also in qualitative agreementwith a 10 6.25% - 1.2
recently published experimental study [43], despite the 1.25% - 1.77
-5
10
fact that the experiment has been performed in low
0 50 100 150 200
Reynolds number conditions (up to Re 30). What
has not been reported yet in experimenλta≃l studies is a | x˙i/uη|
quantification of the three-dimensional spatial distribu-
tion of copepods in turbulence. We do this in the next FIG. 4. Probability densityfunction of absolute valueof sin-
sectionbymeansofafractaldimensioncharacterization. gle component velocity |x˙ /u | for the copepods at constant
i η
jumpintensityu /u =100fordifferentshearratethreshold
J η
values. Fitted PDF curves correspond to the percentage of
jump of copepods. (Inset) Deduced percentage of jump as a
(a) 100 function of theshear rate threshold value τ γ˙ .
Tracers η T
Gaussian distribution
10-1 τη γ˙T = 0.7, uJ/uη = 1
0.7, 50
|)η -2 0.7, 100
u 10 0.7, 200
˙x/i B. Correlation Dimension Analysis
F(| 10-3
D
The distribution of the LC particles is illustrated by fig-
P
-4 ure5,whereweshowtheinstantaneousparticlepositions
10
in two-dimensional slices of thickness η, visualising at
∼
-5 the same time the values of shear-rate of the carrying
10
0 50 100 150 200 flow. Contrary to fluid tracers, LC particles are non-
| x˙i/uη| homogeneouslydispersedinregionswhereturbulencein-
tensityisbelowthegivenshear-ratethreshold,according
to the model. In the panels of figure 5, we also highlight
(b) 100 theγ˙ valuesbycontourlines,wenamerespectivelycom-
Tracers T
Gaussian distribution fort and alert regions the locations which are below or
-1
) 10 RanduoJ/muη j u=m 1p0s0, ,u τJη / γ˙uTη == 01.0201 sapboonvdestthoesγ˙e fi=xe0d.3γ˙5Tτ−va1l,utehse.aIlnertFrigeg.io5n(ai)s,twhehidcohmcionrarnet-
|η -2 100, 1.77 T η
u 10 one. In this situation the great majority of LC particles
DF(| ˙x/i 10-3 110000,, 03..79 aavreailjaubmlepicnogmbfourtt tishleayndms.anTihfeisstlmyafyaibletodureeabcohththteoftehwe
P fact that islands are small and that they are short lived:
-4
10 one shall bear in mind the interplay between space and
time in this problem. The panel (b) shows a condition
-5
10 wherecomfortandalertregionsareequallyprobable. We
0 50 100 150 200
notice a pronounced aggregationof particles in the alert
| x˙i/uη| areas surrounding the comfort regions, while the latter
areefficientlyevacuated. Finally,thepanel(c)illustrates
whathappens when the alertbehaviouris triggeredonly
FIG. 3. Probability density function of absolute valueof sin-
by few extreme shear rate filamentary regions. The LC
gle component velocity |x˙ /u | for the copepods. (a) at con-
i η particlesmanagetoavoidthemquiteefficientlybutinthe
stant threshold value τ γ˙ = 0.7 and different jump inten-
η T
overall picture they seems to be mostly homogeneously
sities. Gaussian distribution is a statistic distribution here
with themeasured root mean squarevelocity of theEulerian distributed. (See Supplemental Material [64] for the 2D
fieldasthestandarddeviation. (b)atconstantjumpintensity visualisation of copepod’s motion in turbulent flow at
u /u = 100 for different shear rate threshold values. Ran- τ γ˙ =0.92.)
J η η T
dom jumps correspond to the expected velocity distribution InordertobetterquantifythepatchinessoftheLCparti-
when randomly oriented jumps occur uniformly in time on cleswecomputetheirCorrelationDimension(D ),which
2
top of theturbulent velocity field. is a measure ofthe dimensionality ofa set ofpoints. Ac-
7
cordingtotheGrassbergerandProcacciaalgorithm[65],
4.5 theD isdefinedasthescalingexponentoftheprobabil-
2
ityoffindingapairofparticleswithaseparationdistance
4
lessthanr :P (X X < r) rD2 asr 0. Inother
2 2 1
3.5 words if | − | ∝ →
3
2
C(r)= H(r X X ) (8)
2.5 γ ˙η N(N 1)X −| i− j|
τ − i<j
2
1.5 decreases like a power law, then D2 = rli→m0lolgoCgr(r). Fig-
1 ure 6 shows the D2 value in the two dimensional param-
eter space composed by the intensity of the jump and
0.5
the shear rate threshold value. The clustering (D < 3)
2
0 is discernible when the prescribed shear rate threshold
(a) τη γ˙T = 0.35 vaatlaureouisnldes0s.5thτa−n1.2O.8nτηt−h1e,oatnhderithaisndmwaxeimobasle,rDve2a≃sa2t.u3-,
η
rationofclustering as u is increased. In orderto better
4.5 J
appreciate these two features, i.e. the minimum with re-
4
spect to τ γ˙ and a saturation as a function of u /u ,
η T J η
3.5 two two-dimensional cuts of the D2(γ˙T,uJ) surface are
shown in Fig. 7.
3
2.5 γ ˙η
τ 3.1
2
3 3.1
1.5 2.9 3
2.8 2.9
1 DD22 2.7 2.8
2.6 2.7
0.5 2.5 2.6
2.5
2.4
2.4
0 2.3 2.3
(b) τη γ˙T = 0.92 050
100
150
200
4.5 uJ/uη 250
300
4 3 540 04050 4 3.5 3 2.5 τ2η γ˙T 1.5 1 0.5 0
3.5
3 FIG. 6. Correlation dimension D2 of copepods as a function
of jump intensity u /u and threshold valueτ γ˙ .
2.5 γ ˙η J η η T
τ
2
One may wonder why there is an optimum and what is
1.5
its physical meaning. Copepods are prone to jump in
1 order to escape from regions of alert (γ˙ > γ˙T), to reach
regionswhereγ˙ <γ˙ ,thereforethechanceforajumpto
0.5 T
be successful (assuming it to be randomly oriented) de-
0 pendsonthesizeofthecomfortregion,inotherwordsto
(c) τη γ˙T = 1.77 the volume, Vγ˙<γ˙T. Onthe other hand,clustering would
be maximum if we have numerous successful jumps, and
FIG. 5. Patchiness of the copepods from the simulations. obviously the number of jumps depends on . This
Vγ˙>γ˙T
Shading shows the instantaneous field of the absolute value implies that copepods clustering is expected to be pro-
of the shear rate of the Eulerian field, (a) distribution of the portional to . Now, substituting the vol-
copepodswithuJ/uη =250andτηγ˙T =0.35,(b)distribution ume of comfoVrγ˙t<abγ˙Tle·rVeγg˙>ioγ˙nTs with leads to
of the copepods with u /u = 250 and τ γ˙ = 0.92, (c) Vtot −Vγ˙>γ˙T
J η η T ( ). One direct consequence is that
distribution of the copepods with uJ/uη = 250 and τηγ˙T = Vthγ˙e>cγ˙lTus·teVrtiontg−wVouγ˙l>dγ˙Tbe maximumwhen = /2.
1.77. Contour lines are traced at the corresponding value of Vγ˙>γ˙T Vtot
This can explain the existence of the optimum of D as
γ˙ on each panel. 2
T
a function of τ γ˙ as shown in figure 7(a) as well as its
η T
trend as a function of γ˙ . Note however that this argu-
T
ment is based on the simplifying assumption that there
8
is no correlation between the orientation of a LC parti- tion and shows a peak for the same value of γ˙ (in-
T
cle at jump and its position respect to the comfort area, set of Fig.8). The prediction that clustering varies as
andalsoitneglectsthespatialstructureoftheshearrate ( ) is in qualitative agreement with
Vγ˙>γ˙T · Vtot −Vγ˙>γ˙T
field. the observed trend, it is in quite good agreement in the
large γ˙ regime, however it fails to capture the correct
T
(a) 3.2 valueatwhichthemaximumappears,givingτηγ˙T =0.85
instead of 0.5. Finally, we note that the case of maximal
3 clustering at D 2.3 corresponds to a condition where
2
≃
the LC particles concentrate in nearly two-dimensional
2.8
uJ/uη = 1 sheets which envelop the alert regions (as can be also
D2 2.6 100 inferred from the visualisation in Fig. 5(b)).
250
2.4
1
2.2 LC
Tracer
2 0.8 1
0 0.5 1 1.5 2 2.5 3 3.5 4
τη γ˙T > / Ttot 0.6 000...468 TraceTrr*a(c1e-rt r-a cLeCr)
T
(b) 3 .32 τη γ˙T = 001.3..745 <Tγγ˙ > ˙ 0.4 0 .02 0 0.5 1 1.5 2 2.5 3 3.5 4
0.2
2.8 τη γ˙T
D2 2.6 0
0 0.5 1 1.5 2 2.5 3 3.5 4
2.4 τη γ˙T
2.2
FIG. 8. Time fraction spent in alert regions by tracer and
2
0 50 100 150 200 250 300 350 400 450 copepods as a function of the threshold value. Inset: differ-
uJ/uη ence of time fraction between tracers and LC particles, and
prediction based on Vγ˙>γ˙T/Vtot measurement.
FIG. 7. Lateral view of correlation dimension of the cope-
We can offer a qualitative physical explanation for the
podsasafunctionofthejumpintensityu /u andthreshold
J η observedD saturationforhighvalues ofu atfixedγ˙ ,
valueτ γ˙ . Errorbarsindicatetherangeofvariabilityofthe 2 J T
η T
(figure 7(b)). The argument is as follow: one may ex-
measurements from 10 independentparticle snapshots.
pectthatthereisclusteringifthetimetoescapefroman
alert region is less than the lifetime of such a region:
How can we determine the value of γ˙ for which the con-
τ < τ . The former time can be estimated as
dition of = /2 occurs? One possibility is to escape γ˙T
performaVnγ˙E>uγ˙lTerianVmtoetasurementoftheγ˙(x,t)fieldover τescape =lγ˙T/uJ,wherelγ˙T isthetypicalsizeofthealert
region characterised by a shear-rate γ˙ > γ˙ . This im-
T
spaceandtime. Anotheroptionistolookatthefraction
pliesthatLC particlesformclustersandthe D measure
2
oftimespentbytracersinalertregions,T /T (with
T thetotaltimeofthemeasurement). Sγ˙i>nγc˙Tetratoctersex- is lead to saturate to a constant value if uJ > lγ˙T/τγ˙T.
tot Thislatterratiocanbethoughtasathresholddependent
plore evenly all the region of the flow this is equivalent
escape velocity u = l /τ . From the correlation di-
to measure the volume ratio / . In particular γ˙T γ˙T γ˙T
Vγ˙>γ˙T Vtot mension measurement this escape velocity is estimated
in order to increase the statistical sampling we look at
to be of the order of 100u , i.e. of the order of the large
η
the globalmean value T /T where the averageis
h γ˙>γ˙Ti tot scalevelocity,withaweakdecreasingtrendatincreasing
over the total number of particles (N ):
tot γ˙ .
T
We finally observe that when the flow field associatedto
1 Ntot Ttot
the Lagrangian particles, v = x˙, displays a weak com-
T = H(γ˙ (t) γ˙ ) dt. (9)
h γ˙>γ˙Ti Ntot Xi=1 Z0 i − T pressibility, it can be shown [15, 66] that D2 depends on
the flow divergence by the relation D = 3 c v 2
2
− h∇· i
The plot in figure 8 shows the trend of T /T as with c a proportionality constant and angular brack-
h γ˙>γ˙Ti tot
function of γ˙ both for tracers and LC particles. It con- ets denoting time and space average. If this argument
T
firms thatcopepods resideless inalertregionscompared is applied to the LC model we observe that the diver-
to tracers. Moreover the difference among the two time gence can be different form zero only at the interface
fractions can be used as an alternative clustering indica- between comfort and alert regions. This is because in
tor. It has in fact a similar trend as the D (γ˙ ) func- comfort regions ( v = u = 0) and in alert regions
2 T
∇· ∇·
9
( v = u+ J = 0, as we can safely assume
(a) 3.2
∇ · ∇ · ∇·
thejumptermtobespatiallyconstant). Attheinterface
however,thechangefromthefluidvelocityintensityuto 3
u+u hasaspatialtransitionscaleroughlyproportional
J 2.8
to u log(102)τ which leads to a non-null divergence.
J J
·
This explains the LC accumulation that we observe in D2 2.6
correspondenceof the alert/comfortinterfaces, whichef-
2.4
fectively acts as sink or source term of the LC velocity
Random direction
field (see in particular the central panel of Fig. 5). By 2.2 Jeffery, α=1
following this line of reasoning, one can guess that the Jeffery, α=3
minimum value of D will correspond to the case where 2
2 0 0.5 1 1.5 2 2.5 3 3.5 4
the surface of alert/comfort interface is maximum (and
not of volumes, as stated above). This has clearly a de- τη γ˙T
pendence on the threshold γ˙ and much less, if any, on
T (b) 100
uJ. Despitethequalitativeagreementofthisobservation Tracers
wcointhfiromuritnuqmuaenritcitailvreelysuilnts,thweewheaavkelynoctombepernesasbibleleyleitmtiot 10-1 RandJoemff edriyr,e cαt i=on 3, ,ττηη γγ˙˙TT == 00..8855
of the LC model. DF(| ˙x/u|)ηi 1100--32 JefferGya, uαs s=i a1n, τdηi sγ˙tTri =bu 0ti.o8n5
P
C. Particle Orientational Dynamics
10-4
What is the importance of particle orientation for the 10-5
non homogenous distribution of particles? The effect of 0 50 100 150 200 250 300
the geometrical aspect ratio of the particles, together | x˙i/uη|
with the direction of their jump on the fractal dimen-
sionareaddressedhere. Thefluiddeformationratesym-
FIG. 9. (a) Effect of the aspect ratio and direction of the
metrictensorS comesintoplaybymodellingcopepods
ij jumponthefractaldimension,(red)Copepodassolidsphere
as elongated particles with aspect ratio equal to 3 (e.g.
particles, their direction of the jump is random in the solid
the relevant aspect ratio for Eurytemora affinis cope- angle; (blue) Copepod as solid sphere particles with an ori-
pod). Its effect on the jump direction selection leads entation and (green) as elongated particles. Both jump in a
toenhancedclusteringofthe particlesforjump intensity directionfollowing theJeffery’sequation. (b)PDFoftheab-
u /u =250. Copepods can also jump in random direc- solutevalueofthesinglecomponentvelocityfor(red)random
J η
tion in the solid angle independently from the rotation direction case, and (green) Jeffery’s case with aspect ratio of
rateanddeformationrateoftheEulerianfield. Lessclus- 3. Allcases are computed at uJ/uη =250.
tering in this case is logical since the jumping direction
hasnorelationwiththefluidflow. Thesebehaviourscan
befoundinmoredetailsinfigure9,whereweaddressthe tononhomogeneousspatialdistributions ofactiveparti-
influence of jump direction on the PDF of the copepods cles. Thisclusteringmechanismshoweveriseffectiveonly
velocity. whenthereactionthresholdisclosetovaluesoftheorder
of τ−1 in a very narrow range. The fact that the range
η
is narrow is ultimately linked to the intermittent distri-
bution of the turbulence dissipation rate [67]. We have
IV. CONCLUSIONS AND PERSPECTIVES
shownthatclusteringapproachesitsmaximumwhenthe
thresholdratevalueγ˙ roughlydivide the shear-rate(γ˙)
T
Inthis study we haveconsideredaLagrangianmodelfor spatial field in equal volume regions. Since this mech-
activeparticles. Themodelistrimmedinawaytorepro- anisms mainly depends on the average value of small-
ducesomedynamicalfeaturesexperimentallyobservedin turbulence scales rather than on their fluctuations we
themotionofcopepodsinstillwater. Itsmaincharacter- expect it to have a weak dependence on the Reynolds
isticsisthepossibilitytolocallyacquireanextra-velocity number of the turbulent flow. A second implication of
(jump) in response to a variation of the fluid flow condi- the model is that for any given shear-ratereaction value
tionssurroundingtheparticle. Thedirectionofthejump γ˙ there is a maximal intensity jump velocity beyond
T
isruledbythehydrodynamicsofsmallneutrally-buoyant which clustering can not be further increased. Finally,
particles. The Lagrangian model has been coupled to a theanalysisofthecorrelationdimensionsuggeststhefor-
turbulentdevelopedflowdescribedbytheincompressible mationoflocalquasi-bidimensionalclustersenclosingthe
Navier-Stokes equations. non-permitted flow regions. From a physicist viewpoint
We have shown that jump escape reaction from spatio- weremarkthattheclusteringmechanismatworkintur-
temporal events characterised by high shear-rate leads bulencefortheLCmodelisdifferentformtheoneshown
10
in other model systems of particulate active matter. For enhancement of the mating contact rate. Of course such
instance the clustering observed for motile algal cells in copepodconcentrationcouldalsoattractpredators. Due
turbulence is given by the gyrotactic effect, which is a to different trade-offs, each copepod species may have
non-isotropiceffectinducedbythepresenceofthetheex- an optimal jump behaviour in response to turbulence.
ternal gravity field [22]. On the opposite, the LC model Forexamplethe copepodEurytemora affinis usedinour
discussed here is isotropic but it is non-homogeneous in experimental section is an estuarine species adapted to
space (it depends on the local value of the shear-rate). maintainthe bulk ofits populationin a salinity gradient
Wehavetestedthefactthatclusteringalsoappearswhen in highly turbulent conditions [76, 77]. By using high
LC particles are made sensitive to other flow quantities frequency sampling data of all life stages of E. affinis,
such as enstrophy or fluid acceleration. The minimal Schmittetal. [77]confirmedthatthelatedevelopmental
fractal dimension we observed is always above the value stages(mainlyadults)exhibitedactiveverticalmigration
of2,confirmingthe factthatparticlesinthis caseaggre- during the flood. Consequently the population was not
gateinordertocoverthesurfaceoftheforbiddenregions. homogeneouslydistributedinthewatercolumn,asdense
Based on these observations we do not expect that such patchesareobservedduring shorttime windowandnear
clustering processes could lead to filamentary like clus- the bottom [76].
ters, D 1, as the ones observed for microbubbles in Our model can be improved in the future to test such
2
≃
turbulent flows. Another notable result is the negligible situation with tidally induced turbulence in shallow es-
impactoftheparticleorientationaldynamicsontheclus- tuaries where copepods can use their jump abilities to
tering. This is likely to be linkedto the limited duration simplyavoidtobeflushedouttheiroptimalhabitat. This
ofjumps(notethathereτ τ ),butmightbecomeim- couldleadtotheidentificationofsomeoptimalclustering
J η
≪
portantforlongerjumps,particularlyinthemodellingof strategy that may be in relation with the dome-shapes
largermotile plankton. The negligibleimpactoforienta- proposedearlier,onpurelyspeculativeintuitions[78,79].
tion for the case examined here, suggests the possibility The presented LC model can also be improved by refin-
toformulateaccurateeulerianmean-fieldparticlemodels ingthejumpingprotocolinordertotakeintoaccountthe
based on the introduction of a space-dependent effective factthatthetemporalsequenceofjumpsincopepodsoc-
diffusivity(κ)whoseamplitudemaybelinkeddirectlyto cursinfastsequences(bursts)interposedtoinactivemo-
jumpshapeparameters,viaadimensionalrelationofthe ments. Another possible direction of research concerns
type κ u2τ . the investigationofthe impactof aspatialradiusof per-
∝ J J
Fromamorebiologicalperspective,althoughbehavioural ception for copepods to react to turbulent shear. This
mechanisms leading to clustering had been already sug- may produce a smoothing or a delay in the perceived
gested in the past, such as the formation of patches turbulent signal.
through swimming against the flow [68], the possibil- Aknowledgments: The authors acknowledge support
ity of cluster formation by escape jumps in a no-mean from European COST Action MP1305 “Flowing Mat-
flow situation was never reported before. As discussed ter” and the “GIP Seine-Aval ZOOGLOBAL project”.
in [38], clustering of copepods has an ecological impor- DominiqueMenuisacknowledgedfortechnicalhelpcon-
tance: an effect may be to strongly increase the contact cerning the copepod experiment devices, and Regis Sion
rate with mates, and hence improve the reproduction. foradvicesandhelpontheuseofthefastcamera. Marion
Indeed several models have been proposed to express Roussin and Gu´enol´e Alizard are thanked for help with
copepod contact rates in turbulence [69–71], reviewedin copepod sorting. We acknowledge discussion with Mas-
[72]. In case of clustering, the contact rate is strongly simo Cencini and Guido Boffetta during the workshop
increased [38, 73–75]. The clustering which would re- Micoroorganisms in Turbulent Flows, recently held at
sult from a behaviour of predator avoidance (a reaction LorentzCenterinLeiden. H.A.issupportedbythePhD
to turbulent shears similar to predator’s signals) would grantforinterdisciplinaryresearch“AllocationPresident
have as side-effect a positive consequence with a strong 2013” of the University of Lille 1.
[1] H.C.BergandL.Turner,“Movementofmicroorganisms [4] E.Lauga,“Bacterialhydrodynamics,”AnnualReviewof
in viscous environments,” Nature 278, 349–351 (1979). Fluid Mechanics 48, 105–130 (2016).
[2] R. Di Leonardo, L. Angelani, D. Dell’Arciprete, [5] J. Teran, L. Fauci, and M. Shelley, “Viscoelastic fluid
G. Ruocco, V. Iebba, S. Schippa, M. P. Conte, response can increase the speed and efficiency of a free
F. Mecarini, F. De Angelis, and E. Di Fabrizio, swimmer,” Physical Review Letters 104, 038101 (2010).
“Bacterial ratchet motors,” Proceeding of the National [6] J. P. Hernandez-Ortiz, C. G. Stoltz, and M. D. Gra-
Academy of Science 107, 9541–9545 (2010). ham, “Transport and collective dynamics in suspensions
[3] A. W. Visser and T. Kiørboe, “Plankton motility pat- ofconfinedswimmingparticles,”PhysicalReviewLetters
terns and encounter rates,” Oecologia 148, 538–546 95, 204501 (2005).
(2006). [7] A.BaskaranandM.C.Marchetti,“Statisticalmechanics
andhydrodynamicsofbacterialsuspensions,”Proceeding
