Next Article in Journal
Comparison of Empirical and Analytical Solutions for Open-Channel Flow Velocity with Common Grass Species in Taiwan
Next Article in Special Issue
Silting in the Grand Canal in the Domain of Chantilly (Oise, France)—Catchment-Scale Hydrogeomorphological Reconnaissance and Local-Scale Hydro-Sedimentary Transport Modelling
Previous Article in Journal
Isotope Signs (234U/238U, 2H, 18O) of Groundwater: An Investigation of the Existence of Paleo-Permafrost in European Russia (Pre-Volga Region)
Previous Article in Special Issue
Turbulent Flow Structure in a Confluence: Influence of Tributaries Width and Discharge Ratios
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Comparative Analysis of HLLC- and Roe-Based Models for the Simulation of a Dam-Break Flow in an Erodible Channel with a 90 Bend

1
Fluid Mechanics-I3A, University of Zaragoza, 50009 Zaragoza, Spain
2
Insitute of Mechanics, Materials and Civil Engineering, Université Catholique de Louvain, 1348 Louvain-la-Neuve, Belgium
*
Author to whom correspondence should be addressed.
Water 2021, 13(13), 1840; https://doi.org/10.3390/w13131840
Submission received: 28 May 2021 / Revised: 25 June 2021 / Accepted: 28 June 2021 / Published: 1 July 2021
(This article belongs to the Special Issue Fluvial Hydraulics and Applications)

Abstract

:
In geophysical surface flows, the sediment particles can be transported under capacity (equilibrium) conditions or noncapacity (nonequilibrium) conditions. On the one hand, the equilibrium approach for the bedload transport assumes that the actual transport rate instantaneously adapts to the local flow features. The resulting system of equations, composed of the shallow water equations for the flow (SWE) and the Exner equation for the bed evolution, has been widely used to simulate bedload processes. These capacity SWE + Exner models are highly dependent on the setup parameters, so that the calibration procedure often disguises the advantages and flaws of the numerical method. On the other hand, noncapacity approaches account for the temporal and spatial delay of the actual sediment transport rate with respect to the capacity of the flow. The importance of assuming nonequilibrium conditions in bedload numerical models remains uncertain however. In this work, we compared the performances of three different strategies for the resolution of the SWE + Exner system under capacity and noncapacity conditions to approximate a set of experimental data with fixed setup parameters. The results indicate that the discrete strategy used to compute the intercell fluxes significantly affected the solution. Furthermore, the noncapacity approach can improve the model prediction in regions with complex transient processes, but it requires a careful calibration of the nonequilibrium parameters.

1. Introduction

Bedload sediment transport is an important process in natural water bodies such as rivers, reservoirs, and estuaries. Bedload refers to the solid particles being transported within a more or less thick layer near the bed surface. The sediment particles can be transported under capacity (equilibrium) conditions or noncapacity (nonequilibrium) conditions. The equilibrium approach assumes that the actual sediment transport rates are equal to the capacity of the flow to carry solid weight, which is only determined by instantaneous local flow features and can be formulated by different empirical closure relations found in the literature [1]. Most of the numerical models proposed for bedload transport are based on this capacity approach [2,3,4,5,6], which leads to the system composed by the shallow water equations (SWE) for the flow and associated with the Exner equation [7] to estimate the bed evolution. Multiple numerical schemes have been proposed for the resolution of the SWE + Exner system, involving different levels of coupling between the hydrodynamic and morphodynamic components of the system [8]. Among them, totally decoupled [9,10,11], weakly coupled [3,6,12], and fully coupled strategies [4,5,13,14] have been proposed and used to predict laboratory experiments and real-scale field cases. Nevertheless, qualitative and quantitative comparisons of different numerical approaches for the resolution of the same system of equations are rare, even though they could provide helpful information for developers of efficient simulation tools, despite consisting of an important joint effort [15].
These kinds of capacity models for the SWE+Exner system are highly dependent on the model parameters, such as the bed interface roughness, the bulk porosity of the bed layer, the critical shear stress at the bed interface for the initiation of the sediment movement, or the closure relation for the capacity transport rate. Therefore, most of the time, when these models are used to predict experimental results or real-scale field data, these factors are used as tuning parameters and calibrated to approximate the observations [16]. This is a widely accepted procedure, provided that the calibration values remain in the physical range of validity, since these parameters include complex physical processes that must be modeled. However, this calibration procedure blurs the advantages and flaws of the numerical method.
On the other hand, in noncapacity bedload models, the actual transport rates are computed through advection and mass exchange with the static erodible bed [11]. Natural morphodynamical systems are always changing in time and space, and hence, absolute equilibrium states rarely exist in natural conditions. Therefore, intuitively, noncapacity approaches appear to be more suitable than models based on the equilibrium assumption since they account for the temporal and spatial delay of the actual sediment transport rate with respect to its potential capacity. Nonetheless, unlike suspended load models, for which it has been demonstrated that the nonequilibrium assumption helps to properly estimate the solid suspended concentration and the bed evolution [17], the way to cope with inertia effects in the case of bedload transport remains uncertain [18].
The assumption of nonequilibrium conditions in a bedload mathematical model requires computing the actual transport rate, as well as the net exchange flux between static and moving bed layers [19]. The physical interaction between flow and sediment particles at the static–moving bed layer interface was studied at the grain scale by [20,21], who proposed a reformulation of the SWE + Exner system based on the thickness of the bedload transport layer. Following this basis, the work in [22] proposed a 1D generalized model for bedload transport able to evolve from noncapacity to capacity states. This model helps to improve the results in some widespread benchmark tests [2,9,23] without the need to recalibrate the setup parameters.
In this work, we compared the performances of three strategies that have proven efficient for the resolution of the SWE + Exner system under capacity conditions, with application to an experimental transient case [24] (submitted) and imposing a fixed set of model parameters. The selected strategies involve the widespread HLLC [6] and Roe [14] solvers for the intercell numerical flux, as well as coupled [15,25] and weakly coupled [3,6] resolution procedures for the hydrodynamic and morphodynamic components of the system. Furthermore, the experiment was also simulated using the generalized bedload transport model proposed by Martínez-Aranda et al. [22] in order to analyze the need for including transient nonequilibrium states in the final bed elevation. The paper is structured as follows: the governing equations are described in Section 2; then, Section 3 is devoted to the numerical methods, and the experimental dataset is also briefly introduced; in Section 4, the numerical results are reported and discussed; finally, the main conclusions are drawn in Section 5.

2. Governing Equations

