## Abstract

Dam construction continues its rapid expansion around the world primarily for the purpose of hydropower generation. One important consequence of such projects is local scour at the downstream of the dam caused by outflow of excess reservoir water through spillways or bottom outlets that is associated with high velocities. The scour development endangers the dam foundation and river banks and undermines the stability of the hydraulic structures. In this study, a detailed three-dimensional (3D) flow simulation is conducted to investigate the complex fluid–sediment interactions leading to the formation of the scour hole and ridge systems downstream of a near-bottom jet. Three different bed-load equations, including Meyer-Peter–Müller, Nielsen, and Van Rijn formulas, are applied for calculating the bed-load transport rate. Comparison with a series of available experimental data shows that the Meyer*-*Peter–Müller equation results in better predictions than the two other relations. The performance of different turbulence models to reproduce vertical profiles of velocity and scour characteristic against the experimental data were evaluated. The vertical and horizontal profiles of the scour hole-ridge system are also compared with the corresponding experimental ones*.* The numerical model satisfactorily reproduces the geometric parameters representing the scour hole*.* However, the model overestimates the length of the scour hole.

## INTRODUCTION

Flow discharge systems in a dam reservoir are commonly designed as a jet either on the spillway or along the dam body. A jet is generated when a fluid with a high momentum enters a quiescent ambient fluid or a fluid moving with a lower velocity (Boroomand *et al.* 2007; Nyantekyi-Kwakye 2016). Depending on the location of the jet specified by the tailwater level, it can be a free or submerged jet (Faruque *et al.* 2006). When the ejected flow is bounded by a wall, it is called a wall jet (Nyantekyi-Kwakye 2016). The scouring caused by a water jet impinging on erodible beds can endanger the stability of structures located in its path. One example is the scour hole developed along the right bank of the plunge pool of Tarbela Dam in Pakistan. The structure of the flip bucket was somewhat displaced, which eventually led to the displacement of joints and the development of cracks in the drainage passage. More cases can be found in Yildiz & Üzücek (1994) and Annandale (2006). Accordingly, safe design of hydraulic structures against scouring is of high importance. To investigate scouring mechanisms, influencing factors, flow structure, and local scour dimensions, a vast number of experimental and numerical studies have been conducted.

