## Abstract

In this work, a novel 2D depth-integrated numerical model for highly sediment-laden shallow flows over non-uniform erodible beds is presented, including variable density and exchange between the bed layer and the water–sediment mixture flow. 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 flow 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 flow and transported as suspended materials. The rheological behaviour of the flow is directly controlled by the properties of the mixture. A new x-split augmented Roe (xA-Roe) scheme is derived to solve the coupled flow and suspended solid-phase equations in both structured and unstructured meshes. The numerical scheme is defined to properly include density variations and momentum source terms, retaining a well-balanced flux 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 flows, demonstrating its accuracy, stability and robustness in the presence of complex bed topography, wetting–drying fronts and rapid morphological changes.

## HIGHLIGHTS

Numerical modelling of highly sediment-laden flows, such as mud or debris flows, over non-uniform erodible beds is a current challenge due to the complexity of the physical processes involved.

Coupling bewteen flow depth and mixture density requires robust, accurate and efficient numerical methods able to be applied to realistic flows.

Realistic highly sediment-laden flow models require incorporation of complex non-Newtonian friction formulations and robust treatments for the wet–dry fronts.

A new x-split augmented Roe (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.

The new xA-Row scheme is demonstrated to be robust, accurate and efficient for a wide range of unsteady problems involving variable-density flows 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 kinds of flows are usually classified as hyperconcentrated flows or mud/debris flows, mainly depending on the size of the solid particles being transported.

Hyperconcentrated flow lies between clear water and mud/debris flow. A clear water flow transitions into a hyperconcentrated flow when particles on the bed begin to move together and coarse sediment becomes suspended in the flow. The water–sediment flow begins to be affected by the suspended sediment when particle concentration reaches about 4% by volume (Pierson 2005). For higher concentrations, the mixture starts to show a non-Newtonian behaviour. Furthermore, fine sediment fractions have a greater impact on the mixture rheology than coarse sediments. A hyperconcentrated flow transitions into a mud/debris flow when rising concentrations of the sediment generate a critical yield stress in the fluid which allows coarse particles to be suspended indefinitely in the mixture flow (Calhoun & Clague 2018). Mud/debris flows are characterized by high sediment volume concentrations, often greater than 60%. In debris flows, sand/gravel and coarser sediment fractions predominate in the solid phase, whereas dominant fine fractions (silt and clay) are typical for mud flows. The sediment size distribution of the solid phase causes changes in the characteristics of the flow: mud flows show high Darcy numbers and reduced values of the modified Reynolds number, while debris flows are characterized by high values of modified Reynolds numbers and smaller Darcy numbers (Iverson 1997). Nevertheless, these transitional processes in the flow behaviour are extremely complex and continue to be debated. Furthermore, mathematical modelling of hyperconcentrated, mud and debris flows and their numerical resolution is still a challenging topic, especially when dealing with realistic applications.

Generally, natural hyperconcentrated and mud/debris flows consist of highly unsteady shallow flows running over non-uniform erodible beds, i.e. beds composed by different sediment size-classes. Most of the 1D and 2D numerical models recently reported for highly sediment-laden flows are based on quasi-single-phase assumption (Armanini *et al.* 2009; Cao *et al.* 2017; Xia *et al.* 2018), solving the mass and momentum equations for the mixture flow and the continuity equation for the suspended solid phase. Also, some two-phase models have been reported in recent years (Li *et al.* 2018; Greco *et al.* 2019), solving the mass and momentum equations for the liquid and solid phases separately. Although, theoretically, the two-phase mathematical model describes more properly the complex interaction between fluid and sediment particles, the high uncertainty involved in the equations and the difficulty to implement efficient and robust numerical schemes have hindered its application to realistic geophysical problems.

Usually, sediment-laden flows involve high levels of energy which are related to strong flow–bed interaction and important morphological changes at the bed. It has been demonstrated by large-scale laboratory experiments that mud/debris flows gain much of their mass and momentum as they flow over steep slopes as a consequence of the material entrainment from the erodible bed, before deposition begins on flatter terrain downstream (Iverson *et al.* 2010, 2011). For these kinds of highly unsteady flows, the coupling between the flow depth and the mixture density is one of the challenging key points in realistic environmental processes computation. The complexity of the numerical resolution and the computational cost of the solvers increase exponentially with the number of equations involved, and the coupling between flow variables adds special features to the mathematical model (Cao *et al.* 2017; Chen *et al.* 2018).

A novel 2D numerical model for depth-integrated hyperconcentrated and mud/debris shallow flows, including variable-density and non-uniform solid-phase transport, is presented in this work. The resulting system of equations is formed by the 2D conservation equations for the mass and momentum of the mixture, supplemented by the mass conservation equation for the different sediment size-classes suspended in the flow (Murillo *et al.* 2012) and the bed variation 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 mixture flow and transported as suspended materials. A pioneering first-order x-split Augmented Roe (xA-Roe) scheme is derived to solve the flow and the suspended solid-phase equations in both structured and unstructured two-dimensional meshes. The numerical scheme is defined to properly include the density variations and to ensure a well-balanced flux formulation in steady states. The global time step is dynamically controlled by the wave celerities of the coupled system of equations, preserving the scheme stability even for high flow–bed interaction and reducing the computational effort required by the model.

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 size-classes 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

*p*th sediment size-class, with and

*N*the number of sediment size-classes transported. The mixture density is given by , where 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 is the relative buoyant density of the

*p*th solid phase, being the density of the sediment particles.

*et al.*2013), as:where is the vector of conserved variables, and are the convective fluxes along the

*x*and

*y*global coordinates, respectively, is the momentum source term associated with the variation of the pressure force on the bottom, is the momentum dissipation due to the boundary shear stress between the mixture flow and the bed layer and accounts for the mass net exchange flux between the mixture flow and the bed layer:where

*h*represents the mixture flow depth, ( are the components of the depth-averaged flow velocity vector along the

*x*and

*y*coordinates, respectively, is the bed elevation and is the bed-normal projection of the gravity in the GCS,

*g*being the gravitational acceleration and the direction cosine of the bed normal with respect to the vertical axis (Juez

*et al.*2013).

*x*and

*y*coordinates are expressed as:and those of the frictional dissipation at the bottom boundary can be expressed as:being the components of the boundary shear stress vector between the flow and the bed. To date, there is not a universal closure relation for representing the shear stress in hyperconcentrated and mud/debris flows. The formulation selected to model the tangential forces generated by the boundary stresses incorporates the rheological behaviour of the water–sediment mixture in motion. Different kinds of shear stresses can determinate this complex rheology: turbulent stress , yield stress , viscous stress or Coulomb frictional stress .

- •
*n*being the Manning roughness coefficient.

- •
being the dynamic viscosity.

- •
- •
- •

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. (Murillo & García-Navarro 2011). The different friction formulations considered in this work have been summarized in Table 1.

. | Formulation . | Flow resistance relation . |
---|---|---|

1 | Pure turbulent | |

2 | Turbulent & Coulomb | |

3 | Turbulent & yield | |

4 | Simplified Bingham | |

5 | Full Bingham |

. | Formulation . | Flow resistance relation . |
---|---|---|

1 | Pure turbulent | |

2 | Turbulent & Coulomb | |

3 | Turbulent & yield | |

4 | Simplified Bingham | |

5 | Full Bingham |

*p*th sediment size-class and , being the porosity of the

*p*th sediment size-class (Wu 2007). The terms and are the deposition and entrainment vertical fluxes, respectively, for the

*p*th sediment size-class, and is the relative fraction of the

*p*th 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 and the capacity volumetric concentration for each sediment size-class: being an empirical parameter representing the difference between the near-bed sediment concentration and the depth-averaged sediment concentration for the

*p*th sediment class, the size-specific settling velocity of the sediment particles in clear water and the global volumetric concentration in the mixture. Water and suspended sediment are generally well mixed along the flow column in highly sediment-laden flows; hence, it is common to adopt . The parameters and are two empirical exponents accounting for the hindering effect on the settling velocity due to high suspended concentrations. Different empirical relationships can be found in the literature to estimate the hindering exponents and , the settling velocity , as well as the capacity volumetric concentration . For the sake of clarity, only the Zhang & Xie (1993) and Richardson & Zaki (1997) relations have been used for the estimation of the settling velocity and the empirical hindering exponents, respectively: and being the characteristic sediment diameter and particle Reynolds number for the

*p*th sediment size-class, and the median diameter of the suspended sediment and the corresponding particle Reynolds number, and the kinematic viscosity of water. The Zhang & Xie (1993) formula has been used to calculate the capacity suspended concentration 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.

## XA-ROE SOLVER FOR 2D SEDIMENT-LADEN FLOWS

*i*cell, the value of the normal flux through each edge, the length of the edge and the area of the

*i*cell associated with the

*k*th edge.

*k*th cell edge, it is possible to define a local 1D RP along the direction normal to the edge including the momentum source terms: being the coordinate normal to the

*k*th cell edge. The solution of (21) provides the variation of the conserved variables in time and space, and , in time and space. Note that the net exchange term between the mixture flow and the bed, , 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 cell-centred contribution.

*k*th cell edge incorporates the non-conservative contribution of the momentum source term into the convective fluxes . In order to integrate (21) over a suitable control volume , the momentum source terms are involved in the Riemann solver as a singular source at the discontinuity and linearized in time so that:where the time step and is a suitable numerical momentum source vector along the normal direction to the cell edge, which can be expressed in the 2D framework as: and being the bed slope and friction momentum contributions, respectively, spatially integrated in the control volume corresponding to each cell edge.

*et al.*2009) and hence can be expressed as a plane Riemann problem in the local framework , corresponding to normal and tangential directions to each cell edge, respectively. Defining a rotation matrix , with an inverse matrix , as:which satisfy the condition:

*k*th edge in the local framework .

*k*th 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 2010):

The resolution procedure is divided into two steps: first, the augmented fluxes at the intercell edges are determined using a new upwind augmented Roe scheme based on the x-split transformation (xA-Roe); second, the bed mass exchange source terms are later incorporated into the updated solution as a cell-centred contribution to the conserved variables.

### xA-Roe solver for the mixture flow equations

*k*th edge, separating the left

*i*cell and the right

*j*cell, is reduced to the conservation laws for the mixture flow and can be approximated by using the following constant coefficient linear RP (Toro 2009):where is a constant coefficient matrix yet to be defined. For the sake of clarity, the superscript will be removed from now on. Integrating (38) over a suitable control volume between and leads to the following constraint involving conservation across discontinuities:where is the conserved variable jump and is defined as:

One result of Roe's linearization is that the approximate Riemann solution consists of only discontinuities and is constructed as a sum of jumps or shocks. The solution for is governed by the celerities in and consists of four regions connected by five waves, one of them a contact wave with null celerity accounting for the integrated source term.

*k*th cell edge can be estimated using an approximate flux function : and being the value of the approximate intermediate states of the solution at the corresponding side of the

*k*th edge:

*k*th cell edge:where the subscripts and 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 for the equivalent 1D Riemann problem associated with the *k*th intercell edge ensures the well-balanced property of the xA-Roe scheme and avoids numerical oscillations in the solution when large momentum sources appear, especially associated with the friction term.

*k*intercell edge is the key point of the Godunov-type methods, the whole approximated solution participates in the updating of the values at the cells. In order to ensure the stability of the explicitly computed numerical solution, the time step should be small enough to avoid the interaction of waves from neighboring Riemann problems. The dynamical limitation of the time step at each

*k*edge is addressed here using the Courant–Friedrichs–Lewy (CFL) condition and assuming that the fastest wave celerity corresponds to the absolute maximum of the eigenvalues of the mixture flux Jacobian matrix (40). The limiting time steps at

*k*th edge are computed using:and the global time step is limited using the CFL condition as:with CFL .

### Momentum source term integration and friction fix

A correct integration of both the bed slope and the bed shear stress momentum source terms is a key point in augmented upwind schemes (Murillo & García-Navarro 2010) for shallow flows and ensures the proper balance between convective fluxes and source terms in steady-state cases. Furthermore, friction terms can generate numerical instabilities if they are not carefully integrated (Burguete *et al.* 2008) and have been reported to require additional time step restrictions over the classical CFL condition (Murillo *et al.* 2008). These additional time step restrictions can lead to a marked increase of the computational time required by the model. The consequence is a reduction of the efficiency, regardless of how the scheme is implemented (programming language, parallel computing, and available hardware). This drawback has been correctly addressed for constant density shallow flows (Murillo & García-Navarro 2011), but is still unresolved for flows including temporal and spatial density variations. This section is devoted to the integration procedure for the momentum source terms, paying special attention to the friction term in order to avoid additional time step restrictions.

*k*th cell edge:where is the integrated bed slope contribution and is the integrated shear force. Considering (52), the source strengths should agree and , 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

*k*th edge:

*k*th cell edge along the normal direction.

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 at the right side of the edge, including the bed slope term as:

The friction fix proposed can be straightforward extended to the case with negative velocity () at the *k*th cell edge.

## NUMERICAL RESULTS

### Well-balanced equilibrium states

*et al.*(2010) for ensuring the well-balanced character of variable-density shallow flow in the presence of bed-level variations. For a pure one-dimensional flow, under quiescent equilibrium (null velocity), frictionless conditions and null net exchange between the bed and the flow, the mixture mass and the solid-phase mass temporal variations reduce to zero, i.e. the equilibrium state must be maintained along time, and the 1D momentum equation becomes:

*et al.*(2010), the bed-level 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: and being the reference values for the mixture relative density and the flow depth, respectively.

A channel of is discretized using a 1D mesh of square cells with . Values , , and are set. The CFL is 0.5, and the final simulation time is . The exact depth-variable and density-variable solutions are imposed as initial conditions for the flow depth and the suspended concentrations.

Figures 4 and 5 show the comparison of the exact depth-variable and density-variable solutions with the corresponding numerical results at , respectively. The exact quiescent equilibrium is maintained along the time for both cases, demonstrating the well-balanced property of the proposed xA-Roe scheme for the simulation of density shallow flows involving topography variations.

### Large-scale and long-term dambreak

This idealized test was first proposed by Cao *et al.* (2004) and consists of a large-scale and log-term one-dimensional dambreak. Initially, the fluid is clear water. The channel length is set to , the dam is initially located at the middle of the channel, , and the initial water surface elevation is and at the left and right sides, respectively. The movable flat bed is made of an uniform non-cohesive sediment of diameter . The friction term is modelled using the turbulent Manning's relation and assuming a constant value of the roughness coefficient . The aim of this test is to study the influence of the suspended solid phase incorporated into the flow from the erodible bed in the dynamics of the dambreak wave in a relatively long channel and over a comparatively long period rather than under the typical laboratory scales. A 1D mesh of squared cells is used for the simulation with and a CFL = 0.5 is set.

Figure 6 (left) shows the flow free surface evolution and bed evolution at , and after the dambreak starting. For comparison with the proposed VD model, results obtained with a fixed-bed model (FB model) and with a passive-transport model (PT model) have also been plotted. The PT model is a particular case of the proposed sediment-laden model but setting the density of the mixture equals the water density regardless of the sediment concentration in the flow column. Sediments can be exchanged with the bed and transported as passive solutes, i.e. without influencing the hydrodynamics of the flow. This corresponds to (2) with . The FB model is an even simpler particular case where the net exchange bed/flow and suspended sediment transport do not exist. In order to allow suitable comparisons, the three models are solved with the same numerical scheme but incorporating the required restrictions to each model. This corresponds to (2) with , and .

The bed mobility considerably affects the free surface evolution compared with the fixed bed case. This can be significant for flooding prediction as the dambreak wavefront progresses faster when bed mobility is considered. Figure 6 (right) shows the volumetric concentration profiles of the solid phase in the flow at , and after the dambreak starting. As the dambreak wave progresses downstream, the mass exchange term with the bed incorporates into the flow a high quantity of sediment, leading to important depth-averaged sediment concentrations in the flow column. The solid-phase volumetric concentration shows a sharp increment at the dambreak wavefront, with values higher than compared with volumetric concentrations lower than at the central reach of the dambreak.

This suspended concentration change, hence a mixture density change, between the central reach and the wavefront leads to the appearance of an intermediate shock wave in the flow surface upstream the dambreak wavefront due to the incorporation of the mixture density into the mathematical model. However, mathematical models, which incorporate the suspended sediment conservation and bed evolution equations but do not take into account density changes, are not able to predict this intermediate shock wave. Incorporating the bed mobility but ignoring the influence of flow properties change due to the sediment entrainment (PT approach) leads to an overestimation of the wavefront velocity propagation (Figure 7 (left)). Moreover, the flow-level increment caused by the intermediate density shock wave is also unpredicted by this passive transport model, with an important drawback for the hazards determination against severe flooding (see Figure 7 (right)).

Finally, Figure 8 depicts the dynamic computational time-step evolution along the simulation for the FB model, the PT model and the proposed VD model. The time step is closely related to the computational time required by the model to perform simulations and hence to its efficiency. For the first stages after the dambreak starting, the FB model shows higher wavefront propagation velocities, leading to smaller time steps, than the models considering bed mobility. However, for the long-term stages, the models incorporating the bed evolution into the equations requires smaller time steps to ensure the computational stability due to the higher wavefront velocity. Furthermore, as the PT model overestimates the progression of the wavefront, 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 (), medium sand () and fine material () weighted to maintain a constant medium diameter , 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. Non-uniform 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.

. | Bed A . | Bed B . |
---|---|---|

Fines (Fraction 1) | ||

Sand (Fraction 2) | ||

Gravel (Fraction 3) | ||

Medium diameter |

. | Bed A . | Bed B . |
---|---|---|

Fines (Fraction 1) | ||

Sand (Fraction 2) | ||

Gravel (Fraction 3) | ||

Medium diameter |

All the other parameters in simulations are set with the same values as in the above section.

Figure 9 (left) shows the flow free surface evolution and bed evolution at , and after the dambreak starting, with the non-uniform beds A and B. The results obtained in the above section with uniform bed configuration are also depicted as references. The non-uniform composition of the bed slightly affects the free surface evolution compared with the uniform bed. The more marked differences are detected in the density-wave region at long-term stages for the non-uniform bed B. The entrainment of fine materials from the bed to the flow makes the jump in density accompanying the wavefront smoother, which finally causes the appearance of the density wave. Figure 9 (right) shows the total volumetric concentration of the solid phase in the flow. The fine material is incorporated into the mixture in the upstream region of the dambreak wave, where the erosive flow energy is lower, and hence, it reduces the density difference between the wavefront and the upstream region. Furthermore, the density peak associated with the dambreak wavefront is reduced as the presence of fine materials in the bed increases. The main consequence is that the density-wave becomes smeared.

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.

Figure 10 depicts the volumetric concentration of each sediment size-class in the mixture flow for both non-uniform bed configurations at , and after the dambreak starting. The free surface level has also been plotted for each case. On one hand, the coarser non-uniform bed A shows higher concentrations of the gravel fraction at the wavefront during all the stages of the dambreak flow, which causes the density jump between the wavefront and the upstream region. On the other hand, the finer non-uniform bed configuration B presents a higher volumetric concentration of sand material at the wavefront during the early stages of the dambreak flow, whereas the fine material fraction prevails along the whole dambreak flow at the long-term stages. Moreover, the volumetric concentration of the finer fraction shows a progressive transition between the upstream region and the wavefront region, avoiding the appearance of the marked density jump detected in the uniform bed case.

### Circular mud dambreak over positive and negative slopes

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. The domain is and initially, a sediment–water mixture column is considered in the centre of the domain, with height over the bed at and radius . Two different cases are tested with uniform bed slopes and , respectively. An initial volumetric concentration of an uniform sediment (grain diameter , porosity and ) is considered in the mixture column, leading to bulk density . The domain is discretized using an unstructured triangular mesh of cells, the CFL is set at 0.5 and the final simulation time is .

The (1) pure turbulent and (2) turbulent & Coulomb friction formulations have been used in this test for the determination of the flow resistance term (see 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.

. | (1) Pure turbulent . | (2) Turbulent & Coulomb . | (4) Simplified Bingham . |
---|---|---|---|

Manning coefficient | – | ||

Internal stability angle | – | – | |

Plastic viscosity | – | – | |

Yield stress | – | – |

. | (1) Pure turbulent . | (2) Turbulent & Coulomb . | (4) Simplified Bingham . |
---|---|---|---|

Manning coefficient | – | ||

Internal stability angle | – | – | |

Plastic viscosity | – | – | |

Yield stress | – | – |

Figure 11 depicts the depth distribution after flow stops for the positive slope case (left) and the negative slope case (right). The flow stops between and for the turbulent & Coulomb and simplified Bingham resistance formulations in both cases, whereas for the pure turbulent friction the flow does not stop and reaches the domain boundaries at for both cases. Therefore, the depth distribution as shown in Figure 11 for the pure turbulent formulation corresponds to . Regardless of the flow resistance formulation chosen, the scheme is able to maintain a concentric depth distribution, providing a correct treatment of the wet–dry dambreak wavefront without numerical instabilities.

The positive slope case shows a larger advance of the wavefront before the flow stops for the positive slope case than for the negative slope case , regardless of the friction formulation. Figure 12 (left) plots the final flow depth radial profiles for the simulation assuming the turbulent & Coulomb and the simplified Bingham friction formulations, whereas for the pure turbulent friction case the depth profile corresponds to the time after the initial time. Differences in the depth distribution arise directly from the flow resistance formulation considered. Furthermore, Figure 12 (right) shows the accumulated bed exchange mass between the flow and the bed, being deposition negative. The accumulated net exchange for the pure turbulent case is only representative at times before the flow reaches the domain boundaries . The accumulated deposition mass depends inversely on the flow velocity and directly on the wetted area. The pure turbulent case shows a lower deposited mass due to the higher velocity of the dambreak wavefront. Furthermore, the turbulent & Coulomb friction case generates higher deposition rates than the simplified Bingham case since the wetted area after the flow stops is slightly larger.

The mass exchange between the bed and the flow generates a change in the solid-phase mass of the flow and hence changes in the mixture density. Figure 13 depicts the flow density distribution at with both the pure turbulent and the turbulent & Coulomb friction formulations. Densities remain higher for the turbulent case than for the turbulent & Coulomb case as the suspended solid-phase loss due to bed exchange is lower.

Finally, in order to demonstrate the performance of the friction term integration proposed in the section “Momentum source term integration and friction fix”, the results obtained with the turbulent & Coulomb resistance formulation with and without including the fix proposed for the friction term are compared. Figure 14 (left) shows the flow surface level along the *x*-axis at time , whereas Figure 14 (right) depicts the *u*-velocity *x*-axis profile at the same time. When the proposed friction fix is applied, the flow reaches the quiescent equilibrium state (null velocity) with free surface slope angles closed to the internal stability angle of the mixture.

Furthermore, if the friction fix is not considered, residual velocities remain even after the wavefront of the dambreak stops, leading to transient oscillations of the flow free surface. 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.

## CONCLUSIONS

In the proposed 2D model for highly sediment-laden unsteady flows over non-uniform erodible beds, two main novelties have been presented. First, the mathematical 2D system of conservation equations includes a new description of the coupling between the mass and momentum of the water–sediment mixture flow and the mass conservation equation for the different sediment size-classes suspended in the flow. Second, the 2D system is solved using a pioneering first-order xA-Roe scheme for sediment-laden flows, which simplifies considerably the correct integration of the bed slope and friction source terms. A new integration procedure is proposed to avoid numerical drawbacks and loss of efficiency when complex non-Newtonian friction terms are considered.

The scheme is able to deal with structured and unstructured square and triangular two-dimensional meshes. The proposed resolution strategy ensures the well-balanced property for steady states and the correct estimation of the momentum source terms at wet–dry fronts. The stability of the scheme is ensured by a CFL condition based on the eigenvalues of the Jacobian matrix of the coupled system of equations. Furthermore, the procedure proposed for integrating the momentum source terms avoids the necessity of additional time step restrictions to maintain the stability of the solution, even when complex friction terms or wet–dry conditions are considered. Finally, the proposed numerical model has been evaluated with four benchmarking tests, including steady-state cases with exact solution and idealized highly transient cases with wet–dry fronts. The proposed scheme has demonstrated its accuracy, efficiency and robustness for all the cases tested.

## ACKNOWLEDGEMENTS

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.