Abstract

A new coupled finite-volume scheme based on the Augmented Roe solver adapted to simulate morphological evolution of arbitrary cross-sections is presented. In pure hydrodynamic conditions, the Augmented Roe scheme has proven to provide accurate results and a constant discharge in steady-flow conditions. Here, this scheme is extended to solve the one-dimensional Saint-Venant–Exner system of equations written for arbitrary cross-sections. Therefore, new eigenvalues and source-term calculations are proposed to account for the irregular shape of the cross-sections. The performances of the proposed scheme are assessed by comparison with three different one-dimensional numerical models aimed at simulating morphological changes, with coupled or uncoupled approaches, and based on HLL or Roe-based flux calculations. Numerous test cases were examined, including water at rest, steady flows and transient flows for which experimental results exist. The results show that the proposed scheme provides stable and accurate results for a wider range of situations than is available with other classical models.

INTRODUCTION

Floods can be highly damaging natural disasters and their impact is often increased by the presence of sediments. In such cases, the strong morphological changes following an important flood completely reshape the riverbed and the landscape of the surrounding area. Predicting such floods and their morphological consequences is thus important to prevent human losses and damage. In solving this problem, 1D models could play a role as they have the advantages of needing less data and being faster than 2D models. Of course, they are not able to capture as many complex patterns as the 2D models. However, during a flood crisis, time is a key factor, if not the most important factor. In this respect, 1D models are still much faster than 2D models even with the latest developments in the use of GPU to increase the computation speed (Lacasta et al. 2014; Juez et al. 2016a).

In a depth-averaged framework, the shallow-water equations describing the water flow are used together with additional equations describing the morphological evolution of the bed following erosion and deposition of sediments. The links between these sets of equations, and consequently the impact of the presence of sediments on the water flow, are described in the literature using various approaches: (a) only mass conservation is considered for the sediments, and their effect on the flow density is neglected (e.g. Kassem & Chaudry 1998; Goutière et al. 2008); (b) the effect of the sediment concentration on the flow is considered by providing some additional coupling between the equations, either weak (e.g. Wu & Wang 2007; Garegnani et al. 2011) or strong (Juez et al. 2016b); (c) the moving sediments are represented using a transport layer, leading to fully coupled equations, which consider momentum exchanges between layers (e.g. Savary & Zech 2007; Zech et al. 2008); and (d) the moving sediments are considered as a separate and immiscible phase in the water flow with exchanges between the two phases providing a significant level of coupling between the equations (Greco et al. 2012).

In the present research, where only bed load transport is considered, the first approach (a) has been favored because in such cases, the concentration of sediment in the water column is low and can be neglected. The considered sediments consist of non-cohesive particles, such as sand, that are more often transported by bed load rather than suspension load. In this approach, the flow is considered as consisting of a layer of clear water, without sediment, flowing on top of a movable bed. The three main equations in this model are the two shallow-water equations describing mass and momentum conservations of the water layer, and the Exner equation describing the mass conservation of the sediment bed. Two types of numerical strategies can be considered to solve these equations. In the first, corresponding to uncoupled (or weakly coupled) schemes (Juez et al. 2014), the equations are solved using a two-step method. First, the shallow-water equations are solved without considering the sediment transport. Then, the Exner equation is solved without changing the hydraulic variables. In the second approach, the three equations are solved simultaneously, leading to a coupled solver in 1D (Goutière et al. 2008) and also in 2D (Murillo & Garcia-Navarro 2010).

The proposed scheme is a coupled solver based on the Augmented Roe (A-Roe) approach developed by Murillo & Garcia-Navarro (2014) that was extended to sediment transport and morphological evolution for arbitrary cross-sections. Particular attention is paid to the treatment of cross-sections of arbitrary shape which requires adaptations of the fluxes and source-terms calculations. The A-Roe approach has the advantage of allowing for accurate steady-state simulation where the computed discharge remains constant even in highly irregular topographies (Franzini & Soares-Frazão 2016). As sediment transport calculations depend on the discharge, this approach allows for the avoidance of excessive erosion or deposition on abruptly changing beds. In addition, the proposed scheme was complemented with a bank failure mechanism adapted from Swartenbroekx et al. (2010) in order to represent local mass failures and thus more accurate evolution of the cross-sections.

This paper is divided as follows. First, the governing equations used to represent the flow and sediment transport in a one-dimensional framework with arbitrary cross-sections are presented. Then, the proposed new A-Roe solver is presented together with three other more classical finite-volume-based approaches used for morphological evolution. In the next section, a series of tests are conducted, including water at rest, steady flows and transient flows for which experimental results exist to assess the performances of the proposed solver in comparison with the results provided by the other approaches. The applicability of the proposed A-Roe solver is also tested in cases where the cross-sections present initial steep banks, highlighting the importance of including a bank failure mechanism as proposed here. Finally, the performances of the proposed A-Roe solver for unsteady flows involving morphological evolution are discussed and conclusions are drawn.

GOVERNING EQUATIONS

The models used in the present study are based on three equations: the one-dimensional shallow-water equations or Saint-Venant equations for the hydrodynamics in arbitrary cross-sections, and the Exner equation for the mass conservation of the sediments.

The main assumptions for the Saint-Venant equations are an incompressible fluid, a hydrostatic pressure distribution, and, as stated before, a negligible concentration of sediments in the water. For a cross-section of arbitrary shape, they read 
formula
(1)
 
formula
(2)
where Q is the discharge, A the area of water, and S0 the bed slope defined by 
formula
(3)
with zb0 the bed level, more precisely the thalweg level. The integral gI1 represents the hydrostatic pressure thrust while the integral gI2 represents the longitudinal component of the lateral pressure due to the longitudinal width changes. They are integrals written as 
formula
(4)
 
formula
(5)
with η a local variable for the depth integration and b(x, η) the width of the channel at a determined depth (Figure 1). It is important to note that these two relations can be linked and this allows another way to write the source terms (Franzini & Soares-Frazão 2016). Indeed, the derivative of I1 with respect to x can be developed, using the Leibniz rule, as 
formula
(6)
Figure 1

Definition of the hydrostatic pressure.

Figure 1

Definition of the hydrostatic pressure.

This results in the following expressions linking I1, I2 and A: 
formula
(7)
 
formula
(8)
 
formula
(9)
In particular, considering a constant level , Equation (8) becomes 
formula
(10)
The friction is modeled using the Manning equation 
formula
(11)
with the roughness coefficient n (also called Manning coefficient) and the wetted perimeter P.
The Exner equation represents the conservation of the mass of sediments in the bed of the river. In arbitrary topographies, it reads 
formula
(12)
where Ab is the area of the bed of sediments (Figure 2), Qs the sediment transport rate and ɛ0 the bed porosity.
Figure 2

Definition of the area of erodible bed.

Figure 2

Definition of the area of erodible bed.

In this research, the bed porosity is considered constant in both time and space. The Exner equation can thus be rewritten as 
formula
(13)
Finally, it is necessary to add a closure equation to determine the sediment transport (Qs). Several relationships exist in the literature such as Meyer-Peter and Müller, Smart and Jäggi or Camenen and Larson, as summarized e.g. by El Kadi Abderrezzak & Paquier (2011) or García (2008). In the present paper, it has been decided to use the Meyer-Peter and Müller formula. This equation is only valid for non-cohesive particles moving as bedload, without suspension or wash load, and provides the unit-width sediment transport rate as 
formula
(14)
where s is the specific gravity of the sediments, d50 the median bed sediment diameter, R the hydraulic radius (defined as R=A/P), and θcr the non-dimensional shear stress for initial sediment motion (the value is set here to 0.047). According to García (2008), the sediment transport over the entire cross-section is obtained by multiplying the unit-width sediment transport by the width B of the channel at the level of the free surface. The final expression for the sediment transport is thus 
formula
(15)
The final system of equations can be written in vector form as 
formula
(16)
with U the vector of conserved variables, F the vector of fluxes and S the source terms (divided into topographical Sg and friction Sτ source terms). 
formula
(17)
 
formula
(18)
 
formula
(19)
 
formula
(20)
 
formula
(21)

NUMERICAL DISCRETIZATION

Finite-volume schemes are used to solve system (15). The discretized equations are written as 
formula
(22)
where the superscript n refers to the time step and the subscript i to the computational cell.

The details of the uncoupled and coupled solvers are presented below. In uncoupled solvers, the hydrodynamic equations are solved first, considering no evolution of the bed elevation due to erosion or deposition of sediment. Then the morphological changes are computed considering no changes in the hydraulic variables (wetted area A and discharge Q). In coupled solvers, all three equations are solved simultaneously using appropriate flux expressions and wave propagation speeds. The proposed scheme belongs to this latter category, with adapted wave-speeds and source-term treatments. Finally, all solvers presented below are subject to the CFL condition for the time step Δt, without any additional restriction due to the morphological evolution.

Uncoupled models

In this approach, the system is rewritten to separate the hydrodynamic (subscript h) and the sediment parts (subscript s). It, now, reads 
formula
(23)
 
formula
(24)
with 
formula
(25)
 
formula
(26)
 
formula
(27)
 
formula
(28)
 
formula
(29)
 
formula
(30)

The finite-volume scheme (22) is applied to the hydrodynamic system (23) only while Equation (24) is solved by simple finite-differences. For the flux calculations of the hydrodynamic system (22), several flux calculation schemes can be used. Two different schemes are presented here and used in the comparisons: (i) the Lateralized HLL scheme (Petaccia et al. 2013) and (ii) the Augmented Roe solver (Murillo & Garcia-Navarro 2014). Considering only pure hydrodynamic flows, the performances of these schemes were compared by Franzini & Soares-Frazão (2016) on a selection of test cases. The conclusion was that the Augmented Roe scheme is more accurate in terms of water level and discharge in the computational cells, but less robust than the Lateralized HLL (LHLL) scheme in the sense that the computed mass flux between cells is not more constant in the LHLL approach. The key features of these solvers are briefly recalled below.

LHLL solver

This solver is a variation of the HLL solver (Harten et al. 1983) developed for cross-sections of arbitrary shape (Petaccia et al. 2013). To compute more accurately the impact of irregular topographies, the topographical sources terms Sgh are included in the fluxes computations while the friction source terms Sτh remain separated, leading to a two-step scheme as 
formula
(31)
 
