## 2 Materials and Methods  
A three-dimensional _CFD_ with a non-adaptive grid that remained fixed during the calculation process was used to simulate the case in hand. The grid consisted of orthogonal rectangular elements and volumetric hexahedral elements, which eases the mesh generation, enhances the calculation accuracy, and reduces the required memory resources. Flow domain boundaries were defined by calculating the obstacle occupied percentage in the grid elements by applying the Fractional Area Volume Obstacle Representation technique (_FAVOR_) [44]. The flow continuity (1) and momentum Equations (2)-(4) were modified to comply with the _FAVOR_ technique. \[\frac{\partial\rho}{\partial t}+\frac{1}{V_{F}}\rho\
abla\left(\overset{ \rightarrow}{q}\right)\ =\ 0 \tag{1}\]  
\[\frac{\partial u}{\partial t}+\frac{1}{V_{F}}\left(\overset{ \rightarrow}{q}\
abla u\right)\ =\ -\frac{1}{\rho}\frac{\partial P}{\partial x}+G_{x}+f_{x} \tag{2}\]  
\[\frac{\partial v}{\partial t}+\frac{1}{V_{F}}\left(\overset{ \rightarrow}{q}\
abla v\right)\ =\ -\frac{1}{\rho}\frac{\partial P}{\partial y}+G_{y}+f_{y} \tag{3}\]  
\[\frac{\partial w}{\partial t}+\frac{1}{V_{F}}\left(\overset{ \rightarrow}{q}\
abla w\right)\ =\ -\frac{1}{\rho}\frac{\partial P}{\partial z}+G_{z}+f_{z} \tag{4}\]  
where (\(V_{F}\)) is the opened fraction volume to flow, (\(x\),\(y\),\(z\)) are the rectangular cartesian coordinates, \(\overset{\rightarrow}{A}\) is the area fraction vector (\(A_{x}\),\(A_{y}\),\(A_{z}\)) opened to flow, \(\overset{\rightarrow}{V}\) is the velocity vector and its flow components, (\(u\),\(v\),\(w\)) are the instantaneous velocity components in rectangular Cartesian coordinates, \(\overset{\rightarrow}{q}\) is the flux vector and its components (\(u\).\(A_{x}\),\(v\),\(A_{y}\),\(w\).\(A_{z}\)), (\(\rho\)) is the fluid density, (\(G_{x}\),\(G_{y}\),\(G_{z}\)) are the water body acceleration in (\(x\),\(y\),\(z\)) directions,(\(f_{x}\),\(f_{y}\),\(f_{z}\)) are the viscous accelerations in (\(x\),\(y\),\(z\)) directions, (\(t\)) is the time, and (\(P\)) is the pressure. The Volume of Fluid algorithm (_VOF_) developed by the authors of [45] was used to capture the fluid free surface. The function (\(F\)) represents the percentage of fluid occupancy percentage for each grid cell, with values that range between 0 and 100%. The perpendicular line to the rapid change direction of the (\(F\)) value represents the free surface inside the grid cell as shown in Figure 1.  
\[\frac{\partial F}{\partial t}+\frac{1}{V_{F}}F\ \
abla\left(\overset{ \rightarrow}{q}\right)\ =\ 0 \tag{5}\]  
The available _CFD_ techniques have a limited ability in simulating the scour and erosion processes, thus limiting their trusted applicability in scour depth estimation [43]. In the present study, to investigate the flow structure during the scouring process, the flow conditions and the bathym changes of the scour hole during the development of the scouring process were obtained from the experimental work found in [46]. The experimental work was conducted at the Darmstadt University of Technology, the Institute of Hydraulic and Water Resources Engineering, Germany. The scour development around a 20 cm diameter (\(D\)), sand-embedded Plexiglas pier was investigated at different time steps (960, 3840, 10,620, 19,620, 34,140, and 70,320 s). The bathymetry of the scour hole at each time step was measured using a Laser Distance Sensor (\(LDS\)) with \(\pm\)4 mm accuracy located at the pier center and vertically driven with \(\pm\)0.02 mm precision. Sand of \(D_{50}\) = 0.6 mm, with the standard deviation of particle size (\(\sigma_{g}\) = 1.37), and an angle of repose = 30\({}^{\circ}\) was used to form a 50 cm depth movable bed material in the working section of a glass-sided rectangular flume that was 26.0 m long with a 2.0 m width, and a 1.0 m depth.