Scale separation and dependence of entrainment bubble-size distribution in free-surface turbulence

We consider the size spectrum of entrained bubbles under strong free-surface turbulence (SFST). We investigate the entrainment bubble-size spectrum per unit (mean) interface area, ${\mathcal{N}}_{a}^{E}(r)$, with dimension length$^{-3}$, and develop a physical/mechanistic model for ${\mathcal{N}}_{a}^{E}(r)$ through energy arguments. The model obtains two distinct regimes of ${\mathcal{N}}_{a}^{E}(r)$, separated by bubble-size scale $r_{0}$. For bubble radius $r>r_{0}$, the effects of gravity $g$ dominate those of the surface tension force $\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C}$, and ${\mathcal{N}}_{a}^{E}(r)\propto g^{-1}\unicode[STIX]{x1D716}^{2/3}r^{-10/3}$, where $\unicode[STIX]{x1D716}$ is the turbulence dissipation rate. For $r<r_{0}$, surface tension is more important and ${\mathcal{N}}_{a}^{E}(r)\propto (\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C})^{-1}\unicode[STIX]{x1D716}^{2/3}r^{-4/3}$. From the model, we show that $r_{0}\approx r_{c}=1/2\sqrt{\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C}g}$, the capillary length scale, and not the generally assumed Hinze scale $r_{H}$. For an air–water interface and Earth gravity, $r_{c}\approx$ 1.5 mm. The model provides an $\unicode[STIX]{x1D716}$–$r$ entrainment regime map that identifies a critical dissipation rate $\unicode[STIX]{x1D716}_{cr}$ (constant for given $g$ and $\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C}$) above which there is appreciable air entrainment, thus separating SFST and weak FST. We confirm the theoretical model and its predictions using two-phase, high-fidelity direct numerical simulations of a canonical FST flow using the conservative volume-of-fluid method: the respective power laws of ${\mathcal{N}}_{a}^{E}(r)\propto r^{-10/3}$ and $r^{-4/3}$ for $r>r_{0}$ and $r<r_{0}$; the value $r_{0}=r_{c}$; the scaling ${\mathcal{N}}_{a}^{E}(r)\propto \unicode[STIX]{x1D716}^{2/3}$; and the predictions of the $\unicode[STIX]{x1D716}$–$r$ entrainment regime map.