Dey & Sarkar (2006), Sui *et al.* (2008), Mehraein *et al.* (2011), and Melville & Lim (2013) experimentally investigated the scour of wall jets. An extensive review of such studies is presented by Aamir & Ahmad (2016). They studied the effects of different parameters such as densimetric Froude number, tailwater depth ratio, sediment bed characteristics, and channel side walls on scour hole characteristics (e.g., scour depth, scour length, scour width, and ridge height). Increasing computer power along with the development of robust numerical models has resulted in more efficient approaches for studying scouring through 3D numerical models. Such models offer a more detailed view of flow properties such as velocity and pressure fields. This leads to a better estimation of the physical phenomena associated with complex flow fields and their interaction with structures and sediment particles (Johnson & Savage 2006; Karim & Ali 2000). Salehi Neyshabouri *et al.* (2001) developed a mathematical model to investigate the scouring process by a 2D wall jet on a cohesionless bed*.* The simulation was carried out in three main steps, including simulation of flow, sediment concentration distribution, and bed deformation. To determine the sediment concentration near the bed, the deterministic equation of Van Rijn was used. Qualitative agreements were obtained for scouring and sediment deposition patterns with experimental data. However, comparison of model results with observations for temporal evolution of the maximum scour position showed significant differences. It was found that the model correctly predicts the overall profile of experimental measurement. They concluded that the combination of numerical and experimental models could provide more accurately prediction of time development of the scour hole. Karim & Ali (2000) investigated the ability of the FLUENT CFD package to predict the flow patterns generated by a wall jet on a rigid flat bed as well as a scoured bed. The floor velocity and the bed shear stress for flat and scoured beds were compared with experimental data. The results showed close agreement to measurements. Moreover, the effects of available turbulence models in FLUENT as well as different grid sizes on the bed shear stresses were examined. Based on their results, there were no significant differences in simulated bed shear stresses when different turbulence models were used. Abdelaziz *et al.* (2010) investigated scouring downstream of a rigid apron caused by a 2D submerged horizontal jet ejected from the bottom of a sluice gate using numerical models. They used FLOW3D commercial code along with a developed bed-load sediment transport module incorporated into FLOW3D. Based on their results for using FLOW3D, after *t* = 3 minutes of simulation fair agreement between the simulated maximum scour depth and ridge height with the measurements was achieved. However, after *t* = 8 minutes of simulation, the mentioned parameters deviated from the measurements. It also was found that FLOW3D exhibited insufficient accuracy in estimating the scour hole length. In the case of using the developed model, the scour hole length was simulated with a good accuracy, while the downstream slope of the ridge was overestimated. Boroomand *et al.* (2007) used a 2D numerical model developed in FLUENT to study the effect of an offset jet over an erodible bed. Computed vertical profiles of suspended particle concentration and total load transport value were compared with experimental measurements*.* Moreover, a qualitative comparison of scour hole and ridge forms were conducted between the numerical and typical observed patterns of offset jet*.* The results showed an appropriate agreement with experimental findings. However, the simulated near-bed particle concentration profiles substantially deviated from the experimental values. Epley-Chauvin *et al.* (2014) investigated the development of plunge pool scour hole due to a jet impinging onto a non-cohesive bed. The results indicated that the calculated scour depth and ridge height values were generally in good agreement with the measurements. To evaluate the accuracy of the model in predicting the scour depth and ridge height, the coefficient of determination was used after 10 minutes of simulation. Based on their results, the coefficient of determination for scour depth and ridge height were 0.93 and 0.89, respectively. The study revealed that the impact angle of the jet plays a key role in determining the amount of time needed by the system to reach the morphodynamic equilibrium conditions. Castillo & Carrillo (2016) carried out a numerical investigation of the scouring caused by flows coming out of outlets located halfway to the top of the spillway in a double curvature arch dam employing three complementary, namely, empirical formulae, a semi-empirical methodology, and a FLOW3D CFD model. The study demonstrated that numerical simulations could play an effective role in the evaluation of the safety and stability of the dam structure. Later, Castillo & Carrillo (2017) used four different approaches, including a physical model, empirical formulae, a semi-empirical methodology, and FLOW3D to evaluate the scour downstream of a free surface weir ending in a ski jump. It was found that the choice of the turbulence model and the bed-load coefficient in the Meyer*-*Peter and Müller formula play a substantial role in the accuracy of the numerical estimates. Lee *et al.* (2019) developed a new rheology-based three-phase flow model using the open-source package OpenFOAM for sediment transport problems with a water–air interface. They used their model to study sediment transport under 1D open-channel flow conditions for the purpose of evaluating the performance of two models for particle response time. The model with higher performance was adopted for simulating local scour caused by a 2D horizontal wall jet. Their results showed overall agreement between the simulated and measured bed profiles; however, some discrepancies were observed in the dune height and shape of the dune crest at the later stage of the scouring process. They also demonstrated that the three-phase model is capable of capturing sediment avalanche events on the back face of the sand dune. Fei *et al.* (2017) developed a numerical scour model to simulate the local scour around hydraulic structures using the OpenFOAM software. In their studies, a dynamic mesh method was applied to consider the evolution of the bed. Moreover, a 2D sand slide method was employed for more accurate prediction of the bed bathymetry. The sediment transport process included both bed-load and suspended load. To validate the numerical model, the results were compared with three experimental cases. The results showed that the sand slide model performs well. In the case of suspended load, simulated concentration profiles also showed good agreement with the experimental data. In the wall jet scour case, the authors showed that the computed values of the maximum scour depth and the deposition height were slightly overestimated.

As noted above, there have been many recent research studies focused on scour hole formation by various forms of water jets*.* Most of the recent numerical studies have focused on wall jet scouring using primarily 2D instead of fully 3D approaches. In most cases, where scouring is formed under the influence of the wall jet, the jet is released on a rigid apron upstream of the movable bed. The sharp transition between rigid wall and mobile bed allows for the development of intense turbulent structures at the edge of the apron leading to augmented sediment entrainment and transport capacity. In the present study, the jet is discharged directly over a sediment bed. In such cases, the developed turbulent structures are less intense compared to the case of a jet over a rigid apron, leading to slower formation and evolution of the scour hole. In this research, the capability of the FLOW3D to simulate the scouring process due to a 3D wall jet is investigated for the first time. A comprehensive analysis of various geometric characteristics of scour caused by a 3D wall jet, including maximum scour depth and ridge height along with their vertical (*xz*-plane) and planar (*xy*-plane) characteristics were presented and compared quantitatively with a series of measured data. Most previous studies have focused primarily on the numerical investigation of scour parameters caused by corresponding jet without validating the flow field (e.g., Salehi Neyshabouri *et al.* 2001, 2003; Boroomand *et al.* 2007; Abdelaziz *et al.* 2010; Epley-Chauvin *et al.* 2014; An *et al.* 2015; Lee *et al.* 2019). In the present study, in addition to the scouring parameters, the ability of the model to appropriately reproduce experimental velocity profiles was evaluated. The findings of this study are valuable for practitioners devising mechanisms to mitigate bed scour induced by bottom jets.

### Experimental data

The experimental data of Faruque *et al.* (2006) were used to validate the numerical model. For this experiment, a flume with the width, depth, and length of respectively 1.1 m, 0.92 m, and 8 m was used. The experimental properties are summarized in Table 1. The wall jet is ejected through a nozzle with a square cross-section outlet of width *b _{0}* = 76 mm. The bed is made of uniform sand particles with mean particle diameter of

