Performance Evaluation of Different Stilling Basins Downstream of Barrage Using FLOW-3D Scour Models

Abstract

The present study investigated the performance of three different stilling basins, i.e., modified United State Bureau of Reclamation (USBR) Type III, USBR Type II, and wedge-shaped baffle blocks (WSBB), using FLOW-3D scour models. Field data of the riverbed profile are employed to validate the present models. After comparison, the results of statistical indices, i.e., coefficient of determination (R²) and Nash–Sutcliffe model efficiency coefficient (NSE), indicated that the Renormalization Group (RNG-K-ϵ) showed good agreement with the field data, with R2 and NSE values of 0.9094 and 0.896, respectively. Validated models are used to simulate velocity field and local bed shear stress (BSS) and scour for design and flood discharges of 28.30 m³/s/m and 17.5 m³/s/m, respectively. At 28.30 m³/s/m, the results indicated that the riverbed downstream of the remodeled basin was completely exposed, while, at 17.5 m³/s/m, the net change in bed reached 85%. At 28.30 m³/s/m, the net change at the centerline of models reached 51% and 67% in USBR Type III and WSBB basins, respectively. At 17.5 m³/s/m, compared to Type II and III basins, the WSBB basin indicated less BSS, which significantly reduced the scour. Conclusively, the Type II basin showed less energy dissipation for the studied flows, while the WSBB basin improved flow fields downstream of the barrage.

1. Introduction

1.1. Hydraulics of Old and Remodeled Barrages

Barrages in the plain regions of Pakistan were built about 50 to 100 years ago, which provide the required amount of water. However, due to hydraulic and structural deficiencies, some of the barrages were remodeled, i.e., the Taunsa barrage, Jinnah barrage, and Khanki barrage [1,2]. The old stilling basins of these barrages resembled the United States Bureau of Reclamation (USBR Type III) basin, which included energy dissipators, i.e., friction and baffle trblocks. These basins stabilize the hydraulic jump (HJ) and dissipate the kinetic energy of the upstream supercritical flow. In contrast, under the remodeling process, the modified USBR Type III basins are transformed into USBR Type II basins, in which the old basin’s appurtenances are demolished and replaced with chute blocks and end sill [3,4].

The Taunsa barrage was built on the river Indus in 1958 with a design flow of 24.30 m³/s/m. Soon after the first operation, several issues, i.e., baffle block uprooting, lowering of tailwater level, and excessive scouring, were noticed on the barrage downstream. Until the end of 2003, the aforementioned issues were addressed with minor repairs and, from 2004 to 2008, based on a model study, the barrage basin was remodeled [5]. The model study only investigated the locations of hydraulic jumps (HJs) with different tailwater depths. However, even after the remodeling, in 2010, the barrage was unable to discharge the super flood of 23.22 m³/s/m, which breached its left marginal bund, and an uncontrolled flood of 7000 m³/s drowned the Muzaffargarh district. The riverbed data for 2010 revealed large scour pits on the barrage downstream. In addition, in front of some of the bays, the block filter floor was also found to have disappeared. Furthermore, it is also reported that, in high discharge periods, the remodeled basin also raises water levels upstream, which develops operational issues during flood times [6].

1.2. Stilling Basins and Local Scour

Stilling basins are constructed downstream of the hydraulic structures, i.e., dams, weirs, spillways, and barrages, to dissipate the potential and kinetic energy of the upstream flow [7]. These basins are of various shapes and sizes and contain different structural arrangements, such as chutes, baffle blocks, friction blocks, and end sills. Primarily, such arrangements stabilize the HJs and reduce the overall length of basins. USBR developed Type II and Type III basins to design high- and low-head hydraulic structures, which include different shapes of appurtenances [8,9]. Al Mansori et al. [7], Ali [10], Chaudary [11], and Macián-Pérez et al. [12] employed these basins and studied their effects on the HJ, free surface profiles, sequent depth, turbulent kinetic energy, and energy dissipation. Similarly, Habibzadeh et al. [13,14], Eloubaidy et al. [15], Tiwari et al. [16], Tiwari et al. [17], Widyastuti [18], Goel [19], Goel [20], and Pillai et al. [21] carried out experiments to investigate the impacts of energy dissipators on the HJ and energy dissipation downstream of open channel and pipe outlet basins. However, the above-mentioned studies were carried out on rigid beds, which failed to provide information about erodible beds. In the real world, the performance of stilling basins is usually assessed by the extent of erosion and retrogression downstream. Therefore, due to the significant importance of these hydraulic structures, the study of local scouring is an important research area.

Scouring is a process of sediment removal from the downstream end of graded control structures such as weirs, barrages, and spillways [22,23]. This process starts from the downstream end of the stilling basin when bed shear stresses increase beyond the critical limits. To reduce the scour, ample protection measures, i.e., rigid floor with dissipator, inverted filter blocks, flexible apron, and loose stone, are deployed downstream of hydraulic structures [24]. Balachandar et al. [25], Mohammed et al. [26], Wuthrich et al. [27], Elsayed et al. [28], and Amin [29] conducted experiments on different shapes of stilling basins and energy dissipators to investigate flow behavior and local scour. In these studies, velocity field, scour depth, and scour length were the main investigated parameters that were examined under different hydraulic conditions, i.e., Froude number, pond level, tailwater, gate openings, and discharges.

1.3. Computational Fluid Dynamic and Hydraulic Modeling

In the former Section 1.1 and Section 1.2, the significance and issues related to the different basins are highlighted. These studies are carried out in laboratories and on reduced scale models, which increases both cost and time. In addition, in physical modeling, due to the scaling effects [30], considerable differences in the outputs are found to occur between the prototype and laboratory models. Furthermore, measuring devices also create a hindrance in the flow. Similarly, because of the differences in terrain and concrete roughness, the model result deviates from the prototypes. However, due to the advancement in computer technology and turbulence models [31], the hydraulic and scour investigations on the graded control structures can be assisted by numerical models [32,33]. Upon use of different turbulence schemes, i.e., two equation turbulence models (Renormalization Group (RNG K-ε) [34] and Standard K-ε [35]), and large eddy simulation (LES) [36], the numerical modeling tools such as FLOW-3D [37], Fluent CFD [38], and OpenFOAM [12] are considered to be more suitable for the investigations of HJ, energy dissipation, and local scouring. Jalal and Hassan [39], Alasta et al. [40], Mehnifard et al. [41], Samma et al. [42], and Epely-Chauvin et al. [43] applied FLOW-3D to study the local scour downstream of the open-channel hydraulic structures. In these studies, the effects of pond level, gate openings, tailwater depths, gate operation, and flow variations were investigated, and the results highlighted the implications of these parameters on length and depth of scour.

