Flow pattern study in Beshar River and its two straight and meander reaches using CCHE2D model

In this study, flow pattern in Beshar River as a main branch of Karoon River has been analyzed using CCHE2D. For this purpose, a 12-km reach upstream of Shahmokhtar hydrometric station near Yasuj city was considered. The CCHE2D model was calibrated using different Manning’s roughness coefficients and different turbulent models; for this purpose, numerical results were compared with observation data for three different discharges. The results showed that for the medium and high discharges, less Manning’s roughness coefficients (0.015 ≥ n ≥ 0.025) and for low discharge, higher Manning’s roughness coefficients (0.035 ≤ n ≤ 0.050) are more suitable. Also, k–ε turbulent model is more effective in this study. Besides, variations of hydraulic parameters like water depth, velocity, shear stress and Froude number are calculated and discussed. The analysis of the flow and velocity pattern in the straight and meander reaches of the river shows that the changes trend of the water surface gradient and velocity in the cross sections of this two reaches are different. Due to effect of secondary currents, latitude gradient of the water surface and depth average velocities increase to the outer bank of the bend. But in the straight reach, latitude gradient of the water surface is almost zero and the maximum velocities are in the center-line of flow. The R-squared (RSQ) and linear correlation coefficient (r) factors between velocity and shear stress show that there is linear and direct relationship between these two hydraulic parameters in the entire study reach.


Introduction
A river carries water, sediment and solute, and this is important to hydraulic engineers, geomorphologists and sedimentologists (Johnson 2008). In addition, morphological properties of rivers are considered in river management, and different forms of rivers require different approaches to organizing (Radan and Vaghefi 2016). Therefore, hydrological studies of rivers are of great importance for a variety of construction, housing, etc. Several procedures exist in river hydraulic studies such as field measurements, numerical model, physical model and a combination of these methods, for example, a numerical model along with a field approach. Modern hydrological models are typical of process-based environmental models (Hu and Bian 2009). In some models, such as physical model, a large number of physical parameters are needed for accurate and efficient model simulations (Kumar et al. 2010), and therefore, the most accurate model should be selected.
Modeling is advantageous because: it offers real-scale spatially dense results; it allows for the simulation of past, current and future conditions governed by multiple boundary conditions; it provides significant cost advantages over the production of a physical model; and it provides a tool to address questions that have previously been restricted by time, scale and resources (Garde 2005). Therefore, spatial discretization of the landscape is usually accomplished using a fully or semi-distributed into modeling units (Francke et al. 2008).
In this regards, mathematical models have been widely used in solving complex hydrodynamic topics in a wide range of scales. In a hydraulic structure, i.e., canal, mathematical models can be used to check stability under different flow conditions. The most common type of numerical models used by engineers are hydrodynamic models. These models analyze the flow conditions using numerical methods. The simulation of flow characteristics requires at least a two-dimensional model (Shahidan et al. 2012). Recently, numerous 2D and 3D numerical models have been developed for modeling the flow and processes of sediment transport. Most 2D models have lower computing costs than 3D models. So on many issues, 2D models are used (Nassar 2011). However, analyzing and modeling some issues, such as the simulation of horse shoe vortex, local scour hole around piers and abutment of bridges, need to use 3D models, and the use of 2D numerical models in such cases is not logical (Zhang 2006a). CCHE2D model has been used for simulation at the Nile River (Elbogdady reach), and multiparametric sensitivity with different roughness parameters was evaluated using RSQ and r factors (Nassar 2011). Mohanty et al. (2012) studied wide meandering compound laboratory channel using CCHE2D, and then the results were compared with the observed values for validation purpose. In other study, sediment transport in Karkheh River was simulated by Kamandbedast et al. (2013) and the results showed that erosion and sedimentation were the dominant phenomena for floods 50 and 25 years, respectively. Elyasi and Kamandbedast (2014) studied the numerical modeling of flow in a river with 90-degree angle bend using the CCHE2D model.
Due to flow in a meander reach of the Palmanaki River (sub-Arctic Northern Finland), flow pattern and morphological changes were studied using the new acoustic Doppler current profiler (ADCP) method (Kasvi et al. 2017). According to the results which were done based on hydraulic parameters of velocity and water depth, duration of flood and the rate of increase and decrease in discharge have played an important role in determining channel changes by controlling the velocity and depth of flow. In another research, hydraulic parameters of flow depth and velocity were investigated in the Karoon River (Iran country; Khuzestan province) using CCHE2D model. Finally, it was clear that the changes of depth and flow velocity in the river have a fine harmony, and most changes were in the meander and the bends (Yusefi Haghivar et al. 2017). Yang et al. (2017) developed a 2D model for simulating the braiding phenomena and morphological changes in rivers. Their model uses hydrodynamic rules, sediment transport equations and total variation diminishing (TVD) method to forecast flows and changes in bed morphology. The model was used for simulating the process of bed changes in a physical model of river with bed load transport. Increases in the active and total braiding intensity showed the same trend to those observed in the experimental river. Gharbi et al. (2016) used a 2D model to predict the amount of materials transported by the Medjerda River in Tunisia. The purpose was to investigate Medjerda behavior in severe accidents and morphological changes after passing through spectacular floods in January 2003. In their research, a comparative analysis was carried out between the 1D model (HECRAS), and the 2D model (TELEMAC) coupled with SISYPHE. The results showed that the 2D model can calculate flow changes, morphological changes and sediment transport rates in severe incident for complex natural range with high precision compared to the 1D model. ShahiriParsa et al. (2016) investigated flood zoning using one-dimensional model (Hydrologic Engineering Centers-River Analysis System, HEC-RAS) and two-dimensional model with CCHE2D sofware to simulate the flood pattern in the Sungai Maka district in Kelantan state, Malaysia. The results showed that the maximum difference between the 1D and 2D models is 6% in the meander's part of the river (ShahiriParsa et al. 2016). In addition, the results of flow pattern simulation at a meandering reach with CCHE2D model in the Khoshke-rud River of Iran showed that using computational fluid dynamics for water flow modeling is one step closer to having a universal predictor for processes in meandering rivers (Fathi et al. 2012).
Therefore, the main objective of this research was to obtain suitable roughness coefficient and appropriate turbulent model for Beshar River as a coarse-bed river which is located in Iran. Following this, the study of flow pattern in straight and meander reaches and also determination of different hydraulic parameters such as velocity, shear stress and Froude number in the studied reach of Beshar River are other important goals of this research.

