Journal: Chaos Please provide your responses and any corrections by annotating this PDF and uploading it to AIP’s eProof website as detailed in the Welcome email.
Article Number: 022812CHA
Dear Author,
Below are the queries associated with your article; please answer all of these queries before sending the proof back to AIP.
Article checklist: In order to ensure greater accuracy, please check the following and make all necessary corrections before returning your proof.
1. Is the title of your article accurate and spelled correctly?
2. Please check affiliations including spelling, completeness, and correct linking to authors.
3. Did you remember to include acknowledgment of funding, if required, and is it accurate?
Location in Query / Remark: click on the Q link to navigate to the appropriate article spot in the proof. There, insert your comments as a PDF annotation.
Q1 Please check that the author names are in the proper order and spelled correctly. Also, please ensure that each author’s given and surnames have been correctly identified (given names are highlighted in red and surnames appear in blue).
Q2 In the sentence beginning “In Sec. III B...,” please confirm that “next subsection” refers to Sec. III B.
Q3 We were unable to locate a digital object identifier (doi) for Ref. 25. Please verify and correct author names and journal details (journal title, volume number, page number, and year) as needed and provide the doi. If a doi is not available, no other information is needed from you. For additional information on doi’s, please select this link: http//www.doi.org/.
Thank you for your assistance.
Inhomogeneities and caustics in the sedimentation of noninertial particles in incompressible flows
1 2
GáborDrótos,1,2,a)PedroMonroy,1EmilioHernández-García,1andCristóbalLópez1 Q1 3
1Instituto de Física Interdisciplinar y Sistemas Complejos (IFISC,CSIC-UIB), Campus Universitat de les Illes Balears, E-07122 Palma de Mallorca, Spain
2MTA-ELTE Theoretical Physics Research Group, Pázmany Péter sétany 1/A, H-1117 Budapest, Hungary 4
5 6
(Received 31 January 2018; accepted 30 November 2018; published online xx xx 2018)
7
In an incompressible flow, fluid density remains invariant along fluid element trajectories. This implies that the spatial distribution of non-interacting noninertial particles in such flows cannot develop density inhomogeneities beyond those that are already introduced in the initial condition.
However, in certain practical situations, density is measured or accumulated on (hyper-) surfaces of dimensionality lower than the full dimensionality of the flow in which the particles move. An exam- ple is the observation of particle distributions sedimented on the floor of the ocean. In such cases, even if the initial distribution of noninertial particles is uniform within a finite support in an incom- pressible flow, advection in the flow will give rise to inhomogeneities in the observed density. In this paper, we analytically derive, in the framework of an initially homogeneous particle sheet sediment- ing toward a bottom surface, the relationship between the geometry of the flow and the emerging distribution. From a physical point of view, we identify the two processes that generate inhomo- geneities to be the stretching within the sheet and the projection of the deformed sheet onto the target surface. We point out that an extreme form of inhomogeneity, caustics, can develop for sheets. We exemplify our geometrical results with simulations of particle advection in a simple kinematic flow, study the dependence on various parameters involved, and illustrate that the basic mechanisms work similarly if the initial (homogeneous) distribution occupies a more general region of finite extension rather than a sheet. Published by AIP Publishing.https://doi.org/10.1063/1.5024356
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
Sedimentation of small particles in complex flows is an
25
outstanding problem in science and technology. In par-
26
ticular, the sinking of biogenic particles from the marine
27
surface to the bottom is a fundamental process of the
28
biological carbon pump, playing a key role in the global
29
carbon cycle. A complete understanding of this problem is
30
still lacking. It has been recently shown that despite fluid
31
incompressibility, sedimenting particles moving as nonin-
32
ertial particles in the ocean on large scales show density
33
inhomogeneities when accumulated on some bottom sur-
34
face. Here, we analytically derive the relation between the
35
geometry of the flow and the emerging distribution for an
36
initially homogeneous sheet of particles. From a physical
37
point of view, we identify the two processes that gener-
38
ate inhomogeneities to be the stretching within the sheet
39
and the projection of the deformed sheet onto the target
40
surface. We point out conditions under which an extreme
41
form of inhomogeneity, caustics, can develop.
42
I. INTRODUCTION
43
The sinking of small particles immersed in fluids is a
44
problem of great importance both from theoretical and practi-
45
cal points of view.1,2In an environmental context, the sinking
46
of biogenic particles in the ocean is a fundamental process. It
47
plays a key role in the Earth carbon cycle through the biolog-
48
ical carbon pump, i.e., the sequestration of carbon from the
49
atmosphere performed by phytoplankton via photosynthesis 50
in the surface waters and posterior sedimentation over the 51
oceanic floor.3This is a complex problem, which involves the 52
interplay of biogeochemical processes with oceanic transport 53
phenomena where many open questions remain. In particu- 54
lar, some of these open questions concern the identification of 55
the catchment area (the place near the surface where the par- 56
ticles come from) of a given oceanic floor zone, and which 57
the mechanisms are that lead to the observed inhomogeneous 58
distribution of particles in surface and subsurface oceanic 59
layers4–6 or when collected at a given depth by sediment 60
traps.5,7–10 61
In this paper, we shall describe basic ingredients for 62
the processes that lead to large-scale inhomogeneities in the 63
density of the particles after sedimentation.11These inhomo- 64
geneities emerge as a result of advection of the particles by 65
flows in the ocean. For the range of parameters that is rele- 66
vant for marine biogenic particles, a very good approximation 67
for the equation of motion of the particles,7,9,10,12 as it has 68
been explicitly shown in Ref.11, simply consists of motion 69
following the fluid velocity with an additional settling term. 70
Such an equation of motion, if the fluid flow is incom- 71
pressible, preserves phase-space volume (note that the phase 72
space coincides here with the configuration space). Thus, 73
inertial effects, which have been typically identified as the 74
main causes for particle clustering (also called preferen- 75
tial concentration) in other situations,13–18 cannot explain 76
inhomogeneities in mesoscale oceanic sedimentation. Then 77
the question is which are the mechanisms that lead to 78
1054-1500/2018/28(12)/000000/23/$30.00 28, 000000-1 Published by AIP Publishing.
sedimentation inhomogeneities in the absence of particle
79
inertia.
80
In incompressible flows, density is conserved along tra-
81
jectories, so that inhomogeneities can occur only if they are
82
already present in the initial distribution of the particles, and
83
these initial inhomogeneities are carried over as intact during
84
the entire time evolution, as long as characterizing the concen-
85
tration of particles by a density is an appropriate framework.
86
Note that this fact could already be sufficient for explaining
87
the presence of inhomogeneities: for example, biogenic parti-
88
cles in the ocean are not produced in a uniform distribution,
89
of course.
90
At the same time, one can also argue that parti-
91
cles become uniformly distributed for asymptotically long
92
times in bounded incompressible flows of chaotic nature,19
93
which translates to a uniform particle density, at least when
94
coarse-grained. For localized initial particle distributions in
95
unbounded chaotic systems, a (growing and flattening) Gaus-
96
sian is obtained instead of a uniform density, but such a shape
97
can also be regarded as trivial.
98
However, if the initial distribution is localized, even if
99
being homogeneous within the localized support, it is well-
100
known that complicated structures can be observed before
101
reaching the asymptotic state.20,21 In particular, stretching
102
and folding of the phase-space volume in which the parti-
103
cles are located can, at least when looking at coarse-grained
104
scales, considerably rearrange the density. That is, (coarse-
105
grained) inhomogeneities emerge due to advection, which
106
can be regarded as clustering or preferential concentration. In
107
fact, it is the same process that leads to the above-mentioned
108
asymptotic simplification, but the effect of this process is
109
opposite on non-asymptotic time scales.
110
Preliminary numerical studies in a realistic oceanic set-
111
ting showed that a homogeneous layer of particles (with
112
neglecting the interaction between them) indeed evolves to
113
complicated shapes by stretching and folding while it is sink-
114
ing. As a motivation, Fig.1presents such a direct numerical
115
simulation. It is clear that homogenization or simplification
116
is not reached on the time scale of the sinking process.
117
The example of oceanic sedimentation thus emphasizes the
118
practical importance of the investigation of non-asymptotic
119
time scales in general, which, from a practical point of
120
view, has received relatively little attention in the literature
121
so far (an important exception is the paradigmatic problem of
122
weather forecasting).
123
Beyond stretching and folding during the sinking pro-
124
cess, an important additional ingredient for the emergence
125
of observable inhomogeneities in the density of sedimented
126
particles is the accumulation at the bottom of the domain.
127
This is a time-integration of the density at a two-dimensional
128
subset of the full three-dimensional space, and this results
129
in the translation of the complicated shapes to actual
130
inhomogeneities without any coarse-graining: the conserva-
131
tion of density no longer holds for such time-integrated
132
projections.
133
In this paper, we shall describe in detail how inhomo-
134
geneities in the accumulated density emerge in incompressible
135
flows on non-asymptotic time scales. We will derive the
136
basic mechanisms analytically, and we will investigate the
137
FIG. 1. The positions of particles (projected onto a vertical plane) at dif- ferent times in a realistic ROMS simulation22 of the Benguela zone. The numerical experiment consisted of releasing 6000 particles from initial con- ditions randomly chosen in a square with sides of 1/6◦, centered at 10.0◦E 29.12◦S and 100 m depth. The particles’ trajectories X(t)were calculated from dX/dt=vROMS−Wk, where vˆ ROMSis the velocity from the ROMS model, and W=10 m/day corresponds to the sinking velocity,11pointing in the vertical direction given by the unit vectork.ˆ
properties of these mechanisms in a simplified kinematic flow, 138
in order to focus on the particle dynamics. 139
The main results we achieve are (i) we identify and quan- 140
tify two geometrical mechanisms contributing to the creation 141
of inhomogeneities in the density: the stretching due to the 142
flow and the projection onto the constant depth where the par- 143
ticles accumulate; (ii) we obtain an explicit expression for the 144
density at an arbitrary position of the accumulation level in 145
terms of the trajectories arriving to that particular position; 146
and (iii) in the context of a simplified kinematic flow, we study 147
the dependence on parameters that are generic to the problem: 148
the settling velocity, the depth of the accumulation level, and 149
the amplitude of the fluctuating flow. 150
The paper is organized as follows. In Sec.II, we establish 151
the setup for our analysis. In Sec.III, we obtain the expres- 152
sion for the final density and quantify the above-mentioned 153
two effects leading to inhomogeneities. In Sec.IV, we eval- 154
uate these results in the kinematic flow model. This flow is 155
defined in two dimensions (one horizontal and one vertical), 156
and it may show chaotic behavior. The role of the chaoticity 157
of the flow will be explicitly addressed. In Sec.V, we study 158
the parameter dependence. Finally, in Sec.VI, we summarize 159
and comment on the results. A number of appendices contain 160
the more technical aspects of our paper. 161
II. FORMULATION OF THE SETUP 162
A. Equations of motion 163
In this work, we will consider the motion of particles 164
that follow closely the velocity of the fluid in which they 165
are dispersed, except for the addition of a constant vertical 166
velocity arising from the particle weight. This description is
167
adequate in a variety of circumstances. In particular, it was
168
shown by Monroy et al.11to be the adequate one to describe
169
a wide range of biogenic particles sedimenting in ocean flows
170
with turbulent intensity typical of the open ocean. We revise
171
in the following the arguments leading to that conclusion.
172
The dynamics of spherical particles advected in flu-
173
ids is described, in the small-particle limit, by the
174
Maxey–Riley–Gatignol equation.11,23,24When writing it in a
175
nondimensionalized form that uses the characteristic length L
176
and magnitude U of the fluid velocity field as units of space
177
and velocity, two relevant dimensionless parameters appear.
178
The first one is the Stokes number
179
St= a2U
3νβL, (1)
where a is the radius of the particle,νis the kinematic viscos-
180
ity of the fluid, andβ =3ρf/(2ρp+ρf), withρpandρfbeing
181
the densities of the particle and of the fluid, respectively. This
182
number quantifies the importance of inertia with respect to
183
viscous drag. The second dimensionless quantity is the Froude
184
number, quantifying the importance of inertia with respect to
185
gravity,
186
Fr= U
√gL, (2)
where g is the gravitational acceleration. In terms of these
187
numbers, the dimensionless terminal settling velocity of a
188
particle in still fluid is
189
W =(1−β)St
Fr2. (3)
In complex turbulent flows such as in the ocean, the values
190
of St and Fr vary with scale. Monroy et al.11 showed that
191
for a relevant range of sizes and densities of biogenic parti-
192
cles, St is very small. For example, it takes values11 in the
193
range 10−7–10−1 in wind-driven surface turbulence in the
194
open ocean at the Kolmogorov scale (∼0.3–2 mm), where25
195
typical turbulent velocities are in the range 0.5–3 mm/s. At
196
larger scales, St is still smaller. For example, for mesoscale
197
oceanic motions, Lh=100 km and Uh =0.1 m/s for horizon-
198
tal motion and Lv=100m and Uv=10 m/day≈10−4m/s
199
for vertical motion. This gives St≈10−6for both horizontal
200
and vertical motion. In any case, St is typically very small for
201
the type of particles we are interested in. Under these circum-
202
stances, a standard approach24can be used to approximate the
203
equation of motion for the particle in the limit of vanishing St
204
[seeAppendix Aand Eq.(A1)in particular], provided that the
205
settling velocity W is also small [Eq.(3)].
206
In our ocean situation, the Froude number ranges from
207
10−5 at the mesoscale to a maximum of 10−2 at the Kol-
208
mogorov scale. Thus, the combination St/Fr2, appearing in
209
the settling velocity W , Eq.(3), is within few orders of mag-
210
nitude from 1 and is typically larger than 1. This means, first,
211
that W is always orders of magnitude larger than St, W St,
212
and, second, that W is typically not a small quantity.
213
If W 1 does not hold, the standard approach24 for
214
the small-St approximation is incorrect. In this case, what
215
is appropriate is to take the limit defined by St→0 and
216
Fr→0 with the value of W ∼St/Fr2 remaining constant.
217
Both in this new limit (seeAppendix A) and in the standard 218
approach24with W St, the leading order contribution in St 219
to the equation of motion for the particle is a well-known7,9–12 220
approximation 221
X˙ =v(X, t)≡vfluid(X, t)−Wk,ˆ (4) where we have introduced the notation v for the “velocity 222
field of the particle.” An important feature of the approximate 223
Eq.(4)is the absence of any inertial term. 224
The description(4)would be applicable in other circum- 225
stances beyond the ones described above, but, of course, there 226
would be situations—for example, coastal wave-breaking tur- 227
bulence environments, industrial flows, or (other) cases in 228
which St is not small enough—in which inertial terms will 229
have central importance, with effects that have been studied 230
in recent works.13–18 231
In our paper, we shall restrict our investigations to 232
dynamics of the form of Eq.(4). Additionally, we shall assume 233
|vfluid,z(X, t)|<W for the vertical component of the fluid 234
velocity field, which ensures vz<0 for the vertical com- 235
ponent of the “particle velocity field” v. This assumption 236
excludes the presence of particle trajectories that would be 237
trapped forever to the system, which simplifies the technical 238
treatment of the problem and the interpretation of the phe- 239
nomenology in that the accumulated density at the bottom of 240
the domain is obtained by integrating over finite times. This 241
assumption is reasonable in the above-discussed example of 242
oceanic biogenic particles serving as our motivation.11 243
B. Definitions 244
Let us consider a flow in a d-dimensional space in 245
which we distinguish a “vertical” direction, characterized 246
by the “vertical” coordinate z, and the remaining (d−1)- 247
dimensional subspace, which we call “horizontal,” with the 248
position vector x=(x, y,. . .)=(x1, x2,. . ., xd−1). We ana- 249
lyze the case d =2 in detail, with mentioning d =3 at some 250
points due to its practical relevance, but all results can eas- 251
ily be generalized to higher dimensions, which can be useful 252
when analyzing problems with phase spaces of higher dimen- 253
sionality. The flow is defined by the velocity field v(X, t), 254
X=(x, z)=(x, y,. . ., z)=(x1, x2,. . ., xd)being the position 255
vector in the full space and t being time. vz<0 is assumed for 256
all X and t. 257
We initialize noninertial particles at t=t0 on a given 258
level z=z0 whose density within the so-defined horizontal 259
subspace (a material line and surface for d =2 and d=3, 260
respectively) is described by a “surface” densityσ. We let the 261
particles fall until all of them reach a depth z= −a where 262
they accumulate. We are interested in the resulting horizon- 263
tal “surface” densityσ of the particles measured within the 264
accumulation level. 265
In our notation, a vertical line with a variable in the lower 266
index, |α, corresponds to keeping that particular variable,α, 267
constant, while a vertical line with the declaration of a value, 268
|β=β0, denotes evaluating the preceding expression at the indi- 269
cated value,β0. These two notations can also occur together. 270
As an implicit rule in our notation, when taking derivatives 271
with respect to a horizontal coordinate, all other horizontal
272
coordinates are assumed to be kept constant.
273
III. RELATING THE DENSITY TO PARTICLE
274
TRAJECTORIES
275
The final densityσforming at any position of the accu-
276
mulation level can be related to geometric properties of the
277
flow observable along the trajectory of a particle that was ini-
278
tialized on the initial level z0at t0and that arrives at the given
279
position. If we have more particles, the corresponding densi-
280
ties are to be added. In this section, we first explain that the
281
relation can be given in terms of a special Jacobian, and ana-
282
lyze the formula from some practical aspects. Then we present
283
(for simplicity, taking d=2 in the main text) an intuitive way
284
of building up our formula, which lets us distinguish between
285
the contribution of two simple effects: the stretching within
286
the material line or material surface in which the particles
287
are distributed, and the horizontal kinematic projection (i.e.,
288
a projection that takes the horizontal component of the veloc-
289
ity into account) of the density at the points of arrival at the
290
accumulation level. Each of these two effects is well-defined
291
even in setups in which the other is absent.
292
A. General results
293
Let the endpoint of a trajectory at time t that was ini-
294
tialized at x0 be denoted by f(t; x0)=[ fx(t; x0), fy(t; x0),. . .,
295
fz(t; x0)]=[ f1(t; x0), f2(t; x0),. . ., fd(t; x0)]. The horizontal
296
density at the point where a particular trajectory crosses the
297
accumulation plane z= −a is proportional to the density at
298
the initial position of the given trajectory:
299
σ[t(fz= −a, x0), x0]=σ(t=t0, x0)F[t(fz= −a, x0), x0],
300
(5)
301
where x0 is the d−1 dimensional initial position at t=t0 302
of the particular trajectory within the initial level z=z0,
303
σ(t=t0, x0)is the initial “surface” density at x0, andσ(t, x0)
304
is the horizontal “surface” density at the endpoint, at some
305
time t, of the trajectory that was initialized at x0. The time t
306
of arrival at the accumulation level is unique, since vz<0 is
307
assumed, see Sec.II. This time depends on the vertical posi-
308
tion of the accumulation level, where fz= −a, and also on
309
which trajectory is chosen, which is defined by the initial posi-
310
tion x0. (More generally, an arbitrary time t can be regarded
311
as a function of any final vertical position fzand of the initial
312
position x0, t=t(fz, x0). The relation t(fz, x0)is single-valued
313
because of the assumption vz<0.) In case more than one tra-
314
jectory arrives at the same position within the accumulation
315
level, the corresponding densities are summed up.
316
The total factor,F(t(fz= −a, x0), x0), that multiplies the
317
original density at the starting point of the given trajectory, is
318
the reciprocal of the determinant of a Jacobian:
319
F(t(fz= −a, x0), x0)=det [J(t(fz= −a, x0), x0)]−1, (6) where J is a(d−1)×(d−1)Jacobian:
320
Jij[t(fz, x0), x0]= ∂fj[t(fz, x0), x0]
∂x0i
fz
(7)
for i, j∈ {1,. . ., d−1}. This Jacobian is not a usual one in 321
two aspects. First, it is not a full-dimensional Jacobian, but 322
it is restricted to the horizontal coordinates. In particular, 323
for flows with d=2, it is a scalar. Second, the derivatives 324
with respect to the coordinates of x0 are taken at a constant 325
value of the vertical coordinate fz, and not at a constant time. 326
For this reason, the direct numerical evaluation of Eq. (7) 327
for a given trajectory is not straightforward. Nevertheless, 328
Eqs.(6)–(7)are intuitive in the sense that they give the ratio 329
between the final and the initial values of the “area” of an 330
infinitesimal “surface” element neighboring the given trajec- 331
tory within the material “surface” of particles. For a more 332
rigorous derivation, seeAppendix B. Note that the determi- 333
nant of a full-dimensional Jacobian taken at a constant time 334
is always one for volume-preserving flows. In our setup, the 335
reduced dimensionality and the non-instantaneous accumu- 336
lation process lead to changes in the density, and thus the 337
formation of inhomogeneities becomes possible. 338
We show inAppendix C that the derivatives taken at a 339
constant fz in the Jacobian(7)can be replaced by derivatives 340
taken at a constant time t in the following way: 341
∂fi[t(fz, x0), x0]
∂x0j
fz
= ∂fi(t, x0)
∂x0j
t
− vi[t, f(t; x0)] vz[t, f(t; x0)]
∂fz(t, x0)
∂x0j
t
, 342
(8) 343
for i, j∈ {1,. . ., d−1}. The difference between taking 344
derivatives at constant fz and constant t stems from the fact 345
that different trajectories in the material “surface” reach a 346
given level fz at different times t. From a practical point of 347
view, Eq.(8)can easily be evaluated numerically. 348
Transforming the right-hand side of Eq.(6) in an alter- 349
native way, we learn that it can be obtained from an integral 350
along the given trajectory (as derived inAppendix D): 351
F[t(fz= −a, x0), x0] 352
=exp
⎧⎨
⎩− −a
z0 d−1
i=1
∂
∂fi
vˆi(fz, f) ˆ
vz(fz, f)
fz,f=f(fz,x0)
dfz
⎫⎬
⎭, (9)
353
where f =(f1,. . ., fd−1) denotes the horizontal coordi- 354
nates of the trajectory, andvˆi(fz, f)=vi(t(fz, x0(fz, f)), f(t(fz, 355
x0(fz, f)), x0(fz, f))) for i∈ {1,. . ., d}, i.e., vˆ(fz, f) is the 356
velocity as regarded as a function of the endpoints of the 357
trajectories (instead of the time and the “bare” geometrical 358
coordinates of the domain of the fluid flow). When keeping fz 359
constant, the derivatives taken with respect to the coordinates 360
fi, with i∈ {1,. . ., d−1}, correspond to varying the selected 361
trajectory and also the time t, so that these derivatives are not 362
the instantaneous geometrical derivatives of the velocity field 363
(seeAppendix Dfor a more detailed explanation). By replac- 364
ing the derivatives taken at a constant fzwith those taken at a 365
constant t, we can further transform our formula such that it 366
can be directly evaluated numerically, seeAppendix E. 367
One important aspect of the results presented in this 368
section is that the final density at a given point can be obtained 369
in terms of the initial density at one point (or, at least, a 370
countable number of them) and of the particle trajectory (or 371
trajectories) linking the points: these are all local properties, 372
and no spatially extended information (within the material
373
“surface”) is needed to determine the final density at the given
374
point. In Sec. III B, we rewrite Eqs. (6)–(8) in alternative
375
ways which highlight the contributions from two different and
376
physically intuitive processes.
Q2 377
B. Stretching and projection
378
In this section, we obtain Eqs.(6)–(8)via two physically
379
intuitive steps which correspond to two individual effects that
380
modify the original density. For simplicity, we restrict our-
381
selves to d =2. In order to be able to precisely formulate our
382
considerations, we use a parametric notation for the material
383
line in this section.
384
Let f(t=t0; u)describe a line segment of initial condi-
385
tions at time t=t0 (a material line of particles) embedded
386
in 2 dimensions, parameterized by the arc length u, and let
387
σ(t=t0; u)be the initial density along the line segment at
388
u. Note that the initial line segment need not be horizontal:
389
the results of this section apply for a 1-dimensional initial
390
subset of arbitrary shape, which extends the validity of these
391
considerations to more general setups.
392
Let us denote the image of the initial line segment at time
393
t by f(t; u). The densityσ(t; u)along this image at t in a point
394
whose initial position was characterized by u is given by
395
σ(t; u)=σ(t=t0; u)S(t; u), (10) where
396
S(t; u)= df(t; u)
du
−1. (11) This simply follows from imposing the conservation of mass
397
(i.e., continuity) within the material line of the particles. Note
398
that the density (due to the incompressibility of the fluid) is
399
conserved only in the full space, but not along subsets with
400
lower dimensionality. For a precise derivation based on the
401
full-dimensional density, seeAppendix H. The factorS(t; u),
402
multiplying the original density, describes the stretching along
403
the material line up to time t experienced near a particle
404
initialized at position u.
405
We can obtain the horizontal density σ(t; u) by pro-
406
jecting the instantaneous density σ(t; u), which is measured
407
along the material line, to the horizontal direction taking
408
into account the kinematics of the problem. In particular, we
409
need to take into account the instantaneous orientation of the
410
material line at the position characterized by u, and also the
411
velocity at the same position:
412
σ(t; u)=σ(t; u)P(t; u), (12) where, according to simple geometry relating the pre- and
413
the post-projection length of an infinitesimal segment of the
414
material line around the position characterized by u,
415
P(t; u)= dfx(t; u)
ds −dfz(t; u) ds
vx[f(t; u), t]
vz[f(t; u), t]
−1. (13) Here, s is the arc length along the image of the line segment at
416
t, and u can be regarded as a function of s. The first term holds
417
alone when there is no horizontal velocity at the given time
418
instant at the position of the given particle, and the second
419
term originates from an additional change in the length, which
420
is due to the presence of horizontal motion. It is worth empha- 421
sizing that the presence of the second term is due to projecting 422
the material line to a given depth, instead of taking the pro- 423
jection at a given time, in agreement with Eq.(7). For a more 424
detailed explanation of the formula, seeAppendix I. This rela- 425
tion is valid for any t and u, so that it also applies to the time 426
instant when a given particle arrives at the accumulation level. 427
In total, there are two independent effects modifying the 428
initial densityσ(t=t0; u): the stretching and the projection, 429
and both of them appear as a factor multiplyingσ(t=t0; u) 430 σ(t; u)=σ(t=t0; u)F(t; u)=σ(t=t0; u)S(t; u)P(t; u), 431
(14) 432
whereF(t; u)is the total factor [the same as in(6), for d=2], 433
andS(t; u)andP(t; u)correspond to the stretching and the 434
projection as defined by Eqs.(11)and(13), respectively. 435
We can simplify the total factor to obtain(6)with(8)as 436
follows. Applying the chain rule for the partial derivatives in 437
(13)yields 438
P(t; u)= df(t; u)
du
dfx(t; u)
du −dfz(t; u) du
vx[f(t; u), t]
vz[f(t; u), t]
−1, (15)
where 439
du ds
= df(t; u)
du
−1 (16) has been used [see Eq.(H11)and the preceding discussion in 440
Appendix H]. Note that, according to(11), 441
df(t; u) du
=S(t; u)−1, (17)
the substitution of which into(15)cancels outS(t; u)in(14): 442
F(t; u)= dfx(t; u)
du −dfz(t; u) du
vx[f(t; u), t]
vz[f(t; u), t]
−1, (18) which is equivalent to(6)–(8)for d=2. 443
The first term in Eq.(18), 444
δx(t; u)=dfx(t; u)
du , (19)
is the parametric derivative, with respect to the position along 445
the initial line segment, of the horizontal component of the 446
current position vector, while the second term, 447
δ˜z(t; u)= −δz(t; u)vx[f(t; u), t]
vz[f(t; u), t] = −dfz(t; u) du
vx[f(t; u), t]
vz[f(t; u), t], 448
(20) 449
is its vertical counterpart, but it is weighted by the ratio of the 450
two velocity components. As in Eq.(13), the former one is 451
due to a “static” change in length (i.e., not influenced by any 452
horizontal motion), and the latter one is the “correction” when 453
horizontal motion is present. The possibility of simplifying 454
Eq. (14)[with Eqs. (11) and Eq. (13)] to Eq. (18) is not a 455
surprise: it is only the ratio between the final and the initial 456
length of an infinitesimal line segment that is relevant, which 457
we have already learnt in Sec.III A. 458
Results for d=3 corresponding to those of this section 459
discussed so far are given in Appendix J, and formulae for 460
d>3 can be constructed similarly. 461
TABLE I. The main quantities relevant for changes in the density.
Notation Name Defining formula
F Total factor (14)
S Stretching factor (11)
P Projection factor (13)
δx Parametric derivative of the horizontal position
(19) δz Parametric derivative of the vertical
position
(20) δ˜z Weighted parametric derivative of
the vertical position
(20)
For d=2, we can summarize our final expression as
462
σ(t; u)=σ(t=t0; u)F(t; u)
463
=σ(t=t0; u)S(t; u)P(t; u)
464
=σ(t=t0; u)δx(t; u)+ ˜δz(t; u)−1, (21)
465
with the particular quantities collected in Table I. Note that
466
a special situation may occur for those trajectories for which
467
|δx+ ˜δz| =0 at the accumulation level. In this case, the final
468
horizontal density is unbounded. The corresponding positions
469
within the accumulation level characterize the so-called (den-
470
sity) caustics,26 and they refer to the maximum levels of
471
inhomogeneity in the accumulated density, so that their identi-
472
fication and dependence on parameters is of great relevance in
473
our work. Of course, the integral of the density (with respect
474
to the final horizontal coordinate x) over such caustics remains
475
finite. In particular, the generic form of a density caustic orig-
476
inating from a standard parabolic fold with its vertex located
477
at xcis∼1/√ x−xc.
478
We can give a more intuitive condition for the positions of
479
the caustics. We first recognize a simplification of(13), which
480
is useful in general, too, and reads as
481
P(t; u)=
vz[f(t; u), t]
n(t; u)·v[f(t; u), t]
, (22) where n(t; u)is the normal vector of the line f at time t at a
482
position that is characterized by u. Equation(22)is true, since
483
n is obtained by rotating the tangent vector df/ds of the line
484
byπ/2
485
n(t; u)=
−dfz(t; u)
ds ,dfx(t; u)
ds . (23)
A remarkable property of (22) is that it remains valid for
486
d=3; seeAppendix Kfor the derivation.
487
The presence of caustics actually originates from the
488
projection factor P alone, and Eq.(22) gives a particularly
489
intuitive interpretation by identifying the positions of the
490
caustics as
491
n(t; u)·v[f(t; u), t]=0. (24) That is, caustics appear in the accumulation plane wherever
492
the local normal vector of the line is perpendicular to the local
493
velocity, or, equivalently, where the local tangent of the line
494
coincides with the direction of the local velocity.
495
IV. NUMERICAL EXAMPLES 496
In this section, we present the basic phenomenology of 497
our setup via numerical examples in a 2D model flow. 498
A. Model flow 499
The equation of motion for the particles, Eq.(4), relies 500
on a fluid flow vfluid(X, t). For clarity, we choose this velocity 501
field to have zero mean integrated over space. Note, however, 502
that as long as the spatial distribution of the particles is inho- 503
mogeneous, the vertical velocity averaged over all particles 504
will be different from−W due to the inhomogeneities of the 505
velocity field.27 506
In order to present relevant phenomena in a clear way, 507
we use a d=2 model flow vfluid(X, t) for our numerical 508
examples: we choose a modified version of the paradigmatic 509
double-shear flow.28 In its classical version, it is a periodic 510
velocity field consisting of a horizontal shear during the first 511
half of the temporal period and of a vertical shear during the 512
other half. We modify this in two aspects: First, we smooth 513
the discontinuous transition between the two orientations by 514
introducing a hyperbolic-tangent-type transition.15 Second, 515
we rotate the shear directions by 45◦, to break the coincidence 516
of the two instantaneous velocity directions with the horizon- 517
tal and vertical axes, which in our sedimentation setup have a 518
very specific role. The resulting velocity field is written as 519
vfluid,x(X, t)= 1
√2[vfluid,ξ(X, t)−vfluid,η(X, t)], (25) 520 521
vfluid,z(X, t)= 1
√2[vfluid,ξ(X, t)+vfluid,η(X, t)], (26) 522
where 523
vfluid,ξ(X, t)=A{1+tanh[γsin(2πt)]}sin[√
2π(z−x)], (27)
524 525
vfluid,η(X, t)=A{1−tanh[γsin(2πt)]}sin[√
2π(z+x)].
(28)
526
γ =20/π controls the temporal sharpness of the shear- 527
direction switching, and it is fixed throughout the paper (as 528
well as the temporal period of the fluid, which is set to 1). A 529
is half of the amplitude of each elementary velocity compo- 530
nent (in what follows: the “amplitude”). By increasing A, we 531
increase the strength of the flow and, as a consequence, also 532
its chaoticity, i.e., the (largest positive) Lyapunov exponent, 533
which is associated with the separation with time of fluid par- 534
ticle trajectories. Note that the velocity field(26)–(28)is also 535
periodic in space, with a period of√
2 in both x and z. For the 536
trajectories, at variance with other implementations of flows 537
related to the double shear, we do not impose any periodic 538
boundary conditions, so that the particles’ positions evolve in 539
the unbounded directions x and z. 540
If we regard the accumulation level as the bottom of the 541
domain of a realistic fluid flow, the velocity field vfluid(X, t) 542 would have to fulfill a no-flux or even a no-slip boundary 543
condition at z= −a, which is not satisfied by(26)–(28). 544
As for the no-flux boundary condition, we do not expect 545
to introduce any qualitative difference compared to the results 546
FIG. 2. The positions of the particles of the initially horizontal material line of unit length, at the indicated time instants.
Dashed lines mark the accumulation levels taken for Figs. 3–5. A=0.06, W=0.6.
obtained in our example flow, since in all our theoretical
547
formulae, the relevant quantity at the accumulation level
548
appears to be not the fluid velocity vfluid, but the particle
549
velocity v, which would not have a vanishing vertical compo-
550
nent. Indeed, we carried out our main analyses in a different
551
flow, namely, a spatially periodic sheared vortex flow with
552
temporal modulation,29,30with accumulation levels fulfilling
553
the no-flux boundary condition, and obtained the very same
554
qualitative results.
555
In principle, a viscous boundary layer with a so-slip
556
boundary condition, or any kind of a separate flow regime
557
at the bottom of the fluid with different characteristics com-
558
pared to the bulk (e.g., length and time scales, magnitude of
559
the velocity), cannot be excluded to leave an important, spe-
560
cific imprint on the qualitative properties of the accumulated
561
particle density. However, with our assumptions and param-
562
eters, as well as in oceanic settings, the time that is spent by
563
a particle in a given layer is mainly determined by the set-
564
tling velocity W , independent of the flow; hence, the effects
565
of any boundary layer or separate flow regime are expected
566
to be negligible if the boundary layer is thin compared to the
567
bulk of the fluid (like in the ocean).
568
Beyond all of the above, in experimental set-ups such as
569
in sediment traps, the accumulation points are not at the bot-
570
tom of the sea, but at some intermediate depth at which no
571
boundary conditions apply at all.
572
B. Illustrative results
573
We now present typical examples for the final density
574
within the accumulation level and show how its form emerges
575
from the reshaping of the material line, which gives rise to
576
the different density-modifying contributions that have been
577
introduced in Sec.III. We always initialize, at t0=0, 10 000 578
particles at z0=0 uniformly in a line segment x∈[0, 1]. 579
(Note that any initial length of the order of unity would suf- 580
fice for our examples.) We follow the particles’ trajectories 581
in the double-shear flow [Eq. (26)] and compute the rele- 582
vant quantities numerically (see TableI). When more than one 583
branch of the material line arrives at the same position (as a 584
result of folding), we additionally calculate the sumF
of 585
the total factorsF corresponding to the individual branches. 586
Furthermore, we compareF
to a normalized histogram h 587
calculated directly from the arrival positions of the individual 588
particles. 589
We start with a parameter setting that does not produce 590
noticeable chaos but leads to regular motion: the portrait 591
of the corresponding stroboscopic map consists of slightly 592
undulating quasi-vertical lines. However, the net horizontal 593
displacement of a trajectory after vertically traversing one 594
spatial period of the flow is not zero generally; it is just very 595
small. 596
Snaphots from the time evolution of the line of particles 597
are shown in Fig.2. At the beginning, both horizontal and rel- 598
ative vertical displacements of neighboring particles remain 599
small, and the line becomes slightly undulated [Fig. 2(a)]. 600
Later on, relative vertical displacements become much larger 601
[see Fig.2(b)]. When they become large enough, it can happen 602
that certain, more slowly falling, parts of the line are folded 603
above the faster parts, as can be observed in Fig.2(c). Such 604
folds, together with the nearly vertical velocity vector, result 605
in caustics after accumulation. 606
Figure 3 considers the parameter setting of Fig. 2 and 607
shows the quantities of Table I for an accumulation level 608
placed at z= −a= −2.7 [marked also in Fig.2(a)]. We can 609
see in Fig. 3(a) that the total factor F computed along the 610
FIG. 3. (a) The total factorF, at the accumulation level a=2.7 [marked by a horizontal dashed line in Fig.2(b)], computed along the individual trajectories according to(18)(in black), and the histogram h (with bin size 0.02) obtained from the positions of the trajectories on the accumulation level (in dark yellow), both as a function of the position along the accumulation level. (b) The total factorF(black) compared to the stretching factorS(orange) and to the projection factorP(blue). (c) The reciprocal of the total factorF(black) compared to the parametric derivative of the horizontal positionδx(green), to the parametric derivative of the vertical positionδz(thin magenta), and to the weighted parametric derivative of the vertical positionδ˜z(thick magenta). See TableIto locate the corresponding formulae. A=0.06, W=0.6.
