A robust two-dimensional model for highly sediment-laden unsteady ﬂ ows of variable density over movable beds

In this work, a novel 2D depth-integrated numerical model for highly sediment-laden shallow ﬂ ows over non-uniform erodible beds is presented, including variable density and exchange between the bed layer and the water – sediment mixture ﬂ ow. The system of equations is formed by the 2D conservation equations for the mass and momentum of the mixture, the mass conservation equation for the different sediment size-classes transported in the ﬂ ow and the bed evolution equation. The depth-averaged mixture density varies according to the volumetric concentration of the different sediment size-classes that can be incorporated from the bed to the ﬂ ow and transported as suspended materials. The rheological behaviour of the ﬂ ow is directly controlled by the properties of the mixture. A new x-split augmented Roe (xA-Roe) scheme is derived to solve the coupled ﬂ ow and suspended solid-phase equations in both structured and unstructured meshes. The numerical scheme is de ﬁ ned to properly include density variations and momentum source terms, retaining a well-balanced ﬂ ux formulation in steady states and the correct treatment of the wet – dry fronts. The numerical scheme is assessed with steady and transient cases involving highly sediment-laden ﬂ ows, demonstrating its accuracy, stability and robustness in the presence of complex bed topography, wetting – drying fronts and rapid morphological changes. words (xA-Roe) solver is proposed in this work, including transport and exchange of different sediment size-classes, complex friction term integration and wet – dry front treatment. (cid:129) The new xA-Row scheme is demonstrated to be robust, accurate and ef ﬁ cient for a wide range of unsteady problems involving variable-density ﬂ ows over non-uniform erodible beds.


INTRODUCTION
Water-sediment mixture flows are widely present in environmental and geophysical processes such as rivers and estuaries morphodynamics. When water bodies are well mixed and the transported suspended sediment is distributed uniformly over the water column, assuming a single value of the mixture density along the flow depth allows the use of depth-integrated models. In natural river modelling, the description usually includes the transport of different suspended sediment fractions with a total concentration low enough to assume that the mixture density is equal to the density of the water. This is no longer valid in flows dealing with highly laden water-sediment mixtures, where the density of the mixture can be more than twice the density of the water. These This paper is structured as follows: in the next section, the two-dimensional quasi-single-phase equations for the water-sediment mixture with different suspended sizeclasses over non-uniform erodible bed are presented; The section "xA-Roe Solver for 2D Sediment-Laden Flows" is devoted to describe the proposed xA-Roe scheme for variable-density flows, paying especial attention to the formulation of the numerical fluxes at the cell edges and the correct integration of the momentum source terms; in the section "Numerical results", the numerical scheme is validated against steady-state cases with exact solution and idealized highly transient sediment-laden tests, demonstrating its stability and robustness in the presence of complex bed topography, wetting-drying fronts and rapid morphological changes; finally, the conclusions are drawn in the final section.

2D SEDIMENT-LADEN FLOW, SUSPENDED SOLID SIZE-PHASES AND BED EVOLUTION EQUATIONS
The model presented in this work involves the following assumptions: (i) shallow-water approach: the flow is confined to a layer which is thin compared with the horizontal scale of interest, leading to the hydrostatic pressure assumption; (ii) multicomponent flow: the mixture of water and suspended sediment particles is described by using the continuum approach and assuming the same velocity for the liquid and the suspended solid phase; and (iii) the different sediment size-classes presented in the flow are distributed uniformly in the column. Accordingly, ϕ p represents the scalar depth-averaged volumetric concentration of pth sediment size-class, with p ¼ 1, . . . , N and N the number of sediment size-classes transported. The mixture density is given by ρ ¼ ρ w r, where ρ w is the density of the water and r represents the relative density of the bulk mixture to that of the clean water: where ϕ χ is the global modified sediment concentration in the mixture and χ p ¼ (ρ p À ρ w )=ρ w is the relative buoyant density of the pth solid phase, ρ p being the density of the sediment particles.
The relevant formulation for the two-dimensional sediment-laden flow model ( Figure 1) includes the depthaveraged equations for the mixture mass and momentum conservation, the mass conservation equation for the suspended solid phase and the mass conservation equation for the bed layer. The complete 2D system can be expressed in vector form, following a global coordinate system (GCS) (Juez et al. ), as: where U is the vector of conserved variables, F(U) and G(U) are the convective fluxes along the x and y global coordinates, respectively, S b (U) is the momentum source term associated with the variation of the pressure force on the bottom, S τ (U) is the momentum dissipation due to the boundary shear stress between the mixture flow and the bed layer and E b (U) accounts for the mass net exchange flux between the mixture flow and the bed layer: where h represents the mixture flow depth, (u, v) are the components of the depth-averaged flow velocity vector u along the x and y coordinates, respectively, z b is the bed elevation and g ψ ¼ g cos 2 ψ is the bed-normal projection of the gravity in the GCS, g being the gravitational acceleration and cos ψ the direction cosine of the bed normal with respect to the vertical axis (Juez et al. ).
The components of the bed slope source term vector S b (U) (5) along the x and y coordinates are expressed as: and those of the frictional dissipation at the bottom boundary S τ (U) can be expressed as:  • The turbulent effects near the bed can be expressed as a function of the Manning's roughness coefficient n: n being the Manning roughness coefficient.
• In case of a pure Newtonian fluid, the viscous stress can be estimated as: μ being the dynamic viscosity.
• The Coulomb-type laws for granular material estimated the tangential shear stress as a function of the internal stability angle of the mixture fluid θ b : • Non-Newtonian Bingham-type fluids do not flow until a threshold value of the tangential stress, the yield stress τ y , is reached. During the movement, the boundary shear stress τ b is characterized by means of a cubic equation accounting for the plastic viscosity of the sediment-water mixture: • If the ratio jτ b j=τ y is smaller than 0.5, the full Bingham relation (12) can be reduced to: All these different types of tangential stresses act simultaneously along the mixture column and hence participate in the depth-averaged tangential forces in the mixture momentum equations, i.e. jτ b j ¼ f(τ t , τ y , τ μ , τ f ) (Murillo & García-Navarro ). The different friction formulations considered in this work have been summarized in Table 1.
In the source term E b (U) (6), the global net exchange fluxes for the mixture, solid phase and bed layer mass conservation equations (N r b , N χ b and N ξ b , respectively) can be calculated as: where r b,p ¼ 1 þ χ p (1 À p p ) is the relative density of the bed layer for the pth sediment size-class and ξ p ¼ 1=(1 À p p ), p p being the porosity of the pth sediment size-class (Wu ). The terms D p and E p are the deposition and entrainment vertical fluxes, respectively, for the pth sediment sizeclass, and F b,p is the relative fraction of the pth sediment size-class at the bed layer. The size-specific deposition and entrainment fluxes can be expressed, respectively, as a function of the mixture depth-averaged volumetric concentration ϕ p and the capacity volumetric concentration ϕ Ã p for each sediment size-class: d p and Re p being the characteristic sediment diameter and particle Reynolds number for the pth sediment size-class, =ϕ 0 and Re 0 the median diameter of the suspended sediment and the corresponding particle Reynolds number, and ν the kinematic viscosity of water. The Zhang & Xie () formula has been used to calculate the capacity suspended concentration ϕ Ã p for each sediment size-class: This set of equations modelling variable-density flows over erodible bed will be referred to as the VD model. It is a generalized formulation which can be simplified to model the situation of reduced complexity.  (2) is integrated in each cell Ω i using the Gauss theorem:

XA-ROE SOLVER FOR 2D SEDIMENT-LADEN FLOWS
vector. Assuming a constant piecewise representation of the conserved variable U at the cell Ω i for the time t ¼ t n : where A i is the cell area, (18) can be expressed as: NE being the number of edges for the i cell, (En) k the value of the normal flux through each edge, l k the length of the edge and Ω i,k the area of the i cell associated with the kth edge.
Assuming a first-order reconstruction approach for the conserved variables in (20), the theory of Riemann problems (RPs) can be used to solve the 2D problem. For each kth cell edge, it is possible to define a local 1D RP along the direction normal to the edge including the momentum source terms: x being the coordinate normal to the kth cell edge. The solution of (21) provides the variation of the conserved variables in time and space, and U(x, t), in time and space.
Note that the net exchange term between the mixture flow and the bed, E b , has been dropped in the local RP (21), since it implies a mass source rather than a momentum source and will be incorporated into the solution as a cellcentred contribution.
The augmented value of the fluxes through the kth cell edge incorporates the non-conservative contribution of the momentum source term S into the convective fluxes (En) k .
In order to integrate (21) over a suitable control volume Δx=2], the momentum source terms are involved in the Riemann solver as a singular source at the discontinuityx ¼ 0 and linearized in time so that: where Δt ¼ t nþ1 À t n the time step and S ∨ n is a suitable numerical momentum source vector along the normal direction to the cell edge, which can be expressed in the 2D framework (x, y) as: H andT being the bed slope and friction momentum contributions, respectively, spatially integrated in the control volume corresponding to each cell edge. hence can be expressed as a plane Riemann problem in the local framework (x,ŷ), corresponding to normal and tangential directions to each cell edge, respectively. Defining a rotation matrix T, with an inverse matrix T À1 , as: which satisfy the condition: The local RP (21) at each cell edge can be projected in The new set of projected conserved variablesÛ is defined as: whereû ¼ un x þ vn y andv ¼ Àun y þ vn x are the flow velocities along thex andŷ coordinates, respectively. The projected convective fluxes F(Û) can be expressed as: Following Murillo & García-Navarro () and Murillo & García-Navarro (), the projected source vectorŜ, including both the bed slope and friction terms, can be expressed in the local framework (x,ŷ) as: and, applying the procedure used in (22) leading to express the integrated momentum source vectors on the local framework (x,ŷ) as: For the left-hand side of system (26), it is possible to define a singular Jacobian matrix M(Û) ¼ @F(Û) @Û as follows: Note that the solid lines in M(Û) (32) separate the terms associated with the bed evolution equation from the terms linked to the mixture flow equations. Furthermore, M(Û) has five real eigenvalues: where λ 1,...,4 corresponds to the mixture flow-wave structure and λ 5 relates to the contact wave associated with the bed evolution equation.
Using the rotation invariance property and the solution of the projected plane RP (26), it is possible to rewrite (20) as: where F(Û) k and (Ŝ ∨ n ) k are the numerical flux and the integrated momentum source contribution, respectively, resulting from solving the projected plane RP (26) for the kth edge in the local framework (x,ŷ).
Furthermore, the temporal and spatial integration of the net exchange flux term are approximated by: (19) and (35), it is possible to express (34) as: defining an augmented numerical flux F(Û) ↓ k for the kth cell edge which incorporates the integrated momentum source term into the convective numerical fluxes at each cell edge, ensuring the well-balanced property of steady states (Murillo & García-Navarro ): The resolution procedure is divided into two steps: first, where δÛ k ¼Û j ÀÛ i is the conserved variable jump andJ k is defined as: The approximate Jacobian matrix for the mixture flowJ k is diagonalizable with four approximate real eigenvalues: where the averaged celerityc k is defined as: Therefore, using the associated orthogonal basis of eigenvectorsẽ m,k , a matrix e P k ¼ (ẽ 1 ,ẽ 2 ,ẽ 3 ,ẽ 4 ) k can be built as: with the following property: k being the inverse matrix ofP k andΛ k being a diagonal matrix with approximate eigenvalues in the main diagonal.
One result of Roe's linearization is that the approximate Riemann solution consists of only discontinuities andÛ(x, t) is constructed as a sum of jumps or shocks. The solution for U(x, t) is governed by the celerities inΛ k and consists of four regions connected by five waves, one of them a contact wave with null celerity accounting for the integrated source term.
According to the Godunov-type method, it is sufficient to provide the solution forÛ(x, t) at the intercell position x ¼ 0 in order to obtain the updating numerical fluxes F ↓ k (37). The numerical flux at the left and right side of the kth cell edge can be estimated using an approximate flux function F(x, t): Following Toro (), the conserved variable differences δÛ k and source term spatial integral (Ŝ ∨ n ) k at the intercell edge are projected on the eigenvector basis in order to obtain the wave and source strength vectors,Ã k andB k , respectively.
The reconstruction of the approximated solution at the left and right sides of the intercell edge,Û À i andÛ þ j , respectively, can be expressed as: The relation between the intercell approximate fluxes F ↓À where the subscripts mÀ and mþ under the sums indicate waves travelling inward and outward from the i cell. It is worth mentioning that, due to the presence of the momentum source terms, it is no longer possible to define a general intercell flux function contrary to the homogeneous case. The corresponding intercell flux jump for the approximate solution is given by: The correct integration of the momentum source term and the global time step Δt is limited using the CFL condition as: with CFL <1. The integrated momentum source term (Ŝ ∨ n ) k (31) is projected onto the eigenvectors basis, separating the contribution associated with the bed slope and corresponding to the friction loss. Therefore, the total source strengths can be expressed asB k ¼B b,k þB τ,k , where:

Momentum source term integration and friction fix
(Ŝ ∨ b,n ) k and (Ŝ ∨ τ,n ) k being the integrated bed slope and friction terms, respectively, at the kth cell edge: whereH k ¼ (Àg ψrh δz b ) k is the integrated bed slope contribution andT k is the integrated shear force. Considering (52), the source strengths should agreeβ b1,k ¼ Àβ b3,k and β τ1,k ¼ Àβ τ3,k , and hence, the first component of the numerical flux vector, i.e. the mixture mass flux, including the homogeneous flux augmented with the source terms remains equal at both sides of the kth edge: Considering positive velocity (ũ k > 0) and subcritical flow,λ 1,k < 0 and the rest of the waves are positive (Figure 2).
The homogeneous mixture mass flux augmented with the bed slope term F 1↓ b,k can be computed as: F 1 k being the homogeneous numerical flux for the mixture mass conservation equation. Then, the integrated shear forceT k is defined considering the sign of the augmented mass flux F 1↓ b,k : whereτ b,x is the averaged normal shear stress at the cell edge, estimated using any resistance formulation (see Table 1), and d n,k is the distance between the centre of the neighbouring cells associated with the kth cell edge along the normal direction.
Therefore, the numerical flux for the mixture mass augmented with both the bed slope and the friction source terms F 1↓ k allows us to define a numerical fix to correct overestimated friction terms. The numerical flux is computed as: and the friction fix can be expressed as: In the case of supercritical regime, all the waves are positive and the augmented numerical mass flux at the cell edge does not depend on the value of the momentum source terms (Figure 3). However, in order to avoid overestimated source terms, it is possible to define an averaged intermediate state of the flux function F 1þ b;k at the right side of the edge, including the bed slope term as: F 1 þk being the averaged value of the homogeneous numerical mass flux at the right side of the cell edge. In this case, the integrated shear forceT k is defined considering the augmented mass flux F 1þ b;k : Therefore, the averaged intermediate numerical mass flux augmented with both the bed slope and the friction source terms F 1þ k at the right side allows us to define a numerical fix to correct overestimated friction terms. The numerical flux is computed as:  and the friction fix can be expressed as: The friction fix proposed can be straightforward extended to the case with negative velocity (ũ k < 0) at the kth cell edge.

Well-balanced equilibrium states
The At steady state, (66) is an ordinary differential equation which can be reordered as: In the generic solution of (67), density and flow depth are both variable. Nonetheless, the particular variabledepth solution with fixed density and free-surface level and the particular variable-density solution with fixed depth are interesting for numerical model validation and easy to compute exactly. Following Leighton et al. (), the bedlevel profile is defined as: A being the amplitude and L the length of the 1D channel.

The depth variable equilibrium equation leads to following
conditions for the suspended volumetric concentration and the flow depth: whereas the density variable solution can be expressed as:   it also shows much smaller dynamical time steps than the complete model proposed and probably leads to a loss of computational efficiency for real-scale and long-term morphodynamical computations.
Dambreak flow over non-uniform erodible beds The aim of this idealized test is to assess the capability of the proposed model to deal with non-uniform beds and to study the influence of the mixture composition on the hydrodynamic behaviour. The same large-scale and long-term dambreak test described in the above section is again used here but setting two different non-uniform sediment-size distributions in the bed layer. Both non-uniform beds consist of a mix of gravel (d p > 2 mm), medium sand (d p ¼ 1 mm) and fine material (d p ¼ 100 μm) weighted to maintain a constant medium diameter d m ¼ P p F b,p d p ¼ 4 mm, allowing us to compare the results with those obtained in the above section for uniform bed. Table 2 shows the sediment-size distribution for both non-uniform beds considered. Nonuniform bed A is composed mainly of medium gravel with small fractions of coarse sand and fines. In the non-uniform bed B, fine material and medium sand prevails over a small fraction of coarse gravel.
All the other parameters in simulations are set with the same values as in the above section.     Nevertheless, the bed-level evolution shows more marked differences than the free surface evolution, especially as the dambreak wave progresses downstream (see Figure 9). For the coarser non-uniform bed A, the bed level shows significant deviations with respect to the uniform bed case in the wavefront, an important deposition even appearing in the density-wave region at long-term stages. The finer non-uniform bed B presents a more marked erosion than the uniform bed configuration, especially at long-term stages as a consequence of the higher entrainment of fine material from the bed. This test aims to demonstrate the robustness of the proposed scheme with any kind of boundary shear stress formulation (Table 1) and the correct treatment of the wet-dry fronts. A 2D circular dambreak over a completely dry bed is considered.  Table 1) considering a non-cohesive solid phase. Furthermore, the (4) simplified Bingham friction formulation has also been tested, assuming that the solid phase is composed of cohesive materials. The parameters needed for each friction formulation are presented in Table 3. These residual velocities are caused by the imbalance between the convective fluxes and the overestimated friction terms. Furthermore, the proposed friction fix avoids the necessity of additional time step reductions to guarantee the stability of the scheme (see Figure 15). As was previously stated, these additional time step restrictions can lead to a marked increase of the computational time required by the numerical scheme.
Hence, the proposed source term integration procedure also ensures the model remains computationally efficient when complex friction terms are considered.  and robustness for all the cases tested.