1.4. Research Motives and Problem Statement

The old stilling basin of Taunsa barrage consisted of two rows of baffle and friction blocks. Such types of basins are similar to USBR Type III basin, which stabilize the HJs and increase the energy dissipation. During the year 2008, based on the limited hydraulic study on a scaled model of 1:50, the barrage basin was remodeled, and it was considered that the old basin produced retrogression, which consequently reduced the required tailwater depths for the HJs. In remodeling, the friction and baffle blocks were replaced with chute blocks and end sill. In addition, except tailwater, the model study did not investigate the effects of new basin on the downstream river bed [3,4]. Consequently, after the 2010 flood, several issues were reported on the barrage downstream such as launching of block floor apron, large scour pits, and reduction in the discharging capacity of barrage [6].

On the other hand, after the literature review, it is found that the wedge-shaped baffle blocks (WSBB) are only investigated on the downstream of pipe outlets; thereby, their effects downstream of river diversion barrage need to be explored. In view of the above-mentioned, the present study has developed the FLOW-3D scour models for USBR-Type II, Type III, and WSBB basins. The study investigates velocity field, bed shear stress, scour depth, and scour length at design and flood discharges. The performance of the present scour models is assessed with riverbed profiles downstream of USBR Type II basin (existing basin). After the validation, the results of investigated scour parameters are utilized to draw comparison among the studied basins.

2. Materials and Methods

2.1. Study Area and Energy Dissipating Arrangements

Taunsa barrage is an important diversion structure built across mighty river Indus to supply the required amount of water for the arid zone of Southern Punjab. The barrage was constructed from the year 1954 to 1958. The barrage is located about 18 Km away from Kot Addu Town. For passage of flow, a total of 1177 m is available, while 1325 m is the total width between the barrage’s abutments. For energy dissipation, the old barrage’s basin includes USBR baffle and friction blocks [1,2]. The barrage is more than a traditional diversion structure because it accommodates transportation infrastructure, i.e., railroad crossing, an oil pipeline, gas pipeline, transmission lines, and also serves as an arterial road bridge. Furthermore, the barrage also supplies uninterrupted flows to four canals that take off from its upstream, i.e., Muzaffargarh (command area (2)), DG Khan (Command area (1) and (3)), and Kachhi canal. The Taunsa Punjnad (TP) Link canal (command area (4)) supplies water to the Chenab River to meet the water requirements of command area (5) of Panjnad Barrage [3,4]. The general layout of Taunsa barrage irrigation system is presented in Figure 1.

From the year 1958 to 2003, the old barrage’s basin (hereafter, Type I basin) faced multiple issues, i.e., uprooting of baffle blocks, excessive retrogression, malfunctioning of mechanical gates, and the sweeping of HJ. During these years, the above-mentioned issues were addressed by the partial repairs, which continued to aggravate, and, from the year 2004 to 2008, the barrage was remodeled (hereafter, Type II basin). In the remodeling, the downstream basin’s floor level was raised about 0.31 m and the old basin’s appurtenances were replaced with the chute blocks and end sill [1,2]. To study the flow dynamics on Type I and II basins, the scour models are developed using FLOW-3D. Additionally, to address the uprooting problems of vertical face baffle blocks, this study also employs WSBB (hereafter, Type III basin) on the downstream of Taunsa barrage. For the present study, the WSBB with vertex angle of 150° and cutback angle 90° are placed in the baffle blocks region. Three-dimensional (3D) illustrations of Type I, II, and III basins are presented in Figure 2.

2.2. Numerical Model Implementation

Presently, three different energy dissipating systems are investigated downstream of the studied barrage, and their performances are assessed by FLOW-3D scour models. Before employing the scour models, the hydraulic models are run to confirm the suitability of meshing, turbulence models, and initial and boundary conditions. At present, FLOW-3D is considered an effective software for hydraulic and scour investigations that employs structured rectangular grids to solve the solid and fluid domains. In FLOW-3D, the fluid’s properties are carried out by the finite volume method (FVM). For three-dimensional turbulent flow, Reynold average Navier–Stokes Equation (RANS) models are employed in each computation grid to discretize the Navier–Stokes Equation, i.e., mass, continuity, and momentum [32,33]. For incompressible flow of constant density $\mathsf{\rho }$, the present model employs the mass continuity Equation (1), while in x, y, and z directions, the momentum Equations (2a), (2b), and (2c) are utilized, respectively.

$$\frac{\partial }{\partial x}\left(\mathrm{u}{A}_x\right)+R\frac{\partial }{\partial y}\left(\mathrm{v}{A}_y\right)+\frac{\partial }{\partial z}\left(\mathrm{w}{A}_z\right)+\xi \frac{\partial u{\mathrm{A}}_x}{\mathrm{x}}\left(u{A}_x\right)=\frac{{\mathrm{R}}_{SOR}}{\rho }$$

(1)

$$\frac{\partial u}{\partial t}+\frac{1}{{\mathrm{V}}_{\mathrm{F}}}\left[\mathrm{u}{\mathrm{A}}_{\mathrm{x}}\frac{\partial u}{\partial x}+\mathrm{v}{\mathrm{A}}_{\mathrm{y}}\frac{\partial u}{\partial y}+\mathrm{w}{\mathrm{A}}_{\mathrm{z}}\frac{\partial u}{\partial z}\right]=\frac{1}{\rho }\frac{\partial p}{\partial x}+{\mathrm{G}}_x+{\mathrm{f}}_x$$

(2a)

$$\frac{\partial v}{\partial t}+\frac{1}{{\mathrm{V}}_{\mathrm{F}}}\left[\mathrm{u}{\mathrm{A}}_{\mathrm{x}}\frac{\partial u}{\partial x}+\mathrm{v}{\mathrm{A}}_{\mathrm{y}}\frac{\partial u}{\partial y}+\mathrm{w}{\mathrm{A}}_{\mathrm{z}}\frac{\partial u}{\partial z}\right]=\frac{1}{\rho }\frac{\partial p}{\partial y}+{\mathrm{G}}_y+{\mathrm{f}}_y$$

(2b)

$$\frac{\partial w}{\partial t}+\frac{1}{{\mathrm{V}}_{\mathrm{F}}}\left[\mathrm{u}{\mathrm{A}}_{\mathrm{x}}\frac{\partial u}{\partial x}+\mathrm{v}{\mathrm{A}}_{\mathrm{y}}\frac{\partial u}{\partial y}+\mathrm{w}{\mathrm{A}}_{\mathrm{z}}\frac{\partial u}{\partial z}\right]=\frac{1}{\rho }\frac{\partial p}{\partial z}+{\mathrm{G}}_z+{\mathrm{f}}_z$$

(2c)

where in Equations (1), (2a), (2b), and (2c), u, v, and w are velocity values; A_x, A_y, and A_z are the flow areas; G_x, G_y, and G_z are body accelerations; and f_x, f_y, and f_z are viscous accelerations in x, y, and z directions, respectively; R_SOR, P, and R are mass source, pressure, and model coefficient, respectively.

Turbulence Modeling and Free Surface Tracking

The turbulent flows are associated with the fluctuations (turbulence) of velocity components. In engineering problems, due to the small scale and high frequencies, these fluctuations cannot be solved directly. Therefore, to solve these issues, RANS models are employed in numerical modeling; however, these models generate additional unknowns which can be solved by turbulence models, i.e., Standard K-ε, RNG K-ε, and LES. Out of these models, the two equation turbulence models (K-ε) are extensively applied in engineering applications [32]. The studies [39,44] reported that the RNG K-ε model results are in good agreement with the experimental data. In addition, these two models are also studied for the investigations of HJ characteristics, velocity profiles, bed profiles, and bed shear stresses.

Table 1. RANS model application in hydraulic and energy dissipation.

Reference	3D Software	Turbulence Models	Parameters Studied
Carvalho et al. [33]	FLOW-3D	RNG K-ε	Velocity, roller lengths, free surface profile, TKEs
Bayon et al. [45]	Open Foam, FLOW-3D	RNG K-ε	Velocity, HJ efficiency, roller length, sequent depths
Aydogdu et al. [46]	ANSYS-FLUENT	RNG K-ε	Free surface profile velocity, turbulence kinetic energy (TKEs)
Abd El Azim et al. [47]	FLOW-3D	RNG K-ε	Energy dissipation, relative energy loss
Kosaj et al. [48]	FLOW-3D	RNG K-ε	Scour, velocity, channel bed slope
Mirzaei and Tootoonchi [49]	FLOW-3D	Standard K-ε, LES, RNG K-ε,	Free surface profile, sequent depths, roller, and lengths of HJ
Macián-Pérez et al. [50]	FLOW-3D	Standard K-ε, K-ω, RNG K-ε	Free surface profiles, velocity
Karim and Ali [51]	FLUENT	RSM, Standard K-ε, RNG K-ε	Velocity, bed profiles, bed shear stress

highlights a few of the relevant studies which have employed two equation RANS models.

Hence, based on the results of the above-mentioned studies, this study also employed RNG K-ε model (hereafter, called K-ε) for the turbulence and flow fields on rigid and erodible beds. For computing the turbulent kinetic energy (K) and its dissipation (ɛ), RNG K-ε model employs Equations (3) and (4), respectively.

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho kui\right)=\frac{\partial }{\partial x_j}\left[\mathrm{\mu }+\frac{\mathrm{\mu }t}{\mathsf{\sigma }\mathrm{k}}\left(\rho k\right)\frac{\partial k}{\partial xj}\right]+Pk+\rho \mathrm{\varepsilon }$$

(3)

$$\frac{\partial }{\partial t}\left(\rho \mathrm{\varepsilon }\right)+\frac{\partial }{\partial x_i}\left(\rho \mathrm{\varepsilon }{u}_i\right)=\frac{\partial }{\partial x_j}\left[\mathrm{\mu }+\frac{\mathrm{\mu }t}{{\mathsf{\sigma }}_{\mathrm{\varepsilon }}}\left(\rho k\right)\frac{\partial \mathrm{\varepsilon }}{\partial x_j}\right]+C_{1\epsilon }\frac{\mathrm{\varepsilon }}{k}Pk-C_{2\epsilon }\frac{{\mathrm{\varepsilon }}^2}{k}$$

(4)

where $x_i$, $\mathrm{\mu }$, ${\mathrm{\mu }}_t$, k, and $\mathrm{\varepsilon }$ are the co-ordinate in x direction, dynamic viscosity, turbulent dynamic viscosity, turbulent kinetic energy (TKE), and turbulent dissipation, respectively, whereas $\rho$, $P_k$, $\mathsf{\sigma }\mathrm{k}$, ${\mathsf{\sigma }}_{\mathrm{\varepsilon }}$,$C_{1\epsilon }$, and $C_{2\epsilon }$ are the fluid density, production of TKE, and model parameters, respectively.

In contrast, researchers have also employed large eddy simulation (LES) model to investigate local scour in meandering channels and around the bridge piers. The scour results of these studies agreed well with the experimental data. However, the use of LES model on the downstream of graded control structures at higher discharge, i.e., 24.30 and 18 m³/s/m, is yet to be explored. Therefore, to check the effects of eddies and recirculation in the scoured and retrogressed zones, the study also employs LES turbulence model [52].

Presently, the volume of fluid (VOF) method is used to study the water surface profiles, in which a fraction of the fluid is employed to find the fractional volume (i.e., water or air). In the present models, clear water with free surface is employed, while the other advection schemes are automatically selected by the numerical models.

2.3. Sediment Scour Model

Suspended and bed loads are generally used for scour modeling, which include advection, erosion, settlement, and deposition of the sediment species. FLOW-3D scour models represent the behavior of suspended and packed sediments in three-dimensional flow. A number of essential parameters are applied to determine the scour model, i.e., critical shield number (CSN), bed load transport rate equation, maximum sediment packing fraction, sediment diameter, sediment density, and the bed shear stress.

