A two-dimensional numerical model predicting flow over a mobile bed has been developed. Governing equations consist of the shallow water equations and the Exner equation. The finite volume method on an unstructured triangular grid was deployed to discretize the governing equations. The local Riemann problem is solved by the Harten, Lax and van Leer–contact (HLLC) method in the interface of the cells and the equations are solved using a fully coupled method. Then the flux modelling has been deployed by the total variation diminishing (TVD) version of the weighted average flux (WAF) scheme. The model was verified by comparison of the results and available experimental data for dam-break flow, in a laboratory test, via a channel with sudden enlargement and erodible bed conditions. Comparison of these two sets of results shows that increasing the accuracy of flux modelling caused the model results to have a reasonable agreement with the experimental data.

## INTRODUCTION

Dam-break flows usually propagate along rivers and floodplains, where huge and sudden fluid flow occurs with sediment transport and bed changes (Xia *et al*. 2010); however, a flow model is usually applied for flow over a fixed bed. Nonetheless, remarkable models have been developed which are capable of simulating morphological evolution of the bed.

In recent years, applications of the schemes based on the finite volume method (FVM) in simulating dam-break phenomena have been increased. Especially, the upwind Godunov-type schemes via the FVM are very popular due to their robustness and accuracy in capturing discontinuities. A comparison of the performance of upwind schemes based on FVM, namely the Osher scheme, Roe scheme, Harten, Lax and van Leer (HLL) scheme and Harten, Lax and van Leer–contact (HLLC) scheme, was performed by Erduran *et al*. (2002). They investigated and compared the accuracy, applicability, efficiency and stability of different Riemann solvers in the finite volume shallow flow model. According to the mentioned study, it is recommended that the HLLC Riemann solver is suitable for all kinds of applications. Many researchers have developed and customized schemes based on FVM.

One of the first researchers is Patankar. Patankar & Spalding (1972) used FVM in order to find a numerical solution for the heat transfer equations. Then De Vriend (1987) worked on the theoretical basis of the behaviour of two-dimensional sediment transport and morphological models. Struiksma (1985) developed a two-dimensional sediment transport model to simulate large-scale bed change. Van-Rijn (1987) developed a two-dimensional vertical model for predicting morphological changes via laboratory testing. Hirsch (1988) employed the discretization of the integral form of conservative equations across a discontinuity. Shimizu & Itakura (1989) developed a two-dimensional bed load transport model for alluvial channels. Chaudhry (1993) used shock capturing methods for simulating one-dimensional open-channel flows. Kassem (1996) developed a two-dimensional bed load sediment transport model for straight and meandering channels. Anastasiou & Chan (1997) provided a model based on the shallow water equations using FVM on unstructured triangular meshes. Valiani *et al*. (2002) utilized the conservative form of shallow water equations with a Godunov-type scheme. Yoon & Kang (2004) used the HLL approximate Riemann solver and multidimensional slope-limiting approach in some test cases. Cao *et al*. (2004) presented a one-dimensional dam-break model considering erodible sediment beds. Loukili & Soulaimani (2007) developed a model using weighted average flux (WAF) approximation for modelling fluxes at the cell interfaces on an unstructured grid. Goutiere *et al*. (2008) investigated eigenvalues of a complete system (flow and sediment transport) and eigenvalues of the water phase only in a one-dimensional model. Iervolino *et al*. (2010) surveyed the effect of sudden enlargement of the channel on simulation of dam-break with a mobile bed. Benkhaldoun *et al*. (2012) and also Murillo & Garcia-Navarro (2010) used FVM to solve a coupled model of water flow and sediment transport and bed evolution. Soares-Frazao & Zech (2010) introduced a new approach consisting of estimating wave speeds considering sediment properties in two-dimensional shallow-water flow and erodible bed conditions. Li *et al*. (2013) developed a morphological model using a fully coupled, total variation diminishing upwind-biased centred scheme. Hsu *et al*. (2014) presented a model for surveying dam-break wave propagation and its implications for sediment erosion.

In this study, a two-dimensional flow model over a mobile bed has been developed. The HLLC scheme with novel wave speed proposed by Soares-Frazao & Zech (2010) is selected for flux modelling in the cell interfaces. To achieve second-order accuracy, the WAF method is employed, which was introduced by Toro (1999). According to Godunov's theorem, schemes with high-order accuracy generate spurious oscillations near large gradients of the solution. The WAF method is used with a flux amplifying function for avoiding the mentioned problem and getting a non-linear total variation diminishing (TVD) scheme of second-order accuracy.