*d*

_{50}= 2.46 mm, density of

*ρ*= 2,650 kg/m

_{s}^{3}, and geometric standard deviation of

*σ*= 1.24. The horizontal sand bed surface, 325 mm deep and 3 m long, was positioned at the invert of the nozzle outlet. Figure 1 shows a schematic view of the experimental setup. Experimental data that were measured during the experiment included the dimensions and shape of the scour-ridge system, profile and plan view, and velocity profiles at three different sections, from the upstream end of the bed.

_{g}Nozzle width b_{0} mm
. | Jet expansion ratio (B/b_{0})
. | Velocity U_{0} (m/s)
. | Tailwater ratio (H/b_{0})
. | Grain size to nozzle width ratio (d_{50}/b_{0})
. | Froude number (F)
. | Densimetric Froude number (F_{0})
. | Test duration (h)
. |
---|---|---|---|---|---|---|---|

76 | 14.5 | 1.31 | 6 | 0.032 | 1.5 | 6.6 | 72 |

Nozzle width b_{0} mm
. | Jet expansion ratio (B/b_{0})
. | Velocity U_{0} (m/s)
. | Tailwater ratio (H/b_{0})
. | Grain size to nozzle width ratio (d_{50}/b_{0})
. | Froude number (F)
. | Densimetric Froude number (F_{0})
. | Test duration (h)
. |
---|---|---|---|---|---|---|---|

76 | 14.5 | 1.31 | 6 | 0.032 | 1.5 | 6.6 | 72 |

## NUMERICAL MODEL

FLOW3D is a widely used commercial CFD code in which non-linear, transient, and second-order Navier–Stokes and continuity equations are discretized using finite difference and finite volume approaches to solve fluid motions (Meselhe *et al.* 2012; Flow3D User Manual 2016). FLOW3D uses both the TruVOF, which is an improved form of the original volume of fluid (VOF) technique (Usta 2014) and the fractional area-volume obstacle representation (FAVOR) method (Hirt & Sicilian 1985) to determine the location of the free surface and obstacles (rigid boundaries such as bridges, gate, etc.), respectively. Model convergence is reached when both viscous and pressure iterations converge. The convergence criterion for the pressure iterations, EPSI, is the same as the one for the explicit viscous algorithm. The convergence criterion is evaluated automatically by FLOW-3D as a function of the time step at every time cycle. The maximum pressure residual is controlled to be smaller than for the sake of accuracy. The generalized minimum residual (GMRES) algorithm was used as the pressure solver option in this simulation. This option has been successfully applied by many researchers (e.g., Shahrokhi *et al.* 2013; Usta 2014; Zhang *et al.* 2018), and has been recommended by FLOW3D User Manual for a wide range of problems due to its accuracy and efficiency. Pressure solver scheme is implicit, while viscous stresses and advection are computed using an explicit scheme. Second order monotonicity preserving was selected as a momentum advection approximation as recommended by FLOW3D for studying swirling, free-surface flows. In this study, time step is automatically calculated by the software, and can be changed during the simulation to avoid numerical instability. As follows, model governing equations are briefly described in the form of tensor notations:

**and**

*u*_{i}**are velocity and position vectors,**

*x*_{i}*t*denotes the time, p is pressure, ρ is density and

**refers to the viscous stress tensor and can be expressed by: in which is molecular viscosity and is the strain-rate tensor which can be defined as follows:**

*t*_{ij}Note that ** s_{ij}** =

**, so that**

*s*_{ji}**=**

*t*_{ji}**for simple viscous fluids (Wilcox 2006).**

*t*_{ij}*φ*denotes the scalar quantity such as concentrations, temperatures, etc. is a volumetric source/sink term, and Γ is molecular diffusivity for property (Rodi

*et al.*2013; An & Julien 2014).

## TURBULENCE MODEL

*i,j,k*), is the eddy viscosity, is the Kronecker delta and

*k*is the turbulent kinetic energy (Filonovich 2015). In this study, three different two-equation turbulence closure schemes, including RNG,

*k-ɛ*and

*k-ω*models were employed to determine which one yields better performance for simulating the flow field obtained from experimental data. Moreover, these turbulence models were investigated to evaluate their performance in predicting bed elevation changes in profile view caused by a wall jet. Two-equation turbulence models use eddy viscosity hypothesis and are the most widely used for engineering applications among the turbulence models. These models provide two additional differential equations to calculate turbulence length and time scales. In addition, they can be used to investigate the properties of a turbulent flow without previous knowledge of the flow structure or its geometry.

*k-ɛ* model

This model is a semi-empirical model and consists of two transport equations. One for the specific turbulent kinetic energy (*k*) and one for the turbulent dissipation rate ()*.* Equations and closure coefficients for this model are summarized below (Patursson 2008; Kositgittiwong *et al.* 2013):

