Application of the Mike21C model to simulate flow in the lower Mekong river basin

Numerical models are useful tools that play an important role in many research projects. Mike21C is one of the most well-established models for simulating a variety of processes, including bank erosion, bed level variations, aggradation, and degradation. Such processes are caused by a variety of activities, such as construction and dredging, as well as seasonal flow fluctuations. Mike21C based on an orthogonal curvilinear grid, which enables a computational speed that is faster than that of other grids. The hydrodynamic part of the model is based on solving the Saint-Venant equations. In this research, the Mike21C model was applied to simulate water depth, flow discharge distribution, and suspended transport rate along reaches of the Bassac River the on Vietnam. The data files of the curvilinear grid and bathymetry were used to generate the model. A time series for the discharge and water level of hydrological stations was established. The model was calibrated using the water level and suspended load data collected during high and low flow discharges. The simulation results show that the model can be applied well to other areas.


Mike21C model description
In this paper, the Mike21C model is used to simulate the water depth, flow discharge, suspended load and bed load concentration of the study river reaches in the flood and dry season. Mike21C is one of the most comprehensive and well-established tools for simulating river bed and channel planform development caused by changes in the hydraulic regime. Simulated processes include alluvial resistance, bank erosion, and scouring and shoaling caused by various activities, such as construction and dredging, and seasonal flow fluctuations (Khue 1985;Lai 1987;Jin and Steffler1993). This model is approximated by using FDM in curved coordinates (Ahmadi et al. 2009;Beck and Basson 2007;DHI 1995;Dang and Park 2016;Talmon 1992;Gulkac 2005;McGuirk and Rodi 1978). Structurally, Mike21C has three main modules: the flow module, sediment transport module, and river morphology module. Mike21C model is applied to simulate the water level fluctuation, flow discharge distribution, suspended load transport rate, and bed level variation in the river downstream.

Hydrodynamic module
The flow module based on the three-dimensional (3D) hydrodynamic model is complex. Application of the 3D model for simulating long time scales (e.g., months, season, and years) elevation to river morphology is a complicated process. To overcome this obstacle, scientists have converted the main hydrodynamics module into 2D equations representing the conservation of momentum and mass horizontally (DHI 1995;Ye and McCorquodale 1997).
The hydrodynamic equations are expressed as follows: where s, n are the co-ordinates in the curvilinear co-ordinate system; h is the water depth; p and q are the mass fluxes in the s-and n directions, respectively. H is the water level; C is the Chezy roughness coefficient; g is the gravitational acceleration; R s and R n are the radius of curvatures of s-and n-line, respectively; RHS is the right hand side describing the Reynolds' stresses.

Sediment transport module
In the sediment transport module, the suspended load transport equations under the control of convection and diffusion are expressed as follows (Duc 2004;Meyer and Müller 1948;Galappatti 1983;Jia and Wang 2001): (1) where z is the vertical coordinate; W S is the settling velocity of the sediment particles, c is the suspended load concentration; ε is the eddy viscosity coefficient; u, v, and w are the flow velocity components in the x-, y-, and z-directions, respectively. Ignoring the limited diffusion outside of the vertical diffusion, (4) becomes: Bed load (S bl ) is very closely related to the suspended load. Many formulas of the bed load transport are based on the calibration coefficients k b and k s . Engelund and Hansen (1967) had established the relationship k b + k s = 1, which is used in many models, including Mike21C (DHI 1995).
where k s and k b are the suspended load and the bed load coefficients, respectively S tl is the total volume of sediment transported determined according to the formula: where ρ ds is the relative proportion of sediment (relative density of the sediment); d 50 is the diameter of the sediment particles; Shields parameter θ is determined by; where τ is the flow shear stress; ρ is the water density (density of water); Shear flow is divided into two types: form drag and friction, as estimated based on local flow velocity u and the local Chézy C.

