Dynamical regimes and clustering of small neutrally buoyant inertial particles in stably stratified turbulence

Dispersion of inertial particles by turbulent flows plays a fundamental role in many geophysical systems, from cloud formation and the dispersion of pollutants in the atmosphere, to the dynamics of plankton in the ocean [1, 2, 3, 4]. In spite of their interest, the dynamics of particles in these systems is poorly understood. We know the equations of motion of small particles (i.e., such that the Reynolds number at the particle scale is much smaller than unity [5]), and in recent years experiments [6] and particle-resolved simulations [7] have provided valuable insights into particles’ dynamics in other regimes. However, particle transport problems in geophysics necessarily require reduced models that simplify the physics, as even in the cases in which a simulation can be done resolving a broad range of scales, an ensemble of runs to get statistical information becomes rapidly unfeasible. This has resulted, e.g., for stratified flows as in the oceans and the atmosphere, in the modeling of inertial particles simply as Lagrangian trancers [8] or using simplified models [9, 10].

Modeling geophysical flows pose other challenges, even without particles. Stably stratified turbulence is anisotropic, and as a result it is fundamentally different from homogeneous isotropic turbulence (HIT) [11, 12, 13]. In these flows stratification reduces the vertical velocity, confining the flow into a quasi-horizontal layered motion, also generating vertically sheared horizontal winds (VSHWs) with strong vertical variability [14]. The stratification results in a restoring force, allowing for the excitation of waves that can coexist with the turbulence. The spectral scaling of stably stratified turbulence is also different than in HIT, with a rich behavior depending on the scale considered, and on many dimensionless parameters. In broad terms, stably stratified turbulence displays an anisotropic subrange with a direct energy cascade between the buoyancy and Ozmidov scales [15, 16]. Studies also indicate that larger scale quasi-horizontal motions can be a continuous source of small scale turbulence as long as the local Reynolds number does not drop below a threshold [17].

Many recent studies have considered stably stratified flows from a Lagrangian perspective (i.e., by considering tracers, or particles with inertia, that are transported by such flows). As an example, the Lagrangian transport of tracers in stably stratified turbulence was studied in [18, 19]. However, the general equations of motion for inertial particles submerged in turbulent flows are not clear. For very small particles the Maxey-Riley approximation provides a set of equations for their dynamics [5]. As particles become larger, Basset-Boussinesq and Faxen corrections become relevant, but for even larger particles such perturbative expansion breaks down. The case of stably stratified turbulence is simpler than HIT in some way: most particles and aerosols are much smaller than the dissipation scale, and thus the Maxey-Riley approximation should hold (except for, e.g., large rain droplets in clouds, or snowflakes). But even in this regime, the derivation of the equations requires certain approximations and depend on the form of the equations for the fluid [20, 21].

A fundamental feature of particles in homogeneous and isotropic turbulent flows is their preferential concentration. Turbulence sometimes separates the particles instead of mixing them. The detailed mechanisms by which turbulence affects particle motions are still unclear. In the homogeneous and isotropic case, and for the average concentration, the main mechanisms behind heavy particles clustering are centrifugal expulsion [5] and the sweep-stick mechanism [22]. Evidence of preferential concentration of heavy particles in laboratory experiments and numerical simulations was reported, e.g., in [23]. Multiscale flow effects may be also relevant [24, 25], and the role of other effects in preferential concentration, such as finite particle radius or the effect of large-scale flows, are still unclear [26, 27, 28]. In the particular case of stratified turbulence it has been shown that inertial particles also cluster for a wide range of parameters [21]. Vertical confinement caused by density stratification produces strong fractal clustering at isopycnic surfaces. Clustering was found to depend on a single parameter, the combination of the Stokes time $\tau_{p}$ of the particles and the Brunt-Väisälä frequency of the flow. In the limit of small $\tau_{p}$ (i.e., small inertia), clustering was found to increase monotonically with $\tau_{p}$ [21].

In this work we present a study considering the Maxey-Riley model for small inertial particles [5], from which we derive an equation for the dynamics of inertial particles in stably stratified flows. We then perform direct numerical simulations of the Boussinesq equations for the fluid, together with the Maxey-Riley equation for one million particles. We derive a simple model for the particles vertical displacement, and compare the model with the simulations to show that particles behave as forced damped oscillators with different regimes depending on the Stokes and Froude numbers. We also characterize the dependence of the stratification-induced vertical confinement of the particles on these two parameters. Finally, we study the formation of clusters using Voronoi tessellation, and show that particles in stably stratified flows tend to accumulate in regions with low vorticity, at least for the range of parameters considered in the present study.

