## 3 Numerical Modeling
### Numerical Domain, Boundary Conditions, and Numerical Schemes  
A numerical domain was defined to replicate the flow conditions of the experimental facility, as shown in Figure 3. The domain dimensions are 4.00 m \(\times\) 2.65 m \(\times\) 1.50 m, matching those of the experimental setup. Similar to the experimental facility, a 72 mm circular inlet is located 1.00 m above the bottom of the domain. The only difference between the two setups is how the outlet boundary condition is imposed. Instead of reproducing the flap gates considered in the experimental facility, in the numerical model, the outlet boundary condition was set up as detailed below to allow the simulation of different pool depths without a need for mesh adjustments that the simulation of the flap gates would require.  
The boundary conditions were closely defined to match those of the experimental facility and were consistently applied to allow direct comparison with the results from each solver.  
The top part of the domain was kept open to atmospheric pressure for all variables, while the walls and bottom were treated as natural surfaces with a fixed velocity of zero. Wall functions were applied in the turbulent kinetic energy and dissipation fields, and a roughness wall function was used to simulate wall roughness.  
A fixed velocity corresponding to the specified inlet velocity was applied to the inlet. Since only water entered the domain, the water volume fraction was set to one. Turbulent kinetic energy was defined based on turbulence intensity, which was determined through a sensitivity analysis. The selected range of values was derived from experimental measurements of turbulence intensity [4; 10; 35].  
Two gates represented the model outlets. To model this scenario, a wave absorption boundary condition was implemented to maintain a constant water level in the model of 0.8 m [36].  
The Reynolds-Averaged Navier-Stokes (RANS) k-Omega SST model was used in both solvers to characterize flow turbulence. This two-equation model assumes that Reynolds stresses are proportional to the mean velocity gradients [37; 38]. [PERSON] et al. [19; 21] concluded that the k-Omega SST model provided the best agreement between experimental and numerical results in an impinging jet. Similarly, [PERSON] et al. [17; 18] observed a better agreement between numerical and experimental results when simulating a water jet and a spillway aerator with the k-Omega SST model. The k-Epsilon turbulence model was not used, as it produced unstable results in the preliminary simulations when coupled with the turbulent dispersion model in the twoPhaseEulerFoam solver.  
The twoPhaseEulerFoam solver can couple a turbulent dispersion model with the general flow model. Three different models were evaluated: [PERSON] [39], [PERSON] [40], and [PERSON] [41]. Preliminary tests were conducted to determine the most suitable model. The [PERSON] model was discarded due to numerical instability, while the [PERSON] and [PERSON] models produced similar flow conditions. The [PERSON] model was ultimately selected for its lower computational effort compared to the [PERSON] model.  
For numerical schemes, an Euler implicit time scheme was employed for time discretization, and a Gauss linear upwind scheme was applied for the divergence of velocity, k, and omega. Flow characteristics from the numerical simulations were obtained after an initial simulation period to ensure time independence of the water volume in the domain, and the values of the total kinetic and turbulent kinetic energies were also time independent.  
The simulations were performed on an HPC cluster featuring Intel Xeon E5-2680 (2.70 GHz) processors, each with 16 cores. The number of cells per processor remained constant for the first three meshes. However, for the most refined mesh, the number of cells assigned to each processor increased due to limitations on the number of processors. This constraint led to significantly longer computation times.  
Figure 3: Numerical domain.