A 2D HLL-based weakly coupled model for transient flows on mobile beds

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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

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.

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.

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.

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

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.

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.

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

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.

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.

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.

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.

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.

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

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 cannot be made publicly available; readers should contact the corresponding author for details.