Breakup and coalescence of particle aggregates at the interface of turbulent liquids

Abstract

The breakup and coalescence of particle aggregates confined at the interface of turbulent liquid layers are investigated experimentally and theoretically. In particular, we consider conductive fluid layers driven by Lorentz forces and laden with millimetre-scale floating particles. These form aggregates held together by capillary attraction and disrupted by the turbulent motion. The process is fully characterised by imaging at high spatio-temporal resolution. The breakup frequency $\varOmega$ is proportional to the mean strain rate and follows a power-law scaling $\varOmega \sim D^{3\text{/}2}$, where $D$ is the size of the aggregate, attributed to the juxtaposition of particle-scale strain cells. The daughter aggregate size distribution exhibits a robust U-shape, which implies erosion of small fragments as opposed to even splitting. The coalescence kernel $\varGamma$ between pairs of aggregates of size $D_{1}$ and $D_{2}$ scales as $\varGamma \sim ( D_{1} + D_{2} )^{2}$, which is consistent with gas-kinetic dynamics. These relations, which apply to regimes dominated both by capillary-driven aggregation and by drag-driven breakup, are implemented into the population balance equation for the evolution of the aggregate number density. Comparison with the experiments shows that the framework captures the observed distribution for aggregates smaller than the forcing length scale.

1. Introduction

The breakup and coalescence of large particle aggregates/clusters due to the hydrodynamic stresses in turbulence are highly relevant to a wide range of natural and industrial applications, including the sedimentation of marine snow (Burd & Jackson 2009), wastewater treatment (Samstag et al. 2016), the fragmentation of macro-plastics in the ocean (Van Sebille et al. 2020), the growth of rain droplets in atmospheric clouds (Shaw 2003) and the formation of planets by dust grains around young stars (Blum & Wurm 2008). Both processes are crucial for the formation and temporal evolution of the aggregates/clusters (Liu et al. 2020; Ma et al. 2023) and thus for the physical properties of the system. In this work, we focus on particle ‘aggregates’ which indicates groups of particles that adhere to each other and merge into larger entities, rather than ‘clusters’ which is often used to denote groups of particles in close proximity (Monchaux, Bourgoin & Cartellier 2012).

The breakup and coalescence of aggregates in the bulk of fluid flows have been studied extensively. Attempts have been made to estimate the hydrodynamic stresses sufficient to break aggregates of given properties (Kobayashi, Adachi & Ooi 1999; Harshe, Lattuada & Soos 2011). This approach, however, does not lead to an explicit expression for the breakup frequency, which is crucial to aggregate evolution and size distribution. On the other hand, approaches that compare the aggregate strength and the fluctuating turbulent stresses lead to forms of the breakup frequency varying from power law (Pandya & Spielman 1982) to exponential (Kusters 1991) and multifractal (Bäbler et al. 2008, 2012). These analytical models, and similar ones focused on coalescence, have been compared with experiments in which the size distributions of aggregates are evaluated as a function of particle properties and shear rate. Such comparisons, however, are indirect in that the breakup and coalescence events are not usually observed.

In order to model coalescence in turbulence, Saffman & Turner (1956) proposed a framework for the collision frequency of non-inertial particles. This was later adapted to particles with finite inertia (described by a Stokes number ${\textit{St}}$ that compares their response time with the time scale of the flow) to account for their inhomogeneous spatial distribution and increased relative velocity (Sundaram & Collins 1997; Wang, Wexler & Zhou 2000; Pumir & Wilkinson 2016; Brandt & Coletti 2022; Bec, Gustavsson & Mehlig 2024). While most studies have considered particles as non-interacting point masses, recent numerical simulations have enabled more direct insight into the coalescence and breakup of aggregates by modelling their cohesive forces (De Bona et al. 2014; Bäbler et al. 2015), including their two-way coupling with the fluid and with each other (Yao & Capecelatro 2021), and even resolving the flow at the particle scale (Vowinckel et al. 2019).

A widely adopted framework to model aggregate evolution is the population balance equation. This has found extensive application not only in particulate systems such as gas–particle mixtures (Marchisio & Fox 2013) and cohesive sediment (Vowinckel et al. 2019), but also in suspensions of bubbles and droplets (Martínez-Bazán et al. 1999; Liao & Lucas 2009; Martínez-Bazán et al. 2010; Aiyer et al. 2019; Qi et al. 2022; Ruiz-Rus et al. 2022; Ruth et al. 2022; Ni 2024; Qi et al. 2024). Considering a spatially homogeneous system in which advection is negligible, the rate of change in time $t$ of the number density $n$ of aggregates of size $D$ is set by the balance between breakup and coalescence:

