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

Abstract

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.

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 $\rho ={\rho }_{w}n+{\rho }_s\varphi$ and linear momentum $\rho \mathbf{u}={\rho }_{w}n{\mathbf{u}}_{\mathbf{w}}+{\rho }_s\varphi {\mathbf{u}}_{\mathbf{s}}$ of the mixture, being ${\rho }_w$ and ${\rho }_s$ the density of the fluid and solid phases, respectively, $\varphi$ the volumetric solid-phase concentration and $n=1-\varphi$ the volumetric fluid-phase fraction or mixture porosity, ${\mathbf{u}}_{\mathbf{w}}$ the velocity of the pore-fluid and ${\mathbf{u}}_{\mathbf{s}}$ 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 $\mathbf{X}=\left(x,y,z\right)$ as

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)=0$$E1

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}\right)=\mathbf{F}-\nabla p+\nabla \cdot \tau$$E2

where $\mathbf{u}=\left(u_x,u_y,u_z\right)$ is the bulk velocity in any point of the fluidized material, $\mathbf{F}=\left(F_x,F_y,F_z\right)$ are the external forces, such as gravity and $p$ denotes the pressure of the mixture. The term $\tau ={\tau }_{\mathit{ij}}\left(i,j=x,y,z\right)$ 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 $z_s=z_s\left(t,x,y\right)$ and the bottom surface of the flow column $z_b=z_b\left(t,x,y\right)$, which is also considered a movable interface. The kinematic conditions at these boundaries can be expressed as

$$\frac{\partial z_s}{\partial t}+{\left(u_x\right)}_s\frac{\partial z_s}{\partial x}+{\left(u_y\right)}_s\frac{\partial z_s}{\partial y}={\left(u_z\right)}_s$$E3

$$\frac{\partial z_b}{\partial t}+{\left(u_x\right)}_b\frac{\partial z_b}{\partial x}+{\left(u_y\right)}_b\frac{\partial z_b}{\partial y}={\left(u_z\right)}_b+N_b$$E4

being $N_b$ the net volumetric flux through the bed interface along the $z-$coordinate, being positive for net entrainment conditions and negative for net deposition fluxes. The subscripts ${\left(\cdot \right)}_s$ and ${\left(\cdot \right)}_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

$$\underset{z_b}{\overset{z_s}{\int }}{\partial }_t\left(\rho \right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}{\partial }_x\left(u_x\rho \right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}{\partial }_y\left(u_y\rho \right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}{\partial }_z\left(u_z\rho \right)\mathrm{d}z=0$$E5

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

$$\frac{\partial \left(\overline{\rho }h\right)}{\partial t}+\frac{\partial }{\partial x}\left(\overline{\rho }h\overline{u}\right)+\frac{\partial }{\partial y}\left(\overline{\rho }h\overline{v}\right)=-{\left(\rho \right)}_b{N}_b$$E6

where the depth-averaged of the mixture bulk density $\overline{\rho }$ and velocities are defined as

$$\overline{\rho }h=\underset{z_b}{\overset{z_s}{\int }}\rho \mathrm{d}z\overline{u}=\frac{1}{\overline{\rho }h}\underset{z_b}{\overset{z_s}{\int }}{u}_x\rho \mathrm{d}z\overline{v}=\frac{1}{\overline{\rho }h}\underset{z_b}{\overset{z_s}{\int }}{u}_y\rho \mathrm{d}z$$E7

being $h=z_s-z_b$ the flow depth and ${\left(\rho \right)}_b$ the mixture bulk density at the bed interface $z_b$.

Regarding the momentum depth integration, the volumetric force components in (1) and (2) are considered for the sake of simplicity as $F_x=F_y=0$ and $F_z=-\mathit{\rho 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

$$\begin{array}{l}\underset{z_b}{\overset{z_s}{\int }}{\partial }_t\left(\rho u_x\right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}{\partial }_x\left(u_x\rho u_x\right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}{\partial }_y\left(u_y\rho u_x\right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}{\partial }_z\left(u_z\rho u_x\right)\mathrm{d}z\\ =-\underset{z_b}{\overset{z_s}{\int }}{\partial }_x\left(p\right)\mathrm{d}z+\underset{z_b}{\overset{z_s}{\int }}\left({\partial }_x{\tau }_{\mathit{xx}}+{\partial }_y{\tau }_{\mathit{xy}}+{\partial }_z{\tau }_{\mathit{xz}}\right)\mathrm{d}z\end{array}$$E8

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

$$\frac{\partial \left(\overline{\rho }h\overline{u}\right)}{\partial t}+\frac{\partial }{\partial x}\left(\overline{\rho }h{\overline{u}}^2-\overline{\rho }{\mathit{hD}}_{\mathit{xx}}\right)+\frac{\partial }{\partial y}\left(\overline{\rho }h\overline{u}\overline{v}-\overline{\rho }{\mathit{hD}}_{\mathit{xy}}\right)+{\left(\rho u_x\right)}_b{N}_b$$E9

with $\left(D_{\mathit{xx}},D_{\mathit{xy}}\right)$ accounting for the depth-averaged dispersion momentum transport due to the non-uniformity of the vertical velocity profile, and defined as

$$D_{\mathit{xx}}=-\frac{1}{\overline{\rho }h}\underset{z_b}{\overset{z_s}{\int }}\rho {\left(u_x-\overline{u}\right)}^2\mathrm{d}z$$E10

$$D_{\mathit{xy}}=-\frac{1}{\overline{\rho }h}\underset{z_b}{\overset{z_s}{\int }}\rho \left(u_x-\overline{u}\right)\left(u_y-\overline{v}\right)\mathrm{d}z$$E11

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 $x-$coordinate can be expressed as

$$-\underset{z_b}{\overset{z_s}{\int }}{\partial }_x\left(p\right)\mathrm{d}z=-g\underset{z_b}{\overset{z_s}{\int }}\left(\frac{\partial }{\partial x}\underset{z}{\overset{z_s}{\int }}\rho \mathrm{d}z\right)\mathrm{d}z=\underset{\text{Conservative\hspace{0.17em}term}}{\underbrace{-\frac{\partial }{\partial x}\left(\frac{1}{2}g\overline{\rho }{h}^2\right)}}\underset{\mathrm{Bed}-\text{pressure\hspace{0.17em}term}}{\underbrace{-g\overline{\rho }h\frac{\partial z_b}{\partial x}}}$$E12

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 $x-$coordinate

$$\begin{array}{l}\frac{\partial \left(\overline{\rho }h\overline{u}\right)}{\partial t}+\frac{\partial }{\partial x}\left(\overline{\rho }h{\overline{u}}^2+\frac{1}{2}g\overline{\rho }{h}^2\right)+\frac{\partial }{\partial y}\left(\overline{\rho }h\overline{u}\overline{v}\right)=-g\overline{\rho }h\frac{\partial z_b}{\partial x}+{\tau }_{\mathit{sx}}-{\tau }_{\mathit{bx}}\\ +\frac{\partial }{\partial x}\left(\overline{\rho }h\left(T_{\mathit{xx}}+D_{\mathit{xx}}\right)\right)+\frac{\partial }{\partial y}\left(\overline{\rho }h\left(T_{\mathit{xy}}+D_{\mathit{xy}}\right)\right)\\ -{\left(\rho u_x\right)}_b{N}_b\end{array}$$E13

and, similarly, for the $y-$coordinate. The boundary stress terms in (13) at the free surface $z_s$ and the bed interface $z_b$ are

$${\tau }_{\mathit{sx}}={\left({\tau }_{\mathit{xz}}\right)}_s-{\left({\tau }_{\mathit{xx}}\right)}_s\frac{\partial z_s}{\partial x}-{\left({\tau }_{\mathit{xy}}\right)}_s\frac{\partial z_s}{\partial y}$$E14

$$-{\tau }_{\mathit{bx}}=-{\left({\tau }_{\mathit{xz}}\right)}_b+{\left({\tau }_{\mathit{xx}}\right)}_b\frac{\partial z_b}{\partial x}+{\left({\tau }_{\mathit{xy}}\right)}_b\frac{\partial z_b}{\partial x}$$E15

where ${\tau }_{\mathit{sx}}$ denotes the $x-$coordinate component of the wind action at the free surface, whereas ${\tau }_{\mathit{bx}}$ represents the $x-$coordinate 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 $\left(T_{\mathit{xx}},T_{\mathit{xy}}\right)$ are defined as

$$T_{\mathit{xx}}=\frac{1}{\overline{\rho }h}\underset{z_b}{\overset{z_s}{\int }}{\tau }_{\mathit{xx}}\mathrm{d}z{T}_{\mathit{xy}}=\frac{1}{\overline{\rho }h}\underset{z_b}{\overset{z_s}{\int }}{\tau }_{\mathit{xy}}\mathrm{d}z$$E16

2.2 Depth-integrated solid transport equation

The solid-phase transport process is governed by

$$\begin{array}{l}\frac{\partial \left({\rho }_s\varphi \right)}{\partial t}+\nabla \cdot \left({\rho }_s\varphi {\mathbf{u}}_{\mathbf{s}}\right)=0\end{array}$$E17

where ${\mathbf{u}}_{\mathbf{s}}=\left(u_{\mathit{sx}},u_{\mathit{sy}},u_{\mathit{sz}}\right)$ 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

$$\begin{array}{l}\frac{\partial \varphi }{\partial t}+\nabla \cdot \left(\varphi \mathbf{u}\right)=\nabla \cdot \left[\left(\mathbf{u}-{\mathbf{u}}_{\mathbf{s}}\right)\varphi \right]\end{array}$$E18

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 $\mathbf{u}\approx {\mathbf{u}}_{\mathbf{s}}$ 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

$$\overline{\varphi }h=\underset{z_b}{\overset{z_s}{\int }}\varphi \mathrm{d}z$$E19

where $\overline{\varphi }$ 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

$$\frac{\partial \left(\overline{\varphi }h\right)}{\partial t}+\frac{\partial }{\partial x}\left(\overline{\varphi }h\overline{u}\right)+\frac{\partial }{\partial y}\left(\overline{\varphi }h\overline{v}\right)=-{\left(\varphi \right)}_b{N}_b+\frac{\partial }{\partial x}\left(\overline{\varphi }{\mathit{hD}}_{\mathit{sx}}\right)+\frac{\partial }{\partial y}\left(\overline{\varphi }{\mathit{hD}}_{\mathit{sy}}\right)$$E20

where ${\left(\varphi \right)}_b$ denotes the solid concentration in the bed layer surface $z_b$ and $\left(D_{\mathit{sx}},D_{\mathit{sy}}\right)$ 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

$$D_{\mathit{sx}}=-\frac{1}{\overline{\varphi }h}\underset{z_b}{\overset{z_s}{\int }}\rho \left(u_x-\overline{u}\right)\left(\frac{\varphi }{\rho }-\frac{\overline{\varphi }}{\overline{\rho }}\right)\mathrm{d}z$$E21

$$D_{\mathit{sy}}=-\frac{1}{\overline{\varphi }h}\underset{z_b}{\overset{z_s}{\int }}\rho \left(u_y-\overline{v}\right)\left(\frac{\varphi }{\rho }-\frac{\overline{\varphi }}{\overline{\rho }}\right)\mathrm{d}z$$E22

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

$$\left(1-\xi \right)\frac{\partial z_b}{\partial t}={\left(\varphi \right)}_b{N}_b$$E23

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 $\sigma$ in the material as

$$\tau ={\mathrm{\Phi }}_1\left(I_{2D}\right)\mathbf{D}$$E24

where $\mathbf{D}\equiv D_{\mathit{ij}}=\frac{1}{2}\left({\partial }_j{u}_i+{\partial }_i{u}_j\right)\left(i,j=x,y,z\right)$ is the rate of deformation tensor and ${\mathrm{\Phi }}_1$ is a scalar function of the second invariant $I_{2D}=\frac{1}{2}\mathrm{tr}\left({\mathbf{D}}^2\right)$ of the rate of deformation tensor $\mathbf{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 ${\mathrm{\Phi }}_1$ on the three parameters: ${\tau }_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 ${\mathrm{\Phi }}_1$ is expressed as

$${\mathrm{\Phi }}_1\left(I_{2D}\right)=\frac{{\tau }_0}{\sqrt{I_{2D}}}+2K{\left(4I_{2D}\right)}^{\frac{m-1}{2}}$$E25

Considering simple shear state along the flow direction, the velocity vector throughout the flow column is expressed as $\mathbf{u}=\left(U\left(z\right)n_{\mathit{ux}},U\left(z\right)n_{\mathit{uy}},0\right)$, where $U\left(z\right)$ is the modulus of the bulk mixture velocity and ${\mathbf{n}}_{\mathbf{u}}=\left(n_{\mathit{ux}},n_{\mathit{uy}}\right)$ is the velocity unit vector. Therefore

$$\mathbf{D}=\left(\begin{array}{ccc}0& 0& \frac{1}{2}\frac{\mathrm{d}U}{\mathrm{d}z}{n}_{\mathit{ux}}\\ 0& 0& \frac{1}{2}\frac{\mathrm{d}U}{\mathrm{d}z}{n}_{\mathit{uy}}\\ \frac{1}{2}\frac{\mathrm{d}U}{\mathrm{d}z}{n}_{\mathit{ux}}& \frac{1}{2}\frac{\mathrm{d}U}{\mathrm{d}z}{n}_{\mathit{uy}}& 0\end{array}\right)I_{2D}=\frac{1}{4}{\left(\frac{\mathrm{d}U}{\mathrm{d}z}\right)}^2$$E26

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

$$\tau =\left(\begin{array}{ccc}0& 0& \tau \left(z\right)n_{\mathit{ux}}\\ 0& 0& \tau \left(z\right)n_{\mathit{uy}}\\ \tau \left(z\right)n_{\mathit{ux}}& \tau \left(z\right)n_{\mathit{uy}}& 0\end{array}\right)$$E27

being $\tau \left(z\right)$ the shear stress along the flow direction, which depends on the fluid rheology function (25). For the Newtonian constitutive model, the yield strength ${\tau }_0$ is null, $K\left[\mathit{Pa}\cdot s\right]$ is the dynamic viscosity of the fluid $\mu$ and $m=1$. The generalized model reduces to $\tau \left(z\right)=\mu \frac{\mathrm{d}U}{\mathrm{d}z}$ 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 ${\tau }_0={\tau }_y$ for the flow initiation, $K={\mu }_B\left[\mathit{Pa}\cdot s\right]$ the Bingham plastic viscosity and $m=1$. In this case, the generalized model reduces to $\tau \left(z\right)={\tau }_y+{\mu }_B\frac{\mathrm{d}U}{\mathrm{d}z}$.

The frictional non-linear viscoplastic model considers a Coulomb-Terzaghi linear relation between the effective normal stress ${\sigma }_e\left(z\right)$ and the shear stress, hence

$${\tau }_0={\sigma }_e\left(z\right)tan{\delta }_f=\left[\mathit{\rho g}\left(z_s-z\right)-\mathcal{P}\left(z\right)\right]tan{\delta }_f$$E28

where $\mathcal{P}\left(z\right)$ denotes de pore-fluid pressure and ${\delta }_f$ accounts for the effective friction angle bewteen solid particles [19]. Estimation of the pore pressure distribution $\mathcal{P}\left(z\right)$ 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 $\mathcal{P}\left(z\right)=\left(1+\mathcal{E}\right){\rho }_{w}g\left(z_s-z\right)$, being ${\rho }_w$ the density of the pore-liquid and $\mathcal{E}$ a tunning coefficient which usually takes values from about 0.4 to 0.8.

Using (28), and considering a plastic viscosity $K={\mu }_P\left[\mathit{Pa}\cdot s^m\right]$, the generalized model reduces to

$$\tau \left(z\right)={\sigma }_e\left(z\right)tan{\delta }_f+{\mu }_P{\left(\frac{\mathrm{d}U}{\mathrm{d}z}\right)}^m$$E29

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 $\left(T_{\mathit{xx}},T_{\mathit{xy}}\right)$ and $\left(T_{\mathit{yx}},T_{\mathit{yy}}\right)$(16) are null, and

$${\tau }_{\mathit{bx}}={\left({\tau }_{\mathit{xz}}\right)}_b-{\left({\tau }_{\mathit{xx}}\right)}_b\frac{\partial z_b}{\partial x}-{\left({\tau }_{\mathit{xy}}\right)}_b\frac{\partial z_b}{\partial y}={\tau }_b{n}_{\mathit{ux}}$$E30

$${\tau }_{\mathit{by}}={\left({\tau }_{\mathit{yz}}\right)}_b-{\left({\tau }_{\mathit{yx}}\right)}_b\frac{\partial z_b}{\partial x}-{\left({\tau }_{\mathit{yy}}\right)}_b\frac{\partial z_b}{\partial y}={\tau }_b{n}_{\mathit{uy}}$$E31

where ${\tau }_b=\tau \left(z_b\right)$ 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 ${\tau }_b$, the distributed shear stress function $\tau \left(z\right)$ must be integrated along the flow column. This allows to obtain both the velocity distribution along the vertical direction and the basal resistance ${\tau }_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

$$\tau \left(z\right)={\tau }_f\left(1-\frac{z-z_b}{h}\right)+{\mu }_P{\left(\frac{\mathrm{d}U}{\mathrm{d}z}\right)}^m$$E32

being ${\tau }_f=\left(\mathit{\rho gh}-{\mathcal{P}}_b\right)tan{\delta }_f$ the value of the frictional yield stress at the basal surface and ${\mathcal{P}}_b$ 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

$$\tau \left(z\right)={\tau }_b\left(1-\frac{z-z_b}{h}\right)$$E33

the velocity derivative along the vertical direction can expressed as

$$\frac{\mathrm{d}U}{\mathrm{d}z}={\left[\frac{{\tau }_b-{\tau }_f}{{\mu }_P}\left(1-\frac{z-z_b}{h}\right)\right]}^{1/m}$$E34

and integrating (34) leads to the velocity vertical distribution

$$U\left(z\right)=\frac{m}{m+1}{\left(\frac{{\tau }_b-{\tau }_f}{{\mu }_P}\right)}^{1/m}h\left[1-{\left(1-\frac{z-z_b}{h}\right)}^{\frac{m+1}{m}}\right]$$E35

for the non-linear viscoplastic model. Note that the velocity at the free surface $U\left(z_s\right)\equiv U_h$ and the depth averaged velocity $\overline{U}$ can be expressed as

$$U_h=\frac{m}{m+1}{\left(\frac{{\tau }_b-{\tau }_f}{{\mu }_P}\right)}^{1/m}h\overline{U}=\frac{m+1}{2m+1}{U}_h$$E36

Integrating (35) throughout the flow column allows to obtain the basal shear stress ${\tau }_b$ as a function of the flow depth $h$, the averaged flow velocity $\overline{U}=\sqrt{{\overline{u}}^2+{\overline{v}}^2}$ and the basal frictional yield strength ${\tau }_f$ using

$${\tau }_b={\tau }_f+{\left(\frac{2m+1}{m}\right)}^m{\mu }_P{\left(\frac{\sqrt{{\overline{u}}^2+{\overline{v}}^2}}{h}\right)}^m$$E37

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 $N_b$ between the underlying bed layer and the mixture flow column at the bed surface $z_b\left(t,x,y\right)$, appearing in (6), (20) and (23), is usually modeled as the balance between the entrainment and the deposition vertical solid fluxes, $E_b$ and $D_b$ respectively, leading to

$$N_b=\frac{1}{{\left(\varphi \right)}_b}\left(D_b-E_b\right)$$E38

involving the bed surface solid concentration ${\left(\varphi \right)}_b=1-\xi$.

The deposition rate $D_b$ is commonly related to the solid particles settling velocity in the mixture ${\omega }_{\mathit{sm}}$ and to the near-bed solid concentration in the flow ${\varphi }_{z\to z_b}$, which is usually related to the depth-averaged solid concentration in the flow as ${\varphi }_{z\to z_b}=\alpha \overline{\varphi }$, being $\alpha$ an adaptation or recovery coefficient, leading to $D_b=\alpha {\omega }_{\mathit{sm}}\overline{\varphi }$.

The settling velocity of the solid particles ${\omega }_{\mathit{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 ${\omega }_{\mathit{sm}}={\left(1-\overline{\varphi }\right)}^m{\omega }_s$, being ${\omega }_s$ the settling velocity of the sediment particles in clear water and $m$ a hindering empirical exponent depending on the Reynolds particle number (${\mathit{Re}}_p={\omega }_s{d}_s/\nu$, with $\nu$ the clear water kinematic viscosity) which usually takes values close to $m=4$. Therefore

$$D_b={\alpha }_D{\omega }_s\overline{\varphi }{\left(1-\overline{\varphi }\right)}^4$$E39

where ${\alpha }_D$ is a dimensionless parameter which accounts for mulple factor and requires calibration.

The erosion solid flux $E_b$ 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 ${\omega }_s$ and the near equilibrium concentration ${\varphi }_{z\to z_b}^{\ast }$. The near-bed equilibrium concentration is related to the depth-averaged equilibrium concentration $\overline{{\varphi }^{\ast }}$ as ${\varphi }_{z\to z_b}^{\ast }={\alpha }^{\ast }\overline{{\varphi }^{\ast }}$, being ${\alpha }^{\ast }$ an adaptation coefficient under equilibrium conditions. When equilibrium solid transport states are reached, the adaptation coefficients $\alpha$ and ${\alpha }^{\ast }$ coincide but in non-equilibrium states ${\alpha }^{\ast }\ne \alpha$ generally. Therefore, for the sake of simplicity, we assume that the vertical erosion rate $E_b$ is expressed as

$$E_b={\alpha }_E{\omega }_s\overline{{\varphi }^{\ast }}$$E40

where ${\alpha }_E$ is a dimensionless empirical parameter which requires calibration. The capacity solid concentration is usually computed as $\overline{{\varphi }^{\ast }}=q_s^{\ast }/\left(h\sqrt{{\overline{u}}^2+{\overline{v}}^2}\right)$, where $q_s^{\ast }$ 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 ${\varphi }^{\chi }$, referred to as buoyant solid concentration

$$r=\frac{\rho }{{\rho }_w}=1+{\varphi }^{\chi }\text{with}:{\varphi }^{\chi }=\sum _{p=1}^N\frac{{\rho }_{s,p}-{\rho }_w}{{\rho }_w}{\varphi }_p$$E41

where ${\rho }_{s,p}$ and ${\varphi }_p$ are the density and depth-averaged volumetric concentration of the $p$th solid phase respectively. Using (41), the equations forming the system can be recast as five conservation laws and rewritten in vector form as

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot \mathbf{E}\left(\mathbf{U}\right)={\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)+{\mathbf{S}}_{\tau }\left(\mathbf{U}\right)+{\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)$$E42

where $\mathbf{U}$ is the vector of conservative variables and $\mathbf{E}\left(\mathbf{U}\right)=\left(\mathbf{F}\left(\mathbf{U}\right),\mathbf{G}\left(\mathbf{U}\right)\right)$ are the convective fluxes along the $\mathbf{X}=\left(x,y\right)$ horizontal coordinates, respectively, expressed as

$$\mathbf{U}=\left(\begin{array}{c}\mathit{rh}\\ \mathit{rhu}\\ \mathit{rhv}\\ h{\varphi }^{\chi }\\ z_b\end{array}\right)\mathbf{F}\left(\mathbf{U}\right)=\left(\begin{array}{c}\mathit{rhu}\\ {\mathit{rhu}}^2+\frac{1}{2}{\mathit{grh}}^2\\ \mathit{\text{rhuv}}\\ \mathit{hu}{\varphi }^{\chi }\\ 0\end{array}\right)\mathbf{G}\left(\mathbf{U}\right)=\left(\begin{array}{c}\mathit{rhv}\\ \mathit{\text{rhuv}}\\ {\mathit{rhv}}^2+\frac{1}{2}{\mathit{grh}}^2\\ \mathit{hv}{\varphi }^{\chi }\\ 0\end{array}\right)$$E43

The vector ${\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)$ accounts for the momentum source term associated to the variation of the pressure force on the bed interface, whereas ${\mathbf{S}}_{\tau }\left(\mathbf{U}\right)$ is the momentum dissipation due to the basal resistance. Finally, the source term ${\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)$ accounts for the bulk mass exchange between the mixture flow and the bed layer.

$${\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)=\left(\begin{array}{c}0\\ -\mathit{grh}\frac{\partial z_b}{\partial x}\\ -\mathit{grh}\frac{\partial z_b}{\partial y}\\ 0\\ 0\end{array}\right){\mathbf{S}}_{\tau }\left(\mathbf{U}\right)=\left(\begin{array}{c}0\\ -\frac{{\tau }_b}{{\rho }_w}{n}_{\mathit{ux}}\\ -\frac{{\tau }_b}{{\rho }_w}{n}_{\mathit{uy}}\\ 0\\ 0\end{array}\right){\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)=\left(\begin{array}{c}-\frac{{\rho }_b}{{\rho }_w\left(1-\xi \right)}\sum _{p=1}^N{\left(D_b-E_b\right)}_p\\ 0\\ 0\\ -\sum _{p=1}^N\frac{{\rho }_{s,p}-{\rho }_w}{{\rho }_w}{\left(D_b-E_b\right)}_p\\ \frac{1}{1-\xi }\sum _{p=1}^N{\left(D_b-E_b\right)}_p\end{array}\right)$$E44

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 ${\mathrm{\Omega }}_i$. Applying the Gauss theorem leads to

$$\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}\underset{{\mathrm{\Omega }}_i}{\int }\mathbf{U}\mathrm{d\Omega }+\sum _{k=1}^{\mathrm{NE}}{\left(\mathbf{E}\cdot \mathbf{n}\right)}_k{l}_k=\underset{{\mathrm{\Omega }}_i}{\int }{\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)\mathrm{d\Omega }+\underset{{\mathrm{\Omega }}_i}{\int }{\mathbf{S}}_{\tau }\left(\mathbf{U}\right)\mathrm{d\Omega }+\underset{{\mathrm{\Omega }}_i}{\int }{\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)\mathrm{d\Omega }\end{array}$$E45

being $\mathit{NE}$ the number of edges for the $i$ cell, $\mathbf{n}={\left(n_x,n_y\right)}_k$ the outward unit normal vector, $l_k$ the length of the edge and hence ${\left(\mathbf{E}\cdot \mathbf{n}\right)}_k$ the value of the normal flux through the $k$th edge (Figure 4).

Figure 4.

Computational cells in triangular meshes and local coordinates at the $k$th 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

$${\left(\mathbf{E}\cdot \mathbf{n}\right)}_k={\left[\mathbf{F}\left(\mathbf{U}\right)n_x+\mathbf{G}\left(\mathbf{U}\right)n_y\right]}_k={\mathbf{R}}_k^{-\mathbf{1}}\mathbf{F}\left({\mathbf{R}}_k\mathbf{U}\right)$$E46

being $\mathbf{R}$ a rotation matrix which projects the global orthogonal framework $\mathbf{X}=\left(x,y\right)$ into the local framework $\widehat{\mathbf{X}}=\left(\widehat{x},\widehat{y}\right)$ of the $k$th cell edge (Figure 4), corresponding to the normal and the tangential directions to the edge respectively. The rotation matrix ${\mathbf{R}}_k$ and its inverse ${\mathbf{R}}_k^{-\mathbf{1}}$ are defined as

$${\mathbf{R}}_k={\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& n_x& n_y& 0& 0\\ 0& -n_y& n_x& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right)}_k{\mathbf{R}}_k^{-\mathbf{1}}={\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& n_x& -n_y& 0& 0\\ 0& n_y& n_x& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right)}_k$$E47

and the set of local conservative variables $\widehat{\mathbf{U}}$ and fluxes $\mathbf{F}\left(\widehat{\mathbf{U}}\right)$ is defined as

$$\widehat{\mathbf{U}}\equiv {\mathbf{R}}_k\mathbf{U}=\left(\begin{array}{c}\mathit{rh}\\ \mathit{rh}\widehat{u}\\ \mathit{rh}\widehat{v}\\ h{\varphi }^{\chi }\\ z_b\end{array}\right)\mathbf{F}{\left(\widehat{\mathbf{U}}\right)}_k\equiv \mathbf{F}\left({\mathbf{R}}_k\mathbf{U}\right)=\left(\begin{array}{c}\mathit{rh}\widehat{u}\\ \mathit{rh}{\widehat{u}}^2+\frac{1}{2}{g}_n{\mathit{rh}}^2\\ \mathit{rh}\widehat{u}\widehat{v}\\ h\widehat{u}{\varphi }^{\chi }\\ 0\end{array}\right)$$E48

where $\widehat{u}={\mathit{un}}_x+{\mathit{vn}}_y$ and $\widehat{v}=-{\mathit{un}}_y+{\mathit{vn}}_x$ are the components of the flow velocity $\widehat{\mathbf{u}}$ in the local framework.

The value of the fluxes through the $k$th cell edge can be augmented incorporating the contribution of the momentum source terms ${\mathbf{S}}_{\mathbf{b}}$ and ${\mathbf{S}}_{\tau }$ into the homogeneous normal fluxes $\mathbf{F}\left(\widehat{\mathbf{U}}\right)$ [25]. The momentum terms can be included within the local framework $\left(\widehat{x},\widehat{y}\right)$ using the spatial discretization

$$\underset{{\mathrm{\Omega }}_i}{\int }{\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)\mathrm{d\Omega }=\sum _{k=1}^{\mathrm{NE}}{\mathbf{R}}_k^{-\mathbf{1}}\mathbf{H}{\left(\widehat{\mathbf{U}}\right)}_k{l}_k\underset{{\mathrm{\Omega }}_i}{\int }{\mathbf{S}}_{\tau }\left(\mathbf{U}\right)\mathrm{d\Omega }=\sum _{k=1}^{\mathrm{NE}}{\mathbf{R}}_k^{-\mathbf{1}}\mathbf{T}{\left(\widehat{\mathbf{U}}\right)}_k{l}_k$$E49

where $\mathbf{H}{\left(\widehat{\mathbf{U}}\right)}_k$ and $\mathbf{T}{\left(\widehat{\mathbf{U}}\right)}_k$ are the integrated bed pressure [26] and basal resistance [27] throughout the $k$th cell edge, expressed in the local framework, being

$$\mathbf{H}{\left(\widehat{\mathbf{U}}\right)}_k={\left(0,-\mathit{grh}\mathrm{\Delta }{z}_b,0,0,0\right)}_k^T$$E50

$$\mathbf{T}{\left(\widehat{\mathbf{U}}\right)}_k={\left(0,-\frac{{\tau }_b}{{\rho }_w}\left(n_{\mathit{ux}}\mathrm{\Delta }x+n_{\mathit{uy}}\mathrm{\Delta }y\right),0,0,0\right)}_k^T$$E51

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

$$\mathcal{F}{\left(\widehat{\mathbf{U}}\right)}_k={\left(\mathbf{F}\left(\widehat{\mathbf{U}}\right)-\mathbf{H}\left(\widehat{\mathbf{U}}\right)-\mathbf{T}\left(\widehat{\mathbf{U}}\right)\right)}_k$$E52

The net exchange flux term ${\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)$ is discretized in space as

$${\int }_{{\mathrm{\Omega }}_i}{\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)d\mathrm{\Omega }\approx A_i{\mathbf{E}}_{\mathbf{b}}\left({\mathbf{U}}_i\right)$$E53

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

$$\begin{array}{l}{\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n-\frac{\mathrm{\Delta }t}{A_i}\sum _{k=1}^{\mathrm{NE}}{\mathbf{R}}_k^{-\mathbf{1}}\mathcal{F}{\left({\widehat{\mathbf{U}}}_i^n,{\widehat{\mathbf{U}}}_j^n\right)}_k^{\downarrow }{l}_k+\mathrm{\Delta }t{\mathbf{E}}_{\mathbf{b}}\left({\mathbf{U}}_i^n\right)\end{array}$$E54

being $\mathrm{\Delta }t=t^{n+1}-t^n$ the time step, $A_i$ the discrete cell area and $\mathcal{F}{\left({\widehat{\mathbf{U}}}_i^n,{\widehat{\mathbf{U}}}_j^n\right)}_k^{\downarrow }$ the numerical flux through the $k$th 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 $12\cdot {10}^6{m}^3$ mining waste tailings with a height of $70-80m$ and covering an area of $4.13\cdot {10}^5{m}^2$. 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 $S_0=0.0165m/m$. The area affected by the mud was $3.3\cdot {10}^6{m}^2$.

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 $22-26\mathrm{kN}/{\mathrm{m}}^3$ (Table 1).

Dam capacity	$12\cdot {10}^6{\mathrm{m}}^3$
Dam area	$4.13\cdot {10}^5{\mathrm{m}}^2$
Solid concent. ${\varphi }_0$	50%
Specific weight	$22-26\mathrm{kN}/{\mathrm{m}}^3$
Size distr.	Sand		Fines
Rel. content	38%		62%
Heavy metals	Fe	Al	Mn	Ti
Weight conc.	264.9 $\mathrm{mg}/\mathrm{g}$	10.8 $\mathrm{mg}/\mathrm{g}$	4.8 $\mathrm{mg}/\mathrm{g}$	0.5 $\mathrm{mg}/\mathrm{g}$

Table 1.

Brumadinho’s dam and tailings features.

The spatial domain is discretized using a unstructured triangular mesh with $6\cdot {10}^5$ 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 $\times$ 1 $\mathrm{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 ${\varphi }_0=0.5$ and normalized density $r\approx 2.25$. The deposition porosity $\xi$ 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 ${\theta }_{c,p}$ are estimated using the Egiazaroff formula [31]. Also, the solid transport capacity of the flow $q_s^{\ast }$ 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 ${\mu }_P$ as a function of the Manning’s roughness coefficient of the terrain ($n_b=0.065{\mathit{sm}}^{-1/3}$). The basal stability angle is estimated in ${\delta }_f=5^{°}$. The basal pore pressure is estimated as ${\mathcal{P}}_b=\left(1+{\mathcal{E}}_b\right){\rho }_w\mathit{gh}$, being ${\mathcal{E}}_b$ a basal pressure factor which tends to ${\mathcal{E}}_b\to \infty$ for totally fluidized materials. Therefore, the total frictional-turbulent basal shear stress is estimated as

$$\frac{{\tau }_b}{{\rho }_w}=\left(r-1-{\mathcal{E}}_b\right)\mathit{gh}+\mathit{rg}\frac{n_b^2}{h^{1/3}}\left(u^2+v^2\right)$$E55

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 ${\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)$ in (42) is neglected and the value of ${\mathcal{E}}_b$ is varied from ${\mathcal{E}}_b\to \infty$ (pure Newtonian turbulent behavior) to ${\mathcal{E}}_b=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 ${\mathcal{E}}_b$, the wave front loses velocity progressively as it moves downstream. As ${\mathcal{E}}_b$ 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 ${\mathcal{E}}_b=1.05$ allows to predict the observed arrival time ($32-47min$) of the tailings wave to the Paraopeba river at $8.5\mathrm{km}$ 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 ${\mathcal{E}}_b=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 $20\mathrm{m}$ in some zones during the first minutes. After this initial stage, the numerical results show that the dambreak wave is still higher than $10\mathrm{m}$ 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 ${\mathcal{E}}_b=1.05$ but the entrainment term ${\mathbf{E}}_{\mathbf{b}}\left(\mathbf{U}\right)$ in (42) is now considered. The relation ${\alpha }_E/{\alpha }_D$, which controls the balance between the erosion component $E_b$ and the deposition component $D_b$ of the vertical bed-flow exchange flux, is varied to ${\alpha }_E/{\alpha }_D=1\%$, ${\alpha }_E/{\alpha }_D=5\%$, ${\alpha }_E/{\alpha }_D=10\%$ and ${\alpha }_E/{\alpha }_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 ${\alpha }_E/{\alpha }_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 ${\alpha }_E/{\alpha }_D$ increases. Even, for ${\alpha }_E/{\alpha }_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 ${\alpha }_E/{\alpha }_D$ is increased.

Furthermore, the entrainment of bed materials also leads to the segregation of the solid phases along the flow. As ${\alpha }_E/{\alpha }_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, $E_p$ and $D_p$ respectively, along the tailings flow.

Figure 9.

Bulk solid concentration in the flow at $t=25min$ after the dam collapse for (a) ${\alpha }_E/{\alpha }_D=1\%$, (b) ${\alpha }_E/{\alpha }_D=5\%$, (c) ${\alpha }_E/{\alpha }_D=10\%$ and (d) ${\alpha }_E/{\alpha }_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 ${\alpha }_E/{\alpha }_D$ is mainly appreciated in the increment of the flow rate peak from $7500{\mathrm{m}}^3/\mathrm{s}$ for the fixed-bed simulation to almost $15,000{\mathrm{m}}^3/\mathrm{s}$ for the movable bed case with ${\alpha }_E/{\alpha }_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 ${\alpha }_E/{\alpha }_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 time	Speed-up
CPU Intel i7-10700F	59.43 h	—
Nvidia Tesla K40c	1.79 h	× 33
Nvidia GTX 1080 Ti	58.2 min	× 61
Nvidia Tesla V100	26.5 min	× 135
Nvidia A100 40GB	12.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.

Acknowledgments

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.