An improved air entrainment model for stepped spillways

Numerical modelling of flow in stepped spillways is considered, focusing on a highly economical approach combining interface capturing with explicit modelling of air entrainment. Simulations are performed on spillways at four different Froude numbers, with flow parameters selected to match available experimental data. First, experiments using the model developed by Lopes et al. (Int. J. Nonlin. Sci. Num., 2017) are conducted. An extensive simulation campaign is used to carefully evaluate the predictive accuracy of the model, the influence of various model parameters, and sensitivity to grid resolution. Results reveal that, at least for the case of stepped spillways, the number of parameters governing the model can be reduced. A crucial identified deficiency of the model is its sensitivity to grid resolution. To improve the performance of the model in this respect, modifications are proposed for the interface detection algorithm and the transport equation for the volume fraction of entrained air. Simulations using the improved model formulation demonstrate better agreement with reference data for all considered flow conditions. A parameter-free criterion for predicting the inception point of air entrainment is also tested. Unfortunately, the accuracy of the considered conventional turbulence models proved to be insufficient for the criterion to work reliably.


Introduction
Along with a renewed interest in stepped spillways as a flood overflow structure and energy dissipator in hydraulic engineering, attempts at gaining a better physical description of spillway flows have also intensified. A process that is especially challenging to study by means of both physical and numerical experiments, is the self-aeration of the spillway. Yet, since large quantities of entrained air lead to higher flow depths, release energy, and reduce the potential for damage caused by cavitation, accurate prediction of aeration is crucial for spillway design. In this work, mathematical modelling and simulation of air entertainment is in focus. To put the present contribution into context, a brief review of the physics of air entrainement in spillways is given below, followed by an overview of past attempts of accounting for them in a numerical setting.
Generally, air entrainment is driven by turbulent motion and occurs when the turbulent forces at the free surface overcome the stabilizing effects of surface tension and buoyancy [10]. Applied to spillways, it has since the early work of Straub and Anderson [26] been widely accepted that the onset of self-aeration takes place when the turbulent boundary layer, developed from the crest, reaches the free surface. This location is commonly referred to as the 'inception point'. Several contributions consider the onset of the aeration in detail [6; 33; 28], and empirical relations exist for the distance to the inception point from the spillway crest immiscible. Furthermore, RANS turbulence modelling is adopted, leading to the following set of governing equations.
Here the overbar denotes the Reynolds average, ρ is the density, µ is the dynamic viscosity, u is the velocity vector, p ρgh = p − ρg · x is the dynamic pressure, and f s is the surface tension force. The latter is approximated using the Continuous Force Model, see [3] and also [24] for a detailed discussion in the context of OpenFOAM R . The term ρu ⊗ u represents the Reynolds stresses, which are to be approximated by the turbulence model. An algebraic approach to account for the evolution of α l is adopted, with the associated transport equation originally formulated as in OpenFOAM R 's VoF solvers. The third term in the equation is artificial and is meant to introduce additional compression of the interface to ensure its sharpness. However, here the formulation of this term is modified in order to accommodate it into the air entertainment modelling framework. The details are provided in Section 4. The definition of u c is nevertheless not altered: It is aligned with the interface normal, and its magnitude is computed as C α |ū|, where C α is an adjustable constant, here set to 1. Given α, the material properties of the fluids are readily obtained as ρ = α l ρ l + (1 − α l )ρ air , µ = α l µ l + (1 − α l )µ air .
The indices l and air are used to refer to the water and air, respectively. What remains to be discussed is the choice of turbulence model, which for the case of the stepped spillway is far from trivial. In principle, the model should be able to properly account for the interaction between the turbulent and multiphase structures in order to provide accurate prediction in the aerated region of the flow. None of the closures that have found widespread use were developed with this goal in mind. Nevertheless, it is common for conventional two-equation models to be used for aerated flows. In [16], which is the work this article largely builds upon, the k-ω SST model [19] is used for stepped spillway simulations. In [17], it is employed in simulations of a plunging jet and a hydraulic jump. For the latter, many studies also use the k-model and its variations, a comprehensive list of references can be found in Table 5 in [29]. Qian et al. [23] found the realisable k-model to be favourable for stepped spillway flow. Based on the above, we consider both the k-ω SST and the realisable k-model [25] and test which of them leads to better predictive accuracy.
It is pointed out in [11] that in the implementation of the above (and also other) turbulence models in OpenFOAM R the viscous diffusion terms are not treated consistently in regions with a non-zero density gradient. Since in the simulations of the spillway presented here a density gradient is present across the whole aerated part of the flow, this issue can have a significant effect on the results. The authors of [11] also provide alternative implementations, in which the inconsistency is resolved. Here, we test using both the default and the improved implementations.

