Open access peer-reviewed chapter

Efficient Simulation Tools (EST) for Sediment-Laden Shallow Flows

Written By

Sergio Martínez-Aranda and Pilar García-Navarro

Reviewed: 24 August 2022 Published: 11 October 2022

DOI: 10.5772/intechopen.107345

From the Edited Volume

Modeling of Sediment Transport

Edited by Davide Pasquali

Chapter metrics overview

88 Chapter Downloads

View Full Metrics


Rapid flows of water-sediment mixtures are probably the most challenging and unknown geophysical gravity-driven processes. The fluidized material in motion consists of a mixture of water and multiple solid phases with different specific characteristics. Modeling sediment transport involves an increasing complexity due to the variable bulk properties in the sediment-water mixture, the coupling of physical processes, and the presence of multiple layer phenomena. Two-dimensional shallow-type mathematical models are built in the context of free surface flows and are applicable to most of these geophysical surface processes. Their numerical solution in the finite volume framework is governed by the dynamical properties of the equations, the coupling between flow variables and the computational grid. The complexity of the numerical resolution of these highly unsteady flows and the computational cost of simulation tools increase considerably with the refinement of the non-structured spatial discretization, so that the computational effort required is one of the biggest challenges for the application of depth-averaged 2D models to large-scale long-term flows. Throughout this chapter, the combination of 2D mathematical models, robust numerical methods, and efficient computing kernels is addressed to develop Efficient Simulation Tools (EST’s) for environmental surface processes involving sediment transport with realistic temporal and spatial scales.


  • sediment-laden flows
  • depth-averaged rheological models
  • bed-material entrainment
  • shallow flow models
  • GPU high-performance-computing

1. Introduction

Sediment transport is ubiquitous in environmental water bodies such as rivers, floods, coasts and estuaries, but also is the main process in natural ladslides, debris flows, muddy slurries or mining tailings (Figure 1). These are considered highly solid-laden fluids, where the density of the water-solid mixture can be more than twice or three times the water density and the bulk solid phase represents about 40–80% of the flow volume [1].

Figure 1.

Geophysical surface flows involving sediment transport.

The presence of the solid phases, especially the fine material as silt or clay, affects the rheological behavior of the mixture. The water-sediment mixture rheology begins to be affected by fine solid particle transported in the flow when the volumetric concentration of fine sediment particles reaches about 4% by volume [2], creating a slight shear strength within the fluid. For higher concentrations, the mixture shows a marked non-Newtonian rheology. Mud and debris flows lie between hyperconcentrated flows and wet avalanches [3]. High concentrations of solids generate a critical yield stress which allows that gravel particles can be suspended indefinitely into the fluidized material.

The mass exchange between mud/debris mixtures and erodible beds involves complicated physical processes and the understanding of its theoretical basis remains unclear. Experiments in large-scale channel [4, 5] and field observations in real debris events [6, 7] indicate that the entrainment volume in steep beds can be in the same order of magnitude as the initial volume mobilized. Debris and mud flows gain much of their mass and momentum as they move over steep slopes as a consequence of the material entrainment from the erodible bed, before deposition begins on flatter terrain downstream.

The mathematical modeling of solid-liquid mixture flows and their numerical resolution is still a challenging topic, especially when dealing with realistic applications. When liquid and solid phases are well-mixed, assuming that the solid phase is distributed uniformly over the flow column allows the use of depth-averaged models derived from the vertical integration of the Navier-Stokes equations [8]. Shallow-type mathematical models represent a simplified formulation, derived from the general 3D Navier-Stokes equations, which is applicable to a large number of these geophysical surface processes involving sediment transport. The simplest models, used in river and coastal dynamics, assume small enough sediment concentrations throughout the flow to consider the bulk density constant and uniform. Most of the numerical models reported for highly solid-laden flows also use this one-single-phase approach, neglecting the bulk density in the shallow-flow mass and momentum equations [9, 10]. Nevertheless, even small density gradients influence importantly the mixing dynamics in flow confluences [11] and larger gradients can also generate numerical oscillations and instabilities throughout mixing interfaces.

Modeling sediment transport involves an increasing complexity with respect to rigid-bed shallow water models [12] due to the presence of variable sediment-fluid mixture properties, coupling of physical processes and multiple layers phenomena [13]. One of the biggest challenges for the application of depth-averaged models to realistic large-scale long-term flows is the computational effort required. New strategies to reduce the computational effort have been developed in the last decade through the use of parallelization techniques based on Multiprocessing (OpenMP) or Message Passing Interface (MPI), which allow to run simulations on multi-CPU clusters. Their main drawback is the associated hardware cost and energy requirements, which are directly proportional to the number of CPU-cores available and limit their efficiency. In the last years, the usage of Graphics Processing Units (GPU) hardware accelerators for sequential computation has demonstrated to be an efficient and low cost alternative to the traditional multi-CPU strategies [14]. GPU-accelerated algorithms have been developed for real-time floods forecasting [15], real-scale bedload erosive shallow-flows [16] or tsunami prediction [17]. GPU devices are oriented to perform arithmetical operations on vector-based information. Unlike the conventional shared-memory multi-CPU implementations, the GPU solution must be designed taking into account the fact that the GPU is an independent device with its own RAM memory. This means that the memory transfer between the conventional RAM memory and the GPU device memory plays a key role in the performance of GPU-accelerated software.


2. Depth-integrated equations for shallow flows

2.1 Mass and linear momentum conservation

The flow of a water-sediment mixture can be mathematically described assuming the movement of the solid particles as a diffusion phenomenon into the liquid phase. Then, the continuity and momentum conservation for the mixture, supplemented with the transport equation for the solid phase, can be established for modeling these two-phase flows. Although both solid and liquid phases are incompressible when considered independently, the bulk behavior of the solid-liquid mixture is the same as that of a compressible material depending on the local solid phase volumetric concentration. It possible to define the bulk density ρ=ρwn+ρsϕ and linear momentum ρu=ρwnuw+ρsϕus of the mixture, being ρw and ρs the density of the fluid and solid phases, respectively, ϕ the volumetric solid-phase concentration and n=1ϕ the volumetric fluid-phase fraction or mixture porosity, uw the velocity of the pore-fluid and us the advective sediment particle velocity.

The 3D time-averaged Navier-Stokes equations for mass and momentum conservation of a two-phase mixture can be written in the Cartesian coordinate system X=xyz as


where u=uxuyuz is the bulk velocity in any point of the fluidized material, F=FxFyFz are the external forces, such as gravity and p denotes the pressure of the mixture. The term τ=τijij=xyz is the deviatoric stress tensor. With low solid-phase concentrations, the water-sediment mixture behaves as a Newtonian fluid, with a constitutive relation given by the Navier-Poisson law [8]. Nevertheless, for high sediment concentrations the mixture becomes a kind of non-Newtonian fluid with a complex constitutive law, which depends on multiple factors, relating stresses and deformation rates.

In order to develop a shallow-type depth-averaged mathematical model, the Navier-Stokes system (1) and (2) is integrated between the free surface zs=zstxy and the bottom surface of the flow column zb=zbtxy, which is also considered a movable interface. The kinematic conditions at these boundaries can be expressed as


being Nb the net volumetric flux through the bed interface along the zcoordinate, being positive for net entrainment conditions and negative for net deposition fluxes. The subscripts s and b indicate the value of the corresponding variable at the flow free surface and the bottom bed interface respectively.

The mass conservation Eq. (1) is integrated along the water column as


Applying Leibnitz’s rule to each term of (5) and using (3) and (4) leads to


where the depth-averaged of the mixture bulk density ρ¯ and velocities are defined as


being h=zszb the flow depth and ρb the mixture bulk density at the bed interface zb.

Regarding the momentum depth integration, the volumetric force components in (1) and (2) are considered for the sake of simplicity as Fx=Fy=0 and Fz=ρg. Moreover, the z-momentum equation in (2) can be simplified by assuming shallow-flow scaling and neglecting temporal, convective and stress terms, leading to a hydrostatic pressure distribution along the flow column.

The x-momentum in (2) is integrated throughout the flow depth in the following way


Applying Leibnitz’s rule and considering the kinematic boundary conditions at both the free surface (3) and the bed interface (4), the left hand side of (8) is integrated as


with DxxDxy accounting for the depth-averaged dispersion momentum transport due to the non-uniformity of the vertical velocity profile, and defined as


On the right hand side of (8), using the hydrostatic pressure distribution and assuming that the vertical density gradient is negligible compared with those along the horizontal plane, as occurs in natural debris and mud slurry flows, the integral of the pressure gradient along the xcoordinate can be expressed as

zbzsxpdz=gzbzsxzzsρdzdz=x12gρ¯h2Conservative termgρ¯hzbxBedpressure termE12

separating the pressure gradient term into a conservative component plus a bed-pressure component.

The stress terms on the right hand side of (8) are also integrated along the flow column, leading to the final expresion for the depth-integrated momentum equation along the xcoordinate


and, similarly, for the ycoordinate. The boundary stress terms in (13) at the free surface zs and the bed interface zb are


where τsx denotes the xcoordinate component of the wind action at the free surface, whereas τbx represents the xcoordinate component of the boundary shear stress at the bed interface, opposing to the flow movement. Furthermore, the depth-integrated stress terms along the flow column TxxTxy are defined as


2.2 Depth-integrated solid transport equation

The solid-phase transport process is governed by


where us=usxusyusz is the local velocity of the solid particles. Because the solid-phase velocity is not included in the dependent variables of the fluid dynamic system, and assuming that the sediment particles are incompressible and non-porous, (17) is rewritten as


where the term on the right hand side accounts for the drag effects caused by the liquid phase on the advective solid flux. Due to the definition of the bulk linear momentum of the mixture, this term is usually neglected and the approximation uus is actually assumed.

The transport equation (18) must also be integrated throughout the entire flow column defining the total solid phase being transported in the flow column as


where ϕ¯ is the dept-averaged volumetric solid concentration. Therefore, applying the Leibnitz’s integration rule to (18) and considering the kinematic boundary conditions at both the free surface (3) and the bed interface (4), the depth-integrated equation for the solid phase in the flow column reduces to


where ϕb denotes the solid concentration in the bed layer surface zb and DsxDsy account for the depth-averaged dispersive flux due to the non-uniformity of both the velocity and solid concentration profiles throughout the flow column, defined as


Furthermore, assuming that the volumetric solid concentration in the bed layer is 1ξ with ξ denoting the bulk porosity of the bed layer, the solid mass conservation in the bed layer requires that


2.3 Constitutive model for complex viscoplastic flows

So far, there is not a universal closure relation for representing the viscous terms in complex non-Newtonian flows. Constitutive formulations used for environmental sediment-water mixtures are mainly derived from 3D general rheological models which, assuming isotropic material and isochoric flow, allow to express the deviatoric component of the stress tensor σ in the material as


where DDij=12jui+iujij=xyz is the rate of deformation tensor and Φ1 is a scalar function of the second invariant I2D=12trD2 of the rate of deformation tensor D [18] and depends on multiple factor, such as cohesive stress, pore-fluid pressure or flow initial regime. The generalized viscoplastic model, also called Herschel-Bulkley model, assumes dependence of Φ1 on the three parameters: τ0 the cohesive-frictional yield strength, K a consistency coefficient, and m a parameter characterizing the rheological response of the mixture [13, 18], known as behavior index. It is worth mentioning that the dimensions of consistency coefficient K depends on the behavior index m. Therefore, the function Φ1 is expressed as


Considering simple shear state along the flow direction, the velocity vector throughout the flow column is expressed as u=UznuxUznuy0, where Uz is the modulus of the bulk mixture velocity and nu=nuxnuy is the velocity unit vector. Therefore


Replacing (26) and (25) into (24) leads to


being τz the shear stress along the flow direction, which depends on the fluid rheology function (25). For the Newtonian constitutive model, the yield strength τ0 is null, KPas is the dynamic viscosity of the fluid μ and m=1. The generalized model reduces to τz=μdUdz for viscous Newtonian fluids. A widespread non-Newtonian constitutive relation for geophysical surface flows in laminar regime is the Bingham model, which considers a pure cohesive yield stress τ0=τy for the flow initiation, K=μBPas the Bingham plastic viscosity and m=1. In this case, the generalized model reduces to τz=τy+μBdUdz.

The frictional non-linear viscoplastic model considers a Coulomb-Terzaghi linear relation between the effective normal stress σez and the shear stress, hence


where Pz denotes de pore-fluid pressure and δf accounts for the effective friction angle bewteen solid particles [19]. Estimation of the pore pressure distribution Pz is a challenging task [1, 19, 20], although its effects on the reduction of the intergranular shear stress seem to be demonstrated [5, 21, 22]. The simplest models divide the pore pressure into a hydrostatic component plus a dynamic additive component Pz=1+Eρwgzsz, being ρw the density of the pore-liquid and E a tunning coefficient which usually takes values from about 0.4 to 0.8.

Using (28), and considering a plastic viscosity K=μPPasm, the generalized model reduces to


Formulation of closure models for the bottom shear stress requires to integrate the deviatoric stress tensor (24) throughout the flow column. This is not a trivial problem since the structure of the flow along the vertical direction is lost and only the averaged quantities are available. Assuming the deviatoric stress tensor (27), all the depth-averaged stress terms TxxTxy and TyxTyy(16) are null, and


where τb=τzb is the shear stress at the basal interface along the flow direction. In order to obtain a depth-averaged formulation for the basal shear stress τb, the distributed shear stress function τz must be integrated along the flow column. This allows to obtain both the velocity distribution along the vertical direction and the basal resistance τb expressed as a function of the depth-averaged flow variables. If the pore-pressure excess in (28) is considered linear, the the constitutive eq. (29) can be rewritten as


being τf=ρghPbtanδf the value of the frictional yield stress at the basal surface and Pb is the pore-pressure at the bed surface (Figure 2).

Figure 2.

Velocity and stress distribution for the generalized non-linear frictional model.

Therefore, assuming the induced shear distribution along the flow column follows


the velocity derivative along the vertical direction can expressed as


and integrating (34) leads to the velocity vertical distribution


for the non-linear viscoplastic model. Note that the velocity at the free surface UzsUh and the depth averaged velocity U¯ can be expressed as


Integrating (35) throughout the flow column allows to obtain the basal shear stress τb as a function of the flow depth h, the averaged flow velocity U¯=u¯2+v¯2 and the basal frictional yield strength τf using


It is worth mentioning that (37) represents a generalized depth-integrated formulation for viscoplastic flows (Figure 3) which encompasses: pseudoplastic or shear-thinning behavior for m<1, reducing the apparent viscosity as the induced shear rate increases; linear viscoplastic behavior for m=1, with a linear relation between shear stress and shear rate; and dilatant or shear-thickening behavior for m>1, increasing the apparent viscosity as the induced shear rate grows.

Figure 3.

Basal resistance behavior for the generalized frictional non-linear viscoplastic model.

2.4 Net mass exchange between bed and flow layers

The net volumetric exchange Nb between the underlying bed layer and the mixture flow column at the bed surface zbtxy, appearing in (6), (20) and (23), is usually modeled as the balance between the entrainment and the deposition vertical solid fluxes, Eb and Db respectively, leading to


involving the bed surface solid concentration ϕb=1ξ.

The deposition rate Db is commonly related to the solid particles settling velocity in the mixture ωsm and to the near-bed solid concentration in the flow ϕzzb, which is usually related to the depth-averaged solid concentration in the flow as ϕzzb=αϕ¯, being α an adaptation or recovery coefficient, leading to Db=αωsmϕ¯.

The settling velocity of the solid particles ωsm in highly concentrated mixtures is influenced by the presence of other solid particles. Furthermore, in dense-packed mixtures with moderate plastic fine fractions in the flow column such as muddy slurries, the particle settling velocity can be strongly reduced by the development of internal yield stresses in the pore-fluid. Richardson and Zaki [23] proposed ωsm=1ϕ¯mωs, being ωs the settling velocity of the sediment particles in clear water and m a hindering empirical exponent depending on the Reynolds particle number (Rep=ωsds/ν, with ν the clear water kinematic viscosity) which usually takes values close to m=4. Therefore


where αD is a dimensionless parameter which accounts for mulple factor and requires calibration.

The erosion solid flux Eb is directly related to the turbulent fluctuation of the volumetric solid concentration and flow velocity near the bed surface. We assume that this near-bed erosion rate is at the capacity of the flow to entrain solid material from the underlying bed layer, hence it is related to the settling velocity of the particles in clear water ωs and the near equilibrium concentration ϕzzb. The near-bed equilibrium concentration is related to the depth-averaged equilibrium concentration ϕ¯ as ϕzzb=αϕ¯, being α an adaptation coefficient under equilibrium conditions. When equilibrium solid transport states are reached, the adaptation coefficients α and α coincide but in non-equilibrium states αα generally. Therefore, for the sake of simplicity, we assume that the vertical erosion rate Eb is expressed as


where αE is a dimensionless empirical parameter which requires calibration. The capacity solid concentration is usually computed as ϕ¯=qs/hu¯2+v¯2, where qs accounts for the value of the solid transport throughout the flow column in capacity or equilibrium condition, which can be estimated using the multiple empirical relationships from the local hydrodynamic variables [8].


3. Efficient numerical tools for multi-grain mud/debris flows

Considering a multi-grain mixture flow, the resulting system is composed by 3+N+1 conservation equations accounting for the bulk mass (6) and momentum (13), the transport of the N sediment classes (20) and the bed elevation evolution (23). The dimensionless bulk density r can be expressed by defining a new variable ϕχ, referred to as buoyant solid concentration


where ρs,p and ϕp are the density and depth-averaged volumetric concentration of the pth solid phase respectively. Using (41), the equations forming the system can be recast as five conservation laws and rewritten in vector form as


where U is the vector of conservative variables and EU=FUGU are the convective fluxes along the X=xy horizontal coordinates, respectively, expressed as


The vector SbU accounts for the momentum source term associated to the variation of the pressure force on the bed interface, whereas SτU is the momentum dissipation due to the basal resistance. Finally, the source term EbU accounts for the bulk mass exchange between the mixture flow and the bed layer.


3.1 Finite volume method in unstructured meshes

System (42) is time dependent, non linear and contains mass and momentum source terms. Under the hypothesis of dominant advection it can be classified as belonging to the family of hyperbolic systems. In order to obtain a numerical solution, the spatial domain is divided in computational cells using a fixed-in-time mesh and system (42) is integrated in each cell Ωi. Applying the Gauss theorem leads to


being NE the number of edges for the i cell, n=nxnyk the outward unit normal vector, lk the length of the edge and hence Enk the value of the normal flux through the kth edge (Figure 4).

Figure 4.

Computational cells in triangular meshes and local coordinates at the kth cell edge.

The left hand side of (42), also called homogeneous part, satisfies the rotation invariant property [24] and hence the conservative normal flux vector can be rewritten as


being R a rotation matrix which projects the global orthogonal framework X=xy into the local framework X̂=x̂ŷ of the kth cell edge (Figure 4), corresponding to the normal and the tangential directions to the edge respectively. The rotation matrix Rk and its inverse Rk1 are defined as


and the set of local conservative variables Û and fluxes FÛ is defined as


where û=unx+vny and v̂=uny+vnx are the components of the flow velocity û in the local framework.

The value of the fluxes through the kth cell edge can be augmented incorporating the contribution of the momentum source terms Sb and Sτ into the homogeneous normal fluxes FÛ [25]. The momentum terms can be included within the local framework x̂ŷ using the spatial discretization


where HÛk and TÛk are the integrated bed pressure [26] and basal resistance [27] throughout the kth cell edge, expressed in the local framework, being


Using (46) and (49), the augmented flux at the kth edge can be defined as


The net exchange flux term EbU is discretized in space as


and, assuming a piecewise uniform representation of the variables at the i cell for the time t=tn and using explicit temporal integration for the mass and momentum source terms, the updating formulation for the conservative variables U is expressed as


being Δt=tn+1tn the time step, Ai the discrete cell area and FÛinÛjnk the numerical flux through the kth cell edge computed as a function of the local conservative variables at the neighboring cells i and j. The numerical fluxes are here upwind computed using a fully-coupled Roe-type Riemann solver (RS) adapted to compressible two-phase shallow flows. Note that the flow density and depth remain coupled in both the conservative variables and fluxes on the left hand side of the equations, improving the robustness and accuracy of the solution when large density gradients appear in the flow [28]. The RS formulation is based on the augmented Roe strategy [25], ensuring the well-balance in steady states and the correct treatment of wet-dry fronts without requiring additional time step restrictions. A detailled explanation of this fully-coupled RS can be found in Aranda et al. [29].

3.2 High-performance-computing in GPU’s

Graphics Processing Units (GPU’s), were developed to control and manage the graphic operations for video games but currently they have become a powerful tool for solving engineering problems. The main advantage resides in the high number of computation nodes available for workload distribution, leading to higher speed-ups without an increment of investment in large facilities. The main drawback is that GPU-kernels usually require an intensive programming effort compared with multi-CPU strategies, such as OpenMP or other shared-memory strategies. There exist several platforms for High-Performance-Computing (HPC). NVIDIA developed the CUDA (Compute Unified Device Architecture), a computing platform and programming environment which incorporates directives for massive parallel operations. To ensure the applicability of the model to realistic large-scale mud/debris flows, the updating eq. (54) is solved using a GPU-based algorithm implemented in the NVIDIA CUDA/C++ framework.

The preprocess step and CPU-GPU memory transfer are implemented to run on one CPU core, whereas the time loop computation is accelerated using GPU. However, some tasks inside the time loop are controlled yet by the CPU, such as the time advance control, the boundary conditions application and the output data dump. Therefore, it is necessary to transfer information from/to the GPU at each time-step. While the computational effort required for the time and boundaries transference is considerably smaller than that of each kernel function, in order to dump the intermediate output information, all the variables in the domain must be transferred from GPU device to CPU host.

The CUDA toolkit allows that all the processed elements can be distributed by threads and blocks of threads. Each thread uses its own thread index to identify the element to be processed, launching several execution threads at the same time (parallel computation). As computing GPU devices are well designed to work efficiently with ordered information, the variables needed for computation are stored in the GPU memory as structures of arrays (SoA), improving the spatial locality for memory accesses. Only the kernel functions, which require a higher computational effort, have been implemented to run on the GPU device. Some tasks in the GPU kernel are optimized using the CUBLAS library included in CUDA. The memory transfer between the CPU host and the GPU device has been reduced as much as possible for each time step.


4. Application to the mine tailings dam failure in Brumadinho (Brasil)

In this section, the proposed methodology is applied to the catastrophic large-scale mud flow occurred Brumadinho (Minas Gerais, Brazil, 2019). The sudden failure of a mine containing dike resul in an extremely violent tailings flow which traveled downstream more than 10 km and reached the Paraopeba River. This disaster caused more than 260 deaths and important economic and environmental losses. The dam contained almost 12106m3 mining waste tailings with a height of 7080m and covering an area of 4.13105m2. Figure 5 shows an aerial image of the mine site after the dam collapse. The thalweg elevation along the area covered by the mud varies between 860m.o.s.l at the dam-toe and 720m.o.s.l at the Paraopeba River, with an averaged longitudinal bed slope S0=0.0165m/m. The area affected by the mud was 3.3106m2.

Figure 5.

Aerial image of the area affected by the mud.

The dike material and the tailings rapidly became a heavy liquid that flowed downstream at a high speed. Tailings were composed by a mixture of water, sediments and heavy metals, mainly iron, aluminum, manganese and titanium [30]. The size distribution consisted basically of a mineral sand fraction (38%) and a fines fraction (62%), accounting for mineral silt-clay and metals particles. The water content before the failure was estimated around 50% by volume with a specific weight of 2226kN/m3 (Table 1).

Dam capacity12106m3
Dam area4.13105m2
Solid concent. ϕ050%
Specific weight2226kN/m3
Size distr.SandFines
Rel. content38%62%
Heavy metalsFeAlMnTi
Weight conc.264.9 mg/g10.8 mg/g4.8 mg/g0.5 mg/g

Table 1.

Brumadinho’s dam and tailings features.

The spatial domain is discretized using a unstructured triangular mesh with 6105 cells approximately and four control cross-sections are placed downstream of the dam (Figure 5), corresponding to (CS-1) the mine stockpile area, (CS-2) the railway bridge, (CS-3) the national road and (CS-4) the gauge station in the river. Furthermore, the base regime water depth and velocities in the Paraopeba River before to the mud arrival is assessed by a previous simulation. The initial tailings depth at the dam is estimated by comparing the terrain elevation before and after the dam failure using 1 × 1 m DTM’s. Six different solid phases are set for characterizing the fully saturated tailings material, including mineral sand, mineral silt, iron (Fe), aluminum (Al), manganese (Mn) and titanium (Ti), with initial bulk volumetric concentration ϕ0=0.5 and normalized density r2.25. The deposition porosity ξ for the bed layer is estimated using the Wu relation [8] and the hiding-exposure effects on the critical Shields stress for the incipient motion of each solid phase θc,p are estimated using the Egiazaroff formula [31]. Also, the solid transport capacity of the flow qs is computed using the Wu relation [8]. The simulated time was 3 h from the dam collapse.

The mining tailings in the dam showed a low plasticity and high values of pore-fluid pressure [30] before the collapse. Assuming the turbulent behavior of the flow, the behavior index is set m=2, setting a quadratic relation in (37) and allowing to calculate the apparent consistency parameter μP as a function of the Manning’s roughness coefficient of the terrain (nb=0.065sm1/3). The basal stability angle is estimated in δf=5°. The basal pore pressure is estimated as Pb=1+Ebρwgh, being Eb a basal pressure factor which tends to Eb for totally fluidized materials. Therefore, the total frictional-turbulent basal shear stress is estimated as


where the first term on the right hand side accounts for the frictional yield stress whereas the second term denotes the turbulent shear component.

In order to assess the influence of the basal pressure in the numerical results, the entrainment term EbU in (42) is neglected and the value of Eb is varied from Eb (pure Newtonian turbulent behavior) to Eb=0.9 (medium frictional yield strength). Figure 6 shows the temporal evolution of the wave-front location and the total flow rate at the control section CS-2. Note that, for each value of Eb, the wave front loses velocity progressively as it moves downstream. As Eb decreases, the frictional shear stress is enhanced and the mobility of the flow decreases considerably. That it is clearly shown in the marked reduction of the flow rate at the control section CS-2 as the pore pressure factor is increased. Setting Eb=1.05 allows to predict the observed arrival time (3247min) of the tailings wave to the Paraopeba river at 8.5km downstream the dam (gray rectangle in Figure 6).

Figure 6.

Temporal evolution of (left) the wave-front location and (right) the total flow rate at the control section CS-2.

Figure 7 shows the mud depth at t=0min, t=5min, t=15min and t=45min after the dam collapse for Eb=1.05. The numerical results show that practically the whole initial tailing volume flows out of the dam in the first 5 min after the dam collapse, as it was observed in the available videos. The mud wave moves downstream with a computed height larger than 20m in some zones during the first minutes. After this initial stage, the numerical results show that the dambreak wave is still higher than 10m when the wave-front reaches the railway bridge (CS-2). During the real event, the mud wave impact caused the collapse of the railway bridge structure. At t=45min, the mud wave-front has reached the Paraopeba River and the flow is practically stopped.

Figure 7.

Mud depth at t=0min, t=5min, t=15min and t=45min after the dam collapse.

Furthermore, in order to assess the influence of the solid material entrainment from the erodible bed to the flow layer, the pore-pressure factor is set to Eb=1.05 but the entrainment term EbU in (42) is now considered. The relation αE/αD, which controls the balance between the erosion component Eb and the deposition component Db of the vertical bed-flow exchange flux, is varied to αE/αD=1%, αE/αD=5%, αE/αD=10% and αE/αD=50%. Figure 8 depicts the predicted temporal evolution of both the mud wave-front location and the total volume of the fluidized material in the domain as the entrainment parameter is varied. The higher the relation αE/αD, the larger the fluidized mass involved in the flow as the dambreak wave progresses. Hence, the entrainment of bed material leads to a greater flow mobility, causing the wave-front to reach the Paraopeba River faster as the entrainment parameter αE/αD increases. Even, for αE/αD=50%, the mobilized mass becomes more than 150% of the initial mass and the wave-front reaches the river faster than the fixed-bed pure-turbulent case (dashed gray line).

Figure 8.

Temporal evolution of (left) the wave-front location and (right) the total fluidized volume as the relation αE/αD is increased.

Furthermore, the entrainment of bed materials also leads to the segregation of the solid phases along the flow. As αE/αD is increased, horizontal gradients in the concentration of the different solid classes composing the mixture tend to appear (Figure 9). These gradients depends on the temporal evolution of the class-specific erosion and deposition vertical fluxes, Ep and Dp respectively, along the tailings flow.

Figure 9.

Bulk solid concentration in the flow at t=25min after the dam collapse for (a) αE/αD=1%, (b) αE/αD=5%, (c) αE/αD=10% and (d) αE/αD=50%.

The increment of the flow mobility caused by the bed material entrainment is clearly shown in the temporal evolution of the flow rate at the control sections CS-2 and CS-3 (Figure 10). At CS-2 the increment of the entrainment relation αE/αD is mainly appreciated in the increment of the flow rate peak from 7500m3/s for the fixed-bed simulation to almost 15,000m3/s for the movable bed case with αE/αD=50%. Furthermore, at CS-3, the entrainment of bed material does not only increase the flow rate peak respect to the fixed-bed assumption, but also leads to a noticeable decrease in the arrival time of the mud wave from 40 min to 25 min approximately.

Figure 10.

Temporal evolution of the flow rate at (left) the control sections CS-2 and (right) the control section CS-3 as αE/αD is increased.

Finally, the computational efficiency of the GPU-accelerated kernel is shown in Table 2 for different GPU devices. The performance obtained with the GPU-based code is compared with the single-core CPU simulation time. Note that the speed-up obtained respect to the single-core performance has increased progressively as the devices have been developed in the last 5 years. Currently, the speed-up obtained with the GPU accelerated kernel is about 280 (Nvidia A100), more than 8 times the Nvidia Tesla K40c performance. That means that for achieving the GPU-parallelized code performance with a CPU-based algorithm, a cluster with at least 280 CPU cores is required (see speed-up in Table 2).

Computational timeSpeed-up
CPU Intel i7-10700F59.43 h
Nvidia Tesla K40c1.79 h× 33
Nvidia GTX 1080 Ti58.2 min× 61
Nvidia Tesla V10026.5 min× 135
Nvidia A100 40GB12.8 min× 279

Table 2.

Computational times with GPU-based and CPU-based algorithms.


5. Conclusions

The numerical modeling of geophysical surface flows involving sediment transport must be addressed using a comprehensive strategy: beginning with the derivation of proper shallow-type mathematical models, following by the development of robust and accurate numerical schemes within the Finite Volume (FV) framework and ending with the implementation of efficient HPC algorithms. This integrated approach is required to release Efficient Simulation Tools (EST) for environmental processes involving sediment transport with realistic temporal and spatial scales.

The derivation of the generalized 2D system of depth-averaged conservation laws for environmental surface flows of water-sediment mixtures over movable bed conditions is a fundamental step to understand the physical consequences of the mathematical shallow-type simplification. From the mathematical modeling approach, the fluidized material in motion is contained in a flow layer consisting of a mixture of water and multiple solid phases. This flow layer usually moves rapidly downstream steep channels and involves complex topography. The liquid-solid material can be exchanged throughout the bottom interface with the underlying static bed layer, hence involving also a transient bottom boundary for the flow layer. These features lead to an increasing complexity for the mathematical simplification of sediment transport surface flows.

The 2D shallow-type system of equations for variable-density multi-grain water-sediment flows can be solved using a Finite Volume (FV) methods, supplemented preferable with upwind augmented Riemann solvers which provide a robust and accurate computation of the intercell fluxes even involving highly transient density interfaces and ensure the well-balanced character of the solution in quiescent and steady states. The usage of GPU-accelerated kernels for sequential computation of the FV solution is an efficient and low cost alternative to the traditional multi-CPU strategies. The GPU solution must be designed taking into account the fact that the GPU is an independent device with its own RAM memory. This means that the memory transfer between the conventional RAM memory and the GPU device memory plays a key role in the performance of GPU-accelerated software.

GPU-accelerated codes are mandatory for large-scale sediment-laden flow models since they allows the calibration of the multiple processes involved in the movement of the fluidized water-sediment material without spending weeks or even months waiting results. Probably, the most challenging and unknown process involved in sediment-laden surface flows is the development of pore-fluid pressure, hence its estimation requires careful calibration procedure. This pore pressure affects the frictional shear stress between solid particles, reducing the basal resistance and affecting the flow mobility. Furthermore, special attention has to be paid to the calibration of the entrainment/deposition parameters. The entrainment of material from the underlying movable bed to the mixture layer increases the mass and momentum of the fluidized material in movement, enlarging the mobility of the flow and reducing the arrival time of the mud waves.



This work is part of the research project PGC2018-094341-B-I00 Fast and accurate tools for the simulation and control of environmental flows (FACTOOL) funded by the Ministry of Science and Innovation and FEDER.


  1. 1. Iverson RM. The physics of debris flows. Reviews of Geophysics. 1997;35(3):245-296
  2. 2. Pierson TC. Hyperconcentrated flow—Transitional process between water flow and debris flow. In: Debris-flow Hazards and Related Phenomena. Berlin, Germany: Springer Berlin Heidelberg; 2005
  3. 3. Hungr O, Evans SG, Bovis MJ, Hutchinson JN. A review of the classification of landslides of the flow type. Environmental and Engineering Geoscience. 2001;7(3):221-238
  4. 4. Rickenmann D, Weber D, Stepanov B. Erosion by debris flows in field and laboratory experiments. In: 3rd International Conference on Debris-Flow Hazards: Mitigation, Mechanics, Prediction and Assessment. Rotterdam: Millpress; 2003. pp. 883-894
  5. 5. Iverson RM, Reid ME, Logan M, LaHusen RG, Godt JW, Griswold JP. Positive feedback and momentum growth during debris-flow entrainment of wet bed sediment. Nature Geoscience. 2011;4:116-121
  6. 6. Berger C, McArdell BW, Schlunegger F. Direct measurement of channel erosion by debris flows, illgraben, Switzerland. Journal of Geophysical Research: Earth Surface. 2011;116(F1):93-104
  7. 7. McCoy SW, Kean JW, Coe JA, Tucker GE, Staley DM, Wasklewicz TA. Sediment entrainment by debris flows: In situ measurements from the headwaters of a steep catchment. Journal of Geophysical Research: Earth Surface. 2012;117:F03016
  8. 8. Wu W. Computational River Dynamics. NetLibrary, Inc. London: CRC Press; 2007
  9. 9. Murillo J, García-Navarro P. Wave riemann description of friction terms in unsteady shallow flows: Application to water and mud/debris floods. Journal of Computational Physics. 2012;231:1963-2001
  10. 10. Juez C, Murillo J, García-Navarro P. 2d simulation of granular flow over irregular steep slopes using global and local coordinates. Journal of Computational Physics. 2013;255:166-204
  11. 11. Gualtieri C, Ianniruberto M, Filizola N. On the mixing of rivers with a difference in density: The case of the negro/solimões confluence, Brazil. Journal of Hydrology. 2019;578:124029
  12. 12. Cunge JA, Holly FM, Verwey A. Practical aspects of computational river hydraulics. In: Monographs and Surveys in Water Resources Engineering. London: Pitman Advanced Publishing Program; 1980
  13. 13. Dewals B, Rulot F, Erpicum S, Archambeau P, Pirotton M. Sediment Transport, Chapter Advanced Topics in Sediment Transport Modelling: Non-alluvial Beds and Hyperconcentrated Flows. Rijeka, Croatia: IntechOpen; 2011
  14. 14. Lacasta A, Juez C, Murillo J, García-Navarro P. An efficient solution for hazardous geophysical flows simulation using GPUs. Computers & Geosciences. 2015;78:63-72
  15. 15. Ming X, Liang Q, Xia X, Li D, Fowler HJ. Real-time flood forecasting based on a high-performance 2d hydrodynamic model and numerical weather predictions. Water Resources Research. 2020;56:e2019WR025583
  16. 16. Juez C, Lacasta A, Murillo J, García-Navarro P. An efficient GPU implementation for a faster simulation of unsteady bed-load transport. Journal of Hydraulic Research. 2016;54(3):275-288
  17. 17. de la Asuncion M, Castro MJ. Simulation of tsunamis generated by landslides using adaptive mesh refinement on gpu. Journal of Computational Physics. 2017;345:91-110
  18. 18. Pastor M, Blanc T, Haddad B, Drempetic V, Sanchez-Morles M, Dutto P, et al. Depth averaged models for fast landslide propagation: Mathematical, rheological and numerical aspects. Archives of Computational Methods in Engineering. 2015;22:67-104
  19. 19. Jakob M, Hungr O. Debris-Flow Hazards and Related Phenomena. Springer Praxis Books. Springer Berlin Heidelberg; 2005
  20. 20. Iverson RM, Vallance JW. New views of granular mass flows. Geology. 2001;29(2):115-118
  21. 21. Berti M, Simoni A. Experimental evidences and numerical modeling of debris flow initiated by channel runoff. Landslides. 2005
  22. 22. McArdell BW, Bartelt P, Kowalski J. Field observations of basal forces and fluid pore pressure in a debris flow. Geophysical Research Letters. 2007;34(7):L07406
  23. 23. Richardson JF, Zaki WN. Sedimentation and fluidisation: Part I. Chemical Engineering Research and Design. 1997;75:82-100
  24. 24. Godlewski E, Raviart P-A. Numerical Approximation of Hyperbolic Systems of Conservation Laws. New York: Springer-Verlag; 1996
  25. 25. Murillo J, Navas-Montilla A. A comprehensive explanation and exercise of the source terms in hyperbolic systems using Roe type solutions. Application to the 1D-2D shallow water equations. Advances in Water Resources. 2016;98:70-96
  26. 26. Murillo J, García-Navarro P. Energy balance numerical schemes for shallow water equations with discontinuous topography. Journal of Computational Physics. 2012;236:119-142
  27. 27. Martínez-Aranda S, Murillo J, Morales-Hernández M, García-Navarro P. Novel discretization strategies for the 2D non-Newtonian resistance term in geophysical shallow flows. Engineering Geology. 2022;302:106625
  28. 28. Martínez-Aranda S, Murillo J, García-Navarro P. A robust two-dimensional model for highly sediment-laden unsteady flows of variable density over movable beds. Journal of Hydroinformatics. 2020;22(5):1138-1160
  29. 29. Martínez-Aranda S, Murillo J, García-Navarro P. A GPU-accelerated efficient simulation tool (EST) for 2D variable-density mud/debris flows over non-uniform erodible beds. Engineering Geology. 2021;296:106462
  30. 30. Robertson PK, de Melo L, Williams DJ, Ward Wilson G. Report of the Expert Panel on the Technical Causes of the Failure of Feijão Dam I. Technical Report. Brasil: Vale S.A.; 2019
  31. 31. Egiazaroff IV. Calculation of nonuniform sediment concentrations. Proceedings of the American Society of Civil Engineers. 1965;91:225-247

Written By

Sergio Martínez-Aranda and Pilar García-Navarro

Reviewed: 24 August 2022 Published: 11 October 2022