Coupled method for the numerical simulation of 1D shallow water and Exner transport equations in channels with variable cross-section

. This work is focused on the a numerical ﬁnite volume scheme for the resulting coupled shallow water-Exner system in 1D applications with arbitrary geometry. The mathematical expression modeling the the hydrodynamic and morphodynamic components of the physical phenomenon are treated to deal with cross-section shape variations and empirical solid discharge estima-tions. The resulting coupled system of equations can be rewritten as a non-conservative hyperbolic system with three moving waves and one stationary wave to account for the source terms discretization. But, even for the simplest solid transport models as the Grass law, to ﬁnd a linearized Jacobian matrix of the system can be a challenge if one considers arbitrary shape channels. More-over, the bottom channel slope variations depends on the erosion-deposition mechanism considered to update the channel cross-section proﬁle. In this paper a numerical ﬁnite volume scheme is proposed, based on an augmented Roe solver (ﬁrst order accurate in time and space) and dealing with solid transport ﬂux variations caused by the channel geometry changes. Channel cross-section variations lead to the appearance of a new solid ﬂux source term which should be discretized properly. Comparison of the numerical results for sev-eral analytical and experimental cases demonstrate the e ﬀ ectiveness, exact well-balanceness and accuracy of the scheme.


Introduction
Sediment transport in rivers is usually classified into two different phenomena: bed load movement and suspended material transport. The bed load process is usually modeled by means of set of equations that include hydrodynamic and morphodynamic components. The hydrodynamic part can be described by the shallow water equations, commonly used to study water movement in rivers, coastal and channels. On the other hand, the morphodynamical component is commonly represented by a solid mass continuity equation, depending on the solid transport flux [1].
Two-dimensional models can be used to simulate flow and sediment transport using a refined representation of topography and local hydraulic effects. However, their application to natural river cases is still restricted due to the computational time required and the amount of field data needed for the model calibration. Therefore, reported two-dimensional models were only applied to reach-scale domain cases with short event time duration [2,3]. On the other hand, one-dimensional modes require less field data, offering a higher computational efficiency.
Despite numerical modeling of free-surface flows with bed load transport involves transient flow and movable boundaries, the conventional 1D methods found in literature usually decouple the hydrodynamic part and the solid transport equation [4][5][6]. Coupled model reported were exclusively developed for constant cross-sectional cases [7,8]. Therefore, a numerical coupled model able to simulate complex geometries and demonstrate its performance in 1D realistic applications is still required.

Mathematical model
Free-surface flow movement in one-dimensional practical applications can be modeled by the 1D shallow water or Saint-Venant equations (1) and (2) for the mass and momentum conservation, respectively.
where A(x, t) is the wetted cross-sectional area, Q(x, t) is the flow discharge, S f is the friction slope, I 1 represents the hydrostatic pressure forces and I 2 accounts for the pressure forces on the channel walls. S 0 is the bed slope, which can be expressed as: being z b (x, A b ) the lowest cross-sectional level of the channel bed above a datum and A b (x, t) the conserved solid area. In the momentum equation (2), the conservative flux should be separated into the component due to geometrical variations exclusively and the momentum flux caused by the hydrodynamic features [9], avoiding the computations of the integral I 1 and I 2 .
On the other hand, bed-load mass conservation is usually modeled by means of the Exner continuity equation [10]: where Q s (x, A, Q) is the total solid discharge at the cross-section. The total bed load discharge is defined considering a constant bed load transport rate per unit width q s (A, Q) applied to the whole cross-section [11]: Hence the Exner equation (4) must be rewritten distinguishing the solid flux component due to the flow variations from the solid flux component caused by the cross-section changes, which is included as a sediment mass source term.
The Jacobian matrix for the coupled system composes by (1), (2) and (6) is singular, since it does not depend on the conserved solid area A b . This should create difficulties to implement a numerical scheme for this formulation. To overcome this problem and taking into account (3), we can express the system of conservation laws including a non-conservative flux T b = gA ∂z b ∂A b δA b at of the momentum equation and a modified source terms vector: An estimation of the solid transport rate per unit width q s (A, Q) using the Grass law [12] allow us to obtain a explicit expression for the partial derivatives of the sediment discharge in J(x, U): where A g [T 2 /L] is a constant called Grass coefficient. Following [13], partial derivatives ∂z b /∂x and ∂h/∂x should be avoided in the source terms formulation in order to correctly approach the discrete increments of z b and h. Therefore, taking into account that ∂z b /∂A b = c 2 b /(gA) and ∂h/∂A = 1/B, the modified source terms vector could be expressed in terms of total derivatives.
Finally, by simply combining one can obtain the complete matrix for the coupled system, which should be used to design the numerical method.
3 Numerical scheme The system of equations (7) is solved according to a finite volume method, in which the domain is divided in computational cells Ω i of constant size ∆x = x i+1/2 − x i−1/2 . Assuming a piecewise representation of the conserved variables and fluxes, the first order Godunov method provides a way to update the averaged values of the solution U n i to the next time step t n+1 . The numerical scheme could be rewritten as: being δ E| i+1/2 = E n i+1 − E n i the difference of the physical fluxes at the neighbouring cells. The definition of the numerical scheme in the Godunov method must be completed by the definition of an approximate solver for the local Riemann problem governed by the fluxes E n i+1 and E n i .

Approximate augmented Roe's solver for the Grass model
For the Roe's approximate Riemann solver, at each edge i + 1/2 separating cell i and i + 1, we should define a linearized local Riemann problem which is described by the following hyperbolic system of equations: Therefore, this allows us to express the flux differences and the source terms at the cell edge i+1/2 as their projection onto the right eigenvector base e m for the approximate Jacobian matrix M( U i , U i+1 ) of the hyperbolic system (11): where λ m are the eigenvalues of M, α m the wave strength and β m the source strength. Therefore, it is necessary the construction of two proper linearized matrix J( U i , U i+1 ) and H( U i , U i+1 ), with the following Roe's averaged components: Assuming that the Grass coefficient A g remains constant through the cell interface, the eigenvalues of the hyperbolic system (11) are the roots of the characteristic polynomial of matrix M( U i , U i+1 ), defined as: According to [14], this kind of hyperbolic systems has always two eigenvalues of the same sign than the flow velocity and the other one with opposite direction. Numerically, the roots of P M ( λ) could also be calculated by the Cardano-Vieta formula.
Finally, the augmented Roe's upwind scheme considering A g = cte through the intercell could be written as: Numerical stability region is controlled by a CFL condition, limiting the time step ∆t such that there is no interaction between waves from neighboring Riemann problems.

Approximate augmented Roe's solver for empirical models
In realistic application, A g is evaluated as a function of the bottom shear stress |τ b |. In [15], Juez et al. reported a method to adjust the q s value calculated by the Grass law to other empirical models, leading to different Grass coefficients A g (h, u) at cells i and i + 1.
Evaluating the solid discharge at each side of the cell edge considering an averaged Grass coefficientĀ g | i+1/2 = (A gi + A gi+1 )/2, the total flux differences at the intercell i + 1/2 can be decomposed as follows: whereQ si = Q s (Ā g | i+1/2 , U i ). Therefore, δQ s | i+1/2 is evaluated by the method reported in section 3.1. However, to ensure bed load conservation, the solid flux difference δ(Q s −Q s )| i+1/2 should be included by means of a new cellwise contribution δ F * | i+1/2 that does not require any splitting [16].
From now on, the numerical scheme (17) will be referred as Conservative Coupled Model (CCM) to clearly distinguish it from (15), called NCCM (Non-Conservative Coupled Model).

Numerical tests 4.1 Equilibrium slopes for steady flow regimes
The aim of these test cases is to assess the performance of the proposed method to converge to the exact solution in steady flow regimes. The geometry is a rectangular channel which  Results for the contraction case ( Fig. 1-left) show a bed erosion downstream caused by an increment in the solid transport rate per unit width q s . On the other hand, the width contraction leads to a progressively lower solid discharge, hence the bed level increases downstream ( Fig. 1-right). However, the total solid discharge in the channel remains constant for both cases, indicating that the steady state is reached. Numerical results for the bed level agree perfectly with the exact solution derived from imposing equilibrium conditions (S 0 = S f ).

Dam-break with analytical solution
Three test problem with exact solution presented in [16] are reported in this section. The tests are one-dimensional Riemann problem for movable bed equation, in which the friction shear stress term has been neglected in the momentum equation. All simulation were made with CFL = 1, ∆x = 0.01 m and sediment layer porosity p = 0.4. Results predicted by the CCM model are presented for a final simulation time t = 2 s. Cases DAS1 and DAS2 prove the numerical scheme to converge to the exact solution of both Riemann problems. The CCM model correctly reproduces the rarefaction waves involved in these two problems. The contact and shock waves were also captured accurately, in terms of strength and position, for all the conserved variables (Fig. 2). Case DAS3 clarifies the influence of the cellwise flux term in the numerical formulation when A g is not constant through the local Riemann problem. The CCM model is able to describe properly the solution structure for all the conserved variables and all the waves were captured accurately. Furthermore, with the NCCM model the numerical solution degenerates clearly due to the lack of conservation of the solid bed material volume (Fig. 3).

Dam-break over bed step and cross-section changes
In this test, a dam-break with h L = 2.5 m and h R = 0.7 m is simulated over a bed step (z bL = 1.5 m and z bR = 1.0 m) for different cross-section shapes both upstream and downstream. A progressive linear cross-section transition occurs between x ∈ [−1 m, +1 m]. Upstream and downstream of this transition region, the channel is prismatic with the cross-sections indicated in Table 2. All the cases were simulated with CFL = 1, ∆x = 0.05 m, p = 0.6, n = 0.025sm −1/3 and q s was estimated by the MPM model. The solution was evolved until a final simulation time t = 3.0 s.  The numerical model is able to handle the solid discharge variations caused by sudden expansions and contractions in rectangular channels (Fig. 4). An expansion downstream leads to higher erosion upstream the bed step. On the contrary, a contraction downstream caused a lower erosion upstream the bed step. Finally, the of the scheme to handle cross-section shape changes is evaluated by Test 4 and Test 5. Results for a prismatic rectangular (Test 1) and prismatic trapezoidal (Test 6) channels were selected as reference. Results depicted in Fig. 5 show the influence of the cross-section shape change over the erosion on the bed step, leading to marked differences with respect to the case in which only width varies but the shape remains rectangular.

Conclusions
A new finite volume scheme has been proposed for the coupled system of shallow water and Exner equations which is applicable to 1D channels with arbitrary geometry. The equations were treated to deal with cross-section shape variations by distinguishing the intercell conservative fluxes due to geometry variations from that caused by the flow features. The resulting coupled system of equations was rewritten as a non-conservative hyperbolic system with three non-linear characteristic fields. An upwind augmented Roe's scheme is proposed for both the solid discharge evaluated by the Grass law and by others empirical models. The scheme is able to describe properly the exact solution structure for the conserved variables in all the cases tested. Furthermore, we remark that the ability of the model to handle cross-section shape changes was also proved.