In this work we solve numerically the incompressible Boussinesq equations for the velocity $\bf u$ and for mass density fluctuations $\rho^{\prime}$,

where $p$ is the correction to the hydrostatic pressure, $\nu$ is the kinematic viscosity, ${\bf f}$ is an external mechanical forcing, $N$ is the Brunt-Väisälä frequency (which in this approximation sets the stratification), and $\kappa$ is the diffusivity. In terms of the background density gradient, the Brunt-Väisälä frequency is $N^{2}=-(g/\rho_{0})(d\bar{\rho}/dz)$, with $d\bar{\rho}/dz$ the imposed (linear) background stratification, and $\rho_{0}$ the mean fluid density. We write scaled density fluctuations $\zeta$ in units of velocity by defining $\zeta=g\rho^{\prime}/(\rho_{0}N)$. All quantities are then made dimensionless using a characteristic length $L_{0}$ and a characteristic velocity $U_{0}$ in the domain, resulting in

Inertial particles are modeled using the Maxey-Riley model, but we consider an approximation consistent with those made to obtain the Boussinesq equations, in addition to assuming that the typical length over which the velocity field changes appreciably is much larger than the particle radius $a$. Under the latter hypothesis the Faxén terms are negligible. Under the Boussinesq approximation for a stratified flow, Eqs. (4) and (5) are obtained from the Navier-Stokes equations after neglecting all density fluctuations except for those in the buoyancy force. Thus, for the dynamics of the particles we also consider the density and the mass of the fluid displaced by the particles in terms of their mean values, respectively $\rho_{f}\approx\overline{\rho}_{f}=\rho_{0}$ and $m_{f}\approx\overline{m}_{f}=\rho_{0}V_{p}$ (where $V_{p}$ is the volume of the particles), except in the gravity term. In that term we consider the entire fluid density dependence, $\rho_{f}=\rho_{0}+\partial\bar{\rho}/\partial z(z-z_{0})+\rho^{\prime}$, for a linear background density profile. As the flow is stably stratified, $\partial\bar{\rho}/\partial z<0$. Under these approximations the equation for the particles results in

where ${\bf x}$ is the particle position, ${\bf v}$ is the particle velocity, ${\bf u}({\bf x},t)$ is the fluid velocity at the particle position, $D/Dt$ is the Lagrangian derivative, $d/dt$ is the time derivative following the particle trajectory, and $\rho_{p}$ is the particle mass density (particles are assumed to be spherical). For a fluid at rest, note particles will be at equilibrium (i.e., neutrally buoyant) when $1-\rho_{f}/\rho_{p}=0$, and that there is some freedom on how $\rho_{0}$ and $z_{0}$ are chosen. In particular, without loss of generality we can choose $\rho_{p}=\rho_{0}$, such that particles are neutrally buoyant at $z=z_{0}$ in the absence of density fluctuations.

Multiplying and dividing the buoyancy term in Eq. (6) by $\rho_{0}$ we have,

where the first term inside the brackets in the buoyancy is $-N^{2}$, while the second term is $N\zeta$. Reordering the terms in the equation and using dimensionless units we finally obtain,

where the particle relaxation time is $\tau_{p}=(m_{p}+\bar{m}_{f}/2)/(6\pi a\bar{\rho}_{f}\nu)$. For a spherical particle $\tau_{p}=a^{2}/(3\nu)$, with $\gamma=\bar{m}_{f}/m_{p}=1$. We define the Stokes number as $\textrm{St}=\tau_{p}/\tau_{\eta}$, where $\tau_{\eta}=(\nu/\varepsilon)^{1/2}$ is the Kolmogorov time scale and $\varepsilon$ is the fluid kinetic energy dissipation rate. Note that any other choice for $\gamma=\bar{m}_{f}/m_{p}$ is equivalent to changing the reference value $\rho_{0}$, and results in the particles being neutrally buoyant at a different height (or equivalently, it results in a redefinition of $z_{0}$).

Besides the Stokes number that characterizes the particles, Eqs. (4) and (5) have two controlling parameters for the fluid, the Reynolds and Froude numbers,

