## 3 Numerical Modeling  
### Solver Comparison  
Simulating a two-phase flow by solving the Navier-Stokes equations presents several challenges, including the modeling of air entrainment, coupling equations at the air-water interface, tracking the fluids interface, and addressing turbulence interactions between the phases [27]. The present study used two types of OpenFAOM solvers: an interfacial model, interFoam, and an Euler-Euler Model, twoPhaseEulerFoam [28]. InterFoam is a single-phase Euler model that solves the Navier-Stokes equations for two incompressible, isothermal, and immiscible fluids using the VoF method to track the fluids 1 and 2 interfaces, respectively water and air in this study [18; 29; 30; 31]. The fluid density, \(\rho\), in each computational cell is calculated as \(\rho=\alpha\rho_{1}+(1-\alpha)\rho_{2}\) where \(\alpha\) represents the fluid 1 volume fraction: \(\alpha=1\) if the cell is entirely filled with fluid 1, \(\alpha=0\) if the cell is filled with fluid 2. The governing mass and momentum equations are provided in Equations (1) and (2), respectively. \[\
abla\cdot\stackrel{{\rightarrow}}{{\
abla}}=0 \tag{1}\]  
\[\frac{\partial\rho\stackrel{{\rightarrow}}{{V}}}{\partial t}+ \
abla\cdot\left(\rho\stackrel{{\rightarrow}}{{V}}\stackrel{{ \rightarrow}}{{V}}\right)=-\
abla\mathrm{p}_{-\mathrm{rgh}}- \stackrel{{\rightarrow}}{{g}}\stackrel{{ \rightarrow}}{{X}}\cdot\
abla\rho+\
abla\cdot\left[2(\mu+\mu_{t})S \right]+\stackrel{{\rightarrow}}{{F}}_{\sigma} \tag{2}\]  
where: \(\stackrel{{\rightarrow}}{{V}}\), denotes the velocity vector; \(t\), time; \(\mathrm{p}_{-\mathrm{rgh}}=\mathrm{p}-\rho\mathrm{g}\mathrm{h}\), alternative pressure term; \(\mathrm{h}\), elevation; \(\stackrel{{\rightarrow}}{{g}}\), gravitational acceleration vector; \(\stackrel{{\rightarrow}}{{X}}\), position vector; \(\mu\), molecular dynamic viscosity; \(\mu_{t}\), turbulent dynamic viscosity; \(S\), strain rate tensor, and \(\stackrel{{\rightarrow}}{{F}}_{\sigma}\), surface tension. To prevent numerical diffusion and maintain a sharp interface, interFoam incorporates an artificial compression term in the VoF equation, involving a compression coefficient that can have values, typically, between 0 (no compression, i.e., interface undergoes significant numerical diffusion, becoming more stable) and 4 (maximum compression, sharp interface,risk of numerical instability. This term acts near the interface, counteracting numerical smearing by advecting the phase fraction along the interface normal. The present study uses a compression coefficient of 1, ensuring a good balance between interface sharpness and numerical stability. The twoPhaseEulerFoam solver, an Euler-Euler type solver, was also employed in this study. This solver is suitable for solving two compressible, inter-penetrating phases. The continuity and momentum equations associated with this solver are presented in Equations (3) and (4) [18; 32; 33; 34]. \[\frac{\delta\alpha_{\text{d}}}{\delta t}+\
abla\cdot\left(\alpha_{\text{d}} \alpha_{\text{c}}\overset{\rightarrow}{V}_{r}\right)+\
abla\cdot\left(\alpha_{ \text{d}}\overset{\rightarrow}{V}\right)=\alpha_{\text{d}}\
abla\cdot\overset{ \rightarrow}{V}+\alpha_{\text{c}}\
abla\cdot\overset{\rightarrow}{V}_{d}- \alpha_{\text{d}}\
abla\cdot\overset{\rightarrow}{V}_{c} \tag{3}\]  
\[ \begin{split}\frac{\delta\alpha_{\text{op}}\rho_{\text{op}}}{ \delta t}+\
abla\cdot&\left(\alpha_{\text{op},\text{f}}\rho_{ \text{op},\text{f}}\overset{\rightarrow}{V}_{i}\right)-\epsilon_{\text{cont}, \text{op}}\overset{\rightarrow}{V}+\left(\alpha_{\text{op}}\rho_{\text{op}} +C_{\text{om}}\right)\frac{\delta\overset{\rightarrow}{V}_{\text{om}}}{ \delta t}=\\ &=\
abla\cdot R^{eff}-C_{d}\overset{\rightarrow}{V}-C_{\text{om}} \left(\frac{\overset{\rightarrow}{DV}_{\text{op}}}{Dt}-\frac{\overset{ \rightarrow}{DV}_{\text{f}}}{Dt}\right)\end{split}