formula
(32)
where vectors Uh contain the hydrodynamic conservation variables defined in Equation (25) and Fh are the vectors containing information on the hydrodynamic mass and momentum fluxes as in Equation (26).
In the first step, Equation (31), the mass fluxes are calculated by taking the impact of the topography into account using a method based on the Nujic variation (Nujic 1995) – considering the water level instead of the water depth. In arbitrary topographies, the method is adapted considering defined as the wetted area of cell i at an average level (Figure 3). The resulting mass flux is written as 
formula
(33)
Figure 3

Definition of the area for the Lateralized HLL model.

Figure 3

Definition of the area for the Lateralized HLL model.

The wave celerities used in this model were defined by Davis (1988): 
formula
(34)
 
formula
(35)
with and
As regards the momentum fluxes in the first step (31), where the frictionless system is solved, the topographical source terms are included in the fluxes, leading to the lateralized ‘left’ and ‘right’ fluxes and linked by the relation . The momentum fluxes thus read, with : 
formula
(36)
 
formula
(37)

Then, in the second step described by Equation (32) the effect of the friction source terms is computed using the provisionally updated variables , leading to an implicit treatment of this term.

Augmented Roe solver

The Augmented Roe solver (A-Roe solver) is a modified version of the Roe solver (Roe 1981) in which the source terms are represented using a stationary wave allowing for energy conservation. It was developed for arbitrary sections by Murillo & Garcia-Navarro (2014) who showed that this approach results in a well balanced scheme, satisfying the C-property. In this solver, all the source terms (topographical and friction) are included in the fluxes computation. The finite-volume scheme (22) thus becomes 
formula
(38)
Following Murillo & Garcia-Navarro (2014), the fluxes are treated in a lateralized way and are written as 
formula
(39)
 
formula
(40)
where λ represents the wave celerity, β the source terms, α the wave strength and e the eigenvector of the hydrodynamic system, as detailed below. The wave celerities λ are given by 
formula
(41)
 
formula
(42)
with the overbar denoting the Roe averages defined as 
formula
(43)
 
formula
(44)
The matrix of eigenvectors with the eigenvectors being and . The wave strengths α are obtained using 
formula
(45)
 
formula
(46)
 
formula
(47)
The source terms are included in the flux calculation using a stationary wave representation given by parameter β in Equations (39) and (40), as follows: 
formula
(48)
 
formula
(49)
 
formula
(50)
In Equations (49) and (50), ζ is the expression for the source terms, with ζg representing the topographical part of this source term as follows: 
formula
(51)
The friction terms are written considering the spatial average values denoted by the overbar: 
formula
(52)
 
formula
(53)
To introduce the energy conservation in the model, the topographical source terms ζg are written as a linear combination with a weighting coefficient ω and 
formula
(54)
 
formula
(55)
The weighting coefficient is calculated by imposing the energy conservation 
formula
(56)
By combining Equation (56) with the momentum conservation equation, the weighting coefficient is obtained as 
formula
(57)

In a shock, such as a hydraulic jump, as the energy is not conserved but dissipated, .

Morphological evolution

The Exner equation describing the morphological evolution is solved using an upwind finite-difference scheme with the following discretization: 
formula
(58)

In this approach, the sediment transport at the interface i − 1/2 has the value of either the upstream sediment transport rate Qs,i−1 or the downstream one Qs,i. The choice is made according to the Froude number, as Savary & Zech (2007) showed that the eigenstructure of the Exner equation is such that information comes either from downstream or upstream depending on the Froude number, and on the derivative of the sediment transport with respect to the bed level. Such an approach consists in considering for the sediment transport the celerity of a single equation as follows:

If or : 
formula
(59)
else 
formula
(60)
The sediment discharge Qs in this model can be evaluated using any sediment transport equations using the flow variables at time n in an explicit way. In the applications presented here, the Meyer-Peter and Müller formula was used. Then to calculate the morphological evolution of the cross-section after evaluation of Equation (58), the variation ΔAb has to be distributed over the initial shape of the cross-section. Two possible approaches are described in the following section.

Coupled models

As stated before, in coupled models, all three equations of system (16) are solved simultaneously. The two approaches described above (the LHLL scheme and the A-Roe solver) are extended to include the third equation, representing the morphological evolution. This means that new wave propagation speeds need to be calculated to account for the effect of the sediment movement equation in the system. To be able to calculate the wave propagation speeds as the eigenvalues of the Jacobian matrix, the system of Equation (16) is rewritten considering a relation describing the variation of the bed elevation as a function of the variation of cross-section area as , where Bs is a representative cross-section that depends on the erosion assumptions, as detailed further.

System (16) now reads 
formula
(61)
with 
formula
(62)
 
formula
(63)
 
formula
(64)
 
formula
(65)
The Jacobian matrix of the new system (61) is 
formula
(66)
and the eigenvectors can be found by solving the eigenequation , resulting in 
formula
(67)
with 
formula
(68)
 
formula
(69)
 
formula
(70)

The solution to Equation (67) is not straightforward and the eigenvalues cannot be expressed analytically in a simple form, especially because of the irregular shape of the cross-sections. Novel extensions of the coupled Lateralized HLL and the coupled A-Roe schemes are presented below, as these two schemes require different methodologies to solve this problem.