Numerical methods
The computations are performed using OpenFOAM R version 5, provided by the OpenFOAM Foundation. This CFD tool is based on cell-centered finite volume discretization, which is de facto the industry standard. Two custom solvers are used in the study, implementing the air entraiment modelling presented in Sections 4 and 6. Both of them represent modifications of the solver interFoam, which is distributed with OpenFOAM R . This solver implements the VoF methodology discussed in Section 2.1. The governing equations are solved in segregated manner using a variant of the PISO algorithm [15].
A crucial component of the numerical setup is the selection of the spatial interpolation and time integration schemes. Generally, linear interpolation can be used in space except when considering convective fluxes. In the momentum equation, the latter are interpolated using the limitedLinearV scheme, which is a TVD scheme based on the Sweby limiter. The limiter is computed based on the direction of most rapidly changing gradient and then applied to all three velocity components. This improves stability but at a certain expense in terms of accuracy. For the convection of α in equation (3), a TVD scheme using the SuperBee limiter is used. The van Leer limiter was also considered, but SuperBee led to better results on coarser meshes due to being more compressive. Unfortunately, in a multi-dimensional setting, using a TVD scheme does not guarantee that the values of α will be bounded between 0 and 1. Therefore, OpenFOAM R utilizes an additional flux limiting technique, referred to as MULES. It is based on the Flux Corrected Transport theory developed Zalasak [32], more details can be found in [8]. The convective fluxes in the turbulence equations are discretized using the second order upwind scheme, called linearUpwind. This scheme is unbounded, but no significant effect of parasitic oscillations was observed even on coarse grids. Finally, as discussed in Section 4.2 below, the air entrainment model adds an advection-diffusion equation for the flow variable α g , meant to indicate the distribution of the volume fraction of entrained air , to the system. Here, a TVD scheme using the van Leer limiter is employed.
The first-order implicit Euler scheme was used for time-stepping. The choice is not of particular importance here because the flow eventually arrives to an essentially steady state. Nevertheless, a CFL number ≤ 1 was neccessay to maintain in order to keep the simulations stable. This was achieved using adaptive timestepping.

Simulation cases
This section presents the setup of the stepped spillway simulations used to evaluate the performance of the entrainment modelling. In order to have a reference with respect to which the accuracy of simulation results can be analysed, the flow and spillway parameters are selected to match those in the experiments of Bung [4]. These were performed on four different spillways combining two selections for the angle (θ =18.4, 26.6 • ) with two for the step height (s = 0.03, 0.06 m). For each spillway, measurements were made for three flow discharge values (q = 0.07, 0.09, 0.11 m 2 s −1 ). The parameters θ, s, and q can be used to construct the step Froude number, which can be considered the main controlling parameter of the flow [18; 5], Here K = s cos θ is the step induced macro-roughness and g is the acceleration due to gravity. The experiments of Bung cover twelve different Froude numbers in the range, 2.7 ≤ F s ≤ 13. For the simulations four values fairly evenly distributed across this range have been selected: 2.7, 4.6, 8.3, and 13.0. The values of θ, s, and q in the four simulation cases are provided in Table 1. This table also provides the values of some auxiliary geometrical parameters, the definition of which can be found in Figure 1. The figure also shows the origin and orientation of the employed Cartesian coordinate system.  All the simulations are performed on 2D domains. This is chiefly motivated by the fact that the investigated modelling methodology is low-fidelity and most suitable for quick evaluations of the integral characteristics of the flow. Any 3D effects due to sidewalls are expected to be negligible with respect to the overall accuracy of the flow predictions. Additionally, using 2D domains it seems to be ensured that, even on dense grids, no air entrainment is resolved by the VoF. This is shown by the interFoam computations on grid G4. By contrast, in a 3D setting, it is more likely that some interface perturbations eventually start getting captured. Investigating the performance of the entrainment modelling in such a scenario is of interest, but lies out of scope of the current work.
The same boundary conditions were used for all cases, with the exception of the discharge q prescribed at the water inlet. The height of the inlet was adjusted to ensure sub-critical inflow conditions. A zero gradient condition was used for the pressure, while Dirichlet conditions were applied to k, , and ω. The values were set assuming 5% turbulent intensity and 10% of the critical height as the turbulent length scale.
For the outlet, a zero gradient condition was prescribed for velocity, pressure, α l and α g . For k, , and ω the OpenFOAM inletOutlet boundary condition was applied. It acts as a zero gradient condition in case of outflow, but for backflow a homogeneous Dirichlet condition is applied instead.
No slip conditions were used at the walls, with a zero gradient condition set for α l and α g . The turbulent quantities were estimated by regular wall laws, in OpenFOAM named as kqRWallFunction, epsilonWallFunction, omegaWallFunction and nutkWallFunction, for k, , ω, and ν t respectively.
For the top boundary the total pressure was fixed, and a pressureInletOutletVelocity condition was applied for the velocity. Similar to inletOutlet, this imposes a zero gradient for outflow, whereas for backflow, it assigns a velocity based on the flux in the patch normal direction. The inletOutlet boundary condition was used for α l , α g , k, , and ω.
The material properties of the fluids were set to correspond to air and water. The values are provided in Table 2.
The computational grids were constructed using Pointwise R , and consist of square cells with the exception of a small strip close the top boundary, where unstructured meshing was necessary to account for the slope of the geometry. Four grids with increasing cell density, denoted G1, G2, G3, and G4, were constructed for each of the four spillways. In each consecutive grid the edge length of the square cells is halved. On the coarsest grid G1, the edge length is 5 mm, which corresponds to what was used in the simulations by Lopes et al. [16]. This can be related to the the critical height of the spillway flow, defined as h c = q 2 /g 1/3 .