(1.1) \begin{align} \frac {\partial n(D,t)}{\partial t} & = \int _{D}^{\infty }{m\left ( D_{i} \right )\beta \left ( D,D_{i} \right )\varOmega (D_{i})n(D_{i},t)\text{d}D_{i}} - \varOmega (D)n(D,t) \nonumber \\&\quad + \frac {1}{2}\int _{0}^{D}{\varGamma \Big ( D_{i},\big ( D^{2} - D_{i}^{2} \big )^{1/2} \Big )}n\left ( D_{i},t \right )n\Big ( \big ( D^{2} - D_{i}^{2} \big )^{1/2},t \Big )\text{d}D_{i}\nonumber \\&\quad - n(D,t)\int _{0}^{\infty }{\varGamma \left ( D,D_{i} \right )n\left ( D_{i},t \right )\text{d}D_{i}} , \end{align}

where $D=\sqrt {A}=\sqrt {N}d_p$ is the size of an aggregate of $N$ particles of diameter $d_{p}$ , $A$ is the area of the aggregate, $m$ is the number of daughter aggregates generated in a breakup event, $\beta$ is the daughter aggregate size distribution, $\varOmega$ is the breakup frequency and $\varGamma$ is the coalescence kernel between two aggregates. The four terms on the right-hand side of (1.1) then represent the source term due to breakup of larger aggregates, the loss term due to breakup of aggregates of size $D$ , the source term due to coalescence among smaller aggregates and the loss term due to coalescence of aggregates of size $D$ with other aggregates, respectively.

Here we focus on the formation and evolution of particle aggregates at the interface of turbulent liquids. In such systems, particles behave and interact in profoundly different ways with respect to those in the bulk (Garbin 2019; Magnaudet & Mercier 2020; Protière 2023). We recently examined this class of flows by studying experimentally the interactions of millimetric, monodisperse particles at the interface of thin turbulent liquid layers (Shin & Coletti 2024, 2025; Shin, Stricker & Coletti 2025). This revealed distinct regimes characterised by different propensity to form aggregates of various sizes and mobility. In the present study, we leverage the population balance framework to gain insight into the breakup and coalescence of aggregates in such a system. Results for breakup frequency, daughter aggregate size distribution and coalescence kernel are presented and modelled. These are further implemented in the population balance equation to predict the evolution of the aggregate number density.

2. Experimental method

We analyse data presented in Shin & Coletti (2024) and Shin et al. (2025), obtained with an apparatus described in detail in Shin, Coletti & Conlin (2023); therefore, the set-up and the methodology are only briefly summarised here. A conductive layer of 10 % $\textrm {CuSO}_{4}$ aqueous solution by mass (with density $\rho _{\!f} = 1.08$ g cm⁻³ and kinematic viscosity $\nu = 1.0\ \times 10^{- 6}$ m² s⁻¹) is contained in a $320\,{\rm mm} \times 320$ $\textrm {mm}$ tray above an $8 \times 8$ checkerboard array of permanent magnets with alternating polarity. A direct current driven between copper electrodes immersed in the fluid generates horizontal Lorentz forces driving an approximately two-dimensional turbulent flow, with a forcing scale $L_{\!f}$ = 35 mm given by the spacing between the magnets. We employ single-layer and double-layer fluid configurations: the former featuring a 7 mm deep conductive layer; and the latter with a 2 mm thick mineral oil layer over an 8 mm deep conductive layer. The main difference between the single-layer and double-layer configurations is represented by the lower interfacial tension $\gamma$ in the latter. Moreover, in the double-layer configuration, the presence of the oil layer over the conductive fluid causes a different friction coefficient along the bottom wall below the conductive fluid, as reported, for example, in Shin et al. (2023). The present range of parameters, however, warrants that the turbulence properties, such as velocity structure functions and energy spectra, behave similarly in all considered cases, and closely follow the expected behaviour of two-dimensional turbulence, as shown by Shin et al. (2023, 2025).