The shear stress on the sediment bed is computed by CSN, which is the minimum or critical bed shear stress required to remove the sediment particles from the packed bed. The CSN mainly depends on the sediment size and density. FLOW-3D employs the Soulsby–Whitehouse equation to compute the CSN; however, the user can also apply prescribed values. In the present study, prescribed values are adopted and the models automatically selected the CSN as provided in

Table 2. Physical properties and essential parameters for scouring.

Sr. No.	Sediment Properties	Value Assigned
1	Bed Load species	Fine sand
2	Species diameters (d50)	0.0002 m
3	Sediment density	1692 kg/m3
4	CSN	0.05
5	Entrainment Coefficient (EC)	0.018
6	BLC	0.053
7	Angle of repose (Degree)	32°

.

As the current sediment bed is completely made of sand (noncohesive), thereby, for bed load transfer rate, Van Rijn Equation (5) [44,53] is applied.

$$\mathsf{\Phi }\mathrm{i}={\beta }_{VR,i} \mathrm{d}\ast ,\mathrm{i}-0.3(\frac{\theta i}{\mathsf{\theta }\mathrm{c}\mathrm{r},\mathrm{i}}-1.0)2.1 c_{b,I}$$

(5)

where ${\beta }_{VR,i}$ is the bed load coefficient (BLC); C_b,I is defined as the fractional volume of packed sediment bed of i species; and (θ_cr,i) denotes the critical shield parameter.

The maximum packing fraction is the volume of all sediments in a cell. The grain size analysis from the barrage authority indicated that the sediment bed at the barrage site is composed of sand particles which have mean diameter and density of 0.0002 m and 1692 kg/m³, respectively, as given in Table 2. Due to the uniform composition of sediment particles, i.e., sand, this study has employed only one species for the sediment packed bed and, for that, the maximum default fraction of 0.64 is utilized. Considering the roughness of wall, a standard wall function is employed to calculate the bed shear stress; the bed roughness/d₅₀ is set to 2.5. For the grain lifting velocity, the entrainment coefficient (EC) of 0.018 is employed, while, for the angle of repose, a default value of 32 is employed for the horizontal bed.

2.4. Models’ Layout, Meshing Setup, and Boundary Condition

The solid geometries of Type I, II, and III basins are prepared in AutoCAD and then converted into stereolithography (stl.) files. On the other hand, the erodible beds of 100% sand particles are created in FLOW-3D, which are placed on the downstream of the studied basins. The solid and erodible geometries for the present models are shown in Figure 3.

Two structural mesh blocks are employed to solve the solid and fluid domains. The first block covers solid geometries of the models, which ranges from X_min = 15 m to X_max = 71 m. In total, 1,492,000 mesh cells are utilized to resolve the basins’ geometries and appurtenances. On the other hand, the second block is extended from X_min = 71 m to X_max = 100 m, which also includes 30 m of erodible bed. A total of 222,720 mesh cells are employed for the sediment packed bed, as shown in Figure 4. It is worth mentioning that, in the first block, the HJ and velocity fields are closely assessed, while, in the second block, the retrogression and scour patterns are evaluated. For the first mesh block, nonuniformed cell sizes of 0.15, 0.20, and 0.15 m are applied in x, y, and z directions, respectively, while, for the second block, a uniform cell size of 0.25 m is employed in all the directions. The details of cell sizes are provided in Figure 4.

The FAVORized views of solid and fluid domains are shown in Figure 5. From Figure 5, it can be believed that the present meshing setup has reasonably resolved the baffle blocks and other regions, which allowed the simulations to be run.

For the first mesh block, pressure (P) boundaries are employed both on the upstream and downstream, while wall (W) boundaries are applied to the sides (Y_min, Y_max) and solid floor (Z_min). In contrast, in the second mesh block, symmetry (S) and pressure (P) boundaries are set on the upstream and downstream, respectively, while the rest of the boundaries are set similar to the first mesh block. For all the studied models, the upper boundaries are set to atmospheric pressure (Z_max), which allows the fluid to null von Neumann.

2.5. Models’ Operation and Initial Conditions

During a year, the flow over the diversion barrage changes from one bay to another. In these situations, most of the time, the barrage goes through two different flow conditions, i.e., controlled free (CF) and uncontrolled free (UCF) flows. Under these conditions, the riverbed on the downstream of the barrage experiences different flow patterns, i.e., scour and bed retrogression. At present, the focus is only made to the UCF flows of 24.30 and 17.5 m³/s/m for which the inflow Reynolds number reaches 6,479,978 and 3,980,015, respectively. For creating the UCF flows, the present models are operated at H/H_o = 0.998 and 0.988 [54] under different pond and tailwater levels as provided in

Table 3. Operational condition for the design and free flow models.

Discharge (m3/s/m)	Pond Level (m)	Tailwater Elevation (m)	Barrage Operation
24.30	135.93	134	Free design Flow
17.5	134.60	133.65	Uncontrolled Flood flow

, whereas H and H_o are the design and operation heads, respectively. Formula (6) is employed to calculate the studied discharges.

$$\mathrm{Q}=\frac{2}{3}\ast {\mathrm{C}}_d\ast B \ast H_0^{3/2}\ast \sqrt{2g}$$

(6)

where Q is the discharge, which is measured in m³/s/m; C_d and B (m) are the discharge coefficient and the width of bay, respectively; H₀ is the net energy over the crest, which is measured in m; g (m/s²) is the acceleration due to gravity; and P is the weir height. For the design discharge, the used and obtained values for C_d are reached as 0.640 and 0.626, respectively. In contrast, for the flood discharge, a C_d value of 0.726 is used in the simulations, whereas the model value reaches 0.716. The schematic diagram of UCF flows is presented in Figure 6.

2.6. Simulation Time, Models’ Stability, and Convergence

For the present models, the stability solution and models’ convergence are achieved through Courant number criterion [32,33]. For the present models, the computations are performed on a personal computer with 2.8 GHz core (i7, 11 generation) and 8 GB memory. The average time steps for design and flood discharges varied between 0.003 to 0.0032 and 0.004 to 0.0032, respectively. For higher flows of 24.30 and 17.5 m³/s/m, the models’ stability can only be monitored by fluid kinetic energy (MAFKE_s) at the upstream and downstream boundaries. After the analysis, T_s = 500 s is found to be acceptable because all models have shown stable values of MAFKE_s with little variation from the mean values. However, the average actual time (T_a) for the investigated models is reached at 96 hr. Presently, the steady state of design discharge is discussed here.