Property
Value Liquid density, ρ 1 1000 kg/m 3 Gas density, ρ 2 1 kg/m 3 Liquid kinematic viscosity, ν 1 1 · 10 −6 m 2 /s Gas kinematic viscosity, ν 2 1.48 · 10 −5 m 2 /s Surface tension coefficient, σ 0.07 Depending on the flow case, on the G1 grid, h c is discretised by either 15 or 21 cells. The numbers for the G4 grid are, respectively, 126 and 171. The densities are not adjusted to remain equal with respect to h c across all flow conditions, because experiments showed that the relevant parameter for entrainment modelling is the resolution of the interface. The number of cells in each mesh is given in Table 1.
In conclusion, additional characteristic scales of spillway flow are defined. These will be used for nondimensionalising the results. At a given x, the height h 90 is defined as the y-coordinate of the point where α air = 0.9. The velocity u 90 is defined as the x-component of the mean velocity vector at y = h 90 . Similar scales can be defined with respect to other α air values, e.g. h 50 .

Air entrainment modelling
This section presents the air entrainment model developed by Lopes et al. [16]. One can split the model into three components: an estimator for the flux of entrained air, a transport equation for the volume fraction of entrained air, and a coupling procedure between the model and the VoF framework. Sections 4.1, 4.2, and 4.3 each focus on one of these components. Additionally, for the stepped spillway, estimating the location of the inception point is necessary and this constitutes an additional component of the model, which is treated in Section 4.4.

Estimating the flux of entrained air
A key component of the model is the estimation of the quantity of entrained air carried passed some imaginary surface located below the interface. The estimate was proposed by Ma et al. [17]: where a is a length scale associated with the roughness of the interface due to turbulence, and n is the interface normal defined as Here ε is a small number added for numerical stability. It is assumed that entrainment is confined to a layer of thickness φ ent close to the surface. Therefore, in order to obtain a volumetric air entrainment rate, q can be divided by φ ent . Note, however, that (6) is by definition not restricted to being non-zero only in the vicinity of the interface. Theoretically, entrainment can be incorrectly predicted in regions where it should not take place. For this reason, in [16], q is additionally multiplied by some function δ f s , which is non-zero only close to the interface. The final form of the volumetric air entrainment rate estimate is It remains to define how a, φ ent , and δ f s are computed. The common approach for a is to equate it to the turbulent length scale as predicted by the RANS model. The value of φ ent should be related to some characteristic length scale of the problem.
Within the VoF framework, the α l -field stands out as the natural choice as a basis for the development of an interface indicator function such as δ f s . Typically, the interface is defined as the isosurface α l = 0.5, however this is only accurate when the interface is sharp. In the presence of air entrainment, a more robust metric is the magnitude of the gradient of α, which can be expected to reach its maximum close to the boundary between the continuous air region and the air-water mixture. Hänsch et al. [12] used the gradient of α and a function based on tanh to find the interface as part of their air entrainment model. This function was adopted by Lopes et al. [16] and reads as Here |∇α l | cr is a constant representing the critical value of the gradient that is expected to be reached in the interface cells. Its estimate can be computed based on the size of the grid cell, ∆x: |∇α l | cr = 1/(4∆x). The parameter β can be used to control the extent of the interface region with respect to the chosen |∇α l | cr , and thus provides an opportunity to broaden or restrict the number of cells in which the source term is active.

The α g -equation
The source term (8) is introduced into an additional equation for the modelled volume fraction of entrained air, α g : Here ν t denotes the turbulent viscosity. The velocity of the entrained air, u g , is either set to be equal to u or alternatively modified according to [7]: where the correction velocity u r is calculated based on a bubble radius according to The inclusion of the diffusion term in (10) is considered optional. Additionally, Lopes et al. [16] argue that to properly account for the break-up of bubbles at the free surface, α g should be set to zero when α air exceeds a certain threshold value, referred to as the BBA. The suggested value to use is 0.1.
The exact physical meaning of α g and its relation to α l are somewhat elusive. In [16], the authors discuss the possibility of using the entrainment model without backward coupling to the VoF solver. In this case, the situation is clear: α g shows the modelled distribution of the volume fraction of entrained air, which cannot be captured by the VoF. However, when the coupling is two-way (particulars presented below), the idea is that the entrained air should be captured in the α l field, and α g is essentially reduced to a buffer-field used to propagate the effect of S g onto α l .

