Oxygen transfer from bubble-plumes

Rights This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing.The following article appeared in Physics of Fluids 30, 107104 (2018) and may be found at https://doi.org/10.1063/1.5040819.Copyright 2018 Author(s). This article is distributed under a Creative Commons Attribution (CC BY) License. Rights(URL) https://creativecommons.org/licenses/by/4.0/ Type article


I. INTRODUCTION
Ocean wave breaking produces large numbers of air bubbles in turbulent flow that form bubble plumes beneath the sea surface, 1-4 contributing to air-sea gas exchanges. [5][6][7] Deeper entrainment of bubbles into the plume results in gas dissolution from the bubbles in inner bulk sea water and enhances gas diffusion due to turbulent mixing. Wallace and Wirick 8 found that the saturation of dissolved oxygen (DO) in sea surface layers depends on the ocean wave height, and suggested explicit effects of breaking-wave-induced aeration on gas concentrations in the open ocean. Farmer et al. 5 observed deep penetration of bubble plumes under strong wind conditions, suggesting possible gas dissolution from bubbles at deeper levels during storm events. In this context, a diffusive air-sea gas transport process across the sea surface, which is commonly assumed in predictions of equilibrium air-sea gas exchange, 9 may not be appropriate for predictions of net gas transfer in the ocean surface layer. In shallow water, oxygen transport from bubbles may be an important factor in assessing the coastal environment because coastal habitats and marine ecosystems are maintained by the steady oxygen supply from aerated surf zone flows. 10 a) Electronic mail: niida@criepi.denken.or.jp b) Electronic mail: yasunori@eng.hokudai.ac.jp The local dynamics of bubbles and turbulence influence the oxygen concentration in bulk water flows through their interactions. Turbulent shear deforms the bubbles to induce their breakup into smaller bubbles. 11 The total interface area of the bubbles increases through this breakup process, resulting in larger gas transfer across bubble interfaces. Bubble motion relative to the flow is another important factor defining gas concentration in the bubble flows; it depends on the bubble size and can be characterized by the Stokes number St, defined as the ratio of bubble buoyancy to drag force. Because bubble drag is relatively dominant over buoyancy, for a small St and small bubble size, the bubbles tend to be passively transported by the flow and persist in the water for long periods, causing long-term gas dissolution along their trajectories over a wide area. By contrast, larger bubbles (St > 1) may rapidly ascend through dominant buoyancy, limiting their contribution to gas dissolution within the short term by reaching a free surface.
The oxygen flux, S 0 , across the unit interface area is generally expressed as where k is the gas transfer velocity, C s is the oxygen concentration at the equilibrium state, and C b is the oxygen concentration in the liquid. When gas solubility in water is low, as is the case for oxygen, the gas transfer at the air-water interface is defined by the liquid-side gas transfer velocity k with C s estimated by Henry's law. The gas transfer velocity has been theoretically determined for a single rigid-shaped bubble 12,13 and experimentally modeled in terms of the bubble size and velocity. 14,15 Recently, a laser-induced fluorescence (LIF) technique has been used to measure dissolved gas concentrations. 16, 17 Herlina and Jirka 18,19 and Tsumori and Sugihara 20 investigated the dissolved gas concentration in grid turbulence near an air-water interface. Woodrow and Duke 17 found that a concentration boundary layer becomes thinner on wavy interfaces than on flat ones; larger transfer velocity is thus achieved on the resulting deformed interfaces. Valiorgue et al. 21 observed that concentration boundary layer separation in vortex wakes produced by a free-rising bubble contributes to an increase in dissolved CO 2 concentrations. Turbulence along bubble trajectories is intensified by the active motion of air bubbles, 22,23 and bubble behaviors are recursively affected by this modified turbulence. The mechanical process of bubble-laden turbulence has been numerically investigated in two major computational frameworks: Eulerian-Eulerian and Eulerian-Lagrangian approaches. 24 The former treats each phase as comprising interpenetrating continua that interact through the source terms of the momentum equations. 25,26 The latter approach discretely treats individual bubble motions governed by the Basset-Boussinesq-Oseen (BBO) equation, 27,28,24 which is advantageous to define all bubble locations as the gas source. Although most of the previous computational work in this field [29][30][31][32] has focused on the fluid dynamics of bubbly flows; we still do not fully understand size-dependent gas dissolution due to rising bubbles with turbulent behaviors, which in turn causes advection and diffusion in dissolved gases.
In this study, we focused on bubble size-dependent gas dissolution in bubble plume flows from the perspective of convection-diffusion processes. The oxygen dissolution process in rising bubble plumes, via interactions between rising bubbles and turbulent flow, was experimentally identified with the objective of modeling gas transfer velocity based on LIF measurements of DO concentration. The dissolution model was then introduced into a large-eddy simulation (LES) numerical computation, in which the local bubble motion was fully coupled with turbulence following the Eulerian-Lagrangian approach, to identify the mechanisms of gas concentration evolution.
The remainder of this paper is organized as follows. In Sec. II, we explain the experimental setup and imaging techniques used to measure oxygen transfer from air bubbles to water. The experimental features of the buoyant bubble plume flows with DO dissolution are discussed in Sec. III, introducing the new empirical model for gas transfer velocity in the plume. The bubble turbulence coupled stochastic approach using the DO transport model and LES computation is explained in Sec. IV. The computed DO concentration is discussed in terms of the dynamics of bubble flows through comparisons with the experimental results in Sec. V. All of our findings and conclusions are summarized in Sec. VI.