We use two sizes of monodispersed polyethylene spherical particles of diameter $d_{p} =$ 1.1 and 1.8 mm (comparable to the Kolmogorov scale $\eta$ ) and density $\rho _{p} = 0.96$ and 1 g cm⁻³, respectively, which float at the free surface in the single-layer configuration and reside at the liquid–liquid interface in the double-layer one. The Stokes number ${\textit{St}}$ and Bond number $Bo$ (defined in table 1 along with other important non-dimensional parameters) are much smaller than unity, indicating that both particle inertia and interfacial distortion due to buoyancy are negligible. The particles are imaged by a CMOS camera recording at 100–150 frames per second, ensuring inter-frame displacements of about 1 mm. Lagrangian particle tracking with subpixel accuracy resolves individual trajectories. Separate experiments featuring only $75 {-} 90$ μm tracers characterise the root-mean-square fluid velocity $u_{\textit{rms}}$ that sets the turbulence Reynolds number ${\textit{Re}}$ . We note that the aggregates formed by the spherical particles likely alter the two-dimensional turbulent flow underneath. This two-way coupling, which cannot be directly evaluated in this set-up, is expected to be subdominant at the present relatively dilute concentrations: aggregates collide but are rarely locked into large structures that would significantly backreact on the energy-containing eddies. Thus, in this work, we neglect to a first approximation the two-way coupling and use flow statistics based on single-phase experiments in our aggregate dynamics model. This assumption was also used in Shin & Coletti (2024) who could successfully build a phase diagram describing the aggregate dynamics.

Table 1. The main parameters for the considered cases in this work. Here, $St\equiv \tau _p u_{\textit{rms}}/L_{\!f}$ with $\tau _p=\rho _pd_p^2/(18\rho _f\nu )$ being the particle response time, $Bo\equiv ((\rho _f-\rho _p)gd_p^2)/(4\gamma )$ and ${\textit{Re}}\equiv u_{\textit{rms}}L_{\!f}/\nu$ .

In this work, for both single-layer and double-layer configurations, the particle dynamics is dominated by the balance among three mechanisms, i.e. capillary attraction, viscous drag and lubrication, whose magnitudes depend on the area fraction $\phi$ and the capillary number ${\textit{Ca}}$ (Shin & Coletti 2024; Shin et al. 2025). Here, $\phi$ is defined as the ratio between the footprint of all particles in the domain and the area of the same domain. Ratio ${\textit{Ca}}$ is the ratio of viscous drag pulling particles apart to the capillary attraction acting at contact, ${\textit{Ca}}=6 \sqrt {2} \pi d_{p}^{6} u_{r m s}/ (\varTheta L_{\!f} )$ , where $\varTheta =12 \gamma h_{qp}^{2} d_{\textit{p}}^{3} / (\rho _{\!f} \nu )$ and $h_{pq}$ is the amplitude of the quadrupole mode of interfacial distortion (as measured in separate experiments by Shin & Coletti (2024)). The dynamics was found to be consistent for both single-layer and double-layer configurations. Here, we analyse three cases of nearly identical area fraction $\phi$ and varying capillary number ${\textit{Ca}}$ . The three cases span the transition from a regime of capillary-driven clustering ( ${\textit{Ca}}\lt 1$ ) to one of drag-driven breakup ( ${\textit{Ca}}\gt 1$ ) (Shin & Coletti 2024; Shin et al. 2025). In the former, capillary attractions dominate over viscous drag and hit-and-stick collisions form aggregates of complex shape and wide size distribution. In the latter, the drag exerted by strong turbulent fluctuations overcomes capillary attraction, causing particles to remain isolated or form small aggregates.

The particle aggregates are identified in this study as the particles with minimum centre-to-centre distance smaller than a threshold. This threshold is determined by first plotting the distance from the nearest neighbour sorted in descending order and then identifying the elbow in such a distribution (see details in Shin et al. (2025)). This method leads to thresholds of $1.1d_p$ and $1.4d_p$ for small and large particles, respectively. The thresholds are approximately one order of magnitude larger than the pixel size resolved by the camera. The aggregates are followed over time along with their member particles. The addition/detachment of any member of an aggregate corresponds to a breakup/coalescence event. Aggregates that touch the boundary of the field of view are removed from the statistics as their sizes cannot be reliably determined. Figure 1 shows an example of detected aggregates overlapped with the raw image, with one breakup and two coalescence events of particle aggregates.