Two different turbulence schemes are employed to investigate the scour and bed retrogression. At 24.30 m³/s/m, using LES turbulence scheme, at T_s = 485 s, the simulation indicated 0.20% mean variation in hydraulic parameter values in Type I basin and, at T_s = 500 s (finish time), the maximum MAFKE value reached 8.32 m²/s² in the LES model, while, even at T_s = 500 s, no stability condition appeared in K-ϵ model. However, at T_s = 500 s, the averaged value of MAFKE reached 8.36 m²/s², which showed an increment of 0.50% compared to the LES model. At 24.30 m³/s/m, using LES model, in Type II basin, the model showed minimum variation in hydraulic parameter values at different intervals of time and, at these time instants, the variations from the mean values, i.e., MAFKEs, reached 0.30%. However, using the K-ϵ model in Type II basin, the simulation did not show any stability because the mean variations in parameter values increased, which reached 2%. At T_s = 500 s, the maximum MAFKEs in Type II stilling basin were reached at 8 and 8.56 m²/s² in LES and K-ϵ models, respectively. In comparison to the LES model, 8% higher MAFKs values were observed in K-ϵ models, which indicated large fluctuations in free surface due to the continuous bed retrogression. Using LES turbulence scheme, in Type III basin, at T_s = 336 s, the mean variation in MAFKEs reached 0.60%; however, at the finish time, these values reached 8.46 m²/s². In contrast, no stability condition was noted in K-ϵ model and, at the finish time, the maximum MAFKE value in Type III basin reached 8.64 m²/s², which was 2% higher than that observed in the LES model. From the stability and convergence analysis, it is found that, due to the continuous bed retrogression, the models did not achieve hydraulic and scour stabilities; however, at the assigned finish time, the variations in mean MAFKEs values at different time instants were found to be much less and the differences in values between the two different turbulence models were up to 2%. Therefore, for the present study, the finish time at which the sediment bed is completely exposed is set as a simulation time (T_s = 500 s).

3. Results

3.1. Model Verification and Validation

At present, the stilling basin of Taunsa barrage is a modified form of USBR Type II (Type II basin) which includes chute blocks and end sill. The barrage has 64 gates to discharge the upcoming flows, which are named bays. Downstream of these bays, the field authorities periodically measure the bed profiles. For the validation of present scour models, the bed profile data, i.e., maximum scour depths, of bay 33 of year 2015 is utilized. The reason for choosing year 2015 data is because, just before this year, during June to September 2014, a high flood discharge passed from the barrage. The bed profile data of bay 33 were established during off-peak months, i.e., January to February 2015. Therefore, presently, a full-scaled single bay of the studied barrage is modeled, and the model’s assessment is made with the maximum scour depths at different locations on the downstream of bay 33 of the prototype (Type II basin).

The data in the hard form are digitized in AutoCAD software. It is important to mention that, in the field data, the scour depths are measured at 10 different points in the flow direction [55]. The field staff measure these data points from the centerline of the bays. From the available 10 locations, the maximum scour depth values are taken for comparison. Hence, to compare the bed profiles downstream of Type II basin, the data from the computed scour depths are also obtained from the similar locations. Downstream of the rigid bed, up to X = 6, both LES and RNG K-ϵ turbulence models predicted identical results, which showed good agreement with the field data. However, afterward, due to the increase in recirculation and eddies in the scour hole, LES model showed higher length and depth of scour hole, which deviated from the field data. In the LES model, the overall change in the sediment bed reached 97%, while the results of RNG K-ϵ model agreed well with the field data, as shown in Figure 7.

The performance of present scour models is also assessed by the coefficient of determination (R²) and Nash–Sutcliffe model efficiency coefficient (NSE). The model performance is said to be perfect when the value reaches 1, while it is said to be poor if the values reach 0 and negative for R² and NSE, respectively.

Table 4. Model performance indicators for the compared field data.

Turbulence Modes	Statistical Indices
	R2	NSE
RNG K-ϵ	0.9094	0.896
LES	0.900	0.423

shows the comparison of LES and RNG K-ϵ models. From Table 4, it can be noted that the results of R² and NSE are found to be more promising in RNG K-ϵ model, which reached 0.9094 and 0.896, respectively. In contrast, the LES model showed acceptable values of R², while the NSE (0.423) was found to be less than RNG K-ϵ model. The deviation in scour profile by LES can be due to (1) the LES model being mesh-dependent, which requires finer cell size in the turbulent flow zones to capture the small scale fluctuations, i.e., turbulent eddies, and (2) within the scour hole, strong vertexes are generated, for which the model employed a sub-grid turbulence model, which captured small eddies in the scour hole and consequently increased the scour depth and overall retrogression on the downside of the sediment bed, as stated by Man et al. [56] and Hamid et al. [57]. In contrast, due to the limitation in computing the wake velocities in the shear region, the RNG K-ϵ model predicted less scour in the scour hole and on the downward areas. However, in the real fields, at such higher discharges, i.e., 24.30 and 17.5 m³/s/m, only forward velocity profiles are noticed on the downstream of hydraulic structures, which only retrogresses the riverbed.

Conclusively, based on the results of scour profiles, it is believed that the present FLOW-3D models showed reasonable accuracy of scour profiles. From the comparison with the field data (Figure 7) and the results of statistical indices, as compared to the LES model, the RNG K-ϵ model showed agreement with the field data; thereby, in the proceeding sections, for the studied discharges, the velocity fields, bed shear stress, and scour profiles downstream of Type I, II, and III basins are only focused on RNG K-ϵ model.

3.2. Flow Field

The present numerical models have focused on the effects of different basins’ appurtenances on the downstream riverbed; thereby, the velocity profiles are described both for the rigid and scoured beds. However, the other parameters, such as shear stress and bed profiles, are discussed only for the scoured and retrogressed beds. The results of the employed turbulence model are analyzed at T_S = 500 s, for which 2D and 3D illustrations of velocity, shear stress, and bed profiles are explained for the studied discharges.

3.2.1. Velocity Distribution at Design Discharges

