Abstract

We propose a finite-volume model that aims at improving the ability of 2D numerical models to accurately predict the morphological evolution of sandy beds when subjected to transient flows like dam-breaks. This model solves shallow water and Exner equations with a weakly coupled approach while the fluxes at the interfaces of the cells are calculated thanks to a lateralized HLLC flux scheme. Besides describing the model, we ran it for four different test cases: a steady flow on an inclined bed leading to aggradation or degradation, a dam-break leading to high interaction between the flow and the bed, a dam-break with a symmetrical enlargement close to the gate and a dam-break in a channel with a 90° bend. The gathered results are discussed and compared to an existing fully coupled approach based on HLLC fluxes. Although both models equally perform regarding water levels, the weakly coupled model looks to better predict the bed evolution for the four test cases. In particular, its results are not affected by an excessive numerical diffusion encountered by the coupled model. Moreover, it usually better estimates the amplitudes of the maximum deposits and scours. It is also more stable when subject to high bed–flow interaction.

HIGHLIGHTS

  • New two-dimensional weakly-coupled scheme for fast transient flows over mobile beds.

  • Validation using three experimental data sets.

  • Comparison with classical coupled schemes to demonstrate the improvements of the weakly-coupled scheme.

INTRODUCTION

River floods have always been a major driver of abrupt changes in the environment. Although they constitute a perfectly natural phenomenon, humans fear them for the damage they can occasion to land and infrastructures. Not only can fast flows lead to significant land erosion and deterioration of structures, but the increase in the amount of sediments driven by the river in such situations can lead to significant morphological changes, like those caused by the 1996 Lake Ha! Ha! breakout flood for instance (Brooks & Lawrence 1999). For these reasons, particular attention has been dedicated by many authors to the numerical simulation of transient flows on mobile beds in the past few years (e.g., Wu 2004; Zech et al. 2009; Soares-Frazão & Zech 2010; Roushangar et al. 2011; Juez et al. 2013; Hou et al. 2018; Barzgaran et al. 2019; Di Cristo et al. 2019).

Because of the meandering morphology of rivers and the turbulent behaviour of fast transient flows, 3D finite-volume models are probably the most suited to simulate such flows. However, the high computational costs they imply make them irrelevant to simulate flows accurately on large space scales. Consequently, most modellers have been developing 2D depth-averaged models, using the so-called shallow water equations to describe the hydrodynamic processes involved.

As long as sediments are transported as bed load only, the Exner equation is often considered to describe bed evolution in the model. Regarding its implementation in numerical models, a closure equation is required to predict the solid discharge. A wide range of formulations exist in the literature, for example, Meyer-Peter & Müller (1948), Engelund & Hansen (1967), Wong & Parker (2006) and Camenen & Larson (2008) with, however, great variability in the resulting sediment transport as highlighted by Van Emelen et al. (2015).

Besides the choice of the closure equation, one of the key aspects much discussed by scientists then concerns the coupling between the hydrodynamics, expressed through the shallow water equations, and the bed evolution, described by the Exner equation. Goutière et al. (2008) used a coupled finite-volume scheme with HLLC fluxes in 1D cases to simulate flows over steep sediment slopes. Soares-Frazão & Zech (2010) extended this approach to 2D cases. However, even though good prediction of water levels could be obtained, the error on bed levels often remained unsatisfactory. In particular, Goutière et al. (2008) pointed out a diffusive term in the sediment mass flux as the reason for the flattening of the bed if the simulation was to last forever.

A possible solution to this issue would be to decouple the shallow water equations from the Exner one. Nevertheless, the definition of coupling needs some clarification as it may be different from one author to another. Some would say a system is coupled if the governing equations of flow and sediment continuity are completely solved within the same time step (Kassem & Chaudhry 1998; Cao et al. 2002). For this purpose, we prefer to talk about synchronous/asynchronous resolution. Garegnani et al. (2011) then used a non-dimensional analysis to make the distinction between coupled and uncoupled models, simplifying the system by removing all terms related to long-term morphological processes. In this sense, we will rather talk about simplified systems than decoupled ones. Finally, Franzini & Soares-Frazão (2018) used a third definition, in which decoupling actually means the separate resolution of the shallow water equations on one hand and of the Exner equation on the other hand. This is also the definition that we will use throughout this paper. Nevertheless, this separate resolution can be performed during the same time step (synchronous). In particular, the Jacobian matrix used to compute the hydrodynamic fluxes at the interfaces between cells only considers the purely hydrodynamic flow, while the morphodynamic fluxes are solved with a full upwind scheme depending on the Froude number. Using different flux formulations such as the Augmented Roe approach developed by Murillo & García-Navarro (2010), Franzini & Soares-Frazão (2018) tested both coupled and fully uncoupled schemes, but could not conclude on the clear superiority of any of the tested schemes for all the tested cases.

An alternative weakly coupled approach based on the Augmented Roe scheme was successfully developed by Juez et al. (2014) and applied to 1D steep sloping beds (Juez et al. 2017). Rather than using the Froude number, a sediment-related numerical characteristic was defined. However, this characteristic was not derived from the Jacobian matrix and this is why we speak about weakly coupling rather than decoupling. In this paper, we adopted a synchronous and simplified weakly coupled approach similar to Juez et al. (2014), but adapted to the lateralized HLLC scheme. The paper is organized as follows. First, the governing equations are presented, as well as the coupled and weakly coupled discretization with lateralized HLLC fluxes. Then, the two schemes are applied to four different test cases, of which three are accompanied by experimental data: (i) a steep-slope channel evolving following over- or under-supply of sediments, (ii) a 1D dam-break causing high interaction between the flow and the bed, (iii) a dam-break flow in a wide erodible floodplain, and (iv) a dam-break flow in a channel with a 90° bend and mobile bed. Finally, conclusions about the performances of the two approaches are provided.

GOVERNING EQUATIONS AND NUMERICAL RESOLUTION

