Coupled finite-volume scheme with adapted Augmented Roe scheme for simulating morphological evolution in arbitrary cross-sections

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.

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 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.

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.

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.

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

In this approach, the sediment transport at the interface i − 1/2 has the value of either the upstream sediment transport rate Q_s,i−1 or the downstream one Q_s,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:

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 B_s is a representative cross-section that depends on the erosion assumptions, as detailed further.

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.

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

When the variation of the bed area ΔA_b 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 ΔA_b that is represented by a variation of the bed elevation z_bj of each node describing the cross-section.

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

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).

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 Q_s. 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.

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.

This test compares the models results to the analytical solution of a progressive slope aggradation or degradation in a prismatic channel. Both supercritical (slope S₀ = 5%) and subcritical (slope S₀ = 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.

For the supercritical case, the equilibrium solution for a liquid discharge Q = 100 m³/s and a sediment supply of Q_s = 1.517 m³/s corresponds to a bed slope S₀ = 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⁻⁶ m for UR, 4.55 × 10⁻⁶ m for CR, 9.16 × 10⁻⁶ 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.

For the supercritical case, the equilibrium solution for a liquid discharge Q = 100 m³/s and a sediment supply of Q_s = 0.111 m³/s corresponds to a bed slope S₀ = 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⁻⁶m for UR, 1.81 × 10⁻⁶ m for CR, 4.77 × 10⁻⁵ m for UL and 0.0084 m for CL.

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.

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 d₅₀ = 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.

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).

The results are shown in Figures 12–15 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.

As observed from Figures 12–15, 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⁻⁴ 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).

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 Δz_b over the cross-sections: the uniform erosion, i.e constant bed level variation Δz_b over the entire submerged part of the cross-section, and the non-uniform erosion, i.e. a distribution of Δz_b 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 (d₅₀) 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.

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.

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.

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