Shallow-water models for surface flows are derived from the depth integration of the Navier–Stokes equations along the flow column. These shallow-type models can be divided into nonhydrostatic and hydrostatic, which assume a hydrostatic vertical distribution of the pressure along the flow layer and are more widespread for environmental flows. Using this hydrostatic hypothesis for the flow movement and neglecting the presence of solid particles in the water, the 2D mass and momentum conservation equations for a clear-water shallow flow can be written as:
h t + x ( h u ) + y ( h v ) = 0
( h u ) t + x ( h u 2 + 1 2 g h 2 ) + y ( h u v ) = g h z b x τ b x ρ w
( h v ) t + x ( h u v ) + y ( h v 2 + 1 2 g h 2 ) = g h z b y τ b y ρ w
h being the vertical flow depth, ( u , v ) the components of the depth-averaged flow velocity vector u along the global ( x , y ) horizontal coordinates, respectively, z b the bed layer elevation, g the gravitational acceleration, ( τ b x , τ b y ) the components of the depth-averaged shear resistance at the bed interface τ b for the flow along the global ( x , y ) horizontal coordinates, and ρ w the density of water.
Furthermore, the mass conservation equation for the bedload transport layer is expressed as:
( 1 ξ ) z b t + x ( q b x ) + y ( q b y ) = 0
where ξ denotes the bulk porosity of the bed layer and ( q b x , q b y ) are the components of the bulk bedload rate q b along the global ( x , y ) horizontal coordinates, respectively.
The equations forming the two-dimensional system (1) and (2) can be rewritten in vector form as:
U t + · E ( U ) = S b ( U ) + S τ ( U )
where U is the vector of the conserved variables:
U = h , h u , h v , z b T = h , q x , q y , z b T
and E ( U ) = F ( U ) , G ( U ) are the convective flux vectors along the ( x , y ) horizontal coordinates, respectively, expressed as:
F ( U ) = h u h u 2 + 1 2 g h 2 h u v 1 1 ξ q b x G ( U ) = h v h u v h v 2 + 1 2 g h 2 1 1 ξ q b y
The vector S b ( U ) accounts for the momentum source term associated with the variation of the pressure force on the bed surface:
S b ( U ) = 0 g h z b x g h z b y 0
The source vector S τ ( U ) denotes the momentum dissipation due to the basal resistance:
S τ ( U ) = 0 τ b x / ρ w τ b y / ρ w 0
where the components of the shear resistance at the bed surface τ b = ( τ b x , τ b y ) are commonly estimated using the quadratic pure turbulent relation as:
τ b x = ρ w g h C f | u | u τ b y = ρ w g h C f | u | v
with C f = n b 2 h 4 / 3 being a friction coefficient and n b the bulk Manning’s roughness parameter.
The system of Equation (3), integrated over each computational cell and discretized in time and space, constitutes the basis of all the two-dimensional Finite Volume (FV) numerical models described hereafter for bedload transport. Nevertheless, the coupling between the hydrodynamical and morphodynamical components of the system, the flux resolution method at the edges between cells, or the solid discharge closure formulation might differ from one model to another.

Capacity vs. Noncapacity Approach for Bedload Transport

The bedload transport rate q b = ( q b x , q b y ) accounts for the volumetric solid discharge integrated in the bedload transport layer and needs a closure model for estimation. The solid particles can be transported under equilibrium or nonequilibrium conditions [17]. The equilibrium approach assumes that the actual sediment transport rate is equal to the capacity of the flow to carry solid weight. This equilibrium bedload rate, referred to as q b * , is only determined by instantaneous local flow features and can be formulated by different empirical closure relations found in the literature [1]. Most of these relations can be generally written as:
| q b * | = c θ m 1 ( Δ θ ) m 2 ( r s 1 ) g d s 3
where | q b * | denotes the bedload rate modulus under equilibrium conditions, c is a constant dimensionless coefficient, m 1 and m 2 are two constant exponents, r s = ρ s / ρ w is the solid/liquid density ratio with ρ s being the sediment particles density, and d s is the median grain diameter. The values of c, m 1 , m 2 , and θ c depend on the closure relation and are summarized in Table 1 for some widespread bedload transport rate formulae.
The term Δ θ denotes the positive excess of the Shields stress θ over the critical Shields stress value for the incipient motion θ c , expressed as:
Δ θ = θ θ c if θ > θ c 0 Otherwise
For the bedload transport process, the boundary shear stress at the bed surface z b is assumed fully turbulent and is modeled using the Manning’s resistance Equation (8). Hence, the dimensionless boundary Shields stress at the bed surface θ reads:
θ = n b 2 | u | 2 ( r s 1 ) d s h 1 / 3
Replacing Equation (11) into (9) for all the formulations in Table 1 allows demonstrating that | q b * | h 1 / 2 | u | 3 , and hence, a general formulation for the capacity solid transport rate based on the Grass law expressed as:
| q b * | = G ( h , θ ) | u | 3
has been adopted by other authors [3,5,25]. In Equation (12), the equilibrium bedload sediment discharge is related to the depth-averaged flow velocity by means of a factor G ( h , θ ) [ T 2 L 1 ] , which represents the interaction between the flow and the bed layer and which depends on the flow characteristics, contrary to the constant value initially proposed by [30].
On the other hand, in noncapacity transport, the actual bedload rate is computed through advection and mass exchange with the static erodible bed. Natural morphodynamical systems such as alluvial rivers are always changing in time and space, and hence, absolute equilibrium states rarely exist in natural conditions. Therefore, noncapacity approaches seem more suitable than models based on the equilibrium assumption since they account for the temporal and spatial delay of the actual sediment transport rate with respect to its potential capacity. However, if the adaptation delay is sufficiently small, equilibrium models can also be applied, at least in theory [19,31].
A generalized model for 1D noncapacity bedload transport was recently proposed by [22] and extended to the two-dimensional framework by [32]. The generalized bedload transport rate is expressed following a modified Grass-type law as:
q b x = G ( h , θ , η ) | u | 2 u q b y = G ( h , θ , η ) | u | 2 v
with G ( h , θ , η ) being an interaction factor between the flow and the bed. Therefore, G is calculated as a function of the flow depth h, the dimensionless Shields stress θ , and the bedload transport layer thickness η :
G = Γ 1 ( h ) Γ 2 ( θ ) Γ 3 ( η )
Expressions for functions Γ 1 , Γ 2 , and Γ 3 are reported in Table 2 for different empirical relations found in the literature, where k D and k E are the two positive calibration constants associated with the detention and entrainment exchange rates, respectively, between the bedload transport (moving) layer and the underlying static bed stratum.
The generalized Grass-type model of Equation (14) for the bedload transport has two important advantages compared with classical capacity models. First, it allows considering the capacity or noncapacity hypothesis for the bedload solid transport only by changing the expression to calculate the moving layer thickness η . Assuming the capacity approach leads to an instantaneous adaptation of the bedload transport rate to the flow carrying capacity, the equilibrium transport layer thickness η * can be estimated as:
η * = k E r s k D Δ θ d s
whereas the noncapacity assumption requires estimating the temporal evolution of the bedload transport layer thickness:
( 1 ξ ) η t + q b x x + q b y y = ( 1 ξ ) ( η ˙ D η ˙ E )
where η ˙ D and η ˙ E are the detention and entrainment exchange rates, respectively, between the bedload transport layer and the underlying static stratum, which can be estimated as [20]:
η ˙ D = k D η d s ( r s 1 ) g d s 3 d s η ˙ E = k E Δ θ r s ( r s 1 ) g d s 3 d s
Therefore, when equilibrium conditions are reached, η = η * , leading to η ˙ D = η ˙ E , and hence, Equation (16) reduces to the widespread Exner equation [7].

3. Numerical Models and Simulations

The following section contains the description of the different numerical models considered in this paper. All models have already been presented with more details in other publications, so only the most important features are presented here. After that, the test case used to compare them is described.

3.1. Roe-Based Solvers: R-Cap and R-NCap

Here, we briefly introduce the fully coupled model proposed by [14], based on the augmented Roe approach and using the generalized bedload transport formulation (13). The updating formulation for the conserved variables at the i-th cell U i between times t = t n and t = t n + 1 is expressed as:
U i n + 1 = U i n Δ t A i k = 1 NE R k 1 F k l k
where A i is the i-th cell area, Δ t = t n + 1 t n is the time step, NE is the number of edges of the i cell, R k 1 is the inverse of the rotation matrix R k at the k-th cell edge [33], l k is the edge length, and F k is the numerical normal flux to the k-th cell edge.
The upwind computation of F k is based on the eigenstructure of the pseudo-Jacobian matrix of the coupled system M ˜ k = ( J ˜ H ˜ ) k , at the cell edges. The term J ˜ k denotes the approximated Jacobian matrix of the convective flux J ˜ k = ( E · n ^ ) k / U , whereas H ˜ k is the nonconservative matrix of the bed-pressure source term (6) at the k-th cell edge, S b = H ˜ k U / x ^ n , being n ^ the outward normal unit vector and ( x ^ n , y ^ t ) the normal and transverse axes to the cell edge, respectively [14].
Using the generalized bedload transport Equation (13), the fully coupled pseudo-Jacobian matrix M ˜ k can be expressed as:
M ˜ k = 0 1 0 0 g h ˜ u ˜ n 2 2 u ˜ n 0 g h ˜ u ˜ n v ˜ t v ˜ t u ˜ n 0 ( u ˜ n a ˜ + v ˜ t b ˜ ) a ˜ b ˜ 0 k
where h ˜ , u ˜ n , and v ˜ t are the edge-averaged depth, normal velocity, and tangential velocity to the edge, respectively, whereas a ˜ and b ˜ represent the edge-averaged derivatives of the normal bedload flux 1 1 ξ G | u | 2 u · n ^ with respect to the normal ( q n ) and tangential ( q t ) flow discharges, respectively.
The approximate matrix M ˜ k of Equation (19) is diagonalizable with four approximate real eigenvalues, λ ˜ m , k for m = { 1 , , 4 } , resulting from the resolution of the characteristic polynomial:
P λ = ( u ˜ n λ ˜ ) λ ˜ [ ( u ˜ n λ ˜ ) 2 g h ˜ ] + g h ˜ a ˜ ( λ ˜ u ˜ n ) = 0
The associated orthogonal basis of eigenvectors ( e ˜ m ) k of M ˜ k is used to build the matrix P ˜ k = ( e ˜ 1 , e ˜ 2 , e ˜ 3 , e ˜ 4 ) k , which satisfies:
M ˜ k = ( P ˜ Λ ˜ P ˜ 1 ) k Λ ˜ k = λ ˜ 1 0 0 λ ˜ 4 k
Following Murillo and Navas-Montilla [34], the conserved variable increment δ U ^ k and the integrated resistance source terms ( S ˜ τ ) k (7) at the cell edge are projected on the eigenvector basis in order to obtain the wave and source strength vectors, A ˜ k and B ˜ k , respectively, leading to:
A ˜ k = ( α ˜ 1 , , α ˜ 4 ) k T = P ˜ k 1 δ U ^ k B ˜ k = ( β ˜ 1 , , β ˜ 4 ) k T = P ˜ k 1 ( S ˜ τ ) k
The approximate upwind solution F k for the numerical flux at the k-th cell edge is defined as:
F k = F ( R k U i ) + m ( λ ˜ m γ ˜ m e ˜ m ) k
where γ ˜ m = α ˜ m β ˜ m / λ ˜ m and the subscript m under the sum indicates waves traveling inward the i-th cell.
It is worth mentioning that this approach allows a full coupling of the morphological bed evolution with the hydrodynamical flow, leading to a robust scheme able to deal with capacity and noncapacity formulations for the bedload transport. The updating procedure for transport layer thickness η depends on the assumption made for the bedload transport:
  • On the one hand, assuming the capacity hypothesis, the cell-centered value of the transport thickness at the next time step t n + 1 is directly computed using:
    η i n + 1 = k E r s k D ( Δ θ ) i n + 1 d s
    where ( Δ θ ) i n + 1 is the nondimensional Shields excess (10) at the i-th cell computed with the conserved variables updated to time t n + 1 . This model is referred to as R-Cap;
  • On the other hand, the assumption of the noncapacity approach requires solving Equation (16) at each time step using the updating formula:
    η i n + 1 = η i n Δ t A i k = 1 N E F k η l k Δ t ( η ˙ D η ˙ E ) i n
    with F k η being the numerical flux at the k-th intercell edge for the homogeneous transport equation and ( η ˙ D η ˙ E ) i n denoting the cell-centered balance between the detention and entrainment rates (17) at time t n .
    To compute the numerical flux ( F p η ) k for the p-th sediment class at the k-th cell edge, the left-hand side of the transport Equation (16) is integrated along the normal direction x ^ n to the edge, and the numerical flux at the intercell interface is approximated using the linearized homogeneous Riemann problem [32]:
    η t + λ ˜ η , k η x ^ n = 0 η ( x ^ , 0 ) = η i n if x ^ n < 0 η j n if x ^ n > 0
    where the virtual bedload wave celerity λ ˜ η , k is defined as:
    λ ˜ η , k = δ G | u | 2 u · n ^ ( 1 ξ ) δ η k
    Therefore, the intercell numerical flux for the transport thickness is computed as:
    F k η = 1 1 ξ G | u | 2 u · n ^ i n if λ ˜ η , k > 0 1 1 ξ G | u | 2 u · n ^ j n if λ ˜ η , k < 0
Finally, regardless of which approach is selected, the cell-centered value of the bed-flow Grass interaction factor G at the next time step t n + 1 is computed using Equation (14) as:
G i n + 1 = Γ 1 ( h i n + 1 ) Γ 2 ( θ i n + 1 ) Γ 3 ( η i n + 1 )
with Γ 1 , Γ 2 , and Γ 3 as in Table 2.

3.2. HLLC-Based Solvers

For more information about the two models presented in this section, the reader may refer to [6], which described the models in detail and discussed their performance against four different test cases. The most important features of these models are presented here.

3.2.1. HLLC-CM

Equations (1a) to (11) still remain valid for this model, but the friction source terms S τ are not included in the intercell flux formulation, so that Equation (18) may be rewritten as:
U i n + 1 = U i n Δ t A i k = 1 NE R k 1 F ^ k l k + S τ , i Δ t
with F ^ k being the flux vector at the edge k that encompasses the bed-pressure term S b through a lateralization step as proposed by [35].
No transport layer is considered here, and a capacity approach ( q b = q b * ) using Equation (22) and the MPM coefficients (Table 1) is used. Besides this different formulation of the bedload rate, the HLLC-based Coupled Model (HLLC-CM) uses another flux solver at the interface. Except for the transverse momentum equation, the fluxes at the interfaces associated with the variables of Equation (4) use the HLL formula [36]:
F ^ k = S + F L S F R + S + S ( U R U L ) S + S
where F ^ k stands for the flux at the k-th edge between two cells represented by the subscripts R and L, omitting the lateralization step, which concerns the normal momentum equation only. S + and S then respectively correspond to the most positive and negative wave speed estimates. For the mass and normal momentum equations, they are chosen as follows:
S + = min ( λ L 1 , λ R 1 , 0 ) S = max ( λ L 4 , λ R 4 , 0 )
whereas for the bed level evolution equation, S + is expressed as:
S + = max ( λ L 2 , λ R 2 , 0 )
S keeps the same formulation however. λ used in Equations (32) and (33) are the eigenvalues of the pseudo-Jacobian matrix M derived for the system of Equation (3) in which the bed-pressure terms were incorporated. This matrix is very similar to Equation (19), but is cell-centered instead of being edge-averaged, despite being expressed in terms of the local system of coordinates ( x ^ n , y ^ t ) instead of the global one ( x , y ) , thanks to the rotational invariance property [37]:
M = 0 1 0 0 g h u n 2 2 u n 0 g h u n v t v t u n 0 1 1 ξ q b , n h 1 1 ξ q b , n q n 1 1 ξ q b , n q t 0
The three nonzero terms of the last row of Equation (34) depend on the considered bedload formulation, but only q b , n h will actually be present in the characteristic polynomial, hence being important for the eigenvalues. Since the MPM formula was selected for the HLLC-CM, this term equals:
q b , n h = 12 g ( r s 1 ) d s 3 n b 2 q n q n 2 + q t 2 1 2 ( r s 1 ) d s h 7 3 θ c 1 2 7 3 n b 2 q n q n 2 + q t 2 1 / 2 ( r s 1 ) d s h 10 / 3
The eigenvalues of M can then be approximated as proposed by [13]:
λ 1 = 1 2 ( u n c ( u n c ) 2 4 g h 1 ξ q b , n h 1 u n + c ) λ 2 = 1 2 ( u n c + ( u n c ) 2 4 g h 1 ξ q b , n h 1 u n + c ) λ 3 = u n λ 4 = u n + c
with c = g h .
To remedy the overly diffusive behavior of the HLL scheme, Reference [38] suggested incorporating the contact discontinuity into the wave pattern for the transverse momentum equation, which is equivalent to:
μ ^ n t = h u n v t ^ = v t , L q ^ n if u ˜ n 0 v t , R q ^ n if u ˜ n < 0
The last step of this flux scheme consists of the lateralization of S b . For all the equations, F ^ k = F ^ k , with the exception of the normal momentum equation, for which:
σ ^ k , n = σ ^ k , n g h ˜ ( z b , R z b , L ) S i S R S L
where σ n = h u n 2 + g h 2 2 is the normal momentum flux and S i = S L if the considered cell is left ofthe edge. Otherwise, S i = S R .
Altogether, Equations (31), (37), and (38) constitute the lateralized HLLC flux scheme used hereafter for the HLLC-CM.

3.2.2. HLLC-WCM

In contrast with the HLLC-CM, the HLLC-based Weakly Coupled Model (HLLC-WCM), relies on the hypothesis that the transport of sediment is low enough to overlook its influence on the system eigenvalues. The Exner equation is hence decoupled from the SWE. Consequently, the lateralization of the slope source terms no longer influences the Jacobian matrix. The well-known eigenvalues of the Jacobian matrix of the purely hydraulic system of Equations (1a)–(1c) can now be used for the HLLC scheme:
λ 1 = u n c λ 2 = u n λ 3 = u n + c
Regarding the update of z b , the same finite-volume approach is followed, but a new estimation of q ^ b , k , n must be provided:
z b n + 1 = z b n Δ t A i k = 1 NE 1 1 ξ q ^ b , n l k
q ^ b , k , n = q b , n , L if λ ^ b 0 q b , n , R if λ ^ b < 0
where λ ^ b is computed as proposed by [3] and is quite similar to Equation (27):
λ ^ b = 1 1 ξ Δ q n , b Δ z b Δ z b = Δ z b if Δ z b > d s S f Δ s if Δ z b d s
with Δ s being the distance between the centers of the two adjacent cells.
The flux solver associated with z b is hence not based on the HLLC scheme, but on an upwind approach where the sign of the virtual eigenvalue λ ^ b will determine the cell whose q b , n should correspond to the q ^ b , k , n at the k-th interface. For the hydrodynamic equations, nothing changes from the lateralized HLLC-CM, except for the eigenvalues.
As [6] highlighted, the choice of this solver is not only meant to simplify the implementation of the governing equations, but also to avoid the excessively diffusive behavior of the HLLC scheme applied to the bed elevation that [39] also tried to avoid with a slightly different approach. Moreover, the HLLC-WCM looks to be more stable than the HLLC-CM in the case of the high transport of the sediment [6].

3.3. Test Case

The test case selected to compare the different models presented above consisted of a two-dimensional dam-break occurring in a flume with a 90 bend and with an initial finite sandy layer of 0.075 m. Meurice et al. [24] reproduced this experiment in a laboratory and provided both water- and bed level results with different data-retrieval techniques [24]. Here, these results were used to assess the capabilities and limits of the different models.
The geometry of the flume and the position of the gauges used to record the evolution of the water level are illustrated in Figure 1. The water level of the reservoir was initially maintained 0.26 m above the nonerodible bed level of the channel (reference level). An aluminum gate separated the reservoir from the inlet channel and could be raised manually to simulate a quasi-instantaneous dam-break on an erodible bed. After 3.92 m, the water front would enter the second part of the flume that was perpendicular to the first one. This would lead the flow to show 2D and even 3D features. At the downstream end of this second part, the water flowed freely in a dissipation reservoir.
Water levels were recorded continuously at five different points with ultrasonic gauges, and no less than thirty-four different runs were performed. Hereafter, the numerical results were compared with the aggregated experimental results, rather than with one test only.
Bed levels were recorded with two different techniques. The temporal evolutions of different bed level cross-sections were recorded through laser profilometry, while photogrammetry was used to capture the final topography after the drainage of the channel. Hereafter, only the photogrammetric results were used to compare the experiments with the simulations because of their total spatial coverage of the channel. This last technique was used to reconstruct the topography consecutive to three different experimental runs. The average topography was then calculated and is represented in Figure 2. This was the dataset that was compared with the numerical results as far as the bed level surface was concerned.
On the one hand, large accumulation regions appeared downstream of the inner corner and in the outer corner stagnation zone. On the other hand, marked local erosion was detected downstream of the outer corner stagnation zone, as well as at the end of the channel outlet reach, where the rigid floor of the channel was practically reached, if not completely because of the free outflow. Furthermore, marked one-directional antidunes were found at the beginning of the channel inlet reach. These bed forms were created during the first stages of dam-break wave advance and progressively disappeared as they became closer to the corner region. It is worth mentioning that a slightly eroded zone appeared close to the inner corner, with a maximum erosion lower than 2 cm with respect to the original bed level.
In order to assess the suitability of the above numerical schemes for predicting the experimental observations, all the models were set with the parameters summarized in Table 3. A triangular unstructured mesh with 36,212 cells was used for all models. The resolution varied across the mesh. It was fixed at 0.5 m and 0.05 m in the reservoir and in the channel, respectively, but was increased to 0.005 m in the corner region in order to capture the local transient structures of the flow. All the simulations lasted 180 s, even though the bed evolution occurred mainly during the first 20 s and practically stopped after 120 s from the dam-break initiation.
Finally, particular attention was devoted to the downstream boundary condition. Since the flow conditions varied during the experiment, both supercritical and subcritical flows were obtained at the outlet. Different conditions should thus be applied depending on the Froude number. These conditions were described in detail by [6] and are briefly recalled here.
For supercritical flows, all the hydraulic information came from upstream and F ^ k = F L . However, when the flow was subcritical, an M2-type axis developed and went through the critical depth at the boundary. An imposed water depth boundary condition could then be applied so that:
q ^ k , n = q n , L + ( u n , L c L ) ( h c , L h L ) σ ^ k , n = σ n , L + ( u n , L c L ) 2 ( h c , L h L ) μ ^ k , n = μ n , L + v t , L ( q ^ k , n q n , L )
with h c , L being the critical depth of the cell left of the boundary.
When it comes to the conservation of the sediment however, information must come both from upstream and downstream, as suggested by the eigenvalues of the coupled model (see Equation (36)). Using the Rankine–Hugoniot relation and making the hypothesis that the sediment depth at the boundary was zero, another boundary condition could be developed for both super- and sub-critical flows:
q ^ b , k , n = q b , n , L λ h s , L
where λ = λ 1 and λ = λ ^ b for the coupled models (R-Cap, R-NCap, and HLLC-CM) and the decoupled one (HLLC-WCM), respectively, and where h s , L is the sediment depth of the cell left of the interface.

4. Results and Discussion

This section presents the results of the different models considered for the test case described in Section 3.3. First, the performances of the three capacity models were compared to each other. Secondly, the noncapacity features were applied and discussed through a comparison with the experimental data, the original R-Cap model results, and those acquired by a short sensitivity analysis.

4.1. Comparison between the Roe- and HLLC-Based Capacity Models

Figure 3 shows the temporal evolution of the water surface level ( w s l ) at the gauge points measured during the experiment. The arrival time of the dam-break wave was well predicted at the gauge points G 1 , G 2 , and G 3 , located upstream of the corner region, with the Roe- and HLLC-based models. However, at the gauge points placed downstream of the corner region ( G 4 and G 5 ), all the numerical models showed a shorter arrival time of the dam-break wave than those observed in the laboratory. Furthermore, the R-Cap model showed a slightly higher w s l than the HLLC solvers, especially at the wavefront, for all the gauge points measured. Nevertheless, the transient flow structure was reasonably well predicted by all the numerical models.
In order to assess the performance of the different numerical schemes to predict the bed changes caused by the dam-break wave, the final bed elevation obtained with the R-Cap- and HLLC-based solvers was compared against the photogrammetry measurements. Figure 4 shows the 2D maps of the bed elevation z b obtained with the R-Cap- and HLLC-based models at the final time t = 180 s. Several common aspects should be pointed out for the three models:
  • First, the three models were able to predict the bed degradation close to the outlet boundary reasonably well. However, none of them were able to obtain the bed forms observed in the experimental measurements at the beginning of the inlet reach;
  • Second, the R-Cap and the HLLC-WCM reproduced the main structures observed in the experiments for the final bed elevation well. However, the HLLC-CM led to a highly diffusive estimation of the bedload flux at the intercell edges, and the model was not able to reproduce the main features of the measured topography (see Figure 2). However, none of the schemes were able to accurately predict the absolute accumulation of bed material observed in the experiments downstream of the inner corner, nor the depth in the opposite eroded region;
  • Third, the R-Cap and HLLC-WCM computed a marked eroded zone close to the inner corner. Although slight erosion was observed in this region in the laboratory, both numerical models overestimated the bed degradation. That was in reality due to the formation and development of a 3D vortex downstream of the inner corner. Such a vortex cannot correctly be reproduced by depth-averaged models, which would neglect some vertical recirculation, which would itself hamper the erosion near the inner corner. Relaxing the assumption of a hydrostatic pressure distribution along the flow column might help to improve the accuracy of the numerical model predictions, since the vertical accelerations within the flow play an important role in this region.
Table 4 shows the global RMSE for the computed bed level z b 2D maps with respect to the photogrammetric data in Figure 2. The HLLC-WCM showed the lowest RMSE, followed by the R-Cap and the HLLC-CM. Since only a limited area—i.e., the bend region and roughly two meters downstream of it—was subject to large differences between the models, no model performed strikingly better than any other when looking at the RMSE. Nevertheless, these results were in line with the visual observations made previously.
In order to quantitatively identify the performance of the schemes more locally, final bed profiles taken along x = 6.34 m (A), x = 6.77 m (B), and y = 0.60 m (C) were obtained with the three numerical models and compared against the photogrammetric data, as shown in Figure 5.
Regarding Profile (A), the R-Cap model approximated the maximum of the accumulation region better than the HLLC-based models. However, the depth of the overeroded region close to the inner corner also increased with the R-Cap solver due to its lower numerical diffusion. Eventually, the R-Cap model showed the lowest RMSE for the bed level z b along this profile ( RMSE = 0.95 cm), whereas the HLLC-based models showed higher errors (see Table 5).
The HLLC-WCM performed slightly better than the R-Cap model along Profile (B) (see Figure 5), with RMSEs of 1.39 cm and 1.45 cm, respectively. The main reason for this similar performance was that the HLLC-WCM scheme better approximated the accumulation height in the outer corner upstream region and the bed slope at the outlet reach of the channel, whereas the R-Cap model slightly improved the prediction of the erosion zone downstream of the outer corner. However, the HLLC-CM again performed worse than the others ( RMSE = 1.50 cm), especially in the erosion zone downstream of the outer corner.
Finally, the three model performed quite similarly along Profile (C), with an RMSE below 0.5 cm and were able to predict the general trend of the bed change, as is shown in Figure 5. However, none of them were able to approximate the formation of dunes in the inlet reach of the channel. The formation of this kind of bed form is directly related to the vertical structure of the flow near the bed surface [19] and is hence difficult to mimic using depth-averaged bedload transport models.

4.2. Application of the R-NCap Model

The main feature of the R-NCap model is the progressive adaptation of the bedload transport rate q b to the local flow conditions until the equilibrium transport state is reached, contrary to the capacity models, which assume instantaneous adaptation. The celerity of this adaptation is controlled by the entrainment and detention constants k E and k D , respectively, but also directly depends on the dimensionless Shields stress excess Δ θ at the bed interface [22]. This is one of the main differences between the proposed R-NCap scheme and other nonequilibrium bedload models, which assume a constant value for the adaptation length L b [19,31,40] and compute the entrainment rate as:
η ˙ E = | q b * | ( 1 ξ ) L b
where | q b * | is the modulus of the equilibrium bedload rate (9).
Comparing the proposed R-NCap with former noncapacity models based on Equation (45) for the formulation of Meyer-Peter and Müller [26], it can be easily derived that the equivalent adaptation length applied by the R-NCap models scales with:
L b = 8 r s d s 1 ξ Δ θ k E
and hence, the equivalent adaptation length increases at regions where the boundary Shields stress excess is high. Furthermore, the smaller the entrainment constant k E , the larger the adaptation length L b is.
This property of the R-NCap model was used to improve the numerical prediction in the inner corner region. One of the main flaws in the numerical results obtained with the capacity R-Cap and HLLC-WCM models was the appearance of a marked overeroded region near the inner corner. This overeroded zone was not observed in the experimental measurements. The simulation showed that the marked erosion happened at the first stages of the dam-break progress throughout the corner region, when a vortex was formed downstream of the inner corner. The right panel of Figure 6 shows the velocity vectors at the corner region for the R-Cap simulation at t = 15 s. The formation of the vortex was clearly associated with the appearance of the overeroded zone, but by neglecting vertical features, the sediment was highly sheared and eventually advected downstream.
Moreover, the changes in the flow direction in the inner corner region led to high bed Shields stresses, which contributed to increasing the erosion within this region. The left panel of Figure 6 is a 2D map of the maximum values of bed Shields stress excess Δ θ as computed by the R-Cap model. The maximum Δ θ was around 1.0 in most of the channel, but increased in the inner corner region until reaching a maximum value greater than 2.5 , leading to a high erosion in this zone.
Considering the 2D map of the maximum Δ θ recorded for the R-Cap model and the features of the sediment used in this experiment, we analyzed the sensitivity of the R-NCap model by setting the entrainment and detention constants, k E and k D , respectively, to the values summarized in Table 6. Therefore, four simulations using the R-NCap model were carried out, varying the entrainment constant from k E = 1.60 to k E = 0.05 , but maintaining the k E / k D ratio equal to 20. It is worth mentioning that, considering the characteristic maximum value of the bed Shields stress excess Δ θ recorded for the simulation with the R-Cap model, the ratio k E / k D = 20 was chosen to ensure that the relation between the characteristic thickness of the bedload transport layer under equilibrium conditions and the sediment diameter remained η * / d s 10 (15) for all the simulations from T0 to T4.
Therefore, as k E decreased, the characteristic value of the equivalent adaptation length increased from L b 5 cm for the case T1 to L b 150 cm for the case T4. The increment of the adaptation length means that the spatial and temporal delay between the actual bedload transport rate and the flow carrying capacity became larger, and hence, the nonequilibrium states were enabled. Note that these values for the η * / d s ratio and the adaptation length L b only corresponded to the inner corner region, where the bed Shields stress was higher during the first stages of the dam-break wave. In other regions of the channel, the equivalent L b would be shorter and the η * / d s ratio smaller.
Figure 7 shows the final topography at time t = 180 s. When the R-NCap model was applied, the adaptation of the actual bedload solid rate to the flow capacity in the inner corner region became noninstantaneous. Hence, the appearance of the overeroded zone was not only delayed in time, but also moved further downstream of the inner corner. As k E decreased, the noncapacity state in that region was enabled, until the formation of the overeroded zone was avoided. The other main features of the topography observed in the laboratory were maintained, even if alterations in the bed level z b results also occurred in other regions of the channel.
Table 7 shows the global RMSE for the numerical topographies computed with the R-NCap model with respect to the photogrammetric data (Figure 2). The improvement of the global performance using the R-NCap model was not marked, but an optimal value for the entrainment constant k E = 0.20 could be found, corresponding to Test T2. Once again, the difference between the models may look limited regarding this indicator, because of the close results that they showed before the bend (Figure 5C), but it highlights the importance of properly calibrating the model. Poor calibration choices could indeed lead to worse results with the R-NCap than with the R-Cap model.
The final bed level z b computed with the R-NCap model is plotted along the profiles in Figure 8. The results from the R-Cap model and the photogrammetric data are also depicted for comparison purposes. For Profile (A), the activation of the nonequilibrium states led to the avoidance of the overeroded zone, but a spatial delay of the accumulation zone was also predicted, as well as a reduction of the accumulation height. Furthermore, the prediction of the bed slope at the outlet reach of the channel was increasingly worse as k E decreased. Despite the gain of accuracy allowed by the noncapacity feature near the inner corner, these worse and worse slope predictions led to higher RMSE values for the R-NCap model than for the R-Cap one along that profile (see Table 8).
For Profile (B), the R-NCap model improved the prediction of the bed slope in the channel outlet reach without significantly affecting other regions (see Figure 8). This also improved the RMSE of the bed level z b with respect to the photogrammetric data along this profile, in comparison with the R-Cap model, as highlighted by Table 8. The best agreement was given for k E [ 0.10 , 0.20 ] with a ratio k E / k D = 20 .

5. Conclusions

In this work, we compared the performance of different strategies for the resolution of the SWE + Exner system under capacity and noncapacity conditions. The selected capacity strategies involved the coupled HLLC model (HLLC-CM) [13], the weakly coupled HLLC model (HLLC-WCM) [6], and the fully coupled augmented Roe model (R-Cap) [14]. The three models were used to predict an experimental transient case [24] imposing a fixed set of noncalibrated configuration parameters. The experiment consisted of the propagation of a dam-break wave along a channel with a 90 bend and a uniform erodible bed, where the free surface evolution was measured at five different gauge points, while the final topography was measured after the channel drainage using photogrammetry.
Regarding the free surface evolution, although the R-Cap model showed a slightly higher w s l than the HLLC solvers, especially at the dam-break wavefront, the transient flow structure was reasonably well predicted by all the numerical models (Figure 3). When it comes to the final bed level estimation, the three models were able to predict the bed degradation close to the outlet boundary quite well (Figure 4). However, none of them were able to obtain the bed forms observed in the experimental measurements at the beginning of the inlet reach. The R-Cap and the HLLC-WCM models reproduced the main topographical structures observed during the experiments well, but the HLLC-CM led to an overly diffusive estimation of the bedload flux at the intercell edges, while the model was not able to clearly reproduce these main structures. The lowest global RMSE for the bed elevation z b was obtained with the HLLC-WCM, but the R-Cap model improved the prediction of the final bed level in some important regions (Figure 5). Therefore, the selection of the numerical approach for solving the intercell fluxes is not a trivial choice when a movable bed boundary is involved and should be carefully assessed. The HLLC-CM results indicated that the adaptation of fixed-bed diffusive solvers to erodible bed problems might lead to an erroneous bed evolution prediction. Hence, the use of highly diffusive centered solvers designed for fixed-bed conditions, as the widespread FORCE approach [41,42], should be thoroughly analyzed when movable boundaries are considered.
None of the capacity models were able to accurately predict the absolute accumulation of bed material observed in the experiments downstream of the inner corner, nor the depth in the opposite eroded region. Furthermore, the R-Cap and HLLC-WCM computed a markedly eroded zone close to the inner corner. Although slight erosion was observed in this region in the laboratory, both numerical models overestimated the bed degradation in this zone as a consequence of their inability to accurately capture all the 3D features of the vortex that formed and developed downstream of the inner corner. At this inner corner region, the vertical accelerations in the flow played a key role, especially at the first stages of the dam-break wave arrival, and assuming a nonhydrostatic pressure distribution along the flow column for the shallow-type equations might help to improve the accuracy of the predictions obtained with the bedload capacity formulation. Further research is required in order to apply the nonhydrostatic pressure assumption when a movable bottom boundary is involved.
In order to improve the model predictions near the inner corner, the experiment was also simulated using the generalized bedload transport model proposed by Martínez-Aranda et al. [22] (R-NCap) in order to analyze the influence of the transient nonequilibrium states on the final bed elevation. When the R-NCap model was applied, the adaptation of the actual bedload solid rate to the flow capacity in the inner corner region became noninstantaneous, and the appearance of the overeroded zone was delayed in time and space (Figure 7). As the equivalent adaptation length L b was increased by imposing smaller values of the nonequilibrium parameters k E and k D , the noncapacity state in the inner corner was activated until the formation of the overeroded zone was fully avoided. The other main features of the final topography observed in the laboratory were maintained, but alterations in the bed level z b also occurred in these other regions. The lowest global RMSE for the bed elevation z b was obtained with moderate values of the nonequilibrium parameters—around k E = 0.2 . When the noncapacity states were “forced” by imposing an excessive equivalent adaptation length L b , the global RMSE for the bed level z b increased, and the general performance of the R-NCap model was reduced. Even though the noncapacity approach can improve the model prediction in regions with complex transient processes, it requires a careful calibration of the nonequilibrium parameters, which can be a difficult task due to the lack of physical reference values. Furthermore, detailed experimental studies are required to assess whether the noncapacity states in the bedload transport actually occur in nature and how they affect the bed evolution in highly transient flows.

Author Contributions

Formal analysis, S.M.-A. and R.M.; investigation, S.M.-A. and R.M.; methodology, S.M.-A. and R.M.; project administration, R.M.; software, S.M.-A. and R.M.; supervision, S.M.-A.; validation, S.M.-A.; writing—original draft, S.M.-A. and R.M.; writing—review and editing, S.S.-F. and P.G.-N. All authors read and agreed to the published version of the manuscript.

Funding

This work was partially funded by the MINECO/FEDER under Research Project PGC2018-094341-B-I00 and by Diputacion General de Aragon, DGA, through Fondo Europeo de Desarrollo Regional, FEDER. It was also funded by Robin Meurice’s fellowship at the Fonds de la Recherche Scientifique—F.R.S.-FNRS. This research received no external funding for the APC.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
HLLCHarten–Lax–van Leer with Contact
SWEShallow-Water Equations
MPMMeyer-Peter and Müller
R-CapRoe-based Capacity Model
R-NCapRoe-based Noncapacity Model
HLLC-CMHLLC-based Coupled Model
HLLC-WCMHLLC-based Weakly Coupled Model
wslwater surface level

References

  1. Yang, C. Sediment Transport: Theory and Practice; McGraw-Hill Inc.: New York, NY, USA, 1996. [Google Scholar]
  2. Tingsanchali, T.; Chinnarasri, C. Numerical modelling of dam failure due to flow overtopping. Hydrol. Sci. J. 2001, 46, 113–130. [Google Scholar] [CrossRef]
  3. Juez, C.; Murillo, J.; García-Navarro, P. A 2D weakly coupled and efficient numerical model for transient shallow flow and movable bed. Adv. Water Resour. 2014, 71, 93–109. [Google Scholar] [CrossRef]
  4. Castro-Díaz, M.; Fernández-Nieto, E.; Ferreiro, A. Sediment transport models in Shallow Water equations and numerical approach by high order finite volume methods. Comput. Fluids 2008, 37, 299–316. [Google Scholar] [CrossRef] [Green Version]
  5. Martínez-Aranda, S.; Murillo, J.; García-Navarro, P. A 1D numerical model for the simulation of unsteady and highly erosive flows in rivers. Comput. Fluids 2019, 181, 8–34. [Google Scholar] [CrossRef] [Green Version]
  6. Meurice, R.; Soares-Frazão, S. A 2D HLL-based weakly coupled model for transient flows on mobile beds. J. Hydroinform. 2020, 22, 1351–1369. [Google Scholar] [CrossRef]
  7. Exner, F. Über die Wechselwirkung Zwischen Wasser und Geschiebe in Flüssen; Akademie der Wissenschaften: Wien, Austria, 1925. [Google Scholar]
  8. Hudson, J.; Sweby, P.K. Formulations for Numerically Approximating Hyperbolic Systems Governing Sediment Transport. J. Sci. Comput. 2003, 19, 225–252. [Google Scholar] [CrossRef]
  9. Struiksma, N. Mathematical modelling of bedload transport over nonerodible layers. In Proceedings of the IAHR Symposium on River, Coastal amd Estuarine Morphodynamics, Genoa, Italy, 6–10 September 1999; pp. 89–98. [Google Scholar]
  10. El Kadi Abderrezzak, K.; Paquier, A.; Gay, B. One-dimensional numerical modelling of dam-break waves over movable beds: Application to experimental and field cases. Environ. Fluid Mech. 2008, 8, 169–198. [Google Scholar] [CrossRef]
  11. Cao, Z.; Hu, P.; Pender, G. Multiple Time Scales of Fluvial Processes with Bed Load Sediment and Implications for Mathematical Modeling. J. Hydraul. Eng. 2011, 137, 267–276. [Google Scholar] [CrossRef]
  12. Martínez Aranda, S.; Navas-Montilla, A.; Lozano, A.; García-Navarro, P. Experimental study of resonant shallow flows past a lateral cavity: A benchmark test for high-resolution numerical models. In Proceedings of the EGU General Assembly, Online, 4–8 May 2020. [Google Scholar]
  13. Soares-Frazão, S.; Zech, Y. HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed. Int. J. Numer. Methods Fluids 2010, 66, 1019–1036. [Google Scholar] [CrossRef]
  14. Martínez-Aranda, S.; Murillo, J.; García-Navarro, P. Comparison of new efficient 2D models for the simulation of bedload transport using the augmented Roe approach. Adv. Water Resour. 2021, 153, 103931. [Google Scholar] [CrossRef]
  15. Soares-Frazao, S.; Canelas, R.; Cao, Z.; Cea, L.; Chaudhry, H.M.; Moran, A.D.; Kadi, K.E.; Ferreira, R.; Cadórniga, I.F.; Gonzalez-Ramirez, N.; et al. Dam-break flows over mobile beds: Experiments and benchmark tests for numerical models. J. Hydraul. Res. 2012, 50, 364–375. [Google Scholar] [CrossRef]
  16. Caviedes-Voullième, D.; Morales-Hernández, M.; Juez, C.; Lacasta, A.; García-Navarro, P. Two-Dimensional Numerical Simulation of Bed-Load Transport of a Finite-Depth Sediment Layer: Applications to Channel Flushing. J. Hydraul. Eng. 2017, 143, 04017034. [Google Scholar] [CrossRef]
  17. Cao, Z.; Li, Z.; Pender, G.; Hu, P. Non-capacity or capacity model for fluvial sediment transport. Proc. Inst. Civ. Eng. Water Manag. 2012, 165, 193–211. [Google Scholar] [CrossRef]
  18. Zech, Y.; Soares-Frazão, S.; Spinewine, B.; Savary, C.; Goutière, L. Inertia effects in bed-load transport models. Can. J. Civ. Eng. 2009, 36, 1587–1597. [Google Scholar] [CrossRef]
  19. Wu, W. Computational River Dynamics; Taylor & Francis/Balkema: Leiden, The Netherlands; CRC Press: Boca Raton, FL, USA, 2007. [Google Scholar]
  20. Charru, F. Selection of the ripple length on a granular bed sheared by a liquid flow. Phys. Fluids 2006, 18, 121508. [Google Scholar] [CrossRef]
  21. Fernández-Nieto, E.; Lucas, C.; Morales-de Luna, T.; Cordier, S. On the influence of the thickness of the sediment moving layer in the definition of the bedload transport formula in Exner systems. Comput. Fluids 2014, 91, 87–106. [Google Scholar] [CrossRef] [Green Version]
  22. Martínez-Aranda, S.; Murillo, J.; García-Navarro, P. A comparative analysis of capacity and noncapacity formulations for the simulation of unsteady flows over finite-depth erodible beds. Adv. Water Resour. 2019. [Google Scholar] [CrossRef]
  23. Spinewine, B.; Zech, Y. Small-scale laboratory dam-break waves on movable beds. J. Hydraul. Res. 2007, 45, 73–86. [Google Scholar] [CrossRef]
  24. Meurice, R.; Martínez-Aranda, S.; Ebrahimi, M.; García-Navarro, P.; Soares-Frazão, S. Laser Profilometry to measure the bed evolution in a dam-break flow. 2021; submitted. [Google Scholar]
  25. Murillo, J.; García-Navarro, P. An Exner-based coupled model for two-dimensional transient flow over erodible bed. J. Comput. Phys. 2010, 229, 8704–8732. [Google Scholar] [CrossRef]
  26. Meyer-Peter, E.; Müller, R. Formulas for Bed-Load Transport. In Proceedings of the 2nd Meeting on International Association on Hydraulic Structures Research, Delft, The Netherlands, 7 June 1948; pp. 39–64. [Google Scholar]
  27. Nielsen, P. Coastal Bottom Boundary Layers and Sediment Transport; Advanced Series on Ocean Engineering; World Scientific: Singapore, 1992. [Google Scholar]
  28. Luque, R.F.; Beek, R.V. Erosion And Transport Of Bed-Load Sediment. J. Hydraul. Res. 1976, 14, 127–144. [Google Scholar] [CrossRef] [Green Version]
  29. Wong, M.; Parker, G. Reanalysis and Correction of Bed-Load Relation of Meyer-Peter and Müller Using Their Own Database. J. Hydraul. Eng. 2006, 132, 1159–1168. [Google Scholar] [CrossRef] [Green Version]
  30. Grass, A. Sediments Transport by Waves and Currents; Department of Civil Engineering, University College: London, UK, 1981. [Google Scholar]
  31. Cao, Z.; Hu, P.; Pender, G.; Liu, H.H. Non-capacity transport of nonuniform bed load sediment in alluvial rivers. J. Mt. Sci. 2016, 13, 377–396. [Google Scholar] [CrossRef]
  32. Martínez-Aranda, S.; Murillo, J.; García-Navarro, P. A new 2D bedload transport model based on noncapacity approach to overcome the problems associated with finite-depth sediment layers. In Proceedings of the 10th International Conference on Fluvial Hydraulics (River Flow 2020), Delft, The Netherlands, 7–10 July 2020. [Google Scholar]
  33. Godlewski, E.; Raviart, P.A. Numerical Approximation of Hyperbolic Systems of Conservation Laws; Springer Science & Business Media: Berlin/Heisenberg, Germany, 1996. [Google Scholar]
  34. Murillo, J.; Navas-Montilla, A. A comprehensive explanation and exercise of the source terms in hyperbolic systems using Roe type solutions. Application to the 1D-2D shallow water equations. Adv. Water Resour. 2016, 98, 70–96. [Google Scholar] [CrossRef] [Green Version]
  35. Fraccarollo, L.; Capart, H.; Zech, Y. A Godunov method for the computation of erosional shallow water transients. Int. J. Numer. Methods Fluids 2003, 41, 951–976. [Google Scholar] [CrossRef]
  36. Harten, A.; Lax, P.; van Leer, B. On upstream differencing and Godunov type methods for hyperbolic conservation laws. SIAM Rev. 1983, 25, 35–61. [Google Scholar] [CrossRef]
  37. Toro, E. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction; Springer: Berlin/Heisenberg, Germany, 1997. [Google Scholar]
  38. Toro, E.F.; Spruce, M.; Spears, W. Restoration of the contact surface in the HLL Riemann solver. Shock Waves 1994, 4, 25–34. [Google Scholar] [CrossRef]
  39. Franzini, F.; Soares-Frazão, S. Coupled finite-volum scheme with adapted Augmented Roe scheme for simulating morphological evolution in arbitrary cross-sections. J. Hydroinform. 2018, 20, 1111–1130. [Google Scholar] [CrossRef]
  40. Paquier, A.; Goutal, N. Dam and levee failures: An overview of flood wave propagation modeling. Houille Blanche Rev. Int. L’eau 2016, 5–12. [Google Scholar] [CrossRef] [Green Version]
  41. Dumbser, M.; Hidalgo, A.; Castro, M.; Parés, C.; Toro, E.F. FORCE schemes on unstructured meshes II: Non-conservative hyperbolic systems. Comput. Methods Appl. Mech. Eng. 2010, 199, 625–647. [Google Scholar] [CrossRef]
  42. Canestrelli, A.; Dumbser, M.; Siviglia, A.; Toro, E.F. Well-balanced high-order centered schemes on unstructured meshes for shallow water equations with fixed and mobile bed. Adv. Water Resour. 2010, 33, 291–303. [Google Scholar] [CrossRef]
Figure 1. (a) Plane view of the experimental flume and position of the gauges G1 to G5; (b) vertical cut taken along the longitudinal axis of the first part of the flume.
Figure 1. (a) Plane view of the experimental flume and position of the gauges G1 to G5; (b) vertical cut taken along the longitudinal axis of the first part of the flume.
Water 13 01840 g001
Figure 2. Final topography obtained after the channel drainage using photogrammetry and averaged over the three experimental runs available.
Figure 2. Final topography obtained after the channel drainage using photogrammetry and averaged over the three experimental runs available.
Water 13 01840 g002
Figure 3. w s l (cm) at gauge points G 1 , G 2 , G 3 , G 4 , and G 5 .
Figure 3. w s l (cm) at gauge points G 1 , G 2 , G 3 , G 4 , and G 5 .
Water 13 01840 g003
Figure 4. Bed elevation z b 2D maps obtained with the R-Cap- and HLLC-based models at t = 180 s.
Figure 4. Bed elevation z b 2D maps obtained with the R-Cap- and HLLC-based models at t = 180 s.
Water 13 01840 g004
Figure 5. Final bed level profiles along x = 6.34 m (A), x = 6.77 m (B), and y = 0.60 m (C) with the R-Cap, HLLC-CM, and HLLC-WCM. Experimental photogrammetric data are also plotted.
Figure 5. Final bed level profiles along x = 6.34 m (A), x = 6.77 m (B), and y = 0.60 m (C) with the R-Cap, HLLC-CM, and HLLC-WCM. Experimental photogrammetric data are also plotted.
Water 13 01840 g005
Figure 6. Flow structure with the R-Cap model: (left) 2D map of the maximum Δ θ ; (right) zoom of the inner corner region. The velocity vectors are superimposed with the bed elevation.
Figure 6. Flow structure with the R-Cap model: (left) 2D map of the maximum Δ θ ; (right) zoom of the inner corner region. The velocity vectors are superimposed with the bed elevation.
Water 13 01840 g006
Figure 7. Bed elevation z b 2D maps obtained with the R-NCap model at t = 180 s.
Figure 7. Bed elevation z b 2D maps obtained with the R-NCap model at t = 180 s.
Water 13 01840 g007
Figure 8. Final bed level profiles along x = 6.34 m (A) and x = 6.77 m (B) with the R-NCap model. Experimental photogrammetric data and results from the R-Cap model are also plotted.
Figure 8. Final bed level profiles along x = 6.34 m (A) and x = 6.77 m (B) with the R-NCap model. Experimental photogrammetric data and results from the R-Cap model are also plotted.
Water 13 01840 g008
Table 1. Coefficients c, m 1 , m 2 , and θ c for different capacity solid transport rate formulations (9).
Table 1. Coefficients c, m 1 , m 2 , and θ c for different capacity solid transport rate formulations (9).
Formulationc m 1 m 2 θ c
MPM [26]803/20.047
Nielsen [27]121/210.047
Fernández-Luque [28]5.703/20.037
Wong [29]3.9703/20.0495
Table 2. Noncapacity Grass interaction factor G for different solid transport rate formulations.
Table 2. Noncapacity Grass interaction factor G for different solid transport rate formulations.
Formulation Γ 1 ( h ) Γ 2 ( θ ) Γ 3 ( η ) θ c
MPM n b 3 g ( r s 1 ) h 8 Δ θ θ 3 / 2 r s k D k E η d s 0.047
Nielsen n b 3 g ( r s 1 ) h 12 θ r s k D k E η d s 0.047
Fernández-Luque n b 3 g ( r s 1 ) h 5.7 Δ θ θ 3 / 2 r s k D k E η d s 0.037
Wong n b 3 g ( r s 1 ) h 3.97 Δ θ θ 3 / 2 r s k D k E η d s 0.0495
Table 3. Setup of the simulations.
Table 3. Setup of the simulations.
Water density ρ w 1000 kg/m3
Solid density ρ s 2650 kg/m 3
Solid particles’ diameter d s 1.7 mm
Manning’s roughness coeff. n b 0.0165 sm 1 / 3
Bed porosity ξ 0.44
Bedload formulationMeyer-Peter and Müller [26]
Critical Shields stress θ c 0.047
Table 4. Global RMSE for the bed level z b with the HLLC-CM, HLLC-WCM, and R-Cap.
Table 4. Global RMSE for the bed level z b with the HLLC-CM, HLLC-WCM, and R-Cap.
Global RMSE for z b (cm)
HLLC-CMHLLC-WCMR-Cap
1.070.961.02
Table 5. Bed-level z b RMSE for Profiles x = 6.34 m, x = 6.77 m, and y = 0.60 m with the R-Cap, HLLC-CM and HLLC-WCM.
Table 5. Bed-level z b RMSE for Profiles x = 6.34 m, x = 6.77 m, and y = 0.60 m with the R-Cap, HLLC-CM and HLLC-WCM.
Profile z b RMSE (cm)
HLLC-CMHLLC-WCMR-Cap
x = 6.34 m1.531.470.95
x = 6.77 m1.501.391.45
y = 0.60 m0.500.480.45
Table 6. Noncapacity setup for the analysis of the R-NCap model’s behavior.
Table 6. Noncapacity setup for the analysis of the R-NCap model’s behavior.
TestModel k E / K D K E K D Δ θ η * / d s (15) L b (46)
(-)(-)(-)(-)(-)(cm)
T0R-Cap20--1.49.4-
T1R-NCap201.600.0801.49.44.8
T2R-NCap200.200.0101.49.438.1
T3R-NCap200.100.0051.49.476.1
T4R-NCap200.050.00251.49.4152.3
Table 7. Global RMSE for the bed level z b with R-Cap and R-NCap models.
Table 7. Global RMSE for the bed level z b with R-Cap and R-NCap models.
TestGlobal RMSE for z b (cm)
T0: R-Cap1.02
T1: R-NCap k E = 1.60 1.01
T2: R-NCap k E = 0.20 0.98
T3: R-NCap k E = 0.10 1.02
T4: R-NCap k E = 0.05 1.15
Table 8. Bed level z b RMSE for the profiles x = 6.34 m and x = 6.77 m with the R-Cap and R-NCap models.
Table 8. Bed level z b RMSE for the profiles x = 6.34 m and x = 6.77 m with the R-Cap and R-NCap models.
Test z b RMSE (cm)
x = 6.34 m x = 6.77 m
T0: R-Cap0.951.45
T1: R-NCap k E = 1.60 1.701.32
T2: R-NCap k E = 0.20 1.611.05
T3: R-NCap k E = 0.10 1.511.13
T4: R-NCap k E = 0.05 1.691.34
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Martínez-Aranda, S.; Meurice, R.; Soares-Frazão, S.; García-Navarro, P. Comparative Analysis of HLLC- and Roe-Based Models for the Simulation of a Dam-Break Flow in an Erodible Channel with a 90 Bend. Water 2021, 13, 1840. https://doi.org/10.3390/w13131840

AMA Style

Martínez-Aranda S, Meurice R, Soares-Frazão S, García-Navarro P. Comparative Analysis of HLLC- and Roe-Based Models for the Simulation of a Dam-Break Flow in an Erodible Channel with a 90 Bend. Water. 2021; 13(13):1840. https://doi.org/10.3390/w13131840

Chicago/Turabian Style

Martínez-Aranda, Sergio, Robin Meurice, Sandra Soares-Frazão, and Pilar García-Navarro. 2021. "Comparative Analysis of HLLC- and Roe-Based Models for the Simulation of a Dam-Break Flow in an Erodible Channel with a 90 Bend" Water 13, no. 13: 1840. https://doi.org/10.3390/w13131840

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop