2. Experiments

Characteristic size, shape and motion of bubbles involved in breaking waves vary in time and space through breakup, coalescence and convection in anisotropic turbulence evolving in splashing waves (Reference Deane and StokesDeane & Stokes 2002; Reference Watanabe, Saeki and HoskingWatanabe et al. 2005; Reference Czerski, Brooks, Gunn, Pascal, Matei and BlomquistCzerski et al. 2022; Reference Villermaux, Wang and DeikeVillermaux et al. 2022). In order to avoid direct measurements of such complex bubble-laden turbulent flows, simplified model experiments, focusing on major bubble behaviours in wave plunging and degassing, have been commonly performed to find characteristic aspects of the flows in the breaking process. Air entrainment at plunging planar jets supported by a vertical wall has been experimentally investigated to identify the distributions of void fraction (Reference Cumming and ChansonCumming & Chanson 1997a), bubble diffusion in a turbulent shear layer formed under the jets (Reference Cumming and ChansonCumming & Chanson 1997b; Reference Brattberg and ChansonBrattberg & Chanson 1998) and the mechanisms of entrainment and breakup (Reference Cumming and ChansonCumming & Chanson 1999). Reference ChansonChanson (1995) discussed analogies of turbulent bubble flow under planar jets with that developed in turbulent hydraulic jumps and bores. Bubble plumes generated by consecutive bubble ejection from porous plates (Reference Anagbo and BrimacobeAnagbo & Brimacobe 1990; Reference AbadieAbadie et al. 2022) and needles (Reference RissoRisso 2018) have been extensively studied to identify effects of bubble agitation on rise velocity and bubble drag in terms of void fraction. Reference Néel and DeikeNéel & Deike (2021) used needle-generated bubbles as an experimental model of oceanic bubbles aggregating and bursting on surfactant water surfaces, aiming to find fundamental properties of mass transfer between air and sea. These findings obtained in the model experiments have provided additional interpretations for specific aspects of oceanic bubble flows. In this study, we followed the conventional approach introducing model experiments to generate bubbles, which provided characteristic behaviours of surface bubbles responding in the specific flows explained below.

2.1 Experimental set-ups

We performed laboratory experiments using two bubble generation methods: bubble entrainment by planar plunging jets (S1) and air ejection through a circular porous ceramic plate (S2). In experiment S1 (figure 3a–c), a circulation system consisting of a head tank ($200\,{\rm mm} \times 115\,{\rm mm}\times 100\,{\rm mm}$), transparent test tank ($380\,{\rm mm}{\,\times\,}190\,{\rm mm}\,{\times}\, 215\,{\rm mm}$) and submerged hydraulic pumps produced a steady sheet of falling water that splashed onto still water in the test tank, causing bubble entrainment; this is the simplest model of a plunging breaker jet. The head tank, which was divided into two chambers by submerged walls, had a 2 mm opening over the bottom end of the tank. Liquid was pumped up from both bottom sides of the test tank, then delivered to chamber 1 and overflowed the submerged wall; this cycle prevented disturbed liquid from flowing into the test tank. In chamber 2, the stabilized liquid surface passed down through an opening along one side wall of the test tank as a planar sheet, creating a bubble plume beneath the impacting still water in the test tank. The surface level in the head tank was maintained at a constant water head difference of 100 mm. Reference Cumming and ChansonCumming & Chanson (1997a) have explained evolution of bubble flows under the planar plunging jets, as schematically illustrated in figure 3(c). When the descending planar jet impinges on the still water surface, the resulting velocity gradient beneath the outer jet surface induces a shear stress, which pulls the water near the surface deeper into the bulk for creating a cusp-like cavity close to the jet surface (Reference Kiger and DuncanKiger & Duncan 2012). The cavity elongated downward by the shear becomes unstable and finally breaks up into bubbles to be transported downward with turbulent diffusion in the momentum shear layer expanding towards depths (Reference Cumming and ChansonCumming & Chanson 1997a). While the bubbles are decomposed into smaller ones by turbulent disturbances during descending transport in the shear layer, the diffused bubbles rise up where downward momentum sufficiently attenuates away from the plunging jet. Surface bubbles were illuminated by a red-light-emitting diode panel ($410\,{\rm mm} \times 230\,{\rm mm}$) through the transparent tank bottom. Backlit top-down-view images of the surface bubbles over a field of view of $90\,{\rm mm} \times 70\,{\rm mm}$ were recorded using a high-speed digital video camera (resolution of $1280\,{\times}\, 1024$ pixels) from above the tank by means of a 45$^{\circ }$ angle mirror; the recording frequency was 250 Hz. We defined a coordinate system with the origin at 20 mm from the wall covered by the liquid sheet on the $x$ axis and the centre of the tank on the transverse $y$ axis. Original 8-bit images of backlit bubbles in the domain of $0\le x\le 90$ mm were acquired at a resolution of 0.088 mm pixel$^{-1}$.

Figure 3. Two experimental set-ups for artificial bubble generation. (a–c) Experiment S1: entrainment using planar plunging jets. (d–f) Experiment S2: air ejection through a circular porous ceramic plate. (c,f) Schematic of typical flow patterns formed in each experiment.

In experiment S2, a circular porous plate (diameter of 50 mm) was placed on one bottom end of the test tank (380 mm length$\,{\times}\, 190$ mm width$\,{\times}\, 215$ mm height; figure 3d–f). Air was injected from an air pump to the plate at a constant discharge rate of 500 ml s$^{-1}$ to generate air bubbles in a tank filled with liquid. As described for S1, surface bubbles were illuminated by a light-emitting diode panel through the transparent bottom of the test tank, then recorded by a high-speed video camera from above. Backlit images from the outer edge of the ceramic plate (60 mm from the wall, downstream) over a field of view of $110\,{\rm mm} \times 90\,{\rm mm}$ were acquired at a resolution of 0.126 mm pixel$^{-1}$. The buoyant bubbles produced on the porous plate entrain the bulk water upward, resulting in liquid circulation in a tank with horizontal divergent flow near the surface and convergent flow in the bulk (Reference Anagbo and BrimacobeAnagbo & Brimacobe 1990; Reference Niida and WatanabeNiida & Watanabe 2018) (see figure 3f), which induces active bubble transport on the surface flow.

Two different time scales of measurements were considered in this study. A long time measurement at a low recording frequency (1.0 Hz) for 100.0 s was used for statistical analysis for equilibrium state of bubble population (figures 9 and 12 in § 3, and discussions in § 4). Dynamic behaviours of surface bubbles, including bubble tracking, were analysed by high-speed images recorded at 250 Hz for 5.0 s (1250 sampling frames), which are discussed in §§ 5, 6 and 7.

Aliquots of purified water containing the non-ionic surfactant Triton X-100 at concentrations of $C= 0$, 200, 400 and $800\,{\mathrm {\mu }}$g l$^{-1}$ were used in the experiments. Triton X-100, which mimics the effects of natural seawater, has been used to model surfactant enrichment in seawater (Reference Wurl, Wurl, Miller, Johnson and VagleWurl et al. 2011). The surfactant concentrations used in this study (200, 400 and $800\,{\mathrm {\mu }}$g l$^{-1}$) are within the ranges of mean concentrations for oligotrophic, mesotrophic and eutrophic seawater conditions, respectively (Reference Wurl, Wurl, Miller, Johnson and VagleWurl et al. 2011). For each concentration, the surface tension was given by reference to the experimental work of Reference Wu, Dai and MicaleWu, Dai & Micale (1999): $\gamma \approx 71.0$, 70.2 and 67.6 dyn cm$^{-1}$ for $C=200$, 400 and $800\,{\mathrm {\mu }}$g l$^{-1}$, respectively. Reference Néel and DeikeNéel & Deike (2021) studied influences of surface contamination by sodium dodecylsulphate and Triton X-100 on behaviours of surface bubbles. They used surfactant water at concentrations of 0.25–$50\,{\mathrm {\mu }}$M (corresponding to 162–$32\,350\,{\mathrm {\mu }}$g l$^{-1}$) of Triton X-100 in the experiments, covering the current concentration range. The flow set-ups used in S1 and S2 were controlled to produce identical population densities of surface bubbles at $C=0$, allowing comparisons of surfactant-concentration-dependent effects in both bubble populations with the common population of non-surfactant bubbles.

2.2 Sizes of rising bubbles in bulk liquid

We performed preliminary backlit experiments to measure rising bulk bubbles. Backlit bubbles at medium depth, illuminated by a light-emitting diode panel placed on an adjacent side wall, were recorded from another side of the tank (see superposed images in figure 4). The sizes of the backlit bubbles were estimated using edge detection based on an established level-set method (Reference Watanabe, Oyaizu, Satoh and NiidaWatanabe et al. 2021).

Figure 4. Probability densities of circle-equivalent diameters of submerged bubbles for experiments S1 (a) and S2 (b).

Figure 4 shows the probability distributions of circle-equivalent diameters of submerged bubbles, $d_0$, for both experiments. Bubbles generated in S1 had larger diameters in $d_0 < 1$ mm and were mainly distributed in the class $0 < d_0 < 5$ mm for all surfactant concentrations, featuring polydisperse bubbles. The observed mean diameters were $\overline {d_{01}} = 2.12$ mm ($C= 0\,{\mathrm {\mu }}$g l$^{-1}$), 1.97 mm ($C=200\,{\mathrm {\mu }}$g l$^{-1}$), 1.86 mm ($C= 400\,{\mathrm {\mu }}$g l$^{-1}$) and 1.89 mm ($C= 800\,{\mathrm {\mu }}$g l$^{-1}$). In S2, the probability distribution was narrower (1.0–3.5 mm) and nearly symmetric around the mean diameter. The estimated mean diameters, $\overline {d_{02}}$, were 2.28, 2.10, 2.38 and 2.31 mm for surfactant concentrations of 0, 200, 400 and $800\,{\mathrm {\mu }}$g l$^{-1}$, respectively. Thus, mean bulk bubble sizes were similar between experiments S1 and S2, whereas high polydispersivity was observed in experiment S1. It should be noted that, according to Reference Néel and DeikeNéel & Deike (2021), variations of surfactant concentration have no influence on the sizes of needle-generated bubbles. The current bubble sizes for all concentrations also exhibited an identical distribution with minor fluctuations.

2.3 Locations and sizes of surface bubbles

The edge detection technique is commonly used for estimation of the locations and sizes of sprays (e.g. Reference Watanabe and IngramWatanabe & Ingram 2015) and bubbles (e.g. Reference Watanabe, Oyaizu, Satoh and NiidaWatanabe et al. 2021) on backlit images. In this method, the particle size is determined from the area enclosed by a closed particle edge on an image, and the particle position is given as the centre of the closed area. We observed that the planar shapes of the surface bubbles have circular features, which are typically in contact under aggregation (figure 21a in Appendix A). General edge detection methods may provide erroneous detection of aggregated surface bubbles because the ring-shaped backlit shadows of the bubbles are connected to the shadows of adjacent bubbles, thus resulting in unnecessary areas enclosed by outer edges of the connected ring shadows within the aggregate.

The circle Hough transform (CHT) procedure, a circle detection method, is widely used in image analyses. Multiple improvements have been proposed to detect these inherent locations in parameter space (e.g. Reference Yuen, Princen, Illingworth and KittlerYuen et al. 1990). Reference Atherton and KerbysonAtherton & Kerbyson (1999) combined previous CHT techniques and introduced edge orientation (Reference Kimme, Ballard and SklanskyKimme, Ballard & Sklansky 1975), phase coding (Reference Atherton and KerbysonAtherton & Kerbyson 1993) and a Hough transform filter (Reference Kerbyson and AthertonKerbyson 1995). The resulting phase-coded CHT (PCCHT) was adopted in the present study to estimate the locations and radii of the circular bubbles. In this study, a PCCHT, introducing phase coding and a Hough transform filter (Reference Atherton and KerbysonAtherton & Kerbyson 1999), has been used to detect optimal coordinates and radii for all circular bubbles on the images (Appendix A). Since accurate sizes of very small bubbles with radii unresolved by at lease three pixels were unable to be estimated in this method, sufficiently resolved bubbles with diameters larger than 0.75 mm were used in this analysis.

Using the consecutive bubble locations and sizes, we tracked individual bubbles while the bubble velocity was updated on sequential images, which allowed the estimation of each bubble's lifetime (figure 5). The specific tracking procedure is given in Appendix B.

Figure 5. Example of the estimated bubble velocity (a) and horizontal distributions of the mean horizontal bubble velocity (b) and its standard deviation (c). Dotted and solid lines indicate experiments S1 and S2.

Figure 5(b) shows the horizontal distributions of the mean horizontal bubble velocity in both experiments. We found that surface bubbles in experiment S1 had a stationary feature. Accordingly, while we observed progressive behaviours of young bubbles just emerging on the surface near the source region ($x <2$ cm), the horizontal velocity rapidly became attenuated and the bubbles drifted very slowly; this feature was particularly evident in surfactant bubbles (supplementary movie 1 available at https://doi.org/10.1017/flo.2024.19). In experiment S2, the ascending bubble plume caused by submerged air ejection led to induction of divergent bubble flow on the surface; therefore, the surface bubbles were rapidly transported from the source region outward (supplementary movie 1). The bubble velocity decreased during movement in the divergent flow; thus, younger bubbles with higher velocity at locations closer to the source were able to catch up and collide with preceding older bubbles. These differences in bubble behaviour between experiments resulted in distinct features during coalescence, clustering and bursting, which are discussed in the following sections. When coalescence occurred, the centre location of the coalesced bubble might simultaneously move, while the other adjacent bubbles also rapidly changed their positions (supplementary movie 1). The fluctuations in bubble motion were reflected as standard deviations of local bubble velocity, $\sigma _{u_b}$ (figure 5c). The standard deviations gradually decreased with distance $x$ in both experiments as the frequency of coalescence horizontally attenuates (as discussed in § 5).

3. Behaviours of surface bubbles

3.1 Evolution of bubbles generated by water sheet entry

Figure 6 shows backlit images of non-surfactant bubbles during the coalescence process observed in experiment S1 (supplementary movie 1). Surface bubbles, emerged on the still surface near the source region, formed polydisperse clusters primarily composed of large aggregated bubbles surrounded by smaller satellite bubbles (figure 6a); large bubbles after steep meniscus surfaces tend to attract smaller neighbouring bubbles (figure 2c). We also observed successive coalescence after the production of capillary waves from neck extensions that travelled on the cap films (arrows, figures 6 and 2d), which is similar to the coalescence observed in stationary pair bubbles generated by gas injection (Reference Shaw and DeikeShaw & Deike 2021). However, in the aggregated bubble system, when coalescence-induced waves travelling on the film arrive the cap edge, they disturb and rupture the films in contact with adjacent bubbles (figure 6b). These events generated new capillary waves along the neighbouring bubble cap, which then destabilized the films upon arrival. Thus, through successive coalescence, the initial polydisperse bubble cluster (figure 6a) evolved into bidisperse clusters composed of the ballooning primary bubble surrounded by small satellite bubbles (figure 6c). The growing caps of primary bubbles finally ruptured in a process similar to the bursting of an isolated stationary bubble observed by Reference Lhuissier and VillermauxLhuissier & Villermaux (2011). Specifically, a hole in the film bounded by the circular rim was rapidly retracted on the cap film (figure 2b and arrow, figure 7a), resulting in the formation of an unstable rim with ‘fingers’ at regular intervals along the rim axis (figure 7b).

Figure 6. Backlit images of bubble clusters during the coalescing process at $C=0$ in experiment S1 for (a) $t=t_0$, (b) $t=t_0+\Delta t$ and (c) $t=t_0+2\Delta t$, where $\Delta t=8$ ms. Arrows indicate neck expansion of coalescing bubbles.

Figure 7. Backlit images of bubble clusters during the bursting process at $C=0$ in experiment S1 for (a) $t=t_0$, (b) $t=t_0+\Delta t$ and (c) $t=t_0+2\Delta t$, where $\Delta t= 4$ ms. Arrow indicates the boundary of the cap film.

Figure 8 shows typical surfactant bubbles observed at $C= 400\,{\mathrm {\mu }}$g l$^{-1}$ in experiment S1. Since the surfactant prevented coalescence as the inter-bubble force was reduced by the Marangoni effect (§ 4), the emerged polydisperse bulk bubbles remained on the still surface, resulting in the formation of a dense polydisperse bubble raft that covered a large area of the stagnant surface. The bubbles gradually grew through intermittent coalescence to adjacent ones in the rafts until they burst. We also observed that a single bubble burst in the rafts triggered successive bursting of neighbouring bubbles, which might be also caused by capillary waves, generated by the initial bursting, rapidly propagated to the surrounding bubbles and ruptured their cap films (see arrows in figure 8). Figure 9 shows the spatial distributions of the mean bubble population per unit area, $N (d, x)$, as a function of the bubble diameter in experiment S1. The maximal $N$ in $d<2$ mm extended over the $x$ axis in any cases of $C$, indicating small bubbles uniformly distributed in the whole measurement domain. In the non-surfactant case (figure 9a), there was also a broad size distribution of $N$ in $2\,{\rm mm}< d<10\,{\rm mm}$, indicating bubble growth through successive coalescence over the wide size range (figure 6). As $C$ increased, coalescence was inhibited, thereby reducing the maximum bubble diameter and narrowing the size distribution (figure 9b,c). Because the polydisperse bubble raft covered the most surface at $C = 400$ and $800\,{\mathrm {\mu }}$g l$^{-1}$ (figure 8), the size distributions became horizontally uniform.

Figure 8. Backlit images of bubble clusters at $C= 400\,{\mathrm {\mu }}$g l$^{-1}$ in experiment S1 for (a) $t=t_0$, (b) $t=t_0+\Delta t$ and (c) $t=t_0+2\Delta t$, where $\Delta t=4$ ms. Arrow indicates locations of bubbles that successively burst.

Figure 9. Horizontal variations of bubble size distributions in experiment S1.

3.2 Evolution of bubbles generated by air injection on a porous plate

Figure 10 shows non-surfactant bubble clusters evolving through coalescence in experiment S2. As explained in § 2.3, surface bubbles travel horizontally in the divergent bubble flow (figure 5b). When younger bubbles with higher velocity approached older preceding bubbles, because of attraction on the bubble-induced meniscus surfaces, they were rapidly attracted to each other, leading to collision (see also supplementary movie 1). Thus, bubble sizes increased during travel via consecutive coalescence events after collisions through bubble attraction. For example, the original bubbles, A–E in figure 10(a), repeatedly collided and coalesced with approaching younger bubbles, resulting in rapid growth of the primary bubbles in figure 10(b–e).

Figure 10. Backlit images of bubble clusters in the process of coalescing and bursting at $C= 0$ in experiment S2 for (a) $t=t_0$, (b) $t=t_0+32\Delta t$, (c) $t=t_0+50\Delta t$, (d) $t=t_0+91\Delta t$ and (e) $t=t_0+108\Delta t$, where $\Delta t= 4$ ms. Red letters (A–F) indicate original bubbles before coalescence, for reference.

At a surfactant concentration of $C= 400\,{\mathrm {\mu }}$g l$^{-1}$ in experiment S2 (figure 11), the surfactant suppressed coalescence, which reduced the frequency of collision-induced coalescence events, thereby inhibiting bubble growth through this process. Instead, colliding younger bubbles contacted older bubbles to form patch-like clusters that were transported downstream. These clusters were then disassembled through successive bursting of the component bubbles (arrows, figure 11d–f).

Figure 11. Backlit images of bubble clusters during the processes of coalescence and aggregation and bursting at $C= 400\,{\mathrm {\mu }}$g l$^{-1}$ in experiment S2 for (a) $t=t_0$, (b) $t=t_0+45\Delta t$, (c) $t=t_0+122\Delta t$, (d) $t=t_0+241\Delta t$, (e) $t=t_0+252\Delta t$ and (f) $t=t_0+289\Delta t$, where $\Delta t= 4$ ms. Arrows in (a–c) indicate coalescing bubbles, whereas arrows in (d–f) indicate past sites of vanished bubbles.

The size growth of the travelling non-surfactant bubbles in the collision–coalescence process was identified over a wide size range of $3\,{\rm mm} < d<8\,{\rm mm}$ (see the region surrounded by the dashed curve in figure 12a). With an increase of the surfactant concentration, as the frequency of coalescence decreases, the maximum size achieved through coalescence was reduced (figure 12b,c). We also found that any sizes of bubbles dissipated and vanished within a certain travel distance (see vertical lines in figure 12), unlike uniform distributions in surfactant bubble rafts observed in experiment S1 (figure 9).

Figure 12. Horizontal variations of bubble size distributions in experiment S2. Bubble size growth owing to coalescence during horizontal travel is indicated in the area surrounded by a dashed curve. A vertical line indicates a critical travel distance at which most bubbles have vanished.

3.3 Mobility effect

The dynamic behaviours of the observed bubbles differed between experiments S1 and S2; quasi-stationary bubbles formed clusters and rafts covering the surface in experiment S1, whereas travelling bubbles coalesced and burst through collisions in experiment S2. To understand the effects of bubble mobility on bursting processes, we performed additional experiments to visualize the surface liquid on bubble caps in experiment S2. Immiscible dry powders (DIAION HP20SS; mean diameter, $100\,{\mathrm {\mu }}$m) were sprinkled on the surface as tracers of surface liquid on the caps of travelling bubbles (figure 13). We observed that the particles floating on the still surface near the travelling bubbles moved over the cap film (arrows, figures 13 and 14a); that is, the liquid in the cap film of the moving bubble was replaced with liquid from the downstream surface. Since the film is steadily renewed before undergoing drainage, film thinning, causing a bubble to burst, may be inhibited. However, if a bubble moves towards the disturbed surface, as the disturbances are transported to the cap through the liquid replacement, the disturbed film may become unstable to be ruptured (figure 14b). The observed successive bursting of bubbles travelling in aggregate, as shown in figure 11(d–f), may be triggered by this process; specifically, surface disturbance caused by the bursting of one preceding bubble in the aggregation is transported to the cap of the adjacent bubble, causing that bubble to burst. This repeated process eventually leads to overall collapse of the entire aggregation (figure 11f). Further investigations are required to identify details of the mechanisms, including relative acceleration of fluid on the film and dependencies on surfactants.

Figure 13. Successive backlit images demonstrating that floating particles (arrows) passed over the cap films of bubbles travelling from left to right. The time interval of the images in (a–c) was $\Delta t= 4$ ms.

Figure 14. Schematic illustrations of the top (a) and side (b) views of a surface bubble moving from left to right. (a) Particles floating on the cap surface translating in the opposite direction to the bubble movement. (b) Film liquid on a moving bubble is replaced by liquid from the downstream surface. After a disturbance on the surface, the cap film liquid is replaced with disturbed surface liquid during bubble travel, leading to rupture of the film.

4. Film thinning

When two bubbles approach each other, the fluid between them is displaced and drained in radial flow normal to the approach direction (figure 15a). Coalescence can occur when the fluid layer between them becomes thinner through drainage than the critical value where rupture takes place (Reference Chesters and HofmanChesters & Hofman 1982). For submerged bubbles, coalescence has been modelled by the relative contact time during collision to the film drainage time (Reference Prince and BlanchPrince & Blanch 1990; Reference ChestersChesters 1991; Reference Sungkorn, Derksen and KhinastSungkorn, Derksen & Khinast 2011).

Figure 15. (a) A circular disk-shaped film on the bubble wall of two colliding surface bubbles. (b) Thinning rates of the film thickness, $h/h_0$, estimated by (4.5) (blue), (4.6) (red) and (4.7) (black) at $h_0=100\,{\mathrm {\mu }}$m and $C=0$: — ($R=1.0$ mm), - - - ($R=2.0$ mm), $\cdots$ ($R=3.0$ mm), – $\cdot$ – ($R=4.0$ mm).

The drainage process has often been explained by the Stefan–Reynolds flat-film model (Reference FriedlanderFriedlander 2000; Reference Chan, Klaseboer and ManicaChan et al. 2011), assuming the thinning-film form of a circular flat disk with radius $a$ (figure 15b). Considering axisymmetric viscous flow in the disk film of thickness $h$ between approaching bubble walls under external force $F$, the solution of the Stokes flow with lubrication approximation, integrated over the disk, is given by

(4.1)\begin{equation} \frac{{\rm d} h}{{\rm d} t} = {-}\frac{2}{3{\rm \pi}}\frac{h^3}{a^4 \rho \nu}F, \end{equation}

where $\rho$ is the density of liquid and $\nu$ is the kinematic viscosity.

As explained in § 3, aggregations of surface bubbles were caused by the apparent attraction which is a tangential component of buoyancy on the meniscus surface slopes formed around the adjacent bubbles (see also figure 2 and supplementary movie 1). When we consider a circular disk-shaped film between the bubble walls of two colliding bubbles of identical radius $R$ (figure 15a), neglecting Hamaker forces, the force acting on two colliding surface bubbles (Reference NicolsonNicolson 1949) is given by

(4.2)\begin{equation} F = {\rm \pi}R \gamma q, \end{equation}

where

(4.3)\begin{equation} q = \frac{\beta^3}{\sqrt{4-\beta^2}}\frac{K_1(2\alpha\sqrt{1-a^2/R^2})}{K_1(\alpha \beta)}-2 \frac{a^2}{R^2}-M(C, a), \end{equation}

where the Morton number $\alpha =R\sqrt {\rho g/\gamma }$ and $M(C, a )$ is the Marangoni force describing surfactant effects on drainage in a film ($M$=0 for non-surfactant bubbles) (Reference LangevinLangevin 2019). Reference Radoëv, Dimitrov and IvanovRadoëv, Dimitrov & Ivanov (1974) solved the Reynolds flat-film model including effects of surfactant concentration and found the Marangoni effect reduced drainage to decelerate the film thinning. In the current work, for simplicity, $M$ was dealt with as an undetermined constant, assuming slower temporal variation of $M$ than the thinning velocity. The first and second terms of the right-hand side of (4.3) indicate contributions of the attraction force and the excess pressure, respectively. The radius of a circular edge of a bubble cap relative to the bubble radius, $\beta =b/R$ (figure 15a), is determined by the mechanical balance between the buoyancy and static pressure on the meniscus:

(4.4)\begin{equation} \beta^2-\frac{2}{3}\alpha^2\left(1 + \sqrt{(1-\beta^2)(1 + \beta^2/2)}-\frac{3}{2} \frac{\beta^3}{\alpha\sqrt{4-\beta^2}}\frac{K_0(\alpha \beta)}{K_1(\alpha\beta)}\right) = 0, \end{equation}

where $K_0$ and $K_1$are the modified Bessel functions of the second kind of order zero and one, respectively. Assuming the thinning velocity is much faster than the moving velocity of bubbles (i.e. constant $F$) and integrating (4.1) with (4.2), the thickness of the thinning film between the attracted bubbles is given by

(4.5) \begin{equation} h = \left(\frac{4}{3}\frac{\gamma R q}{\rho \nu a^4}t + \frac{1}{h_0^2}\right)^{{-}1/2}. \end{equation}

For colliding bubbles in bulk liquid, the bubble contact time relative to the film drainage time has modelled the occurrence of coalescence (Reference Prince and BlanchPrince & Blanch 1990; Reference Sungkorn, Derksen and KhinastSungkorn et al. 2011). Once surface bubbles collide, they never detach as the attraction force keeps acting on them (figures 6, 8, 10 and 11). The contact time thus cannot be a parameter modelling the coalescence rate for surface bubbles. From (4.5), it may be appropriate to use the characteristic drainage time $\tau _c\sim \rho \nu a^4/\gamma R q h_0^2$ as a measure of the time of commencing film rupture leading to coalescence. The film thinning becomes slower (with longer $\tau _c$) as $q$ is reduced by the Marangoni effect (equation (4.3)).

Next we consider drainage in a bubble-cap film. Following Reference Couder, Fort, Gautier and BoudaoudCouder et al. (2005), and Reference Lhuissier and VillermauxLhuissier & Villermaux (2011), considering the viscous flow in a curved thin film of the bubble cap, subjected to gravity, the thickness of the thinning film during gravity drainage can be given as a similarity solution for the lubrication equation system of Poiseuille flow:

(4.6)\begin{equation} h = G(\theta) \left(\frac{g}{3\nu R}t + \frac{1}{h_0^2}\right)^{{-}1/2}, \end{equation}

where $G(\theta )=(4/(3 \sin4/3\,\theta)\int _0 ^\theta \sin ^{1/3}\phi \,{\rm d}\phi )^{1/2}$ and $\theta$ is the polar angle. Here $G(\theta )$ is almost uniform between $G(0)=1$ and $G({\rm \pi} /2)\approx 0.76$. The characteristic time $\tau _b\sim \nu R/g h_0^2$ may be used for the reference time required for film rupture resulting in bursting of an isolated surface bubble (without coalescing process). Reference Lhuissier and VillermauxLhuissier & Villermaux (2011) considered capillary effects on drainage velocity and suggested another thinning law:

(4.7) \begin{equation} h\sim \rho g^{1/3} R^{7/3} \gamma^{{-}1} \left(\frac{\nu}{t}\right)^{2/3}. \end{equation}

Figure 15(b) shows the temporal variations of film thickness given by (4.5), (4.6) and (4.7) at $C=0$. Initial thickness ($h_0$) of $100\,{\mathrm {\mu }}$m was assumed for all cases (Reference Sungkorn, Derksen and KhinastSungkorn et al. 2011). Constant $G=1$ was used for (4.6). Equation (4.5) was solved for the bubbles placed 0.1 % closer to the contact central distance ($2R$); i.e. $a=0.0316R$. Because of active drainage in the bubble-wall film between attracted bubbles, the film thickness, described by (4.5), decreased with much faster velocity than that of a cap film owing to natural gravity drainage, (4.6) and (4.7), indicating that rupture of the bubble-wall film, leading to coalescence, occurred earlier than bursting through slower thinning of the cap film (see figures 6, 10 and 11 and also supplementary movie 1). In surfactant cases, as drainage velocity was decreased by the Marangoni force, resulting from surfactant gradients over the bubble-wall film (as $M > 0$ in (4.3), the inter-bubble force $F$ in (4.2) decreased), the film thinning decelerated, resulting in inhibition of coalescence (see figure 12c). When a raft of surfactant bubbles was formed (figures 8 and 9b,c), as omnidirectional attraction forces worked each other from all contacting bubbles, further reduction of the total inter-bubble force resulted in slower drainage, causing further reduction of the frequency of coalescence.

5. Population through bubble burst and coalescence

In this section, we statistically characterize variations in the bubble population through bubble burst and coalescence following the mechanisms discussed in § 4. The volume of air in bubbles is conserved during coalescence. When the number density of bubbles is discretized in terms of volume classes $v_i\ (i=1, 2, \ldots )$, assuming binary coalescence of monodisperse bubbles in the volume class $v_1$, added volumes of coalesced bubbles are given as $v_2 = 2v_1$, $v_3=v_1+v_2=3v_1$, $v_4=v_1+v_3=2v_2=4v_1$, $v_5=v_1+v_4=v_2+v_3=5v_1$, and so on. Accordingly, the coalesced bubble volumes are simply described by multiples of the initial volume, $v_1$. Reference Néel and DeikeNéel & Deike (2021) found the bubble number density exhibited harmonic variations achieving the maximal peaks at multiples of the mean volume of bulk bubbles through bubble coalescence events in non-surfactant bubble rafts.

Figure 16 shows the bubble density per 0.1 mm$^3$ increment per unit area as a function of the spherical bubble volume. We found oscillatory features of the bubble densities with respect to the bubble volume in any surfactant scenarios. In the non-surfactant experiment S2, because of the inter-bubble attraction force (4.2), bubbles successively collided with each other, leading to successive coalescence owing to rapid film thinning in an early phase of the bubble travel (see also figures 10 and 12 $a$), which resulted in harmonic volume variations at $x= 20$ mm (figure 16a-ii). The harmonic peaks with intervals of 3–4 mm$^3$ correspond to the volume of a spherical bubble in $1.8\,{\rm mm}< d<2.0$ mm, approximately corresponding to the mean diameter of bulk bubbles (§ 2.2). For moderate concentration of $C=400\,{\mathrm {\mu }}$g l$^{-1}$, the harmonic signatures were observed at $x=40$ mm (figure 16b-ii), indicating that coalescence occurred in a later phase than in the non-surfactant case, since the Marangoni force $M$ in (4.3) reduced the total inter-bubble force of (4.2) and decelerated the film thinning. Further reduction of $F$ at the highest surfactant concentration $C=800\,{\mathrm {\mu }}$g l$^{-1}$ significantly suppressed the volume transition through coalescence (figure 16c-ii). In all cases, the bubble densities over volume significantly decreased with distance $x$ (see also figure 12) owing to successive bursting of travelling bubbles (§ 3.3).

Figure 16. Size distribution plots as a function of the bubble volume in experiment S1 (a-i,b-i,c-i) and S2 (a-ii,b-ii,c-ii) for (a-i,a-ii) $C=0\,{\mathrm {\mu }}$g l$^{-1}$, (b-i,b-ii) $400\,{\mathrm {\mu }}$g l$^{-1}$ and (c-i,c-ii) $800\,{\mathrm {\mu }}$g l$^{-1}$.

In experiment S1, since coalescence was significantly inhibited in a raft structure, as explained in § 4, uniform density with minor variations increased over volume without distinct harmonic peaks, indicating that polydisperse coalescence governed in these cases.

6. Frequencies of coalescence and bursting

A Boltzmann-type transport equation has been used to describe the evolution of population of bulk bubbles (Reference Carrica, Drew, Bonetto and LaheyCarrica et al. 1999; Reference Martínez-Bazán, Montañés and LasherasMartínez-Bazán, Montañés & Lasheras 1999; Reference Marchisio and FoxMarchisio & Fox 2013; Reference Ruiz-Rus, Ern, Roig and Martínez-BazánRuiz-Rus et al. 2022):

(6.1)\begin{equation} \frac{\partial n}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(n \overline{\boldsymbol{u}}) + \frac{\partial \phi n}{\partial v} = {S}, \end{equation}

where $n (v, {\boldsymbol{x}},t)$ is the number density of bubbles with volume $v$ at location ${x}$ and time $t$, $\overline{\boldsymbol{u}} (v, {\boldsymbol{x}}, t)$ is the mean horizontal velocity of bubbles of volume $v$, $\phi$ is the rate of change of the volume $v$ of a bubble due to thermal effects (condensation, evaporation and dissolution) and ${S}(v, {\boldsymbol{x}},t)$ is the source term.

When taking into account production source, coalescence and bursting, assuming one-dimensional bubble flow and ignoring thermal effects, (6.1) is reduced as

(6.2) \begin{equation} \frac{\partial n}{\partial t} + \frac{\partial}{\partial x} n \bar{u} = {Q_{c}} + {Q_m} + {Q_b} + {S_p} , \end{equation}

where ${Q_{c}} (v, {x},t)$ and ${Q_{m}}(v, {x},t)$ are the birth and death rates of bubbles due to coalescence, respectively, ${Q_b}(v, {x},t)$ is the death rate of bubbles due to bursting and ${S_p}(v, {x},t)$ indicates the source associated with bubble emergence on a surface. We may also define frequencies of coalescence, $q_m (v)$, and bursting, $q_b (v)$, by $Q_m=q_m n$ and $Q_b=q_b n$, respectively (Reference Ruiz-Rus, Ern, Roig and Martínez-BazánRuiz-Rus et al. 2022). Rates $Q_c$, $Q_m$ and $Q_b$ can be directly estimated in the current temporal analysis through bubble tracking on the basis of sequential locations and sizes throughout lifetimes of individual bubbles (Appendix B).

Figure 17 shows $Q_c$, $Q_m$ and $Q_b$ and their sum $Q_s=Q_c+Q_m+Q_b$ against the bubble volume at three different locations $x= 20$, 40 and 60 mm. As the source of bubbles is sufficiently away in $x\ge 20$ mm (i.e. $S_p \approx 0$), ignoring the first term of the left-hand side of (6.2) at a steady state, the total sink $Q_s$ balances with the bubble density transported from upstream $\partial n\bar {u}/\partial x$. In the non-surfactant case of experiment S1, the relative death rate due to bursting with respect to the total sink, $Q_b/ Q_s$, took values in 0.13–0.28; that is, the major cause of depopulation in this case was not bursting but coalescence. While the coalescence rates decreased by reduction of the inter-bubble force in the rafts (§ 4), the observed intermittent bursts, induced by cap-film thinning, triggered successive bursting of surrounding bubbles (figure 8). Rate $Q_b$ thus relatively increased in the rafts: $Q_b/Q_s\sim 0.72$–0.80 ($x=40$ mm), $Q_b/Q_s>0.95$ ($x=60$ mm) for $C=400\,{\mathrm {\mu }}$g l$^{-1}$; and $Q_b/Q_s>0.87$ ($x=40$ mm), $Q_b/Q_s>0.80$ ($x=60$ mm) for $C=800\,{\mathrm {\mu }}$g l$^{-1}$. Therefore the bursting process is more important for predicting the total bubble loss in the surfactant scenarios. In experiment S2, the rates of bubble loss due to both coalescence and bursting increased more than ten times compared with those observed in S1, especially in the small-volume range $v<5$ mm$^3$, causing significant depopulation during horizontal travel of bubbles (see also figure 16). Rate $Q_m$ took the largest negative values in an early stage of bubble travel ($x=20$ mm) for any surfactant concentration, which decreased with surfactant concentration as the inter-bubble force decreased by the Marangoni effect. In contrast, frequencies of bursting, $Q_b$, increased at any locations with surfactant concentration (e.g. maximum $Q_b/Q_s< 0.1$ for $C=0$, while $Q_b/Q_s$ exceeds 0.95 for $C=800\,{\mathrm {\mu }}$g l$^{-1}$). According to Reference Shaw and DeikeShaw & Deike (2021), when two expanding surface bubbles coalesce, the liquid from the meniscus between the parent bubble caps is brought up in a bulge on the top surface of the film and eventually spreads out over the newly formed bubble cap, suggesting that the liquid mass in the bulge on the new cap recovers the thickness of the cap film. In this case, the cap film thinning might be suppressed through successive coalescence, which possibly resulted in the observed opposite trends of frequencies of bubble bursts and coalescence with respect to surfactant concentration. Further research on changes in film thickness through coalescence is necessary for validation.

Figure 17. Rates of birth $Q_c$ and death $Q_m$ of bubbles due to coalescence and rate of bursting $Q_b$ as a function of the bubble volume in experiment S1 (a-i,b-i,c-i) and S2 (a-ii,b-ii,c-ii) for (a-i,a-ii) $C= 0\,{\mathrm {\mu }}g$ l$^{-1}$, (b-i,b-ii) $400\,{\mathrm {\mu }}g$ l$^{-1}$ and (c-i,c-ii) $800\,{\mathrm {\mu }}$g l$^{-1}$.

Reference Néel and DeikeNéel & Deike (2021) assumed the exponential decay of a raft composed of $\hat {n}$ surfactant bubbles:

(6.3)\begin{equation} \hat{n} = \hat{n_0} \exp({-t/\hat{\tau_r}}), \end{equation}

where $\hat {n_0}$ is the initial bubble population in a raft and $\hat {\tau _r}$ is the raft decay time. The rates of coalescence ($\hat {Q_m}$) and bursting ($\hat {Q_b}$) were defined in terms of the global frequencies of coalescence, $\hat {q_m}$, and bursting, $\hat {q_b}$, by

(6.4a,b)\begin{equation} {\hat{Q_m}} = \hat{q_m} \hat{n},\quad {\hat{Q_b}} = \hat{q_b}\hat{n}. \end{equation}

The exponential decay time was then approximated as

(6.5)\begin{equation} \frac{1}{\hat{\tau_r}}\equiv\hat{q_r} = \hat{q_m} + \hat{q_b}, \end{equation}

where $\hat {q_r}$ is the global frequency of total raft decay. It should be noted that $\hat {q_m}$ and $\hat {q_b}$ are constant in each raft, regardless of bubble size. Therefore, the raft decay time $\hat {\tau _r}$ does not coincide with lifetimes of individual bubbles from birth to death. In this study, mean frequencies of coalescence and bursting over the size and space domains, which may be comparable to $\hat {q_m}$ and $\hat {q_b}$, are defined in terms of the measured $Q_m$ and $Q_b$:

(6.6a,b)\begin{equation} \left.\begin{gathered} \tilde{q_m} = \frac{\displaystyle\int_0^\infty\int_0^{\infty}n q_m \,{\rm d} v\,{\rm d} x}{\displaystyle\int_0^\infty \int_0^{\infty}n \,{\rm d} v\,{\rm d} x} = \frac{\displaystyle\int_0^\infty\int_0^{\infty}Q_m \,{\rm d} v\,{\rm d} x}{\displaystyle \int_0^\infty\int_0^{\infty}n \,{\rm d} v\,{\rm d} x},\\ \tilde{q_b} = \frac{\displaystyle\int_0^\infty\int_0^{\infty}n q_b \,{\rm d} v\,{\rm d} x}{\displaystyle \int_0^\infty\int_0^{\infty}n \,{\rm d} v\,{\rm d} x} = \frac{\displaystyle\int_0^\infty \int_0^{\infty}Q_b \,{\rm d} v\,{\rm d} x}{\displaystyle\int_0^\infty\int_0^{\infty}n \,{\rm d} v\,{\rm d} x}. \end{gathered}\right\} \end{equation}

The corresponding decay time is thus given by following (6.5):

(6.7)\begin{equation} \frac{1}{\tilde{\tau_r}}\equiv\tilde{q_r} = \tilde{q_m} + \tilde{q_b}.\end{equation}

Figure 18 compares the estimated $\tilde {q_m}$, $\tilde {q_b}$, $\tilde {q_r}$ and $\tilde {\tau _r}$ with $\hat {q_m}$, $\hat {q_b}$, $\hat {q_r}$ and $\hat {\tau _r}$ given in Reference Néel and DeikeNéel & Deike (2021). While all the frequencies, $\tilde {q_m}$, $\tilde {q_b}$ and $\tilde {q_r}$, in experiment S1 were smaller than those in S2, overall features of decreasing $\tilde {q_m}$ and unchanging $\tilde {q_b}$ with increasing $C$ were consistent with those of Néel and Deike. We note that there is an analogy of stationary feature of the surface bubbles observed in the current experiment S1 and the previous one by Néel and Deike. The current mean decay time $\tilde {\tau _r}$ for experiment S1 was also comparable to $\hat {\tau _r}$, suggesting that an analogous mechanism (through thinning cap film due to drainage) might trigger the bursting process for stationary bubbles in both the experiments. Time $\tilde {\tau _r}$ of experiment S2, where different organizations of bubble aggregations were observed, indicates shorter decay time than that reported by Néel and Deike. It should be noted that characteristic time scales for coalescence and bursting (§ 4), governed by cap-film thinning due to local inter-bubble forces and Marangoni effects in rafts, highly depend on bubble sizes, as observed in figure 17, which cannot be described by the macroscopic decay time (6.7). Statistical properties of the size-dependent bubble life are discussed in the next section.

Figure 18. Current estimates of mean frequencies of coalescence and bursting, defined by (6.6a,b), and global decay time, (6.7), compared with the previous ones, (6.5), by Reference Néel and DeikeNéel & Deike (2021) (ND): (a) frequencies in experiment S1, (b) frequencies in experiment S2 and (c) decay time.

7. Lagrangian bubble lifetime

In this section, we discuss the Lagrangian lifetimes of individual bubbles, which are measured by bubble tracking (Appendix B). Figure 19 shows scatter plots of the measured bubble lifetime, $T_l$, against the bubble diameter at bursting, $d_b$. As we consider here lifetimes of bursting bubbles owing to ruptures of cap films, the lifetimes ended by coalescence are excluded in the analysis. We found that $T_l$ has weak correlations with $d_b$ in both the experiments. As discussed in § 4, the cap-film thinning through gravity drainage is characterized by the time $\tau _b$ that is proportional to the bubble size. Therefore, a linear relationship between $T_l$ and $d_b$ was assumed and fitted to the lifetime plots by the method of least squares (see lines in figure 19). The optimal rates of the regression slopes and the root-mean-square deviations are given in table 1. In experiment S1, we observed dispersion of $T_l$ with respect to $d_b$, which might result from successive bursts, as observed in figure 8, vanishing adjacent bubbles with arbitrary sizes. In experiment S2, we found better approximations of the lifetimes by linear regression with less deviations than those in experiment S1. Long-life bubbles ($T_l >1.0$ s) grown through coalescence ($d_b>4$ mm) were also observed for moderate surfactant concentration (200 and $400\,{\mathrm {\mu }}$g l$^{-1}$), which may be interpreted in terms of inhibition of bubble-cap thinning owing to film renewal through coalescence, as discussed in the previous section. Since the frequency of coalescence decreased by the predominant Marangoni effect in the highest concentration, both the lifetime and bursting diameter significantly decreased.