Introduction
Air entrainment occurs and plays important roles in both natural processes and engineering applications. Quantifying the total volume, size distribution and scale dependence of entrained air in free-surface turbulence is key to modelling the natural air entrainment and gas transfer across the air-sea interface in ocean breaking waves (Deane & Stokes 2002) and the design of aeration cascades to increase oxygen concentration in water treatment plants (Chanson 1996).
To date, much of the understanding of air entrainment is in the context of breaking waves. Deane & Stokes (2002) measured the bulk bubble-size spectra (per unit volume) N(r) inside breaking waves in the laboratory. In the initial acoustic phase, when entrainment primarily occurs, they identified two distinct power-law regimes of N(r), depending on the bubble size r. For r 1 mm, N(r) ∝ r −10/3 , while for r 1 mm, N(r) ∝ r −3/2 . They calculated the Hinze scale r H (Hinze 1955) using r H = 2 −8/5 (We c σ /ρ) 3/5 −2/5 , (1.1) where We c ≈ 4.7 is the critical Weber number for bubble breakup by turbulence predicted from experiments (Lewis 1982;Martínez-Bazán, Montanes & Lasheras 1999). Using the fact that r H in the experiments was close to the separation scale r 0 ∼ 1 mm, they argued that r 0 is r H . In the quiescent phase (∼1.5 s after the acoustic phase) the bubble plume evolved rapidly under the influence of advection, turbulent diffusion and buoyant degassing, with an accompanying (more complex) evolution of N(r). For the air entrainment phase, other researchers (Loewen, O'Dor & Skafel 1996;Rojas & Loewen 2007;Blenkinsopp & Chaplin 2010;Deike, Melville & Popinet 2016;Wang, Yang & Stern 2016) confirm the power-law dependence N(r) ∝ r −10/3 for larger bubbles in breaking waves experiments and simulations. For the regime r r 0 , Rojas & Loewen (2007) and Wang et al. (2016) obtained similar results to those of Deane & Stokes (2002). Several parameterizations have been proposed for the air entrainment phase N(r). Garrett, Li & Farmer (2000) performed a dimensional analysis in which N(r) (dimension L −4 ) is assumed to depend on the turbulence dissipation rate , the average rate of supply of air Q and the bubble radius r, and obtained a parameterization of the form N(r) ∝ Q −1/3 r −10/3 .
(1.2) Deane & Stokes (2002) pointed out that (1.2) is valid only for r r H . For r r H , surface tension becomes important, requiring a different scale dependence. Performing a similar dimensional analysis for r < r H , and assuming N(r) is a multiplicative function of surface tension σ , liquid density ρ, jet velocity v and bubble radius r, they obtain a N(r) scaling for these smaller bubbles: The models (1.2) and (1.3) originate from dimensional analysis and it is not clear if they reflect all of the underlying physical mechanism(s) of air entrainment. In particular, the inverse dependence on in (1.2) and the non-dependence on in (1.3) of N(r) seem non-physical. Since turbulence provides energy for bubble entrainment (Brocchini & Peregrine 2001), a greater should increase the magnitude of N(r). Deike et al. (2016) argued that Q should itself be a function of , and the combination of Q and −1/3 in (1.2) could lead to a positive exponent for . Similarly for (1.3), we can obtain a scaling of 2/3 instead of 0 if ∝ v 3 /D j , with D j the jet diameter. Finally, the identification of the regime separation scale r 0 with the Hinze scale r H , while physically plausible, is not definitive. The possibility that r 0 derives from other relevant physical scales of the problem cannot be ruled out.
In this paper, we investigate the air entrainment across a free surface induced by underlying strong free-surface turbulence (SFST), specifically its bubble-size distribution and scale dependence. We define N E a (r) (dimension L −3 ) as the bubble-size spectrum per unit mean interface area (over a unit time in the quasi-steady air entrainment process). Through energy arguments, we show (in § 2) that N E a (r) follows two power-law regimes, depending on bubble size r, separated by a bubble-size length scale r 0 . For r > r 0 , N E a (r) scales as g −1 2/3 r −10/3 , describing a gravity-dominated entrainment regime. For r < r 0 , N E a (r) scales as (σ /ρ) −1 2/3 r −4/3 , describing a surface-tension-dominated entrainment regime. We show that r 0 is in fact the capillary scale r c = 1/2 √ σ /ρg, independent of flow parameters such as ; and hence not related to the Hinze scale r H as is often assumed (Deane & Stokes 2002;Wang et al. 2016). Inspired by the diagram of the (L, q)-plane in Brocchini & Peregrine (2001), we derive a (r, )-plane prediction of when air entrainment occurs in free-surface turbulence (FST). The resulting regime map predicts a critical dissipation value cr below which air entrainment is suppressed, thus demarcating weak FST (WFST) from SFST in terms of air entrainment. Like r 0 , cr depends only on the physical constants g and σ /ρ. In § 3, we confirm our theoretical predictions using high-resolution direct numerical simulations (DNS) of canonical FST flows, solving three-dimensional two-phase incompressible Navier-Stokes equations with a fully nonlinear interface resolved by the conservative volume-of-fluid method (cVOF) (Weymouth & Yue 2010). We note that the surface entrainment size spectrum N E a (r) has a direct relationship to the (bulk) volume size spectrum N a (r) generally measured in experiments and simulations by labelling and sizing bubbles in the bulk. For a relatively short time after the onset of entrainment, before subsurface mechanisms come into play, N a (r) closely resembles N E a (r) in terms of the power-law regimes and slopes. With increasing time, N a (r) evolves (primarily) as a result of turbulent bubble fragmentation and degassing, and can become qualitatively different from N E a (r) ). We perform DNS over a range of Froude (Fr) and Weber (We) numbers for the FST and obtain N a (r) at a relatively early time. The DNS confirms the key predictions of the scaling model: the gravity-and surface-tension-dominated power-law regimes and slopes; the separation scale r 0 = r c (≈1.5 mm for an air-water interface in Earth gravity); and the scaling with 2/3 (over the range we considered). Plotting the DNS data on the -r regime map also confirms the critical dissipation for air entrainment, cr ∼ O(10 −2 m 2 s −3 ), for an air-water interface and Earth gravity.

