## 3 Numerical Simulations
### Turbulence Model  
Free-surface flows in a stepped spillway feature high turbulence intensity and air entrainment. A proper numerical model is essential to obtain reasonable simulation results. Investigations [12; 13; 14] show that the \(k\)-\(\varepsilon\) turbulence model is suitable for simulations of the spillway flows. [PERSON] et al. [11] used this model in examinations of the flow and pressure field in the conventional stepped spillway (\(\theta\) = 0). Thus, this model is used to reproduce the hydraulic features of the proposed layouts. The governing equations read:  
The turbulence kinetic energy (\(k\)) equation:  
\[\frac{\partial(\rho k)}{\partial t}+\frac{\partial}{\partial x_{i}}(\rho u_{ i}k)=\frac{\partial}{\partial x_{i}}\bigg{[}\bigg{(}\mu+\frac{\mu_{t}}{ \alpha_{k}}\bigg{)}\frac{\partial k}{\partial x_{i}}\bigg{]}+G_{k}+G_{b}-\rho \varepsilon-Y_{M}+S_{k} \tag{1}\]  
The turbulence kinetic energy dissipation rate (\(\varepsilon\)) equation:  
\[\frac{\partial(\rho\varepsilon)}{\partial t}+\frac{\partial(\rho u_{i} \varepsilon)}{\partial x_{i}}=\frac{\partial}{\partial x_{i}}\bigg{[}\bigg{(} \mu+\frac{\mu_{t}}{\alpha_{\varepsilon}}\bigg{)}\frac{\partial\varepsilon}{ \partial x_{i}}\bigg{]}+C_{1\varepsilon}\frac{\varepsilon}{k}(G_{k}+G_{3 \varepsilon}G_{b})-C_{2\varepsilon}\rho\frac{\varepsilon^{2}}{k}+S_{\varepsilon} \tag{2}\]  
where \(\rho\) = flow density, \(u_{i}\) = velocity component, \(t\) = time, \(\mu\) = dynamic viscosity, \(\mu_{t}\) = turbulent dynamic viscosity, \(G_{k}\) = generation of turbulent kinetic energy due to mean velocity gradients, \(G_{b}\) = generation of turbulent kinetic energy due to buoyancy, \(Y_{M}\) = contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, \(\alpha_{k}\), \(\sigma_{\varepsilon}\) = turbulent Prandtl numbers for \(k\) and \(\varepsilon\), and \(S_{k}\) and \(S_{\varepsilon}\) = user-defined source terms. The model constants are \(C_{\mu}\) = 0.09, \(\alpha_{k}\) = 1.0, \(\varepsilon\) = 1.3, \(C_{1\varepsilon}\) = 1.44, \(C_{2\varepsilon}\) = 1.92 and \(C_{3\varepsilon}\) = 1.0. Although [PERSON] et al. [15] and [PERSON] et al. [16] suggested that model sensitivity be tested, the often-used default settings are adopted in the model.  
In combination with the \(k\)-\(\varepsilon\) model, the VOF method is employed to track air-water interfaces with the geometric reconstruction scheme  
\[\frac{\partial\alpha_{w}}{\partial t}+\frac{\partial(\alpha_{w}u_{i})}{ \partial x_{i}}=0 \tag{3}\]  
The continuity and momentum equations are  
\[\frac{\partial\rho}{\partial t}+\frac{\partial\rho u_{i}}{\partial x_{i}}=0 \tag{4}\]  
\[\frac{\partial\rho u_{i}}{\partial t}+\frac{\partial\rho u_{i}u_{j}}{ \partial x_{j}}=-\frac{\partial P}{\partial x_{i}}+\frac{\partial}{\partial x _{j}}(\mu+\mu_{t})\bigg{(}\frac{\partial u_{i}}{\partial x_{j}}+\frac{ \partial u_{j}}{\partial x_{i}}\bigg{)} \tag{5}\]  
where \(u_{j}\) = velocity component, \(\alpha_{w}\) = volume fraction of water and \(P\) = pressure. The max. and min. values of \(\alpha_{w}\) are 1 and 0, corresponding to full water and only air in a cell, respectively.  
The simulations are performed in ANSYS Fluent. Based on the finite-volume method, the governing equations are solved with the SIMPLE algorithm for the coupling of pressure and velocity. The least-square cell-based method models the gradient. For solutions of \(P\) and \(\alpha_{av}\), the PRESTOI and Geo-Reconstruct techniques are used. The second-order upwind scheme is applied for the momentum; the first-order upwind scheme for the turbulent kinetic energy and dissipation rate.