CCHE2D model
CCHE2D is a 2D hydrodynamic and sediment transport model. This model is able to simulate the unsteady and steady open channel flows over mobile bed, and also it can be used in natural river for hydraulic engineers. In the CCHE2D, the efficient element method (a special finite element method) is applied to descretize the governing equations. Both supercritical and subcritical flow states can be simulated (Zhang 2006a). The CCHE2D has two software: CCHE2D Mesh Generator and CCHE2D GUI (Zhang 2006b). The Mesh Generator allows the user to create the structured mesh of the geometric characteristics and existing structures. The method applied in the structured mesh generation is grouped into two categories: algebraic method and numerical methods. The depth-integrated 2D formulas are solved in CCHE2D model (Zhang 2009).
Continuity equation: Momentum equations: where u and v are the velocity components in the x and y directions, respectively; z is the water surface level; g is the gravity acceleration; ρ is density of water; h is the water depth; f cor is the Coriolis factor; τ xy , τ xx , τ yx and τ yy are the Reynolds stresses; and τ bx and τ by are the bed shear stresses.

Simulation steps in CCHE2D
In general, the steps for setting up and running the simulation by the CCHE2D model are as follows: • Production of mesh. • Defining initial and boundary conditions.

Site characteristics and data
In the present study, CCHE2D model is used to simulate 12-km reach of the Beshar River in Iran. The study reach is located around the city of Yasuj, the capital of Kohgilouye and Boyerahmad province. The location of the study reach is shown in Fig. 1. The bed of the Beshar River, according to Fig. 2, is mainly coarse, and average diameter of the particles (d 50 ) is estimated to be 22 mm. Shahmokhtar hydrometric station (Shahmoktar station is located at 38 m upstream of the study reach outlet) is located at the downstream of this reach, which is monitored by the Yasuj regional water company, Ministry of Energy.
The topographic database of study reach, including longitude (UTM; X), latitude (UTM; Y) and elevation (Z), was taken by surveying operations in 2011-2012. Hydraulic parameters including the water surface level and maximum water depth at the cross section of Shahmokhtar station were collected in the low, medium and high discharges according to Table 1.