Figure 1. Detected particle aggregates on the liquid interface overlapped with the raw image for ${\textit{Ca}}=0.16$ . The circles represent the detected particle locations and different colours represent different aggregates. From $t=0$ to 0.07 s, the red aggregate breaks into purple and orange aggregates. From 0.07 s to 0.2 s, the orange aggregate coalesces with the small pink one above it and merges into the light blue aggregate. Simultaneously, the blue and green aggregates merge into the yellow aggregate.

Figure 2. (a) The breakup frequency $\varOmega$ as a function of the aggregate size $D$ for different ${\textit{Ca}}$ . The dashed line denotes $3/2$ power-law scaling. (b) A schematic to illustrate the breakup of an aggregate of size $D$ due to the strain cell. The black boxes mark the potential position of strain cells of size $d_s$ that could break the aggregate. This picture is similar to the definition of the box-counting dimension of the aggregate.

3. Results

3.1. Breakup frequency

In all cases included in this study, more than 85 % of the breakup events produce only two daughter aggregates ( $m$ = 2). Thus, only binary breakups are considered in this work. Figure 2(a) shows the normalised breakup frequency $\varOmega L_{\!f}/u_{\textit{rms}}$ as a function of the aggregate size $D$ . The breakup frequency increases as the aggregate size grows, consistent with a power-law scaling of exponent $3/2$ . Conditioning the data for different ranges of daughter aggregate sizes (not shown) indicates that the power-law scaling remains consistent. In addition, as ${\textit{Ca}}$ increases, so does the breakup frequency.

To interpret the data, a phenomenological model for the breakup frequency is proposed. We assume the rupture of an aggregate is due to a local strain event with a length scale $d_{s}\sim \mathcal{O}(d_{p})$ , acting to separate two adjacent particles. In the following discussion this is referred to as a strain cell, which will break apart adjacent particles when the associated strain rate is sufficiently strong to overcome the capillary attraction. Because such a strain cell might occur randomly within the aggregate, the breakup frequency is expected to be proportional to the number of strain cells needed to fill the aggregate of size $D$ , as illustrated in figure 2(b). This picture is analogous to the definition of the box-counting dimension of an object, i.e. counting the number of boxes required to cover it. As the aggregates in this system have a well-defined box-counting dimension $\alpha$ (Shin et al. 2025), the number of strain cells to cover the aggregate area can be written as $N_{s}(D) = Cd_{s}^{- \alpha }$ , where $C$ is a constant depending on $D$ .

Now consider another aggregate with a different size $cD$ , where $c$ is a scale factor. Similarly, the breakup frequency of this aggregate is still proportional to the number of strain cells (denoted as $N_s (cD)$ ) of size $d_s$ needed to cover it. In order to calculate $N_s (cD)$ , we notice that $N_s (cD)$ is equivalent to the number of boxes of size $d_s/c$ required to cover a hypothetical, scaled aggregate of the same geometry but with a size $D$ . Assuming the fractal dimension remains scale-independent, this leads to $N_s (cD)=C(d_s/c)^{-\alpha }=Cc^\alpha d_s^{-\alpha }$ . We note that $N_s (cD)=c^\alpha N_s (D)$ . This suggests that $N_s (D)$ follows a power-law scaling $N_s (D)\sim D^\alpha$ , leading to the breakup frequency $\varOmega \sim D^{\alpha }$ . As Shin et al. (2025) observed $\alpha \approx 1.5$ for aggregates smaller than $L_{\!f}$ (which account for the vast majority in the considered range of concentration), the model is consistent with the data in figure 2(a).

A few considerations are in order. First, we have considered strain cells of similar size $d_{s}\sim \mathcal{ O}(d_{p})\sim \mathcal{O}(\eta )$ . Significantly larger strain cells might also cause breakup. However, since in two-dimensional turbulence the strain rate is $O(u_{\textit{rms}}/L_{\!f})$ at all scales (Davidson 2015), these are comparable in strength and fewer in number than the smaller ones. Therefore, particle-scale strain cells are more likely to dominate breakup. Second, not all strain cells are sufficiently strong to break inter-particle capillary bonds, especially at low ${\textit{Ca}}$ : the probability of a local strain cell breaking an aggregate is a monotonically growing function $p(Ca)$ approaching unity for large ${\textit{Ca}}$ . This phenomenology implies that the breakup frequency is also proportional to the frequency of occurrence of the strain cells, which is expected to scale with the strain rate $u_{\textit{rms}}/L_{\!f}$ . Overall:

