Numerical Investigation of Critical Hydraulic Parameters Using FLOW-3D: A Case Study of Taunsa Barrage, Pakistan

: Hydraulic structures, such as barrages, play an important role in the sustainable development of several regions worldwide. The aim of this novel study is to identify the critical hydraulic parameters (CHPs) of Taunsa Barrage, built on the Indus River. These CHPs, including free surface proﬁles, ﬂow depths, Froude number, velocity proﬁles, energy dissipation and turbulence kinetic energy, were investigated using simulation via FLOW-3D numerical models. Incompressible Reynolds-averaged Navier–Stokes (RANS) equations on each computational cell were solved using the numerical methods available in FLOW-3D. The simulation results indicated that the locations of hydraulic jumps (HJs) were lower than that were reported in the previous one-dimensional study. Similarly, the distances of the HJs from the downstream toe of the glacis were reached at 2.97 m and 6 m at 129.10 m and 130.30 m tailwater levels, respectively, which deviated from the previous studies. In higher tailwater, the sequent depth ratio also deviated from the previous data. The maximum turbulent kinetic energies were observed in the developing regions of HJs, which were found to be decreased as the distance from the HJ was increased. The results of this research will be highly useful for engineers working in the ﬁeld of design of hydraulic structures.


Introduction 1.Significance of Hydraulic Jumps (HJs)
Due to the sufficient head upstream of hydraulic structures, the outflow water has immense kinetic energy that can damage the downstream structures.Researchers and hydraulic engineers have devised many measures such as baffle blocks, friction blocks, chute blocks, and end sills to control the above-mentioned issues.These devices stabilize the hydraulic jumps (a vital energy dissipation process) formed at the location.Hydraulic jumps occur in natural systems and artificial channels such as rivers, streams, spillways, sluice gates, barrages, and weirs.At the different locations of these hydraulic structures, the flow passes through different conditions, i.e., subcritical (Froude number ((Fr) < 1), critical (Fr = 1), and supercritical (Fr > 1)).In the subcritical condition, the actual flow depth is higher than the critical, whereas in the supercritical condition the flow depth is always found to be less than critical.In the hydraulic jumps, the flow changes suddenly from supercritical to subcritical conditions and vice versa.During this chaotic phenomenon, a rapid rise in the free surface occurs that dissipates a large amount of energy due to the turbulent mixing [1][2][3][4].A hydraulic jump (HJ hereafter) occurs in gravity-driven flows when Fr crosses unity and is defined as the ratio of inertial to gravitational forces.
Numerous researchers have conducted experiments on HJs by employing various geometries and hydraulic conditions, such as vegetated bed (Bai et al. [5]), close conduits Fluids 2023, 8, 310 2 of 26 (Li et al. [6]), stilling basins of open-channel hydraulic structures [7-13]; and horizontal smooth and rough beds [14].In addition, Balachandar et al. [15] investigated the effects of tailwater variation on HJs and studied downstream bed profiles, whereas [16][17][18][19][20] investigated submerged HJs for different beds and flow conditions.Furthermore, a few studies have investigated the transition region of an HJ (Zobeyer et al. [21]; Abbaspour et al. [22]) and the weak HJ (Mignot and Cienfuegos [23]) on rough beds and in small open channels.On the whole, these studies tested the basin's appurtenances to increase the energy dissipation and mainly investigated the gate openings, Fr, Reynolds number, roller lengths of HJs, sequent depths, tailwater, momentum and velocity decay, and turbulent kinetic energy.However, the literature did not reveal any detailed experimental study on the flow patterns downstream of the Taunsa barrage, which is an important research area because of the significant importance of the barrage.

Role of Stilling Basins
Stilling basins are accepted for the dissipation of the surplus kinetic energy of an HJ downstream of spillways, dams, barrages, and pipe outlets.The hydraulic performance of the stilling basin is very much dependent on its shape and size, which affect the flow patterns [24,25].Furthermore, a substantial amount of energy is still to be dissipated downstream of the stilling basin, for which flexible aprons are provided.The risk of scouring is also safeguarded by the flexible aprons, which alleviate the uplift pressure that is left behind.Ali and Mohamed [26] and Mishra [27] conducted experiments to study the effects of stilling basins shapes downstream of radial gates.Alikhani et al. [28] conducted experiments on a single vertical sill in a stilling basin for a forced HJ.Elsaeed et al. [29] investigated the effects of end steps on the length of submerged HJs and measured the energy losses and velocity profile along the stilling basin.Tiwari et al. (a, b) [30,31] investigated the stilling basin shape for pipe outlets and used different shapes of intermediate sills with heights equivalent to the diameter of the pipes.Hager and Li [32] and Herrera-Granados and Kostecki [33] conducted experiments to investigate the effects of different energy dissipators on the characteristics of HJ.Ali and Kaleem [34] investigated the energy dissipation of Taunsa Barrage stilling basins before and after remodeling.The results showed that, due to the remodeling of the basin, the stone apron was launched during the 2010 flood and the river course also shifted towards the left side.Chaudary and Sarwar [35] and Chaudhry [36] analyzed the tailwater effects on different stilling basins of Taunsa Barrage.The results indicated that the existing tailwater levels downstream of the prototype barrage were appropriate for the formation of an HJ.

Computational Fluid Dynamic and Hydraulic Modelling
The previous sections (Sections 1.1 and 1.2) discussed the methods and the significance of physical and experimental investigations of HJs and basin appurtenances, which can be assisted by three dimensional (3D) numerical codes.Additionally, experimental investigations and on-site measurements are usually expensive and time-consuming and due to hindrance of flow devices and scaling effects, the output results are found to deviate from the prototypes.Therefore, the use of numerical modelling to investigate the hydraulic characteristics of grade-control structures is becoming popular.Such modeling tools are helpful, especially when the basic fundamental equations are not adequate to give solutions such as multifaceted geometries [33].Recently, with the development in computer technology, the problems of complicated hydraulic structures can be studied by improved via numerical methods and turbulence models, which further harmonize physical modeling in the design stages.These numerical models also test various feasible flow phenomena, with little change in the input to obtain further data for the computed domains [14,37,38].Many numerical investigations are found on HJs and energy dissipation using various numerical codes; however, Table 1 1, some studies, e.g., [39,40,43,44,47,50,51,56], have focused HJs and their associated parameters such as velocity, free surface profile, sequent depths, roller lengths, and turbulent kinetics energy, whereas other researchers have investigated energy dissipation [41,45,46,48,[53][54][55] and scour [42,49] downstream of different hydraulic structures.However, except for the 1D Hydrologic Engineering Center River Analysis System (HEC-RAS) study on Taunsa barrage, the literature did not reveal any study that studied the effects of tailwater on the HJ and flow pattern in the barrage's basin.

Research Motives and Problem Statement
Soon after the operation of the Taunsa barrage, due to the retrogression of the downstream riverbed, the tailwater levels were lowered, which consequently damaged the basin's floor; during 1959-1962 these regular issues were resolved.Additionally, due to the structural flaws, some of the baffle blocks were also found to be uprooted [55].After so many years of partial repairs, it was declared that the barrage was to be remodeled; it was reported that due to the lowering of the tailwater levels, the HJ was sweeping on the floor.On the basis of a model study report, the basin was remodeled; however, even after remodeling, the data from the 2010 flood revealed damage downstream of barrage [57,58].To understand the tailwater levels and HJ locations, Chaudhry [36] carried out a one-dimensional HEC-RAS study.
In the year 2008, Taunsa barrage was remodeled based on a sectional model study on the rigid bed.Except tailwater, no other hydraulic parameters were thoroughly investigated.Furthermore, these kinds of physical studies are generally associated with scaling effects such as roughness and length.However, due to the advent of computer technology and advanced turbulence models, as compared with HEC-RAS, more dedicated hydraulic modeling tools are available to investigate hydraulic issues such as those found downstream of the Taunsa barrage.Therefore, in this study, firstly, using the frequency of occurrence and relative importance index (RII), the critical hydraulic parameters (CHPs) are identified from the literature.The prime objective of this study is to investigate the tailwater and HJ locations at 44 m 3 /s discharge downstream of the old basin of Taunsa barrage using FLOW-3D.This study further investigates the identified CHPs on the basin before its remodeling and mainly focuses on free surface profiles, flow depths, Froude number, velocity profiles, and turbulent kinetic energy.Additionally, the results of CHPs are also compared with the available data from prototypes, HEC-RAS, and studies from the literature.The specific research goals for the present numerical study are as below:

•
To calibrate and validate a 3D model under the field conditions of a full-scale hydraulic structure using data from Taunsa Barrage.

•
To identify the critical hydraulic parameters (CHPs) that play a crucial role in the design of graded control structures.

•
To investigate the flow patterns of CHPs for stilling basins at various tailwater levels using the data from the old Taunsa Barrage.

•
To study the effects of different tailwater levels on the locations of the HJs and compare them with relevant field and numerical data.

•
To provide an overview of HJ locations in various basins of a barrage for different tailwater and discharge levels (a case study for Taunsa Barrage).

Study Area
Pakistan is an agricultural country; the major source of economy and livelihood depends on the agriculture sector.This sector provides about 25% of GDP and engages 50% of labor from rural areas [56].More than 18 million hectares (ha) of land is irrigated by the Indus River and its branches.In the Indus Basin, barrages are essential components that divert water into the canals and also serve as roads, bridges and power transmission lines [56].Taunsa Barrage was completed in 1958 on the Indus River.The barrage diverts water to Dera Ghazi Khan Division through the Muzaffargarh and Taunsa-Panjnad link canals and irrigates about 809,371 ha of land.The barrage was designed for a discharge capacity of 28,313 m 3 /s.The total length of the barrage is 1325 m, whereas the clear width for flow passage is 1171 m.The maximum upstream and downstream flood levels are 136.94m and 135.33 m, respectively, whereas 136.24 m is the normal pond level for the operation of the canals.The upstream and downstream floor levels are designed at 128.31 m and 126.79 m, respectively.The weir's crest is located at 130.44 m, whereas the waterfall is kept at 3.66 m.It should also be noted that all the elevations given herein are from the mean sea level.Taunsa barrage's location and its typical cross section is shown in Figure 1.
Fluids 2023, 8, x FOR PEER REVIEW 4 of 30 identified from the literature.The prime objective of this study is to investigate the tailwater and HJ locations at 44 m 3 /s discharge downstream of the old basin of Taunsa barrage using FLOW-3D.This study further investigates the identified CHPs on the basin before its remodeling and mainly focuses on free surface profiles, flow depths, Froude number, velocity profiles, and turbulent kinetic energy.Additionally, the results of CHPs are also compared with the available data from prototypes, HEC-RAS, and studies from the literature.The specific research goals for the present numerical study are as below:

•
To calibrate and validate a 3D model under the field conditions of a full-scale hydraulic structure using data from Taunsa Barrage.

•
To identify the critical hydraulic parameters (CHPs) that play a crucial role in the design of graded control structures.

•
To investigate the flow patterns of CHPs for stilling basins at various tailwater levels using the data from the old Taunsa Barrage.

•
To study the effects of different tailwater levels on the locations of the HJs and compare them with relevant field and numerical data.

•
To provide an overview of HJ locations in various basins of a barrage for different tailwater and discharge levels (a case study for Taunsa Barrage).

Study Area
Pakistan is an agricultural country; the major source of economy and livelihood depends on the agriculture sector.This sector provides about 25% of GDP and engages 50% of labor from rural areas [56].More than 18 million hectares (ha) of land is irrigated by the Indus River and its branches.In the Indus Basin, barrages are essential components that divert water into the canals and also serve as roads, bridges and power transmission lines [56].Taunsa Barrage was completed in 1958 on the Indus River.The barrage diverts water to Dera Ghazi Khan Division through the Muzaffargarh and Taunsa-Panjnad link canals and irrigates about 809,371 ha of land.The barrage was designed for a discharge capacity of 28,313 m 3 /s.The total length of the barrage is 1325 m, whereas the clear width for flow passage is 1171 m.The maximum upstream and downstream flood levels are 136.94m and 135.33 m, respectively, whereas 136.24 m is the normal pond level for the operation of the canals.The upstream and downstream floor levels are designed at 128.31 m and 126.79 m, respectively.The weir's crest is located at 130.44 m, whereas the waterfall is kept at 3.66 m.It should also be noted that all the elevations given herein are from the mean sea level.Taunsa barrage's location and its typical cross section is shown in Figure 1.The stilling basin of the barrage is of a modified form similar to the United States Bureau of Reclamation's (USBR's) stilling basin type-III.Within the basin, two rows of baffle and friction blocks facilitate energy dissipation and stabilize the HJ even in cases of minimum tailwater requirements.However, soon after the barrage construction, multiple problems such as the oblique right-sided river approach, heavy siltation in one of the major canals, uprooting of impact baffle blocks, damage to the stilling basin floor, and retrogression of water levels appeared in its downstream areas.During 1959During -1962During , 1966During , and 1973, repairs works were carried out to cater to the problems mentioned above; however, the problems remained persistent.In the reports by Zaidi et al. [58] and the World Bank [57], sweeping of the HJ was believed to be the main reason for the issues highlighted above.To resolve these issues, the Punjab Government constituted committees of experts; however, no specific measures were taken and the issues continued to be aggravated.The typical cross section of the Taunsa barrage stilling basin is shown in Figure 2. The stilling basin of the barrage is of a modified form similar to the United States Bureau of Reclamation's (USBR's) stilling basin type-III.Within the basin, two rows of baffle and friction blocks facilitate energy dissipation and stabilize the HJ even in cases of minimum tailwater requirements.However, soon after the barrage construction, multiple problems such as the oblique right-sided river approach, heavy siltation in one of the major canals, uprooting of impact baffle blocks, damage to the stilling basin floor, and retrogression of water levels appeared in its downstream areas.During 1959During -1962During , 1966During , and 1973, repairs works were carried out to cater to the problems mentioned above; however, the problems remained persistent.In the reports by Zaidi et al. [58] and the World Bank [57], sweeping of the HJ was believed to be the main reason for the issues highlighted above.To resolve these issues, the Punjab Government constituted committees of experts; however, no specific measures were taken and the issues continued to be aggravated.The typical cross section of the Taunsa barrage stilling basin is shown in Figure 2.

Material and Methods
The methodology for the present study has two phases, as shown in Figure 3.

Material and Methods
The methodology for the present study has two phases, as shown in Figure 3.The stilling basin of the barrage is of a modified form similar to the United States Bureau of Reclamation's (USBR's) stilling basin type-III.Within the basin, two rows of baffle and friction blocks facilitate energy dissipation and stabilize the HJ even in cases of minimum tailwater requirements.However, soon after the barrage construction, multiple problems such as the oblique right-sided river approach, heavy siltation in one of the major canals, uprooting of impact baffle blocks, damage to the stilling basin floor, and retrogression of water levels appeared in its downstream areas.During 1959During -1962During , 1966During , and 1973, repairs works were carried out to cater to the problems mentioned above; however, the problems remained persistent.In the reports by Zaidi et al. [58] and the World Bank [57], sweeping of the HJ was believed to be the main reason for the issues highlighted above.To resolve these issues, the Punjab Government constituted committees of experts; however, no specific measures were taken and the issues continued to be aggravated.The typical cross section of the Taunsa barrage stilling basin is shown in Figure 2.

Material and Methods
The methodology for the present study has two phases, as shown in Figure 3.

Phase 1: Identification of the Critical Hydraulic Parameters
In the first phase, the critical hydraulic parameters (CHPs) were identified from indepth review of the published experimental and numerical data, whereas in the second phase the identified CHPs from the first phase were numerically investigated using FLOW-3D within the basin of Taunsa Barrage.For the first phase, a systematic review of the previous studies published from 2000 to 2020 was carried out and all those articles were retrieved that addressed the hydraulic parameters on downstream sides of any hydraulic structures such as spillways, barrages, sluice gates, weirs, and falls.The retrieved articles were divided into two categories, one was experimental studies and the other was numerical studies.By doing so, eighty research articles were retrieved, whereas the unrelated articles were discarded.A detailed review of the retrieved article was carried out and, skimming through these papers, forty-two articles on experimental studies and twenty-four articles on numerical studies were selected for further analysis to identify the CHPs.The statistical methods used for the identification of parameters are explained below.
The ranking of CHPs was performed on the basis of the relative percentage score, which was calculated by Formula (1) [59][60][61].
where R f and RPA are the relative frequency and relative portion of the party affected, respectively, which have been assigned to each parameter.In the present investigation, articles from the experimental and numerical categories were selected as parties [60][61][62][63][64][65].
The relative importance of the parameters can be calculated using the relative importance index (RII).A similar approach was used in the present study.The RII for each parameter was calculated using Formula (2): where RII is the relative importance index and W i and A are the weight and highest weight given to each parameter.In Formula (2), N is the total number of research articles from where these hydraulic parameters were taken.The RII value ranges from 0 to 1; the higher the RII, the more critical is the parameter.

Phase 2: Numerical Model Implementation
Choosing the most suitable CFD codes from many available options is crucial [63][64][65].Still, it is tedious as the investigative parameters [66][67][68][69] are strongly case dependent (Bennett et al. [70]).However, according to Chen et al. [14] and Babaali et al. [24], FLOW-3D has been the most widely used modelling tool for hydraulic investigations.Bayon et al. [50] recommended FLOW-3D for HJ characteristics after comparing its results with similar numerical models.Based on the recommendations of studies mentioned above, the present study implemented FLOW-3D numerical models to investigate the identified CHPs within the stilling basin of the studied barrage.FLOW-3D software (https://www.flow3d.com/) is considered as one of the most potent computer tools for performing three-dimensional (3D) flow analyses and uses multiple techniques to investigate the issue of multi-fluid by solving incompressible Reynolds-averaged Navier-Stokes (RANS) equations on each computational cell.FLOW-3D subdivides the flow and solid domains into structured grid blocks to resolve these domains for obtaining solutions of hydraulic issues.These structured rectangular grids are easy to develop and store essential information on cell faces and nodes.However, non-uniform mesh grid facilitates users to create meshes for complex geometries.Each cell within the grid is identified with a specific number in three dimensions: i in the x-direction, j in the y-direction, and k in the z-direction.The fundamentals of finite difference and finite volume methods were formerly developed on such meshes.These methods are central for the development of FLOW-3D.This 3D-modeling tool applies the finite volume method (FVM) derived directly from the conservation law to hold fluid properties.In FLOW-3D, the fluid-solid interface is tracked using the FAVOR method, whereas the generalized minimum residual method (GMRES) is implemented to solve issues of pressure velocity coupling.The proceeding sections show the equations used for the present models.
FLOW-3D discretizes the governing equations such as continuity and momentum equations.The general form of the mass continuity is described by Equation (3).For incompressible flow simulation, considering ρ as constant, Equation ( 3) is transformed for the incompressibility conditions, as provided in the following Equation ( 4): In Equations ( 3) and ( 4), V F is the partial volume of flow, ρ is the fluid density, R SOR is the mass source, and R DIF is the diffusion term of turbulence.In case of Cartesian coordinates, R is equal to unity and ξ is set as zero.The fluid velocity's components in three dimensions are computed using the following Reynolds-averaged Navier-Stokes (RANS) Equation ( 5): where u, v, and w are the velocity components, A x , A y , and A z are the flow areas, G x , G y , and G z are body accelerations, f x , f y , and f z are viscous accelerations, and ρ is the fluid density.

Model Meshing and the Initial and Boundary Conditions
The solid geometry of the model was prepared in AutoCAD and converted into a stereolithographic file (Stl.).Before importing the geometry into FLOW-3D, it was checked by Netfab software (https://inno-venture.com/netfabb/) to remove errors, holes, and facets.The model geometry and simulation domain were resolved using structured hexahedral mesh.
A single mesh block of 55.47 m long, 20.42 m wide, and 10.06 m high was implemented.Initially, a coarse mesh of 0.50 m cell size was applied to resolve the geometry; however, the stilling basin appurtenances were not fully resolved.Gradually reducing the cell size to 0.16 m, the geometry displayed a more suitable resolution to run the simulations.In total, 2,890,443 cubic mesh cells were used for models.To reduce the simulation time, a domain-removing component was added on the downstream side to deactivate the empty cells.The cell deactivation region of the domain-removing component was defined from the gate to the end of the stilling basin, and it was ensured that the region of domain-removing component did not contain the fluid.Figure 4 shows meshing applied to the models.
Routinely, barrage gates are not opened to the same levels.The openings are set according to the flows in the river.For reproducing the similar conditions for 44 m 3 /s of flow, the models were set for a constant elevation of 136.24 m for the pond levels, whereas five different tailwater levels ranging from 129.10 m to 130.30 m with an equal increment of 0.30 m were implemented on the downstream side.The flux surface (porosity = 1) was set at the end of the stilling basin and a movable probe was assigned to measure the free surface profile and other essential parameters in the stilling basin.The volume flow rate  2 shows the initial conditions for gatedand free-flow analysis.0.30 m were implemented on the downstream side.The flux surface (porosity = 1) was set at the end of the stilling basin and a movable probe was assigned to measure the free surface profile and other essential parameters in the stilling basin.The volume flow rate of single bay was used for the discharge calculation of 64 bays of the barrage, thereby the actual conditions for 44 m 3 /s flow were generated in the models.For all the models, the initial boundary conditions of discharge (44 m 3 /s), upstream pond level (136.24m), and turbulence model (RNG-K-ε) were kept constant, whereas the tailwater level was changed from 129.10 m to 130.30 m.The total simulation length of the model was 55.47 m, of which 38.10 m comprised the downstream area.Table 2 shows the initial conditions for gatedand free-flow analysis.Figure 5 shows that the pressure (P) boundaries were set for upstream (Xmin) and downstream (Xmax), whereas the lateral sides and bed were set as the rigid boundaries (W) and no-slip conditions were imposed which were expressed as zero tangential and normal velocity (u = v = w = 0) on the wall.u, v, and w are the velocity in the x, y, and z directions, respectively.For all variables (except pressure (P), which was set to zero), the upper boundaries (Zmax) were set as atmospheric pressure to allow water to null von Neumann.  Figure 5 shows that the pressure (P) boundaries were set for upstream (X min ) and downstream (X max ), whereas the lateral sides and bed were set as the rigid boundaries (W) and no-slip conditions were imposed which were expressed as zero tangential and normal velocity (u = v = w = 0) on the wall.u, v, and w are the velocity in the x, y, and z directions, respectively.For all variables (except pressure (P), which was set to zero), the upper boundaries (Z max ) were set as atmospheric pressure to allow water to null von Neumann.

Turbulence Modelling and Free Surface Tracking
One of the important aspects of computational fluid dynamic (CFD) models is turbulence closure.These numerical models implement Reynolds-averaged Navier-Stokes equations (RANS) to find the turbulence closure and solve high Reynolds numbers that develop flow instabilities.These models solve Reynolds stresses terms in the Navier-Stokes equation and calculate solutions for the additional equations of the turbulent viscosity and transport variable.Out of the six turbulence models in FLOW-3D, the two equation turbulence models (K-ε), standard K-ε and renormalization group (RNG K-ε), are the most widely used in hydraulic investigations.The above-mentioned turbulence models were applied by Bayon-Barrachina and Lopez-Jimenez [4], Macián-Pérez [9], and Bayon et al. [50] to investigate the HJ characteristics and the results showed that RNG K-ε produced reasonable accuracy for the efficiency of the HJ, sequent depths, roller lengths, and turbulent kinetic energy.Similarly, Nikmehr and Aminpour [71] also used RNG K-ε to investigate the free surface profile, flow rate, and Fr 1 on the corrugated bed and the model results were well in agreement with the compared experiments.Furthermore, studies such as those by Carvalho et [73] investigated HJ characteristics within the stilling basins of spillways and indicated that the RNG K-ε turbulence model produced free surface, velocity, pressure profiles, and turbulent kinetic energy that were well in agreement with the experimental results.Based on the published data from similar studies, the present studies implemented the RNG K-ε turbulence model within the stilling basin of the studied barrage.In the RNG K-ε turbulence model [73][74][75][76][77][78][79][80][81], Equations ( 6) and ( 7) were applied to model the turbulent kinetic energy (k) and its dissipation (ε), respectively.
In Equations ( 6) and (7), x i , µ, µ t , k, ε, ρ, and P k are the coordinate in x direction, dynamic viscosity, turbulent dynamic viscosity, turbulent kinetic energy (TKE), turbulent dissipation, fluid density, and produced TKE, respectively.Finally, the terms σk, σ ε , C 1ε , and C 2ε are model parameters whose values are given in the study by Yakhot et al. [67].The volume of fluid (VOF) method was used to track the free surface in which the fraction of the fluid (F) was implemented to find the fractional volume (i.e., water or air).To track the free surface within the simulation domain, Equation (8) was used.
where, in the modelling domain, the fluid fraction (F) is represented by the below three possibilities.

1.
If F approaches 0, the cell is considered as empty; 2.
When F reaches 1, the cell is believed to be occupied by fluid; 3.
If 0 < F < 1, the cell represents a surface between the two fluids.
Presently, one fluid (water) with free surface is considered, whereas other advection schemes are selected by the models.

Model Verification and Validation
For validation of free designed flow analysis, h e /H d = 0.998 was implemented as previously used in [72, 73,82,83], whereas h e and H d were the effective head and designed heads, respectively, as shown in Figure 6.For the designed flow analysis to run the simulations, pond and tailwater levels of 135.93 m and 133.8 m were used, respectively.

Model Verification and Validation
For validation of free designed flow analysis, he/Hd = 0.998 was implemented as previously used in [72, 73,82,83], whereas he and Hd were the effective head and designed heads, respectively, as shown in Figure 6.For the designed flow analysis to run the simulations, pond and tailwater levels of 135.93 m and 133.8 m were used, respectively.On the other hand, for gated flow, Formula (9) was used to calculate the discharge through the orifice.For 44 m 3 /s discharge, D = 0.280 m and Hd = 136.24m were used to operate the models, where D is the gate opening and Hd is the design head for orifice discharge.
where Q, A, and g are the volume flow rate, area of orifice, and acceleration due to gravity, respectively, and are measured in m 3 /s, m 2 , and m/s 2 , respectively.However, in Equation ( 9), hc is the centerline head, which is calculated using the gate openings and pond levels.
For gated flow, a 0.816 value for the coefficients of discharge (Cd) was used, whereas a Cd value of 0.819 was obtained from the models.
Courant number stability criteria [66,67] were adopted to control the time steps and varied from 0.06 to 0.0023 and 0.015 to 0.0025 for free and gated flow, respectively.The volume flow rates at inlet and outlet boundaries were monitored to check the steady state of the models [83][84][85].For modelling the discharges of 44 m 3 /s and 444 m 3 /s, a simulation finish time of (T = 80 s) was selected, as can be seen in Figure 7; the models achieved steady state at T = 60 s and T = 75 s for free designed and gated flow, respectively.On the other hand, for gated flow, Formula (9) was used to calculate the discharge through the orifice.For 44 m 3 /s discharge, D = 0.280 m and H d = 136.24m were used to operate the models, where D is the gate opening and Hd is the design head for orifice discharge.
where Q, A, and g are the volume flow rate, area of orifice, and acceleration due to gravity, respectively, and are measured in m 3 /s, m 2 , and m/s 2 , respectively.However, in Equation ( 9), h c is the centerline head, which is calculated using the gate openings and pond levels.For gated flow, a 0.816 value for the coefficients of discharge (C d ) was used, whereas a C d value of 0.819 was obtained from the models.
Courant number stability criteria [66,67] were adopted to control the time steps and varied from 0.06 to 0.0023 and 0.015 to 0.0025 for free and gated flow, respectively.The volume flow rates at inlet and outlet boundaries were monitored to check the steady state of the models [83][84][85].For modelling the discharges of 44 m 3 /s and 444 m 3 /s, a simulation finish time of (T = 80 s) was selected, as can be seen in Figure 7; the models achieved steady state at T = 60 s and T = 75 s for free designed and gated flow, respectively.Analysis of free designed and flow analysis showed that at the beginning the free surface on the upstream and downstream of weir was found to be fluctuating; this became stable when the models reached the steady state, as shown in Figure 7. Free stable HJs were observed at T = 60 s and T = 75 s for free and gated flow, respectively.In free flow, the model underestimated the flow and the maximum error reached −5%, whereas a 1.14% error was observed in gated flow.From the gated and free flow analysis, it was found that the present models produced acceptable discharge accuracy, which allowed us to study the CHPs downstream of the studied barrage as described in Section 4.2.

Identification of Critical Hydraulic Parameters (CHPs)
Appendices A and B (Supplementary data) show the sources of hydraulic parameters extracted from the numerical and experimental studies.The parameters in the appendixes Analysis of free designed and flow analysis showed that at the beginning the free surface on the upstream and downstream of weir was found to be fluctuating; this became stable when the models reached the steady state, as shown in Figure 7. Free stable HJs were observed at T = 60 s and T = 75 s for free and gated flow, respectively.In free flow, the model underestimated the flow and the maximum error reached −5%, whereas a 1.14% error was observed in gated flow.From the gated and free flow analysis, it was found that the present models produced acceptable discharge accuracy, which allowed us to study the CHPs downstream of the studied barrage as described in Section 4.2.

Identification of Critical Hydraulic Parameters (CHPs)
Appendices A and B show the sources of hydraulic parameters extracted from the numerical and experimental studies.The parameters in the appendixes were based on their frequency of occurrence.From the initial analysis, twenty-four and twenty-three hydraulic parameters were found from the numerical and experimental studies, respectively, as provided in Tables 3 and 4, respectively.However, fifteen parameters were found to be common to both types of studies (numerical and experimental).By taking the union of both sources (numerical and experimental), a total of thirty-three parameters were identified from the published data; these are listed in Table 5.In Table 3, by applying the relative importance index (RII), the velocity profile (VP), free surface profile (FSP), pressure profile (PP), turbulence kinetic energy (TKE), air volume value (AV), discharge measurement (DM), shape of stilling basin (SS), Fr 1 , and HJ efficiency (ï) were ranked as the most CHPs in numerical studies that were conducted from 2000 to 2020.
Table 4 shows the ranking of the parameters that were computed for experimental studies.The RII showed that VP, Fr 1 , TWL, SS, FSP, BP, SP, and ED were the most CHPs, upon which several studies have been performed during recent years.After adding up the parameters from the experimental and numerical studies based on their RII, the VP, FSP, Fr 1 , SS, TWL, PP, BP, TKE, SP, and ED were found to be the most important CHPs, as shown in Table 5.In Table 6, the relative percentage score of the 33 parameters was calculated; to do so, the Rf and RPA for each parameter were computed and results were presented in three different ranks.Table 7 shows the overall relative position of the CHPs based on frequency analysis, RII, and relative % score.It can be seen from Table 7 that, except for VP, all the other parameters changed their position when employing different statistical methods.From Table 7, it is found that VP, Fr, FSP, SS, TKE, and TWL are the most significant parameters that have been widely investigated in the literature.Hence, for the present study, these hydraulic parameters are focused on and investigated downstream of the studied barrage.

Parameters
Ranking Extracted from Table 2 Ranking Extracted from Table 3 Ranking Extracted from Table 4 Ranking Extracted from Table 5 Overall Occurrence (1) (2)  8 compares the free surface profiles at five different TWLs.The free surface data were taken from the downstream of the gate to the end of the stilling basin.The total upstream ponding length was 16.45 m, which ranged from X min to the gate.At start, with the time rate of change, the free surface profiles were fluctuating [86][87][88][89]; however, they became stable when the models achieved a steady state.For the investigated tailwater levels, HJs were located at the glacis.However, the locations and lengths of the HJs were found to be different as the tailwater levels were changed; the shapes of the HJs were also changed.The results further indicated that, at low TWLs, the undulating free water surface was found in the HJ, which continued to the end of the stilling basin, as shown in Figure 8 in 129.10 m and 129.40 m TWLs.At higher TWLs, the HJ was shifted to the upstream side of the glacis and different free surface profiles were noticed compared with the lower tailwater levels, as shown in Figure 8.
At 129.10 m tailwater, the results showed a 127.90 m elevation for the HJ; this is in good agreement with [82,90], in which the HJ elevation was observed at 128 m.Furthermore, the results of the present models are also compared with the designed and prototype data for the year 2010.In these reports, the distance of the HJ at 130.30 m was 1.22 m from the toe of the glacis, as shown in Table 8.However, the HJ locations were found to be missing on the lower tailwater [91][92][93][94].Upon comparison with the available data point, the present model showed a five times higher HJ distance from the toe of the downstream glacis.In addition, even at a lower TWL of 129.10 m, the distance of the HJ was found to be 2.5 times higher than that observed at the barrage site during 2010.A detailed comparison of the HJ elevation and its location is provided in Table 8.From Table 8, it can be realized that compared with the present models, the HEC-RAS models in [36] overestimated the length and location of the HJs.On the contrary, The results further indicated that, at low TWLs, the undulating free water surface was found in the HJ, which continued to the end of the stilling basin, as shown in Figure 8 in 129.10 m and 129.40 m TWLs.At higher TWLs, the HJ was shifted to the upstream side of the glacis and different free surface profiles were noticed compared with the lower tailwater levels, as shown in Figure 8.
At 129.10 m tailwater, the results showed a 127.90 m elevation for the HJ; this is in good agreement with [82,90], in which the HJ elevation was observed at 128 m.Furthermore, the results of the present models are also compared with the designed and prototype data for the year 2010.In these reports, the distance of the HJ at 130.30 m was 1.22 m from the toe of the glacis, as shown in Table 8.However, the HJ locations were found to be missing on the lower tailwater [91][92][93][94].Upon comparison with the available data point, the present model showed a five times higher HJ distance from the toe of the downstream glacis.In addition, even at a lower TWL of 129.10 m, the distance of the HJ was found to be 2.5 times higher than that observed at the barrage site during 2010.A detailed comparison of the HJ elevation and its location is provided in Table 8.From Table 8, it can be realized that compared with the present models, the HEC-RAS models in [36] overestimated the length and location of the HJs.On the contrary, even after the remodeling of basin, the location and distances of HJs from the downstream toe of the glacis were found to be less than in the old basin, i.e., studied presently.

Froude Number
Figure 9 shows Froude number variations in the stilling basin.In the tested models, large numbers of oscillations were observed from the gate opening to the jump initiating point.The maximum value for Fr 1 was found for 129.10 m tailwater, which reached to 5.87, whereas the minimum value for Fr 1 was 5.30 at 129.70 m.
Fluids 2023, 8, x FOR PEER REVIEW 15 of 30 point.The maximum value for Fr1 was found for 129.10 m tailwater, which reached to 5.87, whereas the minimum value for Fr1 was 5.30 at 129.70 m.As compared with [36] for all the investigated tailwater levels, the present models showed higher values for Fr1, whereas in the subcritical region, the results for the Froude number agreed with the study of Chaudhry [36].A gradual decrease in the Fr1 was observed when the TWLs increased.The minimum value for Fr1 was observed at the maximum TWL.After the HJ, the flow changed into the subcritical state and the maximum value of Fr2 was 0.22 at a 129.10 m TWL.At a constant discharge, pond level, and gate opening, the results of TWLs against Fr1 showed a nonlinear trend.Two-dimensional illustrations of Fr in the stilling basin are shown in Figure 10a-e.As compared with [36] for all the investigated tailwater levels, the present models showed higher values for Fr 1 , whereas in the subcritical region, the results for the Froude number agreed with the study of Chaudhry [36].A gradual decrease in the Fr 1 was observed when the TWLs increased.The minimum value for Fr 1 was observed at the maximum TWL.After the HJ, the flow changed into the subcritical state and the maximum value of Fr 2 was 0.22 at a 129.10 m TWL.At a constant discharge, pond level, and gate opening, the results of TWLs against Fr 1 showed a nonlinear trend.Two-dimensional illustrations of Fr in the stilling basin are shown in Figure 10a-e.As compared with [36] for all the investigated tailwater levels, the present models showed higher values for Fr1, whereas in the subcritical region, the results for the Froude number agreed with the study of Chaudhry [36].A gradual decrease in the Fr1 was observed when the TWLs increased.The minimum value for Fr1 was observed at the maximum TWL.After the HJ, the flow changed into the subcritical state and the maximum value of Fr2 was 0.22 at a 129.10 m TWL.At a constant discharge, pond level, and gate opening, the results of TWLs against Fr1 showed a nonlinear trend.Two-dimensional illustrations of Fr in the stilling basin are shown in Figure 10a-e.

Flow Depths
Figure 11 shows that at all the investigated tailwater levels, the flow depths up to the jump initiating point show identical behavior, whereas fluctuations in the flow depths were observed in the HJ region.Except at a 129.99 m TWL, the results indicate the smooth transition from supercritical (y 1 ) into subcritical flow depths (y 2 ), whereas at a 129.99 m TWL, large fluctuations were seen in the jump.At a 129.10 m TWL, due to the presence of friction blocks, the results indicate small oscillations in the flow depths at the stilling basin's end that deflected the flow towards the free surface.

Flow Depths
Figure 11 shows that at all the investigated tailwater levels, the flow depths up to the jump initiating point show identical behavior, whereas fluctuations in the flow depths were observed in the HJ region.Except at a 129.99 m TWL, the results indicate the smooth transition from supercritical (y1) into subcritical flow depths (y2), whereas at a 129.99 m TWL, large fluctuations were seen in the jump.At a 129.10 m TWL, due to the presence of friction blocks, the results indicate small oscillations in the flow depths at the stilling basin's end that deflected the flow towards the free surface.From Figure 12, it can be observed that as the tailwater was increased, the sequent depth ratio increased, and the Froude number values were found to be decreased.Figure 12 also compares the sequent depths of present study with the previous experimental and numerical studies.It can be seen from Figure 12 that at lower tailwater levels, i.e., 129.10 m, the sequent depths agreed well with the experimental data from Kucukali and Chanson [94] and numerical study by Bayon-Batrachia and Lope-Jimenez [4].However, as the tailwater levels increased, the sequent depths showed deviation from the compared studies.From Figure 12, it can be observed that as the tailwater was increased, the sequent depth ratio increased, and the Froude number values were found to be decreased.Figure 12 also compares the sequent depths of present study with the previous experimental and numerical studies.It can be seen from Figure 12 that at lower tailwater levels, i.e., 129.10 m, the sequent depths agreed well with the experimental data from Kucukali and Chanson [94] and numerical study by Bayon-Batrachia and Lope-Jimenez [4].However, as the tailwater levels increased, the sequent depths showed deviation from the compared studies.

Flow Depths
Figure 11 shows that at all the investigated tailwater levels, the flow depths up jump initiating point show identical behavior, whereas fluctuations in the flow d were observed in the HJ region.Except at a 129.99 m TWL, the results indicate the sm transition from supercritical (y1) into subcritical flow depths (y2), whereas at a 129 TWL, large fluctuations were seen in the jump.At a 129.10 m TWL, due to the prese friction blocks, the results indicate small oscillations in the flow depths at the stillin sin's end that deflected the flow towards the free surface.From Figure 12, it can be observed that as the tailwater was increased, the se depth ratio increased, and the Froude number values were found to be decreased.F 12 also compares the sequent depths of present study with the previous experiment numerical studies.It can be seen from Figure 12 that at lower tailwater levels, i.e., m, the sequent depths agreed well with the experimental data from Kucukali and son [94] and numerical study by Bayon-Batrachia and Lope-Jimenez [4].However, tailwater levels increased, the sequent depths showed deviation from the compared ies.     Figure 14 shows the behavior of velocity vectors in the longitudinal direction at various tailwater levels.Due to the supercritical velocity, the contracted flow jet was impinging at the toe of the glacis, which decreased the velocities at the upper fluid regions.Before the stilling basin appurtenances, reverse flows and eddies can be seen in the HJ at all the investigated TWLs.The flow behaviors of the upper fluid region of the HJ indicated typical backward velocity profiles, as described by Ead and Rajaratnam [11][12][13].As the flow reached a steady state, these reverse fluid circulations stabilized, which showed the stagnation zones.The analysis showed that the recirculation region occurred in the HJ and between and after the baffle blocks, and the maximum backward velocity profiles were found in the developed areas.At all the investigated TWLs, after the HJ the flow recovering zone starts and the negative velocity profile becomes positive.The results showed that as the tailwater levels increased, the velocity values in the HJ and in the subcritical region decreased.Compared with lower tailwater levels, at a constant discharge, the velocity decay in the HJ region was found more at higher tailwater levels.Additionally, as compared with lower tailwater levels, at higher tailwater levels, the values of the velocities were small at the end of the end stilling basin.As compared with the numerical study by Chaudhry [36], at all the studied tailwater levels, the observed velocity values in the present numerical study were found to be higher.
Figure 14b shows depth-averaged velocity profiles at the studied TWLs.These velocity profiles are drawn from the centerline of the bay.At all the TWLs, the maximum velocity was found just before the initial locations of the HJs, which reached 9.49, 9.45, 9.30, 9.38, and 9.28 m/s at TWLs of 129.10, 129.40, 129.70, 129.99, and 130.30, respectively.In the transition regions of the HJ, due to the eddies and fluid recirculation, the velocity rapidly decreased and remained consistent in the baffle block region.After the baffle blocks, a slight increase in the velocity values was noticed, from X = 57 m to 62 m.From the jump initiating location to the start of friction blocks region, the maximum velocity was noticed at the 129.40 m TWL; however, after the friction blocks, the velocity values became equivalent to those that were observed at the 129.10 TWL.In addition, in and after the friction blocks region, the velocity values at the studied TWLs were further reduced.The results further showed that as the TWL was increased, the velocity values in the stilling basin were found to be decreased due to the higher tailwater depths.At the basin's end, the maximum and minimum velocity values reached 1.1 m/s and 0.70 m/s at the 129.10 and 130.30 TWLs, respectively.Additionally, the effect of energy dissipation devices shows that after the baffle blocks, the velocity profiles near the bed decreased and became positive on the free surface.In addition, at all the studied tailwater levels the vertical velocity profiles followed the trend of Ead and Rajaratnam [11][12][13].The maximum velocities were found at the floor level before the baffle blocks and decreased as the flow moved forward to the end of the stilling basin.Figure 14 shows the behavior of velocity vectors in the longitudinal direction at various tailwater levels.Due to the supercritical velocity, the contracted flow jet was impinging at the toe of the glacis, which decreased the velocities at the upper fluid regions.Before the stilling basin appurtenances, reverse flows and eddies can be seen in the HJ at all the investigated TWLs.The flow behaviors of the upper fluid region of the HJ indicated typical backward velocity profiles, as described by Ead and Rajaratnam [11][12][13].As the flow reached a steady state, these reverse fluid circulations stabilized, which showed the stagnation zones.The analysis showed that the recirculation region occurred in the HJ and between and after the baffle blocks, and the maximum backward velocity profiles were found in the developed areas.At all the investigated TWLs, after the HJ the flow recovering zone starts and the negative velocity profile becomes positive.
Additionally, the effect of energy dissipation devices shows that after the baffle blocks, the velocity profiles near the bed decreased and became positive on the free surface.In addition, at all the studied tailwater levels the vertical velocity profiles followed the trend of Ead and Rajaratnam [11][12][13].The maximum velocities were found at the floor level before the baffle blocks and decreased as the flow moved forward to the end of the stilling basin.

Turbulent Kinetic Energies (TKEs)
The root mean square values of velocity fluctuations were used to calculate turbulent kinetic energies (TKE) at various locations.By considering the successive velocity values, the value of root mean square velocity (u rms ) can be obtained using Equation (10).
In the above Equation ( 10), u 1 , u 2 , and u 3 are the successive velocity values in the horizontal direction.Using the velocity values, the TKE can be calculated using Equation (11).
whereas u rms , v rms , and w rms are the root mean square velocities in the x, y, and z directions, respectively.
Figure 15a shows the variation in the TKE in 129.10 m tailwater.The TKE values were computed from HJ initiation to the termination points, such as at the supercritical flow region and within the HJ and subcritical regions.The maximum value for the TKE was found after the HJ and decreased afterward.At a 129.10 m tailwater level, the maximum value for the TKE reached 3.72 m 2 /s 2 at X =10.At all the fluid depths at X = 52, after the HJ the TKE values were found to level off, as shown in Figure 15a.The result further showed that the TKEs were found to be maximal in the middle regions of the flow depths and that that minimum values were found at the free surface.In Figure 15b, as in the lower tailwater levels, the maximum value for the TKE was found after the HJ initiating location, and the value reached 3.12 m 2 /s 2 .However, as compared with lower tailwater levels, the values for TKEs in 129.40 m tailwater were found to be lower and the TKEs levelled off earlier within and after the HJ.In the HJ regions, following the similar trends for the 129.10 and 129.40 m tailwater levels, the TKEs in the 129.70, 129.99, and 130.30 m tailwater levels reached 3.20, 3.15, and 2.90 m 2/ s 2 , respectively, and their values decreased after the HJ, as illustrated in Figure 15c-e.Additionally, the results for the TKEs showed a trend that was noted in the numerical studies (Nikmehr and Aminpour [71]; Soori et al. [76]).Figure 16 shows 2D illustrations of turbulent kinetic energies captured using RNG-K-ε at five different tailwater levels.Figure 16a shows that in 129.10 m tailwater, the maximum turbulent kinetic energy was found before and within the HJ.The rest of the kinetic energy after the jump was dissipated by two rows of baffle blocks.After the baffle blocks, the flow became less turbulent and TKE was reduced and gradually leveled off at the end of stilling basin.Similarly, in the case of a 129.40 m tailwater level, the overall values for the TKE decreased compared with the 129.10 m tailwater.However, the maximum TKE is found before and within the HJ, as shown in Figure 16b.After the baffle blocks, TKEs gradually reduced up to the end of stilling basin.Following the same development, at the 129.70 m, 129.99 m, and 130.30 m TWLs, the maximum TKEs are found before and within the HJ; this can be seen in Figure 16c-e.
The 2D illustrations further indicate that as the tailwater increased downstream of the barrage, the TKEs within the HJ decreased but the turbulent behavior of all the fluid layers was observed differently.At the 130.30 m tailwater level, the TKEs within the upper and lower layers of the HJ decreased compared with the central region, as illustrated in Figure 16e.This different behavior within the HJ could be due to the greater flow depth available for the jump.The 2D illustrations further indicate that as the tailwater increased downstream of the barrage, the TKEs within the HJ decreased but the turbulent behavior of all the fluid layers was observed differently.At the 130.30 m tailwater level, the TKEs within the upper and lower layers of the HJ decreased compared with the central region, as illustrated in Figure 16e.This different behavior within the HJ could be due to the greater flow depth available for the jump.

Discussion and Real-World Implications
The preceding sections, Sections 4.1 and 4.2, focused on CHP identification and the effects of these parameters downstream of the Taunsa barrage, respectively.Frequency

Discussion and Real-World Implications
The preceding sections, Sections 4.1 and 4.2, focused on CHP identification and the effects of these parameters downstream of the Taunsa barrage, respectively.Frequency analysis, RII, and relative % score revealed different rankings for the hydraulic parameters, which were focused in many experimental and numerical studies.However, out of the CHPs, the velocity profile and Froude number were found to be the most investigated hydraulic parameters, as can be seen in Table 6.On the contrary, after extensive analysis of the literature, the investigation of the identified CHPs is found to be lacking on the studied barrage.The literature only revealed one study that employed a one-dimensional (1D) HEC-RAS model to investigate the limited hydraulic parameters downstream of the studied barrage [36].Therefore, the present FLOW-3D model results are compared with 1D and the relevant studies from the literature, i.e., Bayon-Barrachina and Lope-Jimenez [4], Chaudhry [36], Ead and Rajaratnam [11][12][13], Nikmehr and Aminpour [71], Soori et al. [76], and Kucukali and Chanson [94].After comparing with the results of the free surface profiles of HJs at the studied tailwater levels, it is revealed that the previously used 1D model underestimates the locations of the HJs; for which the errors reached 8%.Similarly, as compared with the present 3D models, the values for the Froude number in the supercritical region were found to be lower in 1D model, whereas the models' results for the Froude number in the subcritical region agree with previous data [36].At all the tested tailwater levels, the analysis of sequent depths ratios showed agreement with the previous studies.In comparison with the 1D study [36], the results of velocity profiles within the stilling basin were found to be higher in the present models.However, at all the tested tailwater levels, the results of the vertical velocity profiles in the HJ region showed wall-jet-like profiles that showed agreement with Ead and Rajaratnam [11][12][13].The maximum TKEs were found within the HJ region, and these declined as the distance from the HJ increased.In addition, the maximum amount of TKE was found within the central fluid depth; this started to decline towards the free surface and the basin's bed.
The findings of the present study will facilitate both hydraulic researchers and practitioners.Firstly, the results identified critical hydraulic parameters that should be given significant importance when a new hydraulic intervention is to be carried out, i.e., remodeling of hydraulic structures (Chaudary and Sarwar [35]; World Bank [57]; Zaidi et al. [58]).On the other hand, the use of numerical models is becoming prevalent in hydraulic investigations such as those carried out in the present study.The results showed that previous one-dimensional HEC-RAS studies are limited and unable to describe the flow characteristics spatially; this is why the FLOW-3D models' results are found to be different.In addition, the results from the FLOW-3D model also highlight that tailwater levels before the remodeling of the barrage were appropriate to hold the HJs at the glacis well above its toe.Therefore, it is believed that, in future, these modeling tools will also eliminate physical modeling, because conventionally scaled modeling is usually associated with the difficulties of terrain, concrete roughness, and flow measuring devices.The present study was limited to a single discharge value and turbulence model.Therefore, the model should also be tested and evaluated for higher values of discharge and other turbulence models.

Conclusions
The present study identified critical hydraulic parameters (CHPs) from the literature and studied these parameters downstream of Tuansa barrage using FLOW-3D numerical models.The study also investigated the effects of changes in tailwater levels on the HJ characteristics and compared the results with available previous data for the studied barrage.Following main conclusions drawn are:

•
The literature review outlined thirty-three hydraulic parameters; out of those, the velocity profile, Froude number, free surface profile, shape of stilling basin, tailwater, and turbulent kinetic energy were the highly significant hydraulic parameters in the literature that were studied downstream of hydraulic structures.

•
At all the investigated tailwater levels, no sweeping of the HJ was observed as reported in the previous studies.The location and elevation of HJs were observed to be different compared with a previous HEC-RAS one-dimensional hydraulic study.Upon comparison with the HJ results of the designed and downstream of a prototype barrage (i.e., remodeled basin), the distance of the HJ from the glacis toe was found to be higher, which further revealed the old basin (i.e., studied presently) was efficiently holding the HJ at the investigated discharge and TWLs.

•
Non-linear trends for the Froude number and sequent depths were observed as the tailwater levels varied.On comparison with previous studies, the present models showed higher values for the Froude number and sequent depths, which showed deviation at higher tailwater levels.
• At lower tailwater levels, the vertical velocity profiles in the developing region of the HJ near the floor were found to be higher than the results at higher tailwater levels.At the investigated tailwater levels, jet-like velocity profiles were obtained in the HJ regions that levelled off as the distance from the HJ was increased.

•
The maximum turbulent kinetic energy was found in the developing region of the HJ at the minimum tailwater level.After the impact and baffle blocks, the kinetic energy gradually reduced, and the minimum kinetic energy was observed at higher tailwater levels.
Based on the results of the numerical models, it can be said that the tailwater envelope for the studied discharge (44 m 3 /s) for the Tuansa Barrage was within the acceptable limit, which holds the HJ on the downstream glacis.Using FLOW-3D numerical models, the study further confirmed that the previously used one-dimensional HEC-RAS models produced higher values for the hydraulic parameters.As the present study was limited to one discharge, more comprehensive studies for other discharges and turbulence models must be carried out on the stilling basin of the studied barrage.The study further suggests investigating the retrogression analysis downstream of the Taunsa barrage using threedimensional numerical models.

Figure 1 .
Figure 1.Location of the study area.Figure 1. Location of the study area.

Figure 1 .
Figure 1.Location of the study area.Figure 1. Location of the study area.

Figure 2 .
Figure 2. Typical cross section of the barrage.

Figure 2 .
Figure 2. Typical cross section of the barrage.

Figure 2 .
Figure 2. Typical cross section of the barrage.

Figure 3 .
Figure 3. Methodology of the study.

Fluids
was used for the discharge calculation of 64 bays of the barrage, thereby the actual conditions for 44 m 3 /s flow were generated in the models.For all the models, the initial boundary conditions of discharge (44 m 3 /s), upstream pond level (136.24m), and turbulence model (RNG-K-ε) were kept constant, whereas the tailwater level was changed from 129.10 m to 130.30 m.The total simulation length of the model was 55.47 m, of which 38.10 m comprised the downstream area.Table

Figure 4 .
Figure 4. Meshing setup for the simulation domain.

Figure 4 .
Figure 4. Meshing setup for the simulation domain.

Fluids 30 Figure 5 .
Figure 5. Boundary conditions governing the models.3.2.2.Turbulence Modelling and Free Surface Tracking One of the important aspects of computational fluid dynamic (CFD) models is turbulence closure.These numerical models implement Reynolds-averaged Navier-Stokes equations (RANS) to find the turbulence closure and solve high Reynolds numbers that

Figure 5 .
Figure 5. Boundary conditions governing the models.

Fluids 2023, 8 , 30 Figure 7 .
Figure 7. Time rate of change of flow at inlet and outlet boundaries.(a) Gated flow.(b) Free designed flow.

Figure 7 .
Figure 7. Time rate of change of flow at inlet and outlet boundaries.(a) Gated flow.(b) Free designed flow.

Figure 8 .
Figure 8. Free surface profiles and location of HJs at various tailwater levels.

Figure 9 .Figure 9 .
Figure 9. Variation in Froude number in the basin at different tailwater levels.

Fluids 2023, 8 ,
x FOR PEER REVIEW 15 of 30 point.The maximum value for Fr1 was found for 129.10 m tailwater, which reached to 5.87, whereas the minimum value for Fr1 was 5.30 at 129.70 m.

Figure 9 .Figure 10 .
Figure 9. Variation in Froude number in the basin at different tailwater levels.

Figure 11 .
Figure 11.Flow depth variation at different TWLs in the stilling basin.

Figure 11 .
Figure 11.Flow depth variation at different TWLs in the stilling basin.

Figure 11 .
Figure 11.Flow depth variation at different TWLs in the stilling basin.

Figure 12 .
Figure 12.Comparison of sequent depth with the relevant literature studies [2,94] at various tailwater levels.
Figure13a,b shows the free surface and depth-averaged velocities at various tailwater levels in the studied basin.The maximum values for longitudinal velocity were observed in the supercritical region before the HJ.Due to the recirculation and turbulence in the HJ region, negative velocity was observed.The maximum negative velocities were observed at lower tailwater levels, as shown in Figure13a.

Figure 13 .
Figure 13.Velocity profiles in the stilling basin at the investigated TWLs.(a) Free−surface velocity profiles, (b) depth−averaged velocity.
highlights a few of the most relevant studies dealing with HJs and energy dissipators.

Table 2 .
Initial conditions for the gated and free flow.

Table 2 .
Initial conditions for the gated and free flow.

m 3 /s) Minimum Tailwater Required for HJ (m) Maximum Tailwater Required for HJ (m) Upstream Water Level Maintained (m) Turbulence Model Models Operation
al. [56], Savage and Johnson [72], and Johnson and Savage

Table 3 .
Frequency of occurrence, RII, and ranking of parameters from numerical studies.

Table 4 .
Frequency of occurrence, RII, and ranking of parameters from experimental studies.

Table 5 .
RII and ranking of parameters by combining numerical and experimental studies.

Table 7 .
Relative position of CHPs in the literature with different statistical methods.

Table 8 .
Comparison of the results of the hydraulic jump with HEC-RAS and field data.

Table 8 .
Comparison of the results of the hydraulic jump with HEC-RAS and field data.