1. Introduction
Air-entraining free-surface flows driven by underlying turbulence are present and important in a variety of natural and engineering applications, including gas exchange between the atmosphere and the ocean, water quality in lakes and rivers, and the wake around ships. In many of these flows – e.g. dam spillways, rivers (Chanson Reference Chanson1996) and ship wakes (Hendrickson et al. Reference Hendrickson, Weymouth, Yu and Yue2019) – interaction of the free surface with strong turbulence in the liquid beneath leads to entrainment of air (Brocchini & Peregrine Reference Brocchini and Peregrine2001), as illustrated in figure 1. For this interaction to lead to air entrainment, the turbulence must be strong relative to the restoring effects of gravity and surface tension (see ‘region 2’ from the regime map by Brocchini & Peregrine Reference Brocchini and Peregrine2001, figure 4). Unlike breaking waves, where air entrainment is driven by large displacement of the free surface and bubble populations are often understood in terms of entrainment events tied to individual breakers (Deike Reference Deike2022), the multi-scale nature of underlying turbulence interacting with the free surface means that air-entraining free-surface turbulence (FST) entrains bubbles continuously across a range of scales (Yu, Hendrickson & Yue Reference Yu, Hendrickson and Yue2020). In this work, we focus on these air-entraining FST flows.
Rendering from DNS of an air-entraining FST flow (see § 4 for details) with a sketch of the five physical mechanisms that effect the bubble population, in terms of the bulk bubble size distribution 
$N(a)$ (see (2.1) for definitions), and the total entrained volume 
$V$ (see (2.8) for definitions). We will show that this flow is dominated by entrainment and degassing and the other mechanisms (in grey) are negligible.
The size distribution of bubbles in air-entraining flows is often of critical importance and is the result of five physical mechanisms: entrainment where bubbles are created at the free surface; fragmentation where turbulence breaks a bubble into multiple bubbles; coalescence where multiple bubbles join into one; degassing where a bubble bursts at the free surface; and dissolution where air dissolves into the surrounding water. If we assume conditions such that dissolution and coalescence are negligible, then we are left with entrainment, fragmentation and degassing. All could (simultaneously) interact to determine the bubble population. We define an air-generating flow as one where the total volume of bubbles is increasing. In an air-generating flow, the entrainment mechanism is clearly important. In addition to growth of the bubble population, this leaves two possible mechanics to balance entrainment: fragmentation or degassing. Fragmentation as the other dominant mechanism during air generation has been studied widely, especially for breaking waves. However, there is little understanding of bubble populations when degassing (rather than fragmentation) is dominant during air generation, the configuration illustrated in figure 1.
To quantify bubble populations, we consider the bulk bubble size distribution 
$N(a)$, where 
$N(a)\,\delta a$ is the number of bubbles in the flow of (effective) radius 
$[a,a+\delta a]$. For bubble fragmentation in turbulence, Garrett, Li & Farmer (Reference Garrett, Li and Farmer2000) show that for bubbles of radius larger than the Hinze scale, 
$a_H$, a fragmentation cascade arises, producing a power law 
$N(a)\propto a^{\beta }$, where 
$\beta =-10/3$. Below the Hinze scale, surface tension prevents fragmentation, and different power laws are observed (Deane & Stokes Reference Deane and Stokes2002). Studying breaking waves, Deane & Stokes (Reference Deane and Stokes2002) observe 
$\beta =-10/3$ for super-Hinze bubbles during the air-generating period (referred to by them as the acoustically active period), and conclude that fragmentation and entrainment are the primary mechanisms during this period. Many have since reported 
$\beta \approx -10/3$ for super-Hinze bubbles during the air-generating period of breaking waves (e.g. Deike Reference Deike2022).
After the air-generating period, when the total volume of air begins to decrease, Deane & Stokes (Reference Deane and Stokes2002) observe steepening of the bubble size distribution (
$\beta <-10/3$), which they attribute to degassing and dissolution becoming relevant in the absence of entrainment and fragmentation. They refer to this period as the quiescent period. Previous work on the effects of degassing typically focus on this quiescent period to explain how degassing evolves a bubble population after entrainment and fragmentation have already established the size distribution (e.g. Callaghan, Deane & Stokes Reference Callaghan, Deane and Stokes2013; Deike, Melville & Popinet Reference Deike, Melville and Popinet2016). In this work, we study the effect of degassing during the air-generating period, when degassing acts to balance entrainment. The wide agreement on 
$\beta \approx -10/3$ suggests that this effect may be negligible compared to fragmentation for breaking waves; however, in this work, we will show that for air-generating FST, degassing is in fact the dominant balance to entrainment during the air-generating period.
To investigate the effect of degassing-dominated free-surface flows on bubble populations, particularly the power-law slope 
$\beta$, we consider air-generating flows driven by strong underlying FST. In contrast to breaking waves, where the energy to create bubbles comes from a (downward) mean flow, in FST this energy comes directly from the underlying turbulence. In the absence of a mean flow, we theorize that degassing will be stronger relative to other entrainment-balancing mechanisms, such as fragmentation, and our results confirm that degassing is dominant for FST. In § 2, we review the governing equations for bubble populations in air-generating flow. In § 3, we use these equations along with a simple model of degassing to define degassing-dominated bubble populations and predict the resulting 
$\beta$. In § 4, we perform direct numerical simulations (DNS) of a canonical FST flow (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999; Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019), and use a robust volume-conserving bubble tracking method – Eulerian label advection (ELA) (Gaylo, Hendrickson & Yue Reference Gaylo, Hendrickson and Yue2022) – to directly obtain entrainment and degassing statistics. Studying large Weber numbers 
${{We}}$, where fragmentation is strongest, we find that degassing is dominant over fragmentation, independent of Froude number 
${{Fr}}$. Measurements of 
$N(a)$ from DNS (independent of ELA) agree with our predicted 
$\beta$ for degassing-dominated bubble populations. In § 5, we discuss how these results can be extended to describe degassing-dominated flows at large-scale Reynolds numbers, as well as the effects of including surface tension. Surface tension decreases the strength of fragmentation. Because fragmentation is already negligible compared to degassing, the inclusion of surface tension is not expected to affect our findings regarding degassing dominance for the air-entraining FST flows that we study.
2. Modelling bubble populations in air-generating FST
The evolution of the bulk bubble size distribution 
$N(a)$ is described by a population balance equation (PBE), with source terms describing each physical mechanism (Sporleder et al. Reference Sporleder, Borka, Solsvik and Jakobsen2012). We assume short time periods and large bubbles such that the effect of dissolution is negligible, and small void fractions such that the effect of coalescence is negligible. This gives the PBE
\begin{equation} \partial N(a) / \partial t = S_f(a) + I(a) - D(a) {,} \end{equation}
where 
$S_f(a)$ describe fragmentation, 
$I(a)$ entrainment, and 
$D(a)$ degassing. Below, we discuss models for these three terms in air-entraining flows.
Fragmentation can be split into 
$S_f=S_f^+-S_f^-$ for the creation and destruction of bubbles by fragmentation. The destruction term is
\begin{equation} S^-_f(a) = \varOmega(a)\,N(a){,} \end{equation}
where 
$\varOmega (a)$ is the fragmentation rate (units 
$1/T$). For large 
${{We}}$ where surface tension is negligible,
\begin{equation} \varOmega(a) = C_\varOmega \, \varepsilon^{1/3}\, a^{-2/3} {,} \end{equation}
where 
$\varepsilon$ is the turbulent dissipation rate, and 
$C_\varOmega \approx 0.42$ (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). For moderate 
${{We}}$ where surface tension is relevant, it can be incorporated by making 
$C_\varOmega$ depend on 
${{We}}$ (Martínez-Bazán et al. Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montaes and Lasheras2010). While bubble fragmentation in turbulence, particularly at moderate 
${{We}}$, is still an area of active research (e.g. Qi et al. Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022), many models agree with (2.3) in the limit of large 
${{We}}$ (Martínez-Bazán et al. Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montaes and Lasheras2010), and, relevant to the analysis here, decreasing 
${{We}}$ (i.e. increasing the strength of surface tension relative to the flow) will decrease 
$C_\varOmega$, meaning that surface tension effects make 
$\varOmega (a)$ smaller.
For our analysis, the fragmentation destruction term is most relevant as it has the same sign as degassing in (2.1). Still, it is useful to consider the relative magnitude of the creation term
\begin{equation} S^+_f(a) = \int_{a}^{\infty} m(a')\,f(a; a')\,\varOmega(a')\,N(a') \,\mathrm{d} a' {,} \end{equation}
where 
$m(a')$ is the average number of daughter bubbles created by fragmentation of a parent bubble of radius 
$a'$, and 
$f(a; a')$ is the daughter-size distribution, expressed as a probability distribution function of daughter radius 
$a$ for a given parent radius 
$a'$. The choice of 
$f(a; a')$ has little impact on the equilibrium shape of 
$N(a)$ (Qi, Masuk & Ni Reference Qi, Masuk and Ni2020; Gaylo, Hendrickson & Yue Reference Gaylo, Hendrickson and Yue2021), so for estimation, we will consider a simple daughter-size distribution where 
$m$ is a tunable constant and 
$f(a; a')$ is a Dirac delta function, 
$\delta (a' - m^{-1/3} a)$. With this distribution, (2.4) becomes
\begin{equation} S^+_f(a) \sim m^{4/3}\,\varOmega(m^{1/3} a)\,N(m^{1/3} a) {.} \end{equation}
Using (2.3) and assuming a power law 
$N(a)\propto a^{\beta }$,
\begin{equation} S^+_f(a) \sim S^-_f(a)\, m^{10/9+\beta/3} {.} \end{equation}
The DNS measurement by Gaylo, Hendrickson & Yue (Reference Gaylo, Hendrickson and Yue2023) suggests that 
$m\gtrsim 2$ describes turbulent fragmentation. Relevant to the analysis in § 3, this shows that 
$S^+_f(a) \lesssim S^-_f(a)$ for any 
$\beta \lesssim -10/3$.
For FST entrainment, Yu et al. (Reference Yu, Hendrickson and Yue2020) study an energy balance between turbulence, surface tension and gravity. For large bubbles where gravity dominates surface tension (corresponding to 
$a$ larger than the capillary scale 
$a_c = 1/2 \sqrt {\sigma /\rho g}$), this balance predicts 
$I(a)\propto a^{-10/3}$. This prediction for entrainment is unrelated (apart from dependence on the Kolmogorov spectrum) to 
$\beta ={-10/3}$ for fragmentation-dominated populations; however, this coincidence of power-law slopes makes it difficult to confirm 
$I(a)\propto a^{-10/3}$ through measurements of 
$N(a)$. In § 4, we use ELA to directly measure the entrainment distribution and confirm 
$I(a)\propto a^{-10/3}$.
Degassing depends on the vertical movement of a bubble to the free surface. We expect the movement of a single bubble to be driven by buoyancy and the surrounding turbulence, and (excluding very large void fractions) to be independent of the presence of other bubbles. If the degassing statistics of individual bubbles are independent, then the degassing statistics of the bulk population of bubbles, 
$D(a)$, will simply be the individual statistics weighted by the number of bubbles 
$N(a)$. Thus independence means that the degassing distribution 
$D(a)$ is directly proportional to 
$N(a)$. The constant of proportionality defines the degassing rate 
$\varLambda (a)$ (units 
$1/T$) such that
\begin{equation} D(a) \equiv \varLambda(a)\,N(a) {.} \end{equation}
In § 3.2, we elucidate the power-law scaling of 
$\varLambda (a)$ in FST, and in § 4, we confirm this scaling through direct measurement with ELA.
In addition to 
$N(a)$, it is useful to consider the total entrained volume 
$V=(4 {\rm \pi}/3) \int N(a)\, a^3 \,\mathrm {d} a$. Applying the same integral to (2.1),
\begin{equation} \dot{V} = Q_I - Q_D {,} \end{equation}
where 
$Q_I=(4 {\rm \pi}/3) \int I(a)\,a^3 \,\mathrm {d} a$ is the flux of air volume from above to beneath the free surface (entrainment flux), and 
$Q_D=(4 {\rm \pi}/3) \int D(a)\,a^3 \,\mathrm {d} a$ is the flux from beneath to above the free surface (degassing flux). Fragmentation does not contribute to 
$\dot {V}$. We define 
$\dot {V}>0$ as air-generating flow. To quantify the importance of degassing, we consider the ratio between degassing flux and entrainment flux,
\begin{equation} \mathcal{D} \equiv Q_D/Q_I {,} \end{equation}
such that 
$\mathcal {D}\in [0,1)$ for air-generating flow. For flows with negligible entrainment, such as during the quiescent period of breaking waves, 
$\mathcal {D}\gg 1$. For breaking waves during the air-generating period, Deane & Stokes (Reference Deane and Stokes2002) show that entrainment and fragmentation are the primary mechanisms, implying that degassing is negligible (
$\mathcal {D}\ll 1$). Our interest is in air-generating flows where degassing is important (
$\mathcal {D}= \operatorname {O}(1)$).
3. Degassing-dominated bubble populations in air-generating flows
We now apply the models described in § 2 to obtain the shape of the bulk bubble size distribution 
$N(a)$ for degassing-dominated populations. Equation (2.1) includes fragmentation, degassing and entrainment. Noting that both (2.2) and (2.7) are linear with 
$N(a)$, we separate positive and negative terms,
\begin{equation} \partial N(a) / \partial t = \left[I(a)+S_f^+(a)\right] - \left[\varLambda(a)+\varOmega(a) \right] N(a) {.} \end{equation}
If 
$\varLambda (a)>\varOmega (a)$, then a bubble is more likely to degas than fragment. For 
$\varLambda (a)\gg \varOmega (a)$, we expect fragmentation to have a negligible effect on 
$N(a)$ compared to degassing, defining a degassing-dominated population.
3.1. Degassing-dominated PBE
Recalling (2.6), 
$S^+_f(a)$ is similar in magnitude to 
$S^-_f(a)$, so 
$\varLambda (a)\gg \varOmega (a)$ allows us to remove both the fragmentation terms from (3.1). This gives the degassing-dominated PBE, which, after rearranging terms, reads
\begin{equation} I(a)= \varLambda(a)\,N(a) + \partial N(a) / \partial t {.} \end{equation}
We see that entrainment is balanced by two terms: degassing, and increasing bubble size distribution 
$\partial N(a) / \partial t$. We define the ratio of the first term to the second term:
\begin{equation} \gamma(a) \equiv [\varLambda(a)\,N(a)]\big/[I(a)-\varLambda(a)\,N(a)] {.} \end{equation}
As a measure of this ratio across radii, we integrate (weighted by bubble volume 
$a^3$) the numerator and denominator to define 
$\varGamma \equiv Q_D/(Q_I-Q_D)$. This integral measure can be related directly to 
$\mathcal {D}$ as 
$\varGamma =\mathcal {D}/(1-\mathcal {D})$. For a given radius, 
$\gamma (a) \ll 1$ means that 
$\varLambda (a)\,N(a) \ll \partial N(a) / \partial t$, and (3.2) becomes 
$\partial N(a) / \partial t \approx I(a)$, giving 
$N(a)\sim t\,I(a)$. Because the degassing term grows linearly with 
$N(a)$, 
$\varLambda (a)\,N(a) \ll \partial N(a) / \partial t$ is true only for time scales less than order 
$1/\varLambda (a)$. Apart from these short time scales, we expect 
$\varLambda (a)\,N(a) > \partial N(a) / \partial t$, corresponding to 
$\gamma (a) > 1$, for degassing-dominated populations. In this case, the primary balance is between entrainment and degassing, and (3.2) becomes
\begin{equation} N(a)\approx I(a)/\varLambda(a) {.} \end{equation}
Given an entrainment distribution 
$I(a)\propto a^{-10/3}$ (Yu et al. Reference Yu, Hendrickson and Yue2020), and assuming that degassing follows a power law 
$\varLambda (a)\propto a^{\alpha }$, the power-law slope of 
$N(a)$ is
\begin{equation} \beta=-10/3-\alpha {.} \end{equation}3.2. Power-law scaling of degassing in free-surface turbulence
Our objective is to find 
$\beta$ for degassing-dominated populations using (3.5). To predict 
$\alpha$, we derive a characteristic bubble depth 
$L_\varLambda$ and bubble rise velocity 
$U_\varLambda$ such that 
$\varLambda (a)\propto U_\varLambda /L_\varLambda$. The scaling of 
$U_\varLambda$ with bubble radius is a widely studied problem not unique to air-generating flows. A challenge is the characteristic depth 
$L_\varLambda$, which should describe some initial depth of an entrained bubble. For FST, turbulent eddies drive entrainment, and to derive 
$L_\varLambda$, we start by estimating how deep a single eddy will push a bubble. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) show that air-generating FST creates a region of nearly homogeneous isotropic turbulence (HIT) near the free surface. For HIT, an eddy of size 
$\ell$ has velocity 
$u_\ell \sim \varepsilon ^{1/3} \ell ^{1/3}$ and is coherent over 
$t_\ell \sim \varepsilon ^{-1/3} \ell ^{2/3}$. A bubble of radius 
$a$ has an (added) mass 
$m\sim \rho _w a^3$ and (at large 
${{Re}}$) feels a force from the eddy 
$F\sim \rho _w a^2 u_\ell ^2$, giving an acceleration 
$\dot {w}\sim u_\ell ^2/a$. This moves the bubble to a depth 
$z\sim {t_\ell }^2 \dot {w}=\ell ^2/a$. To apply this single-eddy model to HIT where a range of eddy sizes are present, we assume that the energy-containing eddies 
$\ell \sim u_{rms}^3 \varepsilon ^{-1}$ are primarily responsible for 
$L_\varLambda$, giving 
$L_\varLambda \sim u_{rms}^6 \varepsilon ^{-2} a^{-1}$.
For characteristic velocity 
$U_\varLambda$, buoyant rise or turbulent advection can be relevant. The terminal buoyant rise velocity of a bubble in quiescent flow, 
$W_t$, is characterized by the bubble Reynolds number 
${{Re}}_b \equiv 2 a W_t/\nu _w$. The bubbles in our DNS are in the inertial regime 
$1<{{Re}}_b<100$, where
\begin{equation} W_{t} = 0.144\,{g}^{5/6}\,\nu_w^{2/3}\,{(2a)}^{3/2} {,} \end{equation}
with the prefactor coming from experiments (Wallis Reference Wallis1974; Park et al. Reference Park, Park, Lee and Lee2017). While this 
$W_{t}$ is appropriate for our bubbles, we note that there are different 
${{Re}}_b$ regimes, and that corrections for surface effects are available (Park et al. Reference Park, Park, Lee and Lee2017). For small bubbles (small 
$W_t$), we expect the movement of the bubbles to be dominated by turbulent advection rather than buoyant rise. The dominance of turbulent advection implies a characteristic velocity for small bubbles 
$U_\varLambda \sim u_{rms}$ rather than 
$U_\varLambda \sim W_t$ for large bubbles. Assuming that 
$u_{rms}=W_{t}$ defines the transition between these two regimes, we solve (3.6) for 
$a$ to define the critical radius
\begin{equation} a_{\varLambda} = 1.82\,u_{rms}^{2/3}\,\nu_w^{4/9}\,{g}^{-5/9} {.} \end{equation}
Although more advanced models consider the interaction between buoyancy and turbulence and incorporate surface tension effects (e.g. Salibindla et al. Reference Salibindla, Masuk, Tan and Ni2020; Ruth et al. Reference Ruth, Vernet, Perrard and Deike2021), we find that this simple two-regime model is sufficient to explain the power-law scaling of degassing at large 
${{We}}$.
Using these 
$L_\varLambda$ and 
$U_\varLambda$ with 
$\varLambda (a)\propto U_\varLambda /L_\varLambda$, we obtain a split power law for 
$\alpha$. For large buoyancy-driven bubbles, 
$\alpha =5/2$, and for small turbulence-driven bubbles, 
$\alpha =1$:
\begin{equation} \varLambda(a) \propto \varLambda_0 \begin{cases} (a/a_\varLambda), & a < a_\varLambda, \\ (a/a_\varLambda)^{5/2}, & a > a_\varLambda, \end{cases} \end{equation}
where 
$\varLambda _0 = 1.82\,u_{rms}^{-13/3}\,\varepsilon ^{2}\,\nu _w^{4/9}\,{g}^{-5/9}$. Thus, using (3.5) for degassing-dominated populations, we obtain
\begin{equation} \beta = \begin{cases} -13/3, & a < a_\varLambda, \\ -35/6, & a > a_\varLambda. \end{cases} \end{equation}4. Direct numerical simulation of air-generating FST
In this section, we perform DNS of a canonical FST flow. During the air-generating period, we use ELA (Gaylo et al. Reference Gaylo, Hendrickson and Yue2022) to directly obtain 
$\mathcal {D}$, 
$I(a)$ and 
$\varLambda (a)$. We find that, unlike breaking waves (e.g. Deike Reference Deike2022), the bubble population in this FST is degassing-dominated. Further, our DNS results confirm the theoretical predictions of § 3 for degassing-dominated populations.
4.1. Methodology
We study FST generated by an initial shear profile,
\begin{equation} u(z, t=0)/U = 1 - 0.9988 \operatorname{sech}(0.88137 \,z/L) {,} \end{equation}
characterized by the shear velocity 
$U$ and shear length 
$L$. At sufficient Reynolds number 
${{Re}}= U L/\nu _w$, a small, random initial perturbation causes the flow to become turbulent (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999). Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019, Reference Yu, Hendrickson and Yue2020) show that for Froude number squared (
${{Fr}}^2= U^2/L {g}$) greater than a critical 
${{Fr}}^2_{cr}\approx 5$, this FST entrains bubbles. There have also been extensive DNS characterizing this FST at low 
${{Fr}}^2$ (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999; Shen, Triantafyllou & Yue Reference Shen, Triantafyllou and Yue2000). Following all previous work, we choose 
${{Re}}=1000$, and set the domain size to 
$(10.472 L)^2 \times 6 L$, where initially 
$z/L\in [-4,0]$ is water and 
$z/L\in [0,2]$ is air (density ratio 
$\rho _a/\rho _w=0.00123$, viscosity ratio 
$\mu _a/\mu _w=0.0159$). This gives an initial free-surface area 
$A/L^2=10.472^2$. For conciseness, we set 
$U$ and 
$L$ to unity. Thus all otherwise dimensional values that follow should be interpreted as normalized by 
$U$ and 
$L$ such that they are non-dimensional. The rendering in figure 1 comes from a simulation of 
${{Fr}}^2=15$ at time 
$t=60$, with the shear velocity going from left to right.
For DNS, we use the same scheme as Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019, Reference Yu, Hendrickson and Yue2020). It is a pressure-projecting second-order staggered-grid finite-volume scheme, where air and water are evolved using the conservative volume-of-fluid method (Weymouth & Yue Reference Weymouth and Yue2010). Our interest is 
${{We}}=U^2 L/(\sigma /\rho _w)\gg 1$ where surface tension is negligible, and Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) validate that not modelling surface tension accurately captures this limit. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) also study grid convergence to show that a grid 
$384^2 \times 256$ is sufficient for DNS of this flow, giving a grid size 
$\varDelta \approx 0.027$. We include an additional grid convergence study in the Appendix to confirm that this grid size sufficiently resolves entrainment and degassing. Following Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019), only bubbles of radius larger than 
$a_{res}=1.5 \varDelta$ are considered resolved.
For each 
${{Fr}}^2$ studied (see table 1), we repeat the simulation with different realizations of the random initial perturbation to obtain ensemble statistics. We have at least six simulations at each 
${{Fr}}^2$, and 
$\operatorname {O}(10^3)$ resolved entrainment and degassing events during the measured air-generating period for all but the smallest 
${{Fr}}^2$ (close to 
${{Fr}}^2_{cr}$). Because we normalize 
$U$ and 
$L$ to unity, 
${g}$ is varied to obtain different 
${{Fr}}^2$. A more intuitive interpretation of table 1 would be to consider a shear flow with fixed 
${g}$ and length scale 
$L$, where larger 
${{Fr}}^2$ mean larger 
$U^2$. Recall that reported units are 
$L$ for length and 
$L/U$ for time.
Summary of FST simulations performed and values obtained during 
$t\in [40,70]$, including 
$95\,\%$ confidence intervals. Here, 
${N}_{sim}$ is the number of ensemble simulations, 
${N}_{I}$ is the number of entrainment events, and 
${N}_{D}$ is the number of degassing events. Note that all otherwise dimensional quantities are normalized by the shear velocity 
$U$ and length 
$L$.
For comparing this shear flow FST to general FST, it is useful to use turbulence characteristics. We use a water domain spatial average
\begin{equation} {\langle \cdot \rangle_\delta}{{\equiv} \frac{\iiint \cdot\ f \,\mathrm{d}\kern0.06em \boldsymbol{x}}{\iiint f \,\mathrm{d}\kern0.06em \boldsymbol{x}} }\quad\text{for } z>-\delta {,} \end{equation}
where 
$f$ is the volume of fluid (
$f=1$ in water), and 
$\delta =0.2$ captures the near-surface HIT region (Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019). Figure 2 shows the turbulent dissipation rate 
$\varepsilon$ and turbulent kinetic energy 
$k$. As expected, there is little change in the strength of near-surface turbulence with 
${{Fr}}$, and a quasi-steady period exists when 
$k$ remains approximately constant (Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019). We use 
$u_{rms}\equiv \sqrt {2 {\langle k \rangle _\delta }/3}$ and 
${\langle \varepsilon \rangle _\delta }$ averaged over the quasi-steady period to calculate quantities from § 3.
Ensemble average (a) turbulent dissipation rate and (b) turbulent kinetic energy in the near-surface region, for solid lines: blue, 
${{Fr}}^2=5$; light green, 
${{Fr}}^2=8$; dark green, 
${{Fr}}^2=10$; brown, 
${{Fr}}^2=15$; red, 
${{Fr}}^2=20$. Dashed lines indicate 
$t\in [40,70]$ over which we average to obtain the values in table 1.
While the values that we present are normalized by the shear flow 
$U$ and 
$L$, we note that the measured characteristic scales of the turbulence, 
$U_T\sim u_{rms}$ and 
$L_T\sim u_{rms}^3/\varepsilon$, allow comparison to general FST. For example, we calculate a turbulent Froude number squared, 
${{Fr}}^2_{T} = \varepsilon / u_{rms} \,{g}$. Because 
${\langle k \rangle _\delta }$ and 
${\langle \varepsilon \rangle _\delta }$ do not change significantly with 
${{Fr}}$ (see table 1), normalizing by the turbulence characteristic scales rather than the shear characteristic scales would represent only a change in a constant rather than a change in any scaling with 
${g}$.
4.2. Measuring degassing and entrainment with ELA
The evolution of bubbles can be described using a volume tracking matrix (Gaylo et al. Reference Gaylo, Hendrickson and Yue2022), illustrated in figure 3. Consider bubbles at time 
$t^n$ labelled 
$j=1,\ldots, M^n$, where 
$M^{n}$ is the number of bubbles (including one for the air above the free surface). At time 
$t^{n+1}=t^n+T_s$, where 
$T_s$ is the snapshot interval, we have another set of bubbles labelled 
$i=1,\ldots, M^{n+1}$. The evolution of bubbles over 
$T_s$ is described by a volume tracking matrix, where each matrix element 
$q_{ij}$ gives the volume of air from a bubble 
$j$ at time 
$t^{n}$ that contributes to a bubble 
$i$ at time 
$t^{n+1}$. The volume tracking matrix provides a robust description of bubble evolution, and entrainment and degassing statistics can be extracted easily. We use ELA to track the movement of air volume on the grid level to (explicitly) calculate the volume tracking matrix. This method is proven to be volume-conservative, meaning that all air is tracked accurately (Gaylo et al. Reference Gaylo, Hendrickson and Yue2022; for a recent application, see Gaylo et al. Reference Gaylo, Hendrickson and Yue2023).
(a) Illustration of entrainment (red) and degassing (green). (b) The associated volume tracking matrix. At 
$t^{n}$, we have 
$j=1$ and three bubbles 
$j=2,3,4$. At 
$t^{n+1}$, we have 
$i=1$ and two bubbles 
$i=2,3$.
For ELA, we must choose the snapshot interval 
$T_s$. We want 
$T_s$ small because we miss the entrainment and degassing events of a bubble if they are separated by less than 
$T_s$. However, very small 
$T_s$ will generate spurious events because it is numerically unclear if bubbles (or a bubble and the free surface) are separate if their interfaces are within 
$\lesssim \varDelta$ (Chan et al. Reference Chan, Dodd, Johnson and Moin2021). To decrease spurious events, we use interface reconstruction to avoid connecting closely separated interfaces (Hendrickson, Weymouth & Yue Reference Hendrickson, Weymouth and Yue2020). In general, 
$T_s\, \varOmega (a_{res}) > 0.1$ avoids spurious events (Chan et al. Reference Chan, Dodd, Johnson and Moin2021; Gaylo et al. Reference Gaylo, Hendrickson and Yue2022). For our 
$\varepsilon$, this gives 
$T_s>0.4$. We perform ELA with 
$T_s=0.4$, 
$0.8$ and 
$1.6$. A posteriori, we see slight decreases in 
$Q_I$ and 
$Q_D$ as 
$T_s$ increases, as expected, due to missing temporally close entrainment and degassing events. We find that the shapes of 
$I(a)$ and 
$D(a)$ are independent of 
$T_s$, suggesting that spurious events are negligible. Thus we select the smallest 
$T_s=0.4$.
4.3. Results
Figure 4 shows instantaneous 
$\mathcal {D}$, the ratio of degassing flux to entrainment flux. The flow is consistently air-generating (
$\mathcal {D}<1$) up to 
$t=100$, after which there are some instances of 
$\mathcal {D}>1$. Consistent with the analysis in § 3.1, we see 
$\varGamma <1$ (
$\mathcal {D}<0.5$) for a short time near the start of entertainment (
$t<30$), before 
$\varGamma >1$ is obtained. We will focus on the air-generating period 
$t\in [40,70]$ when, in addition, the turbulence is quasi-steady (see figure 2), enabling us to normalize by turbulent scales. Unlike 
$\mathcal {D}\ll 1$ seen for breaking waves (Deane & Stokes Reference Deane and Stokes2002), our FST shows 
$\mathcal {D}\approx 0.7$ during this air-generating period, which appears quite independent of 
${{Fr}}^2$ (see table 1). In general, we expect gravity to affect degassing through two opposing mechanisms. First, for buoyancy-driven bubbles of a given radius, 
$\varLambda (a)$ decreases with decreasing 
${g}$. Second, the size of the largest bubble increases with decreasing 
${g}$ (Yu et al. Reference Yu, Hendrickson and Yue2020), and 
$\varLambda (a)$ increases with increasing 
$a$. Figure 4 suggests that these opposing mechanisms approximately cancel over a broad 
${{Fr}}^2$ range. In summary, we conclude that in air-generating FST, degassing is an important mechanism balancing entrainment, independent of 
${{Fr}}^2$.
Ratio of degassing flux to entrainment flux for: blue, 
${{Fr}}^2=5$; light green, 
${{Fr}}^2=8$; dark green, 
${{Fr}}^2=10$; brown, 
${{Fr}}^2=15$; red, 
${{Fr}}^2=20$. Dashed lines indicate 
$t\in [40,70]$. For clarity, we apply a top-hat filter (width 
$8$) denoted 
$\tilde {\cdot }$.
We now look at the entrainment and degassing size distributions. For entrainment, ELA provides direct measurement of 
$I(a)$ (see figure 5a). Although the scaling of 
$I(a)$ with 
${g}$ is not known precisely, to try to collapse the data, we normalize 
$I(a)$ by 
${{Fr}}_T^2$ in figure 5(b), as done in Yu et al. (Reference Yu, Hendrickson and Yue2020). Figure 5(b) suggests that the 
${{Fr}}_T^2$ scaling may not fully capture the effect of 
${g}$ on the magnitude of 
$I(a)$. Of interest here, figure 5 shows that 
$I(a)\propto a^{-10/3}$ given in Yu et al. (Reference Yu, Hendrickson and Yue2020) clearly captures the data, confirming (3.5).
Entrainment size distribution during 
$t\in [40,70]$ for: blue, 
${{Fr}}^2=5$; light green, 
${{Fr}}^2=8$; dark green, 
${{Fr}}^2=10$; brown, 
${{Fr}}^2=15$; red, 
${{Fr}}^2=20$. Dashed lines show 
$I(a)\propto a^{-10/3}$ predicted by Yu et al. (Reference Yu, Hendrickson and Yue2020).
For degassing, figure 6(a) shows that 
$\varLambda (a)/\varOmega (a) > 1$ for most bubble radii. This means that degassing is not only an important mechanism, but dominant over fragmentation. For small bubble radii where 
$\varLambda (a)/\varOmega (a)\lesssim 1$, degassing may still be dominant as this ratio as a measure of degassing versus fragmentation does not consider the cancellation between 
$S_f^+$ and 
$S_f^-$. For degassing-dominated bubble populations, our primary interest is the degassing power-law slope 
$\alpha$, used to determine the bubble size distribution. Here, ELA confirms the split power law predicted in § 3.2 (see figure 6b). In addition to capturing the scaling, we highlight that (3.8) predicts 
$\varLambda (a)$ within nearly an order of magnitude, and 
$a_\varLambda$ is within a factor of two of the observed transition radius. We note that we measure the quantity of interest, 
$\varLambda (a)$, rather than the characteristic scales 
$L_\varLambda$ or 
$U_\varLambda$ individually. In flows with concurrent entrainment and degassing over a range of turbulence scales, it is not clear what measures could independently quantify 
$L_\varLambda$ or 
$U_\varLambda$. Instead, the agreement between 
$\alpha$ predicted in § 3.2 and DNS gives us confidence in the proposed model. To quantify this agreement, we perform a least squares regression fit of the points in figure 6(b) to
\begin{equation} \log\left[\varLambda(a)/\varLambda_0\right] = \hat{\alpha} \log\left[a/a_\varLambda\right] + C, \end{equation}
for 
$a/a_\varLambda <1.5$ and 
$a/a_\varLambda >1.5$. Including 
$95\,\%$ confidence intervals, this gives 
$\hat {\alpha }=0.92\pm 0.14$ (versus 
$\alpha =1$) and 
$\hat {\alpha }=1.75\pm 0.38$ (versus 
$\alpha =2.5$), respectively. We note that there is a relatively small range of 
$a/a_\varLambda >1.5$, so regression may not be reliable (also indicated by the large confidence interval). Qualitatively, figure 6(b) is compatible with 
$\alpha (a>a_\varLambda )=2.5$.
Degassing rate during 
$t\in [40,70]$ for: blue, 
${{Fr}}^2=5$; light green, 
${{Fr}}^2=8$; dark green, 
${{Fr}}^2=10$; brown, 
${{Fr}}^2=15$; red, 
${{Fr}}^2=20$. Here, 
$\varOmega (a)$ is given by (2.3). In (b), black lines show 
$\varLambda (a\,{<}\, a_\varLambda )\,{\propto}\, a$ and 
$\varLambda (a>a_\varLambda )\propto a^{5/2}$.
Finally, we show the bubble size distribution 
$N(a)$ during the air-generating period (obtained independently of ELA) shown in figure 7. For consistency, we retain the 
${{Fr}}_T^2$ scaling, although as for 
$I(a)$ (cf. figure 5b), the effect of gravity on the magnitude seems not to be fully captured. In figure 7(a), we compare to 
$\beta =-10/3$ expected for a fragmentation-dominated population. In figure 7(b), we further normalize 
$N(a)$ to compare to 
$\beta$ expected for a degassing-dominated population. Recalling (3.4), degassing-dominated 
$N(a)$ depend on 
$I(a)$ as well as 
$\varLambda (a)$. Noting the similarities between the vertical spread in figures 7(b) and 5(b), we conclude that the weaker collapse of figure 7(b) compared to figure 6(b) is primarily the result of 
${{Fr}}_T^2$ scaling not fully capturing the effect of 
${g}$ on 
$I(a)$. Following an approach similar to that for for 
$\alpha$, we perform least squares regression to quantify 
$\beta$, although the included gravity effects limit the reliability of regression here. For 
$a/a_\varLambda <1.5$ we obtain 
$\hat {\beta }=-5.12\pm 0.17$ (versus 
$\beta =-4.\bar {3}$), and for 
$a/a_\varLambda >1.5$ we obtain 
$\hat {\beta }=-6.48 \pm 0.57$ (versus 
$\beta =-5.8\bar {3}$). Ultimately, comparing figures 7(a) and 7(b), it is clear that, as expected based on the ELA measurements of 
$\mathcal {D}$ and 
$\varLambda (a)/\varOmega (a)$, the power-law slope of 
$N(a)$ is much better described by the degassing-dominated model. In summary, we show that this air-generating FST creates a degassing-dominated bubble population with power-law slope given by (3.9), unique from fragmentation-dominated populations.
Average bulk bubble size distribution 
$N(a)$ during 
$t\in [40,70]$ for: blue, 
${{Fr}}^2=5$; light green, 
${{Fr}}^2=8$; dark green, 
${{Fr}}^2=10$; brown, 
${{Fr}}^2=15$; red, 
${{Fr}}^2=20$. Dashed lines show 
$N(a)\propto a^{-10/3}$ for fragmentation dominance. Solid lines show 
$N(a< a_\varLambda )\propto a^{-13/3}$ and 
$N(a>a_\varLambda )\propto a^{-35/6}$ for degassing dominance.
5. Discussion
The derivation of degassing-dominated 
$\beta$ in § 3.2 relies on the characteristic rise velocity 
$U_\varLambda$, for which we use a relatively simple model of fluctuation velocity 
$u_{rms}$ for small bubbles and the quiescent rise velocity 
$W_t$ for large bubbles. Here, we consider how alternative models for 
$U_\varLambda$ predict 
$\beta$.
For small bubbles in the turbulence-driven degassing regime (
$a \ll a_\varLambda$), Ruth et al. (Reference Ruth, Vernet, Perrard and Deike2021) suggest a terminal rise velocity based on a modification of the quiescent rise velocity, 
$U_\varLambda \sim W_t a^{1/2} {g}^{1/2} u_{rms}^{-1}$. Still using (3.6) for 
$W_t$, we repeat the analysis in § 3.2 with this alternative small-bubble 
$U_\varLambda$ model. This gives 
$U_\varLambda \propto a^2$ for 
$a< a_\varLambda$, and we obtain 
$\alpha (a< a_\varLambda )=3$ and 
$\beta (a< a_\varLambda )=-6.\bar {3}$. The model that we use for small bubbles (
$U_\varLambda \sim u_{rms}$) agrees better with the values 
$\hat {\alpha } = 0.92 \pm 0.14$ and 
$\hat {\beta } = -5.12 \pm 0.17$ measured from DNS.
Finally, we remark that our corroboration between theory and DNS is restricted to bubbles in the inertial regime (
$1<{{Re}}_b<100$) due to the limitations of DNS. The 
$W_t$ that we have used in § 3.2 correspond to this inertial regime. For larger physical-scale free-surface flows where 
${{Re}}$ is much larger than captured by our DNS (e.g. Hendrickson et al. Reference Hendrickson, Weymouth, Yu and Yue2019), we expect bubbles primarily in the 
${{Re}}_b>100$ regime. Using existing knowledge of quiescent terminal bubble rise velocity, the analysis in § 3.2 can be extended readily to higher 
${{Re}}_b>100$, where the quiescent rise velocity is 
$W_{t} = 0.721 \, \sqrt {{g} (2a)}$ (Davies & Taylor Reference Davies and Taylor1950; Park et al. Reference Park, Park, Lee and Lee2017) rather than (3.6). By repeating the analysis in § 3.2 with this large-
${{Re}}_b$ 
$W_t$, we obtain 
$a_\varLambda =0.96 \, u_{rms}^2 {g}^{-1}$. Because we have only modified 
$W_t$, turbulence-driven degassing remains unchanged, 
$\alpha (a< a_\varLambda )=1$, but for buoyancy-driven degassing, 
$\alpha (a>a_\varLambda )=3/2$. Thus compared to the 
$1<{{Re}}_b<100$ regime studied here, for larger physical-scale free-surface flows where 
${{Re}}_b>100$, we expect the same power law 
$\beta (a< a_\varLambda )=-4.\bar {3}$ for smaller bubbles (noting that 
$a_\varLambda$ is different), and a less negative power law 
$\beta (a> a_\varLambda )=-4.8\bar {3}$ for larger bubbles.
5.1. Surface tension effects
For breaking waves, it is observed that surface tension causes a departure from the power law 
$\beta =-10/3$ for bubbles smaller then the Hinze scale 
$a_H\sim (\sigma /\rho _w)^{3/5} \varepsilon ^{-2/5}$ (Deane & Stokes Reference Deane and Stokes2002). For these sub-Hinze scale bubbles, surface tension prevents fragmentation and the fragmentation cascade that would otherwise produce 
$\beta =-10/3$ for all 
$a$ (Garrett et al. Reference Garrett, Li and Farmer2000). Thus in breaking waves, the effect of surface tension on 
$\beta$ is through weakening fragmentation (reducing 
$\varOmega (a)$), creating a change in 
$\beta$ at the Hinze scale 
$a_H$. For FST (at large 
${{We}}$), we observe degassing dominance, 
$\varLambda (a)/\varOmega (a)\,{\gg}\,1$, meaning that the effect of fragmentation in general on 
$\beta$ is negligible (as shown in § 3). Reduction of 
$\varOmega (a)$ by surface tension would only increase 
$\varLambda (a)/\varOmega (a)\gg 1$. First, we note that this means that observing degassing dominance for FST at 
${{We}}\gg 1$ (as we do in § 4) implies degassing dominance independent of 
${{We}}$. More broadly, the Hinze scale is a characteristic scale related to fragmentation. For degassing-dominated flows, fragmentation is negligible, so the Hinze scale is not expected to be a relevant characteristic scale. Unlike breaking waves, we do not expect any change in 
$\beta$ at the Hinze scale.
For degassing-dominated flows, we expect the effects of surface tension on 
$\beta$ to be through the degassing rate 
$\varLambda (a)$ and the entrainment distribution 
$I(a)$, as described by (3.4). For degassing, Park et al. (Reference Park, Park, Lee and Lee2017) provide corrections to 
$W_t$ to account for surface tension, which could be incorporated into the model for degassing rate 
$\varLambda (a)$ in the same way as previously discussed 
$W_t$ for different 
${{Re}}_b$. For entrainment in FST, Yu et al. (Reference Yu, Hendrickson and Yue2020) suggest that surface tension creates a split power law for 
$I(a)$, with the change in power-law slope occurring at the capillary scale, 
$a_c= 1/2 \sqrt {\sigma /\rho _w g}$. For an air–water interface and Earth gravity, this corresponds to 
$a_c \approx 1.5\, \mathrm {mm}$; 
$I(a)\propto a^{-10/3}$, on which (3.5) is based, applies for 
$a>a_c$. To include bubbles smaller than this capillary scale in our model for degassing-dominated FST, a sub-capillary scale entrainment model, such as 
$I(a< a_c)\,{\propto}\,a^{-4/3}$ proposed by Yu et al. (Reference Yu, Hendrickson and Yue2020), would be necessary. Through (3.4), this suggests a corresponding change in 
$\beta$ at the capillary scale.
6. Conclusion
We study the effect of degassing on air-entraining flows during the air-generating period, when the total volume of bubbles is increasing. Understanding the bulk bubble size distribution 
$N(a)$ is of critical importance, and we focus on its power-law slope 
$\beta$. In air-generating flow, entrainment is clearly important, while the other mechanism(s) that balance entrainment determine 
$\beta$. For breaking waves at large 
${{We}}$, fragmentation is the dominant balancing mechanism, giving 
$\beta =-10/3$ (Garrett et al. Reference Garrett, Li and Farmer2000; Deane & Stokes Reference Deane and Stokes2002). In contrast, we consider large-
${{We}}$ flows where degassing is the dominant balancing mechanism – specifically, the canonical problem of air-generating flows due to strong underlying free-surface turbulence (FST) (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999; Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019). We model degassing-dominated bubble populations theoretically using the population balance equation (PBE), and computationally using DNS, in conjunction with Eulerian label advection (Gaylo et al. Reference Gaylo, Hendrickson and Yue2022) which provides direct measurements of entrainment and degassing separately.
From PBE, we derive metrics to determine the importance of degassing in air-generating flows: the ratio of degassing flux to entrainment flux, 
$\mathcal {D}=Q_D/Q_I$, which quantifies the relevance of degassing; and the ratio of degassing rate to fragmentation rate, 
$\varLambda (a)/\varOmega (a)$, which determines whether degassing or fragmentation is dominant. For FST, DNS give 
$\mathcal {D}\approx 0.7$ over a broad range of air-entraining 
${{Fr}}^2$, consistent with degassing being a relevant mechanism; and 
$\varLambda (a)/\varOmega (a)>\operatorname {O}(1)$, indicating that degassing is dominant over fragmentation.
Also from PBE, we show that for degassing-dominated populations, 
$N(a)\approx I(a)/\varLambda (a)$, and DNS confirm 
$I(a)\propto a^{-10/3}$ predicted by Yu et al. (Reference Yu, Hendrickson and Yue2020). Assuming a power-law degassing rate 
$\varLambda (a)\propto a^{\alpha }$, we obtain the slope 
$\alpha$ for FST. Separated by the critical radius 
$a_\varLambda$ (see (3.7)), we find a turbulence-driven regime where 
$\alpha (a< a_\varLambda )=1$, and a buoyancy-driven regime where 
$\alpha (a>a_\varLambda )=5/2$ (for moderate bubble Reynolds number 
${{Re}}_b$). The DNS measurements of 
$\varLambda (a)$ confirm this split power-law scaling, which can be extended theoretically to large 
${{Re}}_b$. Based on this, we show that degassing-dominated 
$N(a)$ follow a split power law 
$\beta (a< a_\varLambda )=-4.\bar {3}$ and 
$\beta (a>a_\varLambda )=-5.8\bar {3}$ (for moderate 
${{Re}}_b$), which is confirmed independently by DNS.
We have identified that it is possible for an air-generating (implying large 
${{Fr}}$) free-surface flow to be degassing-dominated, even at large 
${{We}}$ when fragmentation is strongest. More study of the variety of air-generating free-surface flows is necessary to allow general prediction of degassing versus fragmentation dominance, and test our hypothesis that it depends on the presence of a mean downward flow. It also remains to identify and describe transitional flows where the effects of both fragmentation and degassing are comparable. Given an air-generating free-surface flow, the metrics for the importance of degassing, 
$\mathcal {D}$ and 
$\varLambda /\varOmega$, can be difficult to quantify, but can be inferred from the size distribution power-law slope(s) 
$\beta$. For breaking waves, the wide agreement on 
$\beta =-10/3$ (e.g. Deike Reference Deike2022) suggests that such flows are generally fragmentation-dominated. The FST that we consider is one example of an air-generating flow that isdegassing- rather than fragmentation-dominated. For another possible example, Hendrickson et al. (Reference Hendrickson, Weymouth, Yu and Yue2019) reported 
$-5<\beta <-4$ in the rooster-tail and diverging-wave regions of the flow behind a dry transom stern, indicative of the degassing-dominated flows considered in the present work.
Funding
This work was funded by the US Office of Naval Research grants N00014-20-1-2059 and N00014-24-1-2076 under the guidance of Dr W.-M. Lin. The computational resources were funded through the Department of Defense High Performance Computing Modernization Program.
Declaration of interests
The authors report no conflict of interest.
Appendix. Verification of grid independence
In addition to the grid convergence study performed by Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) to confirm that the turbulence and bubble population are sufficiently resolved by the grid described in § 4.1, we perform a convergence study on the values that we measure for our results to confirm that the DNS sufficiently resolves entrainment and degassing in FST. We perform a set of three simulations at 
${{Fr}}^2=15$ using the same method described in § 4.1, but with a finer grid of 
$576^2 \times 384$. This gives 
$\varDelta _{576} \approx 0.018$, versus 
$\varDelta _{384}\approx 0.027$ from § 4.1. For this simulation, we include all bubbles of radius larger than 
$1.5 \varDelta _{576}$. We measure (with 95 % confidence intervals) 
${\langle k \rangle _\delta } \times 10^3= 8.2 \pm 3.5$ and 
${\langle \varepsilon \rangle _\delta } \times 10^4= 6.1 \pm 2.7$ during 
$t\in [40,70]$, consistent with 
${{Fr}}^2=15$ using 
$\varDelta _{384}$ (see table 1).
Figure 8 compares the measured entertainment size distribution 
$I(a)$, degassing rate 
$\varLambda (a)$, and bulk bubble size distribution 
$N(a)$ between the two grid resolutions. These correspond to figures 5(a), 6(a) and 7(a), and we note that the range of the horizontal axis has been extended in figure 8. We have run three simulations with 
$\varDelta _{576}$ (versus six with 
$\varDelta _{384}$), so there is more statistical variability, particularly for larger bubbles, of which fewer are observed in a given simulation. From figure 8, it is clear, especially for the smallest bubbles where resolution would be a concern, that the results are consistent between the 
$\varDelta _{576}$ and 
$\varDelta _{384}$ grids. We conclude that DNS with the 
$\varDelta _{384}$ grid described in § 4.1 sufficiently resolves the relevant physics.
(a) Entrainment size distribution, (b) degassing rate, and (c) average bulk bubble size distribution during 
$t\in [40,70]$ for 
${{Fr}}^2=15$, with brown dots indicating 
$\varDelta _{384}$, and black squares indicating 
$\varDelta _{576}$. For consistency, turbulence values from 
$\varDelta _{384}$ are used to calculate 
$\varOmega (a)$ and 
${{Fr}}^2_T$.