where $L$ and $U$ are respectively the characteristic Eulerian integral length and the r.m.s. flow velocity. Using these parameters we can also define the buoyancy Reynolds number

which provides an estimation of how turbulent the flow is at the buoyancy scale $L_{b}=U/N$, and plays an important role characterizing the flow dynamics. For $\textrm{Rb}\gg 1$ strong stratified turbulence can develop, while for $\textrm{Rb}\ll 1$ turbulent motions are strongly damped by viscosity. Geophysical flows typically have large Rb; observations in the ocean thermocline yield $\textrm{Rb}\approx 10^{2}$ to $10^{3}$ [29]. Considering numerical limitations, here we study flows with $\textrm{Rb}>10$. The Ozmidov scale, $L_{oz}=2\pi/k_{oz}$ (with $k_{oz}=\sqrt{N^{3}/\varepsilon}$) also plays an important role in the dynamics, as for scales sufficiently small compared with $L_{oz}$ the flow is expected to recover isotropy. For $\textrm{Rb}>1$, $L_{oz}$ is larger than the Kolmogorov dissipation scale $\eta$, and quasi-isotropic turbulent transport can be expected to take place at small scales.

Setting these parameters and choosing the forcing prescribes the numerical simulations. For the forcing we use Taylor-Green forcing [30], which is a two-velocity components forcing that generates pairs of large-scale counter-rotating eddies perpendicular to the stratification, with a shear layer in between. Its expression is given by

where $f_{0}$ is the amplitude of the forcing, $k_{f}=1/L_{0}$ is the forcing wave number, and $L_{0}$ a unit length. This flow has been used before to study stratified turbulence [17, 19]. In the stratified case it generates a large-scale circulation with VSHWs (i.e., with a non-zero mean horizontal velocity) only in the shear layer between the large-scale Taylor-Green vortices (see [31] for a movie of the development of the non-zero mean horizontal wind in this layer).

The Boussinesq fluid equations, Eqs. (4) and (5), were solved in a triply periodic domain using a parallelized and fully dealiased pseudo-spectral method, and a second-order Runge-Kutta scheme for time integration [32]. The equation for the particles, Eq. (8), was solved using third-order spline interpolation to estimate the forces at the particles positions, and with a second-order Runge-Kutta method for the time evolution [33].

We performed several direct numerical simulations of the Boussinesq equations with different Froude numbers, using a spatial resolution of $N_{x}=N_{y}=768$ and $N_{z}=192$ grid points, in a triple periodic domain of length $L_{x}=L_{y}=2\pi L_{0}$ in the horizontal directions and $L_{z}=H=\pi L_{0}/4$ in the vertical direction. Three different Brun-Väisälä frequences are considered (times are measured in units of a unit turnover time $T_{0}=L_{0}/U_{0}$, with $U_{0}$ a unit velocity, see Table 1 for all the relevant fluid parameters). All simulations have a Prandtl number $\textrm{Pr}=\nu/\kappa=1$, with the kinematic viscosity chosen so that the Kolmogorov scale $\eta=(\nu^{3}/\varepsilon)^{1/4}\approx 0.0072L_{0}$ is well resolved, where the kinetic energy dissipation rate is computed as $\varepsilon=\nu\left<|\bm{\omega}|^{2}\right>$ and $\bm{\omega}=\bm{\nabla}\times{\bf u}$ is the vorticity. This results in $\kappa\eta\approx 1.9$, where $\kappa=N_{x}/(3L_{0})$ is the maximum resolved wave number, corresponding to spatially well resolved simulations [34, 35].

