Experimental analysis and numerical simulation of bed elevation change in mountain rivers

Studies of sediment transport problems in mountainous rivers with steep slopes are difficult due to rapid variations in flow regimes, abrupt changes in topography, etc. Sediment transport in mountainous rivers with steep slopes is a complicated subject because bed materials in mountainous rivers are often heterogeneous and contain a wide range of bed material sizes, such as gravel, cobbles, boulders, etc. This paper presents a numerical model that was developed to simulate the river morphology in mountainous rivers where the maximum bed material size is in the range of cobbles. The governing equations were discretized using a finite difference method. In addition, an empirical bed load formula was established to calculate the bed load transport rate. The flow and sediment transport modules were constructed in a decoupled manner. The developed model was tested to simulate the river morphology in an artificial channel and in the Asungjun River section of the mountainous Yangyang Namdae River (South Korea). The simulation results exhibited good agreement with field data.

new 1D numerical model to calculate flow and sediment transport in steep mountain rivers. According to Mosconi (1988), failure to predict bed loads in mountain river flows may be due to most common equations not considering the morphological peculiarities of study areas in conjunction with their limited capability to cover a wide distribution range of bed material sizes. Bed load equations obtained by several authors on low slopes are rarely applicable to mountain rivers, where the river beds contain wide ranges of bed material sizes, large roughness elements, etc. Therefore, these features largely affect the research results (D' Agostino and Lenzi 1999).
In this paper, a 2D numerical model has been developed to simulate the flow and morphological changes in steep channels where bed materials have large size distributions. The model system consists of a flow module and bed load transport module. The flow module is based on the mass and momentum conservation equations in the Cartesian coordinate system. The sediment transport module only comprises empirical bed load formulas. The river morphology module is based on the sediment continuity equation, and a grain material distribution is applied for individual size fractions. The solution method was implemented in a computer source code and written in structured Fortran 90.

Governing equations
By assuming a hydrostatic pressure distribution and neglecting wind shear and Coriolis acceleration, the depth-averaged 2D governing equations are expressed in Cartesian coordinates in the following forms (Ahmadi et al. 2009;Horritt 2004;Lai 2010;Dang and Park 2015).

Sediment transport equations
In the past, most bed load formulae have been developed and widely used based on laboratory investigations with uniform particle size. Unfortunately, when applied to natural rivers with non-uniform particle size, the calculated bed load transport rates often differ by orders of magnitude and do not exhibit high confidence levels. One of the major causes is often the influence of the local conditions, which are very different from the laboratory conditions where the bed load formulae were constructed. As noted above, mountain rivers is a dominant and fundamental process in the hydrodynamic rivers, and often have non-uniform particle size distributions, the particle sizes are divided into several fractions. Therefore, in this study, only bed load is used and no suspended sediment is included. An empirical formula has been developed by the Department of Civil Engineering, Gangneung-Wonju National University (Park et al. 2013), based on investigation data from mountain rivers in South Korea. The bed load transport formula can be applied to mountain rivers in which the channel bed materials are non-uniform and includes, gravel, and cobble. The bed load transport formula is given as follows: where q * sb is the bed load discharge, q * sb is the bed load transport capacity, σ s is the specific gravity, g is the gravitational acceleration, d is the particle size of the bed material, τ * ci is the dimensionless critical Shields stress of incipient motion, and τ b is the bed shear stress. τ b , σ b , and σ s are expressed as follows: where H is the water depth, d is the particle size of the bed material, S is the stream slope, γ s is the specific weight of sediment, and γ is the specific weight of water (Fig. 1).

Bed level variation
By considering only the bed load transport, the 2D sediment continuity equation may be written as follows: where Z b is the local bed elevation, p is the porosity of the bed material, ∂Z b is the change in the local bed level during the time interval ∂t, and q bx and q by are the x-and y-components of the bed load transport per unit width, respectively. The components of the bed load transport in the x-and y-directions are related to the bed load q b as follows: where q sb appears in Eqs. (11) and (12) and is calculated by Eq. (5). In Eqs. (11) and (12), α is calculated by Eq. (13) and corresponds to the directional angle of bed load transport in the x-and y-planes, as follows: where f s is the sediment shape factor.
Several studies have been carried out to propose a formulation for f s , such as Ikeda (1982), Kovacs andParker (1994), andZimmerman andKennedy (1978). Talmon (1992) used the following expression of the sediment shape factor: where d 50 /h is the relative roughness parameter. Note that this formulation controls the effect of gravity on the sediment particles. In Eq. (13), the term β is calculated as follows.
The terms A and r s in Eq. (15) are given as follows.
(11) q bx = q sb cos α (12) q by = q sb sin α Relationship between τ * c and σ (d i /d m ) (Park et al. 2013) Rozovskii (1957) and Engelund (1974) suggested that the value of A was 11 and the value of r s was seven. The transverse bed slope was small and had little effect on the flow calculations. Therefore, in the research, the effect of transverse bed slope is ignored.

Discretization of governing equations
The governing equations are discretized in the computational domain using an finite difference method (FDM) and a staggered grid in Cartesian coordinates. In this explicit difference formulation, a first-order approximation was used for the temporal derivative (∆t). Second-order central difference approximations were used for space discretization (∆x, ∆y), where the water depth (h) is defined at the primary grid center (i, j) and the velocity components (u, v) are defined at the cell faces of the secondary grid (i + 1/2, j + 1/2) ( Fig. 2) (Hung et al. 2009;Jia and Wang 2001;Dang and Park 2015).
Equation (1) is applied at grid point (i, j), yielding the following finite difference equation.
The x-momentum equation is applied to the secondary grid at grid point (i + 1/2, j).
Finally, the y-component is applied at grid point (i, j + 1/2), yielding the following finite difference equation.

Fig. 2 Staggered grid with boundary cells
For this purpose, the following approximations are used.
The sediment continuity equation given by (10) is discretized as follows.
In Eq. (18), the discharge components are specified at time step (T + 1). The momentum equations [Eqs. (19) and (20)] in the x-and y-directions are solved first to provide P I+1/2,J and Q I,J+1/2 , (with P = u·h and Q = v·h) the values of the discharge components at time step (T + 1). The water depth at the time step (T + 1) in Eq. (18) is then solved. Equation (18) is solved iteratively to determine the value of H i,j at time step (T + 1) over the entire domain. Subsequently, Eqs. (19) and (20) are used to compute the velocity components at time step (T + 1) over the entire domain (Horritt 2004;Hung et al. 2009;Matyka 2004;Dang and Park 2015). The steps in the process of solving the flow and bed load equations are listed below and are illustrated in Fig. 3.
Step 1 Initializing all the variables. This step usually corresponds to time T 0 . In this step, the values of the water depth and flow field within the computational domain and at the boundaries are specifically established. It is assumed that the velocity components and water depth are known at time T 0 and that the boundary conditions of the velocity components and water depth are given.
Step 2 Partial differential Eqs. (18), (19), and (20) for the flow and Eq. (21) for sediment continuity are solved with the finite difference code. Discretized equations are obtained for the shallow water and sediment continuity equations using the staggered numerical grid. The initial and boundary conditions used to solve the momentum Eqs. (19) and (20) and the continuity Eq. (18), i.e., the values of u, v, and h at time T + 1, are determined at every interior node (I = 2,…, N). The values of the dependent variables u and v at the boundary nodes 1 and N + 1 are determined using the boundary conditions. The values of the dependent variables that are not specified through boundary conditions can be determined by extrapolation of the interior points or equivalently by approximation of the derivatives at fictitious boundary points. We then obtain the corrected water depth and (u, v) velocity components at every interior node in the computational domain.
Step 3 The velocity components are calculated at time step T = T 0 + ΔT until a converged solution is obtained. In this step, convergence criteria must be checked because the scheme used in this research is an iterative scheme. Then, the velocity components and water depth are updated with their corresponding values.
Step 4 The water depth and velocity components are used to calculate the dimensionless particle diameter, dimensionless Shields stress, dimensionless critical Shields stress, critical shear stress, boundary shear stress, etc. Finally, the dimensionless particle diameter is calculated.
Step 5 The parameters calculated in Step 4 are used to calculate the bed load transport rate.
Step 6 Erosion and deposition are calculated using the sediment continuity Eq. (21) to determine the bed level variation and update the new water depth if the channel bed has changed. Step 7 Return to Step 2 and repeat the preceding calculation until the specified final time. If a steady state solution is required, a specified convergence criterion must be satisfied.
Step 8 The last step in the calculation process involves storing and updating variables at each time step, moving to the next time step, and repeating Step 2 through Step 7.

Stability conditions
In explicit difference schemes such as the MacCormack, Lax-Wendroff, and Marker and Cell schemes, the magnitude of the time step is governed by the CFL stability condition (Bellos and Hrissanthou 2003;Chow and Ben-Zvi 1973;McKee et al. 2004McKee et al. , 2008Paulo et al. 2007;Rao 2003). In this study, the following expression for the CFL stability condition was used: where ∆t is the time increment; ∆x and ∆y are the grid spacings; u and v are the velocity components in the x-and y-directions, respectively; h is the water depth; g is the acceleration of gravity; and α is the coefficient (α ≤ 1).

Application of model verification and discussion
To investigate the applicability of the developed model, the present model has been tested in two experimental cases. The first case was obtained from the Large Scale Hydraulic Models of the University of Calabria, Italy (Bellos and Hrissanthou 2003;Bor 2008;Miglio et al. 2009). The second was obtained from a flood event in the Asungjun River. Numerical models can be calibrated by comparing measured and computed results and adjusting the empirical coefficients in the associated empirical relationships. By a trial and error procedure, the agreement between calculations and measurements can be satisfied. However, this procedure is difficult to apply because of the lack of input data, especially for simulating flow and sediment transport in natural rivers. Several researchers have determined the goodness of fit of hydrodynamic models by computing the root mean square differences (RMSD) and mean absolute errors (MAE) between observed and simulated results.
In this study, the expression for calculating the RMSD (see Eq. 23) was selected to determine the error between the calculated results and the measured data, minimizing the RMSD error between the calculated results and measured data. The RMSD error is defined as follows.
The root mean square deviation provides a measure of variance between the observed and simulated results.
The expression used to calculate MAE is defined as follows: where U M is the measured value, U C is the calculated value, and N is the total number of samples.

Model setup
First, the developed model was tested by simulating the artificial channel. An experimental facility was installed at the Large-Scale Hydraulic Model at the University of Calabria, Italy (Fig. 4). The experiment was conducted and established with the following conditions: • The experimental flume data were measured over a period of 30 min with a time step of 1.0 s.

Results and discussion of the seal aggradation test case
The comparison of the bed level variation showed small differences between the measured and predicted results. The simulation result was in agreement with the measured data during the simulation period (Fig. 5). In general, a good correspondence was Fig. 4 Illustration of the experimental flume at the University of Calabria (Lai 2010) observed between the aggradation and degradation tendencies along the experimental flume, with an RMSD of 0.53. Similarly, the validation of the calculated and measured bed level variations using the Nash-Sutcliffe criterion (Nash and Sutcliffe 1970) corresponds to a value of 0.98. This confirms that the simulation model of bed level variation is quite accurate.

Model setup
Next, the developed model was tested by simulating flood events in the Asungjun River section. Bed loads in mountainous rivers play important roles in the evolution of river beds. During the flood season in the studied river section, the bed material is predominantly sand, gravel, and cobble. The main mode of sediment transport when the unit discharge is more than 2.5 m 3 /s is bed load transport (Fig. 6). Hydrographs of the water level at the upstream boundary (Fig. 7) were provided because the flow in steep mountain rivers is often critical; however, no hydrographs were provided for the water level at the downstream boundary.
Bed topography was collected from the Sokkia-C32 measuring device. All topographic and bathymetric data measured before and after the flood event (Fig. 8) are presented in the form of x-, y-, and z-points, corresponding to the longitude, latitude, and water depth of the computational domain ( Fig. 9). They were then imported into an Excel file and interpolated into mesh points. They are based on the original data and interpolated mesh elevations. The surveys were then merged with archived data to form a single point data file in xyz form. Detailed information on the bed topography at an initial state is given in Fig. 9.
Simulation data were collected from January to November 2012 (Fig. 10). Water level stations were established at the inflow and outflow boundaries, and points inside of the study area were established to collect water level data using pressure sensor gauges. The location of each station was selected such that the study area was sufficiently covered to capture the water surface fluctuations with a high degree of accuracy.
Field surveys were used to collect and analyze bed material sizes during a flood event (Fig. 11). The measured time series of sediment discharge at the inflow and outflow boundaries were established as boundary conditions.  According to Bravo-Espinosa et al. (2003), bed load transport conditions in alluvial channels are dependent on individual particle size fractions rather than on the complete spectrum of particle sizes represented by one characteristic particle size. Based on this viewpoint, to increase the accuracy of the sediment transport module in the numerical model, we divided the mean bed material size into several fractions based on data measured in the study area (Table 1).
In Table 1, d im is the material diameter and d i is the mean material diameter. In this test case, the model was applied to simulate a flood event in the Asungjun River section, which has highly complex geometrical features. Therefore, the use of a single value for the Manning's roughness coefficient is not appropriate because the channel bottom topography consists of widely different features. Thus, the computational domain in this study is divided into several different roughness zones, and Manning's roughness coefficient is divided into several roughness zones. In Fig. 12, Zone 1 represents sand bars, Zones 2 represents gravel, Zone 3 represents boulders with heavy vegetation, and Zone 4 consists of boulders. Each zone was assigned a Manning's coefficient that was determined through a calibration study by comparing the measured and predicted water levels. The calibrated Manning's coefficients are listed in Table 2.

Wet/dry treatment
In shallow water regions of natural rivers where the water depth is small and the channel bed exhibits irregular geometry, the water edges change with time. In those cases, Fig. 11 Illustration of a field survey to collect the bed load and suspended sediment wet and dry treatments in numerical simulations are often used to determine the wet and dry cells. The water depth defined by the user will often depend on the scale of the simulation. In numerical models, the process of drying and wetting is represented by a flow domain that becomes dry when the water depth decreases and wet when the water depth increases. Existing 2D models have taken a number of approaches to solve the problem associated with some areas being wet and others dry, or fluctuations between the two (Bates and Hervouet 1999;Begnudelli and Sanders 2006;DHI 2003). Several models turn cells on and off based on the minimum depth criteria (Delft 2002;King and Roig 1988;Leclerc et al. 1990). Other models change the fluid properties at very small depths so that a very thin layer of fluid is always present. Most approaches attempt to reformulate the flow equations over partially wet elements by introducing a scaling coefficient, representing the true volume of water at each element. This coefficient varies from zero to one as the cells tend from fully dry to fully wet (Bates and Hervouet 1999;Defina 2000).
In this study, a threshold value of the water depth (0.03 m) based on the river bed material size is used to establish drying and wetting. If the water depth in a cell is larger  Table 2 Manning's coefficients in the roughness zones shown in Fig. 13 Zone Zone 1 Zone 2 Zone 3 Zone 4 Manning's n 0.033 0.042 0.052 0.048 than this threshold value, this cell is considered wet, and if the water depth is lower than this threshold value, the cell is dry.

Results and discussion of the Asungjun River section case
The developed model was compared with the Mike1C model and calibrated using measured data. The simulation results of the developed model were slightly lower than those of the Mike21C model and survey data (Fig. 13). Additionally, the model slightly under predicted stages compared to stages predicted by Mike21C model and survey data. The simulation results of water level from the Mike21C model showed that the RMSD was 0.0039 and the ADM was 0.043, while the corresponding values of the developed model were 0.0032 and 0.030 (Table 3). Similarly, the Nash-Sutcliffe values between the measured water level and the Mike21C model and developed model were 0.89 and 0.97, respectively. This result suggests that the simulated model of the water level is very reliable. Figure 14 shows a comparison of the developed model, Mike21C model, and measured bed elevation along several cross sections of the studied river section. The maximum aggradation at cross section No.15 was determined to be approximately 0.60 m compared with the original bed. This value was 0.22 m using the developed model and 0.32 m using the Mike21C model. Similarly, the maximum aggradation at cross section No.37 was calculated as 0.46 m compared to the original channel bed. At cross section No. 37, the value of the Mike21C model was 0.91 m and the value of the developed model was 0.34 m. The Nash-Sutcliffe values corresponding to the results of the Mike21C model and developed model were 0.76 and 0.97, respectively. This result suggests that the developed model is very reliable (Fig. 15). Generally, the simulation results of river morphology at different cross sections showed that the developed model was  in good agreement with the observed data compared to the agreement of the Mike21C model. Figure 16 shows the final bed level configuration after the flood event along the channel, and the differences between the simulation results of the developed model, the Mike21C model, and measured data are compared. The simulation results of the river morphology from Mike21C model exhibited an RMSD of 0.0028 and ADM of 0.035, while the corresponding values of the developed model were 0.0037 and 0.024. Similarly, the Nash-Sutcliffe criterion was used to validate the Mike21C model, developed model, and measured data. The Nash-Sutcliffe values corresponding to the Mike21C model and developed model were 0.81 and 0.96, respectively. These values confirm that the developed river morphology module is quite accurate.

Conclusions
A depth-averaged 2D numerical model was developed for simulating river morphology in mountain rivers. FDM is used to solve the momentum equations and sediment continuity equations. The model system consists of flow and river morphology modules. Both the flow and river morphology modules are solved using an iteration method that constitutes a coupling procedure. The bed material size distribution is separated using a fractional approach. This approach is more complex than the classical method, which only uses a value of particle size diameter (D 50 ).
The simulation results of the river morphology of the flood event in the natural river using the developed model are more accurate than those produced by the Mike21C model. Generally, the simulation results were in good agreement with the measured data compared to the results of the Mike21C model.
The advantages of the developed model used for the simulation of the experimental channel and flood event are as follows: • The robustness of the developed model under the various cases studied, such as lateral water and abrupt cross section variations, division of bed material into a number of size fractions, and division of Manning's roughness coefficient into different values in study zones to fit the real bed topography conditions. • The simple structure of the developed model allows users to easily control the calculation procedure, and it has a relatively fast computational speed.
The disadvantages of the developed model are as follows: • The complicated procedure of constructing the numerical grid; • The sensitivity of the coefficients to the river morphology; • The need to calibrate multiple parameters when constructing the bed load formula.
More testing of the model may be necessary to improve its predictive ability. It is expected that the model will become a useful predictive tool for mountainous river studies.

Authors' contributions
In this study, SDP proposed the research project and outlined the project, designed the research proposal, and wrote the manuscript. He also has participated in field surveys, analyzed output data, and prepared and edited the manuscript. TAD participated in the field surveys, analyzed measurement data, analyzed and established input-output data, wrote source code, implemented the simulation model, analyzed the simulated results, edited the manuscript, and wrote the main paper. Both authors contributed equally to this work. TAD jointly conceived the study with SDP. We participated in field surveys, collected and analyzed input-output data, and administered the experiments. Both authors discussed the results and commented on the manuscript in all stages. Both authors read and approved the final manuscript.