2.1 Operational details of Nortek Vectrino II

The Vectrino II profiling ADV was introduced by Nortek in 2011 (Craig et al. Reference Craig, Loadman, Clement, Rusello and Siegel2011). Our work using the instrument in bubbly flows reveals the following. Pulse interference is common because of a tight bubble core, and we enabled the adaptive ping algorithm in the instrument software to eliminate these weak spots. To ensure the highest possible data quality, two precautionary steps were taken in the experiments. First, the water tank was adequately seeded with neutrally buoyant $50~\unicode[STIX]{x03BC}\text{m}$ polyamide particles such that the background (without bubble plume flow) signal-to-noise ratio (SNR) and correlation coefficient (Corr) were better than 30 and 98 %, and that their profiles were parabolic with a peak at the sweet spot of the measurement strip (Nortek 2013). Second, a probe check was performed with the plume flow turned on before each experiment run. The amplitude profiles of the four receivers should overlap and be smooth across the measuring strip, and again with a distinct maximum at the strip sweet spot. A schematic sketch of the thirty-five sampling disks is shown in figure 3. The probe’s sweet spot is approximately 5 cm from the transmitter head and corresponds to the tenth cell in the measurement strip. In this study, only data from the sixth to sixteenth cells are considered due to the decaying sensitivity, hence accuracy, at both ends of the strip. These eleven cells were repeatedly used to sample the plume flow at different radial positions. The data quality is discussed in detail in the result section § 4.

Figure 3.

Characteristics of the sampling disks of the Vectrino II; showing the varying measurement sensitivity and the positions of data considered in present study.

2.2 Post-processing of raw ADV data

In the ADV community, the term despiking refers to the removal of faulty data in a velocity time series (Goring & Nikora Reference Goring and Nikora2002). Faulty data usually appear as extreme values and are removed by either thresholding on the SNR/Corr values, or more commonly, by statistical methods that better respond to flow physics. In recent literature on despiking, the class of phase-space despiking methods has been suggested as an objective rejection algorithm, first advocated by Goring & Nikora (Reference Goring and Nikora2002). The method is based on a classical statistical result of the normal (Gaussian) distribution: extreme values cannot appear more often than a certain limit. The presumed Gaussian behaviour appears to apply in many turbulent fluid flows although non-Gaussianity of the small dissipative scales is recognized (Pope Reference Pope2000). More importantly, this despiking method assumes that there is only one forcing responsible for the data, which is true for flows characterized by a single length scale, e.g. a single-phase jet. In multiscale flows such as bubble plumes, there can be multiple forcings. Each may have a different PDF and contribute a different data fraction. Blindly despiking such data will lead to the removal of physical data. Cea, Puertas & Pena (Reference Cea, Puertas and Pena2007) recognized the issue and applied a modified phase-space algorithm to a dataset obtained in the highly aerated flow behind a sluice gate. We will show in § 3 that our bubble plume data do not need despiking, and we will demonstrate that the extreme values are signatures of liquid fluctuations created by the rising bubbles.

2.3 Noise-reduced estimates of the Reynolds stress tensor

It is known that raw ADV data obtained from flows without any physical obstructions are accurate for computing the mean flow field while turbulent statistics are highly biased by the Doppler noise which must be estimated and removed from measurements, see e.g. Voulgaris & Trowbridge (Reference Voulgaris and Trowbridge1998) and Garcia et al. (Reference Garcia, Cantero, Nino and Garcia2005). For a four-receiver bistatic probe like the Vectrino II, the redundant (along transmitter direction) radial velocity $\tilde{u} 2$ can be used together with $\tilde{u} 1$ to directly estimate the Doppler noise. Noise-reduced estimates of the Reynolds stress tensor can then be obtained after accounting for probe geometry (Hurther & Lemmin Reference Hurther and Lemmin2001). Details of the correction procedure are given in appendix A.

3.1 Empirical probability density functions

Figure 4 shows a six minute segment (out of a full 15 min record) of the raw streamwise velocity $\tilde{w}$ at the plume core $r=0$ and at $z/D=6.6$ for a single sample volume along the strip for Case A. When viewed over the entire segment, the series is stationary in time but is also interlaced with intermittent large negative and positive values. This means that the series is not locally stationary, i.e. over small time intervals. Fluctuations of much higher frequencies than those present during typical times occur when $\tilde{w}$ attains extreme values.

Figure 4.

Raw data for streamwise velocity $\tilde{w}$ at plume core $r=0$ (Case A, $z/D=6.6$ ).

Figure 5.

Empirical PDF of raw velocities at the plume core $r=0$ (Case A, $z/D=6.6$ ): (a) velocities $\tilde{w}$ and $\tilde{v}$ in the Cartesian reference frame; (b) the two independent measurements $\tilde{u} 1$ and $\tilde{u} 2$ of the radial velocity.

The corresponding empirical PDFs are shown in figure 5. Two characteristics are observed for $\tilde{w}$ : (i) a primary peak at approximately $0.18~\text{m}~\text{s}^{-1}$ , characterized by a distribution of positive velocities that is right skewed; and (ii) a broad distribution of negative velocities (mean of all negative values is $-0.35~\text{m}~\text{s}^{-1}$ ) but having low probability densities. These large negative velocities do not occur due to phase wrapping of the ADV as the same values result when doubling the velocity range of the signal processor. The interpretation of these observations is analysed in detail in § 3.2.

In contrast to the PDF for $\tilde{w}$ , the out-of-plane velocity $\tilde{v}$ has only one peak in its PDF. The PDF is symmetric with respect to zero velocity, as can be expected from the flow symmetry at the plume core. Similar observations can be found from the PDFs of the radial velocities $\tilde{u} 1$ and  $\tilde{u} 2$ . Hence, the rapid vertical rise of air bubbles created definite signatures in $\tilde{w}$ only.

Figure 6.

Empirical PDF of the three velocity components at the bubble plume centreline ( $r=0$ ): Case A at $z/D=6.6$ (a,b) and Case B at $z/D=2.2$ (c,d).

Figure 7.

Empirical PDF of the three velocity components at the edge of the bubble plume ( $r=2b_{g}$ ): Case A at $z/D=6.6$ (a,b) and Case B at $z/D=2.2$ (c,d).

Figures 6 and 7 compare the velocity PDFs of our bubble plumes with those of homogeneous bubble swarms reported in Riboux et al. (Reference Riboux, Risso and Legendre2010), Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016) and Almeras et al. (Reference Almeras, Mathai, Lohse and Sun2017). In these figures, PDFs of the standardized velocity $(\tilde{u} _{i}-U_{i})/u_{i,rms}$ (rms, root-mean-square) are plotted, where $i$ is the coordinate axis in indicial notation. This standardization facilitates the comparison among flows with different void fractions, and hence, different magnitudes of kinetic energy of the induced liquid fluctuations. The figures are also plotted with a semi-log ordinate axis to show the stretched exponential tails of the PDFs (Risso Reference Risso2018); PDFs having the form $p(x)\sim \exp (-|x|)$ will appear as straight lines in such a plot.

At the bubble plume centreline $r=0$ , figure 6 shows that for both Cases A and B the streamwise (vertical) liquid fluctuations have very similar PDFs to those of bubble swarms. All PDFs are positively skewed, and fluctuations over the range $2<(\tilde{w}-W)/w_{rms}<6$ have higher occurrence probabilities than that given by the standard Gaussian PDF (black dashed line). In fact, our data almost overlap exactly with the LDA data reported in Riboux et al. (Reference Riboux, Risso and Legendre2010). For the two horizontal velocity components ( $\tilde{u} ,\tilde{v}$ ), the agreement of the PDFs between the bubble plume data and the bubble swarms data (Riboux et al. Reference Riboux, Risso and Legendre2010) is again very good. Unlike the PDF of $\tilde{w}$ , these PDFs are symmetric. Finally, all velocity components exhibit multiple, stretched exponential tails, an observation already reported for bubble swarms (Risso Reference Risso2018). From the comparisons just made, two conclusions on bubble-induced liquid fluctuations can be drawn: (i) rising bubbles create a specific pattern of liquid fluctuations that are insensitive to the details of the flow (e.g. heterogeneity, void fraction and background flow) as seen in the similarities of the PDFs of standardized velocity for these diverse conditions, and (ii) our temporal sampling frequency of 50 Hz was adequate to capture essentially all of the most energetic liquid fluctuations.