Coupling to the VoF solver
Here, we are interested in applying the model in a two-way coupling regime, meaning that the model's predictions should be propagated into the distribution of α l . The premise is that the VoF simulation by itself does not resolve any entrainment, and therefore all of it is accounted for by a subgrid model based on the α g equation (10). The overall idea is that α l should be reduced in regions where α g is large, and in a manner that does not disrupt mass conservation.
Here this is done through a modification of the artificial compression term introduced into the α lequation (3): The term (1 − α l )α l = α air α l is originally meant to serve as an indicator for cells constituting the interface, in which the compression is to be applied. The key observation is that multiplying u r j by some negative number instead would lead to interface expansion and thus a region occupied by a mixture. The goal is then to correlate α g with the change in sign in the term in front of u r j . The most obvious way to do that is exchange α air α l for (α air − α g )α l . Note that since (13) is a transport term, mass conservation is guaranteed. The modified α l -equation then reads Under the definition above, the model is active only when α g > α air , which is reasonable. It is also worth mentioning that otherwise (13) recovers its original compressive function. This occurs even in the regions occupied by a mixture, which can be called into question. As part of the work on improving the model, some experiments have been conducted in which positive values of α air − α g where cut to 0, however the exhibited results were inaccurate, and introducing such a discontinuity is probably best avoided.

Inception point estimation
As discussed in the introduction, surface aeration initiates when the turbulence perturbations exceed the stabilizing forces of surface tension and buoyancy at the free surface. In the model of Lopes et al. [16] no attempt is made to explicitly compute the force balance. Instead, two model parameters, k c and u c are introduced, where the subscript c stands for critical. The inception is considered to occur when k > k c and u · n > u c and u · g > u c .
Appropriate values for k c and u c are extremely difficult to predict a priori, since the selection clearly depends not only on the flow conditions (see Section 5.1), but also on the turbulence model and its prediction of k. Careful calibration with respect to the selected model is therefore necessary. In [16], the authors nevertheless suggest u c = 0.8 m/s and k c = 0.2 m 2 /s 2 , referring to previous experimental results. Unfortunately, how these values relate to the characteristic length and velocity scales of the flow is not discussed.

Stepped spillway simulations with the original model
This section demonstrates results from simulations performed using the model of Lopes et al. [16] described in the previous section. In the original source the model was reported to reproduce experimental data on a stepped spillway with F s =2.7. However, the simulations were performed on relatively coarse grids, with the grid sensitivity study performed using only the baseline VoF solver. Furthermore, initial testing within this work indicated that its effect was reduced upon grid refinement, which motivates the analyses herein.
The implementation of the corresponding solver, called airInterFoam was kindly provided by P. Lopes via personal communication. Below, we abbreviate airInterFoam to AIF. In the following sections, the solver is evaluated in terms of sensitivity to flow conditions (Section 5.1), grid resolution (Section 5.2), and several parameters of the entrainment model (Section 5.3).

Sensitivity to F s
Here, results from four AIF simulations varying in the prescribed Froude number of the spillway flow are presented. To highlight the effect of the entrainment model, data obtained with the baseline VoF solver interFoam (abbreviated IF henceforth) are also provided. The employed numerical parameters are based on [16], where good results for the F s =2.7 case (simulated in 3D) are presented. In particular, the G1 grid is employed, k-ω SST is used for turbulence modelling, the diffusion term is omitted in the α g equation, and k c = 0.2 m 2 /s 2 , u c = 0, u g = u l , BBA = 0.1, and φ ent = 0.05h c . Figure 2 shows the obtained values of α air in the uniform flow region. Good accuracy is achieved for F s =2.7 and 4.6, but for the two higher F s the model fails to predict the reduced penetration of air into the corners of the steps. As a result, in terms of magnitude, the errors in the IF and AIF simulations are similar, although in the case of IF the diffusion of the interface is a purely numerical effect. It is also interesting to note that for F s =2.7 the accuracy is on par with the 3D simulations using similar model settings conducted in [16].  The evolution of the surface elevation, measured as h 90 , is shown in Figure 3. The elevation's value in uniform conditions is well-predicted for all Froude numbers. However, the location of the inception points are not captured as consistently. The difference in the obtained values with respect to the experimental data of Bung is provided in each plot of the figure: ∆n i stands for the difference in terms the step number, and ∆L i in terms of x. It should be noted that when comparing across different F s , using ∆L i is more appropriate, since for the two higher Froude numbers, the length of the step is halved relatively to the lower Froude number cases. The obtained incetion point locations for F s =2.7 and F s =4.6 are reasonably accurate, but, unfortunately, at higher F s the disagreement with the experiment becomes larger. Furthermore, the predicted inception point for F s =8.3 is further downstream as compared to that for F s =13, whereas the experimental data exhibits the opposite trend. Figure 3 additionally shows the experimental values of h w , which is the equivalent clear water depth, i.e. the surface elevation that should be predicted by IF. In the obtained results, IF somewhat over-predicts h w , the reason being the coarseness of the grid.
The predicted profiles of the streamwise velocity are shown in Figure 4. Remarkably, no effect of air entrainment modelling is visible, and accurate profiles can be predicted with IF. This result was reproduced in all the simulations in this paper, and, for that reason, velocity profiles are not further presented or discussed.
The principle conclusion from the obtained results is that the settings used in [16] for the F s =2.7 case fail to provide consistently accurate results as the Froude number becomes larger. This indicates that some parameter values of the model, for example k c , should be made a function of F s .

Grid sensitivity
To explore the grid sensitivity of AIF, the F s =2.7-case was run on all four grids G1-G4. The resulting α air and h 90 profiles are shown in Figure 5. Clearly, the behaviour obtained on the coarse grid G1 is substantially changed as the grid gets denser. With increasing resolution, less air is distributed towards the   pseudo-bottom, and for the densest grid only a tiny air layer is found close to the surface. The profiles, both h 90 and α air , approach the corresponding solutions obtained with IF. Obviously, this is caused by the fact that less numerical diffusion contribute to the transport of α g as the grid is refined, but inspection shows that this is also caused by the shrinkage of the area in which the source term S g is non-zero. This, in turn, is controlled by δ f s , which makes the definition of this function a contributor to the grid sensitivity. A more elaborate discussion follows in Section 6. On the other hand, with respect to h w , the IF solution consistently improves with grid refinement. On the G4 grid the interface is perfectly sharp and the h w profile is very well-matched.

5.
3. Sensitivity to u r , bubble breakup criterion, and α g -diffusion Here we explore the effects of the entrainment model parameters that could arguably be considered non-essential or optional. First, the impact of the air bubble drift velocity model, as given in eq. (12), is investigated. This is followed by an analysis of the bubble breakup criteria, BBA. Finally, the model is tested in terms of activation of the diffusion term in the α g -equation (10). The rest of the numerical setup is similar to that used in the Froude number sensitivity study, see Section 5.1.
Here we restrict the analysis to α air profiles in the uniform flow region. For the u r model, three values of the bubble radius r b are considered, along with setting u r to 0. The result is shown in Figure 6a. Clearly, including u r leads to reduced air entrainment, and this effect becomes larger when the input bubble radius is increased. This is expected, since u r is, by definition, directed upwards. Since the intensity of air entrainment is already heavily dependent on the formulation of δ f s , having an additional controller in terms of u r introduces unnecessary complication. At least in the case of spillway flow, setting u r = 0 can therefore be recommended.
Three values of the BBA are considered. The first is 0.1, which is recommended by Lopes et al. [16], and the other two are 0.05 and 0. Recall that the chosen value refers to the minimal volume fraction of water for which α g is allowed to be non-zero. The results from the three simulations are shown in Figure 6b. As expected, slightly more air is entrained close to the h 90 when the bba value is reduced. Note that compared to the experimental data, the air fractions predicted by AIF close to the interface (from y/h 90 ≈ 0.8) are too low. Thus, deactivating the BBA-criterion produces a small improvement in the accuracy meaning that this parameter can be safely removed from the formulation of the model.
Finally, the effect of the diffusion term in (10) is investigated. Two simulations, with and without the diffusion term included, have been conducted, see Figure 6c. The effect of the diffusion appears to be completely negligible. Consequently, it is possible to remove it from the model.     In summary, u r , BBA, and α g -diffusion can be excluded from the model, significantly simplifying its formulation.

Sensitivity to β
As shown in Section 5.2, at high grid resolutions the effects of the model becomes negligible, or even deactivated. The most important controller of the breadth of the region (in terms of ∇α l ) where entrainment is introduced is the form of δ f s , see Eq. (9). In particular, the parameter β determines the breadth of the tanh function, meaning that with smaller β the region of the model activation becomes larger. In principle, it may thus be possible to maintain an appropriate level of entrainment at high grid resolutions by adjusting β accordingly. However, for this to be possible in practice, the necessary change to β should be easy to predict a priori.
Multiple simulations across different Froude numbers and grid resolutions have been conducted in an attempt to determine whether clear guidelines for setting β could be established. Unfortunately, these efforts were fruitless, and the results of using a given β value change significantly depending on the flow and parameters of the simulation.
As an illustrative example, Figure 7 shows the α air profiles produced using β = 10 and 25 in simulations at different Froude numbers on the G2 grid, and compares it to the corresponding profile produced by AIF, where a β-value of 100 is used as default. For F s =2.7 the sensitivity to β is rather small, but for F s =4.6 a significant increase in aeration occurs. For larger Froude numbers a lower β mainly leads to more air being present in the corner of the steps. While the change in the α air profiles appears to be rather small, uniform flow conditions are not reached, and its effect continues to grow further downstream. Similar sporadic behaviour was observed with respect to other simulation parameters. For example, contrary to what is observed when using the G2 grid, using G3 instead leads to the F s =2.7 case becoming very sensitive β.

Proposed model developments
In the previous sections, two issues with the entrainment model proposed by Lopes et al. [16] have been identified. Perhaps the most critical one is the successive deactivation of the model upon grid refinement. The other one is the difficulty in prescribing the k c value in order to get a good prediction of the inception point. In this section, improvements to the model are proposed aiming at alleviating these problems. First, an alternative formulation for the δ f s function is introduced in Section 6.1. Afterwards, amplification of the diffusion term in the α g is argued for in Section 6.2. Finally, a different way of predicting the inception point is discussed in Section 6.3.

Free surface detection
Before discussing the formulation of δ f s , it is necessary to consider what kind of restriction of S g in space is needed in the case of spillway flow. A typical distribution of an unrestricted S g is shown in the left plot of Figure 8. The first thick yellow line is located right below the interface and represents the region where the entrainment can be expected to take place. However, a discontinuous strip of non-zero values is also observed close to the psuedobottom along with more or less randomly distributed points of activation in the corners of the steps. Physically, no entrainment can occur in these regions, and the δ f s function should filter them out. Note that the spatial separation between the correct and non-physical regions of source term activation is not large, which explains why defining δ f s in a universal way that fits all flow conditions and numerical settings is not trivial.
To arrive to a better formulation for δ f s , it is important to clearly understand the deficiencies of the original definition, see Eq. (9). With ∇α l,cr set to 1/(4∆x), the distribution of δ f s over ∇α l depends on two quantities: β and ∆x. It is instructive to see how δ f s changes shape when the values of these parameters are changed. In the left plot of Figure 9, δ f s (∇α l ) is shown for ∆x values corresponding to grids G1-G3, and the two values of β considered in the sensitivity study in Section 5.4. The tanh function defining the transition region of δ f s from 0 to 1 is centred at ∇α l,cr , shown in the figure with black vertical lines. As the grid is refined, this location is shifted to the right, and for larger values of β, the tanh only spans a limited range of high ∇α l values. On the other hand, for smaller β, the tanh becomes so wide that δ f s remains non-zero everywhere.
This behaviour of δ f s should be related to how ∇α l is typically distributed across y, see the right plot in Figure 8. Note that the high values of ∇α l are always restricted to to a relatively thin region close to the interface. Consequently, for a large β, for example β = 100 as proposed in [16], δ f s will be restricted to an increasingly smaller region in space when the grid gets refined. This explains why the diminishing effect of the entrainment with grid refinement observed in Section 5.2.
As mentioned above, it is easy to make the δ f s function less restrictive by lowering β. However, the issue here is that the β and the effective cut-off value in terms of ∇α l are not intuitively related. This is problematic given that the margin of error is quite small, as discussed above in relation to the typical distribution of S g . Alternative options for the expression of δ f s were explored, including a parabola based function, and a step-shaped function defined purely based on the distance from a defined interface. For both alternatives, the idea was to set the value of the critical gradient of α l relatively tight, to capture the upper peak in the gradient plot in Figure 8, and then expand its prevalence away from the these locations according to an appropriate function or logic. In the distance based alternative, δ f s was set to 0 or 1 for a particular cell, depending on its distance from the defined interface. If this distance was less then the interface thickness, φ ent , δ f s was set to 1, and otherwise its value was set to 0. However, this led to step-shaped profiles of α l , where too much air was entrained below the interface, and indicated a need to apply some functionality to reduce the effects of the source term as the gradient of α l is reduced below its critical value. Acknowledging the above, a parabola based function was considered. The possibility to define ∇α cr as a top point of the function, and to define a cut-off value for the gradient as an intersection point, was viewed as beneficial features of this function in the current setting. The latter leading to the possibility of avoiding the long tail in the tanh-function, with the corresponding uncontrollable potential for generation of none-zero values of S g within the steps. The parabola-based δ f s formulation is proposed as Here, d refers to the distance from the vertex of the parabola to its focus, which can be computed as where ∇α l,cut is an input parameter explicitly defining the lowest ∇α l for which the source term may assume non-zero values. The proposed δ f s , computed for grids G1-G3 and two different values of ∇α l,cut , is shown in the right plot of Figure 9. The non-zero values of the function are always fixed to the interval [∇α l,cut , ∇α l,cr ], which expands upon grid refinement. Unfortunately, due to this expansion, even with this new δ f s , the region of non-zero S g values shrinks as the grid is refined. However, the process is slowed, since δ f s is left to be non-zero at lower ∇α l . During initial testing, it was observed that the parabola-based δ f s was very effective at filtering out the sporadic source term activation in the corner of the steps, while preserving the region where entrainment is expected. However, the secondary strip of non-zero S g values close to the psuedobottom (see left plot in Figure 8) would sometimes still be left unfiltered. This is likely related to the secondary peak in ∇α l . To rectify this, the parabolic surface indicator function in Eq. (16) is combined with a distance based approach like the one outlined above.
In details, the values of δ f s computed according to Eq. 16 are additionally manipulated as follows. First, the cells in which δ f s ≥ 0.9 are selected. These should lie near the interface and are therefore likely to belong to the region where S g should be activated. For the remaining cells, the distance to the nearest cells with δ f s ≥ 0.9 is computed. If this distance exceeds the interface thickness, φ ent , the value of δ f s in the cell is set to 0. Overall, this combined formulation gave improvements compared to the original formulation based on tanh. It is noted that even without the distance cut-off modification, improved results with the parabolic δ f s could be achieved, and the δ f s ≥ 0.9 criterion can give false positives.
Finally, making both ∇α l,cr and ∇α l,cut grid independent parameters has briefly been tested, but abandoned due to the sensitivity of the results to the chosen values. In general, our extensive efforts and experimenting with δ f s formulations and other simulation settings showed that constructing a robust and accurate δ f s is very challenging. The sensitivity of the results to any changes in the simulation parameters or flow conditions tends to be very strong. Nevertheless, as shown in the Section 7, the proposed δ f s does represent an improvement with respect to prior art.

Modelling air propagation into the corners of the steps
By definition, the employed entrainment model is meant to account for aeration occurring close to the free surface, within some layer of thickness φ ent . However, experimental results clearly show that in the case of the stepped spillway, air penetrates all the way down to the surface of the steps, see e.g. the experimental profiles in Figure 2. The physical mechanism through which this occurs is described by Pfister and Hager [22]. Inspection of video recordings from their experiments reveals the occasional generation of air troughs that extend from the surface into the bulk flow. These troughs penetrate deep enough to hit the step edges, and when they do, the air is distributed into the steps.
Capturing this intrinsically transient process in a steady state model is not straightforward. Here, we consider using a somewhat ad-hoc approach, taking advantage of the diffusion term in the α g equation (10), ∇ · (ν t ∇α g ). Recall that in Section 5.3 it was shown that the effect of this term on the solution is essentially negligible. However, the effect can be easily amplified by pre-multiplying it with some constant C t . The increased diffusion of α g will then lead to air being redistributed more evenly across y, and consequently result in stronger aeration closer to the pseudobottom.
The difficulty lies in the choice of the value of C t since there is no clear physical analogy between the modelled phenomenon and diffusion. In light of this, it was attempted to search for a suitable value through experimentation, to see whether one leading to improved results across all the considered Froude numbers could be found. As a result, C t = 150 was selected.

Inception point prediction
Even if AIF predicts the aeration onset correctly when appropriate values for k c and u c are supplied, the inception point estimation using this method completely depends on user input. As shown in Section 5.1, the appropriate critical value k c depends on the flow conditions. An alternative approach, found in the work of Hirt [13], is to directly consider the balance between the energy of turbulent motion and that of gravity and surface tension. Defining the source term S g is activated only in cells where P t > P d . This way the inception point prediction requires no user input. However, it completely relies on the correct prediction of k and a by the turbulence model.

Simulations with the improved model
This section is dedicated to evaluating the effects of the model improvements proposed above. To that end, a new solver incorporating these changes has been implemented. Reflecting the focus on stepped spillway simulations, the solver is called spillwayFlow 2 , which is abbreviated to SPF below. The robustness of the model with respect to grid resolution is evaluated in Section 7.1. Results from application of the model to spillway flow at all four considered Froude numbers are presented Section 7.2. The proposed criterion for inception point location is tested separately, in Section 7.3.

Grid sensitivity
Here the new model is put to the same grid sensitivity analysis as presented in Section 5.2 for AIF. To simplify the analysis, the new inception point prediction approach is not employed, and the original criterion based on k c is used instead. Simulation results for the F s =2.7 case obtained on grids G1-G4 are shown in Figure 10. Here, C t = 0 and the diffusion term in the α g equation is thus inactive. As anticipated, the results still depend on the grid, and the general trend is convergence towards the IF solution. Note that in the α air profiles, the reduction of aeration manifests itself predominately at y/h 90 0.6. Closer to h 90 results remain acceptable even on the G4 grid. Related to the absence of air in the lower parts, the predictions of h 90 itself are more sensitive, and unfortunately on the finer grids the accuracy is poor. Nevertheless, compared to the original AIF results (see Figure 5) the robustness of the model is improved.   [4]. Figure 10a shows the surface elevation, Figures 10b-10e, the void fraction profiles at different steps. Figure 11 shows the results obtained with C t = 150. As expected, a comparatively more even distribution of α air across y is achieved, in particular for the simulations on denser grids. Furthermore, a very clear improvement in the robustness of the model is evident, with much more similar results obtained on all four grids.  (e) Step 19 Figure 11: SPF simulated with αg-diffusion (Ct = 150) and Fs=2.7 at different grids (G1-G4) compared to physical model results by Bung [4]. Figure 11a shows the surface elevation, Figures 11b-11e, the void fraction profiles at different steps in the developing region.

Results for different Froude numbers
location in the experimental data. For completeness, profiles corresponding to both C t = 0 and C t = 150 are shown in all the figures. For comparison, they also include results from AIF simulations on the G1 grid, and also from IF simulations on the G4 grid.
The α air profiles are discussed first, see Figure 12. Qualitatively, the same behaviour with respect to C t is observed for all F s . With C t = 0 the distribution of α air across y is close to step-wise, with decent agreement with experimental data for y/h 90 0.7. When C t = 150, the profiles are smoothed out, which generally increases the accuracy. The exception is the F s =8.28 case, for which the air volume fraction at low y/h 90 becomes excessive. Yet even for this case the agreement is better than what could be achieved with AIF. (d) Fs=13 The surface elevation plots are shown in Figure 13. Overall, C t = 150 leads to better results, which agree well with the experimental data in the uniform flow region. An interesting exception is the F s =4.6, for which C t barely has influence on h 90 in the uniform flow region, whereas in the developing region C t = 0 leads to very good agreement with the experiment. However, it is unlikely that this is explained by any fundamental property of the flow or the model. Compared to AIF, the accuracy of the new model is generally on par. AIF curves are marginally closer to experimental data for the two lower F s , and the other way around for the two higher F s .

Inception point analysis
In this section, the source term activation criterion presented in Section 6.3 is tested. Recall that the inception point location is determined from the flow, and depends heavily on the employed turbulence model. Here we show results from simulations using four models: the k-ω SST and realisable k-from the standard OpenFOAM library, and their respective counterparts in the library by Fan and Anglart [11], in which the density gradient is properly accounted for in the transport equations. The latter are referred to by adding varRho to the name of the model. Simulations using standard turbulence modelling and the k c -criteria purposed by Lopes et al. [16] are added as reference. The simulations are performed on the G3 grid.
The results are summarized in Table 3, and the corresponding distributions of α air are shown in Figure 14. For the two models from the standard library, the inception of entrainment is triggered immediately (downstream the crest) for all the considered Froude numbers. When the varRho variants are employed instead, the inception point shifts downstream. This reflects the fact that these model predict significantly lower values of k. For k-ω SST the agreement with experimental data is nevertheless poor, but for the realisable k-model the results are more promising.
In the results produced using standard turbulence and k c = 0.2 m 2 /s 2 , the inception point is predicted in relatively good agreement with the experimental results for all cases but the F s =13 case, where the correspondence is rather bad. This differs from the results attained for AIF at the G1 grid (see Section 5.1), where the inception points were poorly predicted for both F s =8.3 and F s =13. Compared to the k c = 0.2 m 2 /s 2 criterion, the automatic activation criterion performs on par using the realisable k-turbulence model and the variable density turbulence framework.
Overall, none of the tested models perform well enough to use the proposed source term activation criterion. The improved performance of the varRho models shows the importance of properly accounting for the density gradient in the transport equations for the modelled flow quantities. Improving turbulence modelling accuracy near interfaces is an active area of research. The results presented here warrant a deeper investigation of what models are appropriate for the stepped-spillway flow. k-ω SST, k c = 0.2 1.05 1.14 -0.09 -0.9 8.3 G3 realisable k-, k c = 0.2 0.86 1.14 -0.28 -2.9 8.3 G3 k-ω SST 0.00 1.14 -1.14 -12.0 8.3 G3 realisable k-0.00 1.14 -1.14 -12.0 8.3 G3 varRho/k-ω SST

Conclusions
This study presents developments in numerical modelling of self-aeration in stepped spillways. The model of Lopes et al. [16] is taken as baseline, and a large simulation campaign is conducted in order to explore its properties: robustness with respect to flow conditions, grid resolution, as well as sensitivity to the model parameters. The simulations were performed for spillway geometries and inflow discharge values used in the experiments of Bung [4], and cover four step Froude numbers in the range from 2.7 to 13. The corresponding experimental data was used as reference.
The results showed that for the case of spillway flows, three of the model parameters could be removed without loss of generality, making simulation setup easier. The main weakness of the model is shown to be its significant sensitivity to the density of the grid. In particular, with increased resolution, the effect of the model diminishes until it is, essentially, no longer active. Nevertheless, at selected grid resolutions, the demonstrated accuracy of the model was acceptable for all considered step Froude numbers. Interestingly, the prediction of the mean velocity profiles was shown to not be affected by air entrainement modelling, and good results could be achieved using only the underlying VoF solver.
The main reason behind the model's deactivation on dense grids has been identified to be the form of δ f s , which is the function used for limiting the activation region of the volumetric air entrainment source term, see Eq. (9).The region of non-zero values of δ f s shrinks as the grid gets refined, and, in the limit, the source term is set to zero in the whole domain, regardless of flow conditions. To address this issue, a new formulation for δ f s is proposed, combining a parabolic profile with distance-based cut-off. Simulations reveal that while fundamentally the results still depend on the grid resolution in the same manner, under the new definition of δ f s the robustness of the model is improved.
As an additional modification, amplifying the diffusion term in the α g -equation (10) is proposed in order to account for the propagation of entrained air into corners of the steps. Results reveal that this leads to a significant improvement in predictive accuracy, the new model performing better than the original [16] across the whole considered range of F s numbers. Furthermore, the robustness of the model with respect to grid resolution improves significantly as well. It should be acknowledged that the selection of the value of the diffusion coefficient, C t = 150, is currently not physically motivated and can, therefore, be called into question. Nevertheless, we believe that the possibility to use the same value across different flow conditions and the clearly demonstrated advantages in terms of the performance of the model are sufficiently strong arguments in favour of adopting the proposed modification. Finding a more rigorous connection between C t and the characteristic scales of the flow remains as a line of future work.
Finally, an algorithm for automatic estimation of the inception point is tested. The criterion for the inception point is based on energy balance, as proposed by [13]. The performance of the algorithm are heavily dependant on the underlying turbulence model. Unfortunately, for the four considered models the predictions were not reliable, highlighting the need for a more careful investigation of what turbulence modelling is appropriate for self-aerating multiphase flows.