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

: In geophysical surface ﬂows, 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 ﬂow features. The resulting system of equations, composed of the shallow water equations for the ﬂow (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 ﬂaws 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 ﬂow. 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 ﬁxed setup parameters. The results indicate that the discrete strategy used to compute the intercell ﬂuxes signiﬁcantly 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.


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.

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 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, (τ bx , τ by ) 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: where ξ denotes the bulk porosity of the bed layer and (q bx , q by ) 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: where U is the vector of the conserved variables: The vector S b (U) accounts for the momentum source term associated with the variation of the pressure force on the bed surface: The source vector S τ (U) denotes the momentum dissipation due to the basal resistance: where the components of the shear resistance at the bed surface τ b = (τ bx , τ by ) are commonly estimated using the quadratic pure turbulent relation as: with C f = n 2 b 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 bx , q by ) 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: 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.  [29] 3.97 0 3/2 0.0495 The term ∆θ denotes the positive excess of the Shields stress θ over the critical Shields stress value for the incipient motion θ c , expressed as: 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: 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: 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: 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 η: 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. Table 2. Noncapacity Grass interaction factor G for different solid transport rate formulations.

Formulation
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: whereas the noncapacity assumption requires estimating the temporal evolution of the bedload transport layer thickness: (1 − ξ) ∂η ∂t 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]:η Therefore, when equilibrium conditions are reached, η = η * , leading toη D =η E , and hence, Equation (16) reduces to the widespread Exner equation [7].

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.

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: 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 −1 k 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 , beingn the outward normal unit vector and (x n ,ŷ 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: 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: 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: Following Murillo and Navas-Montilla [34], the conserved variable increment δÛ 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: The approximate upwind solution F ↓− k for the numerical flux at the k-th cell edge is defined as: 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: where (∆θ) n+1 i 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: with F η↓ k being the numerical flux at the k-th intercell edge for the homogeneous transport equation and (η D −η E ) n i 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 direc-tionx n to the edge, and the numerical flux at the intercell interface is approximated using the linearized homogeneous Riemann problem [32]: where the virtual bedload wave celerity λ η,k is defined as: Therefore, the intercell numerical flux for the transport thickness is computed as: 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: with Γ 1 , Γ 2 , and Γ 3 as in Table 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.

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: withF 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]: whereF 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: whereas for the bed level evolution equation, S + is expressed as: 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 ,ŷ t ) instead of the global one (x, y), thanks to the rotational invariance property [37]: 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: The eigenvalues of M can then be approximated as proposed by [13]: with c = gh.
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: 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: where σ n = hu 2 n + gh 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 .

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: Regarding the update of z b , the same finite-volume approach is followed, but a new estimation ofq b,k,n must be provided: whereλ b is computed as proposed by [3] and is quite similar to Equation (27): 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 theq 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].

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

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. Figure 3 shows the temporal evolution of the water surface level (wsl) at the gauge points measured during the experiment. The arrival time of the dam-break wave was well predicted at the gauge points G1, G2, and G3, located upstream of the corner region, with the Roe-and HLLC-based models. However, at the gauge points placed downstream of the corner region (G4 and G5), 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 wsl 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.

Comparison Between the Roe-and HLLC-Based Capacity 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 HLLCbased 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.

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

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