(3.1) \begin{equation} \varOmega \sim u_{\textit{rms}}/L_{\!f}p(Ca)D^{3/2}d_p^{-3/2}, \end{equation}

which is consistent with the data in figure 2(a). Here and in the following, the aggregate size is normalised by $d_{p}$ . As the considered cases have comparable $L_{\!f}/d_{p}$ , this cannot be meaningfully contrasted with the alternative normalisation by the forcing scale.

3.2. Daughter aggregate size distribution

We now consider $\beta$ , i.e. the probability distribution function of the normalised area of the daughter aggregate, $A_d/A_m$ , for various ${\textit{Ca}}$ and mother aggregate sizes. Here $A = {{Nd}_{p}}^{2} = D^{2}$ is the area of the aggregate, with subscripts indicating daughter and mother aggregates, respectively. As discussed, only binary breakups are considered.

Figure 3 shows the daughter aggregate size distribution, with mother aggregates with ${A_{m}}/{d_{p}^{2}} \lt 10$ excluded from the statistics to ensure a smooth distribution. The results are consistent with a universal U-shape independent of ${\textit{Ca}}$ . This means that, for all cases, aggregate breakup happens by erosion of small fragments rather than by splitting in parts of comparable size. By dividing the data into subgroups with $A_{m}/d_{p}^{2} \lt 50$ and $A_{m}/d_{p}^{2} \geqslant 50$ , the daughter size distribution does not appear affected by the mother size (for the example case ${\textit{Ca}} = 1.2$ , the other cases behaving analogously). These findings are similar to those reported for bubble breakup (Hesketh, Etchells & Russell 1991; Qi et al. 2020) in which a U-shape daughter size distribution and independence of mother size are also observed.

Figure 3. The daughter aggregate size distribution for various ${\textit{Ca}}$ and mother aggregate sizes.

We rationalise this observation by considering a simple scenario in which a mother aggregate of area $A_{m}$ breaks into two daughter aggregates, the smaller having area $A_{d} \ll A_{m}$ . When a strain cell initialises the fracture, the fracture grows in such a way that the area of the smaller daughter aggregate $A_d$ is minimised. This occurs because the smaller daughter aggregate has smaller inertia (proportional to its area) and therefore is easier to push away from the mother aggregate, leading to a higher probability. One can thus hypothesise $\beta \sim (A_d/A_m)^{-1}$ for $A_{d} \ll A_{m}$ . As $\beta$ is symmetric around ${A_{d}}/{A_{m}} = 0.5$ due to mass conservation, we expect

(3.2) \begin{equation} \beta \sim \left ( A_{d}/A_{m} \right )^{- 1} + \left (1 - A_{d}/A_{m} \right )^{- 1}. \end{equation}

The data in figure 3 agree well with this model (plotted with a numerical pre-factor to subtend a unit area).

We note that the U-shape distribution can also be interpreted using simple energy-based arguments. Considering the breakup process is indeed separating the capillary bonds between particles, the energy required for a breakup event is proportional to the number of bonds to be separated. This number of bonds further scales with the size of the smaller fragment. Therefore, the probability for breakup events producing small fragments is higher as the corresponding energy requirement is lower (i.e. fewer bonds need to be separated).

3.3. Coalescence kernel

Similarly to breakup, more than 90 % of coalescence events involve only two aggregates merging into one; thus, we consider only binary coalescence. We also note that for particle aggregates in this work, the coalescence efficiency is unity, implying all collisions lead to coalescence. Therefore, the coalescence kernel $\varGamma$ is equivalent to the collision kernel.

Figure 4(a) shows the contour of the normalised coalescence kernel ${\varGamma }/{ ( L_{\!f}u_{\textit{rms}} )}$ for ${\textit{Ca}} = 1.2$ . The sizes of the two coalescing aggregates are denoted as $D_{1}$ and $D_{2}$ . Results for other ${\textit{Ca}}$ are qualitatively similar. The coalescence frequency increases as the sizes of the coalescing aggregates grow. Moreover, the shape of the contours suggests that $\varGamma$ might simply depend on the sum of the sizes of two coalescing aggregates. Therefore, in figure 4(b) we plot ${\varGamma }/{ ( L_{\!f}u_{\textit{rms}} )}$ versus $D_{1} + D_{2}$ . The data collapse for all ${\textit{Ca}}$ cases. Moreover, beyond some experimental scatter, the trends are consistent with a power-law scaling $\varGamma \sim ( D_{1} + D_{2} )^{2}$ .