To extend our analysis outside the plume core, PDFs of the velocity components at the plume edge $r=2b_{g}$ ( $b_{g}$ is the Gaussian plume radius; details in § 4) are shown in figure 7. At this large radial position, the flow is nearly single phase, except for the rare event of bubbles having drifted away from the core. The flow condition can be likened to bubbly flows with a small bubblance parameter  $b$ . We therefore choose the data from Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016) for homogeneous bubble swarms ( $\unicode[STIX]{x1D6FC}_{g}=0.2\,\%$ ) with the smallest $b=0.03$ and with bubble diameters in the range $d=2{-}4$  mm for comparison with our data. For  $\tilde{w}$ , the comparison is good. The PDF is much less positively skewed than in the core of our bubble plume and only deviates from the Gaussian PDF for probability densities smaller than $10^{-3}$ at the positive tail of the distribution. These deviations with low probabilities of occurrence are associated with the stray bubbles. For the horizontal velocities ( $\tilde{u} ,\tilde{v}$ ), their PDFs conform closely to the Gaussian PDF fit. For all velocity components, the stretched exponential tails observed at the plume centreline no longer exist at the plume edge. The liquid fluctuations are created by the plume entrainment and are modified by the rare passage of rising bubbles.

Figure 8.

Evolution of empirical PDF of raw velocities of Case A at $z/D=6.6$ and for $r=0{-}50~\text{mm}=0{-}0.92b_{g}$ ; $\tilde{w}$ (a), $\tilde{v}$ (b) and  $\tilde{u}$  (c).

Figure 9.

Evolution of empirical PDF of raw velocities of Case A at $z/D=6.6$ and for $r=55{-}105~\text{mm}=1{-}1.92b_{g}$ ; $\tilde{w}$ (a), $\tilde{v}$ (b) and  $\tilde{u}$  (c).

The results at the plume centreline (figure 6) and plume edge (figure 7) represent the two limiting cases in the continuous evolution of velocity PDFs as one travels radially outward from the centreline. To illustrate such evolution of the PDFs across the plume, we plot the data from Case A at $z/D=6.6$ for the range $r=0{-}105~\text{mm}=0{-}1.92b_{g}$ in figures 8 and 9. In these figures, all $\tilde{w}$ -PDFs are positively skewed, and a tertiary peak with positive velocity gradually emerges in the $\tilde{w}$ -PDFs as $r$ increases (see figure 8 a). This tertiary peak in the $\tilde{w}$ PDF decays away as $r$ continues to increase, disappearing at around $r=65$  mm. In contrast, all PDFs of the horizontal velocities ( $\tilde{u} ,\tilde{v}$ ) are symmetric across all  $r$ . These qualitative observations are also found for all other elevations $z/D$ and for Case B. In the next section, we decompose the $\tilde{w}$ -PDFs to better understand the origin of the tertiary peaks and of the negative velocities observed in figures 4, 8 and 9.

3.2 Decomposition of streamwise velocity $\tilde{w}$

To decompose the $\tilde{w}$ -PDFs, we first separate the extreme values of $\tilde{w}$ into two groups using a median filter:

(3.1) $$\begin{eqnarray}C_{2}\unicode[STIX]{x1D706}~\text{median}(|\tilde{w_{i}}-\text{median}(\tilde{w_{i}})|),\end{eqnarray}$$

where $C_{2}=1.35$ , $\unicode[STIX]{x1D706}=\sqrt{2\log (N)}$ and $N$ = number of data points in a time series (Goring & Nikora Reference Goring and Nikora2002; Cea et al. Reference Cea, Puertas and Pena2007). This set of constants is applied to all data.

Referring back to figure 5, the median filter identifies extreme values relative to the central primary peak which is assumed to take on a Gaussian shape. We then label the extreme values as follows: those with positive values of $\tilde{w}$ are termed the wake (tertiary peak) and those with negative $\tilde{w}$ termed the return flow (secondary peak). We also calculate the data fraction $n_{i}$ for each of these extreme value components in the full time series. The results for both Cases A and B are shown in figure 10. The radial distance $r$ (abscissa) has been normalized by the Gaussian plume radius $b_{g}$ (see § 4). For $r/b_{g}<1$ , the return flow accounts for less than ( $n_{1}=\;$ ) 2 % of the data; whereas, the wake contributes a significant portion at ( $n_{3}=\;$ ) 5–25 %. As one moves away from $r=0$ , $n_{3}$ decreases steadily to low values. Previous void fraction measurements for the same diffuser at the same gas flow rates (Seol et al. Reference Seol, Bhaumik, Bergmann and Socolofsky2007) found that the time-averaged $\unicode[STIX]{x1D6FC}_{g}$ is in the range 0.7–1.8 %. This correlates very well with our $n_{1}$ and indicates that large negative values of $\tilde{w}$ are related to the passage of an air bubble through the ADV sampling volume. Further, the wake fraction $n_{3}$ is much larger than that of the return flow because the latter is a very localized process. In their numerical simulations on bubble columns ( $\unicode[STIX]{x1D6FC}_{g}<4\,\%$ ), Riboux et al. (Reference Riboux, Legendre and Risso2013) show that bubble wakes have lengths ten times the bubble diameter. It is therefore reasonable to see a significant wake contribution in excess of the void fraction to the total data signal.

Figure 10.

Data fraction identified as return flow $n_{1}$ or wake $n_{3}$ by the median filter in (3.1).

Through this extreme value analysis, we have shown that the positive tails of the $\tilde{w}$ -PDFs correspond to a signature of the bubble wakes. In fact, there are three physical mechanisms creating the velocity fluctuations in  $\tilde{w}$ : (i) return fluid flow caused by rising bubbles, (ii) bulk entrained fluid in the bubble plume and (iii) bubble wakes. The boundary between contributions from the bulk flow and bubble wakes is not distinct because plume entrainment is a collective result of the bubble wakes. Nonetheless, this classification is useful as the bulk fluid and the bubble wake have different velocities, and we expect that the latter is larger and that the time-averaged velocity difference between these two contributions is close to the slip velocity between the two phases.

For the return flow contribution, we can show that the negative values are a necessity from mass conservation of the liquid fraction. Let us consider the return flow caused by the rise of a single bubble inside a cylinder of diameter $D$ in which the ambient fluid is initially motionless. If $d$ is the bubble diameter and $u_{rise}$ is the bubble rise velocity, then a void volume must be created at a rate equal to $(\unicode[STIX]{x03C0}/4)d^{2}u_{rise}$ along the trajectory of the bubble. As this void moves upward in an Eulerian perspective, the void at a position $\boldsymbol{x}$ immediately after passage of the bubble must be filled by a downflow of ambient water. The average velocity of this downflow is $\overline{u_{d}}=d^{2}u_{rise}/(D^{2}-d^{2})$ which shows that $\overline{u_{d}}$ can attain very large values as $D\rightarrow d$ . In figure 4, we interpret the large negative velocities as the return flow which have a mean of $\overline{u_{d}}=0.35~\text{m}~\text{s}^{-1}$ in this case. The rise velocity $u_{rise}$ is $0.50~\text{m}~\text{s}^{-1}$ which is the sum of the local bulk fluid velocity ( ${\sim}0.30~\text{m}~\text{s}^{-1}$ ) and slip velocity (Clift et al. Reference Clift, Grace and Weber1978, $0.24~\text{m}~\text{s}^{-1}$ ). Using $d=d_{50}=2.4$  mm and solving for  $D$ , we find $D=3.6$  mm. The periphery gap $D-d$ that allows drainage of the return flow is, thus, a 0.6 mm thick annular ring. This scale is of the same order as the 1 mm radial resolution of the ADV, and therefore, one can expect the return flow be captured in our data.

With these heuristic arguments at hand, we expect that $\tilde{w}$ is created by the interaction of all three mechanisms of return flow around bubbles, mean flow of the bubble plume and wake flow behind the bubbles. We hypothesize that each of these features has a Gaussian PDF $p_{i}(\tilde{w})$ for which their sum is equal to the total (raw) PDF $p_{total}(\tilde{w})$ . We perform a nonlinear curve fitting to calculate the mean and variance of each PDF, i.e. ( $\unicode[STIX]{x1D707}_{i},\unicode[STIX]{x1D70E}_{i}$ ). A constraint that must be satisfied is

(3.2) $$\begin{eqnarray}\int _{-\infty }^{\infty }p_{total}(\tilde{w})=\int _{-\infty }^{\infty }(p_{1}(\tilde{w})+p_{2}(\tilde{w})+p_{3}(\tilde{w}))\,\text{d}\tilde{w}=1.\end{eqnarray}$$

To satisfy (3.2), the individual PDF takes the following form

(3.3) $$\begin{eqnarray}p_{i}(\tilde{w})=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\frac{n_{i}}{\unicode[STIX]{x1D70E}_{i}}\text{e}^{-((\tilde{w}-\unicode[STIX]{x1D707}_{i})/\unicode[STIX]{x1D70E}_{i})^{2}},\end{eqnarray}$$

where $n_{i}$ is the data fraction corresponding to mechanism $i$ and $\sum _{i=1}^{3}n_{i}=1$ . $n_{i}$ is calculated from the experimental data by the median filter introduced at the beginning of this section (see also figure 10).