Once the flows in these simulations have reached a turbulent steady state, we randomly distributed particles in a horizontal strip of width $H/5$ centered around $z_{0}=H/2$ (i.e., at the height of the shear layer of the Taylor-Green flow), with initial velocities equal to the fluid velocity at the center of each particle. The dynamics of the particles neglected the last term in Eq. (8) (i.e., the Basset-Boussinesq history term). However, that term was computed a posteriori to estimate its relevance in the dynamics. Particles were one-way coupled, and thus they can be considered as test particles: they do not collide, and their volume fraction is irrelevant for the flow dynamics. Thus, the number of particles should be considered solely as a way to improve the statistics. To each simulation in Table 1 we added four sets with $10^{6}$ particles each, with four different values of $\tau_{p}$ (or equivalently, different Stokes numbers), resulting in a total of 10 data sets of particles with different Fr and St numbers. Table 2 lists the relevant parameters of all the simulations with particles, including the particle Reynolds number $\textrm{Re}_{p}=a|\bf{u}-\bf{v}|/\nu$ in each case. Note that Tables 1 and 2 should be read together, as we can have, e.g., particles with $\textrm{St}=0.3$ in a flow with $N=4/T_{0}$, or the same particles in a flow with $N=8/T_{0}$ or $12/T_{0}$. As a rule, the Basset-Boussinesq force is smaller than the drag force in all simulations except in the cases with $\textrm{St}=3$; for those particles the Basset-Boussinesq history term becomes comparable to the Stokes drag (even though both terms are smaller than buoyancy and added mass forces), and thus studying particles with larger St would require taking this force into account in the dynamics (see, e.g., [20]).

We first study the power spectrum of the particles’ vertical velocity. Figure 1 shows this spectrum for different values of the Froude and Stokes numbers; frequencies are normalized by the Brunt-Väisälä frequency of the carrier flow. A peak is always present at $\omega\approx N$, and for small Fr a second peak at lower frequencies is observed. Its position and amplitude depends on Fr, while its amplitude depends only weakly on St. The peak at $\omega\approx N$ is followed for larger frequencies by a steep spectrum, and decays slowly for smaller frequencies.

The origin of the two peaks in the spectra can be explained by a simple model derived from the equation of motion of the particles. Equation (8) can be rewritten in term of the particles’ vertical position $z$ using that $\dot{z}=v_{z}$ and $\ddot{z}=\dot{v_{z}}$, resulting in

where the Basset-Boussinesq history force was neglected. Rearranging terms in Eq. (12) we arrive at the following expression,

where we assumed that vertical displacements of fluid elements are caused by internal gravity waves, and thus we defined $z_{wav}=z_{0}+\zeta/N$, $\dot{z}_{wav}=u_{z}({\bf x},t)$, and $\ddot{z}_{wav}=Du_{z}/Dt$. Equation (14) is the equation of a driven damped oscillator with system frequency $\sqrt{2/3}N$, damping constant $(2\tau_{p})^{-1}$, and forcing $f_{wav}=\ddot{z}_{wav}+\dot{z}_{wav}/\tau_{p}+2N^{2}z_{wav}/3$. The pulsation of the damped system is $\Omega^{2}=2N^{2}/3-(2\tau_{p})^{-2}$. For particles with small inertia this results in an over-damped system (i.e., $\Omega^{2}<0$) and when perturbed, particles slowly decay to the equilibrium position following fluid elements. Particles with large inertia result instead in weak damping ($\Omega^{2}>0$), and perturbed particles oscillate around the equilibrium as they decay, only weakly following the fluid elements. Indeed, the dependence of the frequency of oscillation with $\tau_{p}$ is in qualitative agreement with the results in Fig. 1(c) and (d); note that as St increases, the main peak of the spectrum moves from $\omega\approx N$ to lower frequencies (as a reference, for $\textrm{St}=3$ Eq. (14) yields a frequency $\Omega\approx 0.82N$).

Equation (14) can be integrated numerically if $f_{wav}$ is prescribed. As we do not know the precise evolution of $u_{z}$ as seen by each particle, we assume $u_{z}$ is a random colored process. The spectrum of the fluid vertical velocity in many stably stratified flows is compatible with the Garrett-Munk spectrum, as observed in oceanic observations [36] and in numerical simulations [37]. This is a flat power spectrum for frequencies $\omega\leq N$, resulting from the a superposition of internal gravity waves, followed by a power law decay for $\omega>N$. Thus, we consider a random superposition of oscillators of the form

where $\phi_{\omega}$ are random phases (note that, as we are interested only in vertical motions, the dependence of traveling waves on $x$ and $y$ can be ignored or absorbed into the random phases), and $u_{0}$ is an amplitude chosen so that $\dot{z}_{wav}$ has the same r.m.s. value as that of $u_{z}$ in the numerical simulations. The power spectrum of $\dot{z}_{wav}$ that results is compatible with oceanic observations of the Garret-Munk spectrum [2, 36]. In other words, this process results in $\dot{z}_{wav}$ being a random variable compatible with that spectrum.

