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

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, represents the scalar depth-averaged volumetric concentration of pth 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:
formula
(1)
where is the global modified sediment concentration in the mixture and is the relative buoyant density of the pth solid phase, being the density of the sediment particles.
The relevant formulation for the two-dimensional sediment-laden flow model (Figure 1) includes the depth-averaged 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. 2013), as:
formula
(2)
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:
formula
(3)
formula
(4)
formula
(5)
formula
(6)
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).
Figure 1

One-dimensional sketch of the depth-averaged components involved in sediment-laden flows over erodible bed.

Figure 1

One-dimensional sketch of the depth-averaged components involved in sediment-laden flows over erodible bed.

The components of the bed slope source term vector (5) along the x and y coordinates are expressed as:
formula
(7)
and those of the frictional dissipation at the bottom boundary can be expressed as:
formula
(8)
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 .
  • The turbulent effects near the bed can be expressed as a function of the Manning's roughness coefficient :
    formula
    (9)
    • n being the Manning roughness coefficient.

  • In case of a pure Newtonian fluid, the viscous stress can be estimated as:
    formula
    (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 :
    formula
    (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:
    formula
    (12)
  • If the ratio is smaller than 0.5, the full Bingham relation (12) can be reduced to:
    formula
    (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.

Table 1

Flow resistance formulations

FormulationFlow resistance relation
Pure turbulent  
Turbulent & Coulomb  
Turbulent & yield  
Simplified Bingham  
Full Bingham  
FormulationFlow resistance relation
Pure turbulent  
Turbulent & Coulomb  
Turbulent & yield  
Simplified Bingham  
Full Bingham  
In the source term (6), the global net exchange fluxes for the mixture, solid phase and bed layer mass conservation equations (, and , respectively) can be calculated as:
formula
(14)
where is the relative density of the bed layer for the pth sediment size-class and , being the porosity of the pth sediment size-class (Wu 2007). The terms and are the deposition and entrainment vertical fluxes, respectively, for the pth sediment size-class, and 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 and the capacity volumetric concentration for each sediment size-class:
formula
(15)
being an empirical parameter representing the difference between the near-bed sediment concentration and the depth-averaged sediment concentration for the pth 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:
formula
(16)
and being the characteristic sediment diameter and particle Reynolds number for the pth 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:
formula
(17)

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

This section is devoted to the derivation of a new numerical scheme for 2D sediment-laden flows considering the net exchange mass flux between the mixture flow and the bed layer. System (2) is time-dependent, nonlinear and contains source terms. Under the hypothesis of dominant advection, it can be classified and numerically dealt with as belonging to the family of hyperbolic systems. In order to obtain a numerical solution for the five equations, the spatial domain is divided in computational cells using a mesh fixed in time and system (2) is integrated in each cell using the Gauss theorem:
formula
(18)
being the flux normal to the cell boundary, the outward unit normal vector and the momentum source term vector. Assuming a constant piecewise representation of the conserved variable at the cell for the time :
formula
(19)
where is the cell area, (18) can be expressed as:
formula
(20)
being the number of edges for the 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 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:
formula
(21)
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 , 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.
The augmented value of the fluxes through the kth 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:
formula
(22)
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:
formula
(23)
and being the bed slope and friction momentum contributions, respectively, spatially integrated in the control volume corresponding to each cell edge.
Furthermore, (21) satisfies the rotation invariance property (Godlewski & Raviart 1996; Castro 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:
formula
(24)
which satisfy the condition:
formula
(25)
The local RP (21) at each cell edge can be projected in the local framework as (Murillo & García-Navarro 2012; Murillo & Navas-Montilla 2016):
formula
(26)
The new set of projected conserved variables is defined as:
formula
(27)
where and are the flow velocities along the and coordinates, respectively. The projected convective fluxes can be expressed as:
formula
(28)
Following Murillo & García-Navarro (2012) and Murillo & García-Navarro (2016), the projected source vector , including both the bed slope and friction terms, can be expressed in the local framework as:
formula
(29)
and, applying the procedure used in (22), it is integrated in the control volume corresponding to each cell edge as:
formula
(30)
leading to express the integrated momentum source vectors on the local framework as:
formula
(31)
For the left-hand side of system (26), it is possible to define a singular Jacobian matrix as follows:
formula
(32)
Note that the solid lines in (32) separate the terms associated with the bed evolution equation from the terms linked to the mixture flow equations. Furthermore, has five real eigenvalues:
formula
(33)
where corresponds to the mixture flow-wave structure and 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:
formula
(34)
where and 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 .
Furthermore, the temporal and spatial integration of the net exchange flux term are approximated by:
formula
(35)
being . Using (19) and (35), it is possible to express (34) as:
formula
(36)
defining an augmented numerical flux 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 2010):
formula
(37)

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

The numerical fluxes at the cell edges for the mixture mass and momentum conservation, , are computed by means of an augmented Roe (A-Roe) Riemann solver. The numerical flux for the bed evolution equation is null, and hence, the bed evolution equation can be discarded for the determination of the mass and momentum fluxes of the mixture at the intercell edges. Therefore, the plane Riemann problem (26) projected in the local framework of the kth 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):
formula
(38)
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:
formula
(39)
where is the conserved variable jump and is defined as:
formula
(40)
formula
(41)
The approximate Jacobian matrix for the mixture flow is diagonalizable with four approximate real eigenvalues:
formula
(42)
where the averaged celerity is defined as:
formula
(43)
being .
Therefore, using the associated orthogonal basis of eigenvectors , a matrix can be built as:
formula
(44)
with the following property:
formula
(45)
being the inverse matrix of and 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 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.

According to the Godunov-type method, it is sufficient to provide the solution for at the intercell position in order to obtain the updating numerical fluxes (37). The numerical flux at the left and right side of the kth cell edge can be estimated using an approximate flux function :
formula
(46)
and being the value of the approximate intermediate states of the solution at the corresponding side of the kth edge:
formula
(47)
Following Toro (1997), the conserved variable differences and source term spatial integral at the intercell edge are projected on the eigenvector basis in order to obtain the wave and source strength vectors, and , respectively.
formula
(48)
The reconstruction of the approximated solution at the left and right sides of the intercell edge, and , respectively, can be expressed as:
formula
(49)
where , and the subscripts and under the sums indicate waves travelling inward and outward from the i cell.
Being the solution defined as a sum of jumps or shocks between the different intermediate states, the solution for the approximate flux function involves the initial unaltered fluxes at the left and right cells, and , and three intermediate states. Each inner constant state involves an intermediate constant flux function. The approximate flux function provides the intercell fluxes at the left and right sides of the initial discontinuity at , labelled as and (46), with:
formula
(50)
The relation between the intercell approximate fluxes and can be analysed using the Rankine–Hugoniot (RH) relation at , which includes a steady contact wave (Li & Chen 2006; Rosatti & Begnudelli 2010) between approximate solutions and . Following the linear case, the approximate solution for the fluxes can be constructed defining the appropriate RH condition across each moving wave. The telescopic property of the linear solutions for the approximate flux function provides the definition of the fluxes at the left and right sides of the kth cell edge:
formula
(51)
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:
formula
(52)

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.

Although the computation of the approximate numerical fluxes at each 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 kth edge are computed using:
formula
(53)
and the global time step is limited using the CFL condition as:
formula
(54)
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.

The integrated momentum source term (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 as , where:
formula
(55)
and being the integrated bed slope and friction terms, respectively, at the kth cell edge:
formula
(56)
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 kth edge:
formula
(57)
Considering positive velocity () and subcritical flow, and the rest of the waves are positive (Figure 2). The homogeneous mixture mass flux augmented with the bed slope term can be computed as:
formula
(58)
being the homogeneous numerical flux for the mixture mass conservation equation. Then, the integrated shear force is defined considering the sign of the augmented mass flux :
formula
(59)
where is the averaged normal shear stress at the cell edge, estimated using any resistance formulation (see Table 1), and is the distance between the centre of the neighbouring cells associated with the kth cell edge along the normal direction.
Figure 2

Approximate mass flux function for the subcritical regime. The blue region indicates the intermediate value of the mass flux considered for the proposed friction fix. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.027.

Figure 2

Approximate mass flux function for the subcritical regime. The blue region indicates the intermediate value of the mass flux considered for the proposed friction fix. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.027.

Therefore, the numerical flux for the mixture mass augmented with both the bed slope and the friction source terms allows us to define a numerical fix to correct overestimated friction terms. The numerical flux is computed as:
formula
(60)
and the friction fix can be expressed as:
formula
(61)

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:

Figure 3

Approximate mass flux function for the supercritical regime. The blue region indicates the intermediate value of the mass flux considered for the proposed friction fix. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.027.

Figure 3

Approximate mass flux function for the supercritical regime. The blue region indicates the intermediate value of the mass flux considered for the proposed friction fix. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.027.

formula
(62)
being the averaged value of the homogeneous numerical mass flux at the right side of the cell edge. In this case, the integrated shear force is defined considering the augmented mass flux :
formula
(63)
Therefore, the averaged intermediate numerical mass flux augmented with both the bed slope and the friction source terms at the right side allows us to define a numerical fix to correct overestimated friction terms. The numerical flux is computed as:
formula
(64)
and the friction fix can be expressed as:
formula
(65)

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

NUMERICAL RESULTS

Well-balanced equilibrium states

The idealized case with the exact solution was initially proposed by Leighton 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:
formula
(66)
At steady state, (66) is an ordinary differential equation which can be reordered as:
formula
(67)
In the generic solution of (67), density and flow depth are both variable. Nonetheless, the particular variable-depth 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. (2010), the bed-level profile is defined as:
formula
(68)
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:
formula
(69)
whereas the density variable solution can be expressed as:
formula
(70)
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.

Figure 4

Exact VD quiescent equilibrium and simulation results at: (left) flow surface-level FSL and (right) suspended volumetric concentration.

Figure 4

Exact VD quiescent equilibrium and simulation results at: (left) flow surface-level FSL and (right) suspended volumetric concentration.

Figure 5

Exact VD quiescent equilibrium and simulation results at : (left) flow surface level FSL and (right) suspended volumetric concentration.

Figure 5

Exact VD quiescent equilibrium and simulation results at : (left) flow surface level FSL and (right) suspended volumetric concentration.

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 .

Figure 6

Dambreak long-term hydraulics over the mobile flat bed: (left) flow free surface and bed surface and (right) depth-averaged solid-phase concentration in the flow. Front top to bottom, , and .

Figure 6

Dambreak long-term hydraulics over the mobile flat bed: (left) flow free surface and bed surface and (right) depth-averaged solid-phase concentration in the flow. Front top to bottom, , 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)).

Figure 7

Temporal evolution of the (left) wavefront position and (right) flow free surface elevation at.

Figure 7

Temporal evolution of the (left) wavefront position and (right) flow free surface elevation at.

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.

Figure 8

Time step evolution.

Figure 8

Time step evolution.

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.

Table 2

Sediment-size distribution for non-uniform beds

Bed ABed B
Fines (Fraction 1)   
Sand (Fraction 2)   
Gravel (Fraction 3)   
Medium diameter   
Bed ABed 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.

Figure 9

Dambreak hydraulics over non-uniform bed: (left) flow free surface and bed surface and (right) total solid concentration. Front top to bottom, , and.

Figure 9

Dambreak hydraulics over non-uniform bed: (left) flow free surface and bed surface and (right) total solid concentration. Front top to bottom, , and.

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.

Figure 10

Volumetric concentration of each sediment size-class in the flow: (left) non-uniform bed A and (right) non-uniform bed B. Front top to bottom, , and.

Figure 10

Volumetric concentration of each sediment size-class in the flow: (left) non-uniform bed A and (right) non-uniform bed B. Front top to bottom, , and.

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.

Table 3

Parameters for the friction relationships used in the test “Circular mud dambreak over positive and negative slopes”

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

Figure 11

Flow depth distribution at for (left) positive bed slope and (right) negative bed slope. From top to bottom: (1) pure turbulent, (2) turbulent & Coulomb and (3) simplified Bingham friction formulations.

Figure 11

Flow depth distribution at for (left) positive bed slope and (right) negative bed slope. From top to bottom: (1) pure turbulent, (2) turbulent & Coulomb and (3) simplified Bingham friction formulations.

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.

Figure 12

(Left) Concentric depth distribution profiles after the flow stop and (right) accumulated deposition mass for all the cases tested.

Figure 12

(Left) Concentric depth distribution profiles after the flow stop and (right) accumulated deposition mass for all the cases tested.

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.

Figure 13

Flow density distribution at with (left) pure turbulent and (right) turbulent & Coulomb resistance formulations.

Figure 13

Flow density distribution at with (left) pure turbulent and (right) turbulent & Coulomb resistance formulations.

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.

Figure 14

X-axis profile for (left) the flow level and (right) the -velocity at time with and without the proposed friction fix. Top: positive slope case and bottom: negative slope case.

Figure 14

X-axis profile for (left) the flow level and (right) the -velocity at time with and without the proposed friction fix. Top: positive slope case and bottom: negative slope case.

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.

Figure 15

Time step evolution with and without friction fix: (left) positive slope case and (right) negative slope case.

Figure 15

Time step evolution with and without friction fix: (left) positive slope case and (right) negative slope case.

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.

REFERENCES

REFERENCES
Armanini
A.
Fraccarollo
L.
Rosatti
G.
2009
Two-dimensional simulation of debris flows in erodible channels
.
Computers & Geosciences
35
(
5
),
993
1006
.
Calhoun
N. C.
Clague
J. J.
2018
Distinguishing between debris flows and hyperconcentrated flows: an example from the eastern Swiss Alps
.
Earth Surface Processes and Landforms
43
(
6
),
1280
1294
.
Cao
Z.
Pender
G.
Wallis
S.
Carling
P.
2004
Computational dam-break hydraulics over erodible sediment bed
.
Journal of Hydraulic Engineering
130
(
7
),
689
703
.
Cao
Z.
Xia
C.
Pender
G.
Liu
Q.
2017
Shallow water hydro-sediment-morphodynamic equations for fluvial processes
.
Journal of Hydraulic Engineering
143
(
5
),
02517001
.
Castro
M.
Fernández-Nieto
E.
Ferreiro
A.
García-Rodríguez
J.
Parés
C.
2009
High order extensions of roe schemes for two-dimensional nonconservative hyperbolic systems
.
Journal of Scientific Computing
39
,
67
114
.
Chen
C.-H.
Lin
Y.-T.
Chung
H.-R.
Hsieh
T.-Y.
Yang
J.-C.
Lu
J.-Y.
2018
Modelling of hyperconcentrated flow in steep-sloped channels
.
Journal of Hydraulic Research
56
(
3
),
380
398
.
Godlewski
E.
Raviart
P.-A.
1996
Numerical Approximation of Hyperbolic Systems of Conservation Laws
.
Springer-Verlag
,
New York
.
Greco
M.
Di Cristo
C.
Iervolino
M.
Vacca
A.
2019
Numerical simulation of mud-flows impacting structures
.
Journal of Mountain Science
16
(
2
),
364
382
.
Iverson
R. M.
1997
The physics of debris flows
.
Reviews of Geophysics
35
(
3
),
245
296
.
Iverson
R. M.
Logan
M.
LaHusen
R. G.
Berti
M.
2010
The perfect debris flow? Aggregated results from 28 large-scale experiments
.
Journal of Geophysical Research: Earth Surface
115
,
F03005
.
Iverson
R. M.
Reid
M. E.
Logan
M.
LaHusen
R. G.
Godt
J. W.
Griswold
J. P.
2011
Positive feedback and momentum growth during debris-flow entrainment of wet bed sediment
.
Nature Geoscience
4
,
116
121
.
Juez
C.
Murillo
J.
García-Navarro
P.
2013
2D simulation of granular flow over irregular steep slopes using global and local coordinates
.
Journal of Computational Physics
255
,
166
204
.
Leighton
F. Z.
Borthwick
A. G. L.
Taylor
P. H.
2010
1-D numerical modelling of shallow flows with variable horizontal density
.
International Journal for Numerical Methods in Fluids
62
(
11
),
1209
1231
.
Li
J.
Chen
G.
2006
The generalized Riemann problem method for the shallow water equations with bottom topography
.
International Journal for Numerical Methods in Engineering
65
,
834
862
.
Li
J.
Cao
Z.
Hu
K.
Pender
G.
Liu
Q.
2018
A depth-averaged two-phase model for debris flows over erodible beds
.
Earth Surface Processes and Landforms
43
(
4
),
817
839
.
Murillo
J.
García-Navarro
P.
Burguete
J.
2008
Time step restrictions for well balanced shallow water solutions in non-zero velocity steady states
.
International Journal for Numerical Methods in Fluids
56
,
661
686
.
Murillo
J.
Latorre
B.
García-Navarro
P.
Riemann
A.
2012
A Riemann solver for unsteady computation of 2D shallow flows with variable density
.
Journal of Computational Physics
231
,
4775
4807
.
Pierson
T.
2005
Hyperconcentrated Flow – Transitional Process Between Water Flow and Debris Flow
.
Debris-Flow Hazards and Related Phenomena
.
Springer
,
Berlin
.
Richardson
J.
Zaki
W.
1997
Sedimentation and fluidisation: part I
.
Chemical Engineering Research and Design
75
,
82
100
.
Toro
E.
1997
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
.
Springer
,
Berlin
.
Toro
E. F.
2009
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
.
Springer-Verlag
,
Berlin
,
Heidelberg
.
Wu
W.
2007
Computational River Dynamics
.
NetLibrary, Inc., CRC Press
,
London, UK
.
Xia
C.-c.
Li
J.
Cao
Z.-x.
Liu
Q.-q.
Hu
K.-h.
2018
A quasi single-phase model for debris flows and its comparison with a two-phase model
.
Journal of Mountain Science
15
(
5
),
1071
1089
.
Zhang
R.
Xie
J.
1993
Sedimentation Research in China: Systematic Selections
.
China Water and Power Press
,
Beijing, China
.