Do the Volume-of-Fluid and the Two-Phase Euler Compete for Modeling a Spillway Aerator?

: Spillway design is key to the effective and safe operation of dams. Typically, the ﬂow is characterized by high velocity, high levels of turbulence, and aeration. In the last two decades, advances in computational ﬂuid dynamics (CFD) made available several numerical tools to aid hydraulic structures engineers. The most frequent approach is to solve the Reynolds-averaged Navier–Stokes equations using an Euler type model combined with the volume-of-ﬂuid ( VoF ) method. Regardless of a few applications, the complete two-phase Euler is still considered to demand exorbitant computational resources. An assessment is performed in a spillway offset aerator, comparing the two-phase volume-of-ﬂuid ( TPVoF ) with the complete two-phase Euler ( CTPE ). Both models are included in the OpenFOAM ® toolbox. As expected, the TPVoF results depend highly on the mesh, not showing convergence in the maximum chute bottom pressure and the lower-nappe aeration, tending to null aeration as resolution increases. The CTPE combined with the k – ω SST Sato turbulence model exhibits the most accurate results and mesh convergence in the lower-nappe aeration. Surprisingly, intermediate mesh resolutions are sufﬁcient to surpass the TPVoF performance with reasonable calculation efforts. Moreover, compressibility, ﬂow bulking, and several entrained air effects in the ﬂow are comprehended. Despite not reproducing all aspects of the ﬂow with acceptable accuracy, the complete two-phase Euler demonstrated an efﬁcient cost-beneﬁt performance and high value in spillway aerated ﬂows. Nonetheless, further developments are expected to enhance the efﬁciency and stability of this model.


Introduction
Spillways are essential for dam safety, controlling the reservoir water level and discharge, and preventing dam overtopping, one of the leading causes of structural failure and rupture. Such a crucial hydraulic structure requires careful design validation. Typically, a spillway is constituted by: intake, weir, chute, energy dissipation structure, and river restitution. The water flow from the reservoir approaches the intake in the subcritical regime with low velocity. Next, the flow is guided into the weir, the control section, where the critical regime is attained. Downstream, the chute flow is supercritical, with high velocities (frequently larger than 20 m/s), high levels of turbulence, and significant air entrainment and transport. Frequently, the geometry imposes complex flow patterns, such as crosswaves, significant free-surface variations, and lateral-mixing. At the end of the spillway, a structure dissipates energy, allowing proper restitution to the river stream. Moreover, surrounding and entrained air plays a crucial role in the safe operation of spillways [1] and may significantly alter the flow characteristics. Aeration must be considered due to the effects of flow bulking, drag reduction, prevention or mitigation of cavitation damage, reduction of the breakup length of water jets discharging into atmosphere, interaction with the turbulence field, re-oxygenation, and transference of atmospheric gases that have a vital function on stream ecology [2,3]. Besides free-surface aeration-the natural air entrainment air, followed by the impact zone, i.e., where the water jet re-joins the bottom. The high turbulence levels cause lower-nappe aeration at the air-water interface of the cavity zone.
An unprecedented 3D numerical study is conducted due to the mesh size, the assessment of multiple flow characteristics [44], the physical model with detailed data and high Weber and Reynolds numbers that reduce scale effects [6], and the enormous computational efforts. The paper is organized as follows. After this introduction, the methodology is described, including the laboratory data, the mathematical models, and the numerical setup. Next, the results are analyzed. Finally, the discussion and conclusions are presented.

Methodology
The complete two-phase Euler and the volume-of-fluid models are evaluated against the physical modeling of an offset chute aerator by Bai et al. [11] (Figure 1).
The main flow characteristic to be evaluated is the lower-nappe aeration induced by the offset aerator. This phenomenon is quantified by the flow-rate of lower-nappe entrained air. The bubble diffusion is also assessed through selected water fraction profiles downstream of the impact zone.
Additionally, two characteristics related to the pressure and flow geometry downstream of the aerator are analyzed: the value and location of the maximum pressure increment at the spillway bottom and the length of the cavity zone.
A mesh dependence analysis is performed for six different resolutions. The two most popular Reynolds Averaged Navier-Stokes (RANS) turbulence models are tested in each model: k-ω SST and k-. Hence, the models' performance is analyzed for a total of 24 combinations (2 models × 6 resolutions × 2 turbulence models).

