## 3 Materials and Methods
### Breakup Rate Modeling
#### 3.3.2 Breakup Frequency Modeling  
The model relies on deformation energy that is provided by the turbulent stresses of the carrier phase. This energy stretches and deforms the bubble and opposes the surface restoring pressure. The minimum energy to deform a bubble, \(E_{S}\), with a diameter \(d_{eq}\) depends on the superficial tension \(\sigma\)([PERSON] et al., 2018):  
\[E_{S}(d_{eq})=\pi\sigma d_{eq}^{2} \tag{26}\]  
Hence, the surface restoring pressure reads as ([PERSON] et al., 2018):  
\[\tau_{S}(d_{eq})=\frac{6E_{S}}{\pi d_{eq}^{3}}=6\frac{\sigma}{d_{eq}} \tag{27}\]  
The average deformation force per unit surface caused by the turbulent stresses between two points separated a distance \(d_{eq}\) is computed as follows ([PERSON] et al., 2018):  
\[\tau_{t}(d_{eq})=\frac{1}{2}\overline{\rho\Delta u^{2}}(d_{eq}) \tag{28}\]  
The condition \(\tau_{s}=\tau_{t}\) identifies a critical bubble diameter above which bubbles break up. On the other hand, a bubble with smaller \(d_{eq}\) than the critical value has enough surface energy to overcome the deformation energy and will break up. According to [PERSON]'s theory of isotropic turbulence, the mean value of the velocity fluctuations between two points separated a distance \(d_{eq}\) is estimated by ([PERSON] et al., 2018):  
\[\overline{\Delta u}\left(d_{eq}\right)=\overline{\left|u(x+d_{eq},t)-u(x,t) \right|^{2}}=\beta(\varepsilon d_{eq})^{2/3} \tag{29}\]  
The bubble maximum stable diameter is computed by equaling the surface energy to the deformation energy between two points separated by a distance equal to \(d_{min}\)([PERSON] et al., 2018):  
\[\frac{1}{2}\rho\beta\varepsilon^{2/3}d_{\min}^{2/3}=6\frac{\sigma}{d_{eq}} \tag{30}\]  
Hence, \(d_{min}\) reads as follows:  
\[d_{\min}=\left(\frac{12\sigma}{\beta\rho d_{eq}}\right)^{3/2} \tag{31}\]A larger difference between the restoring and the deformation pressures would break up the bubble at a certain time. On the other hand, the breakup frequency must tend to a zero limit when these two pressures are equal. Therefore, following the proposal of [PERSON] et al. (2019), the bubble breakup frequency is estimated as follows:  
\[\begin{split}\Omega(d_{eq,i})&=K\frac{\sqrt{\frac{ \Delta t^{2}(d_{eq})-\frac{12r}{\mu d_{eq}}}{d_{eq}}}}{d_{eq}}=\\ &=K\frac{\sqrt{\beta\left(\varepsilon_{dd_{eq}}\right)^{2/3}- \frac{12r}{\mu d_{eq}}}}{d_{eq}}\end{split} \tag{32}\]  
where \(\varepsilon\) is computed using Equation (15), \(\beta\) is a constant equal to 8.2 [PERSON] et al. (2019), and \(K\) is a constant equal to 0.25 in the model of [PERSON] et al. (2019) (obtained by using experimental data). In this paper, \(K\) is implemented as a function of \(U_{G}\) (determined as the best fit of the model results and the experimental data):  
\[K=0.85-18.92U_{G}\quad U_{G}\;in\;[m/s] \tag{33}\]  
Indeed, \(K\) is related to the boundary conditions used. The original model of [PERSON] et al. (2019) was proposed with uniform BSDs at the inlet. Contrarily, we used the BSDs sampled at the sparger. We observed that the original value \(K=0.25\) proposed by [PERSON] et al. (2019) is suitable only for uniform BSDs at the inlet. The reader may refer to the sensitivity analysis proposed in Section 4 for a discussion concerning this issue.