An example of this Gaussian decomposition is shown in figure 11. It can be seen that the total PDF is reasonably fitted by the sum of three Gaussians. Notably, the return flow is well separated from the rest of the flow and the upper end of the bulk entrained fluid flow is indistinguishable from the wake. To substantiate our heuristic arguments, the velocity differential $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{3}-\unicode[STIX]{x1D707}_{1}$ , which in principle should correspond to  $u_{s}$ , is plotted in figure 12 across the plume for both Cases A and B. The radial distance (abscissa) has been normalized by the Gaussian plume radius $b_{g}$ (see § 4). At all source elevations $z/D$ , a velocity plateau is observed for $r/b_{g}<0.8{-}0.9$ . For Case A, this plateau region gives $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=0.20~\text{m}~\text{s}^{-1}$ , and in Case B, the corresponding value is $0.20{-}0.25~\text{m}~\text{s}^{-1}$ . The velocity differential then drops off sharply and attains low values beyond $r/b_{g}=1$ . We note that value of $b_{g}$ in each case increases by 50 % across the measured range $z/D$ (see table 4), and the normalization $r/b_{g}$ effectively removes the plume growth and collapses the profiles. The values of the velocity plateaus are comparable to the terminal velocity ( $0.24~\text{m}~\text{s}^{-1}$ ) of a single bubble in a quiescent ambient. The observed variations are caused by group effects of bubbles and local state of fluid turbulence, both of which are known to sometimes increase and sometimes decrease the terminal velocity. From the figure, a possible definition of the bubble core radius is $b_{b}=0.9b_{g}$ , and this will be used to calculate the bubble column spreading ratio $\unicode[STIX]{x1D706}$ in § 4.

Figure 11.

Decomposition of the streamwise velocity $\tilde{w}$ into the return flow of bubbles, entrained fluid and bubble wake (Case A at $z/D=6.6$ , $r=0$ ).

Figure 12.

Velocity differential $\unicode[STIX]{x1D707}_{3}-\unicode[STIX]{x1D707}_{1}$ .

Figure 13.

Time-averaged slip velocity $u_{s}$ between the air- and liquid-phase streamwise velocities reported in the two experimental studies of Deen et al. (Reference Deen, Solberg and Hjertager2001) and Simiano et al. (Reference Simiano, Lakehal, Lance and Yadigaroglu2009). The gas inflow rate in the latter study was $Q_{o}=15~\text{l}~\text{min}^{-1}$ and the data were measured at an elevation 438 mm above the source. The inset shows the time-averaged streamwise velocity $W$ in normalized coordinates; the solid line represents a Gaussian curve fit to the data. The true slip velocity defined in Simiano et al. (Reference Simiano, Lakehal, Lance and Yadigaroglu2009) is $\overline{u_{s}}=\overline{(\tilde{w_{B}}-\tilde{w_{L}})\unicode[STIX]{x1D702}_{B}\unicode[STIX]{x1D702}_{L}}/\overline{\unicode[STIX]{x1D702}_{B}\unicode[STIX]{x1D702}_{L}}$ where $\unicode[STIX]{x1D702}_{B}$ and $\unicode[STIX]{x1D702}_{L}$ are phase-indicator functions that ensure the slip velocity is only calculated when both phases are present.

Finally, we plot in figure 13 the time-averaged slip velocity $u_{s}$ reported in Deen et al. (Reference Deen, Solberg and Hjertager2001) and in Simiano et al. (Reference Simiano, Lakehal, Lance and Yadigaroglu2009). These two experimental studies used phase-discriminating particle image velocimetry (PIV) to measure both air- and liquid-phase velocities. In the figure, we have plotted their data on a normalized $r/b_{g}$ -axis, where $b_{g}$ is determined from the radial profiles of the time-averaged streamwise velocity $W$ (shown in the inset of the figure). Comparing the values of our inferred $u_{s}$ shown in figure 12 to the directly measured $u_{s}$ in these two studies, we observed that the radial profiles have very similar magnitudes and shapes for $r<0.8b_{g}$ . All profiles show the tendency of $u_{s}$ to increase from the plume centreline towards the edge at $r=b_{g}$ . This is most apparent in the data of Simiano et al. (Reference Simiano, Lakehal, Lance and Yadigaroglu2009) and less so for Deen et al. (Reference Deen, Solberg and Hjertager2001) and for the present data. Beyond $r=b_{g}$ , $u_{s}$ attained a constant value in both previous studies. The sharp decrease in the values of $u_{s}$ observed in the present data can be understood by the fact that we did not directly measure the velocity of bubbles and the low void fraction in this region. As mentioned in the introduction (§ 1), the data in Simiano et al. (Reference Simiano, Lakehal, Lance and Yadigaroglu2009) were taken in the plume acceleration phase whilst the present data were in the classical plume regime. Their data have $Q_{o}=15~\text{l}~\text{min}^{-1}$ which is 10- and 30-fold the air inflow rates in Cases B and A, respectively. Although their data were measured at a height ( $=438$  mm) comparable to our first measurement station at 447 mm above source level, their plume was dynamically very different from ours. Their data have a value of $z/D$ of only 0.22 which is much smaller than the $z/D=2.2$ and 6.6 for Case B and Case A, respectively. It is then quite remarkable to see the $u_{s}$ values determined from both studies are so similar to each other, despite the plumes being so dynamically dissimilar. Overall, we conclude that our inferred slip velocities are consistent with direct measurements that were made in bubble plumes generated under very different experimental conditions.

3.3 Wavelet analysis on streamwise velocity $\tilde{w}$

To quantify the characteristic frequency at which rising bubbles inject kinetic energy into the surrounding fluid, a wavelet analysis is carried out on the time series of  $\tilde{w}$ . The wavelet method is widely used to analyse time series with non-stationary statistics (Daubechies Reference Daubechies1990). It is particularly relevant here since the injection of bubble-induced kinetic energy into the bulk plume flow occurs only intermittently in time; the void fraction is at most 2 % inside the plume. Conventional Fourier transform is not applicable for characterizing the rapid and short-lived velocity fluctuations caused by the bubbles because of its poor time localization property (Farge Reference Farge1992). We elect to use the Morlet wavelet in our analysis because of its strengths in decomposing turbulent velocity signals (Farge Reference Farge1992). We employ the continuous wavelet transform and follow the recommendations in Torrence & Compo (Reference Torrence and Compo1998) in choosing the wavelet parameters so that at least 98 % of the energy (variance) in the original time series is captured by the wavelet coefficients. Because the basis functions in the wavelet transform are scaled and shifted versions of a mother wavelet that is not periodic in itself, their scales do not directly correlate with an equivalent Fourier period – which we desire in order to quantify the frequency at which bubbles inject kinetic energy to the flow. It is possible to convert the wavelet scale into its equivalent Fourier period so that one can interpret the results in the usual frequency sense, and this has been reported in Torrence & Compo (Reference Torrence and Compo1998). For a Morlet wavelet, the scale is almost equal to the Fourier period.

Figure 14.

Continuous wavelet transform of streamwise velocity $\tilde{w}$ using the Morlet wavelet; number of scales analysed = 71, smallest scale $s_{o}=2dt=0.04$  s, largest scale $J=655$  s.

Figure 14 shows the wavelet results for Case A ( $z/D=6.6$ ) at two locations: (i)  $r=5~\text{mm}=0.1b_{g}$ (plume interior), and (ii)  $r=65~\text{mm}=1.2b_{g}$ (plume exterior). In the plume interior, there are many needle-like signatures with high frequencies in the time–frequency plot. These signatures are set against the backdrop of a band of energetic low frequencies that persists throughout the entire time record. By contrast, the plume exterior contains only the low-frequency band. From these difference, we can conclude that the high-frequency signatures in the plume core result from the bubble, which are nearly absent in the plum exterior.

A close-up view at $r=5$  mm for times between 230 and 250 s (shown in the same figure) reveals that the needle-like signatures are short-lived patches of high energy content. They span a range of frequencies from 1 to 20 Hz, but many of them have frequencies of approximately 10–15 Hz. The two most energetic patches are indicated by arrows, and these may represent the direct passage of bubbles through the ADV sampling volume, thus leaving behind strong signatures. The Nyquist criterion in our present experiment is $50~\text{Hz}/2=25$  Hz, which is approximately two times the frequency of these energetic signatures; hence, the bubble-induced turbulence observed in these wavelet analyses was well captured in our experiments.

In their model of bubble-induced turbulence based on large-scale wake interactions, Riboux et al. (Reference Riboux, Legendre and Risso2013) obtained the scaling relationship, $f_{cwi}=0.14u_{s}/d$ , for the frequency $f_{cwi}$ at which energy generated by the collective wake instability of a group of rising bubbles is injected into the surrounding fluid. Using present values of $u_{s}$ and  $d$ , $f_{cwi}=14$  Hz by this expression, and this value corresponds to that of the short-lived energetic patches in our wavelet analysis. Although the scaling relationship in Riboux et al. (Reference Riboux, Legendre and Risso2013) was developed from homogeneous bubble swarm simulations, the correspondence of this scaling relationship with our data again suggests that the nature of bubble-induced turbulence is insensitive to different flow conditions.

4.1 Mean flow

The instantaneous velocity vector $\tilde{\boldsymbol{u}}(\boldsymbol{x},t)$ is decomposed into a time-averaged part $\boldsymbol{U}(\boldsymbol{x})$ and a fluctuating part $\boldsymbol{u}(\boldsymbol{x},t)$ , e.g. in the streamwise direction $\tilde{w}=W+w$ . Figure 15 shows $W$ in normalized ordinates. A Gaussian profile was fitted to the data, and the centreline streamwise velocity $W_{c}$ and the Gaussian jet radius $b_{g}$ were obtained as a result of the curve fit. In both Cases A and B, the profile $W/W_{c}=\text{e}^{-(r/b_{g})^{2}}$ fits the data reasonably well; $b_{g}$ is defined as the radial location $r$ where $W=\text{e}^{-1}W_{c}$ . The discrepancies occurring beyond $r\geqslant 1.5b_{g}$ are probably caused by large-scale recirculation cells set up inside the tank by the time-steady plume (Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016).

Figure 15.

Normalized radial profiles of the time-averaged streamwise velocity  $W$ ; a comparison between data in the centre portion of the ADV’s measurement strip (a,b) and the data point at sweet spot (c,d); $Q_{o}=0.5~\text{l}~\text{min}^{-1}$ (a,c), $Q_{o}=1.5~\text{l}~\text{min}^{-1}$ (b,d).

Figure 16.

Comparisons of measured plume centreline velocity $W_{c}$ (a), Gaussian plume radius $b_{g}$ (b) and liquid-phase volume flux $Q$ (c) with published data.

Using dimensional analysis, the time-averaged plume centreline velocity $W_{c}(z)$ , Gaussian plume radius  $b_{g}(z)$ and fluid-phase volume flux $Q(z)$ have the following functional dependence (Bombardelli et al. Reference Bombardelli, Buscaglia, Rehmann, Rincon and Garcia2007),

(4.1a-c ) $$\begin{eqnarray}\frac{b_{g}(z)}{D}=f_{1}\left(\frac{z}{D}\right),\quad \frac{W_{c}(z)}{u_{s}}=f_{2}\left(\frac{z}{D}\right),\quad \frac{Q(z)}{u_{s}D^{2}}=f_{3}\left(\frac{z}{D}\right).\end{eqnarray}$$

Figure 16 shows comparisons of $W_{c}$ , $b_{g}$ and the fluid-phase volume flux $Q=\unicode[STIX]{x03C0}b_{g}^{2}W_{c}$ with previous experimental data in Seol et al. (Reference Seol, Bhaumik, Bergmann and Socolofsky2007) and datasets compiled in Milgram (Reference Milgram1983). The data compiled in the latter reference were obtained from bubble plumes with much larger scales than ours; the water depth was 50 m, 10 m and 3.66 m in Milgram (Reference Milgram1983), Fannelop & Sjoen (Reference Fannelop and Sjoen1980) and Milgram & Van Houten (Reference Milgram and Van Houten1982), respectively. We fit equations using the functional forms given in (4.1) to our data, and these are presented in the figure as the black, solid lines. Remarkable agreement can be seen between these empirical relationships derived from the present small-scale plumes and previous large-scale plumes. For  $W_{c}$ , the PIV data in Seol et al. (Reference Seol, Bhaumik, Bergmann and Socolofsky2007) are approximately 20–30 % lower than the empirical fit and other datasets. We expect this is because the wake velocities immediately behind the bubbles were not captured by their PIV technique, leading to underestimated liquid-phase velocities inside the plume core. Earlier studies like that of Milgram (Reference Milgram1983) employed mechanical current meters that had no phase-discrimination problem and so $W_{c}$ was probably not underestimated. If it is assume that $W$ at the plume edge was only slightly underestimated by Seol et al. (Reference Seol, Bhaumik, Bergmann and Socolofsky2007) due to the low void fraction there, a lower $W_{c}$ would probably lead to a larger  $b_{g}$ . Indeed, $b_{g}$ in Seol et al. (Reference Seol, Bhaumik, Bergmann and Socolofsky2007) is 18 % ( $0.132z$ ) larger than the value derived from the present data ( $0.114z$ ). This interpretation is further supported by  $Q$ ; despite the differences in $W_{c}$ and  $b_{g}$ , the data for $Q$ agree very well with the empirical fit. Since for dilute plumes, $Q=\unicode[STIX]{x03C0}b_{g}^{2}W_{c}$ , if $Q$ is to remain the same when $b_{g}$ is increased by 20 %, $W_{c}$ must decrease by 30 %. This reduction compares very well with their $W_{c}$ data. Finally, the steeper decay slope of $W_{c}$ found by a numerical and analytical approach in Bombardelli et al. (Reference Bombardelli, Buscaglia, Rehmann, Rincon and Garcia2007) (red solid line in figure 16) is not observed in the combined dataset that spans $0.01<z/D<11$ . Further experiments are needed in the range $10<z/D<100$ to explore the $(z/D)^{-1/2}$ decay suggested in Bombardelli et al. (Reference Bombardelli, Buscaglia, Rehmann, Rincon and Garcia2007).

Table 3.

Measured mean flow properties of bubble plumes in present study. C.I. is the 95 % confidence interval of the fitted $W_{c}$ and  $b_{g}$ . The $\sqrt{2}$ appearing in $\unicode[STIX]{x1D706}$ accounts for the difference between top-hat and Gaussian profile formulation of integral plume models.

Table 3 tabulates the numerical values of various mean flow properties. The bubble core radius $b_{b}$ is estimated from figure 12 in the manner described in § 3.2. $b_{b}$ is interpreted as a top-hat radius and, thus, a factor of $\sqrt{2}$ is multiplied to $b_{g}$ when calculating the spread width ratio  $\unicode[STIX]{x1D706}$ . A value of 0.64 is obtained for  $\unicode[STIX]{x1D706}$ , which is somewhat smaller than data in Seol et al. (Reference Seol, Bhaumik, Bergmann and Socolofsky2007). The liquid-phase spreading rate $\unicode[STIX]{x1D6FD}$ appears to decrease with $z/D$ and attains an asymptotic value of 0.11 for $z/D>8$ . Overall, our bubble plume data are consistent with other experimental datasets and exhibit the following dependencies on  $z$ ; $W_{c}\sim z^{-1/3}$ , $b_{g}\sim z$ and $Q\sim z^{5/3}$ . Integral properties of bubble plumes are adequately described by the independent parameters $u_{s}$ and  $D$ . As a further comparison, the entrainment coefficient $\unicode[STIX]{x1D6FC}$ in the asymptotic regime ( $z/D>5$ ) of a bubble plume can be computed from the integral equation of conservation of fluid volume flux, $\text{d}Q/\text{d}z=2\unicode[STIX]{x03C0}b_{g}^{2}W_{c}$ , using the three empirical fits. The coefficient $\unicode[STIX]{x1D6FC}$ can be written as,

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\frac{5}{3}\frac{1}{2\unicode[STIX]{x03C0}}\frac{C_{3}}{C_{2}C_{1}}=0.095,\end{eqnarray}$$

where $C_{1}=0.114$ , $C_{2}=1.6$ and $C_{3}=0.065$ are the empirical-fit constants for $b_{g}(z)$ , $W_{c}(z)$ and $Q(z)$ respectively. This value is somewhat larger than the generally accepted value of 0.083 for a single-phase plume (Fischer et al. Reference Fischer, List, Koh, Brooks and Imberger1979), but it is well within the range of published data. It should be noted that the spreading rate $\unicode[STIX]{x1D6FD}$ for a pure plume is equal to $(6/5)\unicode[STIX]{x1D6FC}_{p}$ (Lee & Chu Reference Lee and Chu2003, p. 95) by the equivalence of spreading hypothesis and entrainment hypothesis. Finally, the asymptotic plume densimetric Froude number $Fr_{p}$ can be computed from $Fr_{p}=(1/\sqrt{2})\sqrt{5/4\unicode[STIX]{x1D6FC}_{p}}\unicode[STIX]{x1D706}$ (Lee & Chu Reference Lee and Chu2003, p. 99); the factor of $1/\sqrt{2}$ is used to convert jet radius into diameter so that the definition of $Fr_{p}$ becomes identical to that in Wuest et al. (Reference Wuest, Brooks and Imboden1992) and Bombardelli et al. (Reference Bombardelli, Buscaglia, Rehmann, Rincon and Garcia2007). Using $\unicode[STIX]{x1D706}=0.64$ from our experiments, $Fr_{p}$ equals to 1.63, which is in agreement with the suggested value of 1.6 in Wuest et al. (Reference Wuest, Brooks and Imboden1992) and the analytical value of 1.7 in Bombardelli et al. (Reference Bombardelli, Buscaglia, Rehmann, Rincon and Garcia2007).

4.2 Turbulent Reynolds stresses

Before analysing the time-averaged turbulent stresses, we present the raw data for $\overline{w^{2}}$ in figure 17. Figure 17(a) shows the raw data at $z/D=8$ for Case A. Inside the plume core ( $r/b_{g}<1$ ), the data for $\overline{w^{2}}$ are discontinuous across overlapping strips of the ADV profile. This is witnessed by the upward trend along each ADV strip as the data approach the plume centreline (i.e. decreasing $r/b_{g}$ ). This inconsistency in the data was caused by the varying sensitivity along the ADV profile described in § 2.1: the ADV probe is more accurate at locations closer to the transmitter head, and these locations were located further away from the plume centre in the experiments. In figure 17(b), data points corresponding to only the sweet spot along each measurement strip are shown, and the turbulence intensity profiles appear continuous across the plume width. This behaviour is representative of all other Reynolds stress data in the present experiments. Hence, only the data obtained at the sweet spot along the ADV profile will be used for analysis of the turbulence quantities in this and later subsections.

Figure 17.

Raw ADV data on turbulent stresses; raw data (a) and data at sweet spot (b).

Figure 18.

Comparison of measured turbulent stresses by ADV and PIV data (Duncan et al. Reference Duncan, Seol and Socolofsky2009) for Case A.

Figure 19.

Comparison of measured turbulent stresses by ADV and PIV data (Duncan et al. Reference Duncan, Seol and Socolofsky2009) for Case B.

Figure 20.

Normalized radial profiles of turbulent stresses and turbulent kinetic energy $k=1/2(\overline{w^{2}}+\overline{u^{2}}+\overline{v^{2}})$ ; a double Gaussian profile (4.3) is fitted to the t.k.e. data.

Figure 18 shows the time-averaged turbulent stresses of Case A. All stresses (streamwise, radial and shear) are normalized using the local streamwise mean velocity at each measurement point $W$ and the Gaussian half-width  $b_{g}$ ; $W$ is used here to facilitate comparison with the PIV data (sampled at 125 Hz) in Duncan et al. (Reference Duncan, Seol and Socolofsky2009). The increasing trend from the plume centre towards the plume edge agrees with the PIV data. In particular, the comparisons on the radial stress $\overline{u^{2}}$ and shear stress $\overline{wu}$ are good. For the streamwise (vertical) stress  $\overline{w^{2}}$ , the ADV data are much higher than the PIV data. This difference is caused by the liquid velocities immediately behind bubble wakes that are not captured in the under-resolved PIV data. Two conclusions can be drawn from this comparison: (i) the 15 min sampling time in this study was adequate to obtain converged and repeatable second-order turbulence statistics and (ii) the plume turbulence is highly anisotropic inside the heterogeneous plume core. This strong anisotropy has also been reported in other experimental studies (Deen et al. Reference Deen, Solberg and Hjertager2001; Simiano et al. Reference Simiano, Zboray, de Cachard, Lakehal and Yadigaroglu2006, Reference Simiano, Lakehal, Lance and Yadigaroglu2009). The results for Case B are very similar to those in Case A, as seen in figure 19.

In figure 20 we plot the same data as in figures 18 and 19, but normalized by $W_{c}$ instead of the local  $W$ ; the out-of-plane stress $\overline{v^{2}}$ is also shown. Inside the bubble plume core, anisotropy of the three stresses is strong (figure 20 a,c). The normalized centreline values of $\sqrt{\overline{w^{2}}}$ , $\sqrt{\overline{v^{2}}}$ and $\sqrt{\overline{u^{2}}}$ are 0.38, 0.23 and 0.18, respectively. When compared to single-phase jets and plumes (Wang & Law Reference Wang and Law2002), the two horizontal intensities are comparable (within 10–15 % of their single-phase counterparts); whereas the streamwise (vertical) intensity in the two-phase plume is 1.5–2 times the single-phase value for $r/b_{g}<1$ . Both $\overline{w^{2}}$ and $\overline{v^{2}}$ show a peak around $r/b_{g}=0.55$ . Conversely, $\overline{u^{2}}$ remains level across the plume core. This feature of $\overline{u^{2}}$ is also observed in a high Reynolds number, single-phase jet ( $Re=1.4\times 10^{5}$ ) (Darisse, Lemay & Benaissa Reference Darisse, Lemay and Benaissa2015). Beyond $r/b_{g}=1.5$ , the liquid turbulence attains isotropy as the flow at the plume edge is due to plume entrainment. As discussed for figure 7, the flow is single phase at $r/b_{g}>1.5$ (except for the rare passage of stray air bubbles), and the PDFs of liquid fluctuations are Gaussian. The off-centre peak in the $\overline{v^{2}}$ -profile and the observation that $\overline{v^{2}}>\overline{u^{2}}$ warrant more discussion here. In Darisse et al. (Reference Darisse, Lemay and Benaissa2015), the authors explained for a single-phase jet in air that the differences between the two horizontal stresses are caused by velocity events with a large incident angle relative to the jet centreplane, where $\overline{u^{2}}$ resides – hence, resulting from the three-dimensionality of turbulence. As a result, the common approximation $\overline{v^{2}}\approx \overline{u^{2}}$ is not true, but is a good approximation in single-phase jets, where their ratio is approximately 1.05. In bubble plumes, the three-dimensionality of the flow is expected to be stronger as the millimetre-sized bubbles have strong wakes and exhibit wobbling trajectories as they rise, leading to a larger $\overline{v^{2}}/\overline{u^{2}}$ ratio. Within the bubbly plume core ( $r<b_{g}$ ), this ratio is between 1.1 and 1.3 compared to the value 1.05 reported in single-phase jets.

The profiles of normalized turbulent kinetic energy $k=(1/2)(\overline{w^{2}}+\overline{u^{2}}+\overline{v^{2}})/W_{c}^{2}$ shown in figure 20 have a peak at $r/b_{g}=0.55$ and can be fitted by a double Gaussian profile:

(4.3) $$\begin{eqnarray}k(r)=0.18W_{c}^{2}(z)\exp \left(-\frac{(r-0.55b_{g})^{2}}{(0.5b_{g})^{2}}\right)+\exp \left(-\frac{(r+0.55b_{g})^{2}}{(0.5b_{g})^{2}}\right).\end{eqnarray}$$

We also use a double Gaussian profile to curve fit the individual components of the turbulent stresses, and the resulting curve fit equations are given in appendix D.

For comparison, we include in figure 20 the predictions (dashed lines) from Reynolds-averaged Navier–Stokes (RANS) simulations of the plume that neglect the production of t.k.e. by bubbles. The computational details are given in Lai (Reference Lai2015), and we briefly describe the settings below. The simulations were done with the commercial computational fluid dynamics (CFD) software FLUENT (Ansys Inc.). A mixture model was used in which a single set of governing equations was solved for the mass-weighted mixture properties. The continuous liquid phase and the dispersed bubble phase were modelled as two interpenetrating Eulerian fluids. The air volume fraction was considered as a conservative scalar and tracked by a transport equation. The momentum coupling between the two phases was modelled using the Schiller and Naumann dynamic drag law. The standard $k{-}\overline{\unicode[STIX]{x1D716}}$ model was used for closure of the liquid phase turbulence. The $k$ and $\overline{\unicode[STIX]{x1D716}}$ equations were not augmented by additional source or sink terms to account for bubble-induced liquid turbulence (Vaidheeswaran & Hibiki Reference Vaidheeswaran and Hibiki2017). The simulations therefore only include the production by fluid mean shear (Pope Reference Pope2000). We simulated only half of the time-steady plume by virtue of symmetry in a computational domain that was half of our experimental water tank, i.e. 1 m high by 1 m across the plume by 0.5 m binormal to the vertical and lateral coordinates. The simulations were initialized using the gas inflow rates in table 2, and the gas was released into the domain through a 10 mm inlet (velocity inflow boundary) at the bottom face of the domain. For simplicity, we considered monodisperse bubbles with a median diameter $d_{50}=2.4$  mm, matching that of our experiments. The $k{-}\overline{\unicode[STIX]{x1D716}}$ model with its standard coefficients do not yield satisfactory predictions of the mean flow quantities $W_{c},\unicode[STIX]{x1D6FD}$ and  $\unicode[STIX]{x1D706}$ . Hence, we adjusted the eddy viscosity coefficient $C_{\unicode[STIX]{x1D707}}$ and the turbulent Schmidt number $\unicode[STIX]{x1D70E}$ from their standard values, 0.09 and 0.75, to 0.15 and 0.1, respectively, in order to match the experimental results in Seol et al. (Reference Seol, Bhaumik, Bergmann and Socolofsky2007). All other model coefficients retained their standard values. Comparing the predicted and measured  $k$ , both exhibit an off-centre peak at $r/b_{g}=0.5$ . The predictions are, however, only 50–70 % of the measured values inside the plume core. Beyond the core, $r/b_{g}>1$ , the two become comparable in magnitude. The comparison, thus, reveals that: (i) while calibration of the closure model coefficients can predict the mean flow, the predictions are correct for the wrong reason – the model neglects the bubble-induced turbulence; and (ii) the t.k.e. injected by bubbles into the liquid is non-negligible.

An example where the t.k.e. production by bubbles is partly accounted for is in the momentum amplification factor $\unicode[STIX]{x1D6FE}$ used in integral plume models (Socolofsky et al. Reference Socolofsky, Bhaumik and Seol2008). We evaluate $\unicode[STIX]{x1D6FE}$ using our measured turbulent stresses; $\unicode[STIX]{x1D6FE}$ accounts for the kinematic momentum flux carried by the plume turbulence in the model formulation, and it is defined by $\unicode[STIX]{x1D6FE}=(\overline{M}+M_{t})/\overline{M}$ , where $\overline{M}=2\unicode[STIX]{x03C0}\int _{0}^{\infty }rW^{2}\,\text{d}r$ is the mean momentum flux and $M_{t}=2\unicode[STIX]{x03C0}\int _{0}^{\infty }r(\overline{w^{2}}-(1/2)(\overline{u^{2}}+\overline{v^{2}}))\,\text{d}r$ is the momentum flux carried by turbulence. These expressions follow from the boundary-layer approximation to the governing equations of the plume (Lee & Chu Reference Lee and Chu2003). For single-phase jets and plumes $\unicode[STIX]{x1D6FE}$ is approximately 1.1 (Wang & Law Reference Wang and Law2002). For bubble plumes, Milgram (Reference Milgram1983) reported values in the range of 1–3. Using the double-Gaussian curve fits for the three normal stresses reported in appendix D, $M_{t}=C[0.22-(1/2)(0.097+0.035)]b_{g}^{2}W_{c}^{2}$ where $C=2\unicode[STIX]{x03C0}\int _{0}^{\infty }(r/b_{g})G(r/b_{g})d(r/b_{g})=6.1262$ is a shape factor. Together with $\overline{M}=(\unicode[STIX]{x03C0}/2)b_{g}^{2}W_{c}^{2}$ , we obtain $\unicode[STIX]{x1D6FE}=(0.943+\unicode[STIX]{x03C0}/2)/(\unicode[STIX]{x03C0}/2)=1.6$ in Case A at $z/D=8$ . If the variability of the data points used in the curve fit of $k$ is accounted for (see figure 20), $\unicode[STIX]{x1D6FE}$ is in the range 1.4–1.6. A similar analysis done for Case B yields $\unicode[STIX]{x1D6FE}$ in the range 1.45–1.55. The fact that two-phase plumes have larger $\unicode[STIX]{x1D6FE}$ values reflects their stronger turbulence anisotropy compared to single-phase jets/plumes.

Finally, some insight on the magnitude of the t.k.e. generated by the bubbles can be extracted from a recent simulation of Case A using a large-eddy simulation (LES) code in Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016). An Eulerian–Lagrangian approach was used to predict the bubble plume dynamics. The Smagorinsky subgrid-scale (SGS) model with a constant coefficient of 0.1 was used to model the unresolved fluid turbulence. The LES code has been previously verified to produce accurate results in complex geometries, including open-channel flows over a rough bed, compound channel flow, vegetated flows and transport and deposition of fine sediments in open channels (Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016). The model performs very well for prediction of $W_{c}$ , $\sqrt{\overline{u^{2}}}$ and $\overline{wu}$ ; however, the model under-predicts $\sqrt{\overline{w^{2}}}$ for $r/b_{g}<1$ by up to a factor of two. This difference is because the thin boundary layers around the bubbles were not resolved in the simulation. This motivates us to quantify the t.k.e. production by bubbles.

4.3 Conservation of the total plume momentum flux

Using the mean streamwise velocity $W$ and the turbulent Reynolds stresses  $R_{ii}$ , we assess in figure 21 the conservation of the total plume kinematic momentum flux  $M$ . $M$ has two contributions and is equal to $\overline{M}+M_{t}$ : $\overline{M}=2\unicode[STIX]{x03C0}\int _{0}^{\infty }rW^{2}\,\text{d}r$ is the mean momentum flux carried by the mean flow; $M_{t}=2\unicode[STIX]{x03C0}\int _{0}^{\infty }r(\overline{w^{2}}-(1/2)(\overline{u^{2}}+\overline{v^{2}}))\,\text{d}r$ is the turbulent momentum flux carried by the turbulence (Wang & Law Reference Wang and Law2002; Lee & Chu Reference Lee and Chu2003). The analytical prediction of $M$ as a function of $z$ for a pure vertical plume is given by (Lee & Chu Reference Lee and Chu2003, p. 83),

(4.4) $$\begin{eqnarray}M=(3\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FD}_{tophat}/4)^{2/3}F_{o}^{2/3}z^{4/3}\approx 0.36F_{o}^{2/3}z^{4/3},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{tophat}=\sqrt{2}\unicode[STIX]{x1D6FD}$ is the top-hat spreading rate, and $F_{o}=Q_{o}g$ is the kinematic buoyancy flux exerted by the air bubbles on the surrounding water. Previous experimental studies on single-phase plumes have found the proportionality constant to be 0.29 (Papanicolaou & List Reference Papanicolaou and List1988), 0.34 (Shabbir & George Reference Shabbir and George1994) and 0.38 (Wang & Law Reference Wang and Law2002).

Using the values of $Q_{o}$ in table 2 and the average (over the four $z/D$ elevations in each case) $\unicode[STIX]{x1D6FD}=0.114$ and 0.117 for Cases A and B, the predicted $M$ for Cases A and B is plotted as a solid line and a dashed line in figure 21(a), respectively. In physical units, the agreement between the analytical predictions and measurements is very good; the deviations are within $\pm 10\,\%$ . The measured data exhibit the expected $4/3$ power law with  $z$ . In figure 21(b), we plot the non-dimensionalized total momentum flux $M/(u_{s}D)^{2}$ as a function of  $z/D$ . The non-dimensionalization effectively collapses the data from both cases onto a single curve. This reinforces the choice of $u_{s}$ and $D$ as the two governing parameters for dilute bubble plumes. We conclude that our bubble plumes obey the scaling law for momentum conservation in single-phase plumes ( $M\sim z^{4/3}$ ) and that our experimental data are in agreement with the law of conservation of momentum over the range of heights measured here ( $2<z/D<12$ ).

Figure 21.

Comparison of measured plume total kinematic momentum flux $M=\overline{M}+M_{t}$ and the analytical prediction $M=(3\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FD}_{tophat}/4)^{2/3}F_{o}^{2/3}z^{4/3}$ (Lee & Chu Reference Lee and Chu2003):  $\unicode[STIX]{x1D6FD}_{tophat}=\sqrt{2}\unicode[STIX]{x1D6FD}$ is the top-hat spreading rate; $F_{o}=Q_{o}g$ is the kinematic buoyancy flux exerted by the air bubbles to the liquid. (a) In physical space. (b) In non-dimensional space; $C_{4}=(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}^{2})^{2/3}(3\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FD}_{tophat}/4)^{2/3}=\unicode[STIX]{x03C0}(3\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}_{tophat})^{2/3}$ .

4.4 Triple velocity correlation – turbulent transport of components of t.k.e.

The normalized triple velocity correlations $\overline{wu_{i}^{2}}/W_{c}^{3}$ and $\overline{uu_{i}^{2}}/W_{c}^{3}$ are shown in figures 22 and 23. In these figures, a correction for Doppler noise is possible in theory, but here this correction was not made because the skewness factor $\overline{\unicode[STIX]{x1D70E}_{D}^{3}}$ could not be estimated (see appendix A). Except for the component $\overline{uv^{2}}$ (figure 23 c,f), all other components exhibit a reasonable collapse of the data to a single curve for each term when they are normalized by $W_{c}^{3}$ and  $b_{g}$ . For the three axial components $\overline{wu_{i}^{2}}/W_{c}^{3}$ , the profile shapes and magnitudes are dissimilar between Cases A and B (see figure 22). This may be explained by their different void fractions: Case B has three times the $Q_{o}$ of Case A. On the contrary, the three radial components have similar profiles in both cases (see figure 23).

Figure 22.

Normalized radial profiles of axial turbulent $\overline{wu_{i}^{2}}$ transport of t.k.e. Different symbols represent data measured at different elevations above source; 44.7 cm, ○ –; 54.7 cm, ▵ –; 64.7 cm, ♢ – and 74.7 cm, ▫ –. The black solid lines are curve fits to the data.

Figure 23.

Normalized radial profiles of radial turbulent $\overline{uu_{i}^{2}}$ transport of t.k.e. Different symbols represent data measured at different elevations above source; 44.7 cm, ○ –; 54.7 cm, ▵ –; 64.7 cm, ♢ – and 74.7 cm, ▫–. The black solid lines are curve fits to the data. The magenta solid lines are data for a high Reynolds number, single-phase pure jet reported in Darisse et al. (Reference Darisse, Lemay and Benaissa2015). Note that the $\overline{uw^{2}}$ jet data have been multiplied by 5 to facilitate the comparison with the corresponding data of the bubble plume.

For comparison, the data of a high Reynolds number, single-phase pure jet reported in Darisse et al. (Reference Darisse, Lemay and Benaissa2015) are included in the plots of figure 23 (solid magenta lines) to highlight this similarity. The $\overline{uw^{2}}/W_{c}^{3}$ from the jet data have been multiplied by five to facilitate a comparison of the profile shape. These data comparisons show that $\overline{uw^{2}}/W_{c}^{3}$ is negative close to the centreline axis of a jet or bubble plume, becoming positive outside the core, and $\overline{u^{3}}/W_{c}^{3}$ is positive across the whole jet or bubble plume. We note that the signs will change on the $-r/b_{g}$ axis because of axisymmetry (Darisse et al. Reference Darisse, Lemay and Benaissa2015). The component $\overline{uv^{2}}$ fails to converge in the present experiments inside the bubble core $r/b_{g}<1$ . This may be caused by the weak correlation between the plume entrainment flow $u$ and the out-of-plane oscillatory liquid motions $v$ caused by the bubbles; apparently a much longer sampling time is required to produce converged statistics. Outside the core, $\overline{uv^{2}}$ is similar to  $\overline{uu^{2}}$ , suggesting a similar radial transport. Overall, the results indicate that transport of t.k.e. by the entrainment velocity fluctuations $u$ is not sensitive to void fraction; whereas, transport by $w$ is sensitive to gas flow rate and, hence, void fraction.

For the purpose of computing the t.k.e. budget, we fit curves to the radial profiles of the turbulent transport components $\overline{wu_{i}^{2}}/W_{c}^{3}$ and $\overline{uu_{i}^{2}}/W_{c}^{3}$ in the form

(4.5) $$\begin{eqnarray}P_{k}(r/b_{g})\exp (-C(r/b_{g})^{2}),\end{eqnarray}$$

where $P_{k}(r/b_{g})$ is a polynomial of degree $k$ , and $C$ is a constant. This choice follows the experimental study on single-phase jets in Hussein, Capp & George (Reference Hussein, Capp and George1994), and the full expressions of each curve fit are given in appendix D. Because the $\overline{uv^{2}}$ data fail to converge, we approximate $\overline{uv^{2}}$ by  $\overline{u^{3}}$ . All profiles are practically zero beyond  $1.2b_{g}$ , indicating that the turbulent transport of t.k.e. is confined within the plume core. This is quite different from a single-phase jet where the profiles extend to $r=2b_{g}$ (Darisse et al. Reference Darisse, Lemay and Benaissa2015).

4.5 Turbulent dissipation rate estimates

We calculate the mean liquid dissipation via the pseudo-dissipation expression $\overline{\unicode[STIX]{x1D716}}=\unicode[STIX]{x1D708}\overline{(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})}$ , where $\unicode[STIX]{x1D708}=0.935\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ is the kinematic viscosity of water at $23\,^{\circ }\text{C}$ (the temperature in the experiments). The difference between pseudo-dissipation and the true dissipation rate (stress power) is small in high Reynolds number turbulent flows (Pope Reference Pope2000). The full expression of $\overline{\unicode[STIX]{x1D716}}$ is:

(4.6) $$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D716}}_{full} & = & \displaystyle \unicode[STIX]{x1D708}\left[\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{2}}{\unicode[STIX]{x2202}x_{2}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{3}}{\unicode[STIX]{x2202}x_{3}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{2}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{3}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{2}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\overline{\left(\frac{\unicode[STIX]{x2202}u_{2}}{\unicode[STIX]{x2202}x_{3}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{3}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}+\overline{\left(\frac{\unicode[STIX]{x2202}u_{3}}{\unicode[STIX]{x2202}x_{2}}\right)^{2}}\right],\end{eqnarray}$$

where $(u_{1},u_{2},u_{3})=(w,u,v)$ and $(x_{1},x_{2},x_{3})=(z,x~\text{or}~r,y)$ . We can directly calculate the mean-square radial gradients $\overline{(\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{2})^{2}}$ , $\overline{(\unicode[STIX]{x2202}u_{2}/\unicode[STIX]{x2202}x_{2})^{2}}$ and $\overline{(\unicode[STIX]{x2202}u_{3}/\unicode[STIX]{x2202}x_{2})^{2}}$ from the data. For the streamwise mean-square gradients $\overline{(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{1})^{2}}$ , we calculate them indirectly from the temporal data via the modified Taylor hypothesis (Heskestad Reference Heskestad1965):

(4.7) $$\begin{eqnarray}\displaystyle \overline{\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}=\frac{1}{W^{2}}\overline{\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}t}\right)^{2}}K^{-1}. & & \displaystyle\end{eqnarray}$$

The scalar coefficient $K\geqslant 1$ is dependent on the local turbulence intensities and the ratios of the mean-square gradients. The original Taylor hypothesis is recovered when $K=1$ . Finally, the unavailable out-of-plane mean-square gradients $\overline{(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{3})^{2}}$ are replaced by $\overline{(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{2})^{2}}$ .

For comparison, we also evaluate $\overline{\unicode[STIX]{x1D716}}$ with the local isotropy assumption (Taylor Reference Taylor1935) and with the local axisymmetry assumption (George & Hussein Reference George and Hussein1991):

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D716}}_{iso}=15\unicode[STIX]{x1D708}\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}=7.5\unicode[STIX]{x1D708}\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{2}}\right)^{2}} & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D716}}_{axi}=\unicode[STIX]{x1D708}\left[\frac{5}{3}\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}+2\overline{\left(\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{2}}\right)^{2}}+2\overline{\left(\frac{\unicode[STIX]{x2202}u_{2}}{\unicode[STIX]{x2202}x_{1}}\right)^{2}}+\frac{8}{3}\overline{\left(\frac{\unicode[STIX]{x2202}u_{3}}{\unicode[STIX]{x2202}x_{2}}\right)^{2}}\right]. & \displaystyle\end{eqnarray}$$

In both (4.8) and (4.9), all terms can be computed from the data either directly, for $\overline{(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{2})^{2}}$ , or indirectly via Taylor’s hypothesis, for $\overline{(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{1})^{2}}$ so that there is no need to approximate unavailable out-of-plane gradients. For $\overline{\unicode[STIX]{x1D716}}_{iso}$ , we also use $\overline{(\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{2})^{2}}$ in the calculation because it can be directly computed from the data. The required gradients are computed by finite difference and the details are explained in appendix B. The effects of finite spatial resolution on the computed gradients are assessed in appendix C.

Figure 24.

Time-averaged normalized mean-square gradients for Case A. Different symbols represent data measured at different elevations above source; 44.7 cm, ○ –; 54.7 cm, ▵ –; 64.7 cm, ♢ – and 74.7 cm, ▫ –.

Figure 24 plots the radial profiles of the mean-square gradients (normalized by $(W_{c}/b_{g})^{2}$ ) for Case A. Inside the bubbly core ( $|r/b_{g}|\leqslant 1$ ), the radial mean-square gradients, except for $\overline{(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}r)^{2}}$ , are an order of magnitude larger than the streamwise mean-square gradients, and the ratios of the gradients satisfy neither local isotropy nor local axisymmetry assumptions, e.g. $2\overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z)^{2}}\neq \overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}r)^{2}}$ and $\overline{(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)^{2}}\neq \overline{(\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}x)^{2}}$ . Of the six measurable gradients, $\overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}r)^{2}}$ is the largest. This can be interpreted as connected to the highly dissipative, thin boundary layer around the bubbles (Santarelli et al. Reference Santarelli, Roussel and Frohlich2016): over a radial distance of the order of the bubble diameter or less, the instantaneous streamwise velocity $\tilde{w}$ has to transit from a value close to the bubble rise velocity to a value corresponding to the much lower background plume flow. The streamwise counterpart to $\overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}r)^{2}}$ , given by $\overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z)^{2}}$ , is small and has only one tenth the magnitude of the radial term. This is expected because the wake region behind each bubble is much longer in length than it is wide (Risso Reference Risso2018). Outside the bubbly core, the mean-square gradients tend to a state of local isotropy, e.g. $2\overline{(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}r)^{2}}\sim \overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}r)^{2}}\approx \overline{(\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}r)^{2}}\approx \overline{(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z)^{2}}$ (Taylor Reference Taylor1935). Together with the t.k.e. (figure 20), this result demonstrates that the boundary-layer flow in the single-phase region of the bubble plume is largely isotropic. The results for Case B (not shown) are similar.

Figure 25.

Time-averaged dissipation rate $\overline{\unicode[STIX]{x1D716}}$ . Different symbols represent data measured at different elevations above source; 44.7 cm, ○ –; 54.7 cm, ▵ –; 64.7 cm, ♢ – and 74.7 cm, ▫ –.

Figure 25 depicts the normalized dissipation rate profiles for both cases. From these data we find that $\overline{\unicode[STIX]{x1D716}}_{axi}\approx \overline{\unicode[STIX]{x1D716}}_{full}$ , that $\overline{\unicode[STIX]{x1D716}}_{iso}=15\unicode[STIX]{x1D708}\overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z)^{2}}$ under-predicts the dissipation inside the bubbly core, while $\overline{\unicode[STIX]{x1D716}}_{iso}=7.5\unicode[STIX]{x1D708}\overline{(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}r)^{2}}$ over-predicts it, and that all expressions for the dissipation produce similar results outside the bubbly core. It is fortuitous that $\overline{\unicode[STIX]{x1D716}}_{axi}\approx \overline{\unicode[STIX]{x1D716}}_{full}$ as we have shown in figure 24 that the mean-square gradients do not satisfy the local axisymmetry relations (George & Hussein Reference George and Hussein1991). Moreover, the behaviour inside the bubble core once again highlights the highly anisotropic nature of the flow in this region.

4.5.1 Conservation of total kinetic energy

Dissipation rate estimates derived from physical experiments commonly suffer from errors due to insufficient spatial resolution and instrument noise; we have implemented the correction method in Hurther & Lemmin (Reference Hurther and Lemmin2001) to reduce the latter’s influence on the results (see appendix A). However, these two sources of error are difficult to disentangle in real experiments in which the true value of $\overline{\unicode[STIX]{x1D716}}$ is unknown and in which the adjustment of experimental parameters can reduce one source of error but augment the other.

Here, we assess the accuracy of our estimates of $\overline{\unicode[STIX]{x1D716}}$ by carrying out a total kinetic energy balance of the bulk plume flow. This balance equation is obtained from the sum of the transport equations for the mean kinetic energy $K=(1/2)(W^{2}+U^{2}+V^{2})\approx (1/2)W^{2}$ (where we have used $W\gg U,V$ in jets and bubble plumes) and the t.k.e. $k=(1/2)(\overline{w^{2}}+\overline{u^{2}}+\overline{v^{2}})$ (Pope Reference Pope2000, pp. 125–126). In Lai (Reference Lai2015), the equation for $K+k$ has been further developed for axisymmetric jets and plumes under the integral framework (Lee & Chu Reference Lee and Chu2003). In essence, the computation of the total kinetic energy balance is to evaluate

(4.10) $$\begin{eqnarray}\frac{\text{D}\langle K+k\rangle }{\text{D}t}=-\langle \overline{\unicode[STIX]{x1D716}}\rangle +\langle Wg^{\prime }\rangle ,\end{eqnarray}$$

where $\langle \cdot \rangle =2\unicode[STIX]{x03C0}\int _{0}^{\infty }r\,\text{d}r$ denotes a volume average, and $g^{\prime }=(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{ref})g\approx g=9.81~\text{m}~\text{s}^{-2}$ is the reduced gravity for air bubbles in water.

By the Gauss divergence theorem, the volume averaging removes the gradient flux terms that only redistribute kinetic energy in space by pressure–velocity correlations and by viscous diffusion. Because we only measured at discrete plume elevations, we integrate the volume-averaged total kinetic energy equation in the streamwise direction  $z$ . This gives the following equation for steady bubble plumes:

(4.11) $$\begin{eqnarray}\overbrace{\left({\textstyle \frac{1}{2}}I_{1}+C_{I_{2}}I_{2}\right)[b_{g}^{2}(z)W_{c}^{3}(z)]_{z_{1}}^{z_{2}}}^{\int (\text{D}\langle K+k\rangle /\text{D}t)\,\text{d}z\approx \int _{z_{1}}^{z_{2}}(\unicode[STIX]{x2202}\langle W(K+k)\rangle /\unicode[STIX]{x2202}z)\,\text{d}z}=\overbrace{-{\textstyle \frac{1}{2}}[\langle \overline{\unicode[STIX]{x1D716}}(z_{1})\rangle +\langle \overline{\unicode[STIX]{x1D716}}(z_{2})\rangle ]\unicode[STIX]{x0394}z}^{-\int _{z_{1}}^{z_{2}}\langle \overline{\unicode[STIX]{x1D716}}\rangle \,\text{d}z}+\overbrace{1.36F_{o}\unicode[STIX]{x0394}z}^{\int _{z_{1}}^{z_{2}}\langle wg\rangle \,\text{d}z}.\end{eqnarray}$$

In (4.11), $z_{1}$ and $z_{2}$ ( $z_{2}>z_{1}$ ) are elevations above the air source, and $\langle \overline{\unicode[STIX]{x1D716}}(z)\rangle =b_{g}(z)W_{c}^{3}(z)(C_{I_{3}}I_{3})$ ; we have approximated the integral $-\int _{z_{1}}^{z_{2}}\langle \overline{\unicode[STIX]{x1D716}}\rangle \,\text{d}z$ by the average of the dissipation measured at $z_{1}$ and  $z_{2}$ : $F_{o}=Q_{o}g$ is the kinematic buoyancy flux due to air injection at the source. The constant factor 1.36 accounts for a narrower bubble core than the streamwise fluid velocity profile. Equation (4.11) is a statement of total kinetic energy balance of a conical control volume bounded by any two elevations. For the results presented below, we consider control volumes formed by any two consecutive elevations, leading to $\unicode[STIX]{x0394}z=10$  cm. For steady homogeneous bubble swarms rising in a quiescent ambient, the advection term (left-hand side of the equation) is zero and the usual estimate $\overline{\unicode[STIX]{x1D716}}=\unicode[STIX]{x1D6FC}_{g}V_{r}g$ is recovered (Lance & Bataille Reference Lance and Bataille1991); hence, for that case, a volume averaging is not needed as the flow is homogeneous.

The shape constants $(I_{1},I_{2},I_{3})$ depend on the form of curve fits applied to  $W^{3}$ , $Wk$ and $\overline{\unicode[STIX]{x1D716}}$ . These are:

1. (i) the mean kinetic energy flux $2\unicode[STIX]{x03C0}\int rW^{3}\,\text{d}r\sim I_{1}=\int _{0}^{\infty }2\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}\text{e}^{-3\unicode[STIX]{x1D702}^{2}}\,\text{d}\unicode[STIX]{x1D702}$ ;

2. (ii) the turbulent kinetic energy flux $2\unicode[STIX]{x03C0}\int rWk\,\text{d}r\sim I_{2}=\int _{0}^{\infty }2\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}\text{e}^{-\unicode[STIX]{x1D702}^{2}}G(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}$ ;

3. (iii) the dissipation $2\unicode[STIX]{x03C0}\int r\overline{\unicode[STIX]{x1D716}}\,\text{d}r\sim I_{3}=\int _{0}^{\infty }2\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}F_{\overline{\unicode[STIX]{x1D716}}}(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}$ .

The shape profiles $G(\unicode[STIX]{x1D702})$ and $F_{\overline{\unicode[STIX]{x1D716}}}(\unicode[STIX]{x1D702})$ are double-Gaussian fits applied to radial profiles of $k$ and $\overline{\unicode[STIX]{x1D716}}$ , respectively (see appendix D). The shape constants $C_{I_{2}},C_{I_{3}}$ are related to fitted profiles of $k$ and $\overline{\unicode[STIX]{x1D716}}$ , respectively. Table 4 tabulates the constants found in each case. The degree to which (4.11) is satisfied is expressed as the ratio right-hand side/left-hand side of the equation (last column in the table). Because the plume is accelerated by buoyancy  $F_{o}$ , the right-hand side is positive. We found the that balance is satisfied to within $\pm 30\,\%$ in both cases with the raw dissipation rate estimates. It is possible to correct for the observed imbalance by enforcing equality in (4.11) by premultiplying the dissipation profile by a constant such that the balance ratio is equal to one. We have chosen not to do this so as to present and interpret the objective experimental results.

Table 4.

Balance of the total kinetic energy equation. The balance is defined as the ratio of the right-hand side to the left-hand side of (4.11). At all elevations, $(I_{1},I_{2},C_{I_{2}})=(1.0472,1.8128,0.18)$ in both cases.