1. Introduction
Breaking waves in oceans generate bubbles of a wide range of sizes (Blanchard & Woodcock Reference Blanchard and Woodcock1957; Medwin Reference Medwin1970; Johnson & Cooke Reference Johnson and Cooke1979; Trevorrow, Vagle & Farmer Reference Trevorrow, Vagle and Farmer1994; Melville Reference Melville1996; Deane Reference Deane1997; Vagle & Farmer Reference Vagle and Farmer1998; Deane & Stokes Reference Deane and Stokes2002; and others). A key driver in establishing this wide range of bubble sizes is the turbulence that emerges as these waves break (Kanwisher Reference Kanwisher1963; Kitaigorodskii et al. Reference Kitaigorodskii, Donelan, Lumley and Terray1983; Rapp & Melville Reference Rapp and Melville1990; Agrawal et al. Reference Agrawal, Terray, Donelan, Hwang, Williams, Drennan, Kahma and Kitaigorodskii1992; Melville Reference Melville1994, Reference Melville1996; Terray et al. Reference Terray, Donelan, Agrawal, Drennan, Kahma, Williams, Hwang and Kitaigorodskii1996; Gemmrich & Farmer Reference Gemmrich and Farmer2004; Drazen, Melville & Lenain Reference Drazen, Melville and Lenain2008; Deane, Stokes & Callaghan Reference Deane, Stokes and Callaghan2016a,Reference Deane, Stokes and Callaghanb; Mortazavi et al. Reference Mortazavi, Le Chenadec, Moin and Mani2016; and others). A travelling water wave carries gravitational potential energy from the variation of its surface height, as well as kinetic energy in the coherent motion of the water beneath its surface. As these waves break, a significant proportion of these potential and kinetic energies is either expended in entraining air beneath the surface or converted to turbulent kinetic energy. This turbulent kinetic energy is then cascaded from large to small scales, establishing a wide range of characteristic scales of turbulent motion before being dissipated as thermal energy (Richardson Reference Richardson1922; Kolmogorov Reference Kolmogorov1941; Onsager Reference Onsager1945). The resulting turbulence breaks up sufficiently large entrained air cavities into bubbles of various sizes (Kolmogorov Reference Kolmogorov1949; Hinze Reference Hinze1955). The smaller the size of the bubbles produced by these break-up events, the slower their rise velocity and the longer their residence time in the ocean (Garrettson Reference Garrettson1973; Thorpe Reference Thorpe1982, Reference Thorpe1992; Trevorrow et al. Reference Trevorrow, Vagle and Farmer1994). As they rise towards the ocean surface, these bubbles are known to scavenge surfactants and other microparticles. This delays their coalescence with the atmosphere and with one another, and stabilizes them for a longer period of time in the ocean (Fox & Herzfeld Reference Fox and Herzfeld1954; Turner Reference Turner1961; Blanchard & Syzdek Reference Blanchard and Syzdek1970; Johnson & Cooke Reference Johnson and Cooke1980; Weber, Blanchard & Syzdek Reference Weber, Blanchard and Syzdek1983; Johnson & Wangersky Reference Johnson and Wangersky1987; Stefan & Szeri Reference Stefan and Szeri1999; Takagi & Matsumoto Reference Takagi and Matsumoto2011; Czerski Reference Czerski2017; and references therein). Some of these bubbles eventually burst at the surface to form film and jet drops (Blanchard & Woodcock Reference Blanchard and Woodcock1957; Thorpe Reference Thorpe1992; Deane & Stokes Reference Deane and Stokes2002; Veron Reference Veron2015; and references therein). On the whole, then, these bubbles impart an enduring influence on various physical phenomena.
Knowledge of the distribution of bubble sizes informs the total interfacial area, as well as size-dependent effects. These are important for quantifying bubble- and drop-mediated mass, momentum and energy transport near the wave surface (Kanwisher Reference Kanwisher1963; Atkinson Reference Atkinson1973; Thorpe Reference Thorpe1982, Reference Thorpe1992, Reference Thorpe1995; Woolf Reference Woolf1993; Melville Reference Melville1996; Wanninkhof et al. Reference Wanninkhof, Asher, Ho, Sweeney and McGillis2009; Emerson & Bushinsky Reference Emerson and Bushinsky2016; and references therein), the reflection and scattering of solar and other electromagnetic radiation by the water-wave surface and surrounding air–water interfaces (Stramski Reference Stramski1994; Zhang, Lewis & Johnson Reference Zhang, Lewis and Johnson1998; Stefan & Szeri Reference Stefan and Szeri1999; Terrill, Melville & Stramski Reference Terrill, Melville and Stramski2001; Reed & Milgram Reference Reed and Milgram2002; Stramski et al. Reference Stramski, Boss, Bogucki and Voss2004; Zhang et al. Reference Zhang, Lewis, Bissett, Johnson and Kohler2004; Seitz Reference Seitz2011; Twardowski et al. Reference Twardowski, Zhang, Vagle, Sullivan, Freeman, Czerski, You, Bi and Kattawar2012; Crook, Jackson & Forster Reference Crook, Jackson and Forster2016; and references therein) as well as the generation, scattering and propagation of acoustic waves beneath the water-wave surface (Strasberg Reference Strasberg1956; Schulkin Reference Schulkin1968, Reference Schulkin1969; Medwin Reference Medwin1970; Farmer & Lemon Reference Farmer and Lemon1984; Prosperetti Reference Prosperetti1988; Hall Reference Hall1989; Melville Reference Melville1996; Deane Reference Deane1997, Reference Deane2016; Duineveld Reference Duineveld1998; Vagle & Farmer Reference Vagle and Farmer1998; Farmer, Deane & Vagle Reference Farmer, Deane and Vagle2001; Stanic et al. Reference Stanic, Caruthers, Goodman, Kennedy and Brown2009; Czerski & Deane Reference Czerski and Deane2010; and references therein). In order to rigorously quantify these effects, one needs to have a good grasp of the initial distribution of bubble sizes generated in the wave-breaking process, as well as how and why this distribution evolves while the resulting turbulence dissipates and the bubbles rise to the surface.
A number of experiments have measured the bubble-size distribution due to breaking waves generated under a variety of conditions. Among these, the experiments by Deane & Stokes (Reference Deane and Stokes2002) and Blenkinsopp & Chaplin (Reference Blenkinsopp and Chaplin2010) were able to resolve bubbles of sizes spanning over two decades (from under $100\ \mathrm {\mu }\textrm {m}$ to over 10 mm) – a formidable size range for a single measurement set-up. Crucially, these experiments could characterize the size distribution early in the wave-breaking process, as the optical methods employed are suitable for these high-void-fraction scenarios. While there are differences in the reported distributions due to the different wave parameters and measurement techniques, as well as the precise stages of wave breaking that were characterized, there appear to be two common themes: first, the distribution of bubble sizes has distinct scalings in different bubble-size subranges; second, the distribution evolves noticeably in time as the wave breaks. As alluded to in Part 1 (Chan, Johnson & Moin Reference Chan, Johnson and Moin2021), these observations imply that different physical mechanisms are at play at different length and time scales in the formation and dynamics of bubbles in these waves. Similar observations were made in the three-dimensional breaking-wave simulations of Wang, Yang & Stern (Reference Wang, Yang and Stern2016) and Deike, Melville & Popinet (Reference Deike, Melville and Popinet2016), which draw heavily on the foundational two-dimensional simulations of Chen et al. (Reference Chen, Kharif, Zaleski and Li1999) and Iafrati (Reference Iafrati2009, Reference Iafrati2011). Numerical simulations have fewer limitations on the access to detailed spatial information than experimental studies, which is crucial for characterizing these physical mechanisms. However, the separation of scales inherent in these oceanic systems makes it prohibitively costly to generate large simulation ensembles that resolve the range of bubble sizes accessed in these seminal experiments. A large number of ensemble realizations are necessary to achieve statistical convergence for various bubble statistics in these statistically unsteady waves. As such, mechanisms governing bubble formation and dynamics have not all been straightforward to isolate and remain a subject of active research. A number of candidate mechanisms have been proposed for various bubble-size subranges and wave-breaking stages, including bubble break-up by turbulence (Kolmogorov Reference Kolmogorov1949; Hinze Reference Hinze1955), air entrainment and microbubble formation due to plunging jets and drops (Deane & Stokes Reference Deane and Stokes2002; Kiger & Duncan Reference Kiger and Duncan2012; Chan, Urzay & Moin Reference Chan, Urzay and Moin2018c; Mirjalili, Chan & Mani Reference Mirjalili, Chan and Mani2018; Chan et al. Reference Chan, Mirjalili, Jain, Urzay, Mani and Moin2019; Mirjalili & Mani Reference Mirjalili and Mani2020), buoyant degassing and dissolution (Garrett, Li & Farmer Reference Garrett, Li and Farmer2000) and intermittency in the energy dissipation rate (Garrett et al. Reference Garrett, Li and Farmer2000; Gemmrich Reference Gemmrich2010), although this intermittency remains a topic of active investigation (Deane Reference Deane2016; Deane et al. Reference Deane, Stokes and Callaghan2016b). Note that microbubble formation due to plunging jets and drops may entail the direct generation of sub-Hinze-scale bubbles with subunity Weber numbers from larger air sheets and films, which may bypass the action of turbulent break-up (Deane & Stokes Reference Deane and Stokes2002; Chan et al. Reference Chan, Urzay and Moin2018c, Reference Chan, Mirjalili, Jain, Urzay, Mani and Moin2019; Mirjalili et al. Reference Mirjalili, Chan and Mani2018; Mirjalili & Mani Reference Mirjalili and Mani2020). As is discussed in § 3, the simulations of this work do not have numerical support for sub-Hinze-scale bubbles, and thus these bubbles are out of the scope of the current work.
Break-up by turbulence is often cited as the dominant mechanism for the generation of bubbles of intermediate sizes during the active wave-breaking phase, when air is entrained beneath the wave surface. The sizes of these fragmenting bubbles are typically comparable to or larger than the global Hinze scale, at which forces of capillary and inertial origins are approximately in balance (Kolmogorov Reference Kolmogorov1949; Hinze Reference Hinze1955). This characteristic length scale was introduced in Part 1 and is dimensionally about a millimetre in most terrestrial oceanic waves (Deane et al. Reference Deane, Stokes and Callaghan2016a). Because these bubbles are smaller than the integral scales of the system, the large-scale flow geometry should not have a strong influence on the break-up dynamics. Garrett et al. (Reference Garrett, Li and Farmer2000) suggested that the break-up of these bubbles by turbulence occurs through a quasi-steady cascade of bubble mass from large- to small-bubble sizes. (This quasi-steadiness is from the point of view of the small- and intermediate-sized bubbles.) As discussed in the introduction of Part 1, prior studies have theorized and observed in various bubbly flows that this mechanism implies a $D^{-2/3}$ power-law scaling for the break-up frequency of bubbles of size
$D$ (Hinze Reference Hinze1955; Martínez-Bazán, Montañés & Lasheras Reference Martínez-Bazán, Montañés and Lasheras1999; Rodríguez-Rodríguez, Gordillo & Martínez-Bazán Reference Rodríguez-Rodríguez, Gordillo and Martínez-Bazán2006; Chan et al. Reference Chan, Dodd, Johnson, Urzay and Moin2018b), and a
$D^{-10/3}$ power-law scaling for the corresponding bubble-size distribution (Filippov Reference Filippov1961; Garrett et al. Reference Garrett, Li and Farmer2000; Deane & Stokes Reference Deane and Stokes2002; Mortazavi Reference Mortazavi2016; Deike et al. Reference Deike, Melville and Popinet2016; Wang et al. Reference Wang, Yang and Stern2016; Chan et al. Reference Chan, Dodd, Johnson, Urzay and Moin2018b,Reference Chan, Urzay and Moinc). The mathematical formalism in Part 1, which quantifies the average rate of bubble-mass transfer from large- to small-bubble sizes due to break-up events, further demonstrates that these power-law scalings are directly compatible with the notions of locality and self-similarity in the bubble-mass-transfer process. This compatibility confirms key physical aspects of the bubble break-up cascade phenomenology, and provides a theoretical basis for the dimensional analysis of Garrett et al. (Reference Garrett, Li and Farmer2000).
This paper (Part 2) aims to determine the extent to which the theoretical results of Part 1 are supported by numerical simulations, by algorithmically embedding this analytical toolkit and directly measuring its key component metrics in the simulations. While previous breaking-wave simulations have observed a $-10/3$ power-law exponent in the bubble-size distribution, this work seeks a more direct observation of the bubble break-up cascade. The bubble-mass-transfer rate introduced in Part 1 achieves this objective, and is itself dependent on the size distribution, the break-up frequency as well as the distribution of child bubble sizes. In Part 2, these bubble statistics are measured by averaging over ensembles of breaking-wave simulations. They were computed using novel post-processing algorithms that identify and track individual bubbles, and record the details of individual break-up events. The evolution of these ensemble-averaged statistics with time is studied, building on the work of Deike et al. (Reference Deike, Melville and Popinet2016), and the effect of time averaging is elucidated with reference to the time-averaged statistics of Deane & Stokes (Reference Deane and Stokes2002), Wang et al. (Reference Wang, Yang and Stern2016) and Deike et al. (Reference Deike, Melville and Popinet2016). The statistics discussed in Part 2 provide direct support for the existence of the quasi-steady bubble break-up cascade proposed by Garrett et al. (Reference Garrett, Li and Farmer2000) and theoretically examined in Part 1, at least during the early wave-breaking stages. This support lends credence to the usage of this analytical toolkit for examining the break-up dynamics of turbulent two-phase flows in general, and provides an avenue for developing and validating break-up kernels typically used for population balance modelling.
The dynamics of bubble break-up appears to be distinct in the early and late wave-breaking stages, suggesting the emergence of distinct bubble generation and evolution mechanisms. The size distribution was observed to deviate from the aforementioned $D^{-10/3}$ power-law scaling by Deane & Stokes (Reference Deane and Stokes2002), Tavakolinejad (Reference Tavakolinejad2010) and Masnadi et al. (Reference Masnadi, Erinin, Washuta, Nasiri, Balaras and Duncan2020) late in the wave-breaking process. Prior breaking-wave simulations either did not investigate the evolution of the ensemble-averaged size distribution as a function of time, or did not conclusively recover these alternative scalings in their ensemble-averaged size distribution in the late wave-breaking stages. The results of this work indicate that the size distribution and other bubble statistics are indeed evolving functions of bubble size and time. This paper systematically identifies bubble-size subranges and times over which the
$D^{-10/3}$ scaling is not fully recovered in the size distribution. The characteristics of the bubble-mass-transfer rate are examined in these bubble-size subranges and times in order to explore candidate mechanisms for the formation and dynamics of bubbles under these conditions.
This paper is organized as follows. In § 2, the key results of Part 1 are recapitulated and reformulated with the inclusion of volume averaging for the ensembles of breaking-wave simulations described in § 3. Features of the algorithms used to identify individual bubbles in these ensembles, and to detect break-up and coalescence events by tracking the lineages of these bubbles, are briefly discussed in § 4. The resulting bubble statistics are examined in § 5 in relation to the theoretical analysis of Part 1 and § 2. Finally, conclusions are drawn in § 6.
2. Analysis of the bubble break-up cascade for breaking-wave simulations
Volume averaging ($\bar {\cdot }$) is commonly used in the computation of bubble statistics in both numerical and experimental studies of breaking waves. As such, this section recasts the essential expressions of the mathematical formulation for bubble statistics in Part 1 using volume averaging in order to facilitate the subsequent analysis of bubble break-up in breaking-wave simulations. References to §§ 4 and 5 are also made to indicate how these statistics are computed in the simulations and where they are reported, respectively. The reader is referred to Part 1 for a more detailed development and interpretation of the analysis tools summarized here, as well as references to relevant prior art.
Volume averaging may be interpreted as an averaging procedure over the small, localized regions in which the turbulent energy and bubble-mass cascades are postulated to coexist, as reviewed in § 2 of Part 1. The volume-averaging procedure is further discussed in appendix A. The order of magnitude of the correspondingly averaged dissipation rate $\bar {\varepsilon }$ may be estimated in a global sense by
$u_L^{3}/L$, where
$L$ is the wavelength of the dominant wave,
$u_L = (gL)^{1/2}/(2{\rm \pi} )^{1/2}$ is the corresponding wave phase velocity and
$g$ is the magnitude of standard gravity. This estimate for the characteristic dissipation rate is addressed in more detail in appendix B. It corresponds to the following dimensional expressions for the global Kolmogorov and Hinze scales:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn1.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn2.png?pub-status=live)
where $\rho _l$ and
$\mu _l$ denote the density and dynamic viscosity of the liquid phase, respectively, and
$\sigma$ denotes the air–water surface tension coefficient.
2.1. The bubble-mass-transfer flux
The population balance equation describing the time evolution of the volume-weighted size distribution, $fD^{3}$, was discussed in § 3.2 of Part 1. Volume averaging the phase-space-based equation (3.4) in Part 1 yields
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn3.png?pub-status=live)
Correspondingly, volume averaging the kernel-based equation (3.7) in Part 1 yields
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn4.png?pub-status=live)
As with the derivation of (3.10) in Part 1, these two equations may be compared and reduced by assuming that size-local break-up dominates in an intermediate subrange of bubble sizes to yield a simplified evolution equation for $\bar {f}D^{3}$:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn5.png?pub-status=live)
The size distribution, $\bar {f}$, is computed using the identification algorithm in § 4.1, and is further discussed in § 5.1. Note that the atmosphere above the wave is treated as a large gaseous reservoir in this work. As such, entrainment events are not included in the tally of break-up events, and degassing events are not included in the tally of coalescence events. The local break-up flux
$\overline{W_b} = -\overline {v_D f} D^{3}$ may be interpreted as the ensemble-averaged and volume-averaged rate of bubble-mass transfer from all bubble sizes larger than
$D$ to all bubble sizes smaller than
$D$. One may derive an expression for the analogous flux arising from
$\overline{T_b}$, which is obtained by volume averaging (3.11) in Part 1 as follows:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn6.png?pub-status=live)
The ensemble-averaged and volume-averaged differential break-up rate $\overline {g_b f}(D_p;t)$ is the expected number of break-up events per unit time, unit domain volume and unit size for bubbles of size
$D_p$ at time
$t$, including all events throughout the characteristic wave volume
$L^{3}$ of each ensemble realization. Then,
$\check {q}_b(D_c;t|D_p)$ describes the probability distribution of sizes of child bubbles in these events. Note that each relevant break-up event is weighted equally within each ensemble realization in the computation of
$\check {q}_b$. As such,
$\check {q}_b$ satisfies the same constraints as the analogous distribution
$q_b$ in Part 1. The quantities
$\overline {g_b f}$ and
$\check {q}_b$ are both computed using the tracking algorithm in § 4.2, and are further discussed in §§ 5.2 and 5.3, respectively. The corresponding bubble-mass-transfer flux,
$\overline{W_b}$, may then be expressed as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn7.png?pub-status=live)
The transfer flux $\overline{W_b}$ may also be directly computed using the tracking algorithm in § 4.2, and is further discussed in § 5.4. In particular, its variations with bubble size and time are discussed in § 5.4.2. In a system where entrainment at large sizes and break-up towards smaller sizes are the dominant physical mechanisms present, the break-up flux
$\overline{W_b}$ may be used as a proxy for the entrainment flux, as suggested in figure 9 of Part 1. Note that a similar procedure may be used to derive the bubble coalescence flux
$\overline{W_c}$ from an appropriate model coalescence kernel
$\overline{T_c}$.
2.2. Assessing locality
The expressions (3.15) and (3.16) derived in Part 1 to quantify infrared and ultraviolet locality may be rewritten as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn8.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn9.png?pub-status=live)
These quantities may also be directly computed using the tracking algorithm in § 4.2, and are further discussed in § 5.4.1.
In this work, the infrared and ultraviolet locality of $\overline{W_b}$ are investigated using ensembles of breaking-wave simulations in two ways. First, the individual scalings of
$\check {q}_b(D_c;t|D_p)$ and
$\overline {g_b f}(D_p;t)$ with
$D_c$ and
$D_p$ are computed and compared against the scalings derived in Part 1. Second,
$\overline{I_p}$ and
$\overline{I_c}$ are directly computed by summing the mass transfers from all relevant break-up events in the simulations, as outlined in appendix B of Part 1. The decay rates of
$\overline{I_p}(D_p;t|D)$ and
$\overline{I_c}(D_c;t|D)$ with increasing
$D_p$ and decreasing
$D_c$, respectively, are then examined. Before these statistics are discussed in § 5, the breaking-wave simulation ensembles are first described in § 3, and the algorithms used to obtain these statistics are introduced in § 4.
3. The breaking-wave simulation ensembles
3.1. Flow solver
In order to investigate the evolution of the bubble-size distribution $\bar {f}$, as well as the other bubble statistics introduced in § 2, ensembles of numerical simulations of breaking waves were generated using an unstructured, collocated, node-centred and unsplit geometric volume-of-fluid-based flow solver developed for incompressible and immiscible two-phase flows (Ham et al. Reference Ham, Kim, Bose, Le and Herrmann2014; Kim et al. Reference Kim, Ham, Bose, Le, Herrmann, Li, Soteriou and Kim2014; Bravo et al. Reference Bravo, Kim, Ham and Su2018b). Among other flow variables, the solver tracks the spatially discretized volume fraction of one of the two phases
$\phi$, noting that the local sum of the volume fractions of the two phases is 1 by definition. Without loss of generality, it is assumed that the solver is used to simulate a liquid–gas mixture, and
$\phi$ denotes the local volume fraction of the liquid phase. The mixture is assumed to be non-reacting and electrically uncharged, and no mass transfer takes place between the phases. In each interfacial computational median dual cell, the interface is represented using the conventional piecewise-linear interface calculation scheme, while the corresponding interface normal vector is computed via a reconstructed distance function (Cummins, Francois & Kothe Reference Cummins, Francois and Kothe2005), and the corresponding interface curvature is estimated using the second-order direct front curvature method (Herrmann Reference Herrmann2008). Note that the piecewise-linear interface calculation scheme is known to minimize jetsam and flotsam, and thus suppresses spurious bubble break-up, compared to the traditional simple line interface calculation scheme with which these structures are commonly associated (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999). (See also a related discussion, and a more detailed definition of jetsam and flotsam, in § 4.1.) The density and viscosity of the fluid in these interfacial cells are defined as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn10.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn11.png?pub-status=live)
where $\rho _g$ and
$\mu _g$ denote the density and dynamic viscosity of the gaseous phase, respectively, and
$0 < \phi < 1$. Mass and momentum are consistently advected at the interface (Mirjalili, Ivey & Mani Reference Mirjalili, Ivey and Mani2019) using a variant of the non-intersecting flux polyhedron advection algorithm (Ivey & Moin Reference Ivey and Moin2017) to maintain numerical stability for large-density-ratio flows, while the viscous terms in the momentum equation are implicitly time-advanced with the second-order Crank–Nicolson scheme. The fractional-step method is used in order to include a pressure-projection step to maintain a divergence-free velocity field. Within this method, the surface tension force is treated using a balanced-force algorithm (Francois et al. Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006; Herrmann Reference Herrmann2008) involving the usual continuum-surface-force representation to minimize spurious currents of capillary origin. Surface tension gradients are assumed to be negligible. The solver has been used to simulate a number of liquid jet break-up problems involving primary atomization, including laminar, transitional and turbulent jets in quiescent gas (Bravo et al. Reference Bravo, Kim, Tess, Kurman, Ham and Kweon2015, Reference Bravo, Kim, Ham, Matusik, Duke, Kastengren, Swantek and Powell2016, Reference Bravo, Kim, Ham and Su2018b, Reference Bravo, Kim, Ham, Powell and Kastengren2019), jets in cross-flows (Bravo et al. Reference Bravo, Kim, Ham and Kerner2018a) as well as jets from swirling injectors (Kim et al. Reference Kim, Ham, Bose, Le, Herrmann, Li, Soteriou and Kim2014; Ham et al. Reference Ham, Kim, Bose, Le and Herrmann2014).
3.2. Description of set-up of ensembles
In this work, a baseline ensemble of 30 numerical simulations of breaking third-order Stokes water waves in air was generated. Specifically, the density and viscosity ratios of the two phases are equal to those of an air–water mixture. These waves have the dimensionless integral-scale parameters $\mathit {We}_L = 1.6 \times 10^{3}$ and
$\mathit {Re}_L = 1.8 \times 10^{5}$, which match those of a 27 cm long water wave at atmospheric conditions, and are similar to those selected by Wang et al. (Reference Wang, Yang and Stern2016). Here,
$\mathit {Re}_L$ is the integral-scale Reynolds number
$\mathit {Re}_L = \rho _l u_L L / \mu _l$ and
$\mathit {We}_L$ is the integral-scale Weber number
$\mathit {We}_L = \rho _l u_L^{2} L / \sigma$. The initial conditions employed are similar to those adopted by Chen et al. (Reference Chen, Kharif, Zaleski and Li1999), Iafrati (Reference Iafrati2009, Reference Iafrati2011) and Wang et al. (Reference Wang, Yang and Stern2016): the initial dimensionless wave-surface height
$\eta$ was initialized in terms of the non-dimensional streamwise coordinate
$x$ using
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn12.png?pub-status=live)
where $S_1=a_1k_1=0.55$ is the slope of the fundamental wave component,
$a_1$ and
$k_1 = 2{\rm \pi} /L$ are, respectively, the corresponding dimensional amplitude and wavenumber and the wave propagates in the
$x$ direction. Note that the lengths in the expression for
$\eta$ have been non-dimensionalized by the wavelength of the fundamental wave component,
$L$. A list of characteristic scales used for non-dimensionalization is provided in table 1.
Table 1. Characteristic scales used for non-dimensionalization. The characteristic time scale was selected to be equal to that used by Wang et al. (Reference Wang, Yang and Stern2016). With this time scale, the first wave-surface impact after the wave overturns occurs shortly after $t = 1$, while a single wave period corresponds to 2.5 characteristic times. The energies are defined in § 3.3.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_tab1.png?pub-status=live)
In order to generate an ensemble of statistically independent but similar realizations (i.e. with the same configuration), the free surface was further perturbed by a set of random displacements $\Delta \eta = \Delta \eta (z)$ smaller than the local grid spacing, such that every interfacial mesh node with the same
$z$ (spanwise) coordinate within the same realization has the same
$\Delta \eta$. While the shift
$\Delta \eta$ is not explicitly resolved by the mesh, it perturbs the volume fraction in these interfacial nodes and provides an implicit disturbance to the original interface. This choice of perturbation preserves the modal content of the wave profile in the streamwise direction since the mesh employed in this work is Cartesian. In other words, the initial conditions are effectively two-dimensional except for a perturbation of spanwise modes. As in the works cited above, the air above the wave surface was initialized at rest, while the water below the surface was initialized with the following dimensionless velocity field (Iafrati Reference Iafrati2009, Reference Iafrati2011):
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn13.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn14.png?pub-status=live)
where $y$ is the non-dimensional coordinate antiparallel to gravity. Periodic boundary conditions were employed in the
$x$ and
$z$ directions, while free-slip boundary conditions were employed on the two remaining boundary faces of the computational domain, which is a cube of length
$L$. In each of the realizations in this ensemble, the computational mesh consists of about 4.2 million mesh nodes with a non-dimensional minimum grid spacing of
$1/216$. This is equivalent to a dimensional grid spacing of 1.25 mm for a water wave in air at atmospheric conditions with the aforementioned dimensionless parameters. To put this in context, the dimensional Hinze scale (2.2) takes a value of the order of
$3$ mm. As evidenced in figures 3, 4, 6 and 7, mesh insensitivity is observed in the energetics and bubble statistics at this mesh resolution. For the latter, this insensitivity is observed over a subrange of super-Hinze-scale bubble sizes where turbulent break-up is expected to be dominant. (See also the description of a more resolved ensemble used in this mesh sensitivity study in the ensuing paragraph.) This resolution was selected to permit more ensemble realizations and enable statistical convergence in a time-resolved sense. The mesh is non-uniform and is finer closer to the central region of the domain, such that a large majority of the generated bubbles are resolved using the minimum grid spacing. Snapshots of the initial and post-breaking waveforms, with the computational mesh overlaid, are illustrated for one of the realizations in figure 1. The non-dimensional time step adopted is
$\Delta t = 6.0 \times 10^{-5}$, leading to a maximum Courant number of about 0.1 throughout the computational domain in the course of the simulations. The relatively small Courant number reduces the shape mismatch between the streak tube and the numerical flux polyhedron in the volume-of-fluid-based advection scheme (Ivey & Moin Reference Ivey and Moin2017), which is important to manage in the case of inadequately resolved mixed-phase regions. The evolution of the waveform for one of these realizations is depicted from two different viewpoints in figure 2 to illustrate the interfacial features generated as the wave breaks.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig1.png?pub-status=live)
Figure 1. Snapshots of the spanwise cross-section of the $\phi =0.5$ isosurface of one of the baseline ensemble realizations corresponding to (a) the initial waveform (
$t = 0$) and (b) the wave sometime after breaking has occurred (
$t=4.03$). The grey lines in the background depict the local mesh configuration. The nodal mesh volumes in the uniform-resolution regions are isotropic.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig2.png?pub-status=live)
Figure 2. Snapshots of (a,c,e,g) an axonometric projection from above the wave and (b,d,f,h) the spanwise cross-section of the $\phi =0.5$ isosurface of one of the baseline ensemble realizations corresponding to various times in the wave-breaking process: (a,b)
$t=2.00$, (c,d)
$t=3.01$, (e,f)
$t=4.03$ and (g,h)
$t=5.04$. The snapshots are sampled at an interval of 1.01 characteristic times. The waves are travelling from left to right, and wrap around the domain due to the periodic boundary conditions in the streamwise direction.
Besides the baseline case, an ensemble of three numerical simulations with the same wave parameters and a higher mesh resolution was also generated for a mesh convergence study of the specifically desired quantities studied in this work, such as the bubble-size distribution. In this ensemble, the computational mesh consists of about 32 million mesh nodes with a non-dimensional minimum grid spacing of $1/432$ (dimensionally equivalent to 0.63 mm) and a non-dimensional time step of
$1.2 \times 10^{-5}$. A more detailed rendering of one of the realizations from this ensemble may be found in the video discussed by Chan et al. (Reference Chan, Mirjalili, Jain, Urzay, Mani and Moin2019). Note that even with this higher spatial resolution, the Kolmogorov length scale (2.1),
$L_{K} \sim 30\ \mathrm {\mu }\text {m}$, remains inaccessible as in many previous breaking-wave simulations. In fact, it may be shown that the smallest scales of turbulent motion near a phase interface could require sub-Kolmogorov resolution (Dodd & Jofre Reference Dodd and Jofre2019) that is rarely attained. As remarked in the introduction, an ensemble of direct numerical simulations of breaking waves that resolves these dynamics with converged statistics remains a formidable undertaking due to the inherent scale separation in these oceanic systems, and in energetic multiphase flows in general. In view of these limitations, only intermediate-scale quantities where sub-Hinze-scale dynamics has a limited influence are discussed in this work, including the bubble statistics introduced in § 2. In other words, the formation of sub-Hinze-scale bubbles, which may be due to mechanisms distinct from a turbulent break-up cascade, is beyond the scope of this work.
3.3. Time evolution of the breaking wave
From the snapshots of the plunging wave-breaking process in figure 2, it is evident that the wave profile and accompanying distribution of bubbles dramatically evolve in the span of about 1 to 2 wave periods, or 3 to 5 characteristic times. A qualitatively similar evolution was observed in the studies by Bonmarin (Reference Bonmarin1989), Chen et al. (Reference Chen, Kharif, Zaleski and Li1999), Deane & Stokes (Reference Deane and Stokes2002), Blenkinsopp & Chaplin (Reference Blenkinsopp and Chaplin2007), Blenkinsopp & Chaplin (Reference Blenkinsopp and Chaplin2010), Iafrati (Reference Iafrati2009), Rojas & Loewen (Reference Rojas and Loewen2010), Kiger & Duncan (Reference Kiger and Duncan2012), Deike, Popinet & Melville (Reference Deike, Popinet and Melville2015), Deike et al. (Reference Deike, Melville and Popinet2016), Lim et al. (Reference Lim, Chang, Huang and Na2015) and Wang et al. (Reference Wang, Yang and Stern2016). In particular, a number of tubular air cavities are formed, and then subsequently deformed and ruptured beneath multiple surface splash-ups. This bubble fragmentation eventually ceases as large bubbles degas more quickly, leaving a plume of smaller, slowly rising bubbles. More details are provided by Chan (Reference Chan2020).
This evolution in the wave dynamics is also manifested in the energetics of the wave, whose time evolution averaged over the wave ensemble is plotted in figure 3. The volume-integrated total energy plotted in figure 3(a) is defined as follows:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn15.png?pub-status=live)
where the domain of integration $\varOmega _d$ is taken here to be the entire computational domain with volume
$\mathcal {V}_d = L^{3}$. In this domain, the range of
$\hat {y}$ is defined as
$\hat {y} \in [-L/2,L/2]$. The energy is computed with respect to the reference state of a quiescent air–water system with a flat interface where the water phase occupies the bottom half of the computational domain
$(y \in [-1/2,0))$ and the air phase occupies the top half
$(y \in (0,1/2])$. It is normalized by the energy
$\hat {E}_n$ of a body of water occupying half the volume of the domain and travelling uniformly at the characteristic speed
$u_L = (gL)^{1/2}/(2{\rm \pi} )^{1/2}$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig3.png?pub-status=live)
Figure 3. (a) Non-dimensional ensemble-averaged total energy ($E_t$) of the waves in the baseline ensemble, as well as the non-dimensional total energy of one of the realizations in the ensemble with increased mesh resolution, as a function of non-dimensional time. (b) Non-dimensional rate of change of the energies in (a), computed by differencing the total energy every 500 time steps in the baseline ensemble and every 2500 time steps in the more resolved ensemble. The horizontal dotted line in (b) denotes the nominal rate of change of energy if all the energy initially present in the wave was to be dissipated in one wave period. The left-hand shaded region spans a third of a wave period between
$t=2.30$ and
$t=3.14$, while the right-hand shaded region spans a fifth of a wave period between
$t=3.45$ and
$t=3.95$. In both panels, the lines corresponding to the baseline ensemble denote ensemble-averaged quantities, while the error bars span, in each direction, twice the standard error of the energies of each realization with respect to the ensemble average, thus representing the variation over the ensemble.
The total energy is observed to decay in time in general agreement with the trends observed by Wang et al. (Reference Wang, Yang and Stern2016) and Deike et al. (Reference Deike, Melville and Popinet2016). Figure 3(a) suggests that the total energy is fairly insensitive to the grid spacing, so the baseline mesh should be considered sufficient for the present study. Figure 3(b) generally supports this claim as well, although some differences in the rate of change of the total energy are visible during the active wave-breaking phase due to the difference in the mesh resolution. Note, however, that the standard errors in the baseline ensemble are also noticeably larger during this phase.
Figure 3(b) suggests that during the large-dissipation phase $t \in (2.3,3.1)$, the system loses energy at a rate comparable to that if the initial energy of the wave was dissipated in one wave period. The extraction of the wave period as a defining time scale is compatible with the expressions for the global Kolmogorov and Hinze scales, (2.1) and (2.2), and lends support to the dissipation rate estimate in appendix B. The rate of energy dissipation also appears to be reasonably approximated as quasi-stationary during this time interval, notwithstanding some temporal fluctuations.
The time evolution of the dissipation rate in figure 3(b) exhibits two intervals of interest. In the first interval $t \in (2.3,3.1)$ discussed earlier, the dissipation rate remains quasi-steady as it attains its peak magnitude. In the second interval
$t \in (3.5,4.0)$, the dissipation rate gradually decays in tandem with the surrounding turbulence. The sequence of events depicted in figure 2 illustrates the importance of bubble break-up in both intervals. In § 5, the evolution of various bubble statistics in these two intervals is analysed, and the characteristics of the bubble-mass-transfer dynamics in these intervals are compared and contrasted.
4. Algorithms
Two algorithms were used to retrieve bubble statistics from the simulation ensembles described in § 3 in order to shed light on bubble-mass transfer between large and small bubble sizes during the active wave-breaking phase. Section 4.1 describes an algorithm used to identify the sizes and locations of bubbles in each simulation snapshot. The algorithm was used to determine the bubble-size distribution and the total interfacial area. The algorithm in § 4.2 tracks these bubbles over successive simulation snapshots in order to detect break-up and coalescence events, which drive the time evolution of the bubble-size distribution as discussed in §§ 2.1 and 2.2. More details of these algorithms are discussed by Chan et al. (Reference Chan, Dodd, Johnson, Urzay and Moin2018a) and Chan (Reference Chan2020) and a brief overview is offered in this section.
4.1. Bubble identification algorithm
An existing bubble identification algorithm (Hebert et al. Reference Hebert, Schmidt, Knaus, Phillips and Magari2008; Herrmann Reference Herrmann2010; Tomar et al. Reference Tomar, Fuster, Zaleski and Popinet2010) was refined in this work to increase the accuracy in determining the volumes of the identified bubbles. The basis of the algorithm is the identification of connected regions of computational nodes, where each region corresponds to an individual bubble. After identifying these regions, one may then integrate the gaseous volume fraction, $1-\phi$, over the nodes in each of these regions to obtain the total volume of each bubble. A region is assembled by linking node pairs that each satisfy a grouping criterion. The traditional criterion requires that
$\phi <1$ in each node of a pair, or that each node contains a non-zero amount of air. This criterion may create bubbles of spurious numerical origin because energetic collisions may create small wisps of air, which are sometimes also referred to as flotsam and jetsam, where
$1-\phi \ll 1$ in each of these nodes. These wisps tend to linger in the flow field due to their low buoyancy and slow degassing rate. Wisps may form node pairs that satisfy the criterion and become spuriously linked together. Clipping small values of
$1-\phi$ mitigates this issue at the expense of the volume accuracy of resolved bubbles. Instead, this work uses the additional criterion that at least one node in a pair satisfies
$1-\phi > \phi _{c,m}$, so that nodes with small
$1-\phi$ are linked only if they neighbour nodes with large
$1-\phi$. In this work,
$\phi _{c,m}=0.5$ was selected.
The algorithm yields a list of bubbles in every simulation snapshot, and can simultaneously yield the volume and centre-of-mass (centroid) location of each identified bubble, among other statistics. Accuracy of the determined volume and centroid location is crucial to the performance of the bubble tracking algorithm to be discussed in § 4.2. Volume accuracy also impacts the trends in the bubble-size distribution to be discussed in § 5.1. A comparison of the errors incurred from various grouping criteria has been detailed in other works, such as Chan et al. (Reference Chan, Dodd, Johnson, Urzay and Moin2018a) and Chan (Reference Chan2020). This involved systematically quantifying the incurred error for a number of simple test cases, as summarized in appendix C.
4.2. An event detection algorithm to track bubbles
A bubble tracking algorithm was developed in this work to detect break-up and coalescence events by maintaining a tally of all bubbles over time and tracing the lineage of each bubble. To construct these lineages, lists of bubbles with their sizes and locations are compared between successive simulation snapshots. These snapshots need not arise from consecutive time steps. This is because the principle of mass conservation and the Courant–Friedrichs–Lewy (CFL) condition are used to determine if a bubble has simply been advected between the two snapshots without any change in topology, or if a break-up or coalescence event has occurred. These two constraints are discussed in appendix D in relation to the errors in the identification algorithm discussed in appendix C. As discussed in § 2.1, disconnections and reconnections with the atmosphere, which is treated as a large gaseous reservoir, are not included in tallies of break-up and coalescence events, respectively. This is to prevent frequent topology changes involving the atmosphere and large, non-spherical bubbles with convoluted shapes near the wave surface from obscuring the remaining break-up and coalescence statistics.
The time interval between successive snapshots is now discussed in relation to other characteristic time scales for the baseline ensemble introduced in § 3. The non-dimensional snapshot interval for a set of snapshots available for all 30 realizations is $\Delta t_{s,30} = 3.0\times 10^{-2}$. A set of more closely spaced snapshots is available for 20 of these realizations with a time separation of
$\Delta t_{s,20} = 6.0\times 10^{-3}$. From (2.2) and the inertial subrange scaling
$u_{L_n} \sim (\bar {\varepsilon } L_n)^{1/3}$, one may estimate the characteristic mean time interval between break-up events for Hinze-scale bubbles to be
$\Delta t_{H} \sim 10^{-1}$. This scaling for
$u_{L_n}$ also suggests that the corresponding time interval for super-Hinze-scale bubbles will exceed
$\Delta t_{H}$. Since the snapshot interval is shorter than the mean super-Hinze-scale break-up time, the resulting statistics should provide a realistic picture of at least super-Hinze-scale turbulent break-up. In particular, there are at least
$O(10)$ snapshots in the 20-realization dataset between two super-Hinze-scale break-up events on average.
Two additional remarks are in order here. First, a snapshot interval that exceeds the simulation time step prevents the algorithm from identifying every break-up event of interest. Second, the discussion in appendix C suggests that the identification algorithm in § 4.1 cannot effectively discern the separation between two bubbles spaced less than a grid cell apart. Bubbles that break up but do not separate quickly enough may then register repeated break-up and coalescence events. An excessively small snapshot interval, such as one that is much smaller than the integral time scale, may capture some of these confounding events. This limitation has also been observed in other event detection algorithms (Rubel & Owkes Reference Rubel and Owkes2019) and appears to be inherent in event detection in turbulent flows. An ideal snapshot interval should resolve the characteristic mean break-up times of most bubbles while being sufficiently long to skip over these confounding events. The impacts of these algorithmic limitations are discussed in § 5.2.
5. Bubble statistics
5.1. The bubble-size distribution
$\bar {f}$
Figure 4 plots the bubble-size distribution $\bar {f}$ at selected time instances for the ensembles introduced in § 3. The bubble size,
$D_i = [(6\mathcal {V}_i)/{\rm \pi} ]^{1/3}$, is the diameter of a sphere with the same volume,
$\mathcal {V}_i$. The size distribution has dimensions
$(\text {length})^{-4}$, is non-dimensionalized by the wavelength (
$L^{4}$), and was computed by histogram binning with
$N_{bin}$ bins of equal logarithmic spacing, where the smallest bin is two orders of magnitude larger than the diameter error in appendix C. The non-dimensional discrete distribution satisfies
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn16.png?pub-status=live)
where $\bar {f}(D_j/L;t)$ includes bubbles of sizes between
$D_j$ and
$D_{j+1}$, and
$\varDelta (D_j/L) = (D_{j+1}-D_j)/L$. The distributions from the two ensembles reasonably overlap at intermediate bubble sizes, suggesting that the distribution is fairly mesh-insensitive at these sizes. While the large-bubble distribution of the more resolved ensemble has significant standard errors due to the small ensemble size and the small number of large bubbles, these large bubbles are well resolved in both ensembles.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig4.png?pub-status=live)
Figure 4. The non-dimensional bubble-size distribution $\bar {f}(D_j/L;t)$ against non-dimensional size
$D/L$ for the baseline and more resolved ensembles at (a)
$t=2.50$, (b)
$t=2.74$, (c)
$t=3.71$ and (d)
$t=3.95$. The dotted vertical lines in the same colour/shade as the respective distributions indicate the mesh resolution of the respective ensembles. The dashed sloped lines indicate the
$D^{-10/3}$ power-law scaling for an idealized turbulent bubbly flow. The error bars span, in each direction, twice the standard error of the distribution in each realization with respect to the ensemble average. The shaded regions represent the bubble-size subrange over which the power-law fit exponents in figures 5 and 6 were obtained. More snapshots in time of the bubble-size distribution are provided by Chan (Reference Chan2020).
One may link the evolution of the size distribution in figure 4 to the wave-breaking process in figure 2. As the cylindrical air cavities deform and rupture between $t=2$ and
$t=3$ (figure 2a–d), the size distribution displays a power-law trend (figure 4a,b) with an exponent near the
$-10/3$ value postulated by Garrett et al. (Reference Garrett, Li and Farmer2000), and observed by Deane & Stokes (Reference Deane and Stokes2002) and Deike et al. (Reference Deike, Melville and Popinet2016), suggesting that the turbulent break-up cascade examined in Part 1 dominates the bubble-mass transfer during these early wave-breaking stages. As evidenced in figure 3(b), the dissipation rate is also large. In addition, the smallest resolvable bubbles rapidly appear, as also observed by Deane & Stokes (Reference Deane and Stokes2002) and revisited in § 5.3. While bubbles continue to deform, fragment and interact between
$t=3$ and
$t=4$ (figure 2c–f), the power-law scaling becomes shallower at intermediate sizes (figure 4c,d). As evidenced in figure 3(b), the dissipation rate also begins to decay. Figure 2(e–h) suggests that large-scale air entrainment has mostly ceased by this time, even as intermediate-scale bubble break-up continues, as revisited in § 5.4.
This work focuses on the statistics of intermediate-sized (super-Hinze-scale) bubbles with non-dimensional diameters larger than $1.1\times 10^{-2}$, or dimensional diameters larger than
$2.9$ mm, for which mesh insensitivity was observed. Figure 5 quantifies the evolution of the exponent of a power-law fit to the size distribution in figure 4 over an intermediate size subrange. Observe its oscillatory variation around
$-10/3$ at early times, which may be attributed to regularity in the successive wave-surface splash-ups and impingements. This variation suggests that the largest scales may be influencing these intermediate-size statistics through analogous fluctuations in the energy and bubble-mass cascade rates. Indeed, fluctuations manifest in the dissipation rate in figure 3(b) and the bubble-mass-transfer rate,
$\overline{W_b}$, in § 5.4. This variation may also reflect the finite convergence time of the break-up dynamics towards quasi-equilibrium (Filippov Reference Filippov1961; Qi, Masuk & Ni Reference Qi, Masuk and Ni2020). Figure 6(a) depicts the size distribution averaged over these early time instances with two power-law fits, which are in general agreement with
$D^{-10/3}$, although the fit exponent exhibits slight sensitivity to the size subrange used for fitting. A similar agreement was obtained in the ensemble-averaged and time-averaged statistics of Deane & Stokes (Reference Deane and Stokes2002) and Deike et al. (Reference Deike, Melville and Popinet2016), and the time-averaged statistics of Wang et al. (Reference Wang, Yang and Stern2016), which builds confidence in the results of this work. Ensemble averaging improves statistical convergence and reduces data scatter such as that observed by Wang et al. (Reference Wang, Yang and Stern2016) in their instantaneous distributions. While the instantaneous distribution in this work oscillates between
$D^{-3}$ and
$D^{-4}$, better agreement with
$D^{-10/3}$ is obtained via time averaging over the entire interval, consistent with the suggestion of Deike et al. (Reference Deike, Melville and Popinet2016) that the scaling of the time-averaged distribution,
$\bar {f}^{T}$, is sensitive to the time-averaging window, as is revisited in § 5.4. A more rigorous test of the
$D^{-10/3}$ scaling is shown in figure 7(a) by premultiplying the time-averaged distribution and using a linear scale on the vertical axis. Time averaging is further explored in appendix A.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig5.png?pub-status=live)
Figure 5. The variation of the exponent of a power-law fit over a subset of the baseline distribution, $\bar {f}$, in figure 4 with non-dimensional time. The fit was performed over bubbles with non-dimensional diameters between
$1.07\times 10^{-2}$ and
$2.23\times 10^{-2}$. This bubble-size subrange is marked in the shaded regions in figure 4. The error bars denote twice the standard error in the fit exponent over the baseline ensemble. The shaded regions span the same time intervals as the corresponding regions in figure 3(b).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig6.png?pub-status=live)
Figure 6. The non-dimensional time-averaged size distribution $\bar {f}^{T}(D_j/L)$ against non-dimensional size
$D/L$ for both ensembles. The time-averaging interval is (a) between
$t=2.30$ and
$t=3.14$, corresponding to the left-hand shaded regions in figures 3(b) and 5, and (b) between
$t=3.45$ and
$t=3.95$, corresponding to the right-hand shaded regions in the same figures. For a description of the vertical lines, dashed sloped line in (a) and error bars, refer to the caption of figure 4. The bubble-size subrange over which the bottom power-law fit in (a) and the power-law fit in (b) were performed is the same subrange highlighted in figure 4 and used in figure 5. The subrange for the top power-law fit in (a) is
$[1.69\times 10^{-2},3.54\times 10^{-2})$. The fit exponent uncertainties denote twice the standard error over the baseline ensemble.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig7.png?pub-status=live)
Figure 7. The compensated time-averaged size distribution using the same time-averaging intervals as in figure 6: (a) $\bar {f}^{T}(D_j/L) D_j^{10/3}$ from the first interval, premultiplied by the inverse of the
$D^{-10/3}$ scaling, and (b)
$\bar {f}^{T}(D_j/L) D_j^{8/3}$ from the second interval, premultiplied by the inverse of the observed
$D^{-8/3}$ scaling, against non-dimensional size
$D/L$ for both simulation ensembles. The compensated distribution is normalized by its value at
$D/L = 1.07\times 10^{-2}$. For a description of the vertical lines and error bars, refer to the caption of figure 4.
Returning to figure 5, the fit exponent deviates from $-10/3$ at later times, about one wave period after the onset of breaking. Interestingly, the exponent does not oscillate, in seeming correlation with the cessation of large-scale entrainment. Figure 6(b) depicts the size distribution averaged over these late time instances. Its increased magnitude at most resolved sizes relative to figure 6(a) (by about 3–5 times; see also figure 15) suggests that there are more small and intermediate-sized bubbles at later times, and the break-up flux from larger sizes consistently outweighs mass loss from degassing. Two distinct power-law scalings emerge in two size subranges in agreement with the measurements of Tavakolinejad (Reference Tavakolinejad2010) and Masnadi et al. (Reference Masnadi, Erinin, Washuta, Nasiri, Balaras and Duncan2020), which builds further confidence in the results of this work (see also Castro, Li & Carrica Reference Castro, Li and Carrica2016). The size distribution is reasonably described by a
$D^{-8/3}$ scaling at intermediate sizes, which is shallower than
$D^{-10/3}$. This is supported by the compensated time-averaged distribution in figure 7(b). At large sizes, the power-law scaling is steeper than
$D^{-8/3}$ by a factor of
$D^{-2}$ in agreement with the increased importance of buoyant degassing postulated by Garrett et al. (Reference Garrett, Li and Farmer2000) and observed by Deike et al. (Reference Deike, Melville and Popinet2016), as well as in figure 2(e–h). For comparison, Tavakolinejad (Reference Tavakolinejad2010) obtained an exponent between
$-2.85$ and
$-2.58$ at small sizes, and between
$-6.57$ and
$-3.85$ at large sizes.
Figures 3(b) and 4–7 reflect the unsteadiness in the wave dynamics and demonstrate the emergence of two time intervals with distinct characteristics, suggesting the presence of distinct bubble generation and evolution mechanisms. In the first interval, the dissipation rate appears to be quasi-steady, and the $D^{-10/3}$ power-law scaling in the time-averaged size distribution supports the presence of a quasi-steady bubble break-up cascade. In the second interval, the dissipation rate is seen to decay, and the time-averaged size distribution deviates from
$D^{-10/3}$. Yet, the robustness of an alternative power-law scaling with a negative exponent is indicative of size-invariant cascade-like behaviour, as Part 1 suggested that size-local break-up may still occur in the absence of the
$D^{-10/3}$ scaling. These observations motivate a closer look at other bubble statistics during these two intervals to test the theoretical analysis developed in Part 1 for the break-up cascade in the first interval, and to extend this analysis to the possible cascade-like behaviour in the second interval. Note that these conclusions were drawn from statistics over a limited size range. Future ensembles with larger scale separation will be valuable for shedding more light on these conclusions. Volume averaging also precludes the reporting of spatial size distributions. The physical parameters and grid size selected in this work reflect the trade-off between the resolved size range and the ensemble size. As noted in §§ 1 and 3.2, a larger ensemble increases statistical convergence in a statistically unsteady flow that does not strictly permit time averaging. However, it also necessitates a smaller resolved range from a practical perspective. The mesh insensitivity in figures 4, 6 and 7, as well as the agreement of their power-law scalings with prior works, supports these choices in light of this trade-off.
5.2. The differential break-up rate
$\overline{g_b f}$
In order to gain more insights into the evolution of the size distribution, $\bar {f}$, as described by (2.4) and (2.6), the differential break-up rate
$\overline {g_b f}$ is now analysed. Break-up events were detected in the baseline ensemble using the algorithm in § 4.2. The ensemble-averaged event rate per unit domain volume and unit size is averaged over the non-dimensional sampling interval
$16 \Delta t_{s,20} = 9.6\times 10^{-2}$ to converge statistics over more events. This interval is kept small enough that temporal variations are not significantly smoothed. The number of histogram bins for
$\overline {g_b f}$ is half that for
$\bar {f}$ to increase the bin width and improve this convergence. To ensure that the snapshot interval discussed in § 4.2 has been reasonably selected, rates obtained using two intervals,
$\Delta t_{s,20}$ and 2
$\Delta t_{s,20}$, are reported to demonstrate snapshot convergence. Since the variation of
$\overline {g_b f}$ with parent bubble size is of interest,
$\overline {g_b f}$ is normalized by its maximum value
$(\widetilde {\overline {g_b f}})$ at every time instance to highlight this convergence. This normalization is also carried out in § 5.4.
Figure 8 shows $\overline {g_b f}$ near the time instances for which
$\bar {f}$ was plotted in figure 4. Signatures of the
$D^{-4}$ scaling derived for a quasi-steady turbulent break-up cascade in Part 1 emerge at intermediate sizes. Note that the computation of
$\overline {g_b f}$ samples energetic regions more frequently since break-up only occurs with sufficient energy to change the bubble topology. Thus,
$\overline {g_b f}$ is susceptible to large-scale inhomogeneities, and is difficult to converge since the scale separation in the simulated waves is not exceedingly large.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig8.png?pub-status=live)
Figure 8. The normalized differential break-up rate $\widetilde {\overline {g_b f}} (D_j/L;t)$ averaged over the non-dimensional sampling interval
$16 \Delta t_{s,20} = 9.6\times 10^{-2}$ plotted against non-dimensional size
$D/L$ for a 20-realization baseline ensemble subset at (a)
$t=2.51$, (b)
$t=2.75$, (c)
$t=3.71$ and (d)
$t=3.95$. The rate is normalized such that the maximum value in each plot is 1, i.e.
$\widetilde {\overline {g_b f}} = \overline {g_b f}/\overline {g_b f}|_{max}$. The selected bin sizes are twice those in figure 4 in logarithmic space. For a description of the vertical lines and error bars, refer to the caption of figure 4. The dashed sloped lines indicate the theoretical
$D^{-4}$ power-law scaling for an idealized turbulent bubbly flow. The shaded regions represent the bubble-size subrange over which the power-law fit exponents in figures 9 and 10 were obtained. More snapshots in time are provided by Chan (Reference Chan2020).
Figure 9 quantifies the evolution of the exponent of a power-law fit on the differential rate in figure 8 over a subrange of intermediate sizes, relative to the $-4$ exponent derived in Part 1 using traditional turbulent-flow scalings. The behaviour of the exponent in figure 9 for
$\overline {g_b f}$ loosely tracks the corresponding behaviour in figure 5 for
$\bar {f}$: in the early wave-breaking stages, the exponent oscillates around the exponent discussed in Part 1; in the late stages, it does not strongly oscillate. As with the size distribution in figure 6, the differential rate may be time-averaged over these two intervals, as depicted in figure 10. In the first interval, it may be described by a
$D^{-4.1\pm 0.3}$ power-law fit, in agreement with
$D^{-4}$. In the second interval, it may be described by a
$D^{-3.3\pm 0.4}$ fit, which is about a factor of
$D^{-2/3}$ steeper than the
$D^{-8/3}$ scaling in the corresponding size distribution, although vestiges of the
$D^{-4}$ scaling remain in the differential rate at larger sizes. These power-law scalings at intermediate sizes are largely consistent with the turbulent scaling of the break-up frequency,
$g_b \sim D^{-2/3}$, in agreement with the references cited in § 1 that were used in the theoretical analysis of Part 1.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig9.png?pub-status=live)
Figure 9. The variation of the exponent of a power-law fit over a subset of the differential break-up rate, $\overline {g_b f}$, in figure 8 with non-dimensional time. The fit was performed on the data with a snapshot interval of
$\Delta t_{s,20}$. For a description of the bubble-size subrange over which the fit was performed, error bars and shaded regions, refer to the caption of figure 5.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig10.png?pub-status=live)
Figure 10. The normalized time-averaged differential break-up rate $\widetilde {\overline {g_b f}}^{T} (D_j/L)$ against non-dimensional size
$D/L$ for the 20-realization baseline ensemble subset. Panels (a,b) correspond to the time-averaging intervals in figure 6. For a description of the vertical and sloped lines, error bars, normalization, bubble-size subrange and dataset over which the fits were performed, and fit-exponent uncertainties, refer to the captions of figures 4, 6, 8 and 9.
In summary, the trends observed in $\overline {g_b f}$ mostly echo those for
$\bar {f}$ in § 5.1. The unusual persistence of the
$D^{-4}$ scaling may be suggestive of a tendency of the bubble-mass flux to remain quasi-local and quasi-self-similar, since
$\overline{W_b} \sim \overline {g_b f} D^{4}$ in this limit as shown in Part 1. These characteristics of the bubble-mass flux are further addressed in § 5.4. Analogous to § 5.1, future ensembles of higher
$\mathit {Re}_L$ and
$\mathit {We}_L$ simulations could bring about more clarity on the robustness of these scalings.
5.3. The break-up probability distribution
$\check {q}_b$
In addition to $\overline {g_b f}$, the probability distribution of child bubble sizes,
$\check {q}_b$, also drives the evolution of
$\bar {f}$, as described by (2.4) and (2.6), and may also be computed by the algorithm in § 4.2. The distribution
$\check {q}_b(D_c^{3};t|D_p^{3})$ is symmetric in bubble-volume space and is presented in this space for a more intuitive interpretation. It satisfies
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn17.png?pub-status=live)
The data for each parent bubble size $D_{p,j}$ are averaged over two histogram bins
$[D_{p,j},D_{p,j+2})$, where the bins are identical to those in § 5.2. The reported
$\check {q}_b$ is time-averaged over the two time intervals identified in §§ 3.3, 5.1 and 5.2.
Figure 11 plots the time-averaged break-up probability distribution, $\overline{\check{q}_b}^{T}$, for three distinct parent bubble sizes and the two aforementioned time intervals. Each row corresponds to a single parent bubble size, while each column corresponds to a single time interval. Some child bubble sizes are smaller than the mesh resolution and have to be excluded. Under the influence of this exclusion,
$\overline{\check{q}_b}^{T}$ falls off towards large and small child bubble sizes for small parent bubbles.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig11.png?pub-status=live)
Figure 11. The time-averaged break-up probability $\overline{\check {q}_b}^{T}(D_{c,j}^{3}|D_p^{3})$ against normalized child bubble volume
$D_c^{3}/D_p^{3}$ for the 20-realization baseline ensemble subset during the (a,c,e) early (
$t=2.30$–
$3.14$) and (b,d,f) late (
$t=3.45$–
$3.95$) wave-breaking stages. The non-dimensional parent bubble sizes considered are (a,b) between
$1.86 \times 10^{-2}$
$(\mathit {We}_{D_p} = 18)$ and
$2.23 \times 10^{-2}$
$(\mathit {We}_{D_p} = 25)$, (c,d) between
$2.68 \times 10^{-2}$
$(\mathit {We}_{D_p} = 33)$ and
$3.23 \times 10^{-2}$
$(\mathit {We}_{D_p} = 45)$ and (e,f) between
$3.54 \times 10^{-2}$
$(\mathit {We}_{D_p} = 53)$ and
$4.25 \times 10^{-2}$
$(\mathit {We}_{D_p} = 72)$. For the definition of
$\mathit {We}_{D_p}$, refer to (2.4) in Part 1. Only child bubbles of radii larger than the mesh resolution are considered. The vertical dashed lines demarcate the child bubble volumes in the break-up event where one of the child bubble radii corresponds to this resolution limit. The error bars denote one standard error. The time intervals in the legend denote the snapshot intervals in the detection algorithm. Fits to the beta distribution have been added in (d,e). These fits were performed on the data obtained by taking the snapshot interval to be
$\Delta t_{s,20}$, which was also used to plot the histogram bars. The two numbers in the corresponding legend entries are the shape parameters characterizing each fit, which were obtained via the method of moments.
In general, the child bubble volume distribution is relatively close to uniform for events within the intermediate bubble-size subrange where power-law fits were earlier obtained for $\bar {f}$ and
$\overline {g_b f}$ in §§ 5.1 and 5.2. Hence, the uniform distribution, which was analysed in Part 1, appears to be a suitable, albeit rudimentary, surrogate model for
$\overline{\check {q}_b}^{T}$ in most cases. One exception occurs at large parent bubble sizes outside the aforementioned size subrange in the early wave-breaking stages. Here, large-size-ratio events involving one large child bubble and one small child bubble do frequently occur as evidenced in figure 11(e), bypassing the break-up cascade and generating small bubbles as observed in § 5.1. These individual non-local contributions to bubble-mass transfer are small in volume and may not necessarily influence the locality of the break-up flux,
$\overline{W_b}$, as noted in appendix B of Part 1. Their relative influence is analysed in more detail in § 5.4.1. Another exception occurs in the late wave-breaking stages in figure 11(d). In these two cases, the break-up probability is reasonably described by a beta-distribution fit, which was also discussed in Part 1 as a potential surrogate model. These exceptions indicate that it does not seem appropriate to characterize a turbulent bubbly flow with a single distribution shape without giving more consideration to the bubble-size subrange or flow stage of interest, as many previous studies have done. For example, experimental studies typically report a single distribution of child bubble volumes, agglomerating data from different parent bubble sizes (see figure 1 of Qi et al. (Reference Qi, Masuk and Ni2020), and references therein). This has contributed to a multitude of existing models for the child bubble volume distribution, as noted in Part 1. Interestingly, the results of this work are in qualitative agreement with the model of Lehr, Millies & Mewes (Reference Lehr, Millies and Mewes2002), where ‘equal sized breakage is preferred for small bubbles …as the size of the parent bubbles increases, unequal breakage is preferred’.
5.4. The break-up flux
$\overline{W_b}$
The results of § 5.1 revealed that the bubble-size distribution dramatically evolves between $t=2$ and
$t=4$. Indirect evidence of a quasi-steady bubble break-up cascade (Garrett et al. Reference Garrett, Li and Farmer2000) was observed in the early wave-breaking stages, although multiple mechanisms could explain the observed
$D^{-10/3}$ power-law scaling in the size distribution (e.g. Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019; Yu, Hendrickson & Yue Reference Yu, Hendrickson and Yue2020). The present simulations, with the accompanying detection algorithm for break-up events, also provide detailed break-up statistics, which are analysed in §§ 5.2 and 5.3. These statistics contribute to the break-up flux,
$\overline{W_b}(D;t)$, as discussed in Part 1 and summarized in § 2.1, and provide further indirect evidence of a break-up cascade. The break-up flux across size
$D$ can be directly decomposed into incoming contributions
$\overline{I_p}(D_p;t|D)$ from parent bubbles of sizes
$D_p > D$ and outgoing contributions
$\overline{I_c}(D_c;t|D)$ to child bubbles of sizes
$D_c < D$, as shown in § 2.2, as well as figures 4 and 5 of Part 1. These complementary decompositions respectively provide information on infrared and ultraviolet locality in the flux. Satisfying these measures of locality is a necessary condition for the presence of a break-up cascade, with
$\overline{I_p}$ and
$\overline{I_c}$ decaying rapidly away from
$D$. Cascade-like behaviour may also be directly characterized by size invariance and quasi-stationarity of the break-up flux at intermediate sizes. These observations are direct consequences of the theoretical analysis performed in Part 1. In this subsection, the simulation results are further leveraged by examining these quantities to obtain more direct observations of cascade-like behaviour.
5.4.1. Infrared and ultraviolet locality
Figure 12 plots the time-averaged differential incoming $(\text {parent}, \overline{I_p}^{T})$ and outgoing
$(\text {child}, \overline{I_c}^{T})$ contributions to the break-up flux for the two time intervals identified in §§ 3.3, 5.1 and 5.2, and an intermediate cutoff bubble size near
$D/L = 2\times 10^{-2}$. This choice of
$D/L$ is representative of a bubble size within the intermediate-size subrange in figures 4–10 in which the break-up cascade scalings were recovered during the first time interval of interest. Other choices of
$D/L$, which are not shown here, display a similar behaviour. The differential incoming contributions from parent bubbles strongly decay with increasing parent bubble size, while the differential outgoing contributions to child bubbles strongly decay with decreasing child bubble size, indicating that the break-up flux is strongly infrared and ultraviolet local, as suggested in Part 1. Power-law fits of these decay rates exhibit exponents near
$-7$ and
$5$ for
$\overline{I_p}^{T}$ and
$\overline{I_c}^{T}$, respectively. The theory in Part 1 predicts these exponents for a uniform child bubble volume distribution,
$q_b$, that remains invariant with parent bubble size. Slight deviations from these idealized values may be accounted for by the mild deviations of
$\overline{\check {q}_b}^{T}$ from uniformity and parent-bubble-size invariance observed in § 5.3. Note that reasonable agreement with the exponents from Part 1 is obtained for both time intervals even though the
$D^{-10/3}$ scaling is absent from the size distribution in the second interval. This suggests that cascade-like behaviour persists in this interval, even as the dissipation rate decays. The recovery of the exponents in the second interval is consistent with
$\overline {g_b f} \sim D^{-4}$, whose persistence in both intervals was noted at the end of § 5.2.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig12.png?pub-status=live)
Figure 12. (a,c) The normalized time-averaged differential outgoing transfer rate $\widetilde {\overline{I_c}}^{T}(D_{c,j}|D)$ to child bubbles as a function of non-dimensional child bubble size
$D_c/L$ and (b,d) the corresponding differential incoming transfer rate
$\widetilde {\overline{I_p}}^{T}(D_{p,j}|D)$ from parent bubbles as a function of non-dimensional parent bubble size
$D_p/L$, for the 20-realization baseline ensemble subset at the non-dimensional cutoff size
$D/L = 1.86\times 10^{-2}$ marked by the dashed vertical line. These statistics are time-averaged over the (a,b) early (
$t=2.30$–
$3.14$) and (c,d) late (
$t=3.45$–
$3.95$) wave-breaking stages. The normalization, histogram bins, shaded regions and dotted vertical lines are identical to those in figure 8. Break-up events involving subgrid child bubbles are excluded. The error bars and fit-exponent uncertainties denote twice the standard error over the ensemble. The power-law fits were performed over bubbles with non-dimensional diameters between
$9.74\times 10^{-3}$ and
$1.69\times 10^{-2}$ for
$\widetilde {\overline{I_c}}^{T}$, and between
$2.04\times 10^{-2}$ and
$3.54\times 10^{-2}$ for
$\widetilde {\overline{I_p}}^{T}$, on the dataset with snapshot interval
$\Delta t_{s,20}$.
One might observe that the decay rate of $\overline{I_p}^{T}$ becomes gentler at very large parent bubble sizes, especially during the first time interval. This is reminiscent of the occurrence of large-size-ratio cascade-bypassing break-up events in large parent bubbles in figure 11(e). Two remarks are in order here. First, the contributions from these events to the break-up flux are orders of magnitude smaller than the transfers from events of highest locality. Second, these contributions arise from intermittent break-up events as evidenced in figure 13, which illustrates the instantaneous differential incoming contributions from parent bubbles at two time instances during the first interval. Notwithstanding these caveats, these results suggest that the break-up dynamics is well approximated as local in bubble-size space in both time intervals of interest.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig13.png?pub-status=live)
Figure 13. The normalized differential incoming transfer rate $\widetilde {\overline{I_p}}(D_{p,j};t|D)$ from parent bubbles as a function of non-dimensional parent bubble size
$D_p/L$ for the 20-realization baseline ensemble subset at the non-dimensional cutoff size
$D/L = 1.86\times 10^{-2}$, and at the times (a)
$t=2.51$ and (b)
$t=2.99$. These statistics are averaged over the non-dimensional sampling interval
$16 \Delta t_{s,20} = 9.6\times 10^{-2}$ used in figure 8. For a description of the normalization, shaded regions, vertical lines, excluded events, error bars, bubble-size subrange and dataset for the power-law fits, and the fit-exponent uncertainties, refer to the caption of figure 12.
5.4.2. Size variation and time evolution of the break-up flux
Figure 14 plots the variation of the time-averaged break-up flux $\overline{W_b}^{T}(D)$ with cutoff bubble size
$D$ for the two time intervals identified in §§ 3.3, 5.1 and 5.2. Observe that there is some degree of size invariance in the break-up flux in both time intervals, even though the scale separation in the simulated waves is not exceedingly large as remarked in §§ 5.1 and 5.2. This size invariance is expected in the theory from Part 1 if
$I_p \sim D_p^{-7}$ and
$I_c \sim D_c^{5}$, or more generally if the magnitude of the power-law exponent for
$I_p$ is larger than that for
$I_c$ by 2. These exponents are approximately recovered for both time intervals in § 5.4.1. Note that the size invariance of the break-up flux,
$\overline{W_b}$, is also consistent with approximate stationarity in the size distribution,
$\bar {f}$, as implied by (2.5), as well as the usage of the break-up flux as a proxy for the entrainment flux, as discussed in § 2.1. Further implications of the consistency of quasi-self-similarity and quasi-stationarity in the break-up dynamics are explored in appendix A with reference to the discussion on time averaging in § 5.1. Finally, observe that the ratio of the flux at intermediate sizes to the flux at large sizes differs in the two time intervals of interest. The ratio is close to unity in the first interval, with variations due to large-scale inhomogeneities, but increases by an order of magnitude in the second interval. This suggests that active entrainment at the large scales directly drives bubble-mass transfer from large- to small-bubble sizes in the first interval. The entrainment appears to have ceased in the second interval as mentioned in § 5.1, and bubble-mass transfer seems to be chiefly the result of the fragmentation of already-entrained large cavities and bubbles.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig14.png?pub-status=live)
Figure 14. The normalized time-averaged break-up flux $\widetilde {\overline{W_b}}^{T}(D/L)$ against non-dimensional size
$D/L$ for the 20-realization baseline ensemble subset. These statistics are time-averaged over the (a) early (
$t=2.30$–
$3.14$) and (b) late (
$t=3.45$–
$3.95$) wave-breaking stages. For a description of the normalization, vertical lines, excluded events and error bars, refer to the caption of figure 12. The shaded regions depict the bubble-size subrange over which the size-averaged break-up flux in figure 15 was obtained.
Figure 15 plots the variation of the flux $\overline{W_b}^{D}(t)$ with time, averaged over a narrow intermediate range of bubble sizes enveloping the cutoff size considered in § 5.4.1. Because of the quasi-size-invariance observed in figure 14, the results of figure 15 are representative of the intermediate-size dynamics considered in this work. The flux has distinct characteristics in the two time intervals identified in §§ 3.3, 5.1 and 5.2. The first time interval with a relatively constant rate of energy dissipation appears to be concomitant with an approximately constant flux, closing the loop on the presence of a quasi-steady turbulent bubble break-up cascade in this interval. The oscillations in the flux mirror the oscillations in the dissipation rate in figure 3(b), and appear to be related to the earlier discussion in § 5.1 on the oscillating power-law exponent for the size distribution. The second time interval with a decreasing dissipation rate appears to be associated with a flux that linearly increases with time. As discussed in § 5.4.1, the flux was observed to be strongly size local during this interval, which is indicative of the persistence of cascade-like behaviour. However, the temporally increasing flux does not drive a uniform increase in the magnitude of the size distribution at all sizes, as depicted in figures 4 and 6. Instead, it is observed in § 5.1 that the power-law scaling of the size distribution transitions from
$D^{-10/3}$ to
$D^{-8/3}$. This suggests that small bubbles may be depleted more rapidly at these late times. The continuation of the flux over several characteristic times suggests that the volume of gas present in small bubble sizes gradually increases over time. The assumption in Part 1 that one may neglect this accumulation should eventually break down towards the late wave-breaking stages. The increased importance of coalescence due to this accumulation, and its potential role in the aforementioned depletion of small bubbles, is addressed in appendix E.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig15.png?pub-status=live)
Figure 15. The normalized size-averaged break-up flux $\widetilde {\overline{W_b}}^{D}(t)$ against non-dimensional time
$t$ for the 20-realization baseline ensemble subset. The flux is averaged over the bubble-size subrange spanning non-dimensional diameters between
$1.69\times 10^{-2}$ and
$2.23\times 10^{-2}$, as indicated by the shaded regions in figure 14. For a description of the normalization, excluded events and error bars, refer to the caption of figure 12. The shaded regions span the same time intervals as those in figures 3(b), 5 and 9.
6. Conclusions
This paper investigates the evolution of the bubble-size distribution and other bubble statistics beneath a breaking wave during the active wave-breaking phase using ensembles of numerical simulations. A large ensemble of simulations enhances statistical convergence in these bubble statistics, which is particularly important for the time-dependent behaviour explored in this paper. Two time intervals of interest are identified, with comparison and contrast between them throughout the paper. This analysis goes beyond the size distribution by identifying and characterizing individual break-up events in order to directly observe characteristics associated with the bubble-mass cascade phenomenology that were identified in Part 1. The evolving dynamics in the bubble-mass-transfer process is summarized in table 2. The results of this paper (Part 2) suggest that cascade-like behaviour is present in both time intervals, even when the equilibrium $D^{-10/3}$ power law is not observed in the size distribution at intermediate sizes. They highlight that different bubble generation and evolution mechanisms emerge at different length and time scales during the wave-breaking process, even within the relatively limited range of length scales captured in these simulations. Future ensembles of higher
$\mathit {Re}_L$ and
$\mathit {We}_L$ simulations with resolution of a larger range of scales may bring about more clarity on these mechanisms.
Table 2. Summary of bubble-mass-transfer dynamics in the active wave-breaking phase.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_tab2.png?pub-status=live)
The evolution of the bubble-size distribution and other associated bubble statistics in the present simulation ensembles may be summarized as follows. In the early wave-breaking stages, large air cavities are entrained and successively fragmented. The corresponding rate of bubble-mass transfer from large- to small-bubble sizes is quasi-steady, as is the energy dissipation rate, which also reaches a maximum value in this time interval. These are the ideal characteristics of a forward bubble-mass cascade driven by turbulent fragmentation. As originally proposed by Garrett et al. (Reference Garrett, Li and Farmer2000), this is accompanied by a $D^{-10/3}$ power-law scaling in the size distribution. The break-up statistics of Part 2 provide strong support for the theoretical results of Part 1 in this early time interval. The time-averaged size distribution matches the
$D^{-10/3}$ scaling in agreement with Deane & Stokes (Reference Deane and Stokes2002), Wang et al. (Reference Wang, Yang and Stern2016) and Deike et al. (Reference Deike, Melville and Popinet2016). Deviations of the instantaneous size distribution from this scaling are observed but are reduced by time averaging over the entire interval. The
$D^{-4}$ power-law scaling derived in Part 1 for the differential break-up rate is also recovered. Notwithstanding several intermittent non-local events for large parent bubbles, the probability distribution of child bubble volumes is roughly uniform in the range of parent bubble sizes where the
$D^{-10/3}$ scaling is observed in the time-averaged size distribution. The resulting break-up flux remains strongly infrared and ultraviolet local. The incoming differential flux from parent bubbles approaches a
$D_p^{-7}$ power-law scaling over a range of parent bubble sizes, while the outgoing differential flux to child bubbles approaches a
$D_c^{5}$ power-law scaling over a range of child bubble sizes. These scaling exponents were predicted in Part 1 for a uniform child bubble volume distribution. The observed scalings suggest that the underlying bubble-mass-transfer process is indeed size local and self-similar, although a more robust demonstration of self-similarity may only be evident in larger waves of higher
$\mathit {Re}_L$ and
$\mathit {We}_L$. Together, these observations provide strong support for the presence of a quasi-steady turbulent bubble break-up cascade at intermediate bubble sizes during the early wave-breaking stages.
In the late wave-breaking stages, large entrained bubbles continue to break up into smaller bubbles. The corresponding rate of bubble-mass transfer from large to small sizes accelerates while the dissipation rate begins to decay. Two distinct power-law scalings appear in the size distribution in agreement with Tavakolinejad (Reference Tavakolinejad2010) and Masnadi et al. (Reference Masnadi, Erinin, Washuta, Nasiri, Balaras and Duncan2020) but in deviation from the $D^{-10/3}$ scaling. The increased magnitude of the size distribution suggests that small and intermediate bubbles are more populous in this time interval. Vestiges of locality, stationarity and self-similarity are nonetheless present in the corresponding break-up statistics, indicating that the break-up dynamics still exhibits cascade-like behaviour. The steepening of the size distribution at large bubble sizes is indicative of buoyant degassing. A shallower size distribution appears at intermediate bubble sizes. The data appear to be consistent with a more rapid depletion of small bubbles, potentially via coalescence.
Several concluding remarks are made in view of these observations. The theoretical framework developed in Part 1 enables one to go above and beyond the traditional $D^{-10/3}$ power-law scaling in the size distribution in diagnosing the presence of a break-up cascade. The results of Part 2 demonstrate that the access to detailed spatial information provided by numerical simulations allows a more precise evaluation of various components of the theoretical framework. The multifaceted nature of this analysis enables a more nuanced description of the underlying break-up dynamics even when the
$D^{-10/3}$ scaling is absent. For example, a number of key physical aspects of the cascade phenomenology, such as locality, were determined to be present even in the late wave-breaking stages. The elements of the formalism in Part 1 could thus be seen as an analytical toolkit that may be used to selectively probe the characteristics of break-up and coalescence dynamics in multiscale two-phase flows. This toolkit may be used to advance understanding of cascading processes, and more generally to inform the validation and development of model kernels for the population balance equation.
Knowledge of the physical mechanisms underlying bubble generation enables the quantification of interfacial fluxes and radiation scattering, which are of practical importance in maritime and climate studies, such as carbon sequestration and the persistent wake signatures of seafaring vessels. The rich temporal evolution of the bubble-mass-transfer process observed in this work reflects the complexity in modelling these mechanisms. Elucidation of sustained cascade-like behaviour in the bubble-mass-transfer dynamics presents opportunities for modelling of subgrid bubble break-up in large-eddy simulations.
Acknowledgements
With a deep sense of gratitude, we dedicate this paper and its companion (Part 1) to the memory of Professor Juan Lasheras. These manuscripts build on fundamental concepts introduced in his early work on turbulent two-phase flows. His untimely passing cut short a brilliant career, but his seminal contributions in the fields of fluid mechanics and biomedical engineering will prove a lasting legacy. The authors acknowledge computational resources from the US Department of Energy's INCITE Program. The authors are grateful to A. Mani for fruitful discussions on the bubble-size distribution and to M.S. Dodd for joint work on algorithms to retrieve statistics from numerical simulations that inspired this work.
Funding
This investigation was funded by the Office of Naval Research, Grant no. N00014-15-1-2726, and is also supported by the Advanced Simulation and Computing programme of the US Department of Energy's National Nuclear Security Administration via the PSAAP-II Center at Stanford University, Grant no. DE-NA0002373. W.H.R.C. is also funded by a National Science Scholarship from the Agency for Science, Technology and Research in Singapore.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Averaging operations for computing bubble statistics
The ensemble-averaged bubble-size distribution $f(\boldsymbol {x},D;t)$ was defined in Part 1 at every location
$\boldsymbol {x}$, for every bubble size
$D$ and at some time
$t$ by ensemble averaging
$(\langle \cdot \rangle )$ the number density function of a bubble population, which is constructed by summing a contribution from each bubble in the population
$i=1,\ldots ,N_b(t)$ having a centroid location
$\boldsymbol {x}_i$ and an equivalent size
$D_i$. Then
$f$ may be expressed as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn18.png?pub-status=live)
It was alluded to in § 2 that the volume-averaged and ensemble-averaged size distribution $\bar {f}(D;t)$ is typically reported in breaking-wave simulations and experiments as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn19.png?pub-status=live)
where $\mathcal {V} = \int _\varOmega \mathrm {d} \boldsymbol {x}$ is the volume of the domain
$\varOmega$ over which the bubble population is integrated. Domain
$\varOmega$ should be selected such that it always contains all
$N_b(t)$ bubbles. In this work,
$\mathcal {V}$ is chosen to be the characteristic wave volume,
$L^{3}$. Since this coincides with the computational domain volume, all
$N_b(t)$ bubbles are thereby always included. The volume-averaging operation is analogously defined for other bubble statistics like
$\overline {g_b f}$. The ensemble-averaging and volume-averaging operations commute when
$\mathcal {V}$ and
$\varOmega$ are identical across all the ensemble realizations. Thus,
$\bar {f}$ satisfies the normalization conditions
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn20.png?pub-status=live)
where $n_b(t)$ is the number of bubbles in the system per unit domain volume for an ensemble realization, or the global bubble number density in the realization. The population balance equation for
$f(\boldsymbol {x},D;t) D^{3}$ was written in Part 1 as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn21.png?pub-status=live)
where $v_i$ and
$v_D$ are the velocities of
$fD^{3}$ in the
$\boldsymbol {x}$–
$D$ phase space along the spatial and bubble-size dimensions, respectively, and
$H$ is a model term that includes source, sink and non-local transfer terms for the transport of
$fD^{3}$. Volume averaging the equation over a domain that always includes all
$N_b(t)$ bubbles yields (2.3). A similar procedure may be used to obtain (2.4). Note that the spatial and bubble-size dimensions are orthogonal in the
$\boldsymbol {x}$–
$D$ phase space. Thus, the expected
$D$-scaling of a quantity, such as the theoretical variations of
$f(\boldsymbol {x},D;t)$ and
$\bar {f}(D;t)$ with
$D$, should remain invariant under volume averaging provided there is sufficient statistical convergence in the quantity.
Volume averaging aggregates statistics at different locations, and provides a convenient way of achieving statistical convergence in a limited number of realizations, as well as a more concise description of the underlying dynamics (Hinze Reference Hinze1955). It is argued here that volume averaging also achieves a good approximation for low-order bubble statistics, including the size distribution, as if the underlying flow was statistically homogeneous (Deike et al. Reference Deike, Melville and Popinet2016). More specifically, it averages these statistics over small, localized regions that may each be treated as statistically homogeneous in the spirit of local isotropy. This may be motivated by the arguments behind Kolmogorov's refinements to his original energy cascade hypotheses (Kolmogorov Reference Kolmogorov1962). The original hypotheses (Kolmogorov Reference Kolmogorov1941) are best applied in regions where the assumption of spatial homogeneity is robust. In a statistically inhomogeneous flow, diffusion fluxes across large-eddy length and time scales may become important, thereby manifesting as large-scale variations in the turbulent kinetic energy, as well as intermittency in the localized production and dissipation rates. The refinements of Kolmogorov (Reference Kolmogorov1962) attempt to characterize these variations systematically through an assumption on their statistical distribution. To enable this, volume averaging is performed in small, localized regions that are each approximately statistically homogeneous. The original hypotheses are then recovered for locally volume-averaged flow statistics, with constants of proportionality that account for large-scale inhomogeneities. The original and refined hypotheses do not yield significantly different results when applied to low-order flow statistics (Pope Reference Pope2000, § 6.7.4), such as the second-order velocity structure functions. The global volume-averaging procedure adopted in this work eliminates these diffusion fluxes by averaging flow and bubble statistics over a volume where no net transport in or out is expected. This assumes each small, localized region in the global averaging volume is subject to the same characteristic large-eddy length and time scales, in the spirit of the original hypotheses of Kolmogorov (Reference Kolmogorov1941). This is thus likely to be reasonable for low-order flow and bubble statistics, with additional corrections necessary for higher-order statistics when, for example, $\bar{\varepsilon} ^{n}$ is no longer well approximated by
$(\bar {\varepsilon })^{n}$ for large
$n$. Hence, volume averaging is a physically relevant way of computing low-order bubble statistics in addition to its expedience. Nevertheless, knowledge of the spatial distribution of bubbles may at times be important, and two representative snapshots of the ensemble-averaged and spanwise-averaged liquid volume fraction are provided in appendix F to give a general sense of this spatial variation in the context of breaking waves.
As an endnote, recall the remark in § 5.1 of Part 1 that the analysis for the bubble-mass cascade resembles the monofractal analysis of Eyink (Reference Eyink2005) for the energy cascade. The identification of a single power-law scaling in the bubble-size distribution under a given set of conditions also suggests monofractality in the bubble break-up dynamics (Turcotte Reference Turcotte1986). For example, the $D^{-10/3}$ power-law scaling corresponds to a fractal dimension of
$7/3$, or that of a convoluted surface (see also Sreenivasan & Meneveau Reference Sreenivasan and Meneveau1986; Sreenivasan, Ramshankar & Meneveau Reference Sreenivasan, Ramshankar and Meneveau1989; Sreenivasan Reference Sreenivasan1991; Vassilicos & Hunt Reference Vassilicos and Hunt1991). It appears that the act of averaging out the previously discussed large-scale fluctuations yields monofractal dynamics that provides an average sense of any underlying multifractal dynamics. This applies to both the volume-averaging procedure discussed above and the time-averaging procedure introduced in § 5.1. The latter implicitly assumes that the dynamics in these small, localized and quasi-homogeneous regions are also quasi-stationary. The additional observation of quasi-self-similarity in these break-up dynamics in an intermediate range of bubble sizes, as discussed in § 5.4.2, brings to mind the five-dimensional turbulent cascade in the three spatial dimensions, time and scale discussed by Cardesa, Vela-Martín & Jimenéz (Reference Cardesa, Vela-Martín and Jimenéz2017).
Appendix B. Estimating the energy dissipation rate
$\bar {\varepsilon }$
Recall from Part 1 and § 2 that the wavelength and deep-water-wave phase velocity of a characteristic wave are respectively used to estimate $L$ and
$u_L$ in the context of breaking waves. These characteristic scales may then be used to estimate the characteristic energy dissipation rate
$\bar {\varepsilon } \sim u_L^{3}/L$, and thereafter the global Kolmogorov and Hinze scales in (2.1) and (2.2). However, some experiments (Rapp & Melville Reference Rapp and Melville1990; Na et al. Reference Na, Chang, Huang and Lim2016) have suggested that the largest eddies and gaseous cavities have characteristic sizes and velocities that are an order of magnitude smaller. Other experiments have further suggested that most of the turbulence that originates from wave breaking initially resides in a region whose thickness is of the order of one wave height from the surface (Rapp & Melville Reference Rapp and Melville1990; Terray et al. Reference Terray, Donelan, Agrawal, Drennan, Kahma, Williams, Hwang and Kitaigorodskii1996; Thomson et al. Reference Thomson, Schwendeman, Zippel, Moghimi, Gemmrich and Rogers2016). These have motivated the choice of other quantities for the estimation of characteristic wave scales. For example, typical surface trajectories during wave breaking may have motivated the development of a ballistic scaling by Drazen et al. (Reference Drazen, Melville and Lenain2008) that uses the wave height to estimate
$L$ and the corresponding ballistic velocity to estimate
$u_L$. The corresponding ballistic time can be smaller than the wave period by an order of magnitude depending on the wave slope. The ballistic scaling was demonstrated to be adequate for single breaking events with moderate wave slope
$S$, where
$S=\sum _{i=1}^{N_m} a_i k_i$ may be interpreted as the maximum slope of the initial waveform. The waveform is assumed to be a wave packet with
$N_m$ component modes, and
$a_i$ and
$k_i$ are the initial amplitude and wavenumber, respectively, of the
$i$th mode. However, it has been observed (Loewen & Melville Reference Loewen and Melville1991; Melville Reference Melville1994; Loewen, O'Dor & Skafel Reference Loewen, O'Dor and Skafel1996; Drazen et al. Reference Drazen, Melville and Lenain2008; Deane et al. Reference Deane, Stokes and Callaghan2016a,Reference Deane, Stokes and Callaghanb) that the dissipation rate eventually saturates with increasing wave slope due to the presence of multiple breaking events. This necessitates a different choice of the characteristic scales for steeper waves or more general breaking patterns. In this work, for example, the characteristic maximum slope of the initial waveform is
$S = \sum _{i=1}^{3} a_i k_i = 0.76$, where
$i=2$ and
$i=3$ refer, respectively, to the first and second harmonics. This value of
$S$ is in the saturated regime identified by Drazen et al. (Reference Drazen, Melville and Lenain2008). The forthcoming discussion suggests that in the absence of additional information, the wavelength and phase velocity remain the most generic choices for the characteristic length and velocity scales, respectively, in breaking waves.
Two canonical scenarios are considered for energy conversion during wave breaking. First, consider the average energy transfer rate from coherent to turbulent motion in a single travelling wave of height $h$ and wavelength
$\lambda$. Before breaking occurs, the wave energy per unit width and wavelength
$\frac {1}{4}\rho _l g h^{2}$ is carried forward at the group velocity
$\sqrt {(g\lambda )/(8{\rm \pi} )}$. Suppose that the action of breaking is to transfer this energy into a volume with cross-sectional area
$h$ by
$h$ (Lamarre & Melville Reference Lamarre and Melville1994; Loewen & Melville Reference Loewen and Melville1994; Deane & Stokes Reference Deane and Stokes2002; Drazen et al. Reference Drazen, Melville and Lenain2008; Rojas & Loewen Reference Rojas and Loewen2010). The resulting average energy transfer rate per unit mass is
$\bar {\varepsilon } \sim 0.8c_p^{3}/\lambda$, where
$c_p$ is the wave phase velocity. Second, consider the average energy transfer rate from large to small waves in a travelling wavepacket with central wavelength
$\lambda$ and a corresponding central phase velocity
$c_p$. Gemmrich & Farmer (Reference Gemmrich and Farmer2004) and Babanin et al. (Reference Babanin, Chalikov, Young and Savelyev2010) reported wave and velocity spectra from individual wave-breaking events indicating the momentary upshifting of the peak wavenumber towards smaller scales before breaking, and then downshifting towards larger scales after breaking. This hints at the accumulation and then release of energy at some small limiting scale. A similar focusing process was observed in physical space during the prebreaking phase in the two-dimensional numerical study of a spilling breaker by Iafrati (Reference Iafrati2011), and in the experimental study of a plunging breaker by Bonmarin (Reference Bonmarin1989). Inspired by these observations, as well as the arguments by Kitaigorodskii (Reference Kitaigorodskii1983), one may estimate
$\bar {\varepsilon }$ as a ratio of the characteristic energy transported by the wave
$c_p^{2}$ to the time scale associated with the limiting breaking scale
$t_b \sim L_b / u_b$. Assuming
$L_b$ is a function of only
$\bar {\varepsilon }$ and
$g$, and taking
$u_b$ to be the phase velocity of a wave of length
$L_b$, one eventually obtains the estimate
$\bar {\varepsilon } \sim 4c_p^{3}/\lambda$. Both estimates for
$\bar {\varepsilon }$ are of the order
$O(c_p^{3}/\lambda )$. In the absence of additional constraints to restrict this estimate further, it is proposed that
$\bar {\varepsilon } \sim u_L^{3} / L$ be estimated by assuming
$L \sim \lambda$ and
$u_L \sim c_p$.
Appendix C. Bubble identification
The error incurred by the bubble identification algorithm introduced in § 4.1 is briefly discussed for a number of simple test cases. The reader is referred to the discussion of Chan et al. (Reference Chan, Dodd, Johnson, Urzay and Moin2018a) and Chan (Reference Chan2020) for a more detailed comparison of the employed grouping criterion with other criteria. First, the error in computing the volume of a single drop was determined by ballistically advecting the drop in a quiescent gas and by advecting it with a constant imposed velocity. In both cases, a uniform mesh was adopted, and there were $O(10)$ grid cells spanning the drop diameter. The volume error normalized by the grid-cell volume
$(\Delta D)^{3}$ was of the order of
$O(10^{-2}$–
$10^{-1})$, while the centroid error normalized by the grid-cell spacing
$\Delta D$ was of the order of
$O(10^{-4})$. Next, these errors were used to estimate the ability of the algorithm to distinguish between a closely spaced large drop–small drop pair and fluctuations in the volume of the large drop from one time step to another. The performance of the algorithm on this task is crucial in determining the accuracy of the tracking algorithm in § 4.2. It is first assumed that the normalized volume error of a drop,
$\Delta \mathcal {V}_{err}/(\Delta D)^{3}$, is proportional to its normalized surface area,
${\rm \pi} D^{2}/(\Delta D)^{2}$:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn22.png?pub-status=live)
The proportionality constant, $M$, is assumed to depend only on the grouping criterion. Next, suppose
$d$ and
$D$ are the diameters of the small and large drops, respectively. Then, in the limit that the two scenarios cannot be distinguished, one may write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn23.png?pub-status=live)
This yields the limiting size ratio
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn24.png?pub-status=live)
For the grouping criterion in this work, $M = O(10^{-4})$ and
$r/\sqrt {N} \simeq 20$–
$40$ for a given limiting bubble resolution
$N = d/(\Delta D)$, so the algorithm can distinguish between these two scenarios for size ratios of the order of
$O(10^{2})$ even when the small bubble is not well resolved. Note that this estimate for
$M$ leads to the non-dimensional bubble diameter error for the breaking-wave baseline ensemble
$\Delta D_{err}/L = O(10^{-6})$.
It is worth noting that many identification schemes are afflicted by the corner case of distinguishing two bubbles spaced a grid cell apart from a dumbbell-shaped bubble where two gaseous masses are connected by a thin gaseous bridge with the dimensions of a grid cell. The occurrence of this corner case does not necessarily suggest a deficiency in any of these schemes. Instead, it is representative of an inherent limitation of a sharp-interface field with finite resolution: when a thin gaseous bridge is numerically indistinguishable from a small under-resolved bubble or a small gaseous protrusion on the surface of a larger bubble, none of the geometries necessarily represent reality more accurately in the absence of additional information. Consequently, any decision by any scheme to favour any of the geometries is essentially arbitrary to a certain degree. The grouping criterion introduced in this work collectively identifies the gaseous masses and the gaseous bridge as a single bubble if the bridge is connected to cells in both gaseous masses with sufficiently large $1-\phi$. Otherwise, it returns two large bubbles and a small under-resolved bubble.
Appendix D. Event detection
The constraints required to track a bubble from one simulation snapshot to another in the event detection algorithm introduced in § 4.2 are discussed here. As a bubble is advected, its mass – and volume in an incompressible setting – remains constant between the snapshots even if it deforms, within the error identified in appendix C. The principle of mass conservation is necessarily satisfied even if one bubble breaks up into two or if two bubbles coalesce into one. Suppose bubble 0 breaks up into bubbles 1 and 2, or 1 and 2 coalesce to form 0. In the case of break-up, one may write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn25.png?pub-status=live)
while in the case of coalescence, one may write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn26.png?pub-status=live)
where $\mathcal {V}_i^{j}$ is the volume of bubble
$i$ in a flow snapshot
$j$,
$n$ denotes a generic flow snapshot and
$\Delta \mathcal {V}_{err}$ is the volume error discussed in appendix C. Also, if the simulation satisfies the CFL condition, then each fluid–fluid interface cannot traverse more than a single cell width every time step, multiplied by the maximum permissible Courant number
$C'$ for the employed advection scheme. Thus, the centroid of a bubble remains stationary between snapshots insofar as the permissible distance error is the product of the local grid spacing
$\Delta x$, the number of time steps between the snapshots
$N_t$ and
$C'$. This condition is also necessarily satisfied during break-up and coalescence. Denote the centroid of a bubble
$i$ in snapshot
$j$ as
$\boldsymbol {x}_i^{j}$. In the case of break-up, one may write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn27.png?pub-status=live)
while in the case of coalescence, one may write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_eqn28.png?pub-status=live)
Note that the identification algorithm also generates spatial errors in the centroid. However, as noted in appendix C, these errors are typically much smaller than the grid spacing, and will be neglected. The constraints (D1)–(D4) are sufficient for the identification of break-up and coalescence events, as well as continuing bubbles, using the bubble volumes and centroids. These constraints assume that all break-up and coalescence events are binary. Note also that (D1) and (D2) imply a critical size ratio $r$ above which events involving a small bubble and a large bubble cannot be distinguished from fluctuations in the volume of the large bubble between snapshots. This ratio is discussed in appendix C. Choosing an identification scheme with a lower volume error increases the critical size ratio keeping the resolution of the smallest bubble constant, allowing more relevant events to be captured accurately.
Appendix E. The role of coalescence
Given that the importance of coalescence events is expected to increase with the passage of time and the increase in the number density of small bubbles, coalescence statistics were investigated to determine the role of coalescence in the late wave-breaking stages. Note that coalescence events driven by large-scale dynamics such as buoyancy and turbulence are governed by the corresponding large-scale characteristic time scales (see e.g. Rodríguez-Rodríguez et al. Reference Rodríguez-Rodríguez, Gordillo and Martínez-Bazán2006). This means that while the precise moment of coalescence may not be captured accurately since most full-scale numerical simulations cannot feasibly resolve the film drainage that precedes coalescence, macroscopic coalescence statistics remain reasonably reliable unless the system is devoid of these large-scale external influences. There is, however, a possibility that numerical coalescence elevates the total observed coalescence flux. The coalescence flux $\overline{W_c}(D;t)$ across a cutoff size
$D$ was computed in an analogous fashion to
$\overline{W_b}(D;t)$ using the event detection algorithm in § 4.2. Figure 16 plots the time evolution of
$\overline{W_c}^{D}/\overline{W_b}^{D}(t)$ where size averaging is performed over the same bubble-size subrange used in figure 15. The plot indicates that
$\overline{W_c}^{D}$ is consistently smaller than
$\overline{W_b}^{D}$, indicating that coalescence is possibly not a major player in the dynamics even towards
$t=4$. A closer look at the statistics suggests that the complete story may be more nuanced. Figure 17 plots the time-averaged coalescence probability
$\overline{\check {q}_c}^{T}(D_{cp}^{3}|D_{cc}^{3})$, which is the probability that a coalescence event generating a child (descendant) bubble of size
$D_{cc}$ involves a parent (ancestral) bubble of size
$D_{cp}$ and another parent bubble such that the total gaseous volume remains constant, for the two time intervals identified in §§ 3.3, 5.1 and 5.2. The distributions reveal that for most child bubbles, large-size-ratio events involving one small parent bubble and one large parent bubble dominate the dynamics. The dominance of large-size-ratio coalescence events may be analytically supported from kinetic theory arguments, such as the one posited by Chan et al. (Reference Chan, Dodd, Johnson, Urzay and Moin2018b), which argues that large-size-ratio collisions are favoured when the size distribution decays with increasing bubble size. In short, the dynamics of coalescence appears to be dominated by parent bubbles that are beyond the reach of the current simulations. It is possible, then, that
$\overline{W_c}$ has been underestimated in these simulations, and the jury is still out as to the importance of coalescence in the late wave-breaking stages. A study of coalescence is deferred to later work where the capability to computationally capture these small bubbles is included. One possible approach is to replace bubbles under-resolved by the current Eulerian treatment with Lagrangian point particles, as alluded to in Part 1 and to be discussed in forthcoming work.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig16.png?pub-status=live)
Figure 16. The ratio of the size-averaged coalescence flux to the size-averaged break-up flux, $\overline{W_c}^{D}/\overline{W_b}^{D}(t)$, against non-dimensional time
$t$ for the 20-realization baseline ensemble subset. Coalescence events involving subgrid parent (ancestral) bubbles are excluded. For a description of the shaded regions, excluded break-up events, error bars and the bubble-size subrange over which the fluxes were averaged, refer to the caption of figure 15.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig17.png?pub-status=live)
Figure 17. The time-averaged coalescence probability $\overline{\check {q}_c}^{T}(D_{cp,j}^{3}|D_{cc}^{3})$ against normalized parent (ancestral) bubble volume
$D_{cp}^{3}/D_{cc}^{3}$ for the 20-realization baseline ensemble subset during the (a,c,e) early (
$t=2.30$–
$3.14$) and (b,d,f) late (
$t=3.45$–
$3.95$) wave-breaking stages. The non-dimensional child (descendant) bubble sizes considered are (a,b) between
$1.86 \times 10^{-2}$ and
$2.23 \times 10^{-2}$, (c,d) between
$2.68 \times 10^{-2}$ and
$3.23 \times 10^{-2}$ and (e,f) between
$3.54 \times 10^{-2}$ and
$4.25 \times 10^{-2}$. Only parent bubbles of radii larger than the mesh resolution are considered. The vertical dashed lines demarcate the parent bubble volumes in the coalescence event where one of the parent bubble radii corresponds to this resolution limit. The error bars denote one standard error over the baseline ensemble.
Appendix F. Spatial distribution of the liquid volume fraction
In order to provide a general sense of how the bubbles are spatially distributed during the wave-breaking process, figure 18 plots the ensemble-averaged and spanwise-averaged liquid volume fraction, $\phi$, at two representative time instances, one just before the first characteristic time interval and one just after the second characteristic time interval identified in table 2. Note that the void fraction may be correspondingly computed as
$(1-\phi )$. In general agreement with the other works discussed in §§ 3.3 and 5.1, high-void-fraction mixed-phase regions are initially concentrated near the areas with rapid topology change as evidenced in figure 18(a), and then gradually evolve into a number of diffuse bubble clouds near the wave surface as evidenced in figure 18(b).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20210215161514726-0248:S0022112020010848:S0022112020010848_fig18.png?pub-status=live)
Figure 18. The ensemble-averaged and spanwise-averaged liquid volume fraction $\phi$ for the 30-realization baseline ensemble, computed at (a)
$t=2.23$ and (b)
$t=4.04$. The solid line denotes the
$\phi =0.5$ isocontour.