As previously stated, the 2D shallow water equations describe the water flow and lead to three equations: one for the conservation of mass and two others for the conservation of momentum in two directions:
formula
(1)
formula
(2a)
formula
(2b)
Then, the mass conservation of sediments is ensured through the Exner equation:
formula
(3)
with h being the water depth, and the unit discharge components in the x and y directions, respectively, g the gravitational acceleration, and the normal momentum fluxes in the x and y directions, respectively, the transverse momentum flux, and the bed slope in the x and y directions, respectively, u and v the depth-averaged velocity components in the x and y directions, respectively, and the friction slope components explicitly calculated with the Manning coefficient n, the bed level, and the unit solid discharge components and the bed porosity.
Equations (1) to (3) can be grouped together under the following vector form, with the detailed expressions of σ and μ the momentum fluxes given in (5):
formula
(4)
where,
formula
(5a,b)
formula
(5c,d)
The solid discharge can be expressed through a lot of different closure equations, but we used the well-known formulation of Meyer-Peter and Müller (MPM) as it gave rather good results in comparison with the experimental tests described hereafter. The two components of the sediment transport rate are thus expressed as follows:
formula
(6)
formula
(7)
where s is the ratio of sediment density to water density, is a representative diameter of the grains, considered to be homogeneous, n is the Manning friction coefficient and is the non-dimensional critical bed shear stress beyond which sediments get transported by the flow. The MPM formula sets this threshold at 0.047.

Coupled model (CM)

The equations in system (4) are solved simultaneously. The system can be rewritten in a lateralized form as proposed by Fraccarollo et al. (2003). Costanzo et al. (2002), indeed, showed that this formulation gives more accurate results when the topography is irregular, which is what is expected with a mobile bed. A first step in this lateralization procedure is to include the topographical source term in the left-hand side of the equation:
formula
(8)
where
formula
(9a,b,c)
These equations are then discretized according to a first-order finite-volume scheme that is run on a 2D mesh composed of unstructured triangular cells, as generated, for example, by the software Gmsh (Geuzaine & Remacle 2009). In a finite-volume scheme, fluxes are calculated at each interface between two cells (Figure 1) in a local system of coordinates attached to the interface (,) so that the local variation of F along the edge may be considered as zero. In this way, the use of the rotational invariance property (Toro 1999; Guinot 2008) with a proper rotational matrix T, as defined in (10), can lead to the consideration of an augmented 1D problem for each edge with:
formula
(10)
so that, with the aim of computing the Jacobian matrix of our homogeneous system, we disregard the friction source terms and get:
formula
(11)
with
formula
(12a,b)
Figure 1

Local coordinates system.

Figure 1

Local coordinates system.

where and being, respectively, the normal and transverse components of the outward unit vector to the edge, and the depth-averaged velocity components along and , respectively, and the unit flow discharge components along and , respectively, the unit solid discharge component along while and are the momentum fluxes normal and transverse to the edge, respectively.

We can then finally express the Jacobian matrix of our system as so that (11) becomes:
formula
(13)
This last expression can then be used to define a pseudo-Jacobian matrix that includes the slope terms to derive the eigenvalues of the system and thus, the characteristics along which the information propagates:
formula
(14)
with being the water wave celerity.
As shown by Soares-Frazão & Zech (2010), considering as positive, the eigenvalues of this pseudo-Jacobian matrix can be approximated and sorted in ascending order as follows:
formula
(15a)
formula
(15b)
formula
(15c)
formula
(15d)
where depends on the considered bedload formulation. As stated before, the MPM formulation only was used throughout this paper so that:
formula
(16)

In Equation (15)a), only and are influenced by the sediments. Cordier et al. (2011) showed in a 1D analysis that even the largest eigenvalue of the pseudo-Jacobian matrix could actually be affected by the evolution of the bed level. Nevertheless, the use of approximated eigenvalues as in Equation (15)a) is computationally far cheaper than the analytical resolution of a third-order equation for each cell, and consequently constitutes the set of equations that we will use hereafter.

Here, we will solve system (11) thanks to a first-order finite-volume scheme. Its discretization reads:
formula
(17)
with being a variable time step that verifies the stability criterion (CFL) given by Courant et al. (1967), the cell-base area, j the considered cell interface, the -interface length, the number of cell interfaces.

Such a first-order finite-volume method relies on the definition of constant states U defining Riemann problems. For the pseudo-Jacobian matrix (14), the wave structure is as illustrated in Figure 2, where one intermediate wave corresponds to the sediment propagation, and the other one is a contact discontinuity. As described in Soares-Frazão & Zech (2010), an adapted HLLC approach is used for the fluxes, taking into account these intermediate waves.

Figure 2

Riemann problem considering the presence of sediments: (a) subcritical flow, (b) supercritical flow; after Soares-Frazão & Zech (2010).

Figure 2

Riemann problem considering the presence of sediments: (a) subcritical flow, (b) supercritical flow; after Soares-Frazão & Zech (2010).

According to Harten et al. (1983), mass and normal momentum fluxes can then be calculated thanks to the extreme hydrodynamic characteristics, i.e., and , assuming only one intermediate state . The presence of sediments is not completely disregarded since does contain sediment-related information. These fluxes can be computed as follows:
formula
(18)
formula
(19)
with and being the water levels of the right and left cells, respectively, and with
formula
(20)
formula
(21)
The calculation of the transverse momentum flux accounting for the contact discontinuity can be done as suggested by Toro et al. (1994):
formula
(22)
Finally, the sediment mass flux can be obtained by writing the jump relations between the intermediate states (Soares-Frazão & Zech 2010). Eventually, by assuming that the sediment-related information is only contained in the two first characteristics, as suggested by Lyn & Altinakar (2002), both in subcritical and supercritical regimes, we can express the flux as follows:
formula
(23)
with
formula
(24)
formula
(25)

As observed by Goutière et al. (2008), the last term of (23) will eventually lead to a horizontal topography if the discharge is kept non-zero, or at least to excessive diffusion and thus inaccurate representation of erosion and deposition processes.

Weakly coupled model (WCM)