Figure 2 shows the power spectrum obtained after integrating Eq. (14) using this random process as the forcing, for different values of the Brunt-Väisälä frequency and the Stokes number. Spectra are qualitatively similar to those shown in Fig. 1. The spectra display two peaks, the main one close to $\omega\approx N$. At fixed Fr, increasing St results in a broadening of this peak towards smaller frequencies. The second peak at lower frequencies has increasing amplitude with decreasing Fr, and appears at similar frequencies as those in Fig. 1. It is natural to ask whether these peaks are caused by the forcing or by the damped oscillations of the particles. Changing the forcing while still maintaining the Garret-Munk spectrum for $u_{z}$ (e.g., setting $f_{wav}=\ddot{z}_{wav}$) yields the same qualitative results, which indicates the peaks in the spectra are partially associated to the damped dynamics of the particles.

Equation (14) can be also rewritten in terms of the vertical slip velocity, $v_{slip}=u_{z}-v_{z}$. Taking $y=z-z_{wav}$, Eq. (14) results in the homogeneous damped harmonic oscillator equation. As before, the resulting oscillation frequency is $\Omega^{2}=2N^{2}/3-(2\tau_{p})^{-2}$, with exponential decay rate $(2\tau_{p})^{-1}$. Noting that $\dot{y}=v_{slip}$, we can expect the vertical slip velocity of the particles to display overdamped or underdamped oscillations depending on the sign of $\Omega^{2}$. Figure 3 shows $v_{slip}$ for particles in the numerical simulations with $\textrm{St}=3$ ($\tau_{p}=0.235T_{0}$) in an stratified fluid with $N=12/T_{0}$ and $8/T_{0}$. Both cases have $\Omega^{2}>0$, and dynamics reminiscent of underdamped oscillations can be identified in the time series. The power spectrum of $v_{slip}$ for multiple simulations, also shown in Fig. 3, shows peaks at the expected value of $\Omega$ in these two simulations, and no peaks in the other simulations with $\Omega^{2}<0$. Thus, the dynamics of the individual particles is compatible with randomly forced damped oscillators, with both $\tau_{p}$ and $N$ controlling the particles’ dynamical regime.

Stratification limits vertical motions of the particles, strongly impairing vertical dispersion, and resulting in saturation of the mean squared vertical displacements of the particles with time. Linear models predict this saturation to take place after $t\approx 2\pi/N$, as particle displacements get constrained vertically by stratification, resulting in oscillatory motions around the neutrally buoyant equilibrium [38]. This was confirmed in numerical simulations with moderate Rb [11]. Later, studies of vertical dispersion of tracers in stably stratified flows [18, 37], and of small neutrally-buoyant particles with small St [20], explicitly confirmed the saturation of the mean squared vertical dispersion. For neutrally-buoyant inertial particles the saturation was found to be faster and stronger than for Lagrangian tracers [20].

Figure 4 shows the particles’ mean squared vertical displacement in the simulations, $\left<\delta z^{2}(t)\right>=\left<[z_{i}(t)-z_{i}(0)]^{2}\right>_{i}$ (where the subindex $i$ indicates the average is computed over all particle labels), for different values of Fr and St. Time is normalized by $2\pi/N$ and $\left<\delta z^{2}(t)\right>$ is normalized by $(U_{z}/N)^{2}$, where $U_{z}$ is the Eulerian r.m.s. fluid vertical velocity in the turbulent steady state. With this normalization curves collapse from $t=0$ to $t\lesssim 2\pi/N$, in a time interval with ballistic behavior. The end of this regime at a time proportional to the wave period $2\pi/N$, instead of the Lagrangian eddy turnover time, indicates that the rapid early vertical displacements are caused by the inertial particles following the inertial waves. Note also that there is more overshooting in the vertical displacements (i.e., $\langle\delta z^{2}\rangle$ reaches larger values in its maximum at the end of this ballistic stage) as St (and particle inertia) increases. After this maximum, inertial particles oscillate around their equilibrium position, displaying a plateau in the mean squared vertical displacements, as also reported in [20]. The amplitude of the plateau is weakly dependent on St, and depends strongly on Fr (see [31] for a movie showing the vertical displacements of the three different types of particles when $N=4/T_{0}$, illustrating the confinement in vertical layers). This is different from the case of tracers in stratified flows, which for sufficiently large Rb display some slow vertical dispersion at late times caused by turbulent eddies or by diffusion [37].

The confinement of particles around a layer can be also characterized using the probability density function (PDF) of finding a particle at a given height, either in terms of $z$, or of the density at each particle position (i.e., of how far the particle is from the equilibrium isopycnal). Figure 5 shows the PDF of $z$ for different Fr and St, centered by the mean value and normalized by the dispersion. Remarkably the PDFs have asymmetric tails, the more stronger as Fr is increased. This can be associated to the occurrence of extreme vertical drafts in stably stratified flows for values of Fr in the range $\approx 0.05$ to $0.30$ [39], which can cause more frequent and larger vertical wanderings of the particles. Figure 6 shows the same PDFs but in terms of the rescaled density fluctuations $\zeta$ at the particles positions, also centered by the mean value and normalized by the dispersion. These PDFs are closer to Gaussian and less sensitive to St, but still display asymmetric tails for $\textrm{Fr}=0.09$.

From Fig. 4 it seems apparent that particles are confined in a narrower layer as Fr decreases, but this is not evident from Figs. 5 and 6 as the PDFs in those figures are normalized by their standard deviations. Figure 7 show the standard deviations in $z$ and $\zeta$ of the particles, $\sigma_{z}$ and $\sigma_{\zeta}$ respectively, as a function of $\textrm{Fr}^{-1}$ for all St considered. Note that both deviations (which can be considered as a measure of the height of the confinement layer) are a fraction of $L_{b}$, and decrease with decreasing Fr. The behavior of $\sigma_{z}$ depends also on St for weak stratification. Sozza et al. [40] observed that for large values of Fr, $\sigma_{z}$ was larger for larger $\tau_{p}$ (i.e., larger St), the opposite behavior of what is found here. The effect in [40] resulted from particles with more inertia being suspended from the equilibrium position for longer. Here, the mean winds in the shear layer of the Taylor-Green flow result in a different effect. Particles with more inertia are less affected by rapid vertical motions, following instead the slower horizontal motions and averaging over the vertical fluctuations as they move (see the movie in [31]). This is evident in Fig. 1, where it is observed that for the case with $\textrm{St}=0.3$ the main peak of the power spectrum of the particles’ vertical velocity is above unity, while for the case with $\textrm{St}=3$ the peak is below unity. This confirms that particles with lower inertia are more affected by the fluid vertical displacements, while particles with larger inertia are less affected by them. Thus, the behavior of $\sigma_{z}$ with St for larger Fr is not universal, and probably dependent on the flow. However, the situation is different for $\sigma_{\zeta}$, as shown in Fig. 4(b): $\sigma_{z}/(L_{b}N)=\sigma_{z}/U$ seems to depend linearly on $\textrm{Fr}^{-1}$, and is independent of St, at least in the range of parameters considered. Note this amounts to the dispersion of the particles around the isopycnal decreasing linearly with increasing Brunt-Väisälä frequency.

The vertical confinement of particles have consequences for cluster formation. To quantify it we use Voronoï tessellation. Tessellations have been shown to be useful to characterize preferential concentration of particles, see, e.g., [41, 42, 43, 6, 44, 45], with the standard deviation of the Voronoï cell volumes or areas being associated to the amount of clustering [42, 43, 23]. For heavy particles in stratified turbulence, clustering has also been studied using radial distribution functions [46], which give the ratio of the number of particle pairs found at a given separation to the expected number of pairs if particles are uniformly distributed. A Voronoï tessellation assigns a cell to each particle, so that each point in that cell is closer to that particle than to any other particle. Large tessellation cells correspond to voids (i.e., regions with far apart particles), while small cells correspond to clustered particles. While in the case of homogeneous and isotropic turbulence both three-dimensional (3D) and two-dimensional (2D) tesselations have been used, here we restrict ourselves to 2D tessellation as most of the particles remain in the thin layers discussed in the previous section. To that end, we project all particles into a plane, and consider only their $x$ and $y$ coordinates.

Figure 8 shows the PDFs of the normalized areas of the Voronoï cells, $\altmathcal{A}=A/\langle A\rangle$, where $A$ is the area of each cell. The figure also shows as a reference the PDF of a random Poisson process (RPP), which corresponds to particles randomly distributed in space [47, 48]. The first crossing from the left between the PDFs and the RPP is often used to define clusters: an excess of smaller cells are an indication of a spatial accumulation of particles in certain regions of the flow. Note that particles in the flow with $\textrm{Fr}=0.2$ are closer to the RPP. This is to be expected, as neutrally buoyant small particles do not cluster in the limit of homogeneous and isotropic turbulence (i.e., for large enough Fr) [49]. For fixed Fr, clustering increases with St. But more importantly, clustering increases rapidly as Fr is decreased. The strongest clustering is obtained for intermediate stratification at $\textrm{Fr}=0.12$ and $\textrm{St}=1$. This indicates that although the increase in stratification is favorable for cluster formation, its effect is not monotonous with Fr.

The clusters in this case seem to form in regions of the flow with low vertical vorticity, thus resulting from centrifugal vortex expulsion [5]. To illustrate this Fig. 9 shows the Voronoï areas of a random subset of $10^{4}$ particles with $\textrm{St}=0.3$, in the three simulations with $\textrm{Fr}=0.20$, $0.12$, and $0.09$. Red regions in panels (a), (b), and (c) correspond to cells larger than the average (voids), and blue areas to cells smaller than the average (clusters). A movie with the time evolution of the particles in the case with $\textrm{Fr}=0.12$ can be seen in [31]. Panels (d), (e), and (f) show the squared vertical vorticity $\omega_{z}^{2}$ averaged in each Voronoï cell, and normalized by its mean value (as a reference, the bottom panels show the same vorticity at full resolution, i.e., not coarse-grained). Note there is some correlation between these panels: regions of low vorticity seem to correspond to smaller Voronoï areas, specially for small Fr. A similar correlation between clusters and low vorticity regions was reported before for the case of heavy particles in stratified turbulence in [46].

Figure 10 further confirms this correlation by showing joint PDFs of $\mathcal{A}$ vs ${\omega_{z}}^{2}$ for all simulations and particles. In the panels the vertical dashed lines indicate, from left to right, the first and second crossings of the PDFs of $\mathcal{A}$ with the RPP (i.e., the values of $\mathcal{A}$ below and above which cells correspond respectively to clustered particles, and to voids). Particles in voids tend to be in regions of larger vorticity, and the correlation is more clear as Fr decreases. As a reference, we indicate different slopes with straight lines. Note that for strong stratification ($\textrm{Fr}=0.12$), the shape of the PDFs becomes almost insensitive to the value of St. Overall, a correlation between large Voronoï areas and large vorticity appears independently of the Stokes number in the strongly stratified cases.

We presented a numerical study of the transport and spatial accumulation of light neutrally-buoyant inertial particles in stably stratified turbulent flows, using the Maxey-Riley equation for small particles. We showed that in the stratified case, the equation can be written as the equation of a driven damped oscillator, with two regimes controlled by the inverse squared particle response time, $\tau_{p}^{-2}$, and the flow Brunt-Väisälä frequency, $N$. When the former is larger particles are overdamped, while when the latter is larger particles are underdamped. This results in the appearence of two peaks in the power spectrum of the particles’ vertical velocity, the main peak with frequency $\Omega\approx[2N^{2}/3-(2\tau_{p})^{-2}]^{-1/2}$.

As observed in previous studies of light and heavy particles in stably stratified turbulence [46, 20], the vertical dispersion of particles is strongly confined in layers. The width of this layer depends on the Stokes and Froude numbers. However, when studied in terms of density isopycnals, the width becomes independent of the particles’ Stokes number (at least, in the range of parameters considered in this study), and varies with the Brunt-Väisälä frequency.

This vertical confinement also has a strong impact in the clustering of particles, and in the physical mechanism behind cluster formation. We showed that a two-dimensional Voronoï tesselation can be used to study clusters; previous studies using other methods for vertically confined particles can be seen in [46, 50]. Our analysis indicates that in sufficiently stratified flows the formation of clusters is governed by centrifugal vortex expulsion, independently of whether the Stokes number is smaller or larger than unity. Moreover, clustering is strongly enhanced as stratification is increased, and this enhancement takes place also when only considering the two-dimensional positions of the particles, i.e., independently of the particles’ vertical confinement. This result can be important to compute particle collisions and particle-turbulence interactions in atmospheric problems [51], and in oceanic flows where patches of phytoplankton and of nutrients are commonly observed [52, 53, 54]. The result is also reminiscent of observations of clustering in floaters in free surface flows, where the large scale flow circulation can play an important role in particle accumulation [55].