*μ*is the dynamic viscosity, is the generation of turbulent kinetic energy caused by the average velocity gradient, is the generation of turbulent kinetic energy caused by buoyancy, and are the turbulent Prandtl numbers for

*k*and , respectively;

*S*and

_{k}*S*are source terms; and

_{ɛ}*Y*represents the contribution of the fluctuating dilatation in turbulence. The constants are chosen based on the classical model of Launder & Spalding (1974):

_{M}### Renormalized group (RNG) model

*k-*model are derived explicitly in the RNG model. This model is different from the standard

*k*- model in terms of constants and some additional terms and functions in the transport equations for

*k*and

*.*The RNG

*k-*model has several advantages over the standard

*k-*model. It is more accurate for rapidly strained flows and swirling flows and for lower Reynolds numbers (

*Re*). The transport equations for

*k*and in this model are given as follows (Patursson 2008; Kositgittiwong

*et al.*2013): where , , , and are as previously described. and are inverse effective Prandtl numbers for k and , respectively. is the effective viscosity = . The term is given as = . The constant values for this model are (Yakhot & Orszag 1986):

### The standard k-*ω* model

*k-ω*model is a semi-empirical model like the standard

*k-*model. It is based on the model proposed by Wilcox (1998). It has the same definition for

*k*as in the

*k-*model. However, it employs the specific dissipation rate,

*ω,*instead of turbulent energy dissipation

*,*which can be written as the ratio of to

*k*. The eddy viscosity is calculated as . The turbulence kinetic energy,

*k*, and the specific dissipation rate,

*ω*, are obtained from the following transport equations (Patursson 2008; Kositgittiwong

*et al.*2013): where is the generation of

*ω*; and are turbulent dissipation terms for

*k*and

*ω*; and is the source term for

*ω*. The model constants given by Wilcox (1998) are:

. More details can be found in Wilcox (1998).

It also should be noted that the above models constants have been used with success by a large number of studies (e.g., Savage & Johnson 2001; Johnson & Savage 2006; Khosronejad 2010; Li *et al.* 2011; Sharifipour *et al.* 2015; Rezaei & Safarzade 2016; Rezaei & Amiri 2018).

## SEDIMENT SCOUR MODEL

In the FLOW3D model, the calculation of sediment movement includes erosion, advection, and deposition processes. The total sediment transport rate has two states of bed-load and suspended load rates. Suspended sediments usually have low concentrations and are transported by fluid flow. Packed sediments can be defined by the user through a maximum packing fraction as a solid component (equal to 1 − porosity). The surface of packed sediment particles can be moved in the form of bed-load transport through rolling, hopping, or sliding along the bed (Wei *et al.* 2014; FLOW3D User Manual 2016).

**is the concentration of the suspended sediment, expressed in terms of sediment mass per volume of fluid-sediment mixture,**

*C*_{s,i}*D*is the directional diffusion coefficient, and

*u*is the suspended sediment velocity. Sediment is entrained by picking up and re-suspension due to local shearing and small eddies at the packed sediment–fluid interface (FLOW3D User Manual 2016). In FLOW3D, an empirical bed erosion model is used based on Mastbergen & Von den Berg (2003).

_{s,i}For non-cohesive sediments, particle entrainment occurs when the Shields parameter exceeds its critical value (Crosato 2008). The critical Shields parameter can be defined by the user and in the present study it was considered equal to 0.04 as estimated by the Soulsby–Whitehouse equation (Soulsby 1997). In the presence of a sloping bed, the critical Shields parameter can be corrected for the angle of repose by another relationship proposed by Soulsby (1997). The local Shields parameter is calculated based on the local bed shear stress, *τ*. A standard wall function is employed to determine the bed shear stress, taking into account the bed surface roughness. The bed surface roughness is estimated by Nikuradse's roughness *k _{s}* =

*C*

_{rough}* d_{50}, where

*C*and

_{rough}*d*

_{50}are a user-defined coefficient and the median grain diameter of the bed material, respectively. In the present study, a user-defined constant coefficient of

*C*= 2.5 is considered (Wei

_{rough}*et al.*2014).

The processes of entrainment and deposition of sediment particles act against each other and often occur at the same time. They are combined to obtain the net exchange rate between packed and suspended sediments. The entrainment lift velocity that causes the particle to leave the packed bed and to become suspended into the flow can be calculated based on the expression presented in Mastbergen & Von den Berg (2003). Subsequently, deposition refers to the tendency of suspended particles to fall under their own weight onto the packed bed and stop their movement as bed-load (Wei *et al.* 2014; FLOW3D User Manual 2016).

