## 1. Introduction

The dynamics of gas bubbles in liquids drives a wide variety of operations in the chemical process industry, mineral processing and food industry, among many other examples. It also leads the mass exchange between ocean and atmosphere and the generation of aerosols (Deane & Stokes Reference Deane and Stokes2002). Breakup and coalescence of bubbles is commonly present in most applications, and equipment is designed based on the understanding of the gas–liquid interaction phenomena and on the bubble's behaviour. Thus, a profound knowledge of interaction between bubbles as well as between the gas and the liquid phases is crucial to understand the processes and optimize their performance. Bubble fragmentation has been extensively investigated under different flow configurations, and different models have been proposed (Tsouris & Tavlarides Reference Tsouris and Tavlarides1994; Martínez-Bazán, Montañés & Lasheras Reference Martínez-Bazán, Montañés and Lasheras1999*a*; Wang, Wang & Jin Reference Wang, Wang and Jin2003; Qi, Masuk & Ni Reference Qi, Masuk and Ni2020, among many others); however, there is still a large amount of scientific effort devoted to this topic. Coalescence is also essential to describe the dynamics and evolution of a population of bubbles. It is usually defined as a three-step process involving three different mechanisms (Prince & Blanch Reference Prince and Blanch1990; Chesters Reference Chesters1991). The first step consists of bringing close together the bubbles involved in the process, typically two of them. This step is controlled by the external liquid flow that induces the bubbles motion and causes them to collide. Once the bubbles are in contact, in order to coalesce, it is necessary to drain the thin liquid film separating them. The last stage initiates when the film becomes thin enough so that inter-molecular forces, such as van der Waals ones for pure fluids, become dominant, breaking the liquid film and thus making the bubbles coalesce. Considerable effort has been also devoted to better understand the underlying physics that characterizes the growth of the neck which forms just after a hole appears on the drained liquid film between coalescing bubbles or drops (Eggers, Lister & Stone Reference Eggers, Lister and Stone1999; Paulsen *et al.* Reference Paulsen, Carmigniani, Kannan, Burton and Nagel2014; Anthony *et al.* Reference Anthony, Kamat, Thete, Munro, Lister, Harris and Basaran2017; Moreno Soto *et al.* Reference Moreno Soto, Maddalena, Fraters, Van Der Meer and Lohse2018). It should be noted that the last stage of the drainage is very fast compared with the previous ones. Due to marked contrasts of coalescence efficiency observed when physico-chemical properties of the fluids vary, most of the coalescence studies concentrate on the drainage of the liquid film once the bubbles are sufficiently close (Marrucci Reference Marrucci1969; Chesters & Hofman Reference Chesters and Hofman1982; Oolman & Blanch Reference Oolman and Blanch1986; Zhang & Thoroddsen Reference Zhang and Thoroddsen2008; Ghosh Reference Ghosh2009; Huisman, Ern & Roig Reference Huisman, Ern and Roig2012). This analysis of coalescence as a three-step process has been fruitful to build several global models combining knowledge obtained from different studies. It remains that, in a given flow configuration, steps come one after another without real discontinuity, thus, the ratios of their respective lifetimes may vary depending on their ambiguous definition. In this sense, it can thus be interesting to analyse the coalescence process as a whole, using phenomenological models that avoid the complexity of describing in detail these stages. In the present work we focus on the global process, analysing the hydrodynamics controlling the bubble coalescence in a high Reynolds number confined bubble swarm. In order to explain precisely the aim of our study we first present the modelling formalism that we adopt and the closure laws that we discuss.

The evolution of sizes of a population of bubbles is often modelled by means of a Boltzmann-type partial integro-differential conservation equation (Williams Reference Williams1985),

where $n(v,\boldsymbol {x},t) \, \mathrm {d}v \, \mathrm {d}\kern0.06em \boldsymbol {x}$ is the probable number of bubbles with volume in the range $\mathrm {d}v$ about $v$ in the spatial range $\mathrm {d}\kern0.06em \boldsymbol {x}$ about $\boldsymbol {x}$ at time $t$, $\bar {\boldsymbol {u}}$ is the mean velocity of bubbles of volume $v$ at location $\boldsymbol {x}$ at time $t$, $\dot {Q}_{r}$ and $\dot {Q}_{d}$ are the birth and death rates of change of the number of bubbles due to breakup and coalescence, and ${\mathcal {R}}$ is the rate of change of the volume $v$ of a bubble which, for flows with no thermal effects, may be due to mass dissolution. In the present work we do not include thermal effects and dissolution can be neglected (${\mathcal {R}} = 0$) since the dissolution times are much larger than the characteristics residence time of bubbles for the bubble sizes considered. Equation (1.1) may depend on space and time if the problem is non-homogeneous and unsteady. A dependence on the velocities of the bubbles can also be introduced if, for a given size, possible velocities are distributed in a large range where coalescence and breakup dominant mechanisms may vary. For simpler presentation, we do not incorporate this effect in the following equations, as the flow regime that we consider does not require this supplementary complexity to be represented. When neglecting changes of volume due to a thermodynamical phase change, taking into account breakup as well as coalescence, this equation writes as (Coulaloglou & Tavlarides Reference Coulaloglou and Tavlarides1977; Martínez-Bazán Reference Martínez-Bazán1999; Marchisio & Fox Reference Marchisio and Fox2013)

where $\dot Q_c$ represents the sink or source terms of $n(v,\boldsymbol {x}, t)$ due to coalescence and $\dot Q_b$ due to breakup. Thus, the equation that determines the transport and the evolution of $n(v,\boldsymbol {x}, t)$ is the Liouville–Boltzmann's equation. It is a generalization of Smoluchowski's equation established for particle coagulation (Smoluchowski Reference Smoluchowski1917), usually called the population balance equation (PBE) (Williams Reference Williams1985). Moments of order 0 and 1 of $n(v,\boldsymbol {x},t)$ are respectively the total number of bubbles per unit volume, $N_{\infty }(\boldsymbol {x},t)$, whatever their sizes, and the volume fraction of the dispersed phase $\alpha (\boldsymbol {x},t)$ (Ramkrishna Reference Ramkrishna2000; Marchisio & Fox Reference Marchisio and Fox2013),

The coalescence and breakup rates in (1.2) read as

and

where $h (v, v')$ is the collision frequency between bubbles of volumes $v$ and $v'$; $\lambda (v, v')$ is the collision efficiency between bubbles of volumes $v$ and $v'$; $g_c(v)$ is the coalescence rate of bubbles of volume $v$ with any other bubble; $g_b(v)$ is the breakup or fragmentation frequency of bubbles of volume $v$; $m(v)$ is the number of bubbles resulting from the fragmentation of bubbles of volume $v$; and $f(v' , v)$ is the bubble size distribution of daughter bubbles resulting from the fragmentation of a mother bubble of volume $v'$. In (1.5) the first integral term on the right-hand side is a source term and the second one a sink term, both due to coalescence. Similarly, in (1.6) source and sink terms due to fragmentation also contribute to the evolution of the population of bubbles. As a complement to (1.5), the coalescence rate of bubbles of volume $v$ with any other bubble is defined by

Several closure laws and models for each of these terms have been proposed in the past. Their validity is most often limited to given hydrodynamical regimes of breakup and coalescence enforced by turbulent agitation or by mean shear flow. Sometimes they also include the influence of physico-chemical properties of the liquid or of the interface. A large amount of information on the adopted models can be found in references such as Coulaloglou & Tavlarides (Reference Coulaloglou and Tavlarides1977), Prince & Blanch (Reference Prince and Blanch1990), Chesters (Reference Chesters1991) and Martínez-Bazán *et al.* (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañés and Lasheras2010) or in literature reviews such as Kolev (Reference Kolev1993), Lasheras *et al.* (Reference Lasheras, Eastwood, Martínez-Bazán and Montañés2002), Liao & Lucas (Reference Liao and Lucas2009, Reference Liao and Lucas2010). However, the present work is focused on the coalescence processes of bubbles rising in liquid initially at rest, leaving breakage out of its scope. Thus, without including bubble breakup, (1.2) reduces to

In the literature there are two types of models that have been proposed for $h (v, v')$, or for the coalescence time, including phenomenological and physical models. The physical models are mainly focused on the description of the drainage process of the liquid film separating two bubbles when they get close enough (Chesters & Hofman Reference Chesters and Hofman1982; Oolman & Blanch Reference Oolman and Blanch1986; Ghosh Reference Ghosh2009; Huisman *et al.* Reference Huisman, Ern and Roig2012). They are thus limited to coalescence in stagnant liquids, unperturbed by the bubbles in motion. In contrast, the phenomenological models are introduced for moving bubbles and are based on models of collision of molecules applied in physical gas dynamics (Coulaloglou & Tavlarides Reference Coulaloglou and Tavlarides1977; Sovova Reference Sovova1981; Prince & Blanch Reference Prince and Blanch1990; Tsouris & Tavlarides Reference Tsouris and Tavlarides1994). In these models, the coalescence rate of bubbles of volume $v$ with any other bubble, defined by (1.7), is commonly reduced to the product between a collision frequency times a coalescence efficiency.

The objective of this work is to determine $\lambda (v, v')$ and $h (v, v')$ from bubble coalescence experiments performed in a high-Reynolds-number swarm of bubbles injected at the bottom of a planar vertical thin-gap cell filled with liquid at rest, and to analyse their contribution to $g_{c}(v)$. This confined configuration favours observation of coalescence. Starting from injection, the bubbles are greater than the gap thickness, favouring their interaction since the bubbles cannot escape out of the plane. Thus, coalescence is enhanced, being the coalescence rate larger than in the three-dimensional configuration (Lundin & McCready Reference Lundin and McCready2009).

In the regime explored here, there is no dewetting of the liquid films between the bubbles and the walls, and bubbles move at large Bond and Archimedes (or Reynolds) numbers. Thus, the cascade of sizes generated by coalescence is expected to create a complex self-induced gravity-driven agitation in the swarm. Indeed, for isolated bubbles, it is already known that bubbles whose sizes vary in a range similar to the one that we observe, exhibit contrasted oscillating paths and shapes (Kelley & Wu Reference Kelley and Wu1997; Roig *et al.* Reference Roig, Roudet, Risso and Billet2012; Filella, Ern & Roig Reference Filella, Ern and Roig2015; Piedra, Ramos & Herrera Reference Piedra, Ramos and Herrera2015; Hashida, Hayashi & Tomiyama Reference Hashida, Hayashi and Tomiyama2019; Pavlov *et al.* Reference Pavlov, D'Angelo, Cachile, Roig and Ern2021). In the present work we do not study the statistical properties of bubble agitation, but it has to be kept in mind that the self-induced agitation results from wake interactions, as already discussed in a homogeneous swarm of confined bubbles where coalescence was inhibited (Bouche *et al.* Reference Bouche, Roig, Risso and Billet2012, Reference Bouche, Roig, Risso and Billet2014).

It is worth pointing out that this flow configuration finds promising applications in chemical engineering since it is expected to be an alternative reactor of intermediate size that takes advantage of the confinement to enhance mass transfer, as in monolith reactors, and of the intense bubble-induced agitation to develop satisfactory in-plane mixing (Roudet *et al.* Reference Roudet, Billet, Cazin, Risso and Roig2017; Alméras *et al.* Reference Alméras, Risso, Roig, Plais and Augier2018). Some recent applications have been developed concerning light-activated reactions or cultivation of micro-algae in photo-reactors that need narrow geometries due to light absorption and attenuation, while keeping efficient mixing and mass transfer requirements (Oelgemoller Reference Oelgemoller2016; Pruvost *et al.* Reference Pruvost, Le Borgne, Artu and Legrand2017; Thobie *et al.* Reference Thobie, Gadoin, Blel, Pruvost and Gentric2017).

The paper is organized as follows. The experimental facility and the techniques developed to describe the evolution of the population of bubbles are presented in § 2. A careful examination of the performances of the bubble tracking algorithm is also presented in this section. The statistical properties of the gas flow and their evolution with the bubble population, for different values of the void fraction at injection, are reported in § 3. In particular, in § 3.1 we will introduce a parameter to properly characterize the evolution of bubble sizes. Then, the results of the rate of change of the population of bubbles and the measurements of the bubble collision frequency are summarized and discussed in § 3.2. The model associated with the collision frequency is introduced in § 4. Finally, § 5 is devoted to conclusions.

## 2. Experimental facility and techniques

The experimental facility used to characterize the evolution of the bubble size distribution in a high Reynolds number confined bubble swarm is presented in § 2.1. This particular flow configuration allowed a direct analysis of the whole bubble population using the shadowgraphy technique, since the planar motion prevented bubble overlap in the recording plane. In addition to the description of the operating conditions and of the shadowgraphy technique used, an overview of the image processing algorithm developed to detect, classify and track the bubbles in the swarm is given in Appendix A.

### 2.1. Experimental facility and operating conditions

The confined bubble swarm was generated within a quasi-bidimensional vertical cell filled with distilled water at ambient temperature, with the top section open to atmospheric pressure (Bouche *et al.* Reference Bouche, Roig, Risso and Billet2012, Reference Bouche, Roig, Risso and Billet2014). The cell consists of two parallel glass plates ($800$ mm high and $400$ mm wide) separated by a thin gap of width $w = 1$ mm (figure 1*a*). Unlike in previous works related to confined bubble swarms (Bouche *et al.* Reference Bouche, Roig, Risso and Billet2012, Reference Bouche, Cazin, Roig and Risso2013, Reference Bouche, Roig, Risso and Billet2014; Alméras *et al.* Reference Alméras, Cazin, Roig, Risso, Augier and Plais2016, Reference Alméras, Risso, Roig, Plais and Augier2018), no electrolyte was added to the liquid so that bubble coalescence was not inhibited. In addition, the distilled water was regularly renewed to prevent interface contamination. Air bubbles of uniform size were periodically injected from the bottom of the cell through an array of 16 capillary tubes of 0.6 mm inner diameter and 0.8 mm outer diameter (figure 1*a*). The tubes were equally distributed along the bottom of the cell and connected to a controlled pressure air feeding chamber. The pressure drop along the air injection tubes was sufficiently large to ensure a constant air flow rate along the tubes (Gordillo, Sevilla & Martínez-Bazán Reference Gordillo, Sevilla and Martínez-Bazán2007). The bubble detachment frequency, $f_b$, was accurately selected firstly setting a certain pressure in the air feeding chamber, $p_g$, and, secondly, controlling the air flow rate through each tube using individual valves. This frequency was checked for each injector using a stroboscopic light before running each experiment.

The volume of the injected bubbles ensured that their sizes were always larger than the width of the gap, $w$. Therefore, the bubbles were flattened between the cell walls, forming a thin liquid film between the bubble interface and the wall (see zoomed area in figure 1*a*). In this configuration, independently of the bubble size distribution of the swarm, no dewetting at the glass plates was observed and the bubbles degrees of freedom were bounded. The flow above the bubbles goes around sideways, and does not enter the thin liquid films at rest as in Pavlov *et al.* (Reference Pavlov, D'Angelo, Cachile, Roig and Ern2021). Therefore, these thin films are not expected to play any role in the bubble coalescence processes described in the present work. Hence, a two-dimensional description of the motion can be adopted as in Roig *et al.* (Reference Roig, Roudet, Risso and Billet2012) and Filella *et al.* (Reference Filella, Ern and Roig2015). In that sense, every bubble in the swarm is described by an equivalent diameter which is defined as $D=(4A_b/{\rm \pi} )^{1/2}$, with $A_b$ the projected area of the bubble on the cell plane. The injected gas volume fraction, $\alpha _0=\varSigma _{i} A^i_b w / (L_x L_z w)$, was determined from the total volume occupied by all the bubbles in a measuring window a few millimetres above the capillary tubes ($W1$ in figure 2), where $L_z=50.83$ mm and $L_x=328.67$ mm are respectively the height and the width of the window. In the current experiments, $\alpha _0$ was varied from $2.4\,\%$ to $6.7\,\%$ by adjusting the bubble generation frequency. Table 1 shows the experimental conditions of the four experimental sets considered in the present work, including the mean equivalent diameter of the injected bubbles, $D_0$. The variations in the sizes of the bubbles generated by each tube was negligible and a monodispersed bubble swarm was initially formed, as in Bouche *et al.* (Reference Bouche, Roig, Risso and Billet2012, Reference Bouche, Roig, Risso and Billet2014). Additionally, the total air flow rate was estimated as $Q_g = 4 {\rm \pi}{D_0}^2 w f_b$, showing a linear increase with $\alpha _0$.

The bubble swarm at the bottom of the cell is characterized by $\alpha _0$ and by the non-dimensional parameters governing the motion of an isolated bubble of equivalent diameter at injection, $D_0$. These include the Archimedes number, $Ar_0 = \sqrt {gD_0} D_0 / \nu$, the confinement ratio, $\delta _0 = w / D_0$, and the Bond number $Bo_0 = \rho g D_0^2 / \sigma$, where $g$ is the gravity, $\nu$ and $\rho$ the liquid kinematic viscosity and density, $w$ the thickness of the cell and $\sigma$ the surface tension. Under the conditions reported here, these parameters lie in the following ranges: $630 \leq Ar_0 \leq 850$, $0.24 \leq \delta _0 \leq 0.29$, and $1.79 \leq Bo_0 \leq 2.11$. For confined bubbles of equivalent diameter $D \geq D_0$, the mean rise velocity of an isolated bubble can be estimated by (Filella *et al.* Reference Filella, Ern and Roig2015)

which in dimensionless form can be expressed as,

where $Re= U_{b}D/\nu$ and $Ar = \sqrt {gD} D / \nu$ are the Reynolds and Archimedes number of a bubble of size $D$. The bubble Reynolds number, $Re_0 = U_{b}D_0/\nu$, can then be defined, and ranges here from $380$ to $500$. The gap Reynolds number, $Re_0{\delta _0}^2$, can then also be introduced to assess the inertial regime, as it varies between $28$ and $32$. Thus, the flow can be considered to be dominated by inertia (Bush & Eames Reference Bush and Eames1998) for all bubble sizes involved in the swarm. Bubbles at injection initially behave as isolated bubbles exhibiting oscillatory paths coupled to their unsteady wakes that generate periodic vortex shedding. Their in-plane projected shape is deformed and can be considered as an ellipse of moderate eccentricity (Roig *et al.* Reference Roig, Roudet, Risso and Billet2012; Filella *et al.* Reference Filella, Ern and Roig2015). For further discussion, general ideas can be retained. First, the velocity disturbances induced by bubbles in the liquid are mainly parallel to the plates, except in the close vicinity of the bubbles. Then, the order of magnitude of the liquid velocity in the bubble's wake is $\sqrt {gD}$. The wake is nevertheless strongly attenuated by shear stress at the walls within a characteristic time scale proportional to the viscous one, $\tau _{\nu } = w^2 / (4\nu )$. Therefore, in the swarm the agitation in the liquid results from direct interactions of localized random flow disturbances of various sizes as in the homogeneous swarm studied by Bouche *et al.* (Reference Bouche, Roig, Risso and Billet2014).

The swarm of bubbles generated was recorded using shadowgraphy in a measurement region that spanned almost the entire horizontal width of the cell (figure 1). To that aim, the cell was illuminated from behind with an uniform, constant and diffused white light perpendicular to the cell plane (figure 1*a*). Placing the light source at one side of the cell, a camera (Photron APX) equipped with a $85$ mm lens was used to take images of $2^{10}$ levels of grey and of size $1024 \times 512$ pixels, with an exposure time of $1/2000$ s, from the other side. In order to analyse the evolution of the population of bubbles as they rise, while maintaining the desired resolution, the backlight and the camera were placed at three different positions (figure 2). Transverse homogeneity of the flow was observed. Therefore, the downstream evolution of the bubble swarm was described by a statistical analysis of the bubble population characteristic parameters averaged over the horizontal width of the cell. In order to avoid possible errors due to border effects, the analysis was performed in a recording region, placed in the middle of each image, which consisted of a rectangle of size $328.67\,{\rm mm} \times 152.49\,{\rm mm}$ with a pixel-size resolution of $350\,\mathrm {\mu }{\rm m}$ (figure 1*b*). Moreover, to increase the spatial resolution of the measurements, the recording region was divided into five windows of horizontal length $L_x = 328.67$ mm and vertical length $L_z = 50.83$ mm, with $50\,\%$ of overlapping, as indicated in figure 2. Two types of measurements were performed, depending on the image acquisition frequency. First, a series of experiments taking video images of the swarm at 250 f.p.s. were conducted. This recording rate was high enough to track the bubbles as they rose along the field of view. Thus, bubble collisions could be detected and tracked to determine whether the colliding bubbles coalesced, generating new larger bubbles, or eventually separated, which allowed us to determine the bubble collision frequency, $h$, as well as the efficiency, $\lambda$. Measurements also allowed us to detect breakup events, giving birth to smaller daughter bubbles, although this phenomenon was rare in this study. For this analysis, high-speed movies of 25 s duration were recorded at the three positions. The total duration of the recordings included between $20\,000$ and $75\,000$ bubble records per position, depending on the injection conditions and on the measuring location. In addition to the experiments recorded at high frequency, at each recording position, sets of around 3000 uncorrelated images were recorded at a frame rate of 1/2 f.p.s. to ensure that the bubbles in one image were not recorded in the following one. Thus, the total number of bubbles detected at each position varied between $25\,000$ and $250\,000$, depending on the injected air flow rate and on the recording position. This ensured a statistically converged and unbiased measurement of parameters of the bubble population that can be obtained from low frequency experiments such as the volume probability density function. Satisfactory comparison of the information that could be obtained from both types of measurements indicated that the statistical parameters extracted only from the high-frequency records were indeed robust and meaningful.

Details of the digital image processing methods specifically developed in this work to detect and classify the bubbles in the swarm are given in § A.1. In addition, the techniques designed to track the bubbles and detect the collisions, as well as the coalescence and breakup events are described in § A.2. Additional information can be found in Ruiz-Rus (Reference Ruiz-Rus2019).

### 2.2. Performance of the bubble tracking algorithm (BTA)

The results obtained with the bubble tracking algorithm (BTA) described in Appendix A, consist of the record of the bubbles along the field of view for each position. The information stored for each bubble includes the projected area (bubble volume), the centroid location, the bubble velocity components, as well as the bubble lifespan and the types of birth and death events. In addition, family trees are established for each newly generated bubble, including the parents in a birth from coalescence, or the mother and the sibling in a birth from breakage. The performance of the BTA algorithm can be estimated, first, from the fact that more than $99 \, \%$ of the detected bubbles can be tracked, the remaining ones being associated to specific events with simultaneous coalescence and breakup in agglomerates of bubbles.

As an example, figure 3 shows a set of bubble trajectories for each injection condition at the first recording position. In this figure, regardless of where the bubble trajectories begin, they have been displaced to the same origin, with $x_o$ and $z_o$ being the initial positions of the bubbles. It can be observed that the horizontal dispersion of the bubbles increases with the injected air volume fraction and with the vertical position, showing the effects of the liquid velocities induced by the wakes of the population of bubbles. The bubble lifespan is typically larger for the lowest values of $\alpha _0$ (figure 3*a*,*b*), since the number of coalescence events is still low at this recording position. The trajectories show the characteristic path oscillations described by Roig *et al.* (Reference Roig, Roudet, Risso and Billet2012) for isolated bubbles, indicating the weak effect of the hydrodynamic interactions at these low void fractions. However, the degree of coalescence substantially increases with $\alpha _0$, resulting in much shorter trajectories (figure 3*c*,*d*).

The first quantitative measurement extracted from the BTA is the collision efficiency. It represents the ratio between the number of coalescence and that of collisions. As mentioned before, the confinement of the bubbles imposed in this configuration highly increases the efficiency of the collisions. The mean collision efficiency, $\lambda _\infty$, obtained considering all the collisions detected at each measuring window, is plotted in figure 4 for the four values of $\alpha _0$ tested in this work. Our results indicate that the collision efficiency barely depends on the size of the colliding bubbles. In fact, considering the collision efficiency of different pairs of bubbles, differences lower than $5 \, \%$ were found with respect to $\lambda _\infty$. In addition to being a very high efficiency, it should be noted that $\lambda _\infty$ does not vary with the concentration of bubbles, nor with $z$. In fact, the figure shows that its value can be considered constant and approximately equal to $80 \, \%$.

Moreover, the tracking method allows a direct analysis of the coalescence events. In that sense, any detected collision is individually tracked (see § A.2), obtaining the position at which it initially occurs as well as the corresponding information for the colliding bubbles. In addition, the BTA can be used to obtain the number of bubbles of volume $v$ that die due to coalescence, and the mean coalescence frequency of all bubbles at each position,

that represents the frequency at which a bubble of any size coalesces with other bubbles. In the experiments performed in this work, most of the daughter bubbles which are born due to coalescence come from previous binary collisions. Thus, since the total number of collisions which end up in coalescence, $\gamma _\infty$, during the measuring time, $T$, in a population of $N_\infty$ bubbles, is half of the bubbles dying due to coalescence, $N_\infty {\langle g_c \rangle }_\infty /2$, the mean bubble coalescence frequency can be directly obtained from the experimental measurements as

Similar to the mean coalescence frequency, ${\langle g_c \rangle }_\infty$, the mean breakup frequency of bubbles at each position, given by ${\langle g_b \rangle }_\infty = \int ^{\infty }_{0} n(v) g_b(v)\,{\rm d}v/ \int ^{\infty }_{0} n(v)\,{\rm d}v$, can be directly obtained as

where $\psi _\infty$ is the number of breakup events observed during the time $T$.

Thus, the accuracy and convergence of the tracking analysis of the experiments performed at high rates of acquisition can be checked by comparing both sides of (1.2) averaged over all bubble sizes, assuming that both coalescence and breakup are binary processes. To that aim, the right-hand side of the averaged PBE can be achieved by integration of (1.2) over the whole range of bubbles (Friedlander Reference Friedlander1977). Thus, applying the Leibnitz rule for integration and substituting (1.5) and (1.6) into (1.2), in the one-dimensional, steady state situation of interest here, the averaged PBE simplifies to (Kocamustafaogullari & Ishii Reference Kocamustafaogullari and Ishii1995; Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019)

where $\bar {U}_z$ is the mean rising velocity of the bubbles in the measuring window. Figure 5 shows the different terms on both sides of (2.6), evaluated at various measuring locations, for all the experimental injection conditions. In that case, the frequency terms have been made dimensionless making use of $\sqrt {g/D_0}$, while the downstream locations are normalized with the injection diameter $D_0$, corresponding to each experimental injection condition. As expected, a fairly good agreement is observed between both sides of the equation (black circles and grey diamonds, respectively). This result confirms the validity of the experimental procedure, as well as the effectiveness of the bubble tracking method and the convergence of the results obtained. In addition, it can be noticed that some breakup events (hollow squares) also take place in the swarm, especially for the higher values of $\alpha _0$ (figure 5*c*,*d*).

## 3. Description and discussion of bubble coalescence in the evolving swarm

An overview of the results obtained by the BTA has shown that the evolution of the population of bubbles is mainly governed by bubble coalescence, with a weak contribution of bubble breakup in some cases (figure 5). These processes, that lead to variations in the bubble size distribution, are driven by the liquid agitation in the swarm, which in turn is induced by the interaction of the wakes of bubbles of different sizes that constitute the swarm. Consequently, it is necessary to characterize the different stages of the evolution of the bubble population to adequately elucidate the mechanisms that govern the coalescence process.

### 3.1. Spatial evolution of the bubble population

In the present configuration there is no external liquid flow that carries the bubbles and they rise due to buoyancy effects. The downstream evolution of the population of bubbles can be described in terms of the flux of bubbles crossing each $z$ position. An estimation of the averaged flux is given by the local net number of bubbles detected in each measuring window, $N^*_\infty$, multiplied by their corresponding mean velocity, $\bar {U}_z$, obtained by the BTA. Figure 6(*a*) shows the downstream evolution of the bubbles flux for each injection condition. The fact that the flux of bubbles decreases with the dimensionless downstream location $z/D_0$, indicates that coalescence leads the evolution of the distribution of bubbles, even in the cases where there are some breakup events, as previously shown in figure 5. The rate of change of the population depends on the initial number of bubbles, $N^*_\infty (0)$, which is directly related to the selected bubble generation frequency, $f_b$, and thus, to the injected air volume fraction, $\alpha _0$. For high void fractions (i.e. $\alpha _0 = 4.9$ and $6.7\,\%$, respectively), near the bottom of the cell, the amount of bubbles of injection size $D_0$ quickly decreases as the bubbles rise. In fact, strong bubble–bubble interactions occur in the regions close to the injectors, giving rise to collisions and coalescence events. Once the bubbles start coalescing, larger bubbles are generated leading to a coalescence cascade which rapidly involve pairs of bubbles of wider ranges of sizes. As $z/D_0$ increases, the flux decreases less rapidly. In fact, the reduction in the total number of bubbles composing the population causes the net amount of coalescence events to decrease too. Although the rate of change of the bubbles flux decreases with the downstream distance, there is no evidence of reaching a final frozen state where coalescence would become negligible. It has to be pointed out that at higher positions of the cell, as previously noticed in figure 5(*c*,*d*), bubble breakup occurs, competing with coalescence. The unstable nature of the largest bubbles generated in these cases as well as the interaction with stronger bubble-induced liquid velocities, increases the relevance of breakage far from the injection point. On the other hand, for lower void fractions (i.e. $\alpha _0 = 2.4\,\%$ and $3.2\,\%$), figure 6(*a*) shows that the bubbles take longer to begin to coalesce. In these cases, the bubble flux initially remains constant up to a certain position where it starts to decay at a rate that decreases as $\alpha _0$ decreases. This is related with the lower amount of bubbles generated under these conditions and the larger distances between bubbles at the initial positions. However, as they rise, the injected bubbles loose memory of the injection conditions and begin to adopt the oscillatory motion which characterizes these ellipsoidal bubbles within the confined cell (Roig *et al.* Reference Roig, Roudet, Risso and Billet2012; Filella *et al.* Reference Filella, Ern and Roig2015). At a given height that depends on the generation frequency (Sanada *et al.* Reference Sanada, Watanabe, Fukano and Kariyasaki2005), the trajectories of the bubbles scatter under the effect of the perturbations in the liquid, giving rise to bubble–bubble interactions and to the subsequent collisions and coalescence events (see $W4$ and $W5$ in figure 2).

Figure 6(*b*) shows, for each injection condition, the downstream evolution of the local air volume fraction $\alpha (z)$, defined as the volume occupied by the entire population of bubbles present at each measuring window, divided by the volume of the window. At variance with homogeneous monodispersed bubble swarms where the gas volume fraction remains constant under constant injection conditions (Martinez Mercado *et al.* Reference Martinez Mercado, Chehata, Van Gils, Sun and Lohse2010; Bouche *et al.* Reference Bouche, Roig, Risso and Billet2012; Colombet *et al.* Reference Colombet, Legendre, Risso, Cockx and Guiraud2015), in the present swarm, $\alpha (z)$ also varies with $z/D_0$, due to the evolution of the distribution of sizes and to the fact that the velocity of the bubbles varies with their sizes. Indeed, the absence of an external liquid flow implies that the mean rising velocity of the gas phase is mainly imposed by the buoyancy exerted on the different bubble sizes coupled to the underlying fluid motion generated by hydrodynamic interactions. Therefore, the local volume fraction $\alpha (z)$ is affected not only by the injected air flow rate, but also by the evolution of the distribution of bubbles that induces buoyancy-driven variations on the mean rising velocity of the gas phase.

In order to determine the distribution of bubble sizes at each measuring window, the equivalent diameters of the bubbles were obtained from image processing (see § A.1) to compute the bubble volume probability density function (v.p.d.f.) (Martínez-Bazán *et al.* Reference Martínez-Bazán, Montañés and Lasheras1999*a*),

which represents the volume occupied by bubbles of size $D$ compared with that of the entire distribution. In (3.1), $D_{min}$ is the smallest bubble size of the distribution and $D_{max}$ the largest one. The downstream evolution of the v.p.d.f. resulting from coalescence, and eventually breakup, is shown in figure 7 for the four experimental conditions reported in table 1. For the sake of clarity, we have only plotted six measuring locations. Qualitatively, similar downstream evolutions are observed for the four cases: the nearly monodispersed distribution of bubbles observed close to the bottom of the cell progressively widens further downstream due to bubble coalescence. Since bubbles of constant size, $D_0$, are periodically injected at the bottom of the cell, for low values of $\alpha _0$, the initial v.p.d.f. is a narrow distribution around $D_0$ (see the distribution at $z = 40.92$ mm in figure 7*a*,*b*). In fact, in these cases, coalescence is not observed until $z = 142.58$ mm (figure 7*a*,*b*), as previously noted from the evolution of the total flux of bubbles shown in figure 6(*a*). However, for larger values of $\alpha _0$, at the first measuring window, a secondary peak is already observed at $D \simeq \sqrt {2}D_0$, indicating the existence of some coalescence events of bubbles of size $D_0$ ($z = 40.92$ mm in figure 7*c*,*d*). It is worth noting the existence of additional peaks at $\sqrt {3}D_0, \sqrt {4}D_0, \ldots$, corresponding to added volumes of the injection bubbles resulting from successive coalescence events (Néel & Deike Reference Néel and Deike2021). Nevertheless, as the coalescence process evolves, such peaks (associated with classes of finite extension) attenuate, generating broader and smoother distributions far from the injection point. The v.p.d.f.s exhibit large tails as the coalescence cascade progresses. In fact, it can be observed that the size of the largest bubbles found at a certain distance increases with $\alpha _0$. This fact is clearly illustrated in the images displayed as insets in figure 7, which show characteristic snapshots of the swarm for each experimental condition at $z = 412.92$ mm.

For all values of $\alpha _0$, despite the differences in the downstream evolution of the rate of change of the number of bubbles, similar shapes of the distribution are found at different positions, which can be understood as equivalent stages of the coalescence cascade process. For example, the distribution at $z=324.58$ mm in figure 7(*a*) is similar to that at $z=222.92$ mm in figure 7(*b*), the distribution at $z=324.58$ mm in figure 7(*b*) matches that at $z=142.58$ mm in figure 7(*c*) and the distribution at $z=514.58$ mm in figure 7(*c*) resembles that at $z=412.92$ mm in figure 7(*d*). The similarity in the distribution evolutions reveals that the various swarms considered follow the same coalescence cascade, independently of the value of $\alpha _0$. However, the rate of change of the swarm strongly depends on the void fraction, which is directly related to the amount of bubbles forming the population, as can be inferred from (1.5). It will take longer for low values of $\alpha _0$ (figure 7*a*,*b*) than for larger ones (figure 7*c*,*d*) to reach a given stage of the distribution of sizes (i.e. a given shape), although each size will undergo the same coalescence cascade independently of $\alpha _0$. Taking this into account, the evolution of the bubble size distribution is characterized by a diameter representative of the population of bubbles. For this purpose, we used the statistically robust parameter, $D_{V90}$ (Hinze Reference Hinze1955; Martínez-Bazán, Montañés & Lasheras Reference Martínez-Bazán, Montañés and Lasheras1999*b*), defined as the diameter of a bubble such that 90 % of the total volume of the population of bubbles is contained within bubbles smaller than $D_{V90}$. This characteristic diameter represents the size of the largest bubbles in the distribution which, as they rise, will induce the largest velocity fluctuations in the liquid.

Figure 8(*a*) shows the downstream evolution of $D_{V90}$ for the four injection conditions, reflecting the increase of its rate of change with $\alpha _0$. Note that, for $\alpha _0= 2.4\,\%$ and $3.2\,\%$, $D_{V90}$ barely changes for $z< 150$ mm, indicating that coalescence does not start in those cases until $z > 150$ mm. However, for $\alpha _0= 4.9\,\%$ and $6.7\,\%$, coalescence is already observed near the injection position. In addition, figure 8(*b*) shows the evolution of the flux of bubbles normalized by that of the first measuring window, $N^*_{\infty } \bar {U}_z(0)$, as a function of the characteristic diameter $D_{V90}$ normalized by the diameter of the injected bubbles, $D_0$. It can be observed that the four plots collapse on the same curve, corroborating that the population of bubbles follow similar coalescence cascade processes and that $D_{V90}$ is the proper variable to describe the evolution of the population of bubbles of the different swarms. In that sense, $D_{V90}/D_0$ can be seen as a variable similar to the dimensionless time, $t/\tau$, proposed by Smoluchowski (Reference Smoluchowski1917) to describe the Brownian coagulation of particles within nearly monodispersed systems (see Chandrasekhar Reference Chandrasekhar1943; Friedlander Reference Friedlander1977, for reports of this work in English), being $\tau$ the characteristic coagulation or decay time, usually called half-life of the population. So, in the following, $D_{V90}$ will be considered as the characteristic parameter to determine the evolution of the swarm. This will allow us to consider the coalescence frequency as an unknown that varies with the self-induced evolution of the population, which includes the influence of the *a priori* unknown velocities of each size, as well as the $\alpha _0$ dependence.

As already indicated above, the aim of the present work is to experimentally determine the coalescence frequency of a pair of bubbles of sizes $D$ and $D'$, respectively, in the bubble swarm. To be able to do this, large enough size ranges have to be defined to ensure that the number of bubbles included in each range is sufficiently high to have an adequate number of colliding bubbles and, thus, obtain statistically converged results. For this purpose, we discretized the population of bubbles obtained from the experiments in different size classes, denoted as class 0, 1, 2 and so on. Every class is represented by a diameter $D_k$ and includes bubble sizes in the range $D_k -\varDelta _k/2$ $\le D \le$ $D_k +\varDelta _k/2$, within a size bin of width $\varDelta _k$. The diameter $D_k$ represents the middle size of the bin and corresponds to integer values of the injection bubble volume, accounting for volume-conservative coalescence events. Thus, the initial class corresponds to the injection bubbles, $D_0$. The following class is associated to the coalescence of two bubbles of size $D_0$, being $D_1^2 = 2D_0^2$. The rest of the classes, $D_k$, are defined as the diameter of the bubble formed from the coalescence of two bubbles belonging to the two preceding classes, $D_k^2 = D^2_{k-1} + D^2_{k-2}$. Therefore, the different classes represent added volumes of the injection bubbles resulting in the following diameters: $D_1=\sqrt {2} D_0$, $D_2=\sqrt {3} D_0$, $D_3=\sqrt {5} D_0$, $D_4=\sqrt {8} D_0$, $D_5=\sqrt {13} D_0$, $D_6=\sqrt {21} D_0$, $D_7=\sqrt {34} D_0$ and $D_8=\sqrt {55} D_0$. The different sizes of the bins, $\varDelta _k$, were chosen looking for substantial but smooth variations of the bubble characteristics. Such definition of bubble classes and their corresponding limits allowed us to have an amount of bubbles sufficiently large in each class to observe enough collision events among bubbles of different classes, ensuring a good statistical convergence of the data. Details of the bubble classes defined for the different injection conditions are listed in table 2.

Considering the bubble classes defined in table 2, the number of bubbles of a certain class per unit volume, $N_k$, is obtained by integrating (1.3) between the size limits of the corresponding bin, $D_k -\varDelta _k/2$ and $D_k +\varDelta _k/2$. The evolution of the fraction of bubbles of each class, $N_k / N_{\infty }$, is represented in figure 9, as a function of $D_{V90} / D_0$ for $\alpha _0= 6.7\,\%$, showing the contribution of the amount of bubbles of each class to the whole population in every stage of the coalescence process. Note that, as shown in the inset of figure 9, the fraction of injection bubbles (class $k=0$) is always larger than those of the other classes, indicating that there is still a considerable amount of bubbles of size $D_0$ remaining even far from the injection point. This fact reveals the appreciable presence of these small bubbles even at stages in which the coalescence cascade has accounted nearly $8^2$ events. After a rapid decrease during the early stages of the evolution, the fraction of injection bubbles almost stabilizes around half of the total number of bubbles. However, the evolution of the fraction of all the other classes remarkably differs from that of the initial bubbles. In fact, a certain class of bubbles does not begin to form until a pair of smaller bubbles, whose sum of volumes is equal to the volume of the forming bubble, coalesce. First, the number of bubbles of classes $k>0$ increases due to coalescence of smaller bubbles. Afterwards, when the number of bubbles of a given class becomes sufficiently large, they start to coalesce forming larger bubbles. Eventually, when the number of coalescing bubbles of a class $k$ is larger than the number of bubbles that are generated from the coalescence of smaller ones, $N_k/N_\infty$ decreases. Indeed, the evolution of the concentration of a certain bubble class within the swarm is driven by the balance between the coalescence rate of smaller bubbles that form bubbles of this class, and their rate of coalescence with all the others, as is clearly established by (1.5). The coalescence rate in the balance equation (1.8) is determined from the collision frequencies of pairs of bubbles, $h(D_k,D_{k'})$, whose experimental values will be obtained in § 3.2 and modelled in § 4.

### 3.2. Determination and analysis of the bubble pair collision frequency, $h(D_k, D_{k'}$)

Considering the discretization in bubble classes given in table 2, the number of bubbles per unit volume of a given class $k$ that die per unit time due to coalescence represents the coalescence death rate, commonly expressed as

In (3.2), $h(D_k, D_{k'})$ is the collision frequency of bubbles of class $D_k$ with bubbles of size $D_{k'}$ and has units of ${\rm m}^3 \,{\rm s}^{-1}$, and $n(D_k)$ is the number density of bubbles with units of ${\rm m}^{-3}$, thus, $\dot D_{ec}(D_k)$ has units of ${\rm m}^{-3}\,{\rm s}^{-1}$. Note that, since $\lambda (D_k, D_{k'})=\lambda _\infty$ is constant in the present case, $h(D_k, D_{k'})$ can be treated indistinctly as the collision or the coalescence frequency of the pair of bubbles. Having this in mind, $h(D_k, D_{k'})$ was obtained from the experiments performed at high acquisition rates, applying the bubble tracking and coalescence detection algorithm described in Appendix A, as

In (3.3), $\varGamma _{kk'}$ is the number of bubbles of class $k$ colliding with bubbles of class $k'$ in the measuring window, of volume $L_z \times L_x \times w$ (see figure 2), per unit time, and $N_k$ and $N_{k'}$ are respectively the number of bubbles of class $k$ and $k'$ accounted in the volume of the measuring window. Thus, the rate of loss of bubbles of class $k$, $\dot D_{ec}(D_k)$, corresponds to the frequency at which effective collisions of bubbles of size $D_k$ with the rest of the bubbles take place, which, assuming that $\lambda (D_k, D_{k'})=\lambda _\infty$ is constant (see figure 4), reduces to ${\dot D_{ec}(D_k)=\lambda _{\infty }}\varSigma {_{k'=0}^{\infty } \varGamma _{kk'}}$. Note that, considering all the bubbles of the distribution which collide per unit time, $\varGamma _\infty =\varSigma _{k=0}^{\infty }\varSigma _{k'=0}^{\infty } \varGamma _{kk'}$, the total number of coalescence events in expression (2.4), can also be obtained as $\gamma _\infty / T = \lambda _\infty \varGamma _\infty /2$. It is worth indicating that the models for $h$ have been commonly derived by making an analogy with the kinetic theory of gases (Vincenti & Kruger Reference Vincenti and Kruger1965). These models (see, e.g. the review in Liao & Lucas Reference Liao and Lucas2010) generally consider $h$ as the volume swept per unit time by the colliding bubbles. It is usually referred to as the collision kernel, or the coalescence kernel if the efficiency is also included (Marchisio & Fox Reference Marchisio and Fox2013). However, taking into account that $h$ is symmetric, resulting in $h(D_k, D_{k'}) = h(D_{k'}, D_k)$, in the present work it will be referred to as the bubble pair collision frequency.

The evolution of $h(D_k, D_{k'})$ is described hereafter for the experimental case corresponding to $\alpha _0 = 4.9\,\%$, as a representative example. Figure 10 shows the results of the bubble pair collision frequency for various bubble pairs at a stage of the evolution of the swarm where $D_{V90} = 18.25$ mm ($D_{V90}/D_0 = 4.74$). The adopted representation displays the evolution of $h(D_k,D_{k'})$ with $D_{k'}$ for the different classes of bubbles (constant values of $D_k$). We should remember that, according to figure 8, since the bubble size distributions evolve in a similar way for all values of $\alpha _0$, $D_{V90}/D_0$ can be used as a parameter to describe the coalescence cascade process. Considering the results corresponding to the smallest bubble size present in the distribution, $D_k = 3.85$ mm ($\CIRCLE$), a monotonous increment of the collision frequency is observed when the size of the other colliding bubble, $D_{k'}$, increases. This behaviour is also initially observed for larger bubbles, i.e. larger values of $D_k$, up to a certain value of $D_{k'}$ beyond which the frequency begins to decrease. Taking into account the interchangeability of $D_k$ and $D_{k'}$, it can be stated that the slope of the curve $h(D_k, D_{k'})-D_{k'}$ increases with $D_k$. This monotonically increasing behaviour of the bubble pair collision frequency, characterized by the fact that at least one of the bubbles involved in the collision is relatively small, will be referred to as the first regime (represented by solid symbols in figure 10). In this regime, as soon as the bubbles collide, they coalesce or (in a very few cases) bounce back, but their surfaces do not deform during a certain time until they coalesce. On the other hand, as observed for $D_k \geq 8.62$ mm ($\blacksquare, \blacktriangledown, \blacklozenge$), the evolution of the bubble pair collision frequency shows a local maximum at certain values of $D_{k'}$, which decreases as $D_k$ increases. After this maximum a different behaviour appears, which defines a second collision/coalescence regime (represented by hollow symbols in figure 10), where the two involved bubbles are large enough to be able to deform during the coalescence process. In this regime, once the bubbles get in touch, their interfaces deform and flatten forming a thin liquid film between the two bubbles, which mainly drains sideways (Huisman *et al.* Reference Huisman, Ern and Roig2012; Pavlov *et al.* Reference Pavlov, D'Angelo, Cachile, Roig and Ern2021). The time spent on the bubbles surface deformation and on the liquid film drainage, before coalescence, increases the time during which the bubbles interact, leading to a reduction of the corresponding collision/coalescence frequency.

Figure 11(*a*–*d*) shows the contours of $h(D_k, D_{k'})$ at four different instants of the coalescence cascade process corresponding to $\alpha _0=$ 4.9 %. In this case, bubbles of diameter $D_0 = 3.85$ mm are initially injected and, as they coalesce while they rise, a population of bubbles of increasing diameters is generated. The values of $D_{V90}$ of the bubble size distributions corresponding to panels (*a*–*d*), indicated with horizontal and vertical dashed lines, are 10.93 mm, 13.96 mm, 18.25 mm and 20.79 mm, respectively ($D_{V90}/D_0 = 2.84$; $3.62$; $4.74$ and $5.40$). The plots clearly show that $h(D_k, D_{k'})$ is a symmetric function that increases with $D_{V90}$. In fact, it can be observed that the coalescence frequencies between two small or two large bubbles are small compared with the collision frequency between bubbles of intermediate sizes. Note, for instance, that the maximum coalescence frequency in figure 11(*d*) falls in a fringe corresponding to the coalescence of bubbles of $D_k=8$ mm ($D_k=20$ mm) with bubbles of $D_{k'}= 20$ mm ($D_{k'}= 8$ mm), or to the coalescence between two bubbles of diameter $D_k=D_{k'} \approx 12$ mm. Interestingly, a close inspection of the contour plots indicates that the colour levels nearly follow lines of constant values of

drawn with solid lines in figure 11(*a*–*d*). The diameter $D_p$ represents the diameter of a bubble whose radius of curvature is equal to the mean radius of curvature of the pair of interacting bubbles, and will be further called the reduced diameter. This diameter has been commonly used to describe the deformation of interacting bubbles and in models of drops/bubbles coalescence (Chesters & Hofman Reference Chesters and Hofman1982; Kamp *et al.* Reference Kamp, Chesters, Colin and Fabre2001; Neitzel & Dell'Aversana Reference Neitzel and Dell'Aversana2002). In addition to $h(D_k, D_{k'})$ shown in figure 11(*a*–*d*), figure 11(*e*–*h*) displays the bubble pair collision frequency weighted by the number of bubbles of size $D_{k'}$. This quantity represents the frequency of collision of a bubble of size $D_k$ with bubbles of size $D_{k'}$ and obviously depends on the number of bubbles of size $D_{k'}$ in the distribution, $n(D_{k'})$. It represents the contribution of bubbles of class $k'$ to the coalescence frequency of bubbles of class $k$, $g_c(D_k)$, as established in (1.7). It can be seen that $h(D_k, D_{k'}) N_{k'}$, obtained experimentally as $\varGamma _{kk'}/ N_k$ is no longer symmetric and it depends on the bubbles size distribution, having the maximum values for $D_{k'}=D_0$ in our particular case, since the number of bubbles of size $D_0$ is larger than the number of bubbles of other sizes, as shown in figure 9.

A description of the bubble pair collision frequency, $h(D_k,D_{k'})$, is therefore sought after in terms of the reduced diameter, $D_p$. The experimental values of $h(D_k,D_{k'})$ displayed in figure 10 are represented as a function of the corresponding $D_p$ in figure 12(*a*). They fall on a single curve, represented by the thick solid line, which corresponds to the averaged bubble pair collision frequency along lines of constant $D_p$ in figure 11(*c*). This curve clearly indicates the existence of the two different collision regimes commented above, properly distinguished now as a function of $D_p$. Initially, in the first regime (solid symbols), $h(D_p)$ increases with $D_p$ until it reaches a maximum value beyond which it begins to decrease with $D_p$ (hollow symbols). In this figure, as in figure 10, the cases indicated by arrows correspond to the time series of experimental images shown in figure 13. Equivalent results are found for different values of $D_{V90}$, as shown in figure 12(*a*) (different lines) for the instants of the coalescence cascade presented in figure 11(*a*–*d*). This result corroborates that $D_p$ properly captures the dependence of the collision frequency on the bubble sizes. The values of the reduced diameter at which the maximum of $h(D_p)$ occurs, i.e. the change from the first to the second collision regime, denoted as $\tilde {D}_p$, are displayed in figure 12(*b*) for the different experimental cases tested. For the smallest value of the injected volume fraction, $\alpha _0 = 2.4\,\%$, the entire coalescence cascade took place in the first regime, without transitioning to the second one, and no data are represented in this case. It can be observed that $\tilde {D}_p$ increases with the concentration of bubbles and with $D_{V90}$, following a power law given by $\tilde {D}_p/w\propto (\alpha D_{V90}/w)^{1/2}$ as it will be commented later on in § 4.

In order to better illustrate the main characteristics of the two regimes mentioned above, different coalescence events are shown in the time series of snapshots displayed in figure 13. The time intervals between images are the same for the four series. The black dots inside the bubbles denote the instantaneous position of their centroids, while the coloured ones indicate their previous positions, describing the bubbles trajectories. The cases shown are representative of the collisions taking place during the evolution of the swarm. The first three series (figure 13*a*–*c*) have been selected for the same injection condition and stage of the swarm ($\alpha _0 = 4.9\,\%$ and $D_{V90}/D_0 = 4.74$), and correspond to the points indicated with arrows in figures 10 and 12. The processes involve pairs of bubbles where the diameter of one of them is $D_k=13.90$ mm ($\blacktriangledown$), whose trajectory is represented with red dots in figure 13. The other bubbles that form the pairs have different diameters $D_{k'}$, as indicated in figure 10, and their trajectories are represented by blue dots in the corresponding panels of figure 13. On the other hand, figure 13(*d*) represents a coalescence event in a swarm also at $D_{V90}/D_0 \approx 4.74$, where the pair of bubbles are similar to those in figure 13(*c*), but at a higher void fraction, $\alpha _0 = 6.7\,\%$. As specifically indicated in the figure, coalescence events belonging to both regimes have been represented. The processes shown in figure 13(*a*,*b*) represent typical collisions taking place in the first regime. It can be seen that in both cases the coalescence happens as soon as the two bubbles collide, around $t=-20$ ms in both series. In these cases, the bubbles do not significantly change their shape when they interact and rapidly contact at a point before coalescing. It could be said that they meet and kiss each other. This type of coalescence was reported for unconfined configurations by Howarth (Reference Howarth1964) and later on modelled as the ratio between the interfacial energy and the energy of collision (Sovova Reference Sovova1981; Tsouris & Tavlarides Reference Tsouris and Tavlarides1994). However, the collision shown in figure 13(*c*) is markedly different and belongs to the second regime. In this case, when the bubbles get close enough (after $t=-100$ ms) they deform and flatten, with the upper bubble surrounding the lower one. They move together for a while before being able to open a hole in the liquid film that separates them, and coalesce. During this period, the lower bubble, $D_k$ (red dots), decelerates before contacting the upper one, $D_{k'}$ (blue dots), using its kinetic energy to deform the interface (Chesters Reference Chesters1991; Kamp *et al.* Reference Kamp, Chesters, Colin and Fabre2001) and to increase the pressure in the liquid film that forms between the two bubbles (Duineveld Reference Duineveld1998). Thus, unlike in figure 13(*a*,*b*), the bubbles meet and dance for a while before kissing in figure 13(*c*). The dancing time increases with the size of the bubbles since large bubbles are able to deform more easily, what makes $h(D_k, D_{k'})$ decrease in the second regime. It is worth noting that in the earlier stages of the cascade process, i.e. for low values of $D_{V90}$, the swarm is formed by relatively small bubbles which collide exclusively in the first regime. The time the bubbles take to deform and adapt their shapes before coalescing not only depends on their sizes but also on the liquid velocity fluctuations, induced by the motion of the bubbles in the swarm. Larger velocity fluctuations favour the destabilization of the liquid film separating the contacting bubbles and, thus, their coalescence. To illustrate this effect, figure 13(*d*) shows a pair of bubbles similar to those in figure 13(*c*) in a bubble swarm also characterized by $D_{V90}\approx 18.25$ mm ($D_{V90}/D_0 \approx 4.74$), but with a higher concentration of bubbles, $\alpha _0 = 6.7\,\%$. In this case, although the interacting bubbles are similar to those displayed in panel (*c*), the coalescing time is shorter because the liquid fluctuation velocities are larger for $\alpha _0 = 6.7\,\%$ than for $\alpha _0 = 4.9\,\%$. In fact, it can be observed in figure 13 that at $t=-160$ ms the distance between the centroids of the two bubbles is larger in (*d*) than in (*c*). Therefore, it can be asserted that the global motion of the swarm is driven by buoyancy and, thus, governed mainly by the evolving bubble size distribution, characterized by $D_{V90}$, and the concentration of bubbles, $\alpha _0$ (see discussion of figure 6*b* in § 3.1). Furthermore, the bubble collision rate increases with the liquid perturbation velocity which increases as the bubbles coalesce and get larger, i.e. as $D_{V90}$ increases.

The role played by hydrodynamic interactions on the coalescence process is also illustrated in figure 14, now considering different events where the smallest bubble is above the largest one. In the figure, typical situations of bubbles interacting in the first regime are presented for different values of $D_p$, in a bubble swarm corresponding to $\alpha _0 = 4.9\,\%$ and $D_{V90} = 23.17$ mm ($D_{V90}/D_0 = 6.02$). The relative location of the smallest bubble of the pair, $D_{k'}$, is plotted as it interacts with a larger bubble, $D_k$. The origin of the system of coordinates ($x_o, z_o$) corresponds to the centroid of bubble $D_k$ at each instant. The temporal evolution of the relative position of the centroid of bubble $D_{k'}$, between the represented stages, is indicated by dots. Bubble pairs for increasing values of $D_p$ (and of $h(D_k,D_{k'})$) are presented in figure 14(*a*–*c*). For a nearly constant value of $D_k$ (largest bubble), the size $D_{k'}$ (smallest bubble) is progressively increased from left to right. In figure 14(*a*), illustrating the interaction of a small bubble with a big one, the small bubble barely deforms, and is ejected away from the larger one due to the overpressure established around its stagnation point. Consequently, this case is not considered as a collision in our analysis (see § A.1) and the collision frequency of this type of bubble is low. As $D_{k'}$ increases (figure 14*b*,*c*), in addition to increasing the possible length of interaction with bubble $D_{k}$, the capability of $D_{k'}$ to be deformed also increases, favouring the collision between both bubbles. Besides the hydrodynamic mechanisms governing the response of bubble $D_{k'}$ to the flow around the top of bubble $D_k$, additional bubble deformation effects have been also observed when the two bubbles are sufficiently large (figure 14*d*). In this kind of collision, the liquid velocity field brings the bubbles sufficiently close so that bubble $D_{k'}$ deforms towards the low pressure area of the wake of bubble $D_{k}$. Indeed, these final instants of the interaction process, based on the wake entrainment mechanism, lead to bubble collision when the interacting bubbles are large enough to deform (Miyahara, Tsuchiya & Fan Reference Miyahara, Tsuchiya and Fan1991). In contrast, smaller bubbles that barely deform and respond faster to the liquid velocity fluctuations can eventually avoid the collision (Filella, Ern & Roig Reference Filella, Ern and Roig2020).

For figure 13(*a*,*b*), the discussion was focused primarily on the short time interaction, as well as on the wake effects of the bubbles once they were carried out close enough by the liquid motion. In contrast, the configurations shown in figure 14(*a*–*d*) focuses on the response of $D_{k'}$ to the hydrodynamic effects caused by its interaction with bubble $D_k$. They illustrate the influence of the bubble sizes, and of their capability to be deformed, on the increase of $h$ with $D_p$ in the first regime. The cases shown in figures 13 and 14 are representative of the collision phenomenology observed in the coalescence cascade process and contribute equally to the bubble pair collision frequency obtained from the experiments.

## 4. Formulation of the coalescence model

To model the rate of change of the population of bubbles due to coalescence, $\dot Q_c$, in (1.2), it is necessary to implement models for the integral kernels involved in the coalescence rate in (1.5), which consider the interaction between bubbles of different sizes, namely $D_k$ and $D_{k'}$, respectively. In this sense, it has to be taken into account that bubble coalescence includes a first approaching step, starting at distances typically larger than the bubble sizes, driven by the local transport mechanisms. Subsequently, a second step, characterized by the short-distance bubble-bubble interaction, takes place (Marchisio & Fox Reference Marchisio and Fox2013). The whole process is commonly modelled considering both stages separately, the first accounting for the bubble collision frequency $h(D_k,D_{k'})$ and the second for the coalescence efficiency term $\lambda (D_k,D_{k'})$. The efficiency is usually expressed comparing the time required to drain and disrupt the liquid film formed between the two bubbles when they collide, $t_d$, with the time the bubbles remain in contact under the influence of the external flow, i.e. the residence time $\bar {t}_c$.

As mentioned above, in the two-dimensional confined configuration of interest here, the coalescence efficiency $\lambda$ has been defined as the fraction of pairs of colliding bubbles which end up coalescing. In this case, once the bubbles collide, they cannot move away in the direction perpendicular to the walls and, in most of the cases, eventually coalesce. Thus, the coalescence efficiency is high and nearly constant, being thus independent of the size of the colliding bubbles and of the external flow conditions, such as $D_{V90}$ or $\alpha _0$ (figure 4). Consequently, the bubble collision frequency $h(D_k,D_{k'})$ will be considered as the bubble coalescence frequency since both functions behave in the same way and, in the following, we will focus on developing a model for $h(D_k,D_{k'})$, which jointly accounts for the local transport phenomena and for the effects derived from the bubbles deformations and their hydrodynamics interaction when they approach.

In general, $h(D_k,D_{k'})$ has been modelled as a characteristic relative velocity between the two colliding bubbles, $\bar {u}_r(D_k,D_{k'})$, multiplied by a characteristic cross-sectional collision area $S_c(D_k,D_{k'})$, which usually depends on the size of the colliding bubbles,

In particular, Coulaloglou & Tavlarides (Reference Coulaloglou and Tavlarides1977) proposed an expression for $h(D_k,D_{k'})$ based on the collision of molecules in kinetic theory of gases, given by

where $u(D_k)$ and $u(D_{k'})$ are the root mean square of the velocity fluctuations of bubbles of sizes $D_k$ and $D_{k'}$, respectively. The relative velocity between the two bubbles can be established by different mechanisms, depending on the liquid field where the bubbles are immersed, i.e. turbulent fluctuations of the carrier fluid, size-dependent differences in the bubble rising velocities, wake entrainment or shear-layer induced velocity differences, among others. However, in the present flow where the liquid velocity is not externally imposed, the mean motion of the swarm and its agitation are driven by a gravity effect and, thus, the mechanisms controlling the relative approaching velocity must be related to the distribution of sizes in the population of bubbles. Under the present conditions, the liquid fluctuating velocities are generated by the interaction of all the bubbles, whose sizes range from $D_0$ to the largest ones, represented by the characteristic diameter $D_{V90}$. Experimental observations suggest that the swarm-induced agitation (represented by the white eddies in the sketch shown in figure 15) is the mechanism that controls the approach of a pair of bubbles. This mechanism is dominant compared with any other mechanism, such as wake entrainment or buoyancy-induced velocity difference between bubbles of different sizes. The relative bubble approaching velocity, $\bar {u}_r$, will be thus assumed to be the standard deviation of the bubbles velocities, which in the present case is proportional to the liquid velocity fluctuations induced by the largest bubbles in the swarm,

The effect of the concentration of bubbles, $\alpha$, has been included in (4.3), according to Bouche *et al.* (Reference Bouche, Roig, Risso and Billet2012) for a monodispersed confined swarm of bubbles. Such relative velocity can be understood as the fluctuating velocity associated with the integral scale of the self-induced agitation, as it is produced by the largest and more intense bubbles.

Concerning the cross-sectional collision area, experimental observations of the collision phenomena indicate that, for a couple of bubbles to collide, they must be within a certain maximum distance. This characteristic length, $\ell _{c}$, defines a region of the cell plane where the liquid velocity fluctuations generated by the bubble swarm are correlated and are able to bring the bubbles together (represented by a dashed circle in the sketch of figure 15). Under these considerations, $S_c(D_k, D_{k'}, \ell _{c})$ can be expressed as

where $w$ is included for dimensional consistency. Here, $\hat {F}$ is a dimensionless function that depends on the reduced diameter of the pair of bubbles, and accounts for the capability of the bubbles to be deformed. In (4.4) the characteristic length of influence of the external mechanisms transporting the bubbles, $\ell _{c}$, is proportional to the integral scale of the liquid flow, $D_{V90}$. Therefore, this influence length varies as the distribution of sizes evolves with the vertical position. The local concentration of bubbles, however, is expected to affect this correlation length, for instance, by inducing a screening effect due to successive passages of bubbles. In fact, Alméras *et al.* (Reference Alméras, Risso, Roig, Plais and Augier2018) reported that, for void fractions $\alpha \geq 5\unicode{x2013}6\,\%$ within monodispersed swarms of bubbles of size $d \sim D_0$, the correlation length for dye transport scaled with the mean distance between two bubbles, which is given by $d / \alpha$. It can therefore be presumed that, for polydispersed swarms, the wide distribution of sizes may, in fact, decrease the values of the gas volume fraction from which the correlation length shortening occurs. Therefore, it can be considered that a pair of bubbles can be transported close to each other if they are initially separated by a distance shorter than

being its variation with the position determined in turn by the changes of the population of bubbles. Substituting (4.5) into (4.4), one gets

In order to determine $\hat {F}$, the values of $h(D_k, D_{k'})$ divided by $\alpha ^{-0.54} \sqrt {g D_{V90}} D_{V90} w$ have been represented in figure 16 as a function of the normalized reduced diameter, $D_p/w$. It is shown that, for the collision events falling in the first regime (solid symbols), the values of the normalized frequency collapse on a single curve, following a linear dependence with $D_p/w$. Therefore, it can be concluded that $\hat {F}$ is a linear function of $D_p/w$ for the pairs of bubbles colliding in the first regime, leading to the following scaling law for the collision frequency:

Only some moderate data scattering can be noticed for the first regime around this linear curve (solid line in figure 16), revealing that the main physics controlling the process is captured by (4.7). In this sense, (4.3) and (4.5) properly scale the magnitude of the self-induced liquid velocity fluctuations and the interaction region for a pair of bubbles, driving the transport of the bubbles during their approach. In addition, the function $\hat {F} \propto D_p/w$ includes the bubble size dependence and accounts for the bubbles deformation process as they interact (see figures 13 and 14).

On the other hand, as expected from the experimental results described in § 3.2, a clear deviation is observed for the collision events taking place in the second regime (hollow symbols in figure 16). These events happen for values of $D_p/w$ larger than the corresponding $\tilde {D}_p/w$, established at each stage of the evolution of the swarm for the different $\alpha _0$ (figure 12*b*) and, thus, when both bubbles are relatively large, i.e. close to $D_{V90}$. In this second regime, most of the coalescence time is spent on the bubbles deformation process, with the approaching time of the bubbles being a small fraction of the total time. We further assume that the bubbles deform and adapt their shapes due to the mean local strain (see figure 13*c*),

acting in a volume enclosing the pair of bubbles $V_c\sim D_{V90}^2 w$. Thus, the bubble pair collision frequency can be estimated by $h(D_k,D_{k'}) \sim s \times V_c$, providing

Equation (4.9) indicates that, in the second regime, $h(D_k,D_{k'})$ decreases as $D_p/w$ increases as shown in figure 16. Figure 17(*a*) shows the experimental values of the bubble pair coalescence frequency that fall in the second regime vs the model given by (4.9), exhibiting an excellent agreement. In fact, matching equations (4.7) and (4.9), it can be inferred that

as shown in figure 12(*b*). Furthermore, since the bubble pair Weber number, ${We= \rho \bar {u}^2_r D_p / \sigma }$, has been typically used to describe the drainage and rupture of the liquid film in drops/bubbles coalescence models (Chesters & Hofman Reference Chesters and Hofman1982), (4.9) can also be expressed in dimensionless form in terms of the Weber number as

where $We= \rho \alpha ^{0.92} (g D_{V90}) D_p/\sigma$. The experimental values of $\hat {h}(We)$ are shown in figure 17(*b*) vs the Weber number for $\alpha _0=$ 3.2 %, 4.9 % and 6.7 %, respectively. Again, the agreement between the experimental measurements and the model proposed by (4.11) is excellent, following a linear dependence of $\hat {h}$ with $We$. In addition to the different mechanisms driving the coalescence of bubbles in both regimes, resulting in the models proposed in (4.7) and (4.9), it is worth noting that the contribution of both types of collisions to the global evolution of the swarm is also quite different. In fact, given that the number of large bubbles in the swarm is considerably low, the contribution to bubble coalescence within the second regime is less significant than that in the first one (figure 11*e*–*h*). However, the coalescence events in the second regime lead the evolution of the tails of the distribution of sizes. This type of coalescence generates even larger bubbles, which in turn eventually establish the velocity and length scales governing the coalescence cascade process.

Finally, the effect of the bubble pair collision frequency on the evolution of the coalescence cascade process can be generally assessed examining the results of the global evolution of the swarm displayed in figure 5. Indeed, the mean collision frequency ${\langle h \rangle }_\infty$ in the swarm can be obtained considering the mean coalescence frequency ${\langle g_c \rangle }_\infty$ and the total number of bubbles in the population $N_\infty$, which assuming a constant efficiency $\lambda _\infty$, reduces to ${\langle h \rangle }_\infty = {\langle g_c \rangle }_\infty / (N_\infty \lambda _\infty )$. This mean bubble pair collision frequency represents the averaged frequency at which a bubble collides with other bubbles of any size, that according to the PBE can be determined as

where $D$ and $D'$ represent the diameters of bubbles whose corresponding volumes are $v$ and $v'$, respectively. Figure 18 shows that ${\langle h \rangle }_\infty$ follows a $3/2$ power law with $D_{V90}$, as inferred from (4.7), for the collisions of bubbles in the first regime. This result suggests that most of the collisions controlling the evolution of the number of bubbles in the swarm take place in the first regime, according to the large amount of small bubbles in the distribution (figure 9). However, it has been observed that ${\langle h \rangle }_\infty$ slightly increases with $\alpha _0$, most likely due to a combined, integral effect of the volume fraction of bubbles on the bubble size distribution in (4.12), $n(D)$, and of the increased relevance of the coalescence events occurring in the second regime with $\alpha _0$. Thus, in figure 18, ${\langle h \rangle }_\infty$ has been divided by a constant, $C_\infty (\alpha _0)$, that slightly increases with $\alpha _0$.

## 5. Conclusion

The bubble coalescence cascade has been analysed for a high Reynolds number confined swarm. Bubble-induced agitation is the main mechanism driving the process, which in turn is highly dependent on the bubble size distribution. In this configuration, a significant evolution of the bubble population is observed, with sizes that grow from the injected bubble size up to 10 times larger, even for the moderate gas volume fractions tested (up to 6.7 %).

A detailed characterization of the coalescence process has been performed by direct observation of the coalescence events taking place in the swarm. In particular, the bubble pair collision frequency and the coalescence efficiency have been measured. Due to the confinement, the coalescence efficiency is considerably high and nearly constant in the present configuration. Thus, the bubble pair coalescence frequency is proportional to the collision frequency and we can talk indistinctly of any of them. This result motivated us to focus our attention on the collision frequency in order to elucidate the underlying physics of coalescence as a global process. Such a process includes a first agitation-driven approach of a pair of bubbles and a subsequent drainage and rupture of the liquid film separating the two bubbles before they coalesce.

Comparing the downstream evolution of initially monodispersed bubble populations, obtained for various injected gas volume fractions, it has been shown that the distributions of sizes evolve following a similar cascade of coalesce events which can be characterized, independently of the concentration of bubbles, by $D_{V90}$ representing a typical size of the largest bubbles in the distribution. However, the rate of change of the size distribution depends on the bubble concentration.

We provide experimental evidence that three parameters control the bubble pair collision/coalescence frequency, $h(D_k,D_{k'})$. These are the bubble pair reduced diameter, $D_{p}$, which accounts for the bubble pair deformation, the diameter $D_{V90}$ and the local gas volume fraction $\alpha$, both characterizing the population of bubbles and the resulting agitation. The interplay of these parameters is however complex and two different regimes of coalescence have been identified. For low values of $D_{p}$, the bubbles are first transported close to each other by the agitation induced in the swarm and as soon as they get close, they collide and eventually coalesce. In this first regime, $h(D_k,D_{k'})$ increases linearly with $D_p$. In contrast, for larger values of $D_p$, pairs of relatively large bubbles interact in a regime mainly controlled by the bubble deformation dynamics. In this case, once the bubbles are close, they elongate and deform adapting their shape as they move together before being able to break the liquid film that separates them, and coalesce. In this second regime, $h(D_k,D_{k'})$ decreases with $D_p$. The characteristic velocity governing the coalescence process is considered to be associated with the integral scale of the liquid motion induced by the largest bubbles in the swarm. This velocity has been estimated as $\alpha ^{0.46} \sqrt {g D_{V90}}$, which includes the effect of the local concentration of bubbles according to Bouche *et al.* (Reference Bouche, Roig, Risso and Billet2012). We conjectured that the coalescence interaction surface is associated with a correlation length of the swarm-induced agitation defined by the successive passages of bubbles across otherwise correlated motions as $\ell _c \sim D_{V90}/\alpha$. Thus, the bubble pair collision frequency in the first regime scales as ${h(D_k,D_{k'}) \sim \alpha ^{-0.54} \sqrt {g D_{V90}} D_{V90} D_p}$ for ${D_p/w\ge 3.52\approx D_0/w}$. However, in the second regime, observed for bubble pairs of reduced diameter greater than $\tilde {D}_p$ provided in relation (4.10), most of the coalescence time is dedicated to deform the bubbles and adapt their surface to each other due to the strain induced by the liquid field, given by ${\alpha ^{0.46} \sqrt {g D_{V90}}/D_p}$. Consequently, in the second regime, ${h(D_k,D_{k'}) \sim \alpha ^{0.46} \sqrt {g D_{V90}} D^2_{V90} w/D_p}$, characterized by the product of the strain rate and the characteristic interaction volume, ${V_c\sim D^2_{V90} w}$. The dependence of the collision frequency in both regimes on the diameter of the largest bubbles, $D_{V90}$, and on the gas volume fraction, $\alpha$, strongly supports the idea that the overall excitation of collisions is a consequence of swarm-induced agitation and not of the relative terminal velocities of the bubbles. The velocity fluctuations of the liquid agitate the bubbles, and depending on their respective sizes, and consequently, on $D_p$, they collide in one or the other regime.

Our results indicate that $D_{V90}$ is essential to understand the bubble cascade process in the present study. In fact, it characterizes the evolution of the population of sizes, but it also mainly drives the self-induced agitation in the swarm that controls the evolution of the bubble population. Regardless, as these confined flows are characterized by lightly interacting wakes of bubbles, it would be interesting to explore if $D_{V90}$ also plays such a relevant role in a swarm of bubbles free to move in an unconfined volume of liquid when the Reynolds number of their relative motion is moderate. In that sense, it is expected that the behaviour of confined swarms differ from that of unconfined ones, not only because of the reduced mobility of the bubbles one around the other in a pair, but also because their hydrodynamics are quite different. However, considering that the characteristic fluctuating velocity is that induced by the largest bubbles in the swarm, as proposed in expression (4.3), we believe that the proposed models may remain partially valid for three-dimensional inertial swarms of bubbles. Also the role played by $D_p$ will be relevant in three-dimensional flows, although the dependency may be different due to the additional degree of freedom of the bubble motion. In this regard, the cross-sectional collision area given in (4.6) or the enclosing volume defined for the collisions in the second regime, included in (4.9), should be redefined in the three-dimensional configuration. Nonetheless, the two-dimensional configuration analysed here has very promising applications, such those concerning light-activated reactions or cultivation of micro-algae, among others. In a future work we plan to investigate the agitation in the swarm, which plays a crucial role in the coalescence cascade. We also plan to apply the models developed here to describe the evolution of the bubble size distribution in confined swarms with different distributions of sizes at the injection point.

## Funding

This work has been partially supported by the Spanish MINECO and European Funds under projects DPI2017-88201-C3-2-R and by the Programa Operativo FEDER Andalucía 2014–2020 and Consejería de Economía y Conocimiento de la Junta de Andalucía under project 1263528. This work has also been financially supported by the French Agence Nationale de la Recherche (ANR) under reference ANR-19-CE43-0002-02 (ALLIGATOR project). J.R.R. wants to acknowledge the Spanish MINECO for the financial support provided by the fellowship BES-2015-071329. Support from the Red Nacional para el Desarrollo de la Microfluídica, RED2018-102829-T, is also acknowledged.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A

#### A.1. Image processing methods for bubble detection and classification

Since even the smallest bubbles in the swarm are already larger than the thickness of the cell, the surface of each bubble is mainly perpendicular to the direction of the light, showing up in the images as a region of connected pixels with a grey level similar to that of the liquid background, enclosed by a thin dark stripe representing the air–water interface not aligned with the cell plane (figure 19*a*). The thickness of the dark line that delimits the perimeter of the bubble is almost constant, regardless of its size, because the curvature of the bubble in a plane perpendicular to the field of view is constant (Bongiovanni, Dominguez & Chevaillier Reference Bongiovanni, Dominguez and Chevaillier2000). Detailed images of isolated bubbles within the same experimental facility and for the same range of sizes reached in this work can be found in Roig *et al.* (Reference Roig, Roudet, Risso and Billet2012). Using this property, each image is analysed making use of a specifically developed image processing algorithm, followed by a bubble detection and classification method. The processing algorithm involves a first pre-processing step, followed by a second one where the image is binarized. Once the bubbles present in each image are detected, their centroid positions as well as their projected areas are measured. Figure 19 shows an example of the processing steps in one of the measuring windows.

The first step of the process implies the improvement of the contrast of the original grey scale image (figure 19*a*). It involves the subtraction of a background reference image without bubbles and the normalization of the image brightness by correcting each value of the pixel intensity matrix. More detailed information regarding the brightness correction method can be found in Fu & Liu (Reference Fu and Liu2016). This brightness normalization reduces the uncertainties due to the variation of the illumination conditions and facilitates the next binarization step. Figure 19(*b*) shows an inversion of the resulting corrected image, where any physical noise (e.g. glass wall scratches and background noise) has been removed while the grey-level gradient between the bubbles edges and the background has been enhanced. Afterwards, the well-known Otsu's method (Otsu Reference Otsu1979) is used in the second step as an automatic, robust, global binarization-threshold selection technique.

After binarization, single bubbles can be detected as blobs of connected low-level pixels enclosed by unique edges of high-level pixels which are totally surrounded by background. Bubbles involved in a collision share the same edge of connected high-level pixels, however, there exists an independent blob of low-level pixels for each bubble. Unlike in other works which deal with bubble collisions or formation of clusters through a separation distance criteria (see Figueroa-Espinoza & Zenit Reference Figueroa-Espinoza and Zenit2005; Figueroa-Espinoza *et al.* Reference Figueroa-Espinoza, Mena, Aguilar-Corona and Zenit2018, for example), the contact of at least one edge pixel of each bubble is required in the present work to define a collision. This definition is crucial, since it determines the measurement of the collision rate and the efficiency of coalescence. Figure 19(*c*) shows detected single bubbles as filled objects and colliding bubbles as hollow ones. This procedure to detect the bubble collisions allows us to clearly distinguish between two independent bubbles involved in a collision (see the example highlighted by an arrow in figure 19) and the newborn ones, just formed due to the coalescence of two colliding bubbles (boxed by dashed lines in figure 19).

Once the bubbles have been detected and classified as single or in-collision bubbles, their instantaneous characteristics are determined. The bubble position is obtained as the centroid position of the in-side blob of low-level pixels. The projected area is defined as the area occupied by the pixels belonging to the in-side blob plus those belonging to the edge. A difficulty arises obtaining the area of the bubbles involved in collisions, since the pixels composing the edge of the agglomerate are shared among the colliding bubbles. To deal with this issue, the total area of the agglomerate edge is distributed among the bubbles according to the ratio between the number of pixels enclosing the inner perimeter of each bubble and those composing the outer perimeter of the agglomerate. Since the edge width remains constant independently of the size of the colliding bubbles, this simple procedure allows us to avoid any further computation devoted to the separation of the bubbles, as occurs in more complicated three-dimensional bubbly flows (see, e.g. Rueda Villegas *et al.* Reference Rueda Villegas, Colombet, Guiraud, Legendre, Cazin and Cockx2019).

#### A.2. Bubble tracking and coalescence/breakage detection algorithm (BTA)

The time evolution of each bubble in the swarm was obtained using a BTA specifically designed and developed for this work, which includes a coalescence/breakage detection algorithm. It consists of the detection of the bubbles in a frame, $j$, followed by the search and identification of the same bubbles in the previous one, $j-1$. For new bubbles, born in frame $j$ either due to coalescence or breakup, family trees are established between the daughter (in frame $j$) and the parents (in frame $j-1$).

More precisely, the algorithm involves a first step where an image (frame $j$) is processed using the digital analysis described in § A.1. As a result, the bubbles in the frame are detected, obtaining the positions of their centroids as well as their projected areas. Afterwards, each bubble is classified as single or in-collision bubble. In addition to the bubbles, the detected collisions, defined as agglomerates of two or more bubbles in direct contact (see § A.1), are treated as independent entities and, thus, their characteristic parameters are also calculated, including the total number of bubbles in collision. In order to facilitate the search of corresponding objects in two consecutive frames, a bounding box containing the target object is defined for each bubble or agglomerate of bubbles (colliding bubbles). Figure 20 shows the BTA performance superimposed on the original grey scale image at various instants, being the reference time, $t=0$, the frame where a coalescence event takes place (figure 20*g*). The trajectories of the bubbles are represented as a sequence of dots, which correspond to the locations of their centroids in the previous frames. The solid circles on each bubble denote the position of the centroid in the previous frame, $j-1$. The above mentioned bounding boxes are plotted in figure 20(*a*–*d*) for two different bubbles in red and blue, respectively, while those corresponding to a collision event are shown in figure 20(*e*,*f*) in green. Notice that the bounding box enclosing the collision becomes that of the newborn coalesced bubble in figure 20(*g*) since the collision event ends up with the coalescence of the two bubbles, as described in detail below.

Once the bubbles in the images (frame $j$) have been detected, the key point of the tracking algorithm is to search for the corresponding ones in frame $j-1$. In that sense, every single and in-collision bubble must be related, at least, to one object in the previous frame. Moreover, every detected collision must be related to a previous collision or, at least, to two previous bubbles. The procedure works sequentially identifying objects, taking into account the continuity of the bubbles trajectories and the conservation of volume (projected areas) of the objects. Initially, only single bubbles in both frames are considered. Therefore, a single bubble in frame $j$ is related to the single one in frame $j-1$ whose centroid falls inside the bounding box of the bubble in frame $j$. Given the experimental acquisition rate, this bounding-box criteria is highly effective, even for the smallest bubbles that are accelerated when they are trapped in the wake of larger ones (figure 20*c*,*d*). However, when bubbles of very different sizes get closer (without being in contact), more than a centroid detected in frame $j-1$ can be inside the bounding box of the larger bubble in frame $j$. To avoid possible errors, an additional criterion based on the conservation of volume is imposed (Rodríguez-Rodríguez, Martínez-Bazán & Montañés Reference Rodríguez-Rodríguez, Martínez-Bazán and Montañés2003). Therefore, the volume of the corresponding bubbles in both frames must be equal. Bubbles located near the bottom edge of the field of view in frame $j$, which cannot be associated with any object in frame $j-1$ (see, e.g. the bubble highlighted by an arrow in figure 20*b*), are directly classified as new bubbles just entering the analysis region. Any other single bubble that cannot be related to a previous one is taken out for further analysis. The single bubbles properly tracked are stored in the data base and removed from both frames to facilitate the analysis of the colliding ones.

The following step is devoted to the analysis of the agglomerates, which define collision events, detected in frame $j$. When several bubbles collide, the process can end with the bubbles coalescing or bouncing off each other. Thus, every collision detected in frame $j$, considered as a unique object, is analysed searching for the corresponding collision object in frame $j-1$, applying the algorithm described above to track individual bubbles. In addition to the volume of the agglomerate, the number of bubbles involved in the collision must also be conserved in both frames. If the corresponding collision is found in the previous frame, the bubbles that form the agglomerate (bubbles marked with coloured stars in figure 20*e*,*f*) are identified as well using the criteria used for single bubbles. On the other hand, if no corresponding collision is found in frame $j-1$, it is assumed that the collision detected in frame $j$ is a new one occurring because several bubbles (usually two) have been brought together (figure 20*e*). In that case, the bubbles involved in the collision are analysed searching for the corresponding previous single bubbles leading to the collision. Once they have been processed, the agglomerates, as well as all the involved bubbles, are stored and no further action is performed with them in the current frame.

At this point, the remaining bubbles in frame $j$ are those emerging either from the breakup of a mother bubble in frame $j-1$ or from the death of a previous agglomerate. The latter gives rise to two different situations: (i) death by coalescence of the bubbles which form the agglomerate, creating a new larger bubble, or (ii) death by the separation of the involved bubbles, leading to different single bubbles in frame $j$. These two possibilities represent the basis of the collision efficiency concept. Situation (i) results in a efficient collision, giving rise to a coalescence event. On the other hand, (ii) indicates an *inefficient* collision, where the involved bubbles continue living without changes in the population. To determine this efficiency, it has to be pointed out that both situations respectively arise from a collision detected in frame $j-1$ which does not have a corresponding agglomerate in frame $j$. Therefore, a forward analysis, from frame $j-1$ to frame $j$, is applied to the remaining agglomerates in frame $j-1$. In this case, the correlation method used to track a single bubble is applied here for each bubble involved in the collision detected in frame $j-1$, searching for the corresponding bubble in frame $j$. Only the remaining bubbles in frame $j$ whose centroid falls inside the bounding box of the analysed agglomerate in frame $j-1$ are considered. If the death of the collision is due to separation, any bubble forming the agglomerate in frame $j-1$ will be related to a corresponding single bubble in frame $j$, satisfying both the bubble bounding box as well as the bubble volume conservation criteria. However, if the collision event ends up in a coalescence, the bubble emerging from the coalescence is determined as the single bubble in frame $j$ whose centroid falls inside the bounding box of the analysed agglomerate in frame $j-1$, being its projected area equal to the sum of those of the parent bubbles. A typical coalescence event is shown in figure 20(*f*,*g*). The parent bubbles involved in the collision can be seen in frame $j-1$ (figure 20*f*), being their centroids indicated by coloured stars, while the newborn bubble is shown in frame $j$ (figure 20*g*), with its centroid marked with a green diamond. From this point, the new bubble generated by coalescence is tracked as a single bubble (figure 20*h*).

Finally, the daughter bubbles remaining in frame $j$, which are generated due to the breakup of a mother bubble in frame $j-1$, are identified through a backward–forward implementation of the correlation method, following the ideas proposed by Rodríguez-Rodríguez *et al.* (Reference Rodríguez-Rodríguez, Martínez-Bazán and Montañés2003). For a potential daughter bubble in frame $j$, the corresponding mother bubble is searched in frame $j-1$ as the larger bubble whose bounding box includes the centroid of the analysed daughter one. Then, the second daughter bubble is additionally searched in frame $j$ as that whose centroid falls inside the bounding box of the mother one and whose volume corresponds to the volume of the mother minus that of the sibling one. Therefore, both daughter bubbles in frame $j$ are identified as new single bubbles appearing due to breakup, while the corresponding mother bubble in frame $j-1$ is defined as death due to breakup.