Laboratory Data
Bai et al. [11] conducted unmatched physical modeling of a large spillway offset aerator with high velocity (4 to 9 m s −1 ) and high inflow water depth (0.15 m). Hence, presenting high values of the Reynolds (5.5 × 10 5 < Re < 1.2 × 10 6 ) and Weber numbers (180 < We 0.5 < 405) which comply with Pfister [6] criteria to mitigate aeration related scale effects: Re > 2.2 × 10 5 and We 0.5 > 140. Moreover, the extent number of tests and a complete set of measured flow properties (i.e., pressure at the bottom, velocity profiles, air-concentration, turbulence, and bubble sizes) provide relevant data to validate numerical models.
The rectangular channel is 5 m long and 0.25 m wide (Figure 1a). The aerator offset height (h s ) varies from 0.025 m to 0.045 m. The upstream emergence angle (θ 0 ) inclination is variable from 0°to 14.1°and the channel bottom angle (θ b ) ranges from 5.7°to 14.1°. The channel bed and side-walls are made of smooth polyethylene with a roughness height of 1 × 10 −5 m.
A case with inlet velocity (V 0 ) of 9 m s −1 , water depth (h 0 ) of 0.15 m, equal emergence and bottom angle of 14.1°and offset height of 0.045 m are selected to mitigate the scale effects. Froude number (Fr) is 7.4, Re = 1.2 × 10 6 and We 0.5 = 405 .

Mathematical Models
The CFD solvers and the turbulence models used are included in the OpenFOAM ® toolbox version 2012 [46]. The volume-of-fluid solver is named interFoam. The complete two-phase Euler solver is named twoPhaseEulerFoam.

Two-Phase Volume-of-Fluid
The two-phase volume-of-fluid solves the Reynolds-averaged Navier-Stokes equations for two incompressible, isothermal, and immiscible fluids. The interface capture is based on a volume-of-fluid (VoF) method that incorporates an interfacial compression flux term [47,48]. Hereinafter, this model is mentioned as TPVoF.
The mass Equation (1) and a single-set of momentum conservation equations (2) are solved. Hence, the two fluids-water and air-share the same velocity field.
where V is the RANS velocity vector, ρ is the density, t is the time, p * is the pseudo-dynamic pressure, g is the gravitational acceleration vector, X is the position vector, µ is the molecular dynamic viscosity, µ t is the eddy viscosity coefficient (i.e., turbulent dynamic viscosity) and S is the strain rate tensor. F σ is the surface tension force term for the momentum Equation (3). σ is the surface tension, κ is the curvature of the interface and α is the water volume fraction. Any VoF phase property Φ (e.g., ρ, µ) is a volume-average of the intrinsic fluid property of water (subscript w) and air (subscript a) (4).
Additionally, according to the VoF method implementation in OpenFOAM ® , the fluids volume fraction in each cell is defined by a scalar function (α) ranging from 0 to 1 that allows the interface tracking: α = 1 is a water cell and α = 0 is an air cell. Other α A compression velocity ( V ic ) (6) proportional to the local velocity field magnitude is applied perpendicular to the interface. The interface compression is triggered by the C α coefficient that usually ranges from 0 to 1.
The solver uses a segregated solution approach. Each time step starts an update of the interface, followed by a prediction of the velocity. A Pressure Implicit with Splitting Operators (PISO) type algorithm corrects velocity and implicitly the pressure, advancing both in time. Finally, the turbulence is calculated.

Complete Two-Phase Euler
The complete two-phase Euler solves the Reynolds-averaged Navier-Stokes equations for two interpenetrating and compressible fluid phases-a continuous (c) and a dispersed (d)-including heat transfer [49,50]. Hereinafter, this model is mentioned as CTPE.
The continuity equation is solved for the dispersed phase (7), where V r is the phases relative RANS velocity vector, T D,d is the dispersed phase turbulent dispersion term and T D,c is the continuous phase turbulent dispersion term. The momentum Equation (9)-written for a generic phase iis solved for each phase but not directly. Instead, following Issa [51] methodology, the face flux is initially predicted and afterward corrected by the pressure, which is shared between the two phases, and solved iteratively.
where φ is the face-flux vector, ε cont is the continuity error, C d is the drag coeficient and C vm is the virtual mass coeficient. R e f f is the stress rate tensor (10), ν e f f = ν + ν t is the effective viscosity and k is the turbulent kinetic energy. The energy conservation equation is solved for each phase (see Equation (11)) in the form of internal energy or enthalpy. Therefore variable he is a hybrid variable that represents one or the other.
where K is the kinetic energy, T is the fluid temperature, α e f f is the effective thermal diffusivity, k h is the convective heat transfer coefficient, k e f f h is the effective volumetric convective heat transfer coefficient and C pv is the specific heat capacity. The IF superscript refers to an interfacial property.
Several closures for heat transfer, drag, lift, and turbulent dispersion are available. Only a single-size air bubble can be adopted. The CTPE allows considering both fluids as continuous or dispersed, depending on the respective volume-fraction value. Using a blended interfacial model, for each cell, the solver identifies one of the following scenarios: phase 1 is dispersed and phase 2 is continuous, phase 2 is dispersed and phase 1 is continuous, or it is a cell with no obvious dispersed phase. Several blending options are available: constant, linear, and hyperbolic.
The turbulent dispersion force is defined by Burns et al. [52] and Otromke [53], as follows: where C D is the non-dimensional drag coefficient and σ α is the turbulent Prandtl number for interfacial area density. The turbulent dispersion force is a function of the blending interfacial model chosen and its settings. The solver calculates the turbulent dispersion at each cell for one of the previous three blending model scenarios. Hence, the turbulence dispersion action depends totally on the blending model's chosen parameters, especially on the maximum volume-fraction value to a phase be considered as dispersed. For example, the user may define this value as 0.6, i.e., if the volume-fraction values of a phase range from 0 to 0.6, that phase is considered dispersed. Higher values, consider the phase continuous or with no obvious dispersed phase, depending on the model applied. Each time step starts by solving the phase continuity, followed by discretization and linearization of the momentum equations that predict the flux. Next, energy conservation is solved. At this moment, the pressure is solved, and the flux and velocity are corrected. Finally, the turbulence is calculated.