Morphological module
In the river morphology module, the hydrodynamic solution must first be obtained before solving the sediment transport equation. Next, the river bed and hydrodynamic model are applied (Vriend and Struiksma 1983;Koch 1980;Mosselman 1992;Odgaard 1983;Olesen 1987;DHI 1995). The equation describing bank erosion is given as follows: where E b is the bank erosion rate in m/s; S is the near bank sediment transport; z is the local bed level; α, β, and γ are the calibration coefficients specified in the model. Calculation of the bed level variation is based on the sediment continuity equation in the Cartesian coordinate system (Vanoni 1975;Galappatti 1983;DHI 1995): where S x , S y are the total volume of sediment transported along x and y, respectively; p is the porosity of the bed; ∆S e is the excess sediment supply from erosion of the bed (Fig. 1).

Study area
In this study, the Mike21C model is applied to the Hau River section belonging to the lower Mekong River Delta. The river reaches under study has a length of 20 km ( Fig. 2) from the area adjacent to the Chau Thanh District downstream to the area adjacent to Can Tho City. The left side is Cho Moi District, and the right side is Long Xuyen City in An Giang province. The flood discharge in the main river at the right bank can reach up to 13,000 m 3 /s.

Initial conditions
The important steps in the procedure to set up and solve the Mike21C model are described as follows. The first step involves the creation of a suitable curvilinear grid. A curvilinear grid is created by establishing K = 50 m, J = 100 m, and Jn × Km = 157 × 32 cells within the computational domain (Fig. 3). Next, a bathymetry data file with coordinates and river bed level were obtained from Echo-sounding in Jan 2014. Topographic and bathymetric data measured and presented in the form of x, y, and z points, corresponding to the longitude, latitude, and water depth of the computational domain. Next, these data are imported into an excel file and interpolated into the mesh points. These mesh points are based on the original data and the interpolated mesh elevations. Detailed information on the bed topography at the initial state is given in Fig. 4. By analyzing the mean water level for many years at the study area we found that if we select respectively the water level values of 40 cm and 160 cm and flow discharge is 3000 m 3 /s and 6500 m 3 /s in the dry season and flood season. The simulation run time would significantly reduce, because these water level values are very close to the real values of the water level. In the numerical simulation of the open channel with irregular geometry, the water edges change with time are considered with part of the nodes being possibly wet or dry. In shallow water regions, where the water depth has a small value, the momentum terms are often ignored. The Mike21C model requires the smallest water depth to define a wet or dry cell. This dry value is set along the boundaries (Dac 1986;DHI 1995;Dang and Park 2016;Wu 2004). Manning's coefficient and other parameters are selected as model calibration parameters along the river reach during simulation (Table 1).

Boundary conditions
The inflow and outflow boundaries that describe the hourly water level time series were obtained from the hydrology station of Long Xuyen (Fig. 2). The inflow boundary that describes the hourly discharge time series was obtained from the hydrology station of Long Xuyen using Acoustic Doppler Current Profiler (ADCP) measuring device. The calibration process was performed for the time period from 04 to 30 Apr 2014 (dry period) and 10 to 30 Sep 2014 (high flood period).

Flow velocity during dry season and during flood season at the study area
The calculated results of the flow velocity distribution during the driest season from the Mike21C model (Table 2) showed that the flow velocities and ebb velocities at the two ends of the river reaches are the same. The flow velocity at the head and end of the river reach is 1.27 and 1.26 m/s, respectively. The ebb velocity at the head and end of the river reach is 1.03 and 1.01 m/s, respectively. At the branch of My Hoa Hung and Pho Ba, the islet flow and ebb velocities are 1.0 m/s. At the right-hand side of the My Hoa Hung islet, the flow velocities are higher than 1.2 m/s. General simulation results are consistent with the real condition because the bed river topography in this river reach is narrow; thus, the velocities must increase. At the left-hand side of the Pho Ba, My Hoa Hung islet, Tien, Noi sand dune and at the right-hand side of islet Pho Ba and Tien sand dune, the velocities ebb and flow are smaller than 1.0 m/s due to the sediment deposition at the river bed.
According to the simulated results of the flow velocity at the flood season (Table 3) at the head of the left-hand islet of Pho Ba, My Hoa Hung, and Tien, the velocities are less than 2.0 m/s. There is only flow water and no ebb water in such seasons. At the head of the river reach, the head of the right-hand islet of My Hoa Hung, and the left-hand islet of My Hoa Hung, the velocities are over 2.0 m/s. Flow velocities greater than 2.0 m/s is the cause of frequent flooding or river bank landslides during flood season.