Coupled Lateralized HLL

This model was developed by Goutière et al. (2008) for a rectangular channel. An extension to arbitrary topographies is presented here. The general finite-volume discretization (22) is rewritten as follows, with the topographical source terms included in the flux terms and Sτ representing the friction source term only: 
formula
(71)
To find the three wave celerities, i.e. the eigenvalues of (66) and thus the solutions of (67), it is first assumed that the largest eigenvalue is not influenced by the sediment (Lyn & Altinakar 2002). This yields 
formula
(72)
Then, using the definitions in Equations (68)–(70), the relations linking together the roots λ of a polynomial equation of the third degree are written 
formula
(73)
 
formula
(74)
Equations (73) and (74) can be solved for the wave celerities λ1 and λ2 
formula
(75)
 
formula
(76)
With these expressions for the wave celerities, the fluxes can now be written as follows, considering a lateralized treatment of the momentum flux: 
formula
(77)
 
formula
(78)
 
formula
(79)
 
formula
(80)
Using the definitions (72), (75) and (76), the wave celerities are written as 
formula
(81)
 
formula
(82)
 
formula
(83)

It is important to note that, if no sediment transport occurs, i.e. if Qs = 0, the flux expressions (77)­–(79) become the same as those of the LHLL scheme for hydrodynamics as λ1 = 0 according to (75).

Coupled Augmented Roe

This model is based on the 2D model developed by Murillo & Garcia-Navarro (2010). It is here adapted to one-dimensional cross-sections with arbitrary shapes, which implies an adaptation of the concept of energy conservation (Murillo & Garcia Navarro 2014; Franzini & Soares-Frazão 2016) to this particular context. The finite-volume discretization (22) is adapted to include all source terms and is written as 
formula
(84)
where the fluxes require adapted expressions for the wave celerity λ, the wave strength α and the source term stationary wave β 
formula
(85)
 
formula
(86)
The wave celerities λ are the solution of the system eigenequation 
formula
(87)
with coefficients ai given by (68)–(70). The solution requires the calculation of derivatives of the sediment transport rate, which can be evaluated using a wide range of sediment transport formula. To obtain simple expressions for these derivatives, the Grass formulation (Grass 1981) is used: 
formula
(88)
The methodology can however be adapted to any sediment transport formula. In the case of the Grass formulation, considering a constant coefficient G, the coefficients ai, expressed in terms of appropriate Roe averages , , (representative average width) become 
formula
(89)
 
formula
(90)
 
formula
(91)
where the Roe averages between cell i and i + 1 are estimated using 
formula
(92)
 
formula
(93)
 
formula
(94)
 
formula
(95)
In order to apply the energy conservation principles, the area at the interface is written using the weighting coefficient ω defined in Equation (57) as 
formula
(96)
Then, the eigenvalues are obtained directly as the solution of (87) as 
formula
(97)
 
formula
(98)
 
formula
(99)
with 
formula
(100)
 
formula
(101)
 
formula
(102)
It is important to note that the system (61) is hyperbolic only if . This condition is fulfilled in all the types of flows considered here. With the eigenvalues, the matrix of eigenvectors can be determined as 
formula
(103)
 
formula
(104)
and the wave strengths and the source terms are obtained using 
formula
(105)
 
formula
(106)
With the vector of source terms written as 
formula
(107)
As stated before, other sediment transport formulations can be used, such as the Meyer-Peter and Müller formula. In that case, the model has to be slightly modified to ensure the sediment mass conservation. First, the MPM formula has to be expressed in the same form as the Grass formula considering now a variable G coefficient (Murillo & Garcia-Navarro 2010) defined by 
formula
(108)
Even though this formulation of the Grass coefficient is a function of the water velocity, it will be considered constant in each cell at each time step leading, for interface i+1/2, to the following sediment fluxes in cells i and i + 1 
formula
(109)
 
formula
(110)
where Gi is calculated using the variables in cell i at time n. Expressions (109) and (110) are then included in (95) for the flux calculation.

Erosion in arbitrary cross-sections

When the variation of the bed area ΔAb has been obtained as the solution of (58) or (71), it is necessary to translate it into a variation of the cross-section shape (Figure 4). This can be done using different assumptions on the way the computed erosion is distributed over the cross-section, resulting in uniform or non-uniform erosion as detailed below. Consider the case illustrated in Figure 4: a cross-section A undergoes an erosion ΔAb that is represented by a variation of the bed elevation zbj of each node describing the cross-section.

Figure 4

Discretization of the cross-section.

Figure 4

Discretization of the cross-section.

According to Figure 4, the total bed area variation ΔAb is given by 
formula
(111)

The variation Δzbj of the bed elevation of each single node will depend on the assumptions made for the erosion mechanism, as described below. The partial width Bjk is defined as the free-surface width between node j and node k.

Uniform erosion

In this case, the erosion is uniformly distributed over all the immerged bed, yielding 
formula
(112)
Introducing this relation into the total bed area variation relation for any node j yields 
formula
(113)
 
formula
(114)
 
formula
(115)
Finally, the bed level variation of node j is obtained with 
formula
(116)
So, in the case of uniform erosion assumption, the equivalent average width Bs of Equation (61) is equal to width at the free surface B 
formula
(117)