Turbulence Models
The kand k-ω SST RANS turbulence models are among the most used in hydraulic structures, hence are applied in the numerical study of the spillway bottom aerator. The kmodel is conceived for internal flows and is the most common in RANS simulations. The k-ω SST is used due to its advantage for boundary flows. Both models have two transport equations: one for the turbulent kinetic energy (k) and another for the turbulent dissipation rate ( ) or the turbulent specific dissipation rate (ω), which are used to determine the eddy viscosity. In the volume-of-fluid, both models do not include density explicitly. Therefore, instead of the dynamic form (µ t ), is calculated the kinematic eddy viscosity (15) [54,55]. The default coefficients are employed.
In the complete two-phase Euler, instead of the standard kmodel, is applied the k-mixture which is a specific turbulence model for two-phase gas-liquid systems, based on Behzadi et al. [62] and Lahey [63]. This model solves a single set of equations for the airwater mixture-m subscript-to determine k m (20) and m (21). The mixture variables are calculated through an effective density weighted average of air and water properties (19).
where ρ a,e f f = ρ a + C vm ρ w . ν m = Φ m (ν t,w , ν t,a ) is the mixture turbulent kinematic viscosity, ν t,w and ν t,a are water and air tubulent kinematic viscosity. G m is the k m production rate, V m is the mixture velocity vector, G b is the k m production rate by air-bubble, ρ w is the water density, ρ a is the air density and ρ m is the mixture density (22). C vm is the virtual mass coeficient. C 3 = 1.92. The remaining constants share same values with the kmodel. Furthermore, the k-mixture model is improved following Weller et al. [64]. Thus, a phase fraction limiter is applied to the bubble-generated turbulence (G b ), avoiding spurious turbulence generation where bubbles are not present. In the current work, G b is only activated if α d > 0.3, a standard value.
In the complete two-phase Euler, the standard k-ω SST model is modified to include the bubble-induced turbulent viscosity model of Sato et al. [65]. Thus, a term is added to the turbulent viscosity Equation (18), as follows: where y + is the dimensionless distance from the wall, c b = 0.6 and d b is the characteristic bubble size. Hereinafter, this model is mentioned as k-ω SST Sato.

Numerical Setup
The numerical domain ( Figure 2) is defined by a nozzle with 0.5 m of length and a rectangular section with 0.25 m of width per 0.15 m of height (h 0 ), followed by a rectangular channel with 2 m of length, 0.5 m of height, and the exact width of the nozzle. On each side wall, by the bottom and next to the step, is placed an air-vent with 0.02 × 0.02 m 2 . The nozzle and channel bottom have a downward slope angle of 14.1°.
Rules of thumb in hydraulic structures CFD recommend a minimum of 10 to 20 cells per characteristic hydraulic length (e.g., water depth, pipe radius, etc.). Nevertheless, the dependence of the result on the mesh resolution must be analyzed for every flow and numerical setup. Thus, a set of fully orthogonal and hexahedral 3D meshes with 6 different resolutions (R m ) relative to the nozzle height (h 0 ) is considered ( Table 1). The cell edge length (d m ) ranges from 2.5 mm (R m = 60) to 15 mm (R m = 10). Mesh resolution is limited to 60 cells due to the extremely small time-step needed to verify the Courant-Friedrichs-Lewy condition (∆t < d m /V 0 = 2.5 × 10 −3 /9 = 2.8 × 10 −4 s) that leads to impractical calculation times, even in high-performance computing clusters (HPC). To calculate in parallel at the HPC, the mesh was decomposed into 16 to 256 sub-domains, using the Scotch method.  The boundary conditions are presented next. At the inlet, the velocity (V 0 ) is 9 m s −1 ; the pressure is automatically calculated to assure the flux. The characteristic turbulent mixing length is 0.0105 m (0.07h 0 ). At the top, outlet and air vents, a total pressure condition is defined. At the outlet, the velocity condition is zero gradient for outflow, and the inflow is blocked. Both top and air vents have a binary velocity condition: zero gradient for outflow and exclusive pressure-driven normal air inflow. Turbulence for the previous three boundaries is zero gradient for outflow and for inflow k = 3.1 × 10 −5 m 2 s −2 , = 1.6 × 10 −6 m 2 s −3 and ω = 0.58 s −1 , considering a turbulent mixing length equal to the channel width (0.25 m) and an intermediate turbulent intensity of 0.05. Walls have no-slip tangent velocity and a turbulent wall function for the roughness of 1 × 10 −5 m, as described in the physical modeling. Due to the absence of specification of the inflows' turbulence intensity (TI) in the physical modeling, two values were tested: 0.005 and 0.2. In the complete two-phase Euler, the temperature of the fluids is intended to be 300 K (26.85°C). Hence, the initial temperature of the domain and the inflow temperature are set to 300 K, and the walls have a zero-gradient condition.
An arbitrary bubble diameter of 0.1 mm is adopted after observing the multitude of sizes of air-bubbles and air-packets in Figure 1a). Furthermore, to respect the CTPE's conception, the bubble diameter must be inferior to the cells' size.
The following numerical settings are used: Euler time derivative; gradient is Gauss linear with multi-directional cell limiter 1.0; the divergence scheme is linear-upwind for velocity. For turbulence fields, i.e., k, and ω, the divergence scheme is upwind. Regarding the TPVoF model, interface compression is based on a generic limited scheme for the divergence of α and an interface compression coefficient (C α ) of 0.5. For the CTPE model, among the available closures, the following models were used: Schiller and Naumaan [66] for drag, Ranz and Marshall [67] for heat transfer, Tomiyama et al. [68] for lift and the turbulent dispersion method based on Burns et al. [52] and Otromke [53]. This study applies a hyperbolic blending function that considers the fluid continuous if the phase volume-fraction value is larger than 0.6.
The simulations are transient. The PISO pressure-velocity algorithm is used with standard settings. Hence data is sampled every time step (approximately 2000 to 10,000 Hz, depending on the mesh) during 1 s and averaged after a warm-up period where semi-steady flow conditions are reached. Approximate time-independence is attained for lower-nappe aeration (β), maximum pressure at chute bottom (∆p b max ), open-boundaries air and water flow-rate, and several other domain properties: water and air volume, total kinetic energy, total turbulent kinetic energy and minimum, and maximum pressure. The presented profile plots and bottom data are collected in the central plane along the chute, i.e., the mid-plane of the chute.
Calculations were performed in two HPC clusters composed of Intel Xeon E5-2680 or AMD EPYC 7501 CPUs (Central Process Unit). Simulations ranged from 16 to 256 cores and needed 10 to 50,000 CPU.core hours (2000 CPU core.days). Due to the extremely high calculation time, the two highest resolution meshes (R m ≥ 40) are considered unfeasible for engineering studies. Overall, more than 18,000 CPU core.days, or 50 CPU core.years, were utilized to deliver the presented results.

Results
The main purpose of this study is to assess the lower-nappe aeration (i.e., induced by the offset aerator) in the different model combinations. Hence, the results are presented and analyzed as follows. First, a global flow comparison. Second and most important, the air-water mixture and the lower-nappe aeration. Next, the maximum pressure increment at the spillway bottom, which is a frequent design criterion. Finally, a geometric parameter of the flow: the cavity zone length. Globally, Figure 3 shows similar flow depth (solid blue surface) and velocity (darkblue and red streamlines) in all model combinations. The only exception is the CTPE with k-mixture that, downstream of the impact zone, presents a significant increment of the water-depth and a smeared air-water interface. In all model combinations, the offset aerator separates the flow from the spillway bottom, creating a jet and a cavity zone filled with air immediately downstream of the aerator. At the cavity zone, the turbulence of the air-water lower interface of the jet promotes the entrainment of air into the flow. Thus, to feed this process, air enters continuously through the air-vents located laterally at the aerator.

Air-Water Mixture
Regarding aeration, the primary conditioning to the performance of the models is how they deal with the air-water interface and if the air is entrained into the water flow. Figure 4, shows the water fraction along the spillway, in an xz plane, for all model combinations and four mesh resolutions.
The major standout is the very distinct behavior of the CTPE with k-mixture. The lower-nappe aeration is much more intense, inducing an air-water mixing region at the bottom that occupies more than half of the flow depth. Moreover, it shows a significant air-entrainment at the upper water surface, resulting in a diffuse interface, which starts more upstream as the resolution increases. Therefore, the CTPE with k-mixture models combination is not tested with the highest mesh resolutions: 40 and 60 cells per inlet height. These phenomena can be explained by the concept of the k-mixture turbulence model and are addressed in Section 3.5.
The second standout is the absence of an air-water mixture region at the bottom in the TPVoF combinations. Oppositely, this region is identified in the CTPE. Moreover, in the TPVoF, an air-water interface is present next to the spillway bottom, though becoming sharper as resolution increases. Thus, presenting a clear fluid separation and an airflow next to the bottom. This phenomenon is also noticed in the lower part of the water-fraction vertical profiles analyzed in Section 3.2.

Lower-Nappe Aeration
A spillway aerator is designed to entrain air into the water flow, eliminating cavitation damage next to the bottom and walls [5]. Air-concentrations ranging from 1 to 8% are needed to mitigate cavitation damage [8,69]. Further downstream of the aerator, air bubbles tend to rise to the flow surface, reducing the air-concentration next to the bottom. Therefore, it is frequent to find several aerators along a spillway.
Thus, the performance of a spillway aerator is determined by the amount of air entrained at the cavity zone (i.e., lower-nappe aeration) and in the air-concentration profile next to the bottom along the spillway.
Commonly, lower-nappe aeration is characterized by the ratio (β) between the flowrate of air-entrained at the cavity zone (Q al ) and water flow-rate at the inlet (Q w ), see Equation (24). Q al is measured by the air flow-rate entering the cavity zone through both lateral air-vents ( Figure 2). Figure 5 shows the β coefficient for all model combinations and different mesh resolutions. For low mesh resolutions (R m < 20), all model combinations show the same behavior: β is approximately twice the reference value ( β re f = 0.0475 in [11]), decreasing as the resolution increases.
A significant result is that in the TPVoF lower-nappe aeration tends to a null value, constantly reducing with the increase of the mesh resolution.
The CTPE with k-ω SST Sato presents acceptable results for higher mesh resolutions (R m ≥ 30), converging to a lower-nappe aeration slightly inferior to the reference value. Oppositely, the CTPE with k-mixture reveals an unexpected increase of β for R m ≥ 20, becoming even more abrupt for higher mesh resolutions. This behavior can be explained by the concept of the k-mixture turbulent model and is addressed in Section 3.5.

Aerator: inte
x / L [-] chute bottom pressure / L  The distribution of the air entrained due to lower-nappe aeration is analyzed in three vertical profiles of the water volume fraction (α): at the cavity zone, at the impact zone, and downstream of the impact zone, i.e., 0.2, 0.7, and 1.0 m downstream of the step. Figure 6 presents these profiles for an intermediate mesh resolution (R m = 30), which represents the behavior of the models.
For the three profiles, the TPVoF shows no air-water mixing. α is null at the bottom, and there is an abrupt transition to the complete water flow above, which is exacerbated with the mesh resolution increase. Thus, all the flow-rate of air from lower-nappe aeration is located next to the chute invert. The absence of air-water mixing is justified because the TPVoF model equations do not comprehend turbulent dispersion. Hence, any air-water mixing is due to the numerical diffusion in the interface, which is reduced with higher mesh resolution. To respect the TPVoF model conception, the mesh must be fine enough to reproduce the individual air-bubbles, which is unfeasible in engineering applications of spillways.
Better behavior is found in the CTPE with k-ω SST Sato. The mixed flow region is located in the lower half of the flow, sharing some similarities with the water fraction (α) profiles of Bai et al. [11]. However, at the bottom, in the cavity zone α = 0.4, and in the impact zone α = 0.7, which is considerably inferior to the reference profile, where α is 0.8 and 0.92.
At the impact zone (see Figure 6b), the experimental data shows a mixed flow region with half the thickness of the flow and larger air-concentration in the inner part of the flow compared to the numerical models, except the CTPE with k-mixture. Downstream of the impact zone (i.e., 1.0 m downstream of the step, see Figure 6c)), despite the reference profile shows only a point at z = 0.11 m where α ≈ 1, the CTPE with k-ω SST Sato exhibits a thick layer (0.08 ≤ z ≤ 0.13 m) with α = 1. The presented facts demonstrate the imperfect modeling of air-bubble vertical transport and turbulent dispersion, which may be justified by the single-size of the air-bubbles and the absence of significant roughness and oscillations in the air-water interface of the cavity zone due to the RANS framework, which is observed in Figure 1a).
The CTPE with k-mixture shows a low and non-realistic water concentration (α) vertical profiles, exposing the exacerbated lower-nappe aeration and bubble diffusion. This behavior is associated with the k-mixture formulation and is explained in Section 3.5.

Pressure Increment at the Spillway Bottom
The maximum pressure at the spillway bottom resultant from the jet impact is an important design criterion, especially in terms of materials lifetime. Thus, it is evaluated the pressure increment above the atmospheric pressure of 101325 Pa = 1 atm and compared against the laboratory data of Bai et al. [11]. The profiles of the pressure increment along the bottom (∆p b ) are presented in Figure 7. This analysis is focused on two properties: the value of the maximum pressure increment (∆p b max ) and its location (L m ). (d) Regarding the maximum pressure increment (∆p b max , see Figure 8), both the TPVoF and CTPE models yield an increase of ∆p b max as the mesh resolution increases, which is even more noticed in the k-ω SST type models. Roughly, for the model combinations, the ∆p b max is only near to the reference value ( ∆p b max re f = 6.7 × 10 3 Pa in [11]) for the higher mesh resolutions (R m ≥ 40). Although the CTPE with k-mixture shows simillar ∆p b max for the lower mesh resolutions (R m = {10, 15}), at higher resolutions ∆p b max suffers an abrupt drop. For the different model combinations and the tested mesh resolutions, it is impossible to foresee a convergence of the ∆p b max . Besides the flow trajectory and momentum, the pressure next to the bottom is directly linked to local flow density, i.e., the mixture air-concentration. Thus, the results are coherent with the lower-nappe aeration ( Figure 5) and the water volume-fraction profiles ( Figure 6) presented in Section 3.2.
A sensitivity analysis of the mesh refinement near the walls was performed, evidencing that the results are not significantly affected. Hence, it is purposely not presented because it has a much smaller impact on the results than the remaining analyzed variables. Additionally to the mesh resolution near the bottom, the reproduction of the flow boundary layer, and the solver framework, including the pressure-velocity coupling scheme, may have a significant role in determining the pressure next to the bottom.
A y+ analysis was also performed, evidencing that the y+ condition for both models is never achieved (k-: 30 < y+ < 300, k-ω SST: 5 < y+ < 20). Moreover, an estimate based on the freestream flow (9 m/s) indicates that the height of the first cell next to the chute bottom should be approximately h 0 /1000 to respect the k-y+ condition and even more to respect the k-ω SST y+ condition. Therefore, any significant resolution increase is incompatible with practical engineering, even refining only near the walls. Nevertheless, the authors consider it is important to assess the turbulence models that are undoubtedly the most used with these solvers, in mesh resolutions representative of the commonly applied in research and practical engineering of hydraulic structures. The location of the maximum pressure increment L m is the distance along the chute bottom from the offset step to a point where the pressure increment is maximum (∆p b = ∆p b max ), see Figure 9.
L m is similar for all model combinations, except for the CTPE with k-mixture, increasing in higher mesh resolutions. L m starts to converge in intermediate mesh resolutions R m ≥ 30 to a value approximately 15 to 35% larger than the reference (L m re f = 0.62 m in [11]). The CTPE with k-mixture shows abnormally small values for R m ≥ 20, due to the exacerbated jet diffusion. The TPVoF with kand the CTPE with k-ω SST Sato present the peak of pressure nearer the reference value.

Cavity Zone Length
Following Bai et al. [11], the cavity zone length (L, see Figure 10) is the distance along the chute bottom between the offset step and the most upstream point of the impact zone where ∆p b = 0.1 × ∆p b max . L also characterizes the jet trajectory.
In all model combinations using a k-ω SST type turbulence model, L tends to converge to a value 15 to 20% larger than the reference (L re f = 0.52 m in [11]). Comparing Figures 1a and 4, the air-water interface at the cavity zone is absent of significant roughness and oscillations in the air-water interface of the cavity zone due to the RANS framework. Also, the interface sharpness increases with higher resolutions. Thus, the lower interface of the jet is not diffuse enough, and the first impact point is located downstream than what is observed experimentally.
The TPVoF with kseems to converge to a value nearer to the reference. The CTPE with k-mixture displays minimal values for R m = 30, similarly to what is observed in L m (Figure 9). TPVoF, k-ε TPVoF, k-ω SST

