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

(2.3)\begin{equation} \frac{\partial [\bar{f}(D;t) D^{3}]}{\partial t} + \frac{\partial [\overline{v_D f}(D;t) D^{3}]}{\partial D} = \bar{H}(D;t). \end{equation}

Correspondingly, volume averaging the kernel-based equation (3.7) in Part 1 yields

(2.4)\begin{equation} \frac{\partial [\bar{f}(D;t) D^{3}]}{\partial t} = \overline{T_b}(D;t) + \overline{T_c}(D;t). \end{equation}

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}$:

(2.5)\begin{equation} \frac{\partial [\bar{f}(D;t) D^{3}]}{\partial t} = -\frac{\partial[\overline{v_D f}(D;t) D^{3}]}{\partial D} = \overline{T_b}(D;t). \end{equation}

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:

(2.6)\begin{equation} \overline{T_b}(D;t) = \int_D^{\infty} \mathrm{d} D_p \check{q}_b\left(D;t|D_p\right) \overline{g_b f}\left(D_p;t\right) D^{3} - \overline{g_b f}(D;t) D^{3}. \end{equation}

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

(2.7)\begin{equation} \overline{W_b}(D;t) = \int_0^{D} \mathrm{d} D_c D_c^{3} \int_D^{\infty} \mathrm{d} D_p \check{q}_b \left(D_c;t|D_p\right) \overline{g_b f}(D_p;t). \end{equation}

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

(2.8)\begin{gather} \overline{I_p}(D_p;t|D) = \int_0^{D} \mathrm{d} D_c D_c^{3} \check{q}_b \left(D_c;t|D_p\right) \overline{g_b f}\left(D_p;t\right), \end{gather}

(2.9)\begin{gather}\overline{I_c}(D_c;t|D) = \int_D^{\infty} \mathrm{d} D_p D_c^{3} \check{q}_b \left(D_c;t|D_p\right) \overline{g_b f}\left(D_p;t\right). \end{gather}

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.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

(3.1)\begin{gather} \rho = \phi \rho_l + (1-\phi) \rho_g, \end{gather}

(3.2)\begin{gather}\mu = \phi \mu_l + (1-\phi) \mu_g, \end{gather}

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

(3.3)\begin{equation} \eta(x,t=0) = \frac{1}{2{\rm \pi}} \left[ S_1 \cos(2{\rm \pi} x) + \frac{1}{2} S_1^{2} \cos(4{\rm \pi} x) + \frac{3}{8} S_1^{3} \cos(6{\rm \pi} x) \right], \end{equation}

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.

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):

(3.4)\begin{gather} u(x,y,t=0) = S_1\sqrt{1 + S_1^{2}}\cos(2{\rm \pi} x)\exp(2{\rm \pi} y), \end{gather}

(3.5)\begin{gather}v(x,y,t=0) = S_1\sqrt{1 + S_1^{2}}\sin(2{\rm \pi} x)\exp(2{\rm \pi} y), \end{gather}

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.

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.

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:

(3.6)\begin{align} E_t &= \frac{\hat{E}_k + \hat{E}_p + \hat{E}_s}{\hat{E}_n} = \frac{\text{kinetic energy}+\text{potential energy}+\text{surface energy}}{\text{reference energy for normalization}} \nonumber\\ &= \frac{\left[\displaystyle\int_{\varOmega_d} \dfrac{1}{2} \rho \left|{\hat{\boldsymbol{u}}}\right|^{2} \mathrm{d} {\hat{\boldsymbol{x}}}\right] + \left[\displaystyle\int_{\varOmega_d} \rho g \hat{y} \, \mathrm{d} {\hat{\boldsymbol{x}}} - \dfrac{1}{8} \left(\rho_g-\rho_l\right) g L^{4}\right] + \left[\displaystyle\int_{\varOmega_d} \sigma \left| \hat{\boldsymbol{\nabla}}\phi\right| \mathrm{d} {\hat{\boldsymbol{x}}} - \sigma L^{2}\right]}{\left( \dfrac{1}{2}\rho_l u_L^{2}\right)\left(\dfrac{1}{2}L^{3}\right)}, \end{align}

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}$.

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.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.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

(5.1)\begin{equation} \left\langle N_b(t) \right\rangle = \sum_{j=1}^{N_{bin}} \bar{f}(D_j/L;t) \varDelta(D_j/L), \end{equation}

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.

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.

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).

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.

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.

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.

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.

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

(5.2)\begin{equation} 2 = \sum_{j=1}^{N_{bin}} \check{q}_b\left(D_{c,j}^{3};t|D_p^{3}\right) \varDelta\left(D_{c,j}^{3}/D_p^{3}\right). \end{equation}

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.

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.

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.

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.

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.

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.