*q*. In this study, the Meyer-Peter–Müller equation was adopted because among the different models tested, it yields better agreement with the measurements, as discussed later. Moreover, the Meyer-Peter–Müller formula has successfully been used by many researchers (e.g., An

_{b,i}*et al.*2015; Castillo & Carrillo 2016, 2017; Esmaeili

*et al.*2017; Mirzaei

*et al.*2019). This formula is recommended for particles in a range between 0.4 mm and 28.65 mm according to Mohamad Noor

*et al.*(2016). The sediment particle diameter in the present study, which is 2.46 mm, is in the recommended range. The volumetric bed-load transport rate of sediment using Meyer

*-*Peter–Müller relationship can be expressed as: where

**is the volume fraction of species**

*C*_{b,i}*i*-in the bed material,

*is the bed-load coefficient with a default value of 8 (for low and very high transport rate of sand, the values of 5 and 13 have been recommended, respectively), is the magnitude of the acceleration due to gravity, is the mass density of the sediment particles, is the mass density of fluid, is sediment particle diameter, is the local shields parameter, and is the critical shields parameter (Wei*

_{MPM,i}*et al.*2014).

## NUMERICAL DOMAIN

The inlet boundary condition was set as specified velocity (V) equal to the experimental jet exit velocity. The boundary condition at the downstream end of the domain was described by a pressure boundary condition (P) corresponding to the constant water depth in the flume. The sidewalls and bottom boundaries were defined as a wall (W) boundary condition. The top boundary of the domain was determined as a symmetry (S) boundary condition to consider the atmospheric pressure on the free surface. Figure 2 shows the computational domain of the present study and associated boundary conditions. A 3D geometry of the erodible bed and wall jet are plotted in Figure 3.

## NUMERICAL SIMULATION

### Mesh sensitivity analysis

FLOW3D uses a simple grid of rectangular elements to discretize the computational domain. Using a structured mesh offers advantages such as fast mesh generation, maintaining numerical accuracy by changing the size of the elements slowly, reducing the computation time, and minimizing required output storage (Acharya 2011; FLOW3D User Manual 2016). For three different computation directions *x*(*i*), *y*(*j*), *z*(*k*), the mean values of flow parameters for the scalar quantities are calculated at the cell centers and for the vectors and tensors at the cell faces using a staggered grid technique at discrete times (Savage & Johnson 2001; FLOW3D User Manual 2016).

A grid sensitivity study was conducted in order to find the suitable cell size. Four different mesh sizes including 3 cm, 3.5 cm, 4 cm, and 4.5 cm were created. Figure 4 shows the variation of the mean relative error as a function of the cell sizes for vertical velocity profile at *x/b0* = 1.3 and *t* = 3,600 s. According to Figure 4, the simulated vertical velocity profiles exhibit better agreement with the measured velocities for the cell size of 1.5 cm. In addition, the variation of mean relative errors can be neglected by decreasing the cell size from 2 cm to 1.5 cm. Similarly, the variations of relative error of the maximum scour depth against cell sizes are presented in Figure 5. From this figure, it can be clearly seen that by decreasing the cell size from 4.5 cm, the relative error decreases and reaches a minimum value for cell size of 3 cm. Then, it increases and again reaches a maximum value at the cell size of 1.5 cm. A similar behavior in terms of mesh convergence is observed by Lee *et al.* (2019). Based on their results, the scouring parameters tend to be either overestimated or underestimated by changing the mesh sizes from the employed one. In summary, the cell size of 3 cm is selected in this study as an appropriate cell size because of two reasons. First, the relative error for the simulated maximum scour depths (as the most important parameter in scouring studies) decreases by reducing the cell size from 4.5 cm to 3 cm. It is also observed that for cell sizes range from 4 cm to 3 cm, the model outputs are not influenced drastically by the mesh size and are almost in the asymptotic range, as shown in Figure 5. Second, the mean relative error of the longitudinal velocity component is acceptable for the cell size of 3 cm. It is important to note that in Khosronejad & Rennine's (2010) study, the maximum relative error of longitudinal square jet velocity in simulating of wall jet on fixed bed reaches 100%. Hence, one can conclude that in this study, the error associated with the numerical simulation of the wall jet in the presence of live bed is acceptable.

## NUMERICAL RESULTS

The simulated results were compared against laboratory test data to validate the numerical model. Equilibrium morphodynamic conditions correspond to the time after which evolution of the scour hole ceased. Hence, the simulation end time was set to 23,000 s. Figure 6 shows the temporal development of the scour hole plotted as *Z/b _{0}* versus

*X/b*.

_{0}*Z*and

*X*represent the vertical and horizontal distances from the initial bed elevation and nozzle, respectively. Both distances are normalized by the nozzle width. The contour plots for the plan view of the simulated bed elevation are presented in Figure 7 at various instances of time. The dark black line in Figure 7 is the lowest bed elevation where the maximum scour hole occurred, and bright white is the highest bed elevation with maximum deposition. The results at different times indicate that the highest scouring propagation rate occurs during the early stages of the scouring process in the vertical, transversal, and longitudinal directions. Such behavior is in reasonable agreement with the work carried out by Lee

*et al.*(2019). It was observed that nearly 77

*%*of the maximum simulated equilibrium scour depth of 11.7 cm occurs during the first 2,000 s, as can be observed in Figures 6 and 7. The scour hole continues to gradually deepen until reaching near equilibrium at

*t*

*=*12,000 s. Also, it was noticed that almost 88% of the maximum erosion in the transversal direction takes place within the initial 2,000 s while reaching nearly equilibrium conditions at

*t*= 4,000 s, which can be inferred from the sequence of the contours shown in Figure 7,

*.*Examining the erosion in the longitudinal direction showed that approximately 67% of the maximum equilibrium scour length has occurred by

*t*

*=*2,000 s. Afterwards, the eroded sediments tend to move downstream as time progresses until the incoming jet is not capable of pushing the eroded sediments any farther away. As shown in Figure 6, after

*t*

*=*12,000 s, only the downstream end of the scour hole continues to evolve, although the evolution is very slow. The latter is reflected by the subsequent increase in scour length as can be seen in Figure 7. The ridge formed just downstream of the scour hole forms mainly within

*t*= 12,000 s, as the sediment removed from the scour hole is gradually accumulated at the ridge. The ridge and scour hole growth rates decelerate at the same time. The decrease in growth rate of the ridge is due to the fact that the mean flow velocity reduces in the scour hole as time goes on. This is because the growth of scouring size increases the cross-sectional area, eventually leading to reduction in the sediment transport capacity.

## MEAN FLOW VELOCITY FIELD

The ability of the numerical model to accurately reproduce the experimental scour hole and ridge patterns relies first on how well the simulated flow field compares against the experiments. Figure 8(a)–8(c) shows comparison between measured and simulated vertical profiles of velocity magnitude as a function of the non-dimensional vertical distance *z/b _{0}* at three sections from the nozzle outlet. Here,

*z*and

*b*are the distance from the bed surface in

_{0}*z*-direction and the nozzle width, respectively. As can be observed in Figure 8(a), there is a uniform velocity distribution close to the nozzle outlet section. The reverse flow, which is formed due to the flow separation at the top of the deposited sediments, also affects the velocity distribution. The effects of the reverse flow increase not only near the location of the flow separation but also with vertical distance from the nozzle centerline where the positive water velocity is maximum. This explains the larger values of negative velocities in Figure 8(c) in comparison with Figure 8(b) and 8(a). Figure 8(c) clearly shows that the negative flow velocities increase closer to the water surface. Due to the higher importance of near-bed velocities in the scouring process, these velocities were considered as benchmark data for calculating the relative errors between modeled and measured vertical profiles. The difference of maximum velocities between modeled and measured results has a relative error of 3.3% for Figure 8(a), 22.3% for Figure 8(b), and 10.8% for Figure 8(c). It is important to note that the numerical model adequately follows the general trend of measured vertical velocity profiles. However, the obtained discrepancies between the measured and numerical results can be due to the complex interactions between flow field

*,*sediment transport process, and sediment bed morphology that were previously pointed out in this study.

Figure 9 shows the simulated velocity vector field in the *xz-*plane along the centerline of the experimental arrangement, i.e., at *y* = 0. Vertical distance of the jet body from the sediment bed increases as the sediment bed elevation decreases during the scouring process. Hence, the jet is neither a wall jet nor a free jet in the presence of the scour hole. It is believed that although the entrainment of water from the above and below is not symmetric, the jet characteristics depart from the wall jet and enhance the free jet. Note that due to relatively large particle sizes, i.e., *d*_{50} = 2.46 mm, sediment transport process is mostly carried out as bed-load, while a very small amount is transported as the suspended load.

The upstream slope of the ridge increases by accumulation of sediment over time. This results in the flow separation at the top of the ridge as the flow cannot follow the bed surface any longer due to the sudden pressure drop past the crest of the ridge. Among the turbulence models studied, further examinations revealed that only the *k-ω* turbulence closure scheme is able to reproduce the observed flow separation and reverse flow. Therefore, it was employed in this study. The reverse flow shown in Figure 9 is caused by counter-clockwise vorticity that is formed due to flow separation. This reverse flow eventually joins the normal jet flow. After passing over the top of the ridge, the flow splits into two components. One relatively larger part moves toward the water surface while the other advances toward the bed. This flow distribution causes the formation of several other vortices. The unchanged bed elevation downstream of the ridge is likely due to the fact that each pair of symmetric vortices formed along the flume are rotating in opposite directions, as can be seen in Figure 9.

Figure 10(a)–10(c) indicate the comparison of performance accuracy of different turbulence models in simulation of vertical profiles of velocity magnitude at three different sections. The x-axis and y-axis labels are previously described. Three common two-equation turbulence models, namely: the *k-*, RNG, and were used. From this figure, it can be easily seen that the *k-* and RNG turbulence schemes fail in predicting the negative measured velocities. The result from the *k-ω* model demonstrates that this model is capable of capturing the reverse flow caused by flow separation, indicating that the *k-ω* turbulence model has a superior performance for flows with adverse pressure gradients than *k-* models (Wilcox 1998; Bates *et al.* 2005; Tu *et al.* 2005).

## EVALUATION OF TURBULENCE MODELS IN PREDICTING SCOUR PARAMETERS

Three different turbulence models are assessed to suggest the best numerical turbulent modeling scheme for simulating scour under submerged wall jet (Figure 11). As can be seen, there is a remarkable difference among scour hole lengths obtained from different turbulence models compared to the measured ones. In contrast to the *k-ω* model, the RNG and *k-* models predict the scour hole length with sufficient accuracy. The relative error for the RNG and *k-* models are 8.8% and 17%, respectively. However, these models give inaccurate results in predicting maximum of ridge height (the *k-* and RNG turbulence models underestimate the ridge height by 39% and 36%, respectively, compared with the measurements). Regarding scour depth estimation, the best fit to experimental data is obtained by the *k-ω* model predicting the maximum scour depth with a relative underestimation of 11.6% compared to measured data, while the RNG and *k-* models underestimate the scour depth by 19.5% and 20%, respectively. As a result, the *k-ω* turbulence model is chosen in this study due to better prediction of scour depth and ridge height as well as its ability to capture the negative measured velocities.

## SCOUR DEPTH AND RIDGE HEIGHT

Figure 12 compares experimental and simulated bed profiles at *t* = 23,000 s. The three simulated bed profiles were obtained using the bed-load equations of Meyer*-*Peter–Müller, Nielsen, and Van Rijn. In Figure 12, *Z* and *X* are the distance from the bed surface in the z-direction and longitudinal distance from the nozzle, respectively. The x-axis and z-axis are normalized with respect to the nozzle width. As can be seen, the equations of Meyer*-*Peter–Müller and Van Rijn result in very similar scour bed profiles. However, the Meyer*-*Peter–Müller equation exhibits better agreement with the measurements in terms of maximum scour depth and ridge height. Furthermore, Meyer*-*Peter–Müller equation appropriately reproduces the experimental scour hole shape*.* Comparing with the measurements, the maximum computed scour depth and ridge height by Meyer*-*Peter*–*Müller formula are satisfactorily underestimated with relative errors of 11.6% and 2.3%, respectively. Regarding the Van Rijn formula, the maximums of the relative error for scour depth and ridge height are 24.2% to 3.2% showing underestimation and overestimation compared to the observed data, respectively. The Nielsen's relation underestimates the scour depth by 11.6% and overestimates the ridge height by 10.5% with respect to the measured ones. The main difference can be observed between the measured and calculated scour hole lengths. It is also observed that the scour hole length predicted by Van Rijn equation is slightly more accurate than the Meyer*-*Peter–Müller formula. Regarding Nielsen's equation, the corresponding graph indicates a significant overestimation of scour hole length. Moreover, this equation does not capture the shape of scour hole appropriately. Hence, the Meyer*-*Peter–Müller equation is considered for bed-load calculations in this study.

A more quantitative comparison of bed elevation changes between modeled and measured data at multiple instances of time is shown in Figure 13. As illustrated in this figure, the modeled scour hole length is 58% larger than the measured one during the first 3,600 s. From *t* = 3,600 s to *t* = 10,800 s, the changes in simulated scour hole length is relatively large. Beyond this time, there is a slight change between scour hole lengths due to approaching equilibrium conditions. Finally, the predicted scour hole length is 60% larger than the one observed at *t* = equilibrium time. Overall, it can be inferred that the predicted scour hole length tends to be overestimated significantly compared to the experimental findings. A similar discrepancy between the simulation and measurements was previously reported by Abdelaziz *et al.* (2010). It can therefore be concluded that the model performance in reproduction of the spreading coefficient of the wall jet is poor.

## SCOUR HOLE AND RIDGE PERIMETER