The model developed by Juez et al. (2014) decouples hydrodynamic equations from the Exner equation describing the bed morphological evolution. Only hydrodynamic quantities involved in mass and momentum conservation laws are now coupled in system (26), whereas the Exner equation (Equation (3)) is handled separately:
formula
(26)
where,
formula
(27a,b)
formula
(27c,d)
System (26) can then be lateralized and transformed in order to be attached to the local coordinates system of the interfaces in the very same way that we did with the CM. Considering a vector , system (13) remains valid, but the pseudo-Jacobian matrix now reads:
formula
(28)
and actually corresponds to the hydrodynamic Jacobian matrix so that only three pure hydrodynamic characteristics can be derived. These ones correspond to:
formula
(29)
formula
(30)
formula
(31)
Consequently, mass and momentum fluxes, both in normal and transverse directions, can be computed with these characteristics, using the HLL or HLLC approach, as explained for the CM:
formula
(32)
formula
(33)
formula
(34)
where
formula
(35)
formula
(36)
Juez et al. (2014) used a Roe-type solver to calculate fluxes at the interfaces in their weakly coupled model. Here, the method is adapted to the lateralized HLLC solver by deriving a numerical characteristic that will ensure consistency between the integral of the exact solution (37) and the integral of the solution linearized locally at the interface (38):
formula
(37)
formula
(38)
such that
formula
(39)
formula
(40)
formula
(41)
with being the wave speed crossing normally the interface between two cells and guaranteeing the consistency between (37) and (38), its discretized form, the maximum between the local energy slopes of the two adjacent cells calculated with the Manning formula and ds the distance between the centres of the two cells.
In the same manner as with the hydrodynamic laws of conservation, the bed evolution is calculated with a first-order finite-volume method:
formula
(42)
where
formula
(43)
with and being the components of the unit solid discharges normal to the edge of the left and right cells, respectively.

In this way, the sediment flux is calculated in a very simple way and does not rest on a diffusive term that would eventually cause the bed to flatten in spite of the physics.

Boundary conditions

The boundary conditions are calculated directly as a flux at the boundary interfaces of the computational domain, without the use of any ghost cell. The expression for these boundary fluxes depends on the type of boundary condition. For three test cases considered hereafter, the same transmissive boundary condition was applied at the downstream end of the channel. Owing to its striking influence on the bed evolution, its mathematical description should not be left aside.

On the boundaries, no right cell exists. Formulae (18), (19), (22) and (23) can thus no longer be applied for the CM and neither can (32), (33), (34) and (43) for the WCM.

Except for the sediment mass flux, the CM and the WCM will share the same flux formulations. On one hand, for a supercritical flow, the description of the mass and momentum fluxes is rather straightforward. All the information comes from upstream and a full upwind flux scheme can be used:
formula
(44)
On the other hand, a subcritical flow raises some issues. The information cannot come from downstream as there is no right cell. Since the water level suddenly decreases, a M2-type water profile going through the critical depth is expected to develop. Regarding the mass and normal momentum equations, Savary & Zech (2007) suggested using the compatibility equations of the pure hydrodynamic system (26) to derive the boundary conditions:
formula
(45)
Discretizing equations (Equation (45)) and focusing only on the transmission of the information along the positive characteristic (31), for Equation (29) cannot be defined at the boundary, we get:
formula
(46)
These two equations contain three unknowns. If we impose the water depth at the boundary to be the critical depth of the left cell , we can easily derive , and :
formula
(47)
As mentioned before, the sediment mass flux always relies on both upstream and downstream information. For the CM, the Rankine–Hugoniot relation around is used:
formula
(48)
Assuming that at a free boundary, for the sediment can only be evacuated through the spilling section, the expression for the sediment flux becomes:
formula
(49)
For the WCM, the following equations are used:
formula
(50)
with the systems of Equations (49) and (50) being valid for the CM and the WCM, respectively.

RESULTS AND DISCUSSION

The two models were applied to four challenging test cases, compared with each other and with experimental data when available: a theoretical supercritical flow over a steep bed slope made of erodible material as well as dam-break flows over mobile beds with strong erosion and deposition, both in one and two dimensions. All simulations were run on a laptop housing eight Intel Core i5-8350 U 1.7 GHz CPU and 8 GB of RAM.

Equilibrium slopes

As a first application, our two 2D models were tested on a 1D inclined channel, which is 4 m long, 10 m wide and with elements with a maximum edge size of 0.1 m. They are to be compared with the tests carried out by Franzini & Soares-Frazão (2018). Upstream, both water and solid discharges were imposed whereas a transmissive boundary condition at critical depth was applied to the downstream end, so that only the left cell was considered. The bedrock level was set to zero on the whole channel so that only the movable bed could be responsible for the slope. The downstream bed level was also set to zero. Consequently, no erosion could occur there and the downstream end would behave as a rotation point. A water discharge by unit width and a solid discharge by unit width were imposed upstream. The expected flow was thus highly unidirectional. The bed was considered to be constituted by a uniform coarse sand with and a relative specific gravity . The Manning coefficient of the grains was set to and the porosity of the bed to . The corresponding equilibrium slope could thus be derived as . Initially, the channel was also already filled with the corresponding uniform water depth .

Three different cases with different initial slopes were then investigated. Besides the equilibrium slope as an initial value, milder () and steeper () slopes than were considered. Aggradation and degradation were then respectively expected.

Unlike the results obtained for this case by Franzini & Soares-Frazão (2018), the CM (Figure 3(a)) and the WCM (Figure 3(b)) show different equilibrium slopes. While the WCM reaches the correct equilibrium slope, the use of the CM leads to a milder one. Also, the same milder slope is to be found for the three different initial slopes. It is thus inherent to the CM, rather than related to the initial state. Even if the CM does not lead to a complete flat bed, it appears the HLLC flux scheme overestimates the equilibrium sediment flux. The usefulness of the WCM is thus hereby demonstrated for simple and steady applications.

Figure 3

Equilibrium slope reached with (a) the CM and (b) the WCM for a supercritical flow after aggradation (dashed line), degradation (dotted line) and with the expected equilibrium slope as an initial condition (dash-dotted line). The continuous line stands for the expected equilibrium slope.

Figure 3

Equilibrium slope reached with (a) the CM and (b) the WCM for a supercritical flow after aggradation (dashed line), degradation (dotted line) and with the expected equilibrium slope as an initial condition (dash-dotted line). The continuous line stands for the expected equilibrium slope.

1D dam-break over mobile bed

The second test case aims at assessing the limits of our models, by applying them to dam-break flows causing an important interaction between the flow and the bed layer. In order to do so, we confronted our models with the experimental dataset provided by Spinewine & Zech (2007) and carried out at the Hydraulics Unit of the LEMSC (Mechanical and Civil Engineering Laboratory, Université Catholique de Louvain, Belgium). These experiments consist of a 1D dam-break in a 6 m-long flume with sediment both upstream and downstream of a gate with an opening time of roughly 0.1 s. In the investigated configuration, the bed is perfectly flat and no water lies downstream of the gate while the upstream water level reaches above the bed. Two bed materials were used. A first dataset was constituted with coarse sand, characterised by , a relative specific gravity and a porosity . The Manning coefficient of the grains was measured to be . The other one used PVC pellets with , and . Since Spinewine & Zech (2007) did not mention any measure of the Manning coefficient for the PVC pellets, we used the Strickler's equation (Chaudhry 2008) so that .

The models used the bank-failure operator described in Swartenbroekx et al. (2010) and were compared with the experimental data at and for the sand case at and for the PVC case. Since Spinewine & Zech (2007) recorded different profiles for the upper () and lower () bed interfaces, the bed RMSE we will be referring to corresponds to the following formula:
formula
(51)
with being the bed level computed by the numerical model, and the experimental levels of the upper and lower bed interfaces, respectively, and k the number of experimental data points. A minimum RMSE corresponds to a bed level staying between both interfaces and showing a characteristic level of the whole bedload layer.

Regarding the sand test case (Figure 4), the CM and WCM behave similarly, as can be deduced from Table 1. The bed level computed by both models follows the same tendency as the two bed interfaces, i.e., a noticeable erosion around the gate along with a slight increase of the characteristic bed level downstream of it. However, it hardly shows any strong erosion or deposition. Both models also overpredict the arrival time of the wave, especially when . The two models show more agreement with the bed level than with the water level, for which the RMSE is the largest. As regards computational time, the CM lasted longer than the WCM. To summarize, even if the performance of the models is not completely satisfactory, their performance is very similar and we cannot really state the superiority of one model over the other.

Table 1

RMSE between the numerical and experimental results for the sand test case


RMSE (mm)0.250.50.751Mean
WCM  15.5 8.7 8.33 8.2 10.18 
 4.93 4.91 6.11 5.45 5.35 
CM  15.36 8.5 7.92 7.83 9.9 
 4.92 4.92 6.13 5.47 5.36 

RMSE (mm)0.250.50.751Mean
WCM  15.5 8.7 8.33 8.2 10.18 
 4.93 4.91 6.11 5.45 5.35 
CM  15.36 8.5 7.92 7.83 9.9 
 4.92 4.92 6.13 5.47 5.36 
Figure 4

Evolution of bed and water profiles with time with sand as bed material.

Figure 4

Evolution of bed and water profiles with time with sand as bed material.

The use of PVC as bed material enables much more solid transport than with sand, for it is much lighter, and pushes the models towards their limits. In that sense, Figure 5 clearly shows that the CM is not able to behave correctly in this case, especially regarding the water level during the latter stages. This contrasts with the WCM that does not show unphysical irregularities, despite its poor performance (Table 2). Finally, applying the CM to PVC as bed material lasted 41% longer than with the WCM. A look at the time step evolution during the simulation indeed shows that the CM is unstable in these conditions (Figure 6). Figuring out the reason for these instabilities goes beyond the scope of this paper, but might be worth future investigation. The approximation of the eigenvalues might be a lead to start with, as Cordier et al. (2011) suggested. For such an application, the CM must thus be rejected whereas the WCM could be used as a first estimation of the water and bed levels. Two-phase or two-layer models such as described, respectively, by Di Cristo et al. (2016) and Swartenbroekx et al. (2013), would certainly obtain much better performances, despite a higher computational cost. The Roe-based non-capacity bedload approach developed by Martínez-Arranda et al. (2019) could also lead to improved performance of the models.

Table 2

RMSE between the numerical and experimental results for the PVC test case


RMSE (mm)0.250.51.251.5Mean
WCM  14.88 14.2 18.94 18.5 16.63 
 27.55 13.16 16.62 16.5 18.46 
CM  12.31 14.1 27.05 30.43 20.97 
 27.06 13.04 16.57 16.39 18.27 

RMSE (mm)0.250.51.251.5Mean
WCM  14.88 14.2 18.94 18.5 16.63 
 27.55 13.16 16.62 16.5 18.46 
CM  12.31 14.1 27.05 30.43 20.97 
 27.06 13.04 16.57 16.39 18.27 
Figure 5

Evolution of bed and water profiles with time with PVC as bed material.

Figure 5

Evolution of bed and water profiles with time with PVC as bed material.

Figure 6

Time step evolution of the WCM (dotted line) and the CM (continuous line) with PVC as bed material.

Figure 6

Time step evolution of the WCM (dotted line) and the CM (continuous line) with PVC as bed material.

NSF-PIRE benchmark test

This second bunch of simulations is to be compared with the results obtained by Soares-Frazão et al. (2012) for which they also collected experimental measures, both in bed elevation and water level. Two test cases were run in a wide and long flume at the Hydraulics Unit of the LEMSC. The dimensions of this flume are given in Figure 7. For a complete description of the experimental conditions, please refer to Soares-Frazão et al. (2012). Initially, the upstream part of the flume (before the gate) was filled with water, unlike the rest of it. Only of the reservoir and 9.5 m of the downstream part were covered by an 8.5 cm sand layer. This sand was considered to be uniform and characterized by , a relative specific gravity and a porosity . The Manning coefficient was measured to be .

Figure 7

(a) Plane view of the flume. The shaded area stands for the zone covered by sand. (b) Position of the water level gauges.

Figure 7

(a) Plane view of the flume. The shaded area stands for the zone covered by sand. (b) Position of the water level gauges.

As the gate spans only over a limited part of the cross section, the resulting flow is two-dimensional. Two configurations were tested: one with a just saturated sediment layer in the downstream channel, and one with a water level of 0.15 m in the downstream part, resulting in a submerged sediment layer. Water levels were recorded at different locations downstream from the gate, and the final bed elevation was measured using laser profiling.

Configuration 1

The first test case aims at simulating a dam-break on a dry mobile bed. The measured final bed elevation is illustrated in Figure 8. An important scour was noticed just after the enlargement, while the maximum deposits took place in the first metres, along the edges of the flume. These deposits then formed a tongue-like shape until both deposition fronts merged. The focus is only put on the first 8 metres after the gate.

Figure 8

Configuration 1 measured topography after 20 s in metres.

Figure 8

Configuration 1 measured topography after 20 s in metres.

Numerical simulations were run on two meshes. The coarsest one had a spatial resolution (hereafter denoted by r) of while the finest one halved this value. In order to determine which mesh to use, we compared the results obtained by both models in terms of bed elevation. Regarding the WCM, the results of the coarse mesh (Figure 9(a))) are barely less accurate than those of the fine mesh (Figure 9(b)) on average. Nevertheless, the amplitudes of the maximum deposits (Figure 10(a)) and scours (Figure 10(b)) are better estimated on the fine mesh, sometimes in spite of a lag between the numerical and experimental results leading to a slightly larger RMSE for the fine mesh along the profile (Table 3). The largest difference comes with the CM. The coarse mesh (Figure 9(c)) indeed delivers really poor results in comparison with those collected with the fine one (Figure 9(d)). Even if the accuracy of the results delivered by the latter is still inferior to that of the WCM run on the coarse mesh, the increase of the resolution definitely constitutes a way to improve the results of the CM (Figure 10(c) and 10(d)). However, comparing the performances of the two models despite the mesh refinement, it appears the WCM (Figure 10(a) and 10(b)) is better than the CM (Figures 10(c) and10(d)) at reproducing extreme values of bed elevation, in particular when it comes to deposition along the edges. Nonetheless, it is interesting to note that the improvement of the results of the WCM with the resolution is maybe not worth the associated computational cost. While the CM and WCM have a similar computational cost for the same mesh resolution (the WCM lasts only 4% longer than the CM for both resolutions), using the coarse mesh cuts the simulation time by 7.8 times for both models. Hence, the accuracy gain of the WCM by refining the mesh is quite limited. Nevertheless, to maximize the chances of the CM hereafter, we will keep focusing on the results obtained with the fine mesh.

Table 3

RMSE between the numerical and experimental bed levels for configuration 1

WCM
CM
 5.5 6.5 13.2 8.6 
 10.2 10.5 10.3 9.4 
WCM
CM
 5.5 6.5 13.2 8.6 
 10.2 10.5 10.3 9.4 
Figure 9

Bed elevation (m) numerical results obtained for the WCM applied to the coarse mesh (a) and the fine mesh (b) and for the CM applied to the coarse mesh (c) and the fine mesh (d) for configuration 1.

Figure 9

Bed elevation (m) numerical results obtained for the WCM applied to the coarse mesh (a) and the fine mesh (b) and for the CM applied to the coarse mesh (c) and the fine mesh (d) for configuration 1.

Figure 10

Comparison of bed elevation results for the WCM along (a) and (b) and for the CM along (c) and (d). Numerical results obtained on a coarse mesh (dotted lines) and a finer one (dashed lines) are compared with experimental results (continuous lines).

Figure 10

Comparison of bed elevation results for the WCM along (a) and (b) and for the CM along (c) and (d). Numerical results obtained on a coarse mesh (dotted lines) and a finer one (dashed lines) are compared with experimental results (continuous lines).

When it comes to water level results (see Figure 7(b) for the position of the gauges), the WCM does not stand out from the CM as for bed elevation. Although numerical results are really close for the furthest gauges, as indicated by Figure 11(a), a significant difference can be noticed for gauges located closer to the gate. For these gauges, the WCM generally provides worse results than the CM (Figure 11(b)).

Figure 11

Configuration 1 water levels' comparison between experimental results (continuous line), the WCM (circles) and the CM (stars) at gauges no. 6.1 (a) and no. 1 (b).

Figure 11

Configuration 1 water levels' comparison between experimental results (continuous line), the WCM (circles) and the CM (stars) at gauges no. 6.1 (a) and no. 1 (b).

Configuration 2

In this second experiment, the water level downstream of the gate is no longer zero. The final bed elevation (Figure 12) shows that beside the tongue-like shape of the bed that we already had with the first configuration, we now have a second scour downstream of the first one and oscillatory bed forms.

Figure 12

Configuration 2 measured topography after 20 s in metres.

Figure 12

Configuration 2 measured topography after 20 s in metres.

Numerical simulations were run on the same fine mesh as for configuration 1. Once again, the comparison between Figure 13(a) and 13(b) shows that the WCM delivers results closer to the experimental ones, but to a very limited extent in this case. Both models indeed did not reproduce any bed form, nor the second scour. Moreover, the amplitude of the deposits and the first scour are less accurately calculated than in the first case. Although the WCM gives slightly better results than the CM, Figure 14 and Table 4 show that both models fail to deliver satisfactory bed elevation results in this case. Once again, the WCM lasted longer than the CM.

Table 4

RMSE between the numerical and experimental bed levels for configuration 2

WCMCM
 14.1 16.1 
 18.7 19.6 
WCMCM
 14.1 16.1 
 18.7 19.6 
Figure 13

Bed elevation (m) numerical results obtained with the WCM (a) and the CM (b) for configuration 2.

Figure 13

Bed elevation (m) numerical results obtained with the WCM (a) and the CM (b) for configuration 2.

Figure 14

Bed elevation comparison between the results obtained with the WCM (dotted line), the CM (dashed line) and the experimental measures (continuous line) along profile (a) and (b) for configuration 2.

Figure 14

Bed elevation comparison between the results obtained with the WCM (dotted line), the CM (dashed line) and the experimental measures (continuous line) along profile (a) and (b) for configuration 2.

With regard to water levels, the WCM and the CM generally give very close results (Figure 15(a)), even more for the downstream gauges (Figure 15(b)).

Figure 15

Water levels' comparison between experimental results (continuous line), the WCM (circles) and the CM (stars) for configuration 2 at gauge no. 1 (a) and no. 6.2 (b).

Figure 15

Water levels' comparison between experimental results (continuous line), the WCM (circles) and the CM (stars) for configuration 2 at gauge no. 1 (a) and no. 6.2 (b).

Channel with a 90° bend

Our last test case considered a channel with a 90° bend. In order to do so, we used preliminary data gathered by Soares-Frazão et al. (2019) on an experimental flume at the Hydraulics Laboratory of the iMMC. Figure 16 illustrates the dimensions of this flume.

Figure 16

Plane view of the experimental channel with a 90° bend. Dimensions in metres, after Soares-Frazão & Zech (2002).

Figure 16

Plane view of the experimental channel with a 90° bend. Dimensions in metres, after Soares-Frazão & Zech (2002).

The light shaded area of Figure 16 represents the reservoir, whose bed level is below the reference level and is closed by a gate that can be lifted rapidly, manually. For this experiment, the reservoir was filled with of water, i.e., cm above the reference level. Downstream of the gate lies a uniform layer of sand. The final time was considered after opening the gate, as no more significant bed evolution could be noticed. The final topography that Soares-Frazão et al. (2019) measured, thanks to a laser sheet technique, is given in Figure 17.

Figure 17

Final topography as measured by Soares-Frazão et al. (2019) with the focus put only on the second part of the channel that was rotated by 270° for convenience reasons.

Figure 17

Final topography as measured by Soares-Frazão et al. (2019) with the focus put only on the second part of the channel that was rotated by 270° for convenience reasons.

The sand used in this experiment was also considered uniform and could be characterized by , a relative specific gravity , a bed porosity and a Manning friction coefficient , while the friction of the glass or Plexiglas walls of the flumes was disregarded.

Comparison between experimental and numerical results

The simulations were run on an unstructured mesh and used the bank-failure operator described in Swartenbroekx et al. (2010). In the reservoir, the resolution was set to , while the flume itself was meshed with a resolution of (Figure 18). Running the CM on this mesh did not reproduce deposits and scours satisfactorily (Figure 19(a)). These zones were actually barely distinguishable. Although its results were still far away from the expectations, the WCM then highlighted the deposition and erosion maxima a little more (Figure 19(b)). It also led to another scour around the sharp edge of the bend, which was undetected by the CM. However, the poor spatial resolution of the preliminary data does not let us tell whether this scour actually occurs or not. Furthermore, the WCM shows a slight oscillatory behaviour, particularly in the beginning of the simulation, when high Froude numbers are to be dealt with. However, these oscillations look to be self-stabilizing. In terms of computational time, a 115 s simulation with the CM lasted longer than with the WCM. Here again, the difference remains limited.

Figure 18

Unstructured mesh of the 90° bend test case.

Figure 18

Unstructured mesh of the 90° bend test case.

Figure 19

Bed elevation (m) computed after with the CM (a) and the WCM (b).

Figure 19

Bed elevation (m) computed after with the CM (a) and the WCM (b).

Once again, we can use longitudinal profiles to compare both models on zones of interest. To assess the capacity of the models to reproduce the maximum deposit in the inner side of the meander and the maximum erosion on its outer side, we will refer to profiles and , respectively, as indicated by Figure 17. Figure 20(a) shows that the WCM better reproduces the maximum deposit than the CM, even if the error is still significant. Moreover, as mentioned earlier, it computes an important scour around the northwest corner of the bend that a look at the experimental results also confirms. However, the amplitude of the measured scour is not as important as given by the WCM. The accuracy of this measure is to be questioned, though. A local minimum could indeed have been smoothened by the resolution of the laser sheets grid. Regarding this profile, the RMSE between the WCM and the experimental data remains limited to 12.14 mm while that of the CM reaches 17.36 mm. Figure 20(b) then shows that the WCM also better estimates the scour on the outer side of the meander, even though it is far away from the measures too. Here, the RMSEs of the WCM and of the CM are closer to each other (11.5 mm and 13.9 mm, respectively). Regarding the erosion downstream of the bend, the WCM looks to perform much better along both profiles. Nevertheless, the WCM still fails to accurately reproduce deposits and scours caused by the meandering flow. We think this is mainly due to secondary currents that are not considered in our model. It would definitely be worth investigating this hypothesis. Despite the 3D behaviour of these processes, 2D methods, as developed by Kassem & Chaudhry (2002), for example, would let us consider their impact in our model.

Figure 20

Comparison between bed elevation results computed with the WCM (dotted line), the CM (dashed line) and the experimental results (continuous line) after 115 s along longitudinal profiles A (a) and B (b).

Figure 20

Comparison between bed elevation results computed with the WCM (dotted line), the CM (dashed line) and the experimental results (continuous line) after 115 s along longitudinal profiles A (a) and B (b).

CONCLUSION

We presented here two models aimed at simulating transient flows on mobile beds. Both were finite-volume models using a lateralized HLLC flux scheme and were run on unstructured grids. The difference between the two concerned the coupling between the shallow water equations and the Exner equation. On one hand, the coupled model (CM) takes the sediment-related information into account when calculating hydraulic fluxes. On the other hand, the weakly coupled model (WCM), developed after Juez et al. (2014), does not consider sediments at all in hydraulic fluxes. Nevertheless, its sediment mass flux is based on a numerical characteristic linked to the hydraulic state via a solid discharge closure formulation.

Both models were confronted with four test cases: an inclined channel, a 1D dam-break causing a large solid transport, a dam-break with a sudden enlargement and another dam-break, but on a channel with a bend. The first test case showed that the CM was not able to converge towards the right equilibrium slope, because of a diffusive term linked to the bed level difference in its sediment mass flux. The second test case emphasized the instability of the CM when subject to a very high interaction between the flow and the bed and showed its inadequacy for such applications. Nonetheless, the performance of the two models was rather poor when a bed material as light as PVC was considered. The third test case then showed that the WCM was better at reproducing large morphological changes than the CM. Nevertheless, both models were not able to reproduce bed forms in the second configuration and the WCM results were only just better than those of the CM. Regarding water levels, both models approximately equally performed. Finally, the last test case confirmed the ability of the WCM to better reproduce large deposits and scours than the CM. However, the error of both models with the measurements was significant. The integration of secondary currents into our model could certainly help to increase the accuracy of our results. As regards computational time, both models generally almost equally performed, except when the CM became unstable in the second test case.

Despite the issues highlighted by the different test cases, we showed that the WCM performs better than our fully CM to simulate transient flows on mobile beds, provided that the hypotheses do not differ from those used hereabove. Despite its oscillatory behaviour, the WCM also delivers better results than the CM when subjected to a steady flow. Eventually, we showed that the WCM also delivers accurate results with a lateralized HLLC flux scheme, and not only with a Roe-based one, as developed by Juez et al. (2014).

ACKNOWLEDGEMENTS

Besides the Université Catholique de Louvain, Belgium that supported this work through an FSR PhD grand for the first author, we would also like to thank Stéphanie Abbeels and Marie Messens, who did a great job at gathering and treating the data we used for our numerical simulations on the channel with the bend, and Yves Zech for the outside perspective and precious advice he gave us.

DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

REFERENCES

REFERENCES
Barzgaran
M.
Mahdizadeh
H.
Sharifi
S.
2019
Numerical simulation of bedload sediment transport with the ability to model wet/dry interfaces using an augmented Riemann solver
.
Journal of Hydroinformatics
21
(
5
),
834
850
.
doi:10.2166/hydro.2019.046
.
Camenen
B.
Larson
M.
2008
A general formula for noncohesive suspended sediment transport
.
Journal of Coastal Research
24
(
3
),
615
627
.
doi:10.2112/06-0694.1
.
Cao
Z.
Day
R.
Egashira
S.
2002
Coupled and decoupled numerical modeling of flow and morphological evolution in alluvial rivers
.
Journal of Hydraulic Engineering
128
(
3
),
306
321
.
Chaudhry
H.
2008
Open-Channel Flow
, 2nd edn.
Springer
,
New York
,
USA
.
Cordier
S.
Le
M. H.
Morales de Luna
T.
2011
Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help
.
Advances in Water Resources
34
,
980
989
.
doi:10.1016/j.advwatres.2011.05.002
.
Costanzo
C.
Macchione
F.
Viggiani
G.
2002
The influence of source terms treatment in computing two-dimensional flood propagation
. In:
Proceedings of River Flow 2002
(D. Bousmar & Y. Zech, eds), Vol.
1
.
Louvain-la-Neuve
,
Belgium
, pp.
277
282
.
Courant
R.
Friedrichs
K.
Lewy
H.
1967
On the partial difference equations of mathematical physics
.
IBM Journal
11
,
215
234
.
Di Cristo
C.
Greco
M.
Iervolino
M.
Leopardi
A.
Vacca
A.
2016
Two-Dimensional Two-Phase Depth-Integrated Model for Transients over Mobile Bed
.
Journal of Hydraulic Engineering
142
(
2
),
04015043
.
Di Cristo
C.
Greco
M.
Iervolino
M.
Vacca
A.
2019
Interaction of a dam-break wave with an obstacle over an erodible floodplain
.
Journal of Hydroinformatics
22
(
1
),
5
19
.
doi:10.2166/hydro.2019.014
.
Engelund
F.
Hansen
E.
1967
A Monograph on Sediment Transport in Alluvial Streams
.
Teknisk Forlag
,
Sklbrekgade 4, Copenhagen V, Denmark
.
Fraccarollo
L.
Capart
H.
Zech
Y.
2003
A Godunov method for the computation of erosional shallow water transients
.
International Journal for Numerical Methods in Fluids
41
(
9
),
951
976
.
doi:10.1002/fld.475
.
Franzini
F.
Soares-Frazão
S.
2018
Coupled finite-volume scheme with adapted Augmented Roe scheme for simulating morphological evolution in arbitrary cross-sections
.
Journal of Hydroinformatics
20
(
5
),
1111
1130
.
doi:10.2166/hydro.2018.109
.
Garegnani
G.
Rosatti
G.
Bonaventura
L.
2011
Free surface flows over mobile bed: mathematical analysis and numerical modelling of coupled and decoupled approaches
.
Communications in Applied and Industrial Mathematics
2
(
1
).
Geuzaine
C.
Remacle
J.-F.
2009
Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities
.
International Journal for Numerical Methods in Engineering
79
(
11
),
1309
1331
.
doi:10.1002/nme.2579
.
Goutière
L.
Soares-Frazão
S.
Savary
C.
Laraichi
T.
Zech
Y.
2008
One-dimensional model for transient flows involving bed-load sediment transport and changes in flow regimes
.
Journal of Hydraulic Engineering
134
(
6
),
726
735
.
doi:10.1061/(ASCE)0733-9429(2008)134:6(726)
.
Guinot
V.
2008
Wave Propagation in Fluids: Models and Numerical Techniques
.
ISTE
,
London
,
UK
.
Harten
A.
Lax
P. D.
Van Leer
B.
1983
On upstream differencing and Godunov-type schemes for hyperbolic conservation laws
.
SIAM Review
25
(
1
),
35
61
.
doi:10.1137/1025002
.
Hou
J.
Han
H.
Li
Z.
Guo
K.
Qin
Y.
2018
Effects of morphological change on fluvial floods patterns evaluated by a hydro-geomorphological model
.
Journal of Hydroinformatics
20
(
3
),
633
644
.
doi:10.2166/hydro.2018.142
.
Juez
C.
Murillo
J.
García-Navarro
P.
2013
Numerical assessment of bed-load discharge formulations for transient flow in 1D and 2D situations
.
Journal of Hydroinformatics
15
(
4
),
1234
1257
.
doi:10.2166/hydro.2013.153
.
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
.
doi:10.1016/j.advwatres.2014.05.014
.
Juez
C.
Soares-Frazão
S.
Murillo
J.
García-Navarro
P.
2017
Experimental and numerical simulation of bed load transport over steep slopes
.
Journal of Hydraulic Research
55
(
4
),
455
469
.
doi:10.1080/00221686.2017.1288417
.
Kassem
A. A.
Chaudhry
M. H.
1998
Comparison of coupled and semicoupled numerical models for alluvial channels
.
Journal of Hydraulic Engineering
124
(
8
),
794
802
.
Kassem
A. A.
Chaudhry
M. H.
2002
Numerical modeling of bed evolution in channel bends
.
Journal of Hydraulic Engineering
128
,
5
.
doi:10.1061/(ASCE)0733-9429(2002)128:5(507)
.
Lyn
D. A.
Altinakar
M.
2002
St. Venant-Exner equations for near-critical and transcritical flows
.
Journal of Hydraulic Engineering
128
(
6
),
579
587
.
doi:10.1061/(ASCE)0733-9429(2002)128:6(579)
.
Martínez-Arranda
S.
Murillo
J.
García-Navarro
P.
2019
A comparative analysis of capacity and non-capacity formulations for the simulation of unsteady flows over finite-depth erodible beds
.
Advances in Water Resources
130
,
91
112
.
doi:10.1016/j.advwatres.2019.06.001
.
Meyer-Peter
E.
Müller
R.
1948
Formulas for Bed-Load Transport
. In:
Proceedings of the 2nd Meeting of the IAHF
,
Stockholm, Sweden
.
Murillo
J.
García-Navarro
P.
2010
Weak solutions for partial differential equations with source terms: application to the shallow water equations
.
Journal of Computational Physics
229
,
4327
4398
.
doi:10.1016/j.jcp.2010.02.016
.
Roushangar
K.
Hassanzadeh
Y.
Ali Keynejad
M.
Tagi Alami
M.
Nourani
V.
Mouaze
D.
2011
Studying of flow model and bed load transport in a coarse bed river: case study – Aland River, Iran
.
Journal of Hydroinformatics
13
(
4
),
850
866
.
doi:10.2166/hydro.2010.010
.
Savary
C.
Zech
Y.
2007
Boundary conditions in a two-layer geomorphological model. Application to a hydraulic jump over a mobile bed
.
Journal of Hydraulic Research
45
(
3
),
316
332
.
doi:10.1080/00221686.2007.9521766
.
Soares-Frazão
S.
Zech
Y.
2002
Dam break in channels with 90° bend
.
Journal of Hydraulic Engineering
128
(
11
),
956
968
.
Soares-Frazão
S.
Zech
Y.
2010
HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed
.
International Journal for Numerical Methods in Fluids
66
,
1019
1036
.
doi:10.1002/fld.2300
.
Soares-Frazão
S.
Canelas
R.
Cao
Z.
Cea
L.
Chaudhry
H. M.
Die Moran
A.
El Kadi
K.
Ferreira
R.
Fraga Cadórniga
I.
Gonzalez-Ramirez
N.
Greco
M.
Huang
W.
Imran
J.
Le Coz
J.
Marsooli
R.
Paquier
A.
Pender
G.
Pontillo
M.
Puertas
J.
Spinewine
B.
Swartenbroekx
C.
Tsubaki
R.
Villaret
C.
Wu
W.
Yue
Z.
Zech
Y.
2012
Dam-break flows over mobile beds: experiments and benchmark tests for numerical models
.
Journal of Hydraulic Research
50
(
4
),
364
375
.
doi:10.1080/00221686.2012.689682
.
Soares-Frazão
S.
Abbeels
S.
Messens
M.
2019
Dam-break flow in a channel with a 90° bend and mobile bed: experimental and numerical simulations
. In:
38th IAHR World Congress
,
Panama City
.
Spinewine
B.
Zech
Y.
2007
Small-scale laboratory dam-break waves on movable beds
.
Journal of Hydraulic Research
45
(
Suppl. 1
),
73
86
.
doi:10.1080/00221686.2007.9521834
.
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
(
3
),
302
314
.
doi:10.1080/00221686.2010.481856
.
Swartenbroekx
C.
Zech
Y.
Soares-Frazão
S.
2013
Two-dimensional two-layer shallow water model for dam-break flows with significant bed load transport
.
Internation Journal for Numerical Methods in Fluids
73
(
5
),
477
508
.
doi:10.1002/fld.3809
.
Toro
E. F.
1999
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
.
Springer
,
Berlin
,
Germany
.
Toro
E. F.
Spruce
M.
Speares
W.
1994
Restoration of the contact surface in the HLL-Riemann solver
.
Shock Waves
4
,
25
34
.
doi:10.1007/BF01414629
.
Van Emelen
S.
Zech
Y.
Soares-Frazao
S.
2015
Impact of sediment transport formulations on breaching modelling
.
Journal of Hydraulic Research
53
(
1
),
60
72
.
doi:10.1080/00221686.2014.939111
.
Wong
M.
Parker
G.
2006
Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database
.
Journal of Hydraulic Engineering
132
(
11
).
doi:10.1061/(ASCE)0733-9429(2006)132:11(1159)
.
Wu
W.
2004
Depth-averaged two-dimensional numerical modeling of unsteady flow and nonuniform sediment transport in open channels
.
Journal of Hydraulic Engineering
130
(
10
),
1013
1024
.
doi:10.1061/(ASCE)0733-9429(2004)130:10(1013)
.
Zech
Y.
Soares-Frazão
S.
Spinewine
B.
Savary
C.
Goutière
L.
2009
Inertia effects in bed-load transport models
.
Canadian Journal of Civil Engineering
36
(
10
),
1587
1597
.
doi:10.1139/L09-052
.