Clustering, coalescing and bursting processes in surface bubbles of surfactant water flows

Abstract

Breaking waves aerate seawater surfaces and form whitecaps in the open ocean. The aerated surface area, or whitecap coverage, has been used to macroscopically parametrize air–sea momentum and gas exchange. However, the microscopic mechanisms of the generation, evolution and attenuation of surface bubbles in whitecaps remain poorly understood. In this study, we examined the size distributions and size-dependent lifetimes of surface bubbles generated by water sheet entry and air injection on a porous plate during the clustering, coalescing and bursting processes, depending on surfactant concentrations and bubble mobility. Mechanisms of coalescence through film thinning of adjacent bubble walls owing to the inter-bubble attraction and Marangoni forces experimentally described the surfactant-dependent bubble growth, finally achieving bubble bursting, which were statistically characterized in a population balance analysis. Lagrangian bubble lifetimes were described by the Weibull distribution, providing that surfactant in seawater extended the probabilistic survival periods of surface bubbles two times longer than those of clean bubbles.

Impact Statement

The aggregation of breaking-wave-induced surface bubbles on seawater surfaces contributes to the enhancement of air–sea gas transfer, thereby impacting marine ecosystem management. The transient features of oceanic bubble aggregation – dependent on the polydispersivity of bubble size, bubble mobility and surfactant concentration – must be identified to predict their effects on the ocean environment. In this study, we characterized individual bubble behaviours during the coalescing and bursting processes in simplified experiments modelling plunging jets and divergent surface flows through precision image analysis. We observed active coalescence of non-surfactant bubbles colliding through attraction by inter-bubble forces, while aggregating behaviours were predominant at surfactant concentrations that suppress coalescence owing to the Marangoni forces, which modified mechanisms of the total bubble loss and thus probabilistic bubble lifetimes. We also observed another bursting mechanism for travelling bubbles through surface replacement on the bubble cap, which influenced the decay process of surface bubble aggregations. The findings will contribute as the first step in understanding complex oceanic bubble behaviours in surface flows.

1. Introduction

In the surf zone, air bubbles generated by breaking waves are entrained into bulk seawater and entrapped within breaking-wave-induced vortices and turbulence (Watanabe, Saeki & Hosking 2005; Watanabe 2022) (figure 1a), thereby contributing to the dissipation of wave energy (Lamarre & Melville 1991; Callaghan, Deane & Stokes 2016) and dissolution of gases into seawater (Niida & Watanabe 2018). When turbulence weakens, buoyant bubbles are released from the vortices and rise towards the sea surface. Bubbles arriving at the surface may aggregate and remain as a foam structure for long periods, enhancing the chemical and biological activities that support marine ecosystems (Schilling & Zessner 2011). In the open ocean, wind-induced breakers aerate wave faces to create active whitecaps; new foam patches emerge behind the crests through buoyant degassing of entrapped bubble plumes (Deane & Stokes 2002), which exponentially decay over time (Monahan & Muircheartaigh 1980) (figure 1b). The whitecap has roles in exchanging heat and moisture between air and sea (Marmorino 2005; Watanabe & Mori 2008; Wanninkhof et al. 2009) and in producing sea spray aerosols through bubble bursting (Callaghan 2013; Veron 2015), which is significantly influenced by variations of surfactants covering ocean surfaces (Jähne & Haußecker 1998). The area covered by this foam has been used in the parametrization of gas transfer velocity (Monahan & Lu 1990). The extension of the foam area along wave crests, defined as the crest length, has been used as a parameter to estimate wave energy dissipation and flux exchange (Melville 1996). Although macroscopic investigations of photographic whitecap data are commonly performed for oceanic parametrization with wind velocity (Monahan 1971; Callaghan, Deane & Stokes 2008), the fundamental mechanisms of whitecap evolution through the aggregation, coalescence and bursting of bubbles floating on the sea surface remain poorly understood.

Figure 1. Photographs of aerated water generated by shoaling wave splash in the surf zone, taken at Pecatu, South Bali (a) and patches of surface bubbles drifted behind whitecaps, taken at Mombetsu, Hokkaido (b).

Bubble–bubble interactions in bulk liquid flows have been studied both experimentally (Soto et al. 2018) and theoretically (Chesters 1991; Chan, Klaseboer & Manica 2011). When a bubble approaches and collides with a neighbouring bubble, drainage thins the disk-shaped film between the bubbles. If the thinning film is destabilized and ruptured by van der Waals forces, capillary waves generated around the rupture hole rapidly propagate on the coalescing bubble surfaces until the spherical bubble form is recovered (Soto et al. 2018). The bubble may break up under sufficient bulk flow disturbances (Hinze 1955). The mechanisms that underlie the deformation and breakup of a bulk bubble rising in a turbulent flow via wake instability have been examined elsewhere (Ravelet, Colin & Risso 2011). In breaking waves, as surface tension modifies surface forms at the onset of wave breaking (Duncan 2001), subsequent splashing behaviours, characterizing bubble entrainment through interactions between turbulence and surface, are also dependent on seawater surfactants (Liu & Duncan 2006). Deane & Stokes (2002) found that turbulent breakup of the entrained bubbles results in a distinct bubble size distribution of power-law scaling. The void fraction associated with the bulk bubble entrainment is an essential parameter to estimate energy budget in breaking wave turbulence (Na, Chang & Lim 2019). The presence of wave-induced vortices also influences bubble motion; bubbles initially entrained at shallow-level seawater are transported in depths by a Langmuir type of vortical convection (Czerski et al. 2022). In a stage of bubble rise, the bulk vortices modify paths, velocity and drag of buoyant bubbles and enhance aggregation to form bubble clusters through bubble–vorticity interactions (Watanabe et al. 2021).

When rising bulk bubbles arrive at a free surface and float, their local behaviours are governed by distinctive dynamics. Surface bubbles comprise three types of surfaces: a liquid film covering the top of the bubble (i.e. the bubble cap), a meniscus surface surrounding the cap and a submerged bubble wall (Toba 1959). The meniscus forms through a balance of hydrostatic pressure and surface tension on the circular edge of the bubble cap, thereby creating a sloping free surface around the cap (figure 2a). If another bubble is present on the meniscus surface, that bubble must climb the sloping surface because of the tangential component of buoyancy (figure 2c); both bubbles attract each other to the point of aggregation, regardless of whether there is sufficient current to transport them. Nicolson (1949) developed a mechanical model of bubble aggregation based on the apparent attractive force acting between the bubbles. The aggregation of multiple surface bubbles results in the formation of a bubble cluster or raft (Néel & Deike 2021; Watanabe 2022).

Figure 2. Schematic illustrations of (a) bubble-cap film thinning through drainage, (b) retraction of the rim bounding the film as the bubble bursts, (c) bubble attraction on the meniscus surface and (d) neck expansion during coalescence.

The coalescence of contacting surface bubbles may progress under a mechanism similar to that of bulk bubbles; when the film separating adjacent bubbles ruptures, various modes of capillary waves propagate on the cap films of the coalescing bubbles (figure 2d), causing oscillatory deformation until a spherical form is recovered (Shaw & Deike 2021). While the size distribution of clustering bubbles changes during coalescence, cluster attenuation is also caused by bursting of the surface bubble. Lhuissier & Villermaux (2011) described the microscopic bursting mechanism of an isolated stationary surface bubble. The Bénard–Marangoni effect destabilizes the thinning bubble-cap film through gravitational and viscous drainage (figure 2a), resulting in the nucleation of a hole near the edge of the bubble cap. This opening rapidly expands at Taylor–Culick velocity until the cap is eliminated (figure 2b). Néel & Deike (2021) experimentally studied changes in a bubble population caused by coalescing and bursting processes in stationary rafts of monodisperse surfactant bubbles that were generated by air ejection from a circular array of submerged needles.

These findings contribute to the broader understanding of the fundamental coalescing and bursting processes of surface bubbles; however, further understandings of bubble–flow interactions, which may modify the processes of clustering, coalescing and bursting, are required for applications to complex oceanic bubble flows. Additionally, bubble dynamics is strongly dependent on the physiochemical properties of seawater (Slauenwhite & Johnson 1999), which contains numerous surfactants, marine organisms and particulates (Schilling & Zessner 2011) that exhibit seasonal and regional effects according to annual fluctuations of biological activity (Wurl et al. 2011). Because surfactants are responsible for the formation of sea foam and stabilization of bubbles against coalescence and breakup, the surfactant concentration must be considered when predicting the growth and dissipation of whitecaps in the open ocean and in aerated seawater in the surf-zone foam (Callaghan et al. 2012; Callaghan, Deane & Stokes 2013; Villermaux, Wang & Deike 2022) (figure 1). In this study, the mechanisms to induce bubble clustering, colliding, coalescing and bursting are discussed in terms of the experimental parameters of bubble mobility and surfactant concentration; the study objective is to provide fundamental parameters characterizing the transient bubble population, size distribution and bubble lifetime, modified by mobility and organization of surface bubbles, in experiments modelling plunging jets and divergent surface flows, as the first step in understanding complex behaviours of whitecap bubbles in turbulent breaking waves.

The remainder of this paper is organized as follows. Section 2 describes the experimental set-ups, experimental conditions and bulk bubble sizes, as well as image analyses to measure the locations, sizes and travelling velocities of individual bubbles. Section 3 explains behaviours of surface bubbles in coalescing and bursting processes, characterizing the bubble size distributions. Analytical models of the film thinning on adjacent bubble walls, causing coalescence, and the drainage in a bubble-cap film, leading to bubble bursts, are explained in § 4. The regional properties of bubble populations through binary coalescence and their dependency on the surfactant concentration are statistically investigated in § 5. Section 6 discusses contributions of relative frequencies of coalescence and bursting to total bubble loss in population balance analysis. Bubble lifetime probability distributions derived from bubble tracking are examined in § 7. The conclusion is presented in § 8.

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 (Deane & Stokes 2002; Watanabe et al. 2005; Czerski et al. 2022; Villermaux 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 (Cumming & Chanson 1997a), bubble diffusion in a turbulent shear layer formed under the jets (Cumming & Chanson 1997b; Brattberg & Chanson 1998) and the mechanisms of entrainment and breakup (Cumming & Chanson 1999). Chanson (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 (Anagbo & Brimacobe 1990; Abadie et al. 2022) and needles (Risso 2018) have been extensively studied to identify effects of bubble agitation on rise velocity and bubble drag in terms of void fraction. Né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. Cumming & 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 (Kiger & 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 (Cumming & 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 (Anagbo & Brimacobe 1990; Niida & 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 (Wurl 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 (Wurl et al. 2011). For each concentration, the surface tension was given by reference to the experimental work of Wu, 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. Né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 (Watanabe 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 Né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. Watanabe & Ingram 2015) and bubbles (e.g. Watanabe 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. Yuen et al. 1990). Atherton & Kerbyson (1999) combined previous CHT techniques and introduced edge orientation (Kimme, Ballard & Sklansky 1975), phase coding (Atherton & Kerbyson 1993) and a Hough transform filter (Kerbyson 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 (Atherton & 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 (Shaw & 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 Lhuissier & 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 (Chesters & Hofman 1982). For submerged bubbles, coalescence has been modelled by the relative contact time during collision to the film drainage time (Prince & Blanch 1990; Chesters 1991; Sungkorn, 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 (Friedlander 2000; Chan 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 (Nicolson 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) (Langevin 2019). Radoë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 (Prince & Blanch 1990; Sungkorn 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 Couder et al. (2005), and Lhuissier & 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). Lhuissier & 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 (Sungkorn 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$. Né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 (Carrica et al. 1999; Martínez-Bazán, Montañés & Lasheras 1999; Marchisio & Fox 2013; Ruiz-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 (Ruiz-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 Shaw & 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}$.

Né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 Né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 Né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 Lhuissier & 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$. Zheng, Klemas & Hsu (1983) proposed a Rayleigh distribution of lifetimes of stationary surface bubbles generated by air injection through a capillary tube. Lhuissier & Villermaux (2011) suggested $m=4/3$ through reanalysis of the experimental results by Zheng 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.

8. Conclusions

In this study, we experimentally examined the physical properties of clustering, coalescing and bursting surface bubbles generated by water sheet entry (experiment S1) and an air-injected porous plate (experiment S2) based on precise measurements of backlit images. The non-ionic surfactant Triton-X100, which mimics the effects of natural seawater, was used in our experiments to model surfactant enrichment in seawater. The PCCHT method was used to identify the locations and diameters of circular forms representing the surface bubbles, which were then used to quantify the bubble size distributions and track individual bubbles to estimate bubble lifetime.

In the non-surfactant S1 experiment, the largest primary bubble in quasi-stationary, polydisperse surface bubbles attracted and aggregated with nearby smaller bubbles on the meniscus surface, resulting in the formation of bubble clusters where primary bubbles rapidly grew through coalescence with smaller satellite bubbles. The bubble cap of the fully developed primary bubble then ruptured in a process similar to the bursting of a stationary isolated bubble (Lhuissier & Villermaux 2011). In this case, as the frequency of coalescence was much higher than the bursting frequency, total bubble loss was governed by the coalescence process. As surfactant concentration increased, coalescence was suppressed by the Marangoni effect preventing film thinning on bubble walls, resulting in the formation of a bubble raft exhibiting polydisperse distribution spread over the liquid surface. Since bubbles in the raft received omnidirectional inter-bubble forces from neighbouring bubbles, the total force is significantly reduced and thus the frequency of coalescence further decreased. The rates of bubble loss due to bursting in the rafts became predominant, since intermittent bursts emitted capillary waves and disturbed surrounding bubbles, causing successive bursting in the rafts.

In experiment S2, bubbles were transported in a divergent bubble flow, where younger bubbles successively approached older bubbles. When the bubbles came close, attraction force acting on them accelerated younger bubbles towards preceding ones, causing collision. As bubble-wall films of collided bubbles rapidly thinned through drainage by inter-bubble attraction, non-surfactant bubbles successively coalesced with colliding bubbles, resulting in the rapid growth of primary bubbles. Since inter-bubble forces of surfactant bubbles were reduced by the Marangoni force, frequencies of coalescence decreased with surfactant concentration, and instead, patch-like bubble clusters were organized. We detected another potential mechanism for the bursting of travelling bubbles, which constituted cap-film liquid renewal with ambient surface liquid. If disturbances occur on ambient surfaces, the cap film of the travelling bubble is replaced with the disturbed surface liquid, which may destabilize the film. This process might cause successive bubble bursting during bubble drift, resulting in significant depopulation with travel distance. The rates of bubble loss due to both coalescence and bursting were recorded as more than ten times higher than those in experiment S1, resulting in shorter bubble lifetimes than in S1.

Finally, we estimated the Lagrangian bubble lifetime via bubble tracking. We found that the lifetimes had linear correlations with the sizes of bursting bubbles in experiment S2, which was supported by a theoretical time scale of cap-film thinning through drainage. The observed successive bursting of surfactant bubbles in rafts in experiment S1 caused dispersion of their lifetimes from the linear regression. Probability densities of the lifetimes ended by surfactant bubble bursting in experiment S1 have longer tails than those in experiment S2, which may indicate uncertainty associated with successive bursting triggered by intermittent bursts in the rafts. The lifetimes ended by both coalescence and bursting were well approximated by an exponential distribution in both the experiments, while optimal scale parameters of the Weibull distribution of the lifetimes ended only by bursting took values similar to those previously estimated by Lhuissier & Villermaux (2011) and Zheng et al. (1983). Considering the linear correlation between the lifetime and bursting diameter, the lifetime distributions may provide probabilities of occurrence of bursting diameters, and bursting frequencies against bursting sizes. The analysis of the Weibull reliable life led to the conclusion that surfactants in mesotrophic seawater increase bubble lifetime about twice as long as that for clean water.

Through all the analyses of surface bubbles generated for different experimental set-ups and surfactant concentrations, we found that bubble growth and dissipation highly depend on a distinct coalescence process via modification of inter-bubble forces and distinct bursting process including uncertain successive bursting in rafts as well as cap-film replacement, resulting in different behaviours of submerged bubbles in bulk liquid that have been revealed by many previous studies. The current findings provided fundamental parameters characterizing the evolution of surface bubbles, sensitively modified by mobility and organization of surface bubbles, which will be useful for considering suitable parameters required to understand realistic oceanic bubbles in future studies.

All findings obtained in this study relate to a limited range of bubble size, $d>0.75$ mm. Further analysis for smaller bubbles remains to be performed. We also note that the current experiments modelling plunging jets and bubble plumes in a small tank had inherent flows limited by walls, which are different from what is observed in oceans. Further extensions of research are required for applications to whitecap bubbles.