Scaling for bubble-size distribution of air entrainment in SFST
To model the size spectrum of entrainment across a given free-surface area, we consider the entrainment size spectrum N E a (r), of dimension L −3 , where N E a (r )δr is the number of bubbles of (equivalent) radius r < r < r + δr entrained per unit (mean) free-surface area. For simplicity, we assume quasi-steady air entrainment over unit time. A physical/mechanistic derivation for the scaling of N E a (r) through an energy argument requires two assumptions. First, we assume that the only source of bubble generation is the surrounding turbulence (Brocchini & Peregrine 2001). Second, we assume that the void fraction of bubbles is small and does not affect the dynamics of the bulk water turbulence (Garrett et al. 2000).
Consider a simplified entrainment scenario of a static spherical air bubble of radius r formed at a depth 2rα under the free surface (a change in energy between Scenario 1 and 2 sketched in figure 1a). We assume the constant α is of order 1 to account for this simplification. The minimum energy required to entrain such a bubble, E b (r), is the sum of the potential energy and the surface energy: The total number of bubbles of radius r entrained across a (initial) surface area A h is N E a (r)A h dr and the total energy required for entrainment of these bubbles is Using our first assumption, this amount of energy is supplied by turbulence over a near-surface region of depth r (around the entrained bubbles). Defining E(k) as the turbulence energy spectrum, the turbulence energy per unit volume at length scale 1/k is given by |E(k) dk|. We argue that the turbulence of length scale 1/k is primarily responsible for the entrainment of bubbles of radius r ∼ 1/k. The available energy supplied by near-surface turbulence for entrainment within a volume A h r is proportional to For SFST, the characteristic large surface deformations (and large Fr and We) weaken the free-surface boundary conditions, resulting in essentially isotropic near-surface turbulence . Thus, the Kolmogorov spectrum E(r −1 ) ∝ 2/3 r 5/3 provides the near-surface turbulence energy spectrum (this is also confirmed in the DNS, see figure 4). Using the Kolmogorov spectrum, (2.4) becomes For flows with different underlying turbulence spectra (for example, due to large bubble density, Alméras et al. 2017), the predicted bubble-size spectra will be modified accordingly. From (2.5), we observe two regimes of entrainment by considering f (r) = (2/3)grα + (σ /ρ)(1/r). For large r, gravity dominates f (r) as (2/3)grα is large, The scale separating the two regimes r 0 occurs when the two terms in f (r) are comparable at r 0 = √ 1.5σ /αρg. This capillary scale r 0 ≈ r c = 0.5 √ σ /ρg corresponds to a balance between gravity and surface tension forces, or Bond number Bo ≡ ρg(2r c ) 2 /σ = 1. We note that (2.6) and (2.7) are also obtained directly from dimensional analysis using the relevant set of physical parameters. In contrast to Garrett et al. (2000) and Deane & Stokes (2002), for free-surface air entrainment, the size spectrum per unit surface area, N E a (r) (dimension L −3 ), is powered by the FST kinetic energy measured by (L 2 T −3 ). This is balanced by gravity g (LT −2 ) and surface tension σ /ρ (L 3 T −2 ). Matching dimensions, we immediately obtain the scaling results for the respective size regimes.
Several remarks are in order. First, for the r > r 0 regime, the power-law slope r −10/3 happens to coincide with that of Garrett et al. (2000) -for example, (1.2). However, the underlying physical processes in these two models are different in that Garrett et al. (2000) considered a cascade scenario of bubble fragmentation rather than air entrainment driven by surface turbulence as in our model. In fact, the slopes of the two power-law regimes in (2.6) and (2.7) result directly from the power-law slope of the near-surface isotropic SFST energy spectrum. One physical explanation for the coincidental similarity of the power-law slope of our model and Garrett et al. (2000) is to frame the air-water boundary of an initially large cavity immersed in an isotropic turbulent field as a free surface with underlying isotropic near-surface turbulence. In this context, the bubble generation process through fragmentation is similar to turbulence air entrainment . Second, for r < r 0 in the surface-tension-dominated regime, the power-law slope in (2.7) differs (slightly) from the model of Deane & Stokes (2002) in (1.3). Third, the (positive exponent) dependence on turbulence dissipation rate ∼ 2/3 in the current models (2.6) and (2.7) is physically reasonable compared to that of (1.2) (∼ −1/3 ) and (1.3) (∼ 0 ). In the current model, stronger surface turbulence leads to more entrainment. Finally, the separation scale r 0 is independent of flow parameters and not strictly related to the Hinze scale r H as generally assumed (without strict theoretical consideration). The present theory provides a direct argument that r 0 is the physical capillary scale r 0 ≈ r c , which, for an air-water interface under the influence of Earth gravity, is r c ≈ 1.5 mm.
As discussed in § 1, the scale dependence of the surface entrainment size spectrum N E a (r) should be reflected in the underlying volume size spectrum N a (r) providing it is early enough in time, when fragmentation and degassing are not a factor. Figure 2 shows the theoretical prediction for N E a (r) in (2.5) superposed on available experimental and DNS data for N a (r) obtained for both breaking waves and SFST flows. For comparison of the spectral shape/slope, we normalize the magnitude of each dataset by its values at some (arbitrary) radius r = 2 mm. Despite the fact that most of the datasets in figure 2 are from breaking waves, the scale dependence of the data appears to be well described by (2.5) for the entire range of r. We point out that in figure 2 the transition between the two regimes for most of the data occurs near r 0 = 1-2 mm, despite presumably different values of turbulent dissipation rates for the different datasets.
Air entrainment in SFST is a result of two competing effects on the free surface: the disrupting effect of near-surface turbulence and the stabilizing effect from gravity and surface tension. Comparing the two effects provides criteria for the lower and upper bounds r l and r u of bubble size that can be entrained. We argue that the lower bound r l corresponds to a scale where turbulence disturbance is much smaller than surface tension forces -that is, We t = ρv 2 (2r l )/σ = We l ( 1), where v 2 ≈ 2.0( d) 2/3 is the turbulence fluctuation over a distance d = 2r l ; r l can be further written as r l = 2 −8/5 (σ /ρ We l ) 3/5 −2/5 .
(2.8) Equation (2.8) is the same form as the Hinze scale (1.1), with We l replacing We c . Both r l and r H correspond to the critical length scales resulting from the competition between near-surface turbulence and surface tension forces. While the length scale r H is for bubbles with a stable spherical interface, r l is the length scale associated with the breakup of a free surface with a relatively unstable flat interface. As such, we expect We l We c (and consequently r l r H ). We define the upper bound bubble length scale r u similarly from the (turbulence) Froude number Fr 2 t = v 2 /g(2r u ) = Fr 2 u ( 1), which yields r u = 4(Fr 2 u g) −3 2 .
With the expressions (2.8) and (2.9) in terms of , we obtain a critical turbulent dissipation cr when r u ( cr ) = r l ( cr ), below which no entrainment would occur: cr = C g 5/4 (σ /ρ) 1/4 . (2.10) As with the separation scale r 0 for N E a (r), cr is a constant that depends only on the interface media and value of gravity. The upper and lower bounds (2.8) and (2.9) enable us to predict an -r regime map for bubble entrainment in FST (see figure 3, also see Brocchini & Peregrine (2001)). The dependencies of r l and r u on demarcate four regions in the ( -r)-plane. The strong free-surface turbulence regime ( > cr ) is where entrainment occurs within the size range r l < r < r u , which we denote Region 2 in figure 3. The bubble-size range increases with . Turbulence will not entrain bubbles outside of this size range due to dominant gravity or surface tension forces. The weak free-surface turbulence regime ( < cr ) is where no entrainment occurs for three different physical reasons. In Region 1, r < r l and surface tension forces stabilize the free surface. In Region 3, r > r u and strong gravity forces suppress air entrainment. In Region 0, both gravity and surface tension forces suppress entrainment.

Direct numerical simulations of FST
We perform DNS of a canonical problem of three-dimensional incompressible two-phase (air and water) viscous turbulent flow with turbulent kinetic energy supplied by a sheared underlying bulk water flow (see figure 1b). Extensive documentation of the turbulence characteristics of this flow exists for non-entraining WFST at low For the data presented here, the mean shear-flow depth L and velocity deficit U are the characteristic length and velocity scales and normalize all the variables unless otherwise stated, with Fr 2 = U 2 /gL, We = ρU 2 L/σ and Re = UL/ν. The domain size is 6 3 and the grid size ∆ ≈ 0.01. In the analyses below, we report (only) bubbles that are considered resolved in the presence of surface tension in the DNS corresponding to equivalent diameter 2r 7∆. The unreported entrained volume is 1 % for all cases. For SFST, we present results for two high-resolution DNS, both with Bo = ρgL 2 /σ = 100. DNS-1 has parameters Fr 2 = 21, We = 2100 and Re = 1200 and DNS-2 has Fr 2 = 24.64, We = 2464 and Re = 1300. We note that for constant g, σ /ρ and ν, DNS-2 has a larger U and the same L as DNS-1. The DNS results capture the evolution of air-entraining SFST corresponding to large Fr and We . Briefly, starting from a quiescent free surface, the FST grows rapidly due to production by the bulk flow mean shear, followed by a quasi-steady active entrainment SFST period, after which the FST eventually decays (and entrainment ceases) due to turbulence diffusion and dissipation. We report results within the quasi-steady active entrainment SFST period. In this period, turbulence forms vortical structures of different length scales near the free surface and entrains bubbles of different sizes . Figure 4 shows the one-dimensional energy spectra of near-surface turbulence during the SFST period for DNS-1 (DNS-2 is similar). Consistent with the assumption in § 2, the turbulence shows isotropy, with the three components having comparable magnitudes and following the Kolmogorov k −5/3 spectral slope. Figure 5 shows the volume bubble-size spectra N a (r) for DNS-1 during the active entrainment period (r is normalized by the capillary scale r c for clarity). The measured bulk spectrum provides direct support of the scale dependencies of (2.6) and (2.7), respectively, for r > r 0 and r < r 0 with the separation scale r 0 ≈ r c . We estimate r H (shown in figure 5) using the dimensionless form of (1.1), r H = 2 −8/5 (We c /We) 3/5 −2/5 , using the DNS-1 near-surface turbulence dissipation rate ∼ O(10 −3 ) and We c = 4.7. For DNS-1, r H is approximately three times greater than r c , and this scale separation allows us to confirm that the separation scale between the two regimes is indeed r 0 ≈ r c and not r H . For this particular DNS case, turbulence fragmentation of entrained bubbles should not have significant effects on N a (r) within the entrainment period, because r H is at the upper bound of the measured entrainment. We confirm this by comparing the three spectra at different times for r < r H in figure 5. This enables us to estimate N E a (r) using N a (r) with confidence. The inset of figure 5 shows N a (r) for both DNS-1 and DNS-2 averaged over the entrainment period. For DNS-2, is increased by a factor of 1.5. However, r 0 is the same for both DNS datasets and occurs at r c (which is the same for both DNS). This further confirms that r 0 between the two power-law entrainment regimes is in fact the capillary scale r c and not the Hinze scale r H , which changes with .
To confirm the scaling of the entrainment model with turbulent dissipation, we perform a series of lower-resolution DNS with ∆ = 0.027 for FST with Fr 2 = 10, We = 50 000 and Re = 1000, varying the slope of the initial mean shear profile to achieve different values for . Figure 6 shows the linear dependence of the integral of averaged N a (r) (or entrained volume) plotted against 2/3 (averaged during the entrainment period), which confirms the positive exponent scaling of N E a (r) ∝ 2/3 in (2.5). Linear regression obtains a horizontal intercept corresponding to 2/3 ≈ 0.002. When scaled by an air-water interface under Earth gravity, it provides an independent Finally, we confirm the prediction of the -r entrainment regime map, figure 3, using FST DNS over a broad range of Fr 2 , We and (by varying the initial mean shear profiles). We estimate the physical units of the characteristic length scale L and near-surface turbulence dissipation rate using the defined values of Fr 2 , We and Re (and assume Earth gravity and an air-water interface). The measured size range of entrainment as functions of for the different cases are superimposed onto the map in figure 3 (heavy lines). Notice that r l and r u in (2.8) and (2.9) assume an unbounded turbulence inertial range. For each DNS case, the range of the entrainment bubble size has a lower bound of ∆ and an upper bound of L as well. This range of (∆, L) is represented by dashed lines. The selected DNS cases with/without air entrainment correctly fall into the predicted regions of strong/weak FST, substantially confirming the predictions of § 2.

Conclusion
We perform theoretical and numerical investigations of air entrainment, and specifically the entrainment bubble-size distribution under SFST. Using arguments relating the energy required for bubble entrainment/formation and available turbulent kinetic energy in the SFST, we obtain a theoretical model for the surface entrainment bubble-size spectrum N E a (r). The model describes two size regimes separated by a scale r 0 , with N E a (r) ∝ g −1 2/3 r −10/3 for r > r 0 , and N E a (r) ∝ (σ /ρ) −1 2/3 r −4/3 for r < r 0 , where r is the (equivalent) bubble radius, g is gravity, σ is the surface tension, ρ is the water density and is the turbulent dissipation rate of the SFST. From the model, we obtain that r 0 is the capillary scale r c (≈1.5 mm for an air-water interface under Earth gravity), distinct from the Hinze scale r H . We also propose an -r regime map to describe air entrainment in free-surface turbulent flows. The map predicts a critical value of the dissipation rate cr ≈ 0.02 m 2 s −3 above which air entrainment occurs. High-resolution volume-conserving two-phase DNS of canonical shear-flow FST confirm the theoretical model and predictions. Namely, it provides direct evidence of the two power-law regimes and slopes, the separation scale r 0 ≈ r c , the positive exponent of the turbulent dissipation 2/3 , and the -r entrainment regime.
We remark that, while the (shear-flow) FST we consider is canonical and has been studied extensively theoretically and computationally, there exist few detailed experimental investigations of SFST air entrainment. Direct experimental measurements to support the present theoretical and DNS results would be invaluable. In particular, it would be useful to confirm the scale separation we find for additional combinations of ρ, σ and . The energy argument and predictions for SFST air entrainment may provide insights into more complex air-entraining, free-surface problems, such as wave breaking (see figure 2) and ship wakes . These are subjects of our current investigations.