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

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.

• •
  The turbulent effects near the bed can be expressed as a function of the Manning's roughness coefficient ⁠:

  

  (9)

• • n being the Manning roughness coefficient.

• •
  In case of a pure Newtonian fluid, the viscous stress can be estimated as:

  

  (10)

• • 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 ⁠:

  

  (11)
• •
  Non-Newtonian Bingham-type fluids do not flow until a threshold value of the tangential stress, the yield stress ⁠, is reached. During the movement, the boundary shear stress is characterized by means of a cubic equation accounting for the plastic viscosity of the sediment–water mixture:

  

  (12)
• •
  If the ratio is smaller than 0.5, the full Bingham relation (12) can be reduced to:

  

  (13)

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.

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.

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.

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.

The correct integration of the momentum source term for the equivalent 1D Riemann problem associated with the kth 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.

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.

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 kth cell edge.

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.

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.

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.

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.

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.

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.

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.

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.