Production of mesh
Mesh generation is the first step in modeling, and it is important to create a suitable network to increase the accuracy of final results. The first step in the mesh production is loading the topographic database (*. Mesh_xyz). Then, the study reach is defined by the block boundaries, and algebraic mesh is produced using the measurement data points by CCHE2D mesh generator with 24 nodes across and 500 nodes along the river.

Defining initial and boundary conditions
To solve each partial differential equation, it is important to determine the initial and boundary conditions for the functional domain. The initial conditions are containing initial bed elevation and initial water surface. Initial bed elevation is defined by the mesh file and topographic data in CCHE-MESH generator and doesn't need to change in CCHE2D-GUI model, but the initial water surface level should be defined for the model. The initial water surface level according to the measured depth of water at the Shahmokhtar hydrometric station, for discharges of 13.32, 49.22 and 335.38 m 3 /s, for inlet, outlet and whole of the study reach according to Table 2 was done.
In other words, the initial conditions (water surface level) were applied to the inlet and outlet sections of the computational mesh and then were generalized by linear Fig. 2 Particle-size distribution curve of the study reach (Yasuj regional water company, Ministry of Energy)  interpolation to the entire study reach. So to start modeling, the initial water surface level at each point of the study reach is known. In addition to the initial conditions, the boundary conditions should also be applied at the inlet and outlet sections. So at the inlet boundary, the flow discharge can be defined as a constant value or as a discharge hydrograph. There are both options in the model, but the first one was considered. Steady state discharges of 13.32, 49.22 and 335.38 m 3 /s were used in different cases of simulation as the inlet boundary conditions for simulation. At the outlet section, the water surface level should be defined for the model as a boundary condition, which was 1712.46, 1713.15 and 1716.01 m, respectively, for discharges of 13.32, 49.22 and 335.38 m 3 /s.

Parameters setting
Numerical simulation is reproduced true physics by solving mathematic equations; therefore, many physical parameters are needed. Some physical parameters have been provided in the graphic user interface as default, which should be treated as guidance only. Many have to be provided by users for their particular applications (Zhang 2006b).
Users must also provide the parameters that control the simulation process. The roughness factor is an important parameter for the influence of the flow field. In the CCHE2D model, there are two methods for choosing the roughness coefficient: • The first method is to use a specific roughness parameter, so that the user can apply Manning's roughness coefficient (n) or roughness height (Ks). • The second method is to use the formulas of bed roughness.
Different turbulent models such as k-ε model, mixing length model and parabolic eddy viscosity model are also considered for CCHE2D. Selecting the appropriate turbulent model helps to solve the flow simulation and obtain more accurate results.
In this study, the Manning's roughness coefficient and turbulent models mentioned above are used, and finally the most suitable Manning's roughness coefficient and the best turbulent model are used to study the flow pattern.

Executing the model
After all the initial conditions, the boundary conditions and the model parameters are set, the simulation can be performed. The total simulation time and time step (calculated time scale) to achieve a stable state flow are set equal 432,000 s (5 days) and 2 s, respectively. The CCHE2D-GUI model was executed after all the setting and necessary information such as "Depth to consider dry" parameter of "simulation parameters" that it should be between 0.02 and 0.05 (Zhang 2006a). This was considered 0.03 m in this study.

Results and discussion
Changes in the flow and sediment regime of rivers often lead to modifications of their channel geometry (Borgohain et al. 2019). Therefore, verification of cross section profile at Shahmokhtar hydrometric Station is very important before starting flow modeling. Figure 3 shows the evaluation of the Shahmokhtar cross section profile produced by the CCHE2D-Mesh generator and measured. It is clear that both field and model provide the same results. As shown in this figure, the changes in the water profiles are consistent with the original structure and model. In some sections, there are minor differences that arise from the errors involved in data mining, and the software interpolation is defined according to the basic information, and most importantly human errors.
The water depth profiles (water surface level minus the river bed level) derived from CCHE2D numerical model and observational data according to eight roughness coefficients of 0.015, 0.020, 0.025, 0.030, 0.035, 0.040, 0.045 and 0.050 are presented in Fig. 4. The water depth profiles Fig. 9 Cross-profiles of water surface in the straight and sudden contraction reaches of the numerical model and observational data have the same trend for different roughness coefficients. To select the best roughness coefficient, the water depth profile was evaluated at the Shahmokhtar hydrometric station. For this purpose, statistical analysis was carried out between the observed and simulated data using two factors of "coefficient of determination or R 2 " and "linear correlation coefficient or r," as defined in Eqs. (4) and (5), for all tests according to Fig. 5.
where B is the simulated variable, A is the observed variable, Ā is the sample mean of A value, and N the total number of variables. Figure 5 shows that Manning's roughness coefficients in the range of 0.015 ≥ n ≥ 0.025 for both medium and high (5) Moreover, for a low discharge of 13.35 m 3 /s, the range of 0.035 ≤ n ≤ 0.050 offers more accurate results. It can be said that for low discharge, the effect of bed roughness on the flow is greater due to its lower depth, and therefore higher roughness coefficients for modeling have better and closer Fig. 11 Cross-profiles of water surface in the meander reach results to observational data. But, for medium and high discharges, the bed roughness has a lower impact on the upper layers (flow layers farther from the river bed) of the flow due to its greater depth. Therefore, for medium and high discharges, lower roughness coefficients are more appropriate. Using a roughness coefficient of 0.015 and a discharge of 49.22 m 3 /s, different turbulent models of k-ε model, mixing length model and parabolic eddy viscosity model were evaluated. For this purpose, the water depth profiles obtained from CCHE2D model for different turbulence models of k-ε, mixing length and parabolic eddy viscosity were compared with the observation data at Shahmokhtar hydrometric station, as shown in Fig. 6, and analyzed using statistical parameters of RSQ and r. According to Fig. 7, the highest values of RSQ and r (0.88 and 0.94 respectively) are for the k-ε model. So, this turbulent model for modeling is better than mixing length and parabolic eddy viscosity models and offers better results.
After calibration and executing the models, various hydraulic parameters such as velocity, shear stress and Froude number can be studied in different sections. In this study, various parameters for low, medium and high Along their direction, rivers usually have variant forms such as straight, braided and meandering, and the hydraulic flow in each one is different from the other. Flow around the bends and meandering reaches of rivers is one of the most important subjects for engineers, and their study is significant in order to know the hydraulic flow in these areas. However, the study of river bends is not a topic that has only existed in previous researches (De Vriend 1981;Odgaard 1984), but today, new methods such as numerical modeling (Ulke et al. 2017) can be used to analyze the flow in the bends and meander reaches of rivers.
As shown in Fig. 8, in Beshar River a straight reach from the cross sections 1-3, which leads to a sudden contraction, and a meander reach from sections 5 to 10 are considered.
The cross-profiles of the water surface and bed profiles obtained from the CCHE2D numerical model in different cross sections of the straight reach and sudden contraction from study reach are shown in Fig. 9. Due to the uniform distribution of flow in cross sections, and the uniform distribution of secondary currents from the flow center toward the left and right bank of the cross sections of the straight reach of the river, the latitude gradient of the water surface is approximately zero and the water surface is horizontal.
Moreover, the study of the velocity distribution profiles in different cross sections of the straight reach in Fig. 10 shows that, due to the greater resistance of the river bed on the flow at the river sides (left and right banks of the river), the flow velocity on the left and right banks is very low and close to zero. On the other hand, the maximum velocities are located in the middle of the cross sections, and this trend is more evident for higher discharges. Figures 11 and 12 show the cross-profiles of the water surface and the depth average velocity in the meander reach, respectively. In the bend 1, the latitude gradient of the water surface due to the presence of secondary currents increases to the outer bank of the bend (the left bank of the river). Also the largest velocities develop toward the outer bank, where erosion is stronger than the inner bank. The secondary currents are the result of energy difference between the water surface and the bed, and also the reduction in centrifugal force at the bottom. At the point of attachment bends 1 and 2, the direction of latitude slope of the water surface changes between cross sections 7 and 8. Then, in the bend 2, the gradient of water surface and the largest velocities due to secondary currents tend to the outer bank (the right bank of the river).
The results show the water surface level increases in upstream of the contraction region for discharge of 335.38 m 3 /s. This is because the flow in the upstream is subcritical (F r < 1) and then in the sudden contraction region becomes supercritical (F r > 1). The blockage in the contraction area increases the water surface level at the upstream. For medium and high discharges, due to the change in the flow from the supercritical to subcritical, a hydraulic jump occurs in this reach. But this phenomenon is not seen for low discharge. Table 3 shows the latitude gradient of the water surface and Froude number in different sections for both straight and meander reaches.
The critical bed shear stress (τ b,cr ) according to Shields is a function of the particle diameter (d) for a temperature of T = 10°, 20° and 30 °C, (ρ s = 2650 kg/m 3 , ρ = 1000 kg/ m 3 ). The influence of temperature is only significant for a particle diameter larger than 700 µm (7 mm); the τ b,crvalue is independent of temperature (viscosity), and in this case, the critical Shields parameter (θ cr ), defined as Eq. (6), is constant and equal to 0.055 according to the   Bonnefille (1963) and Yalin (1972) definitions (Van Rijn 1993. In this study, the average diameter of the bed material, which is mainly quartz origin (ρ s = 2650 kg/m 3 ), is 22 mm and larger than 7 mm. Therefore, the critical Shields parameter (θ cr ) is 0.055, and the critical bed shear stress (τ b,cr ) is equal to: Figure 13 shows the distribution shear stress in the studied reach from the Beshar River. As it is clear, the shear stress in some places, like around the Beshar bridge, is very high. Therefore, in these locations, necessary measures should be taken to prevent the hazard for these important and communications structures. On the other hand, Fig. 13 shows the flow zoning in all three-low, medium and high-discharges, the brown areas in the figure represent the places where there is no water and they are dry. Table 4 shows the average of some of the flow parameters in the studied reach for the three mentioned discharges. The RSQ and r factors between velocity and shear stress, according to Table 4, show that there is linear and direct relationship between these two hydraulic parameters.

Conclusion
The results show that, for low discharge due to the greater effect of the bed roughness on the flow, the Manning's roughness coefficient ranges of 0.035-0.050 will have better results, and for medium and high discharges, the roughness coefficient ranges of 0.015-0.025 will provide results that are more accurate. Also, using the k-ε turbulent model for this river offers better results.
The analysis suggests that the straight and meander reaches are different in structure and physics of the flow. So the latitude gradient of the water surface and the maximum velocities in the river arch develop toward the outer bank, and this is one of the reasons for erosion in the outer bank and deposition in the inner bank of the arch.
The simulation results also indicate that, on average, the depth average velocity (V ave ) for the discharges of 13.32, 49.22 and 335.38 m 3 /s in the study reach is 1.23, 1.33 and 2.73 m/s, respectively. Moreover, the average shear stress for the threelow, medium and high-discharges mentioned above is 6.69, 7.58 and 21.45 N/m 2 , respectively. Considering that the critical shear stress of the study reach is 19.58 N/m 2 , therefore, in areas of the study reach, where shear stress is much higher than this value, there is a need for protective measures.
The CCHE2D model has the ability to simulate longitudinal and cross-profiles of the water surface, physical phenomena such as hydraulic jump and local contraction, and generally the pattern of flow in the Beshar River. Therefore, in the engineering and operational projects of the Beshar River, the results of this numerical model can be trusted and used.