Figure 8 shows velocity distribution downstream of Type I, II, and III basins at 24.30 m³/s/m discharge using K-ϵ model. Figure 8a shows velocity patterns downstream of Type I basin and the results revealed highly turbulent flow on the fore side of HJ. The maximum forward and backward velocities were found at the downstream glacis and in the HJ regions, which reached 9 m/s and −2.56 m/s, respectively. The results also showed large fluid recirculations at the free surface. After the rigid floor, the velocity field on the scoured bed showed two different flow zones, i.e., slow-moving zone near the scoured bed and high-velocity zone from the central fluid region to free surface. In addition, the results also showed forward velocity contours downstream of Type I basin, which eroded/retrogressed the bed. At the initial length of the bed, the depth of the slow-moving wake zone was lower, which increased as the flow moved downstream.

Using K-ϵ turbulence model, the results indicated that a concentrated supercritical jet was found to be impinging on the floor of Type II basin, as shown in Figure 8b. The maximum velocity was noticed from the basin’s floor to the central fluid depth, which ranged between 6 and 8 m/s. However, at the free surface of HJ, the results showed rollers and recirculations. Soon after the end sill, the results showed wake zones and slow-moving flow, which caused bed erosion and developed a large scour hole. Additionally, the results also indicated that, as the depth of the low-velocity zone increased, the dimensions of scour hole and the overall retrogression was also found to be increased. Using K-ϵ model, Figure 8c shows velocity distribution downstream of Type III basin. Similar to Type I basin, the maximum velocity was observed on the fore side of HJ; however, its pattern was found to be dissimilar to that noticed on the downstream of Type I and II basins. In Type III basin, the results indicated a stable free surface profile on rigid bed, which indicated that, shortly, the flow would achieve hydraulic stability. In Type III basin, due to the development of large wake zones around the appurtenances, the velocity values on the basin’s floor were found to be decreased. Like the rigid bed, as compared to Type I and II basins, a different flow pattern was observed on the scoured bed downstream of Type III basin. Furthermore, the result also indicated that, as the flow travelled to the end of the sediment bed, due to the increasing depth of the wake zone, the bed retrogression was also found to be increased at those locations.

3.2.2. Velocity Field at High Flood Discharge

In Section 3.2.1., the velocity fields are analyzed downstream of the studied basins at the design discharge of 24.30 m³/s/m. The results have indicated unstable free surface downstream of Type I and II basins, while a relative stable profile is noticed downstream of Type III basin. Additionally, such high discharges usually do not occur in the rivers; therefore, this study has also investigated velocity fields on 75% of the design discharge (17.5 m³/s/m).

Using K-ϵ model, at 17.5 m³/s/m discharge, large fluid recirculations were observed at the free surface of HJ, while, after the HJ, the free surface was found to be stable on Type I basin, as shown in Figure 9a. In comparison to the design discharge, a large wake region was noticed after the first row of baffle blocks. Furthermore, downstream of Type I basin, the depth of the high-velocity zone was found to be reduced, which decayed before the end of the sediment bed. On the retrogressed bed, the depth of slow-moving flow was less than that observed at the design discharge; therefore, the results indicated less bed retrogression at 17.5 m³/s/m discharge. Figure 9b illustrates the velocity field on Type II basin using K-ϵ model at 17.5 m³/s/m discharge. After the downstream glacis, the supercritical flow was striking on the basin’s floor, which travelled to the end of the rigid bed. On the rigid bed of Type II basin, three different velocity zones were observed: (1) higher velocity zones from floor to central fluid depth, (2) intermittent velocity zones, and a lower velocity zone near the free surface. Out of the three zones, the results indicated the maximum velocity values in the first zone, which reached 6.5 m/s. On the scoured bed, a large fluid recirculation and wake zones were observed, which developed a scour hole; however, the net length and depth of the wake zone was found to be less than that observed at the design discharge. After the scour hole, as compared to the design discharge, the depth of the slow-moving zone on the scoured bed was also found to be reduced, which consequently decreased the overall bed retrogression.

On the other hand, at 17.5 m³/s/m discharge, the velocity pattern in Type III basin was found to be dissimilar to that observed in Type I and II basins, as illustrated in Figure 9c. The results showed higher fluid recirculation in the roller region of HJ and indicated a stable HJ. The results further indicated that, after impacting with WSBB, the supercritical flow was directed towards the free surface. Furthermore, after the WSBB, a large wake zone region was noticed, which further reduced the magnitude of velocities near the basin’s floor. Due to the reduced velocity at the rigid floor, the depth of the slow-moving zone on the erodible bed was also reduced on the downside of Type III basin. At the start of retrogressed bed, the depth of slow-moving fluid was lower, whereas, as the fluid moved towards the end of the sediment bed, the depth of slow-moving fluid also increased, which consequently increased the erosion process in those regions. However, as compared to the Type I and II basin, the depth of central and higher velocity zones were found to be lower in Type III basin, which indicated that the high-velocity currents were moving above the central fluid depth and at the free surface.

3.3. Bed Shear Stress (BSS)

Another parameter that influences local scour is the local bed shear stress (BSS) on the bed, which is calculated by wall function for the three-dimensional turbulent flow. Therefore, the contours are drawn to illustrate the BSS on different bed locations of the studied discharges.

3.3.1. Bed Shear Stress (BSS) at Design Discharge

Figure 10 shows the BSS contours downstream of Type I, II, and III basins at 24.30 m³/s/m using K-ϵ turbulence model. From Figure 10a, the maximum magnitude of BSS downstream of Type I basin was found on the left side of the sediment bed, which increased the scour depths in those regions. Figure 10b shows the BSS downside of Type II basin using K-ϵ model. From Figure 10b, the results indicated maximum BSS within the scoured hole and on the downside of the scour hole. Figure 10c shows shear stress contours downstream of Type III basin and the results indicated that the maximum magnitude of BSS occurred on the right side of the sediment bed.

Figure 11 shows BSS downstream of Type I, II, and III basins at the design discharge. The plots are drawn from the centerline (Y/2) of the bed. Downstream of Type I basin, the BSS gradually increased as the flow moved toward the end of the sediment bed, whereas, downstream of Type II basin, the maximum BSS was found in the scour hole region. After the scour hole, the shear stress rapidly decreased due to the large fluid recirculation, which excessively retrogressed the bed. It is important to mention that, as the bed was eroded, the shear stress was found to be decreased. On the other hand, downstream of Type III basin, the BSS gradually increased towards the end of the sediment bed. Out of the studied basins, the maximum BSS was noticed downstream of Type III, which indicated the sediment particles were still intact and required higher critical shear stress to remove the particles from the bed surface. In conclusion, in comparison to Type II and Type III basins, the minimum shear stress was observed downstream of Type I basin, which reached 3 P_a.

