Numerical Simulation of Micro-Bubbles Dispersion by Surface Waves

: This paper presents an algorithm for numerical modeling of bubble dispersion occurring in the near-surface layer of the upper ocean under the action of non-breaking two-dimensional (2D) surface waves. The algorithm is based on a Eulerian-Lagrangian approach where full, 3D Navier-Stokes equations for the carrier ﬂow induced by a waved water surface are solved in a Eulerian frame, and the trajectories of individual bubbles are simultaneously tracked in a Lagrangian frame, taking into account the impact of the bubbles on the carrier ﬂow. The bubbles diameters are considered in the range from 200 to 400 microns (thus, micro-bubbles), and the effects related to the bubbles deformation and dissolution in water are neglected. The algorithm allows evaluation of the instantaneous as well as statistically stationary, phase-averaged proﬁles of the carrier-ﬂow turbulence, bubble concentration (void fraction) and void-fraction ﬂuxes for different ﬂow regimes, both with and without wind-induced surface drift. The simulations results show that bubbles are capable of enhancing the carrier-ﬂow turbulence, as compared to the bubble-free ﬂow, and that the vertical water velocity ﬂuctuations are mostly augmented, and increasingly so by larger bubbles. The results also show that the bubbles dynamics are governed by buoyancy, the surrounding ﬂuid acceleration force and the drag force whereas the impact of the lift force remains negligible.