In the following part, the governing shallow water equations are described, and then the details of the development of the numerical model are presented. The model is applied for dam-break flow in a channel with a sudden enlarged section and erodible bed conditions.

## METHODS

### Governing equations

*t*is the time,

*x*and

*y*the space coordinates,

*g*the gravitational acceleration,

*h*the water depth,

*u*and

*v*the depth-averaged velocity components, respectively, in the

*x*and

*y*directions, and the unit discharge components, the erodible bed elevation, and the components of the sediment transport rate per unit width and the bed porosity.

*S*, and are the bed slopes in the

*x*and

*y*directions, respectively, and the friction slope components. The friction slope is calculated using the Manning formula. The sediment transport rate is calculated using a suitable relation. The Meyer-Peter and Müller-type formula with ad hoc coefficients is used (Ribberink 1998) for sediment transport calculation because it appeared well suited and reasonable to the experimental test with an enlargement section in our study. The Ribberink formula is as follows: where

*s*is the relative density of sediment, is a representative grain diameter, denotes the non-dimensional bed shear stress and the non-dimensional critical bed shear stress. It is thought that sediment starts to move when . The non-dimensional bed shear stress is calculated as follows: where and

*n*is the Manning friction coefficient.

In schemes based on FVM, the fluxes are calculated across the cell interfaces; so, using the rotational invariance property in the direction normal to the considered interface, the variables are described in a local coordinate system related to the interface. The local coordinate system for an unstructured triangle mesh is shown in Figure 1.

*x, y*) to the system , the system of Equation (1) becomes (Soares-Frazao & Zech 2010): where As was mentioned before, according to the basic principle of the finite volume schemes, only the terms in the normal direction are considered. For analysis of the wave-propagation properties of the system, Equation (5) is transformed to the following equation: Equation (7) corresponds to the so-called augmented one-dimensional problem (Toro 1999) that is identical to a split multidimensional system in Cartesian coordinates. The eigenvalue analysis was carried out completely from Equation (7) and the validation process of the approximate analytical expression for the eigenvalues was performed accurately (Soares-Frazao & Zech 2010). Some corrections are implemented because of using Ribberink's formula (Ribberink 1998) instead of the Meyer-Peter and Müller original formula. Then the eigenvalues resulting from the mentioned process are in two types as follows (Soares-Frazao & Zech 2010): where is the sediment transport rate in the normal direction and calculated by the following relations: As was mentioned before, the validity of two types of eigenvalues has been checked by comparison with the exact eigenvalues (Soares-Frazao & Zech 2010). But expressions in Equation (8) appear less convenient than those in Equation (9), because there is a singular point when

*u*+

_{n}*c*= 0, although both types in the present research, with some programming consideration, are used.

### The HLL approximate Riemann solver

### The HLLC approximate Riemann solver

The HLLC scheme is a modification of the HLL scheme described in the previous section, whereby the missing contact and shear waves are restored. The wave structure in the HLLC approach is shown in Figure 3.

Then the HLL fluxes are suitable for the first and second components of the flux.

In the HLLC solver, the wave speed was estimated without considering the properties of bed material. Therefore, these properties need to be considered to improve the wave speed estimation.

### Modified HLLC scheme for the fluxes in the case with sediments

*λ*

_{2}in the subcritical case, and

*λ*

_{1}in the supercritical case), the classical HLLC flux formulation has been used for the first three components of the intercell flux . According to the research of Soares-Frazao & Zech (2010) the HLLC intercell flux in the case with sediments is calculated by the following relation: where and are mass flux and momentum flux in the normal direction and is momentum flux in the transverse direction. Then the flux expressions are evaluated as follows (Soares-Frazao & Zech 2010): where is the water level and Finally, the sediment flux in the normal direction is (Soares-Frazao & Zech 2010):

The waves speed and are related to the minimum and maximum wave speeds corresponding to the sediment movement at the interface (Soares-Frazao & Zech 2010).

### Formulations of the WAF method

The WAF approach for obtaining second-order extensions of the Godunov first-order upwind method has its origins in the random flux method (Toro 1999). WAF is a deterministic approach and has two versions. In the first version the intercell flux is an integral average of the physical flux across the full structure of the solution of a local Riemann problem. In the second version, an average state is computed, and the intercell flux, , resulted from evaluating the physical flux F at this state (Toro 1999).

