Brumadinho’s dam and tailings features.

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

### Keywords

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

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

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

where

In order to develop a shallow-type depth-averaged mathematical model, the Navier-Stokes system (1) and (2) is integrated between the free surface

being

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

being

Regarding the momentum depth integration, the volumetric force components in (1) and (2) are considered for the sake of simplicity as

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

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

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

and, similarly, for the

where

### 2.2 Depth-integrated solid transport equation

The solid-phase transport process is governed by

where

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

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

where

Furthermore, assuming that the volumetric solid concentration in the bed layer is

### 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

where

Considering simple shear state along the flow direction, the velocity vector throughout the flow column is expressed as

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

being

The frictional non-linear viscoplastic model considers a Coulomb-Terzaghi linear relation between the effective normal stress

where

Using (28), and considering a plastic viscosity

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

where

being

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

Integrating (35) throughout the flow column allows to obtain the basal shear stress

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

### 2.4 Net mass exchange between bed and flow layers

The net volumetric exchange

involving the bed surface solid concentration

The deposition rate

The settling velocity of the solid particles

where

The erosion solid flux

where

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

Considering a multi-grain mixture flow, the resulting system is composed by

where

where

The vector

### 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

being

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

and the set of local conservative variables

where

The value of the fluxes through the

where

Using (46) and (49), the augmented flux at the

The net exchange flux term

and, assuming a piecewise uniform representation of the variables at the

being

### 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

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

Dam capacity | ||||
---|---|---|---|---|

Dam area | ||||

Solid concent. | 50% | |||

Specific weight | ||||

Size distr. | Sand | Fines | ||

Rel. content | 38% | 62% | ||

Heavy metals | Fe | Al | Mn | Ti |

Weight conc. | 264.9 | 10.8 | 4.8 | 0.5 |

The spatial domain is discretized using a unstructured triangular mesh with

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

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

Figure 7 shows the mud depth at

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

Furthermore, the entrainment of bed materials also leads to the segregation of the solid phases along the flow. As

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

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 |

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

## References

- 1.
Iverson RM. The physics of debris flows. Reviews of Geophysics. 1997; 35 (3):245-296 - 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.
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.
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.
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.
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.
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.
Wu W. Computational River Dynamics. NetLibrary, Inc. London: CRC Press; 2007 - 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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Jakob M, Hungr O. Debris-Flow Hazards and Related Phenomena. Springer Praxis Books. Springer Berlin Heidelberg; 2005 - 20.
Iverson RM, Vallance JW. New views of granular mass flows. Geology. 2001; 29 (2):115-118 - 21.
Berti M, Simoni A. Experimental evidences and numerical modeling of debris flow initiated by channel runoff. Landslides. 2005 - 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.
Richardson JF, Zaki WN. Sedimentation and fluidisation: Part I. Chemical Engineering Research and Design. 1997; 75 :82-100 - 24.
Godlewski E, Raviart P-A. Numerical Approximation of Hyperbolic Systems of Conservation Laws. New York: Springer-Verlag; 1996 - 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.
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.
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.
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.
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.
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.
Egiazaroff IV. Calculation of nonuniform sediment concentrations. Proceedings of the American Society of Civil Engineers. 1965; 91 :225-247