Figure 14 depicts the plan view of the perimeter of the scour hole and the ridge obtained from simulations using the bed-load equations of Meyer*-*Peter–Müller, Nielsen, and Van Rijn. The y- and x-coordinates are transverse and longitudinal distances from the nozzle centerline and the nozzle location, respectively. Both distances are normalized by the nozzle width. The jet ejected from the nozzle interacts with the sediment bed; subsequently, the jet expands and scouring occurs in the transversal direction symmetrically on both sides of the jet. It should be noted that the scour hole length and the ridge height increase until the equilibrium state is reached. The simulated results from the FLOW3D model are in reasonable agreement with the experimental observations, as discussed later in this section. Particularly, the overall planar shape of the scour hole and ridge along their corresponding maximum values are relatively well estimated by the model. The model overestimates the scour hole length, leading to a hole and ridge characterized by milder longitudinal slopes than the experimental ones. Referring to Figure 14, it is observed that the three different bed-load equations estimate the maximum perimeter of scour hole with similar accuracies (a relative error of 18.2%) against the measured one. From this figure, it can be observed that the predicted results of three different bed-load equations are different in terms of maximum perimeter of ridge. The relative errors in predicting maximum perimeter of ridge are 13.3% for Meyer*-*Peter–Müller formula, 7.2% for Van Rijn formula, and 1.2% for Nielsen formula.

The perimeter of scour hole and ridge as simulated using Meyer*-*Peter–Müller formula at multiple instances of time can be observed in Figure 15. During the first 3,600 s, the relative errors in predicting maximum perimeter of scour hole and ridge with respect to the measured ones show overestimations and underestimations by 5.5% and 12.9%, respectively. It is to be noted that the perimeter of the scour hole and ridge reach nearly equilibrium conditions at *t* = 4,000 s and by approaching the final simulation time (*t* = 21,000 s), respectively.

## CONCLUSIONS

Three main objectives were considered in implementing the present study: first, to test the ability of the FLOW3D model to reproduce experimental flow velocities in the scour hole induced by a submerged wall jet; second, to evaluate the ability of FLOW3D to accurately reproduce scour hole and ridge formation and evolution; and, third, to assess the accuracy of two equation turbulence models in estimating the velocity field as well as scour parameters. The main findings of this study can be summarized as follows.

This study revealed that FLOW3D is capable of reasonably reproducing measured flow velocities. Three different turbulence models were employed to investigate the impact of turbulence models on flow field properties. In general, all of the turbulence models used in this study adequately follow the trend of measured vertical velocity profiles. However, in contrast to the *k-* and RNG turbulence models, the *k-ω* turbulence model successfully simulated the reverse flow (at farther distances from the bed) caused by flow separation at the top of the deposited sediments.

Investigations indicated that among the turbulence models studied, the k-*ω* model presents the best agreement with the observed scour depth and ridge height. The scour depth and ridge height were underestimated by 19.5% and 36% respectively using the RNG turbulence model. The underestimation was 20% and 39% for the *k-* turbulence model and 11.6% and 2.3% for the *k-ω* model, showing higher accuracies for the *k-ω* model. However, it was found that the *k-* and RNG turbulence models provide more accurate predictions of scour hole length than the *k-ω* model.

Three different bed-load equations of Meyer*-*Peter–Müller, Nielsen, and Van Rijn were evaluated in this study. Comparing with the measurements, the results showed that maximum predicted scour depth and ridge height by Meyer*-*Peter–Müller formula are underestimated by 11.6% and 2.3%, respectively. Regarding Van Rijn formula, the maximums of scour depth and ridge height are underestimated by 24.2% and overestimated by 3.2%. The Nielsen's relation underestimates the scour depth by 11.6% and overestimates the ridge height by 10.5%. However, all formulas significantly overestimated the length of the scour hole. Also, it was found that the Meyer*-*Peter–Müller equation reproduces the experimental scour hole shape better than the other two equations.

The simulated perimeters of scour hole and ridge, i.e., at *y**=* 0 using three different bed-load equations were also compared to their experimental counterparts. In general, the maximum simulated perimeter of scour hole and ridge demonstrated good agreement against experimental data. All formulas predicted the maximum perimeter of scour hole with almost the same accuracy, while the results indicated difference performances in predicting the maximum perimeter of the ridge. The relative errors in predicting the maximum perimeter of ridge are 13.3% for Meyer*-*Peter–Müller, 7.2% for Van Rijn, and 1.2% for Nielsen formulation.

The temporal evolution of bed erosion and deposition indicates that a significant rate of scour development occurs during the early stages of the time evolution in all directions. Additionally, it was observed that the height of the ridge and the length of the scour hole keep increasing continuously until *t**=* 23,000 s. After this time, the bed morphological changes are negligible. Nevertheless, the erosion in the vertical and lateral directions reach equilibrium state at about *t**=* 12,000 s and *t**=* 4,000 s, respectively.

The simulation results presented in this paper showed that Flow3D model can appropriately reproduce planar and vertical geometry of the scour hole and ridge, while the length of the scour hole is significantly overestimated. The reason for the lower accuracy in prediction of the scour hole length could be the fact that the model is not capable of properly estimating the spreading coefficient of the wall jet. Overall, the model results provide estimations of the different scour parameters with acceptable accuracy. This can be exceedingly valuable from an engineering standpoint to assist in the evaluation of alternative scenarios to mitigate jet-generated scouring downstream of dams during planning and design phases.

## REFERENCES

**FLOW-3D**