3.3.2. Bed Shear Stress at High Flood Discharge

Figure 12 shows BSS on the sediment bed of Type I, II, and III basins at 17.5 m³/s/m discharge. As compared to the design discharge, at 17.5 m³/s/m discharge, different patterns of the BSS were noticed on the downstream of Type I basin, which indicated a uniform distribution of shear stress from X = 10 to 30 m. However, at the end of sediment (from X = 25 to 28 m), the BSS on the right side of the bed was found to be higher than the rest of the locations, as presented in Figure 12a. At 17.5 m³/s/m discharge, Figure 12b shows BSS on the scoured bed of Type II basin using K-ϵ model. It can be seen from Figure 12b that the maximum BSS occurred in the regions of scour hole and at the end of the retrogressed bed. In contrast, at the design discharge, the maximum BSS was found at the downstream end of the scour hole. At 17.5 m³/s/m discharge, after the scour hole, the distribution of BSS was found to be uniform. Conclusively, as compared to the design discharge, at 17.5 m³/s/m discharge, the bed shear stress was found to be dissimilar on the downstream of Type II basin.

Using K-ϵ model, Figure 12c shows the distribution of BSS on retrogressed bed of Type III basin at 17.5 m³/s/m discharge. The results showed a uniform distribution of BSS on the retrogressed bed downstream of Type III basin. The maximum shear stress was found from X = 18 m to 30 m. As compared to the design discharge, at 17.5 m³/s/m discharge, the patterns of BSS were found to be regular. However, as compared to Type I basin, at 17.5 m³/s/m discharge, fewer sediment regions were found to be affected by the BSS downstream of Type III basin.

Figure 13 shows BSS profiles downstream of Type I, II, and III basins at flood discharge. As compared to the design discharge, at flood discharge, the trends of BSS profiles are found to be different downstream of Type II and Type III basins, while the results indicated similar trends on the downside of Type I basin. In Type I basin, the amount of BSS was found to be increased as the fluid moved towards the end of the sediment bed. On the retrogressed bed of Type I basin, the maximum BSS reached 3.49 P_a. Downstream of Type II basin, up to 10 m, a nonlinear trend was noticed, while, in the scour hole, the shear stress abruptly decreased, which reached 0.43 P_a. After the downside end of the scour hole, the shear stress started rising till the end of the sediment bed, as shown in Figure 13. In contrast, as compared to the design discharge, a dissimilar BSS profile was noticed on the downstream of Type III basin. The results further indicated that, at the flood discharge, the flow produced less shear stress; thereby, less bed was eroded on the downstream of Type III basin, as can be seen from Figure 13.

3.4. Local Scour and Bed Retrogression

Section 3.2 and Section 3.3 described the velocity and BSS patterns on the downstream of Type I, II, and III basins using RNG K-ϵ turbulence models. Following the flow fields, this result section focuses on the local scour and bed retrogression downstream of the studied basins.

3.4.1. Scour and Retrogression Patterns at Design Discharge

Figure 14 shows scoured bed downstream of different basins using K-ϵ model. The 3D plots are drawn at the finish time of T_s = 500 s. Figure 14a indicates scoured bed downstream of Type I, and the results showed maximum bed retrogression at the end of sediment, for which the maximum scour depth reached 2.44 m. Near the rigid floor, the sediment bed was less retrogressed; however, as the flow moved towards the downside, the retrogression was found to be increased.

Using K-ϵ model, soon after the rigid floor, the results revealed a scour hole on the sediment bed downstream of Type II basin, as shown in Figure 14b, and, at T_s = 500 s, the overall length and depth of the scour hole reached 14 m and 3.57 m, respectively. In addition, after the scour hole, the riverbed continued to retrogress and, at T_s = 500 s, the riverbed was found to be completely exposed, as shown in Figure 14b. Figure 14c shows the scour and retrogression patterns downstream of Type III basin using K-ϵ model. The results showed only bed retrogression; however, the retrogression patterns were found to be dissimilar to those observed in Type I and II basins.

3.4.2. Scour and Retrogression Patterns at High Flood Discharge

As compared to the design discharge, at 17.5 m³/s/m, the overall scour patterns are found to be different downstream of studied basins. Figure 15 describes the scour and retrogression patterns at the high flood discharge of 17.5 m³/s/m.

Figure 15a shows retrogression patterns downstream of Type I basin at 17.5 m³/s/m discharge. Using K-ϵ model, the maximum scour depth reached 2.87 m at the end of the sediment bed, as illustrated in Figure 15a. As compared to the design discharge, at 17.5 m³/s/m discharge, 14% higher scour depth was observed. However, the pattern of bed retrogression was found to be identical as witnessed at the design discharge. Figure 15b displays the scour pattern downstream of Type II basin at 17.5 m³/s/m discharge. Similar to the design discharge, a scour hole was noticed near the rigid floor for which the maximum scour depth reached 2.93 m, which was about 6% less than was noticed at the design discharge. Figure 15c shows bed retrogression patterns downstream of Type III basin. As compared to Type I and II basins, the results indicated dissimilar scour patterns downstream of Type III basin, as presented in Figure 15c. In addition, on the downstream of Type III basin, the results also showed ripple and wavy bed surface on the entire retrogressed bed. At the flood discharge, the overall change downstream of Type III basin reduced to 10% and 37% compared to Type I and II basins, respectively. In conclusion, as compared to the design discharge, the overall change in the sediment bed downstream of Type III basin was about 19% lower.

4. Discussion

The former Section 3 has presented the results of the models’ validation, flow fields, i.e., velocity profiles, and bed shear stress for the investigated basins at the design and high flood flows. However, keeping in view the scour and bed retrogression, this section mainly compares the performance of different stilling basins at the two studied flows. The discussion is further extended by comparing the results of net scour and retrogression on the downstream of different basins.

At the design flow, using K-ϵ turbulence model, due to the presence of baffle and friction blocks, less flow velocity was observed near the rigid floor of Type I basin. Near the sediment bed, a slow-moving fluid zone was noticed, which triggered the bed erosion. However, except bed retrogression, the results did not indicate any scour hole on the sediment bed. At the design discharge, the overall change in the sediment bed of Type I basin reached 51%. In contrast, due to the absence of basin appurtenances, a high-velocity flow was found to be impinging on the floor of Type II basin, which travelled to the entire length of the rigid floor. Soon after the rigid floor, a large wake zone was found on the sediment bed, which developed a large scour hole near the rigid floor. The analysis further revealed that, as the depth of the slow-moving zone increased, the sediment bed was also found to be eroded at higher rate. Downstream of Type II basin, the overall bed change reached 75% and 100% in the scour hole and at the end of the sediment bed, respectively. In contrast, compared to the Type I and Type II basins, dissimilar flow fields were noticed on the downstream Type III basin. In addition, due to the inclusion of WSBB, considerable wake regions were developed near the bed of Type III basin and, compared to Type I basin, the free surface profile was also found to be stable. Furthermore, it is important to mention that, downstream of Type III basin, only bed retrogression was noticed as observed in Type I basin. Hence, from the present results, it can be said that the WSBB basins can be used downstream of diversion barrage up to the investigated discharges. However, at the design discharge, the maximum bed retrogression on the downstream of Type III reached 67%, which was about 16% higher than the Type I basin.

At the flood flow, using K-ϵ model, the results showed a stable free water profile on the rigid bed of Type I basin and, compared to the design discharge, the results indicated lower depth of the wake zone near the scoured bed. Due to the hydraulic stability at free surface and lesser wake zones near the scoured bed, the overall retrogression was found to be lower at the flood discharge. However, the results showed higher bed retrogression at $\frac{y}{4}$ and $\frac{3y}{4}$ on the downstream of Type I basin, as given in

Table 5. The overall change in beds downstream of the Type I, II, and III basins.

Stilling Basins	Flow (m3/s/m)	Sediment Bed Retrogressed (%)
		RNG K-ϵ Model
		$\frac{\mathit{y}}{4}$	$\frac{\mathit{y}}{2}$	$\frac{3\mathit{y}}{4}$
Type I	23.30	69%	51%	64%
Type II		100%	100%	97%
Type III		60%	67%	69%
Type I	17.50	60%	58%	54%
Type II		87%	85%	82%
Type III		46%	48%	48%

. At the flood flow, the overall change in the bed of Type I basin reached 67%. Similar to Type I basin, the results showed lower depth of the slow-moving zone near the scoured bed and, due to this reason, the overall change in the sediment bed was also found to be less. In addition, the wake zone in the scour hole and fluid recirculation at the scoured bed was also found to be less; as a result, in comparison to the design flow, a lower depth of scour hole was noticed on the downstream of Type II basin, which reached 84%.

In contrast, at the high flood discharge, the water surface profile on Type III basin was found to be stable compared to the Type I and II basins. Additionally, compared to the Type I basin, in WSBB block basin, the depth of the high-velocity zone was found to be lower on the rigid bed and WSBB blocks also generated large wake zones on the rigid bed. In conclusion, on the downstream of Type III basin, the overall change in the sediment bed was found to be less than that noticed at the design discharge. In comparison to the Type I basin, the total change in the bed was about 20% less on the downstream of Type III basin. From the results, it can be said that, at the design discharge, out of the employed basins, the Type II basin produces higher bed retrogression. Similar to the design discharge, at the flood discharge, the maximum bed change was also noted on the downstream of Type II basin.

In conclusion, out of the studied basins, the scour results downstream of Type I basin are found to be lower at the design discharge. However, due to the complex scour phenomena, at $\frac{y}{4}$, the Type III basin showed better results than the Type I basin, as shown in Table 5. Based on the results and the bed profile analysis, the study revealed that the WSBB basin downstream of the studied barrage has reduced the scour and bed retrogression up to 17.5 m³/s/m discharge, while the efficiency of WSBB basin is found to be slightly reduced at the design discharge. Additionally, from the results of velocity field, bed shear stress, and scour patterns, it is believed that the remodeled basin (Type II) of Taunsa barrage is dissipating less energy and producing high-velocity currents on the riverbed.

5. Conclusions

The present study developed numerical models to study local scour downstream of Taunsa barrage for design and flood discharges. Field data of riverbed profile are employed to validate the present scour models. Based on the validation, this study examined velocity pattern, bed shear stress, and scour patterns downstream of three different stilling basins. The following conclusions are drawn from the current study:

• Out of the tested turbulence models, RNG K-ϵ showed agreement with the field data, for which the values of R² and NSE reached 0.9094 and 0.896, respectively.
• From the results, the study revealed a slow-moving velocity region near the beds, which produced scour and bed retrogression. Additionally, the analysis also revealed a forward velocity profile downstream of Type I and III basins. On the other hand, large recirculating regions were noted downstream of Type II basin, which developed a scour hole in those regions.
• At the design and flood discharges, Type I and Type III basins showed only bed retrogression, while large scour hole and excessive bed retrogression were observed downstream of Type II basin. In addition, up to 17.5 m³/s/m discharge, the stilling basin with WSBB (Type III basin) controlled the bed retrogression, while, at the design discharge, as compared to Type I basin, higher bed retrogression was noticed at the centerline $\left(\frac{\mathrm{y}}{2}\right)$ of the model. In contrast, at the design discharge, the riverbed was found to be completely exposed downstream of Type II basin, which revealed that Type II basin developed high-velocity flow on the basin’s floor, which consequently eroded the riverbed. Similar to the design flow, at flood discharge, 85% of the riverbed was found to be eroded downstream of Type II basin.

Based on the models’ results, it is concluded that, out of the tested basins, Type II basin generates turbulent flows on the riverbed, which has developed a scour hole near the basin’s floor. Conversely, the newly proposed WSBB basin has improved the energy dissipation downstream of the studied barrage up to 75% of the design discharge. Additionally, from the results of investigated parameters, it is believed that WSBB improves the flow near the basin’s floor. However, as the present study is limited to Van Rijn transport rate equation, this study recommends investigating the riverbed profile downstream of the studied barrage with Meyer-Peter Müller and Nielsen equations. In addition, the study also suggests examining the effects of multiple gates openings on the local scour downside of the investigated basins.