### The original version of WAF

*N*is the number of waves in the solution of the Riemann problem and ,

*k*= 1, … , 4 are the normalized lengths of the segments In terms of the wave speeds , the weights are (Toro 1999): where is the Courant number for wave

*k*of speed . The alternative form of WAF flux is obtained by the substitution of from the Equation (26) into the Equation (24). where is the flux jump across wave

*k*of the Courant–Friedrichs–Lewy (CFL) number . In the CFL condition, is taken to be the minimum distance between two volume centres, but at boundaries, the distance between the centre of the element and the volume centre is used.

### TVD version of WAF schemes

*r*is calculated as follows (Toro 1999):

For the shallow water equations, *Q = h* for the acoustic waves, *Q = v* for the contact wave and *Q =* for the sediment are selected. The notation means the jump in *Q* across the wave *k* in the Riemann problem with data . Since the waves in the Riemann problem generally travel in different directions one must construct a function for each wave. For instance, in the triangular mesh, the left and right upwind states are computed as the mean value of states in , , , weighted by their relative areas according to Figure 6.

### Combining HLLC and WAF schemes (HLLC-WAF) in mobile bed conditions

### Boundary conditions

For the WAF scheme as the second-order method the application of boundary conditions is fundamentally the same as the Godunov method. The computational domain is assumed which discretized by *N* cells, so that cells *i* = *1*, … , *N* lie within the computational domain. Also the *M*th cell is assumed to be the boundary cell. For applying boundary conditions two fictitious cells next to each boundary are needed. For the left boundary the fictitious cells are denoted by *i = M−2* and *i = M−1* and for the right boundary they are denoted by *i = M + 1* and *i = M + 2*. Two types of boundary conditions are discussed.

*h*is the water depth,

*u*is the velocity component in the

*x*or

*y*direction considering boundary direction

**,**is the sediment depth (Toro 1999). In the domain of the test case all of the boundary forms using for the right upwind and for the left upwind side are depicted in Figure 7.

### Finite volume numerical model

*n*is the outward-directed unit vector normal to the boundary and

*F*is the flux vector. The integrand

*F*·

*n*is the normal flux across a surface. Equation (35) implies that the time rate of change in

*U*inside a control volume depends on the total flux through the surface of the volume plus the sum of sources within the volume. To facilitate solution of Equation (35) in this paper the approach of Zoppou & Roberts (2000) is followed so the rotation matrix

*T*is as below: Using the rotational invariance property of the two-dimensional shallow water equations, then Conversely, by applying Equation (37), Equation (35) is transformed into: where is the cell-base area, the

*j*-interface length and

*nb*the number of cell interfaces (three for triangular cells). As in Equation (5), vectors and express

*U*and

*F*in terms of normal and tangential velocity components and , attached to the considered interface. The transformation of coordinates is obtained via rotational matrix

*T*, where and are the components of the outward unit vector normal to the interface: For the source term , the friction slopes and are calculated using the Manning formula as follows:

Finally, due to the explicit nature of the scheme, stability is dictated by the CFL condition on the time step (Courant *et al*. 1967).

### Test case: dam-break flow in an erodible open channel with a sudden enlarged section

The present model is applied to model dam-break flow in an erodible channel. Palumbo *et al*. (2008) reported an experimental investigation on a model of dam-break flow in an erodible channel with a sudden enlargement. The flow domain consists of an upstream reservoir which has a length of 3 m, and there is a widening from 0.25 m to 0.50 m at 1 m downstream of the gate.

As can be seen in Figure 8, there is an initial sediment layer of 0.1 m all over the flume and the initial water depth is 0.25 m in the reservoir. In the downstream part of the flume, the sediment layer is saturated with water, but there is no water depth above the sediment layer. The sediments consist of coarse sand with median diameter *d*_{50} = 1.65 mm, deposited with a porosity . At the downstream end of the flume, there is a weir where the crest is at the top level of the sediment layer, so a free outflow occurs. The water level evolution and the final bed topography were recorded at specific points in the flume (Figure 8).

In the present model, the computational domain is divided into 8,035 triangular cells. The triangular grids are shown in Figure 9.

Comparisons between the measurements of Palumbo *et al*. (2008) and the present model results for the water level at different stations and the bed elevation in different sections are shown in Figures 10 and 11 respectively.

## RESULTS AND DISCUSSION

Comparing the numerical modelling results at the points U2, U3, U6 and U7 with the experimental test results (EXTR) shows very good convergence at the point U2. The start of the increase in water level, the maximum water level , the time related to the maximum water level, and the water level related to the end time are generally well reproduced by the results of the present model using the HLLC-WAF method.

According to the water level increase after t = 8 sec at the point U2 (Figure 10(a)), the present model result by the TVD-WAF method shows that kind of increase very well. At the point U3 (Figure 10(b)), the increase time of the water level results in the TVD-WAF method is different by about 0.3 sec comparing with EXTR. The maximum water level is a bit higher than EXTR but lowering the water level has been well reproduced as well as the water level related to the end time.

At the point U6 (Figure 10(c)), the shape of the water level increase between *t* = 1 sec and *t* = 2 sec has been well reproduced as well as the maximum water level in the model results by the HLLC-WAF method. Also the shape of the water level increase after *t* = 6 sec is simulated at t = 7 sec due to the influence of the second-order accuracy upwind method which has been employed in the present model.

At the point U7 (Figure 10(d)), the maximum water level has been well simulated. Probably it is resulted from some deviation in predicting bed level. Furthermore, the water level increase shape from the results of the TVD-WAF method is similar to EXTR.

Comparisons between measurements of Palumbo *et al*. (2008) and the present model results for bed elevation at different sections are shown in Figure 11. As is shown in Figure 11, the morphological change on the bed is rather well simulated at sections S3, S5 and S9.

At the section S3 (Figure 11(a)), the shape of the bed profile change in the model results is rather near to EXTR and the maximum and minimum level in the profile is well reproduced at *y* = 0.2 m and *y* = 0.42 m, considering the fact that the same coordinates in the test are *y* = 0.22 m and *y* = 0.38 m approximately. Only in two parts (*y* = 0.13 m and *y* = 0.25 m) are there little unexpected bed increases that are not considerable because of the good general shape of the bed profile in the numerical simulation at the section S3.

Also at the section S5 (Figure 11(b)) the shape of the bed profile is very well reproduced and compared to S3 and S9. The results of the present model in this section are the best ones comparing to the maximum and minimum bed level in the EXTR. So the values of the maximum and minimum bed level are approximately equal to the EXTR and they are located at *y* = 0.23 m and *y* = 0.44 m for the both numerical and experimental results.

At the section S9 (Figure 11(c)) general convergence of the shape of the bed profile between the numerical and experimental results is valid and discrepancies between *y* = 0.15 m and *y* = 0.3 m may result from the other effects of the flow that are not considered in the present model, such as the effect of secondary currents and changing of flow velocity in the vertical profile of the channel. Nonetheless, relatively acceptable predictions of flow in the present model results compared to the experimental results are provided.

For further evaluating use of the second-order accuracy TVD-WAF scheme in conjunction with the HLLC scheme, comparison of the present model results with the numerical results of Soares-Frazao & Zech (2010) are presented for the water level and sediment bed level in Figures 10 and 11.

According to the present model results, the complicated pattern of water level has been predicted better by the TVD-WAF method than by the HLLC method. The comparison shows that the results of the HLLC scheme are diffusive and smooth compared to the present model results. Because the HLLC method has first-order accuracy, the results of the HLLC method underestimate the water level and also the bed elevation. The HLLC-WAF scheme indicates better results in view of the convergence between numerical and experimental results. There is an improvement in the results of the HLLC-WAF scheme; however, there are discrepancies between the present model results and EXTR.

## CONCLUSIONS

A two-dimensional flow model over a mobile bed has been developed. The governing equations are solved in a coupled approach. The model is based on the finite volume method using a second-order approximate Riemann solver. The TVD version of the WAF scheme in conjunction with the HLLC scheme is capable of simulating changes in the water surface and morphology in the present model. Also the novel analysis of eigenvalues proposed by Soares-Frazao and Zech has been applied. It was successful because of utilization of sediment properties aimed at calculating the eigenvalues and also maximum wave speeds. Furthermore, application of the Ribberink formula for the sediment bed load presents suitable results for the flow over a mobile bed. Finally, the present model is applied for modelling dam-break flow over a mobile bed in a channel with sudden enlargement. Comparison of the model results to the experimental data and the result of the HLLC model shows that the present model can predict better the complicated water level and bed elevation than can the HLLC model. However, there are discrepancies between the present model results and experimental data especially in the bed elevation. The model results might improve by using other formulas for bed changes and considering the effect of three-dimensional flow.