II. EXPERIMENTAL METHODS AND PROCEDURES
We performed laboratory experiments to simultaneously measure bubble motion and DO concentration under bubble plumes in a water tank (see Fig. 1). Oxygen gas was supplied to a needle fixed to the bottom of a transparent rectangular tank (length: 150 mm, width: 170 mm, height: 200 mm) and was ejected from the needle to generate bubble plumes. The flow rate was controlled by a discharge controller. The oxygen discharge frequency and needle outflow diameter were manipulated to produce five cases of bubble conditions (see Table I). The coordinate system for the experiment was defined with the bubble release point (needle tip 65 mm above the bottom) at the origin of the horizontal x and vertical z axes (see Fig. 1).
DO concentration in the bubble plume flows was measured using the LIF technique with fluorescence emission from pyrene butyric acid excited by ultraviolet (UV) light. Since oxygen quenches the fluorescence intensity of pyrene butyric acid, 33 fluorescent intensity, F, decreases with DO concentration, C, following the Stern-Volmer equation; where F 0 is the reference fluorescence intensity in the absence of oxygen and a is a local constant determined by a calibration procedure. Pyrene butyric acid was initially dissolved in warm 0.05-mol/l NaOH to form a 0.004-mol/l pyrene butyric acid solution, 17   to 0.000 04 mol/l for use as a DO indicator liquid to fill the experimental tank to a water depth of 170 mm. The liquid around the plume was illuminated by UV light emitted from the transparent tank bottom (see Fig. 1). Bubble forms and locations were measured using a backlit imaging method. 34 Shadow images of bubbles illuminated by a red light-emitting diode (LED) light panel (length: 80 mm, height: 70 mm) across the transparent rear wall of the tank were recorded to detect the bubble interfaces (see Fig. 2) in the manner described by Watanabe and Ingram. 34 Simultaneous fluorescence images (indicating DO concentration) and backlit images (indicating bubble interfaces) were recorded separately by two cameras via a light ray splitter (half mirror). The fluorescence emitted from the pyrene butyric acid, whose light spectrum is between 375 and 425 nm, was recorded by a low-light camera equipped with an optical bandpass filter (wavelength: 410-470 nm) at 30 Hz [see the light path (1) in Fig. 1], and the backlit images were recorded by a high-speed camera equipped with a high-pass optical filter (cut-off wavelength: 630 nm) at 500 Hz [see the light path (2) in Fig. 1]. Because light refraction at the bubble interfaces might cause uncertainties in the DO concentration estimated at bubble locations, estimates of DO concentration within the detected bubble interfaces were eliminated from the analysis to ensure the reliability of the current results. A rectangle measurement area of 40 mm × 40 mm, with a center level of z = 45 mm, was used in analyses of both image types. Image noise in both image types was reduced using a median digital filter. The experimental cases are summarized in Table I. The initial DO concentration was controlled at approximately 2.5 ppm by nitrogen bubbling for all experiments.
Additionally, we performed an experiment to measure the fluid flow in the plumes by super-resolution particle imaging velocimetry (SRPIV). 35 In these measurements, neutral buoyant fluorescent particles were seeded into the water tank, and an yttrium aluminum garnet (YAG) laser light sheet  (wavelength: 532 nm) was installed to illuminate the fluorescent particles (wavelength: 630 nm) on the x-z plane. A high-speed camera recorded the fluorescent-particle displacement at 500 Hz. Laser light reflected from the bubble surfaces was filtered using a high-pass optical filter (cut-off wavelength: 590 nm) attached to the camera. Figure 2 displays typical backlit bubble images with detected interfaces in the measurement area where the bubbles ascended at a terminal velocity. We observed changes in bubble shapes over successive bubble ascents, especially in larger bubbles (cases [3][4][5]. Because deformation results in fluctuations in the ascent velocity due to changes in the drag on the bubble, larger ascent velocity dispersion occurs in cases involving larger bubbles (see Table I). According to Woodrow and Duke, 17 the deformation of the air-water interface increases the gas transfer velocity k [see Eq. (1)], such that size-dependent deformation locally modifies the gas dissolution process.

A. Bubble behaviors
Bubble shapes were classified using the well-known Grace diagram 36 in terms of the Morton M = gµ 4 /ρ f γ 3 , Eotvos Eo = gd 2 ρ f /γ , and Reynolds numbers Re = ρ f du r /µ (see Fig. 3), where g, µ, ρ f , γ, d, and u r are the gravitational acceleration, molecular viscosity, liquid density, surface tension, bubble diameter, and relative bubble velocity, respectively. The experimental results of the current study fall into the ellipsoidal shape and wobbling regimes of the Grace diagram. Figure 4 shows the time records of the aspect ratio of the major to minor axes of the ellipse-shaped bubbles (A). When the value of A was in the range 0.5-0.6 for small bubbles (cases 1 and 2), the fluctuation range increased significantly with the mean bubble diameter; A ≈ 0.4-0.8 for cases 3 and 4, 0.3-0.97 for case 5. Bubble deformation may also modify local turbulence and the pressure field around the bubble to deform the  bubble shape again, resulting in changes in the orientation of bubble motion along typical zigzag trajectories to extend the area in which oxygen can be dissolved, which is discussed in Secs. III B and III C. Figure 5 shows the measured concentration, together with the bubble locations estimated from the sequential high-speed backlit images. As the concentration gradually increased along the zigzag bubble trajectories, concentration increased within the range covering lateral bubble displacements. Because oxygen dissolved in the plume may be transported vertically by upward bubble-driven flows to be circulated in the tank, the background concentration gradually increased in time as well as within the plume range. For the largest bubble (case 5), the concentrations in both the bubble plume and the background rapidly increased compared with those of case 1, due to the larger interface area where oxygen could be dissolved into bulk water. The separation of a concentration boundary layer, formed around large bubbles, 21 may also be relevant to the increased concentration rate, which will be discussed in Sec. III C. Figure 6 shows the horizontal distribution of concentrations at the center level of the measured plane (z = 45 mm). We found that the horizontal concentration was well approximated by a Gaussian distribution with increasing background concentration during the measurement period, Fig. 6). The three parameters, reference concentration at the center C p , background concentration C b , and standard deviation σ were determined by least-squares approximation, are shown in Fig. 7. Because the three parameters depend on the bubble size and gas discharge, C p and σ were almost constant in time (C p ≈ 4.28 ± 0.35 ppm and σ ≈ 2.76 ± 0.18 ppm throughout the measurement period in case 5), that is, the horizontal profile of the concentration was temporarily unchanged, whereas C b monotonically increased in time, approaching the saturation concentration (about 39 ppm).

B. DO concentration around bubbles
The oxygen flux due to transport from a bubble to bulk water was estimated by the temporal increment of the approximated maximum concentration, Figure 8 shows the time record of the increasing rate of C max for a single bubble, n −1 b dC max /dt. The rate increases with the bubble size because larger bubbles have larger interface areas. Steady gas flux was observed in smaller bubbles (cases 1 and 2), whereas temporal decreases in the rates of increase were observed in larger bubbles (case 3, 4, and 5). Although the gas flux per unit interface area is generally described in proportion to the relative concentration C s − C b [see Eq. (1)], the background concentration C b temporally increased to approach the equilibrium state; thus, C s − C b decreased, especially in larger bubbles with higher C b , resulting in a temporal decrease in gas flux under a constant gas transfer velocity k.

C. Oxygen transfer velocity from a bubble to water
Gas transfer across bubble interfaces has been used to derive the gas transfer velocity, k, in previous studies. Lewis and Whitman 37 proposed the well-known "film model" where D is the molecular diffusivity of oxygen in water and δ c is the concentration boundary layer thickness. Higbie 38 developed the "penetration model," which included the effect of the surface renewal process, leading to the following relationship: 15 where Sc is the Schmidt number, defined as the ratio of the kinematic viscosity of the liquid ν to the molecular diffusivity D. Levich 12 theoretically derived the mass transfer from a bubble for creeping flow Finally, Frössling's equation 13 for a rigid sphere, often used for small bubbles, is given by In the current study, the oxygen transfer velocity k c was explicitly estimated on the basis of Eq. (1) with the measured oxygen flux (Fig. 8), mean bubble surface area, and relative concentration C s − C b . The least-squares approximation for the current results yields where the size-dependent constant α ≈ 1.13 for d < 3.6 mm, and 2.21 for d ≥ 3.6 mm. Figure 9 shows the estimated k c as a function of D/d with comparisons to the previous models, k H , k L , and k Fr , assuming that the bubble ascends at a terminal velocity. We found that k L and k F were not acceptable over the whole range of D/d for the current plume flow. For smaller bubbles (cases 1 and 2), the current model, k c , was consistent with k H , whereas there was a significant deviation between k c and k H when D/d < 6. In the current study, a gap was observed in the transfer velocity around D/d ≈ 6, and the transfer velocity dispersion at D/d < 6 was much larger than those for larger D/d values. This result may have been caused by the dependency of the deformation process on the bubble size. As discussed above, in larger bubbles (D/d < 6), the amplitude and velocity of the interface displacement were much larger than those of smaller bubbles (see Fig. 4). Valiorgue et al. 21 observed that an ascending CO 2 bubble leaves trails of high-CO 2 concentration from the bubble bottom due to vortex wakes. Thus, high concentrations persist along the bubble trajectory. Figure 10 shows measured distributions of turbulent kinetic energy, E k = u · u /2, where the fluctuation velocity u = u − u, u is instantaneous velocity, and u is the temporal mean velocity. We observed higher turbulence energy in larger bubbles; in particular, turbulence was significantly intensified when a bubble ascended along the plume axis at D/d < 6 (cases 3-5), which may have been the result of vortex wakes induced by the ascent of the large bubble. Accordingly, the concentration boundary layer formed at the deformed bubble surface may have trailed along the lee side of the bubble and been shed by vortex wakes in the turbulent flow. The separated boundary layer within the wake remained on the plume axis, and a boundary layer is renewed on the bubble bottom. This repeated separation and renewal of the concentration boundary layer on the bubbles can significantly increase the gas transfer velocity in larger bubbles. Although previous analytical models, which assumed laminar flow, failed to predict distinct gas transport for large-bubble plumes, the current model empirically includes this effect through the bubble-turbulence interaction (see Fig. 9). This effect may also contribute to maintaining the high gas flux observed in Fig. 8 for larger bubbles.
While the current paper focuses on the mass transfer of oxygen, Eq. (7) may be applicable to dissolution processes for the other gases with different S c and D. Introducing the Prandtl number and thermal diffusivity instead of S c and D, the model may be extended to heat transfer problems.

IV. NUMERICAL MODEL
The current empirical model of bubble-size-dependent oxygen transfer was coupled with stochastic Eulerian-Lagrangian computation to comprehensively predict oxygen concentration via dissolution, diffusion, and convection processes in the bubble plume. Turbulent flow was computed in the LES framework coupled with the individual bubble motion description from the BBO equations. A generalized Langevin equation was used to define the sub-grid liquid velocity at the bubble location to predict the local liquid force acting on the bubble.

A. Liquid phase
Assuming dilute bubble flow, the filtered Navier-Stokes equation for incompressible fluid was used as a governing equation where u f ,i is the liquid velocity, ρ f is the liquid density, p f is the fluid pressure, g i is the gravity vector, τ ij is the subgridscale (SGS) stress, and Q i is the bubble stress. The overline indicates grid-scale (GS) variables. The GS bubble force, Q i , which is determined using a top-hat filtering operation over a computational cell, may be identical to the mean force due to bubbles included in the cell where ∆ is the grid spacing, N is the number of bubble centers contained in the cell volume, and F D,i (n) and F A,i (n) are the drag and added mass forces of the nth bubble. Eddy viscosity SGS stress was defined as where the GS eddy viscosity ν e = c ν ∆ √ q sgs , δ ij is the Kronecker delta, and the SGS kinetic energy is represented by A well-known Smagorinski model is based on the assumption of the balance between the SGS energy production and dissipation (local equilibrium hypothesis). This energy balance holds in isotropic turbulence, but may break down in turbulence involving jets and wakes (Ref. 39). Yoshizawa and Horiuti 39 statistically derived a SGS kinetic energy LES model on the basis of two-scale direct-interaction approximation, which re-sorts the energy balance assumed in the Smagorinski model, for modeling anisotropy of turbulent intensity. We used the SGS kinetic energy model by Yoshizawa and Horiuti 39 to derive the SGS kinetic energy equation for bubble-laden flow governed by Eq. (8) where the dissipation rate = c q 3/2 sgs /∆, c = 1.8, and c ν = 0.043.

B. Bubble dispersion phase
The motion of a bubble with equivalent volume πd 3 /6 is governed by the BBO equation 27,40,41 where u b,i is the instantaneous bubble velocity at the bubble location x b,i . The drag force F D,i , added mass F A,i , Lagrangian acceleration F L ,i , and gravity (and buoyancy) F G,i are given by and C D is the drag coefficient, C A is the added mass coefficient, ρ b is the air density, and u s,i is the instantaneous velocity of the fluid seen by the bubble. 42 The derivatives du b,i /dt and Du s,i /Dt represent the accelerations of the bubble and fluid, respectively. The so-called Faxen terms, accounting for nonuniformity of the flow field, 27 were neglected in F A,i , assuming the relative small bubble to the length scale of the ambient flow. 41 The drag force depends on the bubble shape, in addition to fluid viscosity, bubble diameter, and surface tension, which has been modeled in terms of the Reynolds number (Re) 12,43 and Eotvos number (E o ). 44 Tomiyama et al. 45 combined the previous C D models to derive the universal empirical model covering the ranges of 10 −2 < E o < 10 3 and 10 −3 < Re < 10 5 C D = max min 16 Re This model has been confirmed to be consistent with the previous drag coefficient 44 and terminal velocity 46,47 for a nonspherical bubble, which ensures that Eq. (17) is acceptable for bubbles with ellipsoidal and wobbling shapes observed in the current experiment (see Fig. 3). The effect of the Basset history term has been investigated by Hjelfelt and Mockros, 48 Coimbra and Rangel, 49 and Vojir and Michaelides. 50 According to Vojir and Michaelides, 50 the effect of the history term is negligibly small when the dimensionless frequency of velocity fluctuations is less than 0.5. Assuming that the fluid velocity fluctuation at the bubble location synchronizes with the bubble shape oscillation causing wobbling motion, the normalized characteristic frequency, f, with respect to the bubble relaxation time, τ b = ρ b d 2 /18µ, was estimated from variations of the bubble aspect ratios, observed in Fig. 4, for evaluating effectiveness of the history term. The estimated f for experimental case 1-5 are f 1 = 0.0148, f 2 = 0.0285, f 3 = 0.0180, f 4 = 0.0245, and f 5 = 0.1300.
In all cases, the characteristic frequency is less than 0.5, so the history term is unimportant in the BBO equation and is not considered in the computation.
According to Bataille, 51 the added mass coefficient C A = 0.5 is approximately valid even for moderately ellipsoidal bubbles with diameters up to 3.5 mm. Therefore we used C A = 0.5 for the computational cases with the bubble diameter less than 3.5 mm.
It should be noted that Eq. (12) requires the assumption that the relative velocity between the bubble and fluid is sufficiently small; that is, the bubble Stokes number St = τ b /τ f should be smaller than one where τ f is the time scale of fluid flow. St for all experimental cases are summarized in Table I. We find that St for all cases are less than one, and Eq. (12) is thus applicable to all experimental bubble regimes.
The bubble location x b,i was updated with the computed u b,i , given by The fluid velocity u s,i at an arbitrary sub-grid location is uncomputable in the standard LES framework; therefore, a stochastic representation of u s,i from the generalized Langevin equation [52][53][54] was used in this study where the constant C 0 = 6.2 (Ref. 55) and W i represents the vector-valued Wiener process. represents the mean quantity at the bubble location, which is estimated by quadric interpolations of the filtered variables. The Lagrangian integral time scales T L and T * L,i are given by Csanady. 42 The turbulent flow computation for bubble-fluid twoway interactions is achieved by the simultaneous solution of Eqs (8), (11), (12), and (19). Since the drag coefficient (17) is dynamically changed with local Re and E o according to the fluid turbulence described by Eq. (19), probabilistic disturbances of the drag force cause lateral dispersions of the ascending bubbles with typical zigzag motions, which will be explained in Appendix B.

C. DO concentration
The turbulence diffusion and advection processes of passive concentration are governed by the filtered advectiondiffusion equation in the LES framework where S c is the Schmidt number and S c,sgs is the turbulent Schmidt number. S is the filtered source term of the oxygen dissolved from bubbles contained in a computational cell where the oxygen transfer from a bubble, S n , is defined by S n = πd 2 n k c (C s − C). (22) k c is the oxygen transfer velocity proposed in the current experiment, Eq. (7), and C s is the DO concentration for the equilibrium state, defined by where T w is the water temperature, R is the ideal gas constant, α is the Bunsen solubility coefficient given by Weiss, 56 and p o is the oxygen partial pressure. In the current case, where oxygen occupies the bubbles, the oxygen partial pressure p o is identical to the total pressure within the bubble p b , expressed by where p a is the atmospheric pressure, h is the water depth, and γ is the surface tension.
In this study, we assumed no gas exchange except across the bubble interfaces because the gas transfer velocity at the still water surface was much smaller than that at the bubble interface. 17

D. Numerical methods
The numerical methods used in the current computation are briefly explained in this section. Full details may be found in Watanabe et al. 3,57 Within the framework of a fractional step method, discretized forms of Eqs. (8) and (11) were decomposed into advection and non-advection equations. The cubic interpolation polynomials method 58 was used to solve the advection equation, and a predictor-corrector method was used to update the non-advection equation. A Poisson equation for pressure was iteratively solved using a multi-grid method. A level-set method 59,60 was used to compute the deformation of the free surface in the water tank, and spherical SGS bubbles were tracked by solving Eq. (18) without computing the local bubble interfaces. The current numerical techniques used to simulate liquid turbulence were validated in Watanabe et al. 3,57 The computation was performed in a Cartesian staggered grid system with a dimensional resolution of 5.0 mm per cell, in a rectangular domain with non-slip side walls and bottom. The measures of the domain were identical to those of the experimental water tank (see Fig. 1). The dimensional time step interval was 5.27 × 10 −4 s.
The Eulerian-Lagrangian approach has the fundamental assumptions that the bubble is isolated and the bubble-bubble interaction is excluded; that is, the dilute flow condition should be fulfilled. The bubble volume contained in a computation cell needs to be sufficiently smaller than the cell volume to achieve the dilute condition, while high resolution is required for resolving smallest vortices formed in the plume. The smallest bubble case (with 3.4% volume fraction of a bubble in a computational cell) was selected to be computed in this study, as small-scale turbulence may be unable to computed for larger bubble cases with coarser cells. Numerical bubbles 2.0 mm in diameter were ejected at a constant rate (50 bubbles per s) in experimental case 1 from the same level as in the experiment. Ascending bubbles arriving at the free surface were removed and excluded from further fluid computation.

V. COMPUTATIONAL RESULTS
In this section, the major features of the computed bubble plume flows are explained, and the computed DO concentration, based on the current oxygen transfer model, is then discussed through comparisons with the experimental results. The computed turbulent flows are validated in Appendix A. Figure 11 shows the evolution of the fluid velocity driven by the bubble plume. Upward liquid flow is induced by the ascent of successive bubbles along the plume axis during the early stage [see Fig. 11(a)]. The bubble-driven liquid lifts the free surface above the plume to produce divergent surface flows toward the side walls [see Fig. 11 . The major quasi-steady circulation system is observed above the bubble release level (z/d = 0), which may contribute to the convection of oxygen dissolved in the plume. It should be noted that the current computation provided consistent estimates of the experimental bubble ascent velocity and liquid velocity within the tank, validating the current model (see Appendix A). Figure 12 shows the iso-surface of SGS turbulent kinetic energy together with the bubble locations. We found that turbulent energy is produced above the level where the vertically accelerated bubbles achieved terminal velocity in the plume (z/d = 6.0) and intensified along with the ascending bubbles. Bubble-induced turbulence was transported to the free surface and radically spread through surface divergent flows.

A. Bubble turbulence flows
We discuss how the bubble plume flow contributes to the oxygen transport process in Sec. V B. Figure 13 shows the evolution of the gas concentration in the tank. The concentration was observed to increase along the plume axis, where the bubbles successively supplied oxygen during the early stage (a). Oxygen-rich liquid was transported by bubble-induced vertical flow and then radically spread through divergent surface flow [ Fig. 13(b)]. Vortices formed near the side walls convected the liquid, increasing the background concentration in the circulation flow system [ Fig. 13(c)]. The concentration along the bubble plume was temporally increased by the cumulative oxygen supply in the plume, which is identical to the experimental observation (see Fig. 5). Figure 14 shows the time records of the computed concentration at z/d = 22.5, which corresponds to the center level of the experimental measurement area (see Fig. 5), in comparison with the experimental plots for case 1. The computed concentrations within the plume, x/d = 2.5, and outside of the plume, x/d = 5.0 and 7.5, monotonically increased with some fluctuation after t ≈ 50 s once a quasi-steady state of circulation flow was attained (Figs. 11 and 13). The concentration at x/d = 2.5 was consistently higher than that outside of the plume, x/d = 5.0 and 7.5 (see also Fig. 6), as the oxygenrich liquid from the plume was delivered to the outside of the plume at a lag time of 50 s via circulation flows.

B. DO concentration in the bubble plume
The instantaneous concentration may be decomposed into slowly varying and rapidly varying components as C = C +C .
The slow term C may be evaluated as the temporal movingaverage of C, and the rapid C indicates fluctuations associated with bubble ascent and the resulting turbulence. Assuming that a uniform concentration in the tank with water volume V w is achieved by a steady supply of oxygen dissolved from n b bubbles of diameter d, the increasing rate of C(t) may be described by The solution for the above equation is where C 0 is the initial concentration. Equation (26) demonstrates the slowly increasing background concentration (see Fig. 14). Fluctuations in concentration, C , are highly dependent on the horizontal location; that is, high-frequency energy fluctuations are involved in the bubble plume (x/d = 2.5 in Fig. 14) and slower variations are observed outside of the plume (x/d = 5.0, 7.5). Figure 15  dissolved oxygen in the plume, transported by circulating flow, returned to the plume location after about 50 s (see Fig. 13). Identical long-term oscillations were observed in u/U T and w/U T , which in turn modified the transfer velocity k c through Eq. (7). As previously discussed, the increasing concentration rate was larger in large-bubble cases because concentration boundary layers around large bubbles became separated during their higher-velocity ascent due to lee vortex wakes shed from the bubbles. The boundary layer then remained along the plume axis. The gas transfer velocity k c , defined in Eq. (7), takes this effect into account empirically and describes the evolution of the mean (background) concentration in the tank consistently. Figure 16  velocity k c , Eq. (7), Higbie's k H , Levich's k L , and Frössling's k F . The Higbie model and the current model well predicted the slowly increasing background concentrations in small bubbles (cases 1 and 2), and the other models k L and k F underestimated the measured concentration. For larger bubbles, all previous models underpredicted the experimental results, and the proposed model provided a reasonable prediction.

VI. CONCLUSION
In this study, size-dependent gas transfer from bubbles into bulk water in bubble plume flows was identified through an image analysis of air-water two-phase flow and DO concentration. The temporal increasing rate of the DO concentration in the plume was highly dependent on the bubble size, which defined the terminal ascent velocity and the intensity of bubbleinduced turbulence, as well as the interface area. For flows with large Reynolds numbers (indicating large bubble size and ascent velocity), the concentration boundary layer for the oxygen-rich liquid, which formed around bubble interfaces, was separated behind the ascending bubbles to remain on the plume axis, resulting in high concentration and high turbulent energy along the plume axis. This effect causes significant deviation in the oxygen transfer velocity from those of previous models in larger bubbles (D/d < 6), whereas the transfer velocity for small bubbles was consistent with those determined by the Higbie 38 model. This study therefore proposes a new empirical model of size-dependent transfer velocity in bubble flows.
The current model was introduced into an LES-stochastic, bubble-turbulence-coupled computation to determine the dissolution, diffusion, and advection processes in an experimental turbulent bubble plume. The concentration in the plume was observed to vary on three major time scales. The fastest variations were correlated with the local fluid velocity and turbulent energy, indicating that bubble-induced turbulence governed fluctuations in concentration within the plume. Long-term oscillations were associated with the transport of oxygen-rich liquid in the circulation flow system, where bubble-driven divergent surface flows changed their orientation relative to the tank walls, to return to the plume. The slowest variations, indicating net transport of oxygen from bubbles to the tank water, were observed to monotonically increase toward the saturation concentration.
The mass conservation of DO with the proposed transfer velocity can be used to estimate the monotonic increase in the slowly varying background concentration. The current and Higbie models predicted the experimental concentrations for small bubbles (D/d > 6). In this study, the current model provided reasonable predictions regardless of the bubble size, whereas previous models underestimated the measured concentrations for plumes with large bubbles (D/d < 6).
Limitations of the Eulerian-Lagrangian approach for the bubble plume computation are associated with a dilute flow condition that constrains the bubble size and grid width.
Further investigation is required to determine practical applications for plumes with the wide range of bubble sizes observed in sea surface layers.

ACKNOWLEDGMENTS
This research was supported by the Japan Society for the Promotion of Science (JSPS) Research Fellowships for Young Scientists (No. 11J04354) and a JSPS Grant-in-Aid for Scientific Research (No. 18H03791).

APPENDIX A: VALIDATION OF THE COMPUTATIONAL MODEL
In this section, the computed bubble ascent velocity, liquid velocity, and liquid turbulent energy are compared to the experimental results to validate the LES stochastic model presented in this study.
First, because the bubble motion was stochastically determined in the computation, the stochastic representations of computed results are compared with the experimental results.  Figure 17 shows the probability density distributions of the computed and experimental bubble ascent velocities W b , normalized by the terminal ascent velocity for a single bubble U T (case 1). We found reasonable agreement between the computed and experimental distribution of the ascent velocities, indicating the validity of the current computation. We also found maximum provability at W b /U T ≈ 1.1, indicating that bubble ascent in the plume W b was faster than that of a single bubble in still water U T because bubble ascent was followed by bubble-driven vertical flows. Figure 18 shows the horizontal distributions of the ensemble-averaged horizontal liquid velocity u, vertical liquid velocity w, and turbulent energy E k at z/d = 22.5 (case 1). The computed water velocity reasonably described the experimental bubble-driven plume flow in the circulation flow system, and the computed turbulent energy underestimated the experimental turbulent energy around the plume axis.

APPENDIX B: BUBBLE TRAJECTORIES
The stochastic bubble-turbulence coupled computation has been performed only for the smallest bubble case (case 1) since the computation for larger bubbles requires coarser computing cells to achieve the dilute flow condition required by the Eulerian-Lagrangian approach. The turbulence defined by coarse grids in the LES framework may be inconsistent with turbulence statistics. In this section, instead of the constrained LES, the experimental flow field measured by SRPIV was used for computing the Lagrangian bubble dispersion phase for all experimental cases. While steady mean velocity was given for computing Eqs. (12), (18), and (19), the stochastic turbulence component estimated by (19) modifies the fluid forces (13)-(15) as well as the drag coefficient (17), causing lateral fluctuations of the ascending bubbles. Figure 19 compares the experimental and computational trajectories of ascending bubbles. We find the lateral dispersions of the ascending bubbles with wobbling behaviors in the computed trajectories. The computed width of the lateral dispersions is consistent with the experimental one especially in larger bubble case 3-5, while significant difference is found in the dispersion width between the experimental and computational trajectories in smaller case 1 and 2. Tomiyama et al. 61 found that the lateral displacement of the bubble trajectories highly depends on the inner diameter of the needle to create bubbles as gas pressure in the needle affects the initial bubble forms. The wide dispersion observed in the experimental case 1 and 2 may be caused by the initial bubble deformation induced by this effect as fine needles were used for these two cases (see Table I). Overall features of the computed dispersion process agree with the experimental ones, suggesting our valid model for the Lagrangian bubble computation.