Figure 19. Scatter plots of bubble lifetime $T_l$ as a function of the diameter of the bursting bubble $d_b$ in experiments S1 (a) and S2 (b).

Table 1. Optimal proportional constants fitted with a linear equation of $d_b$ and the root-mean-square deviation from the linear approximation.

The size-dependent bubble lifetimes, shown in figure 19, indicate that coalescence, causing bubble growth, significantly affects the bursting process ending the lives of individual bubbles, which depends on local inter-bubble forces, organization of clusters and rafts and bubble mobility. In terms of prediction of total bubble attenuation via these local processes, it is also important to understand probabilistic features of the lifetime distributions. We consider here probabilities of ($a$) bubble lifetimes ended only by bursting and ($b$) lifetimes ended by bursting or coalescence.

The Weibull distribution was previously used to describe lifetime probability for an isolated stationary bubble by Reference Lhuissier and VillermauxLhuissier & Villermaux (2011):

(7.1)\begin{equation} f(T_l) = \frac{m}{\tau_0}\left(\frac{T_l}{\tau_0}\right)^{m-1}\exp({-(T_l/\tau_0)^m}), \end{equation}

where $\tau _0$ is the characteristic time and $m$ is the scale parameter. It is known that (7.1) can describe characteristics of different types of distributions with the value of $m$; i.e. an exponential distribution with $m=1$ and a Rayleigh distribution with $m=2$. Reference Zheng, Klemas and HsuZheng, Klemas & Hsu (1983) proposed a Rayleigh distribution of lifetimes of stationary surface bubbles generated by air injection through a capillary tube. Reference Lhuissier and VillermauxLhuissier & Villermaux (2011) suggested $m=4/3$ through reanalysis of the experimental results by Reference Zheng, Klemas and HsuZheng et al. (1983). While the previous studies used the Weibull lifetime probability for an isolated stationary bubble, the effects of bubble aggregation and mobility on the bubble lifetimes were examined through comparisons with optimal Weibull parameters of the measured bubble lifetimes estimated by fitting the distribution (7.1).

Figure 20 shows the Weibull distributions fitted to the experimental bubble lifetimes ended by bursting. The lifetime distributions in experiment S1 have a longer tail than those in S2 for any surfactant concentration, which may be interpreted by shorter bubble lives ended by cascade bursts of travelling bubbles observed in S2 (§ 3.3). As discussed in § 6, a possible mechanism of cap-film recovery through intermittent coalescence, reducing the frequency of bursting for moderate surfactant concentration, may result in a slightly wider distribution at $C = 400\,{\mathrm {\mu }}$g l$^{-1}$ in experiment S2. The optimal parameters $\tau _0$ and $m$ for the cases of both (a) bubble loss by bursting and (b) total bubble loss are specifically given in table 2. The characteristic time $\tau _0$ in any cases of experiment S1 were larger than those in S2, which was consistent with the relative values of macroscopic decay time $\tilde {\tau _r}$ (figure 18). The values of $\tau _0$ for (b) were smaller than those for (a) in all cases, suggesting that frequent coalescence of short-life bubbles predominated the total bubble loss. We also found the shape parameter $m$ took values close to 1 in (b), indicating that an exponential distribution well approximated the probabilistic lifetimes of total bubble loss. We note that the linear relation of the lifetime and bursting diameter, $T_l\propto d_b$, observed in figure 19, may be able to specify the probability of occurrence of the bursting diameter, as well as the bursting frequency against the bubble volume (figure 17).

Figure 20. Probability densities of bubble lifetimes in experiments S1 (a) and S2 (b). Solid lines indicate Weibull distributions (7.1), with optimal parameters given in table 1.

Table 2. Parameters of a Weibull distribution, and reliable lifetimes.

The Weibull reliable life for a specified reliability $L$ is defined by

(7.2)\begin{equation} \tau_L = \tau_0(-\ln(L))^{1/m}.\end{equation}

If $L=0.5$, (7.2) indicates the median lifetime (half of the bubbles can survive by $\tau _{0.5}$). The estimated $\tau _{0.5}$ and $\tau _{0.1}$ are provided in table 2. Considering mesotrophic seawater, corresponding to $C = 400\,{\mathrm {\mu }}$g l$^{-1}$, both $\tau _{0.5}$ and $\tau _{0.1}$ increase 1.3–1.5 times in (a) and 1.9–2.2 times in (b) compared with non-surfactant cases. Accordingly, surfactant in seawater contributes to suppress the dissipation of bubbles and increase bubble lifetime about twice as long as that for clean water.

The Weibull parameters $\tau _0$ and $m$ should be generalized according to bubble mobility in a future study. The additional effects of turbulence on the parameters should also be determined to predict sea-surface coverage by any surface bubbles that arise from breaking waves.