## 3 Materials and Methods
### Breakup Rate Modeling
#### 3.3.1 Modeling Approaches  
One of the most common breakup models is the one proposed by [PERSON] and [PERSON] [9]. In this model, the interactions between bubbles and eddies are modeled using an analogy with the molecular gas collisions in an ideal gas. Hence, the breakup rate is computed by multiplying the collision frequency with breakup efficiency. The breakup efficiency is the probability that turbulent eddies have enough energy to break up a bubble. The collision frequency is defined as the volume swept in unit time by the approaching bubble and turbulent eddy; it is computed as the product of the cross-sectional area, the relative bubble-eddy velocity and the density number of the eddies. This last parameter is the weakest point of the model since it requires some assumptions and the model performance depends highly on its value. [PERSON] and [PERSON] [9] expressed the number of eddies as a function of wave number using the equation proposed by [PERSON] and [PERSON] [21], which reads as follows:  
\[\frac{dn_{e}(k_{b})}{dk_{b}}=0.1{k_{b}}^{2} \tag{23}\]  
where \(n_{e}(k)\) is the number of eddies of wave number \(k\):  
\[k_{b}=2/d_{eq} \tag{24}\]  
Equation (24) provides an infinite number of eddies as the wave number \(k\) goes to infinity (i.e., for very small-sized eddies, Equation (25)). For this reason, [PERSON] and [PERSON] [9] stated that this expression should not be integrated over the eddies that lie in the viscous dissipation range. They argued that the eddies with a characteristic size smaller than the 0.2 \(d_{eq}\) are unlikely to contribute significantly to the overall breakup rate. Indeed, such small eddies have only 0.5% of the kinetic energy associated with an eddy, having a size equivalent to the bubble one. [PERSON] et al. (2018) stated that this model is very sensitive to both the upper and lower limit of integration and that they both must always be specified ([PERSON] and [PERSON] (2018) did not define the lower limit of integration of Equation (25)). For this reason, due to the uncertainties in the evaluation of the limits of integration, we implemented the breakup source term as the phenomenological model of [PERSON] et al. (2018). In particular, in this approach, the breakup rate is computed as a function of the breakup frequency only:  
\[S_{b,i}=\Omega(d_{eq,i}) \tag{25}\]  
In this approach, the breakup rate is modeled taking into account the turbulent fluctuations of the carrier phase. This model is based on a one-way coupling between the two phases: the motion of bubbles does not affect the velocity field of the carrier phase. This is, of course, a main drawback of the present paper as the motion of the liquid is caused by the moving of bubbles due to buoyancy, mainly through the drag force between the two phases. Also, the bubble-induced turbulence contributed much to the liquid turbulence. This assumption is valid for small void fractions (low \(U_{G}\)). Furthermore, it is assumed that the bubble sizes are in the inertial subrange, meaning that there is an equilibrium between the energy transferred from the larger turbulent eddies toward the smallest dissipative eddies.