Introduction
Understanding the dynamics of small-scale processes occurring in the vicinity of the sea surface is important for modeling the exchange of mass, heat and momentum between the atmosphere and the ocean [1,2]. In many practical situations, clouds of bubbles, produced mainly by breaking surface waves, may also affect the state of the near-surface water layer [3,4]. Laboratory and field observations [3,[5][6][7] as well as recent direct numerical simulations of breaking waves [8] indicate that whereas comparatively large bubbles (with diameters d > 1 mm) quickly rise to the surface and burst, smaller (or micro-) bubbles (d~100 µm) are suspended in water for a considerably longer times and thus contribute mostly to the void fraction observed in the near-surface oceanic layer.
According to field observations the concentration of microbubbles (and thus the void fraction) in the upper ocean layer can be considerable even at relatively low winds and affects gas exchange between the air and water [9,10], production of sea-salt aerosols [11], and the propagation of sound in the upper ocean [12,13]. Therefore, modeling microbubble dynamics is important for many practical applications.
Numerical modeling of the dispersion of microbubbles by surface waves and currents taking into account their impact on the near-surface turbulence represents a challenging problem, and various models are employed in numerical investigations of bubbly flows [14]. Known numerical studies of microbubble dispersion are mainly restricted to either bubbleladen isotropic turbulence (with periodic boundary conditions) [15], or boundary-layer flows in the vicinity of a fixed, flat boundary [16]. However, performing numerical simulation of a flow in the vicinity of a waved boundary involves additional efforts required to resolve a strong geometric nonlinearity. There are mainly two methods employed to cope with this problem in numerical studies. One approach is related to using the Volume of Fluid method where the air and water phases are directly resolved (cf., e.g., [8]). However, in order to model a sufficiently high void fraction (on the order of 10 −5 [12]), the number of microbubbles to be considered in the framework of a fully-resolved numerical experiment with an adaptive mesh refinement may become prohibitively large. Another approach (also used and discussed in the present study) employs a mapping of the physical coordinates onto the curvilinear coordinates where the waved surface is reduced to a flat surface.
The objective of the present paper is to present a numerical algorithm for evaluation of the dispersion of micro-bubbles by progressive, non-breaking surface (Stokes) waves. The wave shape is prescribed and assumed to be stationary, and unaffected by either bubbles or induced turbulent motions. Thus it is assumed that for typical void fractions observed in the upper ocean layer (on the order of O(10 −5 ) or less), the impact of bubbles on the carrier surface wave flow remains negligible. Results of numerical experiments also reveal that typical amplitudes of turbulent wave-induced fluctuations are by orders of magnitude smaller as compared to the energy-containing mother wave [17]. However, in the present study, the impact of bubbles on the induced turbulence, although not so significant at the considered void fractions, is accounted for. The algorithm also allows to investigate how wind-induced surface drift affects bubbles dynamics employing a parameterization of the drift velocity at the water surface in terms of the air friction velocity (determined by the velocity scale tied to the surface wave celerity) [18]. The numerical method is based on the algorithms previously developed for modeling atmospheric boundary layer over progressive surface waves under various conditions (stable air-stratification, parasitic capillaries, droplets) [19][20][21][22][23] and adapted here for modeling a bubble-laden near-surface water layer. boundary-layer flows in the vicinity of a fixed, flat boundary [16]. However, performing numerical simulation of a flow in the vicinity of a waved boundary involves additional efforts required to resolve a strong geometric nonlinearity. There are mainly two methods employed to cope with this problem in numerical studies. One approach is related to using the Volume of Fluid method where the air and water phases are directly resolved (cf., e.g., [8]). However, in order to model a sufficiently high void fraction (on the order of 10 −5 [12]), the number of microbubbles to be considered in the framework of a fully-resolved numerical experiment with an adaptive mesh refinement may become prohibitively large. Another approach (also used and discussed in the present study) employs a mapping of the physical coordinates onto the curvilinear coordinates where the waved surface is reduced to a flat surface.

Governing Equations
The objective of the present paper is to present a numerical algorithm for evaluation of the dispersion of micro-bubbles by progressive, non-breaking surface (Stokes) waves. The wave shape is prescribed and assumed to be stationary, and unaffected by either bubbles or induced turbulent motions. Thus it is assumed that for typical void fractions observed in the upper ocean layer (on the order of O(10 −5 ) or less), the impact of bubbles on the carrier surface wave flow remains negligible. Results of numerical experiments also reveal that typical amplitudes of turbulent wave-induced fluctuations are by orders of magnitude smaller as compared to the energy-containing mother wave [17]. However, in the present study, the impact of bubbles on the induced turbulence, although not so significant at the considered void fractions, is accounted for. The algorithm also allows to investigate how wind-induced surface drift affects bubbles dynamics employing a parameterization of the drift velocity at the water surface in terms of the air friction velocity (determined by the velocity scale tied to the surface wave celerity) [18]. The numerical method is based on the algorithms previously developed for modeling atmospheric boundary layer over progressive surface waves under various conditions (stable air-stratification, parasitic capillaries, droplets) [19][20][21][22][23] and adapted here for modeling a bubble-laden near-surface water layer.  A domain with sizes 6 , 4 ,  A domain with sizes L x = 6λ, L y = 4λ, L z = λ, with periodic side boundaries, a solid bottom and a waved upper boundary, z s (x, t), is considered, and the Cartesian coordinates are employed, 0 ≤ x ≤ L x ; −L y /2 ≤ y ≤ L y /2; −L z ≤ z ≤ z s . The motion of the fluid (water) is driven by the upper boundary where a progressive, two-dimensional (2D) stationary wave of amplitude a, celerity c and length λ propagating in the positive x-direction is prescribed. A Eulerian-Lagrangian approach is adopted where the Navier-Stokes equations of the water motion are solved in a Eulerian frame, and the bubbles are tracked simultaneously by solving their respective equations of motion in a Lagrangian frame.

Governing Equations
The Navier-Stokes equations for the carrier fluid are written in the dimensionless form [24]: where x i ≡ (x, y, z), U j = (U x , U y , U z ) are water velocity components, P is the pressure, and S n i is a momentum source term (cf. Equation (8) below) contributed by the n-th bubble, and n = 1, . . . , N b , the latter being the total, constant, number of tracked bubbles. Variables in Equation (1) are normalized with velocity scale, U 0 , set equal the surface wave phase speed, c, and length scale, L 0 , equal to the wave length, λ. The pressure is normalized with ρU 2 0 where ρ is the water density (≈1 g/cm 3 ). The carrier flow Reynolds number, Re, measures the ratio of the inertial vs. viscous effects and is defined as: where v is the water kinematic viscosity (≈0.01 cm 2 /s). Equation (1) is supplemented by the incompressibility condition: which implicitly defines the pressure field, P (cf. Equation (18) below). The bubbles equations of motion are written as [15]: In Equations (4) and (5), r n i , V n i (i = x, y, z) are the n-th bubble coordinate and velocity components, U n i is the fluid velocity at the bubble location, d/dt = ∂/∂t + V n j ∂ j and D/Dt = ∂/∂t + U n j ∂ j are the material derivatives along the bubble trajectories and the surrounding fluid Lagrangian paths, respectively; g is the dimensionless gravitational acceleration, and ω n i = ε ijk ∂ j U n k is the surrounding fluid vorticity; δ ij and ε ijk are the Kronecker and Levi-Civita tensors. The forces acting on the bubble (in the order as in the right hand side of Equation (5)) are the fluid acceleration, viscous drag, lift, and buoyancy. The correction in the viscous drag force (accounted for by factor f ) is caused by a finite Reynolds number of the bubble and defined as: where [15], Both the momentum source term in the right hand side of Equation (1), S n i , and Equation (5) are formulated under an assumption of a negligible mass of the bubble as compared to water. A "point force" approximation is adopted, whereby the following formulation for the source term, S n i , is employed [23]: where w(r n , r) is a geometrical weight-factor inversely proportional to the distance betweeen the n-th bubble located at r n = (x n , y n , z n ) and the grid node at r = (x, y, z), and Ω g (r) is the volume of the considered grid cell. Thus, for each individual bubble, eight weightfactors are defined (for each of the surrounding grid nodes) and normalized, so that the sum of partial contributions distributed to these nodes exactly equals the respective total source contribution. Therefore, there is no numerically induced loss or gain of momentum in the bubble-water exchange processes (cf., e.g., [23] and references therein for a more detailed discussion).

Numerical Method
In order to avoid coping with a strong geometric nonlinearity caused by the wavy upper boundary during the integration of the governing equations discussed above, a mapping is introduced transforming the domain with a wavy upper boundary into a domain with a flat upper boundary, and relating the Cartesian coordinates (x, z) to curvilinear coordinates (ξ, η) as: Mapping (9) transforms the wavy (upper) boundary at z = z s (x, t) into a flat-plane boundary at η = 0. The water surface elevation, z s (x, t), is defined implicitly by Equation (9) and up to the second order in ka coincides with the Stokes-wave solution [25]: An additional mapping is also employed for the vertical coordinate in the form: where −1 ≤ η ≤ 0, so that the grid nodes are clustered in the vicinity of the upper boundary, at η = η = 0, and stretched with increasing depth. Since the mapping, Equation (9), is conformal, the following relations between the derivatives hold: where the Jacobian of the transformation is: Due to the properties in Equations (12) and (13), the derivatives over the Cartesian coordinates, x and z, can be related to the derivatives over curvilinear coordinates, ξ, η, as, The Laplacian operator is also rewritten as: Equations (1) and (3) for the water velocity are discretized on a staggered grid consisting of 360 × 240 × 180 nodes, in the ξ, y and η coordinate directions using a secondorder-accuracy, finite-difference method. At each time moment, t k , the mesh is redefined and adapted to the shape of the water surface, z s (x, t k ), according to Equation (9), and all The integration of Equation (1) is advanced in time by the second-order-accuracy Adams-Bashforth method in two stages to calculate the water velocity at each new time step, U i (t k+1 ). First, an intermediate velocity, U * i , is computed using the velocity fields at the preceding time steps [26]: where the flux, F i , is evaluated as, Further, the new pressure, P(t k+1 ), is computed by solving its respective Poisson equation in the form, Equation (18) is solved by iterations by performing, at each iteration step (j) the FFT in the horizontal directions and Gaussian elimination in the vertical direction. The iteration procedure stops when the condition P j+1 − P j /P j < 0.1% is satisfied. Usually this condition is met after j = 3-5 iterations. The new velocity at k + 1 time step satisfying the incompressibility condition (2) is then computed as: At the upper boundary, η = η = 0, the no-slip (Dirichlet) condition for the velocity is prescribed: where the surface drift velocity, U d , is expressed as [27]: Two different cases are considered. In one case, parameter q is put to zero, U d = 0, and the wind stress effects are not taken into account. In this case, the water velocity at the boundary coincides with orbital velocities of the fluid particles in the surface wave. In another case, q = 0.05c is prescribed, and thus U d is finite so that the wind-stress effects upon the water surface are accounted for.
Periodic conditions for all fields are prescribed at the side boundaries, and the no-slip (Dirichlet) condition for the water velocity is prescribed at the bottom boundary ( η = −1).
The bubble equation of motion is solved by employing the Adams-Bashforth method: In Equation (23), the surrounding fluid velocity, its acceleration, and the ambient-flow vorticity at the location of each bubble are obtained by a Hermitian 4th-order accuracy interpolation procedure. The dimensionless bubble response time, τ n , (or the Stokes number) is expressed as: The bubble coordinate equation is advanced in time by employing the Adams method as: The water velocity is initialized as a random, zero-mean field with an amplitude of 0.1% (normalized by the velocity scale, U 0 ). After initiation, a transient occurs during which the velocity field adjusts to boundary conditions, and a statistically stationary flow state is reached (at dimensionless time t = 100). The results (not shown) indicate that the the mean, wave velocity field is established comparatively early (~O(10) wave periods) and the transient to the stationary flow state is mainly related to the development of a near-surface, wave-induced turbulent layer (of~O(100) duration, cf., e.g., [17]). After the stationary flow state is reached, the bubbles are injected into the flow at random locations with a concentration (number density) distribution exponentially decreasing with depth with an e-folding scale close to the surface wave length (similar to void-fraction distributions observed in natural, oceanic conditions [12]). Since bubbles rise due to buoyancy, they reach the upper boundary (i.e., the water surface) and thus leave the computational domain. In order to maintain a constant void fraction throughout the simulation, these bubbles are re-injected at random locations with a spatial distribution exponentially decaying with depth and a velocity equal to the surrounding water velocity.
The simulation of the bubble-laden flow continues until the flow again reaches a stationary state (at t = 250) where its statistical characteristics are evaluated. Similar to the previous DNS studies of flows over waved surfaces [19][20][21][22][23]27,28]), in the statistical post-processing analysis, phase averaging, equivalent to averaging over an ensemble of turbulent fluctuations, is performed. This averaging (denoted below by angular brackets) is firstly performed over the y-coordinate and time t, and further window-averaged over the ξ-coordinate over six wave lengths as: where F is the averaged field, N y = 240, and N t = 50. Time averaging is preformed over interval 250 < t k ≤ 300, so that the phase of the surface wave at consecutive steps, t k , t k+1 , changes by 2π and thus the shape of the upper boundary remains the same for all considered time moments, cf. Equation (10). Further the mean vertical profile is obtained by additional averaging of <F> along the ξ-coordinate as: where N x = 360. Phase-averaged fields of the bubbles concentration (void fraction) and its fluxes are also evaluated as: and the respective vertical mean profiles, [CV i n ], are obtained as in Equation (27). Below the term DNS (i.e., Direct Numerical Simulation) is used, although the surface wave dynamics is prescribed and assumed to be stationary, and the bubbles are consid-Algorithms 2022, 15, 110 7 of 14 ered as non-deformable and spherical. However, these assumptions are justified for the considered case of non-breaking waves (sufficiently small wave-slope ka) and bubble sizes (d less than 1 mm). On the other hand, the primitive, full 3D Navier-Stokes equations, Equations (1), (4) and (5), are integrated without employing any closure assumptions.

Results and Discussion
DNS were performed for the surface wave length λ = 15 cm and wave slope ka = 0.1 and 0.2 (amplitude a ≈ 0.25 cm and 0.5 cm). For the chosen length, the wave celerity computed from the linear dispersion relation for surface gravity waves [1] equals c ≈ 49 cm/s, and the Reynolds number, Equation (2), Re ≈ 73019. Both wind-driven case (with q = 0.05c and finite surface drift velocity, U d , Equation (21)) and zero surface-stress case, U d = 0, were considered. Bubble diameter, d, was varied from 200 to 400 microns, the total number of bubbles in each simulation was maintained constant (up to N d = 10 6 for the smallest, d = 200 µm, bubble-size cases) corresponding to a mean (reference) void fraction (concentration, C 0 ) of about 5 × 10 −5 .  [17]. The figure also shows that rising bubbles create vertical streaks of vorticity in their wakes. Figure 3 presents distributions of phase-averaged water velocity components, U x and U z , obtained using DNS with wave slope ka = 0.1 and 0.2 (left and right panels), both with and without bubbles. The components are compared against an analytical, potential-flow, solution for the velocity field in the surface deep-water gravity wave (in dashed line) [1]. The figure shows that the computed fields agree well with the analytical solution except a near surface layer where the wave-induced turbulence modifies the water velocity. The influence of bubbles on the phase-averaged surface-wave flow field is found to be insignificant (in Figure 3, the full black and blue lines practically coincide).

Carrier Flow Modification
Mean profiles of the fluctuations of x,y,z water-velocity components, [U i ], evaluated as a root mean square deviation, and obtained employing DNS for different bubble diameter, d = 200, 300 and 400 µm, with the same reference void fraction (C 0 ≈ 5 × 10 −5 ), both with and without wind-induced surface drift, are compared in Figure 4 with the corresponding bubble-free flow cases. The comparison shows that bubbles enhance wave-induced turbulence, and, in both cases (with and without surface wind drift), mostly z-component of water velocity is affected (Figure 4c,f). This turbulence enhancement is caused by the presence of wakes induced by rising bubbles (Figure 2f). Figure 4f also indicates that the influence of a wind-driven surface stress synchronizes the mean velocity fluctuations for different bubble sizes. The explanation of this observation requires further research and is to be reported elsewhere.

Bubbles Dispersion
The trajectories of individual bubbles of different sizes (with diameters d = 200 µm and 400 µm) obtained using DNS without and with wind-induced surface drift, are presented in Figures 5 and 6, respectively. The figures also show the depth-dependence of the forces imposed on the bubbles by the surrounding water (drag, fluid acceleration, and lift) along the trajectories.   05 (b,e). Panels (a-c) are for the no-bubbles case, and panels (d-f) are for DNS with bubbles with diameter d = 400 µm. In both cases, surface drift velocity, Equation (21), is added at the waved boundary accounting for the wind stress. Wave slope ka = 0.2. both with and without bubbles. The components are compared against an analytical, potential-flow, solution for the velocity field in the surface deep-water gravity wave (in dashed line) [1]. The figure shows that the computed fields agree well with the analytical solution except a near surface layer where the wave-induced turbulence modifies the water velocity. The influence of bubbles on the phase-averaged surface-wave flow field is found to be insignificant (in Figure 3, the full black and blue lines practically coincide).
and obtained employing DNS for different bubble diameter, d = 200, 300 and 400 µ m, with the same reference void fraction (C0 ≈ 5 × 10 −5 ), both with and without wind-induced surface drift, are compared in Figure 4 with the corresponding bubble-free flow cases. The comparison shows that bubbles enhance wave-induced turbulence, and, in both cases (with and without surface wind drift), mostly z-component of water velocity is affected (Figure 4c,f). This turbulence enhancement is caused by the presence of wakes induced by rising bubbles (Figure 2f). Figure 4f also indicates that the influence of a wind-driven surface stress synchronizes the mean velocity fluctuations for different bubble sizes. The explanation of this observation requires further research and is to be reported elsewhere.

Bubbles Dispersion
The trajectories of individual bubbles of different sizes (with diameters d = 200 µm and 400 µ m) obtained using DNS without and with wind-induced surface drift, are presented in Figures 5 and 6, respectively. The figures also show the depth-dependence of the forces imposed on the bubbles by the surrounding water (drag, fluid acceleration, and lift) along the trajectories. Figures 5 and 6 show that the influence of the surface-wave-induced motions on the bubbles dynamics increases as they rise toward the water surface: the trajectories, being    Figures 5 and 6 show that the influence of the surface-wave-induced motions on the bubbles dynamics increases as they rise toward the water surface: the trajectories, being almost straight-vertical at z/λ < 0.5, oscillate with an increasing amplitude as the depth decreases and start spiraling sufficiently close to the surface (Figures 5a,e and 6a,e). The bubbles experience a mean drift in the direction of the surface wave propagation. This drift is similar to the classical Stokes drift of the Lagrangian fluid particles in a surface wave and caused by the decreasing dependence of the wave motion amplitude on depth [29]. As expected, the drift becomes more pronounced in the presence of the wind-induced surface drift ( Figure 6). In the considered case, however, this drift is also modified by the bubbles' ascent due to buoyancy, so that larger bubbles (d = 400 µm) are less subject to this mean Lagrangian drift since their terminal rise velocity is almost 4 times larger as compared to 200 µmbubbles (cf. Figure 7 below). As a result, their dynamics are less affected by the oscillatory wave-field motion. In this aspect, there is some analogy between the observed reduction of the drift due to bubble rising and the "crossing trajectories" effect governing the dispersion of inertial particles by turbulence, where particles' settling reduces their dispersion by ambient turbulence [30]. The analysis of forces governing bubbles dynamics reveals that the drag force, Fd, and the fluid-acceleration force, Fa, are both of the same order, and much larger as compared to the lift force, FL. The forces, Fa and Fd, increase monotonically as the bubbles rise toward the water surface, whereas the lift force becomes non-zero only in close vicinity of the surface, where the fluid motion is not strictly irrotational due to the presence of wave-induced turbulence [17]. The profiles of the mean bubble concentration (or void which is obtained from the bubble equation of motion, Equation (5), rewritten for a bubble rizing with constant (terminal) velocity, Vt, in a quiescent water. Equation is solved by the Newton's method [26].  Figure 7 shows that in our numerical experiment, the bubble concentration decreases exponentially with depth, as observed in natural oceanic conditions [12]. Thus, the injection mechanism employed in our DNS adequately reproduces the required void-fraction vertical distribution.
The parameterization for the horizontal concentration flux, [CVx], in the absence of the wind-induced surface drift (Ud = 0, (b), Equation (31)) is obtained using the Stokes-drift solution [1] adapted here for a somewhat smaller surface drift velocity (Vsd ≈ 0.02) as compared to the classical surface Stokes-drift velocity (c(ka) 2 ≈ 0.04 for the considered surface-wave slope ka = 0.2). Note that a similar effect of a reduced Stokes drift velocity was observed in DNS of a free-propagating, non-breaking surface waves [31]. In The analysis of forces governing bubbles dynamics reveals that the drag force, F d , and the fluid-acceleration force, F a , are both of the same order, and much larger as compared to the lift force, F L . The forces, F a and F d , increase monotonically as the bubbles rise toward the water surface, whereas the lift force becomes non-zero only in close vicinity of the surface, where the fluid motion is not strictly irrotational due to the presence of wave-induced turbulence [17]. The profiles of the mean bubble concentration (or void fraction), [C], and x and z components of the concentration fluxes, [CV x ] and [CV z ], were evaluated in DNS according to Equations (28) and (29) for different bubble size, and both with and without wind stress at the water surface (Figure 7). Figure 7 also compares numerical solution for the fluxes with their parameterizations in the form: where V sd is the bubble drift velocity at the water surface (V sd ≈ 0.02 and V sd ≈ 0.03 for the zero (Figure 7b) and non-zero (Figure 7e) wind-induced-drift cases, respectively); V t is a terminal velocity determined by numerical solution of the following equation: which is obtained from the bubble equation of motion, Equation (5), rewritten for a bubble rizing with constant (terminal) velocity, V t , in a quiescent water. Equation is solved by the Newton's method [26]. Figure 7 shows that in our numerical experiment, the bubble concentration decreases exponentially with depth, as observed in natural oceanic conditions [12]. Thus, the injection mechanism employed in our DNS adequately reproduces the required void-fraction vertical distribution.
The parameterization for the horizontal concentration flux, [CV x ], in the absence of the wind-induced surface drift (U d = 0, (b), Equation (31)) is obtained using the Stokes-drift solution [1] adapted here for a somewhat smaller surface drift velocity (V sd ≈ 0.02) as compared to the classical surface Stokes-drift velocity (c(ka) 2 ≈ 0.04 for the considered surface-wave slope ka = 0.2). Note that a similar effect of a reduced Stokes drift velocity was observed in DNS of a free-propagating, non-breaking surface waves [31]. In the presence of the wind-induced surface drift (U d ≈ 0.04), the same parameterization is used but with V sd ≈ 0.03. Both parameterizations agree well with DNS results.
The vertical component of the void-fraction flux, [CV z ], is well predicted by Equation (32) in both (zero and non-zero wind induced drift) cases. That means that neither waveinduced irrotational motions nor turbulence affect on average the bubble rising rate, so that the void fraction vertical flux is analogous to that in a quiescent water.

Conclusions
We have presented an algorithm for numerical modeling of microbubble dispersion in the near-surface water layer of the upper ocean, under the action of non-breaking, progressive surface waves. The algorithm is based on a Eulerian-Lagrangian approach where full, 3D Navier-Stokes equations for the carrier flow induced by a waved water surface are solved in a Eulerian frame, and the trajectories of individual bubbles are simultaneously tracked in a Lagrangian frame, taking into account the impact of the bubbles on the carrier flow and the impact of a wind-induced surface drift. The bubbles diameters are considered in the range from 200 to 400 microns (thus, micro-bubbles), and the effects related to the bubbles deformation and dissolution in water are neglected. The wave shape is prescribed and assumed to be unaffected by either bubbles or induced turbulent motions.
The simulations results show that bubbles are capable of enhancing the carrier-flow turbulence, as compared to the bubble-free flow, and that the vertical water velocity fluctuations are mostly augmented, and increasingly so by larger bubbles. The results also show that the bubbles dynamics are governed by buoyancy, the surrounding fluid acceleration force, and the drag force whereas the impact of the lift force on the bubble dynamics remains negligible. On the basis of the simulation results, parameterizations for the void fraction fluxes have been obtained. The results show that the vertical component of the void-fraction flux remains unaffected by either the wave motion or wave-induced turbulence as compared to that in a quiescent water. The horizontal void-fraction flux is produced by a mean drift of the bubbles in the direction of the surface wave propagation and can be regarded as analogous to the Stokes drift of Lagrangian (non-inertial) particles in a 2D surface wave modified by the bubbles' ascent. The developed algorithm and parameterizations for void-fraction fluxes can be used for prediction of microbubble dispersion in the ocean upper layer and further employed in large-scale prognostic models.
As a next step in the development of the algorithm, a more complicated problem is to be considered where the wave motion is fully resolved (i.e., not prescribed) thus allowing modification of the surface wave by either induced turbulence or bubbles. This however remains a subject for future research.