## 3 Numerical Model
### RANS Equations with k-\(\epsilon\) Turbulent Model  
The governing mass and momentum equations for Newtonian, incompressible fluid flow, which constitute the RANS equations, are expressed as below:  
\[\frac{\partial}{\partial x_{i}}(u_{i}A_{i})=0 \tag{1}\]  
\[\frac{\partial u_{ii}}{\partial\mathrm{t}}+\frac{1}{V_{F}}\left(u_{j}A_{j} \frac{\partial u_{i}}{\partial x_{j}}\right)=-\frac{1}{\rho}\frac{\partial p}{ \partial x_{i}}+G_{i}+f_{i} \tag{2}\]  
where \(x_{i}\), \(i=1\),..., \(3\) represents the coordinate along the \(i\)-th direction in space, \(u_{i}\) (\(\kappa_{i}\), \(t_{i}\)) is the velocity, \(A_{i}\) is the fractional area open to flow over the \(i\)-th direction, \(t\) is time, \(V_{F}\) is the fractional volume open to flow, \(p\) is pressure, \(\rho\) is the fluid density, \(G_{i}\) represents the gravity accelerations, and \(f_{i}\) is the \(i\)-th component of the viscous acceleration term, expressed as follows:  
\[rV_{F}f_{i}=\tau_{b,i}-\frac{\partial}{\partial x_{j}}\big{(}A_{j}\tau_{ij} \big{)}, \tag{3}\]  
in which \(\tau_{b,i}\) represents the wall shear stress, and \(\tau\) the shear stress, given by:  
\[\tau_{ii}=2\mu\left\{\frac{\partial u_{i}}{\partial x_{i}}-\frac{1}{3}\left( \frac{\partial u_{j}}{\partial x_{j}}\right)\right\},\quad\text{ and }\tau_{ij}=-\mu\left\{\frac{\partial u_{j}}{ \partial x_{i}}+\frac{\partial u_{i}}{\partial x_{j}}\right\}\,, \tag{4}\]  
with \(\mu\) as the dynamic viscosity. Wall boundary conditions are evaluated differently according to the chosen turbulence method, which uses a law of the wall formulation. In the present study, the standard \(k\)-\(\epsilon\) turbulence closure was adopted, as was done in previous works involving rapidly varied high-Reynolds unsteady flows [16; 25; 26], in order to calculate the turbulence eddy viscosity with the turbulence kinetic energy \(k\) and the turbulent dissipation rate \(\epsilon\) per unit fluid mass. The Volume-of-fluid (VOF) method was used to detect the free-surface water profile. Grid cells are classified as full, empty, or partially filled with fluid, through a number between 0 and 1 for the free water surface in each cell, which is a fraction of water \(F\), equal to 1 when the cell is completely filled with water and to 0 when the cell is completely empty. The VOF function \(F(x_{i}\), \(t)\) is defined using the following transport equation:  
\[\frac{\partial F}{\partial\mathrm{t}}+\frac{1}{V_{F}}\left\{\frac{\partial}{ \partial\mathrm{x}}(FA_{i}u_{i})\right\}=0 \tag{5}\]  
With FLOW-3D, complex geometries can be efficiently modeled by partially blocking each cell in a rectangular mesh using the FAVOR (Fractional Area/Volume Obstacle Representation) technique. However, like other discrete methods, this method is affected by the resolution of the computational grid. With the improvement of the mesh resolution, the solid geometry is better represented. When the mesh cells are partially blocked by a solid object, each computational cell is represented by a single volume fraction (\(VF\)) and three area fractions (\(A_{i}\)). The VF indicates how much volume is occupied by fluid in the cell after subtracting the volume occupied by solids (if there are any). In this case, the fraction of each mesh cell face through which the fluid can flow were determined as area fractions. Volume fraction and the fractional areas for each cell was calculated by the preprocessor. Numerical computations were carried out on a structured finite-difference grid. The grid was staggered so that while scalar quantities of the fluid such as pressure (\(p\)), density (\(\rho\)), viscosity (\(\mu\)), gravitational acceleration (\(G_{i}\)), volume fraction (\(V_{F}\)), and fluid fraction (\(F\)) were computed at the cell centers Only vector and tensor quantities such as velocities (\(u_{i}\)) and areas (\(A_{i}\)) were taken into account at the cell faces. For equations density, dynamic viscosity and gravitational acceleration were known constants at solution initiation.