1. Introduction
Settling of solid particles in a fluid is a fundamental process relevant to a wide range of natural and industrial systems, including sediment transport, pollutant transport in aquatic environments and particle separation in chemical engineering.
The dynamics of individual particle settling in quiescent fluids have been investigated extensively, resulting in a detailed classification of motion regimes based on key non-dimensional parameters. Jenny, Dusek & Bouchet (Reference Jenny, Dusek and Bouchet2004) described particle trajectory regimes using the particle-to-fluid density ratio
$\rho _p/\rho _{\!f}$
and Galileo number,
$\textit{Ga} = ( \lvert \rho _p/\rho _{\!f} - 1 \rvert g d^3 )^{1/2}/\nu$
, which quantifies the relative influence of gravitational to viscous forces. Here,
$\rho _p$
and
$\rho _{\!f}$
are the particle and fluid densities,
$g$
is gravitational acceleration,
$d$
is the particle diameter, and
$\nu$
is the kinematic viscosity of the fluid. For spherical particles, large Galileo numbers (
$\textit{Ga} \gtrsim 200$
) are associated with the onset of wake instability, as the steady axisymmetric wake becomes unstable and gives rise to vortex shedding, generating unsteady lift forces that result in chaotic and non-vertical settling trajectories. The transition to this regime depends sensitively on the particle-to-fluid density ratio and has been well established through a combination of laboratory experiments (Veldhuis & Biesheuvel Reference Veldhuis and Biesheuvel2007; Horowitz & Williamson Reference Horowitz and Williamson2008, Reference Horowitz and Williamson2010) and direct numerical simulations (Zhou & Dušek Reference Zhou and Dušek2015).
While the behaviour of individual particles is well understood, the dynamics becomes significantly more complex when inter-particle interactions are introduced. In particular, the settling of particle pairs can give rise to phenomena such as drafting, kissing and tumbling (DKT), in which the motion of one particle directly alters the trajectory of its neighbour (Fortes, Joseph & Lundgren Reference Fortes, Joseph and Lundgren1987; Feng, Hu & Joseph Reference Feng, Hu and Joseph1994). Building on these foundational observations, subsequent studies have examined the dynamics of two-particle systems across several configurations. Li et al. (Reference Li, Liu, Zhao, Yin and Lu2022) demonstrated numerically that attractive interactions between particles are highly sensitive to both their initial separation and relative orientation. In side-by-side arrangements, experimental and numerical investigations by Liu et al. (Reference Liu, Zhang, Xiao, Wang, Yuan and Tang2021) showed that asymmetry in the induced vortex structures generated repulsive forces, with interaction strength governed by the particle Reynolds number and initial spacing. Ghosh & Kumar (Reference Ghosh and Kumar2020) further explored the influence of differing particle size and density, revealing that both parameters strongly modulate the onset and evolution of DKT.
The settling behaviour of large particle populations is strongly influenced by collective hydrodynamic interactions, which can lead to either enhanced or hindered settling (Penlou et al. Reference Penlou, Roche, Manga and van den Wildenberg2023). When particles form clusters, the ensemble may descend cooperatively, with wake overlap and entrainment generating low-pressure regions that accelerate the group relative to an isolated grain. In dilute regimes (
$\phi \lesssim 0.5\,\%$
) at sufficiently large Galileo numbers (
$\textit{Ga} \gtrsim 150$
), settling velocities have been reported to exceed single-particle values by approximately 20 %–25 %, associated with the formation of vertically aligned particle trains and wake-mediated entrainment (Kajishima Reference Kajishima2004; Huisman et al. Reference Huisman, Barois, Bourgoin, Chouippe, Doychev, Huck, Morales, Uhlmann and Volk2016). In contrast, at higher concentrations where particles are more uniformly distributed, wake interference and direct hydrodynamic hindrance increase the effective drag, reducing mean settling velocities below those of an isolated particle (Baas et al. Reference Baas, Baker, Buffon, Strachan, Bostock, Hodgson, Eggenhuisen and Spychala2022; Penlou et al. Reference Penlou, Roche, Manga and van den Wildenberg2023).
A widely used description of hindered settling is the empirical Richardson–Zaki correlation (Richardson & Zaki Reference Richardson and Zaki1954), which relates mean settling velocity to the solid volume concentration as
$w = w_0(1-\phi )^n$
, where
$w_0$
is the terminal velocity of a single particle and
$n$
is an exponent that depends on the particle Reynolds number. While this relation is widely used, several numerical and experimental studies have shown that it tends to overpredict settling velocities in the dilute inertial regime, where wake interactions lead to anisotropic particle distributions and enhanced fluctuations (Yin & Koch Reference Yin and Koch2007; Zaidi, Tsuji & Tanaka Reference Zaidi, Tsuji and Tanaka2015). Shajahan & Breugem (Reference Shajahan and Breugem2023) performed particle-resolved DNS of monodisperse suspensions for
$\phi = 0.005$
–0.3 at
$\textit{Ga} = 144, 178, 210$
, covering regimes where single particles display axisymmetric, oblique and oscillatory wakes. They identified three regimes: dilute (
$\phi \lesssim 0.02$
) with vertical clustering, moderate (
$0.02 \lesssim \phi \lesssim 0.1$
) dominated by DKT-driven horizontal pairs and dense (
$\phi \gtrsim 0.1$
) with nearly random hard-sphere configurations. The mean settling velocity followed
$w/w_0 \propto (1-\phi )^n$
with
$n \approx 3.1$
–3.5 in dilute-to-moderate regimes, approaching
$n \approx 2.7$
only for
$\phi \gt 0.1$
, consistent with diminished wake-trapping and DKT interactions at higher
$\phi$
. They also reported kinematic waves whose speed agreed with Kynch theory (Kynch Reference Kynch1952), using the Richardson–Zaki relation for the mean sedimentation velocity, and whose frequency coincided with peaks in the particle velocity spectra, linking concentration fluctuations to velocity modulation. Together, these results attribute deviations from classical hindered-settling laws at low-to-moderate
$\phi$
to the competition between vertical clustering and horizontal pair formation, which shape microstructure anisotropy and effective drag. Kang et al. (Reference Kang, Hong, Cheng, Best and Chamorro2023, Reference Kang, Best and Chamorro2024) experimentally studied vertically aligned arrays of settling particles in homogeneous and density-stratified media, highlighting the role of the Galileo number and the presence of fluid density interfaces in modulating inter-particle interactions and overall descent behaviour. Liu et al. (Reference Liu, Xiao, Liang, Zhang, Wang, Liu, Zhang and Zhou2024) examined settling of multiple particles initially arranged side-by-side and observed that spatial trends in settling velocity were strongly dependent on initial inter-particle spacing.
Although previous studies have provided valuable insights into settling dynamics, they have predominantly considered particles of uniform density. However, many natural and industrial suspensions contain particles with varying densities, whose interactions can significantly alter transport behaviour. Although some studies have explored particle density heterogeneity, they have largely been limited to pairwise interactions or small particle ensembles. As a result, the collective dynamics of mixed-density particle settling at larger scales, particularly in configurations where spatial arrangement influences hydrodynamic coupling, remain insufficiently understood.
In many natural and industrial systems, particles with differing densities coexist and interact during sedimentation, leading to complex collective dynamics that influence transport, deposition and mixing. For example, in salt marshes, the interplay between mineral sedimentation and organic matter accumulation governs net carbon sequestration, a process sensitive to climatic drivers such as warming and sea-level rise (Kirwan & Mudd Reference Kirwan and Mudd2012). Understanding how differences in particle density modulate inter-particle interactions and settling behaviour is essential for improving predictions of material flux in heterogeneous suspensions. A critical case involves microplastics in aquatic environments, where particle mixtures span a wide range of densities from 0.88 to 2.2 g
$\textrm {cm}^{-3}$
with grain sizes ranging from 1
$\mu$
m to 5 mm (Baensch-Baltruschat et al. Reference Baensch-Baltruschat, Kocher, Stock and Reifferscheid2020; Goßmann et al. Reference Goßmann, Mattsson, Hassellöv, Crazzolara, Held, Robinson, Wurl and Scholz-Böttcher2023; Poulain-Zarcos et al. Reference Poulain-Zarcos, Pujara, Verhille and Mercier2024; Stride et al. Reference Stride, Abolfathi, Bending and Pearson2024), and can form aggregates of heterogeneous density that may incorporate organic particles. Riverine pollutant transport has been the subject of considerable scientific attention, with particle settling velocity being one of the most sensitive and influential parameters in predictive transport modelling (Unice et al. Reference Unice, Weeber, Abramson, Reid, Van Gils, Markus, Vethaak and Panko2019a
,
Reference Unice, Weeber, Abramson, Reid, Van Gils, Markus, Vethaak and Pankob
; Waldschläger et al. Reference Waldschläger, Lechthaler, Stauch and Schüttrumpf2020; Li et al. Reference Li, Bai, Chen, Li, Hu, Krause and Schneidewind2026). However, the settling velocities of pollutants are typically measured in experiments using particles of a single density (Goral et al. Reference Goral, Guler, Larsen, Carstensen, Christensen, Kerpen, Schlurmann and Fuhrman2023; Dittmar et al. Reference Dittmar, Ruhl, Altmann and Jekel2024). Our current understanding of riverine pollutant transport behaviour is, therefore, incomplete given the complex settling behaviour expected when particle mixtures of varying densities are present. Tyre wear particles (TWP) and tyre and road wear particles (TRWP), formed after combination with road particles, have recently been recognised as a significant source of pollution into the environment (Baensch-Baltruschat et al. Reference Baensch-Baltruschat, Kocher, Stock and Reifferscheid2020). These particles have wide density ranges (TWPs = 1.13–1.16 g
$\textrm {cm}^{-3}$
, TRWPs = 1.5–2.2 g
$\textrm {cm}^{-3}$
) (Rhodes, Ren & Mays Reference Rhodes, Ren and Mays2012; Baensch-Baltruschat et al. Reference Baensch-Baltruschat, Kocher, Stock and Reifferscheid2020). Understanding the settling behaviour of mixtures of different density particles is thus essential in predicting the fate of TRWPs in the environment. Despite substantial progress, the collective settling dynamics of particles with differing densities remain poorly understood. Clarifying how cross-density interactions modify wake coupling, entrainment and relative dispersion is essential for advancing the physical understanding of multi-particle hydrodynamic interactions at finite Reynolds numbers (Nie, Lin & Gao Reference Nie, Lin and Gao2017).
To begin to address this gap in knowledge, we investigate experimentally the multi-particle settling dynamics of a layer of grains, by combining particles of two densities under various initial arrangements. Two types of spherical particles with density ratios
$\rho _p/\rho _{\!f} = 1.14$
and
$1.28$
are used in this study, which are within the range of densities relevant to TWPs and TWRPs. The particle arrangements studied herein enable the identification of configuration-sensitive, hydrodynamic interactions that emerge at intermediate-to-high Galileo numbers (
$\textit{Ga}$
), revealing flow-induced segregation and wake-mediated dispersion not accessible in confined or small-scale systems. The present study aims to understand mesoscale organisation in multiphase flows and provide new insights into how particle heterogeneity governs bulk transport processes. Section 2 describes the experimental set-up, including the particle release mechanism, measurement techniques for particle tracking and flow characterisation, and the configurations tested. Section 3 presents the principal results and § 4 discusses the underlying hydrodynamic mechanisms. Finally, § 5 summarises the main findings and outlines directions for future research on multiphase collective settling.
2. Experimental set-up
The experiments were conducted in a transparent acrylic tank filled with quiescent water, with dimensions of 600 mm in height and a square cross-section of 300 mm
$\times$
300 mm. Arrays of 100 spherical particles of diameter
$d = 4$
mm were arranged in a horizontal 10
$\times$
10 Cartesian grid near the free surface and released simultaneously. The particles in the outermost layer were located approximately 27.5 particle diameters away from the wall. At this separation, the wall-correction factor is close to unity and wall effects can therefore be assumed to be minimal (Miyamura, Iwasaki & Ishii Reference Miyamura, Iwasaki and Ishii1981). The particles were composed of nylon and acetate with density ratios
$\rho _p/\rho _{\!f} = 1.14$
and
$1.28$
, respectively. The corresponding Galileo numbers were
$\textit{Ga} = 297$
for the lighter particles and
$\textit{Ga} = 419$
for the heavier grains, placing them both within the chaotic regime of sphere settling (Jenny et al. Reference Jenny, Dusek and Bouchet2004).
(a) Experimental set-up showing particle tracking velocimetry (PTV)/PIV systems, release chamber and tank. (b) Initial particle array configurations for the five cases; blue and red dots denote
$\rho _p/\rho _{\!f} = 1.14$
and
$1.28$
, respectively. (c) Definition of particle layers
$L_i$
, with
$L_1$
outermost and
$L_5$
innermost.
Five distinct particle array configurations were tested to investigate the effects of density heterogeneity and initial spatial arrangement on collective settling behaviour. Cases 1 and 2 serve as homogeneous reference cases, composed entirely of lighter and heavier particles, respectively. Case 3 introduces lateral heterogeneity by dividing the array into two equal regions, one with lighter particles and the other with heavier particles. Case 4 consists of a 6
$\times$
6 central core of lighter particles surrounded by heavier particles, while Case 5 reverses this arrangement, placing heavier particles at the centre and lighter ones on the periphery. Figure 1 illustrates the experimental set-up and all five particle arrangements. For reference, we define concentric particle layers numbered 1–5, progressing from the outermost to the innermost ring of the array (figure 1
c).
A custom-designed release mechanism was developed to ensure the simultaneous deployment of all 100 particles with minimal disturbance to the fluid surface and negligible initial particle rotation. The system consisted of a rigid rectangular chamber with a
$10 \times 10$
array of precisely machined holes, each slightly smaller than the particle diameter, to ensure stable positioning. A vacuum pump was connected to the chamber to generate a negative pressure that held the particles in place. During set-up, particles were manually positioned over each hole, sealing them against the suction flow. Once all particles were secured, the chamber was carefully positioned 10 mm above the free surface of the water. The release was initiated by turning off the vacuum pump, eliminating the pressure differential and allowing all particles to fall freely into the fluid synchronously. An example of this particle release is illustrated in supplementary movie 1 available at https://doi.org/10.1017/jfm.2026.11570.
The three-dimensional (3-D) trajectories of the settling particles were captured using a multi-view imaging system composed of three 2.3 MP synchronised high-speed cameras, each operating at a sampling frequency of 120 Hz. The cameras were arranged with oblique viewing angles to maximise the reconstruction volume and minimise image overlap, and were focused on a region containing the interrogation depth of particle descent. Each experimental case was repeated 10 times to ensure statistical representativeness and account for inter-run variability. The three-dimensional particle positions were reconstructed through triangulation based on the polynomial Tsai camera model. Calibration was performed by imaging a planar calibration target traversed through multiple known axial positions using a precision translation stage, similar to the procedure described by Kreizer & Liberzon (Reference Kreizer and Liberzon2011). For the 3-D particle reconstruction, only regions with a reconstruction error below 1.5 pixels were considered, corresponding to approximately 0.1 % of the investigation domain. The reconstructed trajectories were post-processed using the trackpy linking algorithm (Allan et al. Reference Allan, Caswell, Keim, van der, Casper and Verweij2021), which links particle positions across frames based on predictive kinematic filtering and proximity constraints. The velocity, estimated from a short history window of up to the most recent eight frames and including directional information, is used to predict particle positions and constrain the search region in the subsequent frame.
The flow field generated by the settling particles was measured using planar particle image velocimetry (PIV) in a vertical plane and stereoscopic planar PIV in a horizontal plane. Vertical PIV planes intersected the centre of the initial particle array and spanned a vertical extent of
$35d$
(142 mm), beginning at a depth of
$2d$
(8 mm) below the initially quiescent free surface and extending
$37.5d$
(150 mm) in the horizontal direction. The horizontal PIV plane was positioned at a fixed depth of
$z/d = 25$
. Silver-coated particles with a diameter of 44
$\unicode{x03BC}$
m and density
$\rho _t=1.1$
g ml-1 acted as tracers and a 532 nm laser was used to illuminate the field of view. Flow fields were obtained at 200 Hz using a 2.3 MP high-speed camera mounted perpendicular to the laser sheets for planar PIV, and two angled cameras for stereoscopic PIV. Velocity fields were obtained through standard two-frame cross-correlation of image pairs using interrogation windows of
$16 \times 16$
pixels with 50 % overlap, resulting in a vector spacing of 1.6 mm (
$0.4d$
).
3. Results
This section presents the main experimental results, while detailed discussion and further perspectives are presented in subsequent sections. Representative three-dimensional particle trajectories from each of the five configurations (figure 2) form the basis for the subsequent analysis of particle motion and kinematics. These trajectories are coloured by particle type, with blue corresponding to particles of density ratio
$\rho _p/\rho _{\!f} = 1.14$
and red to
$\rho _p/\rho _{\!f} = 1.28$
. The inset panel displays four particle trajectories coloured by the local vertical velocity, illustrating their spatial and temporal variability. The tracking of each particle allows for detailed quantification of grain dispersion and interparticle dynamics. All results herein are based on ensembles compiled from 10 repeated experiments per configuration, supporting statistical reliability. While the figure provides a qualitative view, the observed spatial patterns already suggest the emergence of configuration-dependent interparticle interactions and collective settling behaviours.
Example three-dimensional particle trajectories for the five configurations. Each trajectory is coloured by particle type: blue for particles with
$\rho _p/\rho _{\!f} = 1.14$
and red for
$\rho _p/\rho _{\!f} = 1.28$
. Cases 1 and 2 correspond to homogeneous arrays of light and heavy particles. The top-right inset shows representative particle velocities for four grains in Case 4, normalised by the terminal velocity of a solitary grain.
Ensemble-averaged bulk particle velocities in the lateral
$v_L$
(left column) and vertical
$w$
(right column) directions. Solid lines represent particles in the outer two layers and dashed lines correspond to the inner three layers. Shaded regions indicate the standard deviation across ten realisations. Velocities are normalised by the terminal velocity of a single lighter particle (
$\rho _p/\rho _{\!f} = 1.14$
) for Case 1 and by that of the heavier particle (
$\rho _p/\rho _{\!f} = 1.28$
) for Cases 2–5. Blue secondary ordinate axes indicate normalisation with respect to the terminal velocity of the lighter particles
$w_0^L$
to facilitate comparison across density configurations.
To characterise and quantify the settling dynamics, we first examine the averaged particle group velocities in the vertical and lateral directions for each scenario. Particular attention was given to differences in settling behaviour between particles near the periphery of the array and those situated in the central region, following the layer nomenclature defined in figure 1(c). For each experimental realisation, a reference time
$t_0 = 0$
was defined as the instant when the average vertical velocity of the particle ensemble (or of the lighter particles, in mixed-density cases) reached a minimum, marking the end of the initial release phase. Velocities were ensemble-averaged across the ten runs for each case and figure 3 presents the resulting behaviours with shaded bands indicating one standard deviation across experimental repetitions. Quantities were non-dimensionalised using the terminal velocity of an isolated particle, estimated from the semi-empirical expression
$w_0 = [(729 + 3 \textit{Ga}^2)^{1/2} - 27] \nu /d$
(Goossens Reference Goossens2020). The reference velocity for normalisation was
$w_0 = 122$
mm s
$^{-1}$
in Case 1, corresponding to the terminal velocity of a light particle, and
$w_0 = 175$
mm s
$^{-1}$
for Cases 2–5, associated with the denser particle.
Vertical velocity profiles (figure 3) reveal a monotonic sub-linear velocity increase phase until reaching a terminal value, with the group mean value being less than that of the terminal velocity of a single grain, illustrating that all experiments are within the broad regime of hindered settling. Although the particles are initially arranged in a monolayer, the configuration can be interpreted within a Richardson–Zaki hindered settling framework by introducing an effective initial solid volume fraction. For an initially planar release, this quantity is not uniquely defined and is dependent on the effective thickness
$H$
assigned to the monolayer. Here, we adopt
$H = d$
, corresponding to one particle diameter, and define the reference volume as
$V_0 = 20d \times 20d \times d$
. This convention converts the initial areal packing into an equivalent initial volumetric concentration
$\phi _{V_0} = 0.13$
. For this effective volume fraction, the bulk terminal velocity is approximately 80 % of the corresponding single particle settling velocity, resulting in an effective Richardson–Zaki exponent
$n = 1.6$
. This reflects the reduction in mean settling velocity arising from collective wake interactions. It is worthy of note that alternative choices of
$H$
will yield different values of
$n$
. This sensitivity indicates that
$n$
in the present context should not be interpreted as a classical suspension exponent, but rather as an empirical measure of collective drag modification for a finite, initially planar particle ensemble.
In homogeneous systems, the outer-layer (
$L_1+L_2$
, solid lines; figure 3) particles settle faster than inner-layer counterparts, thus forming a concave ‘parachute’ shaped structure that is persistent in time due to the same terminal velocity regardless of the layer set. In mixed-density systems (Cases 3–5), vertical velocities of individual particle types remained consistent with those observed in their corresponding homogeneous cases, reflecting the dominant role of particle density in determining terminal speed.
In the homogeneous Cases 1 and 2, the lateral velocity (figure 3) exhibited a non-monotonic temporal evolution, characterised by an initial decay, followed by a brief increase before decreasing again. This behaviour reflects transient interactions among particles as they adjust from the release condition to the quasi-steady settling regime. Comparison between Case 1 (
$\rho _p/\rho _{\!f} = 1.14$
) and Case 2 (
$\rho _p/\rho _{\!f} = 1.28$
) reveals that heavier particles exhibit a larger variability in lateral velocity throughout settling, indicating greater horizontal dispersion, likely due to stronger wake-induced deflections and enhanced susceptibility to trajectory instabilities at higher Galileo number (
$\textit{Ga} = 419$
). When comparing particle locations within the Case 1 and 2 arrays, outer-layer particles (
$L_1+L_2$
, solid lines) tend to display lower lateral velocities than those in the inner layers (
$L_3+L_4+L_5$
, dashed lines), particularly during the early settling phase. Note that in Case 2, the inner-layer particles exhibited a delayed increase in velocity, indicating that the central particles experience extended mutual interference or sheltering effects before fully interacting with the surrounding fluid. However, lateral velocities exhibited a marked configuration dependence. In Case 3, lateral velocities for both particle types were comparable to those in homogeneous systems, suggesting limited inter-species interaction. In contrast, Cases 4 and 5 reveal that heavier particles, when positioned adjacent to lighter ones, exhibited reduced lateral velocities relative to Case 2. This attenuation points to configuration-sensitive hydrodynamic coupling, where the presence and placement of lower-inertia particles modify wake development and lateral dispersion of heavier neighbours.
Probability density functions (p.d.f.s) of particle vertical velocity fluctuations,
$w' = w - \overline {w}$
, normalised by the standard deviation
$\sigma _w$
, for (a)
$\rho _p/\rho _{\!f} = 1.14$
and (b)
$\rho _p/\rho _{\!f} = 1.28$
. Solid curves denote fitted distributions and coloured vertical dashed lines indicate the single-particle terminal velocity.
(a) Bulk particle configurations at two instants for the five cases (
$tw_0/d \approx 0$
and 20); blue and red markers denote lighter and heavier particles. The shaded region represents the smooth surface fitted to the particles. (b) Temporal evolution of the planar concentration
$\phi$
normalised by its initial value
$\phi _0$
. The inset shows the same data without normalisation.
Figure 4 presents the distributions of particle terminal velocities, normalised by the standard deviation after subtracting the mean, i.e.
$w'/\sigma _w$
with
$w'=w-\overline {w}$
. For the mixed-density systems (Cases 3–5), the distributions of lighter and heavier particles are shown separately in panels (a) and (b). The homogeneous arrays (Cases 1 and 2) and the laterally heterogeneous arrangement (Case 3) show broadly symmetric distributions of particle velocities. In contrast, for Cases 4 and 5, the lighter and heavier particles show slight asymmetries, with velocities shifted towards higher and lower velocities, respectively. Importantly, portions of the distributions extend beyond the dashed line corresponding to the terminal velocity of an isolated particle, indicating that certain grains settle faster than when in isolation, while others are hindered substantially. This asymmetry highlights the role of configuration-sensitive hydrodynamic interactions: lighter particles may be entrained and accelerated within the downdrafts generated by heavier neighbours, whereas heavier particles may be slowed when surrounded by lighter ones. Such non-uniform acceleration and retardation provide a mechanism for the enhanced vertical dispersion observed later in the analysis.
Figure 5(a) illustrates the three-dimensional dispersion of particles at two instants,
$tw_0/d \approx 0$
and
$tw_0/d \approx 20$
. In the homogeneous systems, dispersion occurred primarily in the lateral direction, with Case 2 (heavier particles) spreading more rapidly than Case 1. By contrast, mixed-density systems exhibited pronounced vertical spreading driven by density contrasts: Case 3 produced an oblique asymmetric pattern, while Cases 4 and 5 developed concave and convex fronts, respectively, reflecting the influence of cross-density wake interactions.
Figure 5(b) illustrates the temporal evolution of the planar concentration
$\phi$
, normalised by its value evaluated at
$t$
= 0. The planar concentration is defined as
$\phi = A_{\textit{particle}}/A_S$
, where
$A_S$
is the surface area of the smooth surface
$S$
fitted to the particle distribution using least squares regression (indicated by the shaded region in figure 5
a) and
$A_{\textit{particle}}$
is the total projected area of the particles onto
$S$
. The projected particle area is given by
$A_{\textit{particle}} = N\pi d^2/4$
, where
$N = 100$
is the number of particles. The particles are arranged initially in a 10
$\times$
10 array with an initial planar concentration of
$\phi = 0.2$
and then disperse over time, resulting in a decrease in planar concentration. In homogeneous systems, dispersion is dominated by lateral spreading, with heavier particles exhibiting a faster decay in
$\phi$
due to enhanced lateral dispersion. In contrast, mixed-density systems show additional vertical spreading driven by density heterogeneity. Despite these different dispersion mechanisms, the planar concentration in the mixed-density cases follows a trend similar to that of the homogeneous heavy particle case (Case 2), indicating that vertical separation due to density contrast contributes to planar dilution to a degree comparable to the lateral dispersion observed for heavier particles.
The basic lateral displacement (or projected dispersion) of particles during settling is quantified in figure 6. This metric is defined by projecting the particle positions onto the horizontal plane at a given time and computing the mean radial distance from the centroid of the particle cluster
$L_0 = ({1}/{n}) \sum _{i=1}^{n} \left \| \boldsymbol{x}_i - \boldsymbol{x}_c \right \|_2$
, where
$n$
is the number of particles under interrogation, and
$\boldsymbol{x}_i$
and
$\boldsymbol{x}_c$
denote the in-plane position of particle
$i$
and the instantaneous centroid of the array (figure 6
a).
(a) Diagram of particle positions projected onto the horizontal plane, with the red cross indicating the ensemble-mean particle position,
$\boldsymbol {x}_c$
, and distance,
$l_i$
, of particle
$i$
from
$\boldsymbol {x}_c$
. (b) Time evolution of bulk lateral displacement
$L_0$
for homogeneous arrangements (Cases 1 and 2) and (c)
$L_0$
for mixed-density arrangements (Cases 3–5). (d)
$L_0$
for the inner three layers (
$L_3+L_4+L_5$
) and (e)
$L_0$
for the outer two layers (
$L_1+L_2$
). Note that particles within a given layer have the same density and shaded regions indicate one standard deviation across 10 experimental realisations.
The evolution of
$L_0$
for the homogeneous cases is presented in figure 6(b), showing consistently greater lateral dispersion in Case 2 (heavier particles) relative to Case 1 (lighter particles), in agreement with the enhanced lateral velocities observed in figure 3. In the mixed-density configurations (figure 6
c), Case 3 displays higher lateral dispersion compared with Cases 4 and 5, which in turn exhibit values similar to the lighter homogeneous system (Case 1). These trends suggest that asymmetric arrangements (e.g. Case 3) can enhance lateral spreading via interfacial shear, while more symmetric, radially structured distributions (Cases 4 and 5) suppress outward bulk displacement. To further isolate the influence of particle density and spatial arrangement, we compare lateral displacement within layers (figure 1
c) of various particle groups. Figure 6(d, e) shows the mean displacement of particles in the inner three layers (
$L_3+L_4+L_5$
) and the outer two layers (
$L_1+L_2$
), defined as
$L_{\textit{inner}}$
and
$L_{\textit{outer}}$
. In figure 6(d),
$L_{\textit{inner}}$
is higher in the mixed-density cases relative to their homogeneous counterparts, indicating that particles in the centre undergo greater lateral spreading when surrounded by neighbours with differing densities. Conversely, figure 6(e) shows that
$L_{\textit{outer}}$
is reduced in the mixed particle arrangements.
(a) Raw images of settling particles at mean vertical position of the light particle cluster
$Z_p/d \approx 15$
for Cases 1 and 4. (b) Relative vertical position of each layer
$\Delta Z_{L_i}/d = (Z_{L_i} - \bar {z})/d$
, for the symmetric particle configurations (Cases 1, 2, 4 and 5).
$Z_{L_i}$
and
$\bar {z}$
denote the mean vertical position of particles in the layer
$i$
(
$L_i$
) and group mean, respectively.
The influence of particle arrangement on vertical settling dynamics was examined by analysing the variation of particle descent across the arrays and layers in each array. Particles located near the centre of the array generally settled more slowly than those near the periphery, giving rise to a characteristic concave or ‘parachute-shaped’ settling pattern. This structure is attributed to hydrodynamic interactions that inhibit the descent of centrally clustered particles. Figure 7(a) displays example snapshots of the particle array during settling for Case 1 (homogeneous) and Case 4 (mixed-density with lighter particles in the core). While both configurations exhibit the parachute pattern, the effect is more pronounced in the mixed-density case, where the lighter particles in layer 3 experience an enhanced vertical velocity due to surrounding heavier particles and associated wake interactions. To quantify this behaviour, figure 7(b) presents the relative vertical positions of particles within each layer for symmetric configurations (Cases 1, 2, 4 and 5). The relative vertical position of each concentric layer with respect to the overall mean height of the cluster
$\bar {z}$
is defined as
$\Delta Z_{L_i}/d = (Z_{L_i} - \bar {z})/d$
, where
$Z_{L_i}$
is the mean vertical position of particles within the interrogated concentric layer
$L_i$
. To improve statistical robustness, concentric layers 4 and 5, which contain fewer particles, were combined for the averaging procedure. The evolution of
$\Delta Z_{L_i}/d$
is included at three distinct stages of settling, corresponding to times when the vertical position of the particle cluster
$Z_P/d=5$
, 10 and 15, as calculated using all particles in homogeneous systems and only the lighter particles in mixed-density arrangements.
In the homogeneous particle arrangements (Cases 1 and 2; figure 7), a similar vertical stratification developed, with particles in the inner layers (
$L_3$
and
$L_{4+5}$
) lagging behind those near the periphery, consistent with the formation of a parachute-shaped structure. This trend was more pronounced with heavier particles (Case 2,
$\rho _p/\rho _{\!f} = 1.28$
), where the inner and outer layers exhibited a vertical gap of approximately
$2d$
at
$Z_P/d = 15$
, compared with a smaller gap of
$c$
.
$1d$
with the lighter particle cluster (Case 1,
$\rho _p/\rho _{\!f} = 1.14$
). Once the parachute-like structure forms, the relative arrangement of particles remains largely unchanged, particularly in the case of heavier particles. These heavier particles settle at larger
$\textit{Ga}$
and particle
$Re$
, producing more energetic unsteady wakes with enhanced vortex shedding. These inertial wakes increase wake overlap and entrainment within the finite array, strengthening the peripheral downdrafts and the associated interior recirculation. The resulting wake shielding and increased hydrodynamic drag on the central particles amplify the parachute-like curvature, which exert an influence on the confined central region of the array.
In the mixed-density configurations (Cases 4 and 5; figure 7), layer-wise vertical separation became much more pronounced due to the contrasts in particle density. In Case 4, lighter particles in layers 3–5, located at the centre of the array, settled significantly more slowly than their heavier neighbours in the outer layers. This disparity persisted and grew over time, leading to an increasing separation between the central and peripheral layers. The observed behaviour highlights a fundamental distinction from the central particles in the homogeneous scenarios, where such sustained divergence was absent. At
$Z_P/d = 15$
, the vertical gap between the inner (lighter) and outer (heavier) particles reached approximately
$10d$
, compared with only
$1d$
in the homogeneous light-particle case (Case 1). Case 5, however, exhibited an inversion of the parachute pattern. Here, the heavier particles are concentrated in the central layers (
$L_3$
to
$L_5$
) with the lighter ones at the periphery; the inner particles thus settled more rapidly than the outer ones, effectively reversing the concavity of the settling front.
Figure 8 illustrates a sequence of flow fields in the central vertical plane for homogeneous Case 1 and heterogeneous Case 4, shown at characteristic stages of descent corresponding to array depths
$Z_P/d = 10$
, 20, 30, 40 and 50. Two distinct phenomena are evident from a qualitative comparison. First, a comparison of the light particles present in both scenarios reveals greater dispersion in the central region of the heterogeneous case, which is attributed to the influence of wake interactions generated by the surrounding heavier particles. While this observation is qualitative, it is supported by subsequent quantitative analysis of directional pair dispersion. Second, the velocity field in Case 1 displays relatively homogeneous and symmetric fluctuations, indicative of roughly uniform wake development across the array. In contrast, Case 4 reveals significant flow heterogeneity associated with the imposed particle density gradient. A persistent shear layer forms at the interface between the heavier outer particles and the core of lighter grains, giving rise to spatially localised vortical structures (highlighted by the red arrows in figure 8 for
$Z_p/d = 20$
). These structures propagate inward and impinge upon the central region, thereby increasing mixing and enhancing the lateral and vertical spread of lighter particles. The emergence of these vortical structures, and their coupling with the particle distribution, suggest that the flow in heterogeneous systems reflects complex, configuration-dependent, hydrodynamic interactions that strongly modulate dispersion and particle clustering as shown later.
Flow field sequences for Case 1 (a) and Case 4 (b) at array depths of
$Z_P/d = 10$
, 20, 30, 40 and 50. Red arrows highlight the pronounced shear layer generated by the particle density contrast and the resulting coherent motions.
Instantaneous velocity fields superimposed with vorticity contours in the vertical mid-plane for all five configurations. Schematics denote initial particle densities: blue for
$\rho _p/\rho _{\!f} = 1.14$
and red for
$\rho _p/{}\rho _{\!f} = 1.28$
.
(a) Temporal evolution of the vertical velocity field,
$w/w_0$
, in the mid-plane at
$z/d = 10$
for homogeneous Case 1 and heterogeneous Cases 3 and 5.
$\tau$
is the time starting when
$Z_P/d = 10$
. (b) Iso-contours of vertical velocity (
$w/w_0$
= 0.06) of wake flows over time at
$z/d=25$
for Case 3, illustrating various interaction mechanisms.
To further illustrate flow organisation induced by the settling particles, figure 9 shows instantaneous velocity fields superimposed with out-of-plane vorticity contours in the vertical mid-plane for all five configurations when
$Z_p/d \approx 25$
. Across all cases, the descending particles generate superimposed wake dynamics, with vortex shedding, coherent shear layers and localised regions of comparatively high vorticity. These flow structures are sensitive to the initial particle arrangement and contrast in grain density. The coexistence of light and heavy particles leads to localised intensification of shear and spatially heterogeneous wake interactions, particularly along the interfaces between particle types.
A more detailed examination of the flow evolution is presented in figure 10(a), which shows the temporal progression of vertical velocity at a fixed height
$z/d = 10$
for homogeneous Case 1 (light particles), and heterogeneous Cases 3 and 5. These plots highlight the development and interaction of particle wakes, the emergence of spatial uniformity and variations in wake width over time, with
$\tau$
being referenced to the instant when
$Z_P/d = 10$
. In Case 1, the vertical velocity field becomes relatively uniform across the horizontal span by
$\tau w_0/d \sim 100$
, indicating strong wake interactions and effective lateral mixing within the particle array. A similar trend is observed in Case 5, which features light particles surrounding a heavier particle core; despite the initially non-uniform velocity field at earlier times (
$\tau w_0/d \sim 20$
–50), the flow becomes homogenised at later stages as the particles settle. Note that Case 5 displays a pronounced narrowing of the wake cluster over time, indicating a reduction in lateral spreading, which may be attributed to the influence of the central heavy particles organising the surrounding flow. In contrast, Case 3 exhibits a distinctly different evolution. Here, the velocity field remains horizontally asymmetric throughout the duration shown, reflecting the persistent influence of the sharp lateral interface in particle density. This asymmetry suggests that the presence of a strong density contrast alone does not guarantee enhanced wake mixing or symmetry. In this case, the faster-settling heavy particles leave behind residual wake structures that interact asymmetrically with the wakes of slower-moving lighter particles, thereby limiting the degree of flow development across the interface. Figure 10(b) shows a spatial reconstruction illustrating the time evolution of vertical velocity isosurfaces at the plane
$z/d = 25$
for Case 3, obtained from stereoscopic PIV. This reconstruction captures key flow features such as individual particle wakes, wake merging and the emergence of complex superimposed wake structures. Complementing this, figure 11 provides stereoscopic views in the
$z/d = 25$
plane, where in-plane velocity vectors are overlaid and colour-coded by the out-of-plane velocity component. These snapshots highlight the relatively uniform and well-mixed wake field in homogeneous Case 1 at two times, in contrast to the heterogeneous arrangements (Cases 3 and 5), which exhibit persistent asymmetry and reduced wake homogenisation.
We now quantify the temporal evolution of the bulk mean-squared pair dispersion for all five configurations, decomposing it into total three-dimensional (
$R^2$
), lateral (
$R_L^2$
) and vertical (
$R_z^2$
) components. These calculations consider all initial particle separations to capture the ensemble-averaged collective behaviour. In the mixed-density configurations, pair dispersion is evaluated separately for lighter and heavier particles to isolate species-specific dynamics and cross-interaction effects. The results, summarised in figure 12, reveal distinct dispersion regimes that are governed by particle density and initial spatial arrangement. In all configurations, the early-time lateral spreading exhibits a ballistic regime, with
$R_L^2 \propto t^2$
. This scaling reflects the initial acceleration of particles from rest and the influence of wake-mediated interactions, rather than classical Batchelor scaling in homogeneous turbulence (Batchelor Reference Batchelor1950). For lighter particles, the initial ballistic regime transitions towards a diffusive-like scaling,
$R_L^2 \propto t^1$
, as their reduced inertia limits the persistence of coherent trajectories. In the vertical direction, dispersion grows faster than diffusion, with
$R_z^2 \propto t^{3/2}$
, driven primarily by gravitational settling and wake-induced interactions. The duration and magnitude of this regime depend on particle configuration. In homogeneous Case 1, composed entirely of light particles, vertical dispersion rapidly transitions to a diffusive-like regime due to uniform settling. In contrast, Case 5 maintains faster than diffusive vertical spreading for longer times, as lighter peripheral particles are continuously entrained by the downdraft produced by heavier central particles.
Instantaneous velocity fields in the horizontal plane at
$z/d = 25$
for Cases 1, 3 and 5 at (a)
$\tau w_0/d = 30$
and (b)
$\tau w_0/d =100$
after the particles cross the plane. Arrows depict the horizontal velocity component, while colour shows the vertical velocity
$w$
. The
$w_0$
corresponds to the heavier particle terminal velocity.
Bulk pair dispersion of the particles in 3-D space,
$R^2/d^2$
, lateral,
$R_L^2/d^2$
, and vertical,
$R_z^2/d^2$
, directions. For mixed-density particle systems, pair dispersion was calculated separately for each particle type.
Pair dispersion in the lateral (
$R_L^2/d^2$
) and vertical (
$R_z^2/d^2$
) directions for varying initial separations
$r_0/d$
. Top panels, homogeneous Case 1. Bottom panels, lighter particles in heterogeneous Case 4. Insets show the scaled dispersion using
$\varphi _L$
and
$\varphi _z$
.
Normalised pair dispersion scaling coefficients (a)
$\varphi _L$
and (b)
$\varphi _z$
as a function of initial spacing
$r_0/d$
for homogeneous Case 1 and heterogeneous Case 4.
To further examine the influence of initial particle separation, we analyse the dependence of dispersion on the initial centre-to-centre spacing
$r_0/d$
. For both lateral and vertical directions, dispersion is expressed as
\begin{equation} R_L^2(t, r_0) = \big \langle \left | \delta l(t) - \delta l(0) \right |^2 \big \rangle , \quad R_z^2(t, r_0) = \big \langle \left [ z(t) - z(0) \right ]^2 \big \rangle , \end{equation}
where
$\delta l$
denotes the lateral particle distance in the
$x$
–
$y$
plane from a reference origin, and the averages are conditioned on the initial separation
$\delta l(0)$
and
$z(0)$
. A comparison of homogeneous Case 1 and heterogeneous Case 4 (figure 13) reveals a distinct dependence on initial separation. In Case 1, lateral dispersion increases with
$r_0$
, whereas vertical dispersion remains largely independent of
$r_0$
. In contrast, Case 4 displays a monotonic increase in both lateral and vertical dispersion with increasing
$r_0$
, reflecting stronger sensitivity to initial configuration in the heterogeneous arrangement. To assess these trends, we explore the normalised dispersion ratio
\begin{equation} \varphi _i(r_0) = \frac {R_i^2(r_{0 \textit{min}})}{R_i^2(r_0)}, \quad i = L,\ z, \end{equation}
where
$r_{0, \textit{min}}$
denotes the smallest initial particle separation in the array corresponding to nearest-neighbour pairs. This normalisation is introduced as a reference scale to compare the relative growth of separations within a given finite configuration. The quantities
$\varphi _L$
and
$\varphi _z$
are shown as functions of
$r_0/d$
in figure 14. In the homogeneous configuration,
$\varphi _L$
exhibits a systematic decay with increasing separation over the accessible range, while
$\varphi _z \approx 1$
for all
$r_0$
, indicating weak sensitivity of vertical dispersion to initial spacing. In the heterogeneous configuration, both lateral and vertical components show a stronger, non-trivial dependence on
$r_0$
, reflecting the influence of heterogeneity and interactions between particles of different densities. When the dispersion curves
$R_i^2(t)$
are rescaled by
$\varphi _i(r_0)$
, a collapse is observed across different initial separations (see insets in figure 13). This normalisation removes the leading-order
$r_0$
dependence and isolates the temporal evolution of dispersion, providing a compact and generalisable factor for characterising pairwise separation dynamics in heterogeneous suspensions.
4. Discussion
The experiments herein focus on the early to intermediate settling dynamics of different arrangements of particles of variable densities, and demonstrate that collective settling is strongly influenced by distinct particle-particle interactions and resulting flow dynamics, which are dependent on their initial configuration and density contrast. Homogeneous arrays (Cases 1 and 2) consistently develop a self-organised, concave parachute-like settling pattern, where outer particles descend faster than those at the core due to asymmetric wake interactions, consistent with clustering and preferential path formation observed in dilute suspensions (Nicolai et al. Reference Nicolai, Herzhaft, Hinch, Oger and Guazzelli1995; Guazzelli & Hinch Reference Guazzelli and Hinch2011). In mixed-density arrays (Cases 3–5), this structure is significantly modified by spatial heterogeneity in the arrangement of different density particles. When heavier particles surround a lighter core (Case 4), their rapid descent generates a curtain of downward flow and interfacial shear, increasing the settling velocity of the light particles at the edge of the central core (layer 3), thereby amplifying the curvature of the parachute-shaped settling pattern (figure 7 a). In contrast, when particles form a core of denser grains (Case 5), they produce a downdraft that accelerates the surrounding lighter particles, inverting the vertical profile and causing the core to lead the descent (figure 7 b). These results complement prior studies, which have largely focused on monodisperse suspensions (Uhlmann & Doychev Reference Uhlmann and Doychev2014), by revealing how cross-density wake interactions can either suppress or enhance collective settling depending on the initial arrangement.
Pair-dispersion analysis shows that vertical spreading dominates due to gravity and collective wake interactions, while lateral dispersion transitions from a ballistic (
${\sim} t^2$
) to diffusive-like (
${\sim} t^1$
) behaviour as velocity correlations decay. In homogeneous arrays, vertical dispersion is nearly independent of initial separation
$r_0$
, indicating uniform collective acceleration. In contrast, mixed-density arrays exhibit strong
$r_0$
dependence in both the lateral and vertical directions, with greater separations for initially distant pairs due to enhanced local shear and entrainment.
These phenomena show how interfacial shear and collective entrainment establish the characteristic length and time scales of dispersion in settling particle groups. For example, the thickness of the ‘parachute’ layer and duration of coherent motions are both controlled by the strength of shear at interfaces between the heavy and light particles, and by the rate at which the particle ensemble entrains ambient fluid downward. More broadly, the emergent mesoscale structures observed here, clusters of particles, vertical jets of flow and superimposed wakes, arise from multiscale particle–fluid interactions that are not captured by single-particle or two-particle models. Classical descriptions of isolated settling or pairwise drafting–kissing–tumbling cannot predict phenomena such as a heavy particle ring effectively acting as a hydrodynamic barrier that filters and slows the interior flow, or a core of heavier particles generating a sustained downdraft that accelerates neighbouring grains. Key wake-mediated mechanisms (e.g. wake drafting, vortex-induced repulsion and flow shielding) underlie these collective behaviours, leading directly to outcomes such as enhanced vertical dispersion of lighter particles and suppressed lateral spread of heavier grains in mixed configurations. In quiescent environments, large groups of particles can thus self-organise into clusters (Metzger, Guazzelli & Butler Reference Metzger, Guazzelli and Butler2005), stratified layers and other mesoscale patterns purely through hydrodynamic coupling. Density heterogeneity primarily modifies the transient mesoscale organisation of the settling particle layer through configuration-dependent wake coupling and entrainment.
These measurements also give a contextual comparison with the Richardson–Zaki form by treating the observed bulk settling velocity in terms of an effective volume fraction. Using this representation, an effective exponent
$n \approx 1.6$
is obtained at
$\phi _{V_0} = 0.13$
, providing a phenomenological measure of drag modification by hydrodynamic interactions in a finite, initially planar particle ensemble.
The velocity p.d.f.s (figure 4) show a slight bias toward higher velocities for the lighter particles, while for the heavier particles, the distributions are weighted more towards lower velocities, consistent with configuration-sensitive hydrodynamic coupling: we show that lighter particles can be entrained and accelerated in the downdrafts of their heavier neighbours, whereas heavier particles experience increased drag when embedded among lighter ones. These asymmetries align with the enhanced vertical dispersion of lighter particles and hindered settling of heavier ones observed in figures 5 and 12.
5. Conclusions
This study demonstrates that collective particle settling in a quiescent fluid is governed by configuration-sensitive hydrodynamic interactions, which produce coherent structures and anisotropic dispersion patterns. Using systematic experiments with large arrays of inertial particles, we show that both spatial arrangement and density heterogeneity influence the dynamics of settling significantly, generating distinct mesoscale features that differ from predictions from single-particle or dilute models.
The results reveal the emergence of particle stratification and clustering driven by wake interactions. Homogeneous arrays spontaneously develop organised descent patterns with coherent velocity gradients, while mixed-density systems exhibit symmetry breaking due to wake interference between grains and regions of different particle density. These phenomena result from nonlinear flow–particle coupling and cannot be inferred from linear superposition or isolated-particle behaviour.
The results indicate an effective Richardson–Zaki exponent
$n \approx 1.6$
for an effective volume fraction
$\phi _{V_0} = 0.13$
, reflecting the reduction in mean particle settling velocity resulting from collective wake interactions and particle shielding. Although the present configurations do not represent a classical hindered settling suspension, the results indicate that collective hydrodynamics influence the mean settling rate within a finite, initially planar ensemble.
Pair dispersion analysis reveals directionally dependent regimes: vertical spreading dominates due to gravitational forcing coupled with unsteady wake-induced fluctuations, while lateral dispersion reflects decorrelation dynamics linked to the array geometry. Homogeneous particle systems exhibit weak dependence on initial inter-particle spacing, consistent with coherent collective motion, whereas heterogeneous systems display strong sensitivity to spacing due to local shear and entrainment differences. Normalised dispersion ratios
$\varphi _{L,z}(r_0) = R_{L,z}^2(r_0)/R_{L,z}^2(r_{0, \textit{min}})$
reveal a collapse across different separations as a function of
$r_0/d$
, highlighting the presence of a scaling linked to the array geometry and density distribution. Measurements of the velocity field confirm that coherent flow features, such as jets, wake interactions and entrainment layers, are strongly configuration-dependent and central to the spatio-temporal evolution of the particle distributions. These results provide insight into wake-mediated interactions and density heterogeneity, and offer a benchmark for validating high-fidelity multiphase simulations.
The findings offer a clear benchmark for validating multiphase simulations, and point towards future work involving larger particle systems, additional settling layers and more complex particle or flow conditions. A natural extension of this work involves expanding the particle count and domain size, and with multiple settling layers, which will enable exploration of statistical clustering, network-like connectivity and hydrodynamic percolation at larger scales. Finally, consideration of non-spherical particles, polydisperse systems and imposed flow conditions (e.g. shear, oscillation) will broaden the applicability of these findings to natural and engineered particulate flows.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2026.11570.
Declaration of interests
The authors report no conflict of interest.
Open access