Flow discharge distribution
Flow discharge distribution into river branches at the driest season is Q = 6000 m 3 /s when the lowest tidal water is H = −60 cm. Table 4 shows the flow discharge in the river branches. The water discharge into the left-hand islet My Hoa Hung accounts for 84.5% of the total water discharge at the head of the river reach. At the Pho Ba islet, the water discharge is divided into 74.7% at the left-hand side of Pho Ba inlet and 9.8% at its right-hand side. Water discharge into the left-hand islet My Hoa Hung is 15.5% of the total discharge at the head of the river reach. The water discharge is divided into 8.5% at the right-hand side of the Tien sand bar and 4.6% at its left-hand side.
Flow discharge distribution into river branches at the flood season is Q = 13,000 m 3 /s and corresponds to the highest water level of H = 200 cm. Table 5 shows the flow discharge in the river branches. Water discharge into the left-hand side of islet My Hoa Hung accounts for 83.8% and the left-hand side of the islet accounts for 16.2% of total water discharge at the head of the river. As a result, the left-hand side of islet My Hoa Hung often suffers from more serious river bank landslides compared to other regions during flood season. At the Pho Ba islet, the water discharge is divided into 71.8% at its right-hand side and 12.0% on its left-hand side. At Tien sand bar, the water discharge is divided into 8.6% at its right-hand side and 5.2% on its left-hand side.
The simulation results of the flow discharge at Long Xuyen station 9 km away from the upstream boundary are presented in Fig. 5. In general, the simulation results of the flow Table 4 Flow discharge distribution in dry season at the study area  discharge are in good agreement with the measured data. The difference between two successive troughs of the measured and simulated flow discharge is very small, 1.57%, whereas the difference between two successive peaks is 6.75%. By applying the Nash-Sutcliffe criterion (Nash and Sutcliffe 1970), the validation of calculated and measured flow discharge data with Nash-Sutcliffe is 0.83. This result demonstrates that the simulated model of the flow discharge is of high accuracy.

Time (hour)
Measured data Calculated data Fig. 6 Comparison of the measured and predicted water level at Long Xuyen station (see Fig. 2) during the dry season with NASH = 0.80 Simulation results of the water level during the flood season ( Fig. 7) showed that the water level oscillation at the studied river reach is semi-diurnal. The difference between two successive troughs is less than 10 cm, whereas the difference between two successive peaks is 20 cm; ΔH Ltb = 23 cm, ΔH Lmax = 90 cm, ΔH Xtb = 27 cm, and ΔH Xmax = 76 cm. The average time of flow tide ΔT L and ebb tide ΔT X is 5 h 35 min and 8h21 min, respectively.
The validation of the calculated and measured water level data indicates that the Nash-Sutcliffe index ranges from 0.80 to 0.87. This result implies that the simulated model of the water level is very reliable.

Suspended load and bed load concentration distribution
The calculated results (Table 6) showed that the suspended load is low during in the dry season, approximately 0.01-0.03 kg/m 3 , and is rapidly increasing during the first heavy rains of the season, up to 0.45-0.75 kg/m 3 .

Time (hour)
Measured data Calculated data Fig. 7 Comparison of the measured and predicted water level at Long Xuyen station (see Fig. 2) during the flood season with NASH = 0.87