This trend may be modelled following the classic approach of gas-kinetic theory, extensively applied also in bubble–eddy collision and bubble coalescence problems (Luo & Svendsen 1996; Liao & Lucas 2010). In this picture (schematically depicted in figure 4 c), the coalescence kernel $\varGamma$ is estimated by the product of the collisional cross-section $L_c$ and the approaching velocity between both aggregates $\delta u$ , i.e. $\varGamma \sim L_c\delta u$ . Here, the maximum cross-section can be approximated by the sum of aggregate sizes ${L_{c} \approx D}_{1} + D_{2}$ . The approaching velocity is estimated via the second-order longitudinal structure function of the aggregate velocity $D_{LL}$ , obtained from tracking the centroids. This is shown in the Appendix for the representative case ${\textit{Ca}} = 2.7$ . For separations $r\lt L_{\!f}$ , a scaling $D_{LL}\sim r^2$ is observed, analogous to the power-law scaling of the fluid velocity in the same range (Boffetta & Ecke 2012; Shin et al. 2023). This is consistent with the finding of Shin & Coletti (2024) that the fluctuating kinetic energy of particles in this range of the parameter space is only marginally lower compared with the fluid. Therefore, although the aggregate inertia may not be negligible, we assume for the relative velocity a scaling $\delta u\sim r$ . As the centroid-to-centroid separation of coalescing aggregates is $r\sim D_{1} + D_{2}$ , we expect $\delta u\sim r\sim D_1+D_2$ . Substituting $L_c$ and $\delta u$ into $\varGamma$ leads to $\varGamma \sim ( D_{1} + D_{2} )^{2}$ .

Moreover, we consider the order of magnitude of $\varGamma$ . The aggregates considered in this work have sizes of the order of the forcing length scale $L_{\!f}$ (or somewhat smaller). This leads to the cross-section $L_c$ being of the order of $L_{\!f}$ as well. The approaching velocity $\delta u$ , although depending on the separation distance, is of the order of $u_{\textit{rms}}$ (Shin & Coletti 2024). This finally yields $\varGamma \sim \mathcal{O}(L_{\!f} u_{\textit{rms}})$ . Combining this equation with the aggregate size scaling as derived above, a model for the coalescence kernel is obtained:

(3.3) \begin{equation} \varGamma \sim L_{\!f}u_{\textit{rms}}\left ( D_{1} + D_{2} \right )^{2}d_{p}^{-2}. \end{equation}

Independently from the normalisation of the aggregate size by $d_{p}$ (which, as discussed above, is not demonstrably more appropriate than that by $L_{\!f}$ ), (3.3) is in remarkable agreement with the experimental data in figure 4(b).

We emphasise that this model can only hold for aggregates smaller than $L_{\!f}$ . For $D \gt L_{\!f}$ , (3.3) implies that the coalescence frequency increases indefinitely as the aggregate size grows. This scenario is not realistic as the formation of very large aggregates is likely to be interrupted by the forcing scale and/or by the physical boundaries of the system. Modelling the coalescence of very large aggregates thus requires further considerations outside the scope of the present study.

Figure 4. (a) The contour of the coalescence frequency $\varGamma$ for various coalescing aggregate sizes. (b) The coalescence frequency as a function of the sum of two coalescing aggregate sizes $D_{1} + D_{2}$ for various ${\textit{Ca}}$ . (c) A schematic to illustrate the coalescence process of two aggregates using the frame of reference of the lower aggregate for simplicity. The blue arrow illustrates the approaching velocity vector.

3.4. Number density of aggregates

Given the model forms of the breakup and coalescence terms, the evolution of the number density can be obtained numerically from the population balance equation. This is demonstrated for the example case ${\textit{Ca}} = 1.2$ , the data of which we use to determine the numerical pre-factors in (3.1), (3.2) and (3.3). We insert those in (1.1) and integrate in time via a first-order explicit Euler method, starting from an initially uniform number density $n ( {D}/{d_{p}} )d_{p}^{2} = 10^{- 3}$ . The aggregate size is discretised in bins of width $2d_{p}$ and the time step is $3.86 \times {10^{- 5}L_{\!f}}/{u_{\textit{rms}}}$ , though the outcome is not dependent on the choice of such parameters. The largest aggregate size allowed in the simulation is $50d_{p}$ which is about half of the size of the field of view.

Figure 5(a) shows the steady-state number density of the aggregates measured for ${\textit{Ca}} = 1.2$ as well as the evolution of the modelled number density. The number density $n(D,t)$ increases for small aggregates and decreases rapidly for large ones. The total number of particles $N_p$ reduces over time as illustrated in figure 5(b). Here, $N_{p} = d_{p}^{2}\int n ( {D}/{d_{p}} ) ( {D}/{d_{p}} )^{2}\text{d}({D}/{d_{p}})$ is obtained by multiplying the number of aggregates $n(D/d_p )\text{d}(D/d_p)$ by the number of particles in each aggregate $(D/d_p )^2$ and integrating the product over the considered size range. This continuous loss of particles is due to the coalescence kernel model (3.3) which increases indefinitely for large aggregates as previously discussed, generating aggregates of sizes $D \gt 50d_{p}$ which are discarded in the simulation. To eliminate this issue, coalescence frequency models appropriate for $D \gt L_{\!f}$ are needed, though the present experiment cannot inform such models due to the scarce number of large aggregates in the data.

Although a steady state is not achieved, the rate of particle loss decreases rapidly. This is shown up to ${\textit{tu}_{\textit{rms}}}/{L_{\!f}} = 0.035$ , for which ${{\rm d}N_{p}}/{{\rm d}t} = - 9.1u_{\textit{rms}}/L_{\!f}$ . This is much smaller than the characteristic particle loss rate $- 41u_{\textit{rms}}/L_{\!f}$ (the initial number of particles divided by the time scale $u_{\textit{rms}}/L_{\!f}$ ). Indeed, by that time, the number density of aggregates of size $D \leqslant 5d_{p}$ (accounting for 98 % of aggregates observed) evolves relatively slowly, and the modelled number density is in close agreement with the experiments. Considering the inherent limitations of the present models, the proposed scaling relations for the terms of the population balance equation capture the important dynamics of the system.

Figure 5. (a) The normalised aggregate number density obtained from the experiment (black symbols) and its evolution by solving the population balance equation (solid lines). (b) The temporal evolution of the normalised total number of aggregates per unit area.

4. Conclusion

We have examined the breakup and coalescence of particle aggregates floating at the interface of turbulent liquid layers, exploiting the framework of the population balance equation. While the latter has been at the basis of many studies in dispersed turbulent flows, typically experimental observation cannot directly validate specific parametrisations of its various terms. The two-dimensional nature of the present system allows us to detect virtually all coalescence and breakup events of aggregates, leading to direct measurements of the main elements: the breakup frequency, the daughter aggregate size distribution and the coalescence kernel.

In particular, the frequency of breakup is found to scale with the aggregate size as $\varOmega \sim D^{3/2}$ . This follows from assuming that aggregates are broken by particle-sized strain cells. The daughter aggregates follow a universal U-shaped distribution, which is recovered by assuming that the probability of detachment of a fragment is inversely proportional to its inertia. The coalescence kernel is found to depend on the sum of two aggregate sizes, $\varGamma \sim ( D_{1} + D_{2} )^{2}$ , which is successfully predicted by adopting a gas-kinetic model. These scaling relations apply to regimes dominated both by capillary-driven aggregation ( ${\textit{Ca}}\lt 1$ ) and by drag-driven breakup ( ${\textit{Ca}}\gt 1$ ). They allow us to close the population balance equation, which then is numerically integrated in time to simulate the evolution of the size-specific number density $n(D)$ . The agreement with the experiments confirms that the current models capture the important dynamics of the process.

We note that the framework accounts for the behaviour of aggregates smaller than the forcing length scale. Future work that focuses on larger ones is needed to obtain a complete picture of the dynamics of floating particle aggregates. Moreover, the current study could also be extended to aggregates at the free surface above three-dimensional turbulence, in which the divergence of the velocity field plays an important role in the transport of floating particles (Cressman et al. 2004; Li et al. 2024; Qi et al. 2025a , b ).