Non-uniform erosion

In reality, erosion is not distributed uniformly over the entire cross-section. The amount of erosion is proportional to the bed shear stress exerted on each portion of the bed, which is in turn proportional to the local water depth. So, in the proposed non-uniform erosion model, the bed level variation is distributed according to the local water depth, in such a way that no erosion occurs if the water depth is equal to zero. This last characteristic renders the operator smoother as the transition between submerged and emerged bed is directly taken into account. So, the relation linking the bed level variation of two adjacent points j and k (Figure 4) is 
formula
(118)
Assuming a linear relation (n = 1), and introducing this into the total bed area variation relation (111) for node j yields 
formula
(119)
 
formula
(120)
 
formula
(121)
Finally, the level variation of point j is obtained from (121) as 
formula
(122)
In this method, the equivalent average width Bs is defined as 
formula
(123)
with zw the water level and zb0 the thalweg level. It can be observed that, in rectangular channels, the non-uniform erosion gives the exact same results as the uniform erosion.

Bank failure

In nature, when erosion starts at the toe of the bank of a river, the bank slope can become steeper than the stability limits of the constituting material. This will then result in a failure of the bank. However, in numerical schemes where only bed deepening is represented as a result of the governing equations and using processes such as those described above, global mass failures and the resulting channel widening cannot be represented. To account for such mass failures, a bank failure mechanism has to be included in the model. The bank failure mechanism proposed here is derived from Zech et al. (2008) and Swartenbroekx et al. (2010) and is based on a tilting of the unstable portion of the bank. For irregular cross-sections described by straight segments, the stability of each segment is checked and the tested segment is tilted if its slope exceeds the stability limit. Doing this, an adjacent segment can become unstable. So, the whole process is repeated until complete stability of the cross-section is achieved. The process is described in Figure 5: segment BC is initially unstable and after the tilting, segment B′C′ is stable but segment AB′ has become unstable, and will thus be tilted in the next iteration. The convergence of the procedure was demonstrated by Swartenbroekx et al. (2010).

Figure 5

Bank failure operator.

Figure 5

Bank failure operator.

As only non-cohesive sediments are considered, the slope limit depends solely on the internal friction angle of the material. It is important to add that, as the bank failure operator is operated separately from the rest of the model, it only changes the shape of the section without affecting the sediment transport Qs. The sediment deposited at the toe of the bank after a failure event will be transported by the flow during the next time step. No dynamical effects are thus taken into account.

RESULTS

In this section, the results of the novel coupled A-Roe model are compared with the solutions provided by the three other models. Here, the models are abbreviated as follows: UH for uncoupled LHLL, UR for uncoupled A-Roe, CH for the coupled LHLL and CR for the coupled A-Roe model. All tests were run with a CFL number equal to 0.9.

Water at rest

To assess the ‘well balancedness’ of the models, it is important to check if they are able to keep water at rest. The selected tests consider a trapezoidal channel (Murillo & Garcia-Navarro 2014) with a shape defined by 
formula
(124)
with B0 the width of the channel at h = 0 m, as represented in Figure 6 for the CR model. All models are able to maintain the water at rest with a maximum error on the discharge of 10−13 m³/s and, therefore, respect the C-property (Vazquez-Cendon 1999), so only the results for the CR model are illustrated, as they are very similar for the four models.
Figure 6

Results for the water at rest using the CR model.

Figure 6

Results for the water at rest using the CR model.

Equilibrium slopes

This test compares the models results to the analytical solution of a progressive slope aggradation or degradation in a prismatic channel. Both supercritical (slope S0 = 5%) and subcritical (slope S0 = 0.5%) cases are investigated. For each case, three initial configurations are studied (Figures 7 and 8): (1) starting at the equilibrium slope, with equilibrium upstream sediment supply, (2) starting from a steeper slope (6% and 0.8% for the supercritical and subcritical cases, respectively), with a sediment supply corresponding to the final equilibrium to be reached after erosion and thus slope degradation, (3) starting from a milder slope (4% and 0.25% for the supercritical and subcritical cases, respectively), with a sediment supply corresponding to the final equilibrium to be reached after deposition and thus slope aggradation.

Figure 7

Results for the equilibrium of a supercritical slope: (a) initial conditions and (b) final computed bed elevation. All tested models (CR-CH-UR-UH) provide the same results.

Figure 7

Results for the equilibrium of a supercritical slope: (a) initial conditions and (b) final computed bed elevation. All tested models (CR-CH-UR-UH) provide the same results.

Figure 8

Results for the equilibrium of a subcritical slope: (a) initial and (b) final.

Figure 8

Results for the equilibrium of a subcritical slope: (a) initial and (b) final.

Supercritical

For the supercritical case, the equilibrium solution for a liquid discharge Q = 100 m³/s and a sediment supply of Qs = 1.517 m³/s corresponds to a bed slope S0 = 5%. All models provide the exact solution, regardless of the initial condition. Figure 7(a) shows the three initial conditions: the upper line corresponds to the initial steeper slope, the middle line to the equilibrium slope and the lower line to the initial milder slope. All tested schemes were run for 150,000 time steps with a CFL number of 0.9. The correct equilibrium slope is obtained as illustrated in Figure 7(b) showing the final situation, with the following RMSE on the bed elevation: 4.62 × 10−6 m for UR, 4.55 × 10−6 m for CR, 9.16 × 10−6 m for UL and 0.039 m for CL. These values do not depend on the initial conditions. Moreover, it appears that the largest error is obtained for the CL model.

Subcritical

For the supercritical case, the equilibrium solution for a liquid discharge Q = 100 m³/s and a sediment supply of Qs = 0.111 m³/s corresponds to a bed slope S0 = 0.5%. All models but the coupled LHLL (CH) give the same, and correct, results (Figure 8), regardless of the initial condition. This means that starting from an initial slope steeper or milder than the equilibrium one as illustrated in Figure 8, the correct final equilibrium slope is obtained, except for the CH scheme. For this case, the RMSE are: 1.81 × 10−6m for UR, 1.81 × 10−6 m for CR, 4.77 × 10−5 m for UL and 0.0084 m for CL.

The failure of the CH model can be explained by looking at the sediment mass flux expression: 
formula
(125)

The last term of the equation, which constitutes the diffusive part of the flux expression, only cancels out if the bed elevation is constant. This problem, already pointed out by Goutière et al. (2008) for rectangular channels, can be solved by adding a term, which is a function of the equilibrium slope, in the diffusive term. However, this solution proposed by Goutière et al. (2008) was not adopted here, as it cannot be extended to arbitrary topographies with abruptly changing bed slopes. Therefore, the CH model should not be considered for general applications and the proposed CR model can thus be considered as an improvement of coupled models for this type of problem.

Dam break flow over a flat erodible bed

Experiments of dam-break flow over an initially flat mobile bed were presented by Spinewine & Zech (2007). The channel is 6 m long and 0.25 m wide. Upstream, the first 3 m are initially filled with 0.35 m of water. At , the dam is removed and the water starts flooding the flume. The sediment bed is composed of PVC particles with a median diameter d50 = 3.9 mm and a specific mass s = 1.58. During the experiment, three flow regions were observed: a layer of pure water, a layer of sediment transport composed of a mixture of water and sediments and the motionless bed. Therefore, two interfaces were captured: the interface between the flow of pure water and the mixture, and the interface between the mixture and the motionless bed. However, the models presented in this paper do not compute a layer of moving sediments. They only provide the bed elevation after erosion or deposition. This level is assumed to correspond to an intermediate value between the lower and upper limits of the moving sediment layer. So, numerical results are considered to be in good agreement with the measurements when the computed bed elevation is located between the two recorded interfaces. The results are presented for t= 1 s (Figure 9). They show that all four models give similar erosion, with values between the upper and lower limits of the moving sediment layer (dots in the figure). However, the CR model shows a more regular bed evolution than the other three models that present spurious bed level variations (zoom in Figure 9) and thus appear less stable. These spurious bed level variations are reduced but not eliminated by reducing the time-step.

Figure 9

Dam-break flow over a mobile bed. Comparison simulations (lines) and experimental results (dots).

Figure 9

Dam-break flow over a mobile bed. Comparison simulations (lines) and experimental results (dots).

Dike failure by overtopping

This experiment studying erosion and failure of a sand dike was conducted at the Hydraulics Laboratory of the Université catholique de Louvain, Belgium (Van Emelen et al. 2015). The dike was built in a 10 m long, 0.2 m wide and 0.3 m high horizontal flume. The dike is of trapezoidal shape, 0.20 m high with slopes of 1 V:2H (Figure 10) and is made of sand with a mean diameter of 0.61 mm. During the erosion process, both water and bed level evolutions were recorded (Figure 11).

Figure 10

Dike failure: experimental setup.

Figure 10

Dike failure: experimental setup.

Figure 11

Comparison between measured and computed bed and water levels for the uncoupled Lateralized HLL model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 11

Comparison between measured and computed bed and water levels for the uncoupled Lateralized HLL model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

The results are shown in Figures 1215 for four time steps: 8 s, 24 s, 64 s and 200 s (t= 0 s corresponds to the initial overtopping of the dike) for each model. At each of the considered time instants, the computed water and bed levels (continuous lines) are compared to the measurements (discrete symbols). The results are then summarized in Figure 15 where the evolution of the average absolute error is presented for the bed level, the water level and the water depth. The errors are computed for t = 8, 14, 24, 64 and 200 s. Finally, Figure 16 represents the signed error for the bed level, water level and water depth computed as the difference between the predicted and measured values.

Figure 12

Comparison between measured and computed bed and water levels for the uncoupled Augmented Roe model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 12

Comparison between measured and computed bed and water levels for the uncoupled Augmented Roe model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 13

Comparison between measured and computed bed and water levels for the coupled Lateralized HLL model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 13

Comparison between measured and computed bed and water levels for the coupled Lateralized HLL model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 14

Comparison between measured and computed bed and water levels for the coupled Augmented Roe model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 14

Comparison between measured and computed bed and water levels for the coupled Augmented Roe model at (a) 8 s, (b) 24 s, (c) 64 s and (d) 200 s. Measured bed level , computed bed level , measured water level and computed water level .

Figure 15

Evolution of the average absolute error for the dike failure considering the bed level, the water level and the water depth.

Figure 15

Evolution of the average absolute error for the dike failure considering the bed level, the water level and the water depth.

Figure 16

Average error for the bed level, water level and water depth computed as the difference between the predicted and measured values.

Figure 16

Average error for the bed level, water level and water depth computed as the difference between the predicted and measured values.

As observed from Figures 1215, the results are similar for the coupled and uncoupled version of the A-Roe model (CR and UR) for 8 s, 24 s and 64 s. All models also reach similar results at 64 s. However, the intermediary results (at 24 s) show significant differences. CR erodes the dike faster than the experiment while UH erodes slower. Only CH is able to capture the correct crest level evolution. However, all models fail to represent the correct behavior of the bed and water level at 64 s when antidunes start to form in the experiment.

In Figures 15 and 16, it can be observed that the error remains constant for the coupled A-Roe model (CR) for the bed level prediction, contrarily to the other models where the error progressively decreases. In addition, LHLL-based models (CH and UH) provided better bed evolution than A-Roe-based models (CR and UR) that in turn provided better predictions of the water depth. Figure 16 shows that all models have the same error for the prediction of the dike erosion. First, at 8 s, the models do not transport enough sediments and the erosion is underestimated. Then, at 14 and 24 s, the erosion of the dike is overestimated before being again underestimated at 64 s. Finally, for the prediction of the volume of water, the coupled Roe model (CR) gives the best prediction with an average error of less than 2 × 10−4 m after 24 s.

These conclusions in terms of error can be compared to the conclusions drawn in Franzini & Soares-Frazão (2016), where similar schemes were compared on pure hydrodynamic cases with complex topographies. The lowest errors in terms of discharge and water level predictions were obtained for the schemes using the energy balance (HLLS and A-Roe). In the present cases involving morphological evolution, it appears that the major part of the error is due to the computation of the elevation that is highly dependent on the sediment transport relations, as shown for example by Van Emelen et al. (2015).

Dam-break flow in a trapezoidal channel

Experiments of a dam-break flow in an initially trapezoidal channel made of coarse sand with steep banks were conducted at the Hydraulics Laboratory of the Université catholique de Louvain, Belgium (Soares-Frazão et al. 2007). This test aims at validating the procedures described above for the distribution of bed level variations Δzb over the cross-sections: the uniform erosion, i.e constant bed level variation Δzb over the entire submerged part of the cross-section, and the non-uniform erosion, i.e. a distribution of Δzb over the cross-sections according to the actual bed shear stress. Two types of results are presented: uniform erosion without any bank failure mechanism, denoted UE-NBF, and non-uniform erosion with bank failure mechanism, denoted NUE-BF. The initial conditions consist of a semi-trapezoidal, 5.7 m long channel as sketched in Figure 17, the reservoir being filled with a 15 cm water depth. The sediment used is sand with a median diameter (d50) of 1.8 mm and a specific density s of 2.615. The angles of repose are 37° for submerged sand and 85° for humid sand.

Figure 17

Dam-break on a trapezoidal channel: experimental setup. (a) Top view, (b) cross-section and (c) initial conditions.

Figure 17

Dam-break on a trapezoidal channel: experimental setup. (a) Top view, (b) cross-section and (c) initial conditions.

In Figure 18, the results are shown for cross-section S2 located at x = 0.5 m downstream from the gate using the coupled A-Roe approach. It can be easily observed that for the bank evolution the NUE-BF approach gives better results than considering a classical, simpler, UE-NBF approach. Only the NUE-BF is able to capture the variation in the slope of the banks. It is thus important to consider local 2D phenomena in the models by means of non-uniform erosion and bank failure.

Figure 18

Dam-break flow in a trapezoidal channel: comparison between UE-NBF and NUE-BF approaches using the coupled Roe.

Figure 18

Dam-break flow in a trapezoidal channel: comparison between UE-NBF and NUE-BF approaches using the coupled Roe.

DISCUSSION AND CONCLUSIONS

A new coupled finite-volume scheme based on the Augmented Roe solver with energy balance is presented. It was adapted to simulate morphological evolution of arbitrary cross-sections. The proposed model has also been compared to three more classical models, coupled or uncoupled: a coupled LHLL, an uncoupled LHLL model and an uncoupled Augmented Roe model. First, the distinction between uncoupled (or weakly coupled) and fully coupled models has been described. While the former solve the equations in two steps (shallow water equations then Exner equation), the latter solve them together resulting in a more complex scheme. Secondly, from the point of view of the model formulation, as the HLL method uses approximate value for the celerities, it results in simpler expression than Roe-based models.

The coupled Augmented Roe was compared to the other models using different test cases. These tests were both numerical and experimental. It could be observed that, although much simpler, even the uncoupled models provided reasonably good results on most of the test cases. However, the proposed CR model provided more regular and realistic results than the other models, i.e. without spurious oscillations. The tests also highlighted the fact that the coupled LHLL model failed to reproduce bed aggradation or degradation in subcritical flow conditions.

Finally, the necessity to consider non-uniform erosion distribution and a bank failure model for arbitrary cross-sections was highlighted by studying a dam-break flow in a semi-trapezoidal channel. It was shown that these two parameters influence greatly the morphological changes created by a dam-break flow in arbitrary topographies.

ACKNOWLEDGEMENTS

This work was supported by the Fonds National de la Recherche Scientifique, Belgium.

REFERENCES

REFERENCES
Davis
S. F.
1988
Simplified second-order Godunov-type methods
.
SIAM Journal on Scientific and Statistical Computing
9
,
445
473
.
El Kadi Abderrezzak
K.
&
Paquier
A.
2011
Applicability of sediment transport capacity formulas to dam-break flows over movable beds
.
Journal of Hydraulic Engineering
137
(
2
),
209
221
.
García
M. H.
2008
Sediment transport and morphodynamics
. In:
Sedimentation Engineering Processes, Measurements, Modeling and Practice
(
García
M. H.
, ed.).
ASCE Manuals and Reports on Engineering Practice 110
,
ASCE
,
Reston, VA
, pp.
21
163
.
Garegnani
G.
,
Rosatti
G.
&
Bonaventura
L.
2011
Free surface flows over mobile bed: mathematical analysis and numerical modeling of coupled and decoupled approaches
.
Communcations in Applied and Industrial Mathematics
3
(
1
),
1
22
.
Goutière
L.
,
Soares-Frazão
S.
,
Savary
C.
,
Laraichi
T.
&
Zech
Y.
2008
One-dimensional model for transient flows involving bedload sediment transport and changes in flow regimes
.
Journal of Hydraulic Engineering
134
,
726
735
.
Grass
A. J.
1981
Sediments Transport by Waves and Currents
.
SERC London Centre for Marine Technology, Report No. FL29
.
Greco
M.
,
Iervolino
M.
,
Vacca
A.
&
Leopardi
A.
2012
Two-phase modelling of total sediment load in fast geomorphic transients
. In:
River Flow 2012
Proc. IAHR Int. Conf. on Fluvial Hydraulics
,
Costa Rica
,
5–7 September 2012
.
CRC Press
,
Vol. 1
, pp.
643
648
.
Juez
C.
,
Murillo
J.
&
García-Navarro
P.
2014
A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed
.
Advances in Water Resources
71
,
93
109
.
Juez
C.
,
Lacasta
A.
,
Murillo
J.
&
García-Navarro
P.
2016a
An efficient GPU implementation for a faster simulation of unsteady bed-load transport
.
Journal of Hydraulic Research
54
(
3
),
275
288
.
Juez
C.
,
Murillo
J.
&
García-Navarro
P.
2016b
One-dimensional Riemann solver involving variable horizontal density to compute unsteady sediment transport
.
Journal of Hydraulic Engineering
142
(
3
),
article number 04015056
.
Kassem
A. A.
&
Chaudry
M. H.
1998
Comparison of coupled and semicoupled numerical models for alluvial channels
.
Journal of Hydraulic Research
124
,
794
802
.
Lacasta
A.
,
Morales-Hernández
M.
,
Murillo
J.
&
García-Navarro
P.
2014
An optimized GPU implementation of a 2D free surface simulation model on unstructured meshes
.
Advances in Engineering Software
78
,
1
15
.
Lyn
D. A.
&
Altinakar
M.
2002
St. Venant–Exner equations for near-critical and transcritical flows
.
Journal of Hydraulic Engineering
128
(
6
),
579
587
.
Murillo
J.
&
Garcia-Navarro
P.
2010
An Exner-based coupled model for the two-dimensional transient flow over erodible bed
.
Journal of Computational Physics
229
,
8704
8732
.
Petaccia
G.
,
Natale
L.
,
Savi
F.
,
Velickovic
M.
,
Zech
Y.
&
Soares-Frazão
S.
2013
Flood wave propagation in steep mountain rivers
.
Journal of Hydroinformatics
15
,
120
137
.
Soares-Frazão
S.
,
Le Grelle
N.
,
Spinewine
B.
&
Zech
Y.
2007
Dam-break induced morphological changes in a channel with uniform sediments: measurements by a laser-sheet imaging technique
.
Journal of Hydraulic Research
45
(
Suppl. 1
),
87
95
.
Spinewine
B.
&
Zech
Y.
2007
Small-scale laboratory dam-break waves on movable beds
.
Journal of Hydraulic Research
45
(
Suppl. 1
),
73
86
.
Swartenbroekx
C.
,
Soares-Frazão
S.
,
Staquet
R.
&
Zech
Y.
2010
Two-dimensional operator for bank failures induced by water-level rise in dam-break flows
.
Journal of Hydraulic Research
48
,
302
314
.
Van Emelen
S.
,
Zech
Y.
&
Soares-Frazão
S.
2015
Impact of sediment transport formulations on breaching modelling
.
Journal of Hydraulic Research
53
,
60
72
.
Wu
W.
&
Wang
S.
2007
One-dimensional modeling of dam-break flow over movable beds
.
Journal of Hydraulic Engineering
133
,
48
58
.
Zech
Y.
,
Soares-Frazão
S.
,
Spinewine
B.
&
Le Grelle
N.
2008
Dam-break induced sediment movement: experimental approaches and numerical modelling
.
Journal of Hydraulic Research
46
,
176
190
.