Turbulence
A complete assessment of the turbulence is impossible due to the absence of detailed experimental data of k, , ω, and ν t in [11]. Nevertheless, the four turbulence models are compared.
As presented in Figures 4 and 5, the kand k-mixture present larger lower-nappe aeration than the corresponding k-ω SST and k-ω SST Sato. Analyzing the k of the TPVoF in Figure 11 and the k w of the CTPE in Figure 12, the turbulent kinetic energy (k) is in the same order of magnitude at the air-water interface of the cavity zone. In detail, the TPVoF combinations show larger k at the cavity zone, yet also have the smaller values of lower-nappe aeration. Despite the CTPE with k-mixture show the lower levels of k w , it is the model combination with the highest amount of lower-nappe aeration. However, the turbulent viscosity (ν t ) presented in Figures 13 and 14 expose the lower-nappe aeration is more significant where the ν t is higher, which is displayed clearly by the kand k-mixture models.
These results are justified by the conception of the models and their performance in the cavity zone, as reported by Bardina et al. [70]: "[the k-ω SST includes] a limitation of the growth of the eddy viscosity in rapidly strained flows. [...] The shear stress transport (SST) [...] improves the prediction of flows with strong adverse pressure gradients and separation.".
Additionally, in the CTPE, the turbulent dispersion is proportional to ν t as defined in Equation (12). Oppositely, the TPVoF model equations do not comprehend the turbulent dispersion. Hence, any air-water mixing is due to the numerical diffusion in the interface, which is reduced with higher mesh resolution. The k-mixture model exacerbated lower-nappe aeration is explained by the conception of the model and the problems of the kformulation in rapidly strained flows and with large pressure gradients [70]. The k-mixture model solves a single set of equations for the air-water mixture, and the variables are calculated through an effective density-weighted average (19). Thus, this approach increases the gradients, especially at the air-water interface, which can be aggravated with the increase of the mesh resolution. This problem is not so significant in the k-, because the implementation in the TPVoF does not include density explicitly, which, on the other hand, may lead to several turbulent-related issues at the fluids interface (i.e., spurious velocities and unrealistic turbulence production) [54,55]. Especially the high k generation at the interface is observed in Figure 11. Moreover, the ktype models exhibit minimal turbulent viscosity next to the chute bottom, which also conditions the inner turbulence fields.

Computational Cost
Although the same numerical schemes are employed for all model combinations, for stability reasons, the time-step definition required a maximum Courant number of 0.7 for the TPVoF and between 0.25 and 0.4 for the CTPE model. Nevertheless, the CTPE with k-mixture for the higher mesh resolution is unstable for intermediate and high mesh resolutions. Nonetheless, an analysis is done considering the necessary CPU time to attain a stable solution during one second.
Globally, for the same model and mesh resolution, the calculation time is very similar between the two turbulent models. On average, the CTPE is 4 to 6 times slower than the TPVoF. The CTPE with k-mixture calculation time increased, especially due to spurious air velocities. On the other hand, for the mesh resolutions (R m ) of 10, 15, 20, 30, 40 and 60, there is an average increment of 4.5 times between two sequential resolutions. Hence, from R m = 10 to R m = 60, calculation time increases by a factor of 2 × 10 3 . The CTPE with k-ω SST Sato shows very high levels of turbulence near the interface, which is probably due to the overestimation by two-fluid per-phase turbulence models of the turbulent kinetic energy near large scale interfaces with significant shear [71].

Discussion
The TPVoF and the CTPE Reynolds-averaged Navier-Stokes equations models are not directly comparable. The TPVoF is designed for two incompressible and immiscible fluids that share one set of momentum equations. Despite not including air-water mixing, the TPVoF is still used for engineering purposes. The CTPE is designed for two compressible and interpenetrating fluids, including heat transfer, and is more appropriate for highlyaerated flows. Besides the significant difference in the number and complexity of the mass, momentum, and energy conservation equations, the CTPE comprises several complex closure models (i.e., heat transfer, drag, lift, and turbulent dispersion). Thus, the calculation times and instability of the CTPE model are larger, and a single less accurate definition of the boundary conditions or the mesh may lead to pressure or temperature oscillations and solution divergence.
The interface modeling is also extremely different. The volume-of-fluid method is conceived to separate the fluids. Thus, the TPVoF does not comprise air and water mixing, and any air-water mixture present is due to numerical diffusion, especially in low mesh resolutions. In a RANS approach, with a mesh resolution compatible with most hydraulic engineering applications, the VoF model cannot reproduce the diversity of bubbles and air pockets present in the flow.
Oppositely, the CTPE allows the interpenetrability of both phases, which act as dispersed or continuous in each cell, according to the phase volume fraction. The mixing of the phases is mainly promoted by the turbulent dispersion, depending on the local turbulence conditions. However, the CTPE model has some limitations regarding bubbles and air pockets. First, the bubble size is limited to a single diameter, only varying due to pressure or temperature changes. Therefore, no bubble break-up or coalescence is considered. Secondly, the size of the bubbles must be at least an order of magnitude smaller than the size of the cells. Third, the interface tends to be more diffused than in the TPVoF, especially in the presence of high turbulence levels. Fourth, the RANS approach limits the growth of interface instabilities due to the limitations of the turbulence models and the respective spatial and temporal discretizations. The implementation of an expeditious sub-grid bubble model with a multiple-size bubble population in the CTPE could improve the range of spillway applications without significantly increasing the calculation time.
In the design of spillways, the TPVoF is the primary tool if aeration is nonexistent and in the absence of entrained air, benefiting from the easy model setup and low computational cost. The intake is the part of the spillway most suitable for the application of the TPVoF, considering that the flow is commonly in the subcritical regime, presenting low velocities and no air-entrainment. Downstream of the weir, where the flow changes from sub-critical to the supercritical regime, the TPVoF continues to be the best model until the spillway section where, due to the rise of the boundary layer, free-surface aeration initiates, or if in the presence of an aerator. Downstream, for the two analyzed solvers, the CTPE is the only tool that models the phenomena of aeration and air-bubbles transport and dispersion efficiently, despite showing some limitations due to the single-size bubble formulation. In some spillways, downstream of the chute, there is an energy dissipation basin where the flow regime returns to subcritical, and the majority of the air-bubbles rise to the surface and are detrained. Hence, it may be useful to apply the TPVoF between the basin and the river restitution.
Despite not simulating the air-water mixture, the TPVoF may be an appropriate tool for preliminary analysis of flows with low air-concentrations (<20%), where the tracking of the free-surface is dominant. Further research should compare the efficiency of the CTPE with other successful approaches using the VoF or the mixture models combined with a sub-grid bubble model (e.g., [27,37]), in cases where the air-concentration is lower than 15-20%, or even for larger air concentrations, in light of some recent successful applications of the VoF and mixture models for highly aerated flows [16,18,24].

Conclusions
An assessment performed in a spillway offset aerator compares the two-phase volumeof-fluid (TPVoF) with the complete two-phase Euler (CTPE).
As expected, the TPVoF results depend highly on the mesh resolution, showing no airwater mixing. The accuracy for intermediate mesh resolutions is misleading: the lower-nappe aeration tends to null aeration as resolution increases. This is justified by the fact that the TPVoF model equations do not comprehend turbulent dispersion. Hence, any air-water mixing is due to the numerical diffusion in the interface, which is reduced with higher mesh resolution. To respect the TPVoF model conception, the mesh must be fine enough to reproduce the individual air-bubbles, which is inviable in engineering applications of spillways.
The CTPE combined with the k-ω SST turbulence model exhibits the most accurate results. Surprisingly, intermediate mesh resolutions are sufficient to surpass the TPVoF performance with reasonable calculation efforts, and no significant improvements are found in the highest resolutions, which demand exorbitant computational resources. Moreover, compressibility, flow bulking, and several entrained air effects in the flow are comprehended. Nevertheless, the turbulent dispersion of air-bubbles next to the bottom is not accurately reproduced, possibly due to the limitations of the single-size bubble population and the RANS approach to model the air-water interface at the cavity zone. Due to the more diffusive nature of the Euler-Euler approach, the CTPE may show slower convergence. Nevertheless, the lower-nappe aeration is expected to converge to a value similar to the highest resolution mesh. Contrarily to the TPVoF, the CTPE has a turbulent dispersion model which enhances the mixing of the phases.
The k-mixture turbulence model presents an exacerbated lower-nappe aeration, proving inadequate to simulate aeration in interfaces of rapidly strain flows and with high-pressure gradients.
Both the TPVoF and the CTPE show an increase of the maximum chute bottom pressure with higher mesh resolutions, surpassing the reference value, which is linked to the difficulties to mimic the air-concentration distribution next to the bottom.
Overall, despite not reproducing all aspects of the flow with acceptable accuracy, the complete two-phase Euler surpasses the two-phase volume-of-fluid model, evidencing an efficient cost-benefit performance and significant value in hydraulic engineering applications of spillway aerated flows. Further developments are expected to enhance the tool efficiency and stability. Nevertheless, the two-phase volume-of-fluid is appropriate to model the spillway intake and sometimes even the river restitution. Additionally, a comparison of the efficiency of the CTPE with other successful approaches using the VoF or the mixture models combined with a sub-grid bubble model is of major interest to identify the most appropriate model for specific hydraulic engineering applications of spillway aerated flows.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: