Introduction

According to the Special Report: Global Warming 1.5 °C prepared by the IPCC, 2018[1], the temperature of the Earth's surface could increase more than 2 °C above preindustrial levels by the end of the century, which would cause nonlinear alterations in the patterns of natural hazards in terms of their geographical distribution, frequency and intensity. This is due to changes in the volume, intensity and frequency of rainfall, which will cause droughts, floods, cyclones, tropical storms and increasingly frequent mud and debris flows. As a consequence, there will be an increase in deaths, losses and material damage, which will exceed the capacity of risk mitigation and response mechanisms, especially in poor and developing countries [2].

The Global Assessment Report on Disaster Risk Reduction of the United Nations, 2009 [3] shows that the risk of disasters of meteorological origin is already disproportionately concentrated in developing countries, specifically in their poorest populations. The risk of natural disasters is related not only to the threat of the disasters but also to the vulnerability of the populations; in turn, this is associated with their livelihoods, which are unable to sustain minimum levels of quality of life, especially in rural areas where houses and other infrastructures are built with materials and techniques that are not very resistant to threats; additionally, weak government in urban centers and inadequate territorial planning increases the exposure of populated centers to extreme natural events. Thus, a feedback relationship is generated between losses due to natural disasters associated with extreme events, which are exacerbated by climate change and increased poverty, which in turn leads to greater vulnerability of populations to natural hazards.

To reduce climate risk, as well as its socioeconomic impacts, it is necessary to implement prevention measures such as adaptation strategies such that long-term risk management allows the possibility of losses and damage caused by the effects of climate to be reduced. The evaluation of the threat is the first step to achieve this goal, and it can be carried out from the characterization of its location, intensity or magnitude and its frequency or probability of occurrence. The magnitude can be expressed qualitatively, in terms of scales of high, medium and low, or quantitatively, in terms of the affected area, depth and flow velocity. The frequency can be associated with different return periods or probabilities of occurrence.

Mud and debris flows are sequential and combined natural threats because, on the one hand, they are hydrometeorological threats since they are of atmospheric and hydrological origin and, in turn, are closely related to geological threats such as landslides [4]; however, there is still no consensus on their definition. According to the Peril Classification and Hazard Glossary [5], mud and debris flows are a type of landslide that occurs when heavy rains send large amounts of vegetation, mud, or rocks downhill by gravitational forces, whereas for the Disaster Inventory System belonging to the Sendai Framework for disaster risk reduction, it is a violent flow of water in a basin, sometimes reported as sudden rising, or as a torrent that can transport tree trunks and/or abundant sediments from fine to rock blocks; this flow can be generated by rain, by the rupturing of dams or by abundant landslides in a basin [6]. Although mud and debris flows are increasingly frequent phenomena that have generated great economic and life losses, they are still rarely studied, so threat assessment through modeling, including both climatic and geological conditioning factors in a joint and coupled way, can become a very useful tool for risk management and adaptation to climate change.

Conditioning factors of mud and debris flows

Mud and debris flows typically occur in steep basins with little vegetation and in areas that are highly susceptible to landslides and that are inherently unstable. A common way mud and debris flows form is that during intense rains, usually preceded by dry season, surface faults or landslides are generated due to soil saturation with little cohesion are generated. The subsequent detachment of this soil contributes materials, from fine to thick, to the beds of nearby channels, and these materials flow immediately through the steepest sections and stop in areas where the slope decreases. As the rain continues or becomes more intense, the material can accumulate in the less steep sections, which can generate temporary dams [7], which can cause the accumulation of water in natural reservoirs that break due to a gradual increase in pore pressure or a sudden increase in drainage.

Although a less frequent way that debris and mud flows can also be generated by basal erosion of the riverbed. This process of lateral and bottom scour can also occur during an event caused by another conditioning factor, generating a greater contribution of sediment to the flow and possibly altering the rheological behavior. Finally, some events can be triggered by the melting of snow peaks and pyroclastic flows or by a shallow seismic event of great magnitude that, in combination with a period of intense rains, causes generalized landslides throughout the basin [8].

Classification of mud and debris flows

In general, flows can be classified into 3 types: mud flows, hyperconcentrated flows and debris flows. Mud flows are flows with high concentrations of fine sediments, specifically silt and clay, with high viscosity and shear stress [9]; this type of flow can transport large blocks of rock during flotation. Mud flows behave as a single mass, and when they settle, they do not do so in a stratified way (deposit first the thickest and then the fine sediments), and these sediments form lateral deposits along the main channel [10].

Second, there are hyperconcentrated sediment flows, also called mud floods, which can contain a concentration of sediments of up to 60% between fine and coarse; the main sediments are sand, gravel, pebbles and blocks, and there is little cohesion between them. In these flows, the sediments move partially as a bottom load and are partially suspended, and sedimentation of the flow occurs in a classified way, first the thickest and then the finest ones [8]. Hyperconcentrated flows are turbulent, and flow resistance depends on roughness, as is the case for clear water flows [9].

Finally, debris flows contain sediment concentrations greater than 60%, and the predominant sediment is coarse solids, mostly of the gravel- type. They can transport large blocks in suspension [10], and energy dissipation occurs mainly because of the high collision between sediment particles [9].

Flow classification does not depend only on sediment concentration but also on the predominant type of sediment, so in some cases, its classification is more difficult, and there is no single criterion to perform it. Table 1 shows the references of some proposals made by various authors for the classification of this type of flow.

Rheology

Rheological models describe the behavior of fluids in terms of their viscosity, plasticity and elasticity, which are properties that are described by shear stress. Considering that mud and debris flows do not necessarily behave like Newtonian fluids do, mainly because of their high concentration of sediments, it is convenient to use different rheological relationships for their modeling. Next, the mathematical structure of some of the most commonly used rheologies for the modeling of mud and debris flows is presented.

Newtonian fluids

These fluids follow Newton's law of viscosity, in which there is a linear relationship between the viscosity 𝜇 and the shear stress 𝜏, the latter represented as the variation in the strain rate ν with respect to depth 𝑦.

In some cases, debris flows have been modeled as Newtonian fluids by implementing the Manning equation and assigning a “pseudo-Manning” value n, which is much greater than the values commonly used above, to simulate the dissipation of energy generated by collisions between sediment particles. However, this is a methodology with little physical basis since the mechanisms of energy dissipation by friction of a debris flow are very different from those that act in turbulent water flows for which the Manning equation was established [15].

Non-Newtonian fluids

Non-Newtonian fluids can be represented by adding components of shear forces; the forces included in the equation depend on the type of fluid, be it turbulent, pseudoplastic or viscous, and can include the yield stress 𝜏_𝑐 for viscous fluids, the cohesive limit stress 𝜏 _𝑜 , the Mohr Coulomb stress 𝜏_𝑚𝑐, the viscous shear stress 𝜏_𝑣, the turbulent stress 𝜏_𝑡 or the dispersive shear stress 𝜏_𝑑, among others [9]. The Bingham model and the Herschel Bulkley model are some of the most widely used rheologies to characterize viscous flows such as mud flows, whereas the Voellmy model has been used more for debris flows. These and other rheologies included in the simulation models of flow propagation in the flood valley are briefly described below.

Bingham model

The Bingham model introduced by Bingham & Green, 1919 [16] is used to characterize viscous fluids, which behave as a solid until a minimum shear stress or yield stress is exceeded 𝜏_𝑜. From that point on, the relationship between the deformation dv/dy and the shear stress 𝜏 can be linear or nonlinear. A characteristic of this type of fluid is that it can flow on flat surfaces, unlike nonviscous fluids, which can only flow through a slope.

Herschel Bulkley model

The model of Herschel & Bulkley, 1926 [17], very similar to the Bingham model, is also used to represent the rheological behavior of viscous fluids, including a yield shear stress 𝜏_𝑜, a consistency factor 𝑘, representing the viscosity of the fluid and a flow index 𝑛 that indicates the degree to which the fluid is thinning or thickening, that is, it indicates the linear or nonlinear behavior of the flow deformation.

Quadratic model

The quadratic model introduced by O'Brien et al., 1985 [18] divides the yield stress 𝜏_𝑐 in cohesive yield effort 𝜏_𝑜and the effort of Mohr–Coulomb 𝜏_𝑚𝑐. This model also includes the dispersive and turbulent shear stresses represented by the coefficient of inertia stress 𝐶, which in turn depends on the density of the mixture 𝜌_𝑚, the Prandtl mix length 𝑙_𝑚, the diameter of the sediment 𝑑_𝑠 and a function of the volumetric concentration of sediment 𝐶𝑣 [19].

Bagnold model

The model of Bagnold, 1954 [20] with dilating behavior (the fluid flows more easily at higher applied shear stresses) differentiates three types of flows: the macroviscous regime, transitional regime and inertial granular regime. The regime to which the flow belongs can be determined by the Bagnold number 𝐵_𝑎, which depends on the diameter of the sediment 𝑐, the volumetric concentration of the sediment 𝐶, the maximum sediment concentration 𝐶₀, the dynamic viscosity 𝜇 and the density of the sediment particles 𝜌_𝑠. For Bagnold numbers less than 40, the corresponding regime is the macroviscous regime, whereas for Bagnold numbers greater than 450, the regime corresponds to the granular-inertia regime, and in intermediate values, the transition regime is found.

According to the above, an equation describes the shear stress as a function of the deformation of the fluid for each regime, where 𝑎_𝑣 and 𝑎_𝑖 are experimental constants and where 𝛼₁ is the dynamic friction angle, which is different from the internal friction angle.

Macroviscous regime

Granular-inertial regime

Takahashi model

Takahashi, 1978 [21] developed a model from Bagnold's model, introducing variations to the parameter 𝑎_𝑖, and incorporating an equation that allows calculating the value of the dynamic friction angle 𝛼₁ as a function of the volumetric concentration of sediments 𝐶, the maximum concentration of sediment 𝐶₀ and the angle of internal friction 𝜙.

In 1977, Takahashi developed an equation to determine the equilibrium concentration or maximum concentration of sediments in a debris flow, which depends on the density of the water 𝜌_𝑤, sediment density 𝜌_𝑠, angle of inclination of the terrain 𝜃 and angle of internal friction 𝜙 [22]

Mohr-Coulomb model

The Mohr‒Coulomb model was proposed by Johnson & Kehle, 1972 [23] and is based on the Coulomb model for soil movement and includes parameters of the Bingham model. In this model, yield stress is included as a function of flow cohesion 𝐶 and includes normal effort 𝜎_𝑛 related to the angle of internal friction 𝜙.

Voellmy model

The Voellmy, 1955 [24] model was initially developed for snow avalanches and was later implemented by Körner, 1976 [25] for rock avalanches. This model can be used for a wide range of sediment concentrations and takes into account the Chezy coefficient of friction. 𝐶_𝑧, flux density 𝜌, angle of internal friction 𝜙, angle 𝜃 defined by gravity, flow velocity 𝑉 and pressure gradient.

Methodologies for modeling mud and debris flows

There are multiple approaches for the modeling of mud and debris flows, from the way triggering factors are modeled to the propagation of the flow, including empirical approaches to physically based and highly complex models in their mathematical scheme. Some of the methodological approaches are focused on a simplified regional or local scale, in which hydrological modeling or trigger mechanisms are not necessarily carried out; in contrast, flow volumes are estimated through geomorphological analysis [26], interpretation of satellite images and collection of field information, and modeling focuses solely on flow propagation [27]. Other approaches attempt to simulate in detail all the processes involved from the antecedent rain to the detonation of landslides and all the processes that occurred during propagation until their deposition in the flood valley.

In general, the modeling of mud and debris flows can be divided into 4 phases, some of which can be approached in a simplified way depending on the goal of study: 1. characterization of the basin and triggering factors, 2. hydrological modeling, 3. modeling of triggering mechanisms, and 4. modeling of flow propagation. The first phase is carried out mainly by collecting field information and analyzing satellite and cartographic information.

In cases in which the triggering factors are directly related to rainfall, hydrological analysis is a necessary step for the subsequent modeling of triggering mechanisms such as landslides and dam failure, in addition to providing information on the liquid flows of the avenue. Finally, flow propagation can be carried out through different types of models, including empirical approaches to one-dimensional, two-dimensional models of 1, 2 or 3 phases, and incorporating processes such as watercourse erosion and changes in the rheological characteristics of the flow. Table 2 presents some references to the modeling of mud and debris flows, in which different methodological approaches and computational models are used, and each phase or step in the modeling is detailed later, as well as the different options available for its development.

Hydrological modeling

Hydrological rain-runoff modeling not only allows the determination of liquid flows but also is an essential input for the modeling of triggering mechanisms such as landslides or the failure of natural dams. The hydrological models can be classified by the conceptualization of the modeled processes, the spatial distribution and the temporal variability. According to the conceptualization of the processes, 3 types of models are distinguished: empirical or "black box" models, in which the model is based on an input‒output relationship without describing the behavior of each of the processes (this category includes infiltration models such as the Green Ampt or SCS curve number method); conceptual or "gray box" models, in which some of the internal relationships or processes are known and in which physical laws are considered in a simplified way; and finally, physically based or “white box” models, which are those in which all the processes involved are known and are simulated in a detailed way for which they are supported by physical laws [44].

Furthermore, according to their spatial representation, hydrological models can be aggregated, semidistributed and distributed. In the aggregated models, the basin is represented as a single element with homogeneous and constant characteristics, which are valid for very small basins. For larger basins with soil characteristics, slopes, covers and other variable physiographic characteristics, semidistributed models can be used, in which zones with similar behavior or characteristics are distinguished, either by subbasins or by smaller units called units of hydrological response (HRUs). Finally, for basins with very heterogeneous characteristics, to achieve greater representativeness of the characteristics of the basin, distributed models can be used, in which the basin is divided into much smaller elements of uniform geometry called grids or cells, each of which includes the specific characteristics of the soil, humidity, coverage, and slope, among other features [44].

Hydrological models are classified according to their temporal variability or time scale into event models and continuous time models. Storm histograms hietograms are simulated that can be obtained from pluviographic records, forecast meteorological models [40] or synthetic storms built with methods such as alternate blocks or Basto-Salazar curves [37], and their time intervals are on the order of minutes or hours. These rains can be associated with different return periods, events that have occurred in the past, or thresholds that trigger landslides. With these models, flood hydrographs are obtained that are included as liquid flows in the flow propagation modeling. In continuous time models, a period of time that can be from days to years is simulated, with a time interval usually daily. Continuous time models allow the determination of the characteristics of the antecedent humidity that increases triggering events such as landslides and the generation of natural dams, as well as the determination of rainfall on the day of the event; however, if a continuous time model is used, a synthetic flood hydrograph must also be generated via methods such as the triangular hydrograph of the United States Soil Conservation Service [39] to incorporate it as liquid flow in the flow propagation model.

Within the wide variety of hydrological models, some have the advantage of explicitly calculating the subsurface flow; this is the type of flow that increases the pore pressure, triggering instability in the slopes. With a single model, the value of the rainfall–runoff flow of the event can be obtained, as can the antecedent subsurface water flows that can trigger landslides, without the need to develop separate infiltration models, such as the model of O'Loughlin, 1986 [35] or the Richards model, to determine the level of soil saturation or level of the water table that generates soil saturation, nor is it necessary to later couple these models with the landslide models. The TOPMODEL model developed by Beven & Kirkby, 1979 [45] calculates subsurface flow as a function of hydraulic transmissivity 𝑇_𝑜, a parameter that describes subsurface flow behavior with depth 𝑚, flow accumulation area 𝑎 and slope 𝑆 [46].

Table 3 presents some hydrological models that can be used for the modeling of mud and debris flows, as well as their classification, modeled processes and main parameters. The choice of the model depends on the temporal scale of available climatological information, the characteristics of the event being simulated and the approach of the modeling, considering whether it is necessary to represent the entire process of generating the event from the days prior to the event or if the model is intended only to simulate the effects of the flow in the flood valley.

Landslide modeling

Landslides are the most recurrent triggering mechanism for mud and debris flows, generated by either short and very intense rainfall or by moderate-intensity rainfall that lasts for several days, especially when preceded by dry periods during which the cohesion of the river has decreased. This is why a comprehensive landslide threat analysis is essential for determining the magnitude of a mud and debris flow event, as well as the type of flow that would be generated in future events. Landslide hazard analysis can be carried out in different ways by implementing statistical methods, heuristic methods, deterministic methods, or a combination of these methods. The method and/or model to be selected depends on the quantity and quality of information available, as well as the level of detail required for the analysis.

For its part, the heuristic method performs an analysis of cartographic and geomorphological information and weighs the factors causing instability according to the criteria of experts; the method can be complemented by statistical methods of inventory analysis of mass removal events. This type of methodology has been widely used for risk analysis due to mud and debris flows, as is the case of Sepúlveda et al., 2016 [38], who analyzed susceptibility to landslides from maps of events, lithology, soil cover, slopes, curvature and soil thickness; Kain et al., 2018 [41], who mapped erosion and deposition patterns using remote sensing images; and Bertoldi et al., 2012 [33], who used it in combination with field studies and the developed magnitude‒frequency relationships to evaluate future scenarios from empirical geomorphological relationships.

Furthermore, there are deterministic models that are physically based and generally determine the threat from landslides by calculating the factor of safety. Importantly, although these models can have a high level of detail in terms of the amount of slipped mass (which allows the determination of the concentration of sediments and, in turn, the type of flow and its behavior with greater certainty), they are highly complex. Additionally, some of these models include vertical infiltration models with complicated mathematical solutions that require high computational effort and are highly susceptible to uncertainty, which in turn is a function of the quality of the field information. This is why when a study of this type is carried out, it must be weighted according to the level of detail required if a heuristic–statistical analysis is sufficient or if, on the contrary, it is justified to invest greater computational and economic effort in the gathering of field information for the construction of a deterministic model, which must also be calibrated with historical data. Next, the mathematical structure of some of the deterministic models used in various studies for the simulation of mud and debris flows with landslides as a conditioning factor generating the event is presented.

SHALSTAB model

This model is based on the Mohr‒Coulomb failure law, in which the shear stress necessary for the slope to fail is equal to the resistance generated by the cohesion of the soil 𝐶 and frictional resistance due to normal stress. The model calculates the soil saturation in terms of the relationship between the subsurface flow 𝑞 and the hydraulic transmissivity of the soil 𝑇 required for the ground to fail. This relationship depends on the slope 𝜃, the accumulated drainage area 𝑎 per unit width 𝑏, the angle of internal friction of the material 𝜙, the soil density 𝜌, the density of water 𝜌_𝑤, the thickness of the soil 𝑧 and cohesion 𝐶. Stability can be classified from unconditionally stable saturated (the soil will not fail under any precipitation) to unconditionally unstable dry (the soil will fail even if it is completely dry), which possibly are rocky outcrops [59].

Infinite slope model

The infinite slope model combines the Mohr‒Coulomb failure mechanism with the saturation of the soil generated by the horizontal flow of water or subsurface flow, determining the stability by calculating the factor of safety. The model has as an input parameter the value of the water table 𝐷_𝑤, which can be calculated by means of subsurface flow from a hydrological or infiltration model [45]

Bishop’s method - limit equilibrium

The limit equilibrium method of Bishop, 1955 [60] is used for failures of the circular type, in which the failure surface is divided into a series of vertical slices. The simplified Bishop method assumes that the lateral forces between the slices are horizontal and that the shear forces are neglected; calculating the safety factor iteratively [10] considers the equilibrium of moments in relation to the center of the circle for each one. of the slices. Additionally, this method allows the effects induced by the evolution of the limits where flow propagation occurs, that is, the failure of the riverbed itself [36].

STEP TRAMM model

The STEP TRAM model allows modeling of the chain reactions of small faults that generate landslides, in addition to the mobilization of the material to the drains. By means of the landslide modulus, the table of water in the soil column that affects the resistance of the soil is calculated, considering that an imbalance of driving forces and resistance must exist [61]. The tensile strength 𝜏_𝑡 is calculated using the Mohr−Coulomb criterion; in this case, the total cohesion 𝐶, which differentiates between the cohesion of the roots 𝐶_𝑟 and the cohesion of the soil 𝐶_𝑠, also includes the term capillary pressure ℎ_𝑐 and a coefficient 𝜒, which defines the relationship between the capillary force and capillary pressure.

The compressive strength 𝜏_𝑐𝑜 is calculated by the following expression:

Once the imbalance of forces generated by the small faults is calculated, the model represents the mechanical interactions of the soil with its neighboring columns by means of conceptual mechanical links called FBMs or “fiber bundles” that are composed of numerous mechanical elements called fibers that break to a predefined threshold redistributing the load on the other fibers, triggering local failure of the entire column [61].

TRIGRS model

The transient rainfall infiltration and grid-based regional slope-stability model (TRIGRS) was developed by the United States Geological Survey (USGS) to determine progressive changes in pore pressure and changes in the factor of safety due to infiltration of rain, both for wet initial conditions and for unsaturated initial conditions [62].

The model for wet initial conditions is based on the solution of Iverson, 2000 [63] of the Richards vertical infiltration equation, in which the infiltration depends on the depth of the water table and a constant infiltration rate, with a vertical hydraulic gradient, which is also constant and depends on the infiltration rate and the saturated hydraulic conductivity [62].

Furthermore, for initial unsaturated conditions, the model divides the soil into two layers: the saturated zone up to the water table and an unsaturated upper layer up to the soil surface. A part of the infiltrated water reaches the saturated zone that raises the water table, generating a load that propagates downward as waves of diffusive pressure increase the pore pressure; this process is described by a one-dimensional form of the Richards equation [62].

Dam failure modeling

Another triggering mechanism for recurrent mud and debris flows is the breaking of natural dams, which are generally formed by landslides, and the accumulation of plant material, in which water accumulates for several days before breaking due to excess pore pressure, or by erosion in the dam by tubing or overflow. Importantly, the simulation of the dam failure must be complemented with an adequate identification of the possible dam sites by analyzing the topography of the channel, as well as possible landslides of great magnitude that could generate obstruction of the channel.

For the simulation of dam failure, there are multiple types of software that allow the determination of both the liquid and solid hydrographs of failure, including DEBRIF1D, NWS-DAMBRK, BOSS-DAMBRK, and BREACH, among others. For its part, DEBRIF1D was first used as a model to simulate the dynamics of snow avalanches; later, in 2000, it was used as a model for mud and debris flows detonated by dam failure. Četina et al., 2006 [28] implemented it in the modeling of two consecutive debris flows that occurred due to the failure of a natural dam and the additional contribution of flow due to the obstruction of a bridge in the city of Gorenji Log, Slovenia. This model solves the conservation equations of momentum and the continuity equation by means of a finite difference scheme and can consider the energy loss generated by the wave front. An advantage of the model is that it allows the initial flow hydrograph to be determined just at the instant after the dam collapses. In addition, the BREACH model was implemented by Páez, 2016 [39], to simulate the failure of a natural dam that generated a mud flow in Utica, Cundinamarca; it is a physically based model designed to simulate the breaking of earth dams and has the ability to simulate the outlet flow of the spillway if it exists, the flow rate of the rupture either by tubing or overflow and the transport of sediments in the gap generated by erosion by overflow on the crest of the dam.

Flow propagation modeling

Flow propagation is the final step necessary in determining the threat from mud and debris flows in terms of flooded area, depth and flow velocity. For this purpose, there is a wide variety of options from simplified empirical relationships that calculate the travel distance and flow volume to one- and two-dimensional models that involve 1, 2 or 3 phases (liquid and solid) and include rheological relationships and channel erosion. The choice of the model to be used depends on the level of detail to be obtained in the results and the scale applied, whether regional or local, since the higher the level of detail is, the higher the field information requirements, computational effort and associated costs.

Empirical relationships

There are empirical relationships, such as those presented in Table 4, that allow the determination of different flow parameters, such as the length of the flow, the length of the flood valley or the volume of the flow. These types of relationships are established using large datasets of debris flows that can be analyzed while considering the characteristics of the basin and implemented to generate preliminary hazard maps or complement field analyses [64], and they can also be built from machine learning algorithms [65][66].

Hürlimann et al. (2008) [27] implemented several empirical relationships to generate hazard maps for the Guingueta Basin, Spain, for which they used the longitudinal profile of the flow path and the volume of the event. Furthermore, Bertoldi et al., 2012 [33] applied the methodology of D'Agostino & Marchi, 2003 [67] to determine approximate volumes of debris flows in the Bondes and Tracuit basins in Switzerland, supported by geomorphological data, geological analysis, analysis of aerial photographs and field inspection, this as a complement to a two-dimensional model for the delimitation of threat areas. Tin such a way that those generated by the two-dimensional models and medium and low threat were delimited as areas of high and very high threat those calculated through empirical relationships.

Some empirical relationships have also been incorporated into computational models, such as Bernard et al., 2020 [68] and the DFLOWZ model modified by Berti & Simoni, 2007 [69] from the LAHARZ algorithm developed by Schilling, 1998 [70], which uses empirical relationships to simulate pyroclastic flows. The DFLOWZ model also simulates the unconfined deposition of the flow, which done based on the use of two empirical relationships for the cross-sectional area of the flow 𝐴 and the planimetric flood area 𝐴_𝑓 from the flow volume together with uncertainty adjustment factors 𝑎 and 𝑏 [71]. The above, which is based on a digital elevation model and a flow path, allows the delimitation of deposition areas for a range of probable events.

This type of approximation, although simplified, can include great uncertainty, since the determination of the starting point of the flow and the deposition zone can be difficult; like the intensity of the flow, the depth and speed can be quantified only indirectly, and the rheological behavior of the flow cannot be considered. Therefore, this type of method must be complemented with good field exploration and geomorphological and topographic analysis to determine with some certainty the volume and path of the flow; similarly, these relationships are clearly limited to providing an initial and very general estimate of potential threat areas.

One-dimensional models

The empirical relationships do not consider the rheology of the flow, so a simple way to incorporate it into the propagation analysis is through numerical or analytical one dimensional models, which allow simulation of the total travel distance of the flow, as well as the depth and velocity at each point of the trajectory. Some of them have incorporated rheology through the Voellmy frictional–turbulent model, either with an analytical or numerical solution. The one-dimensional Voellmy analytical model (Equation 20) incorporates the curvilinear distance that the flow travels 𝑠, the slope of the channel 𝜃, and two coefficients of friction; the first is called the coefficient of sliding friction 𝜇_𝑚, and the second turbulence coefficient is called the mass‒drag relationship, which was originally applied to snow avalanches by Perla et al., 1980 [80].

Furthermore, the Voellmy numerical model (Equation 21) can be expressed in terms of the total friction slope 𝑆_𝑓, which is a function of the dry friction coefficient 𝜇_𝑑, the pseudo Chezy coefficient 𝐶 and the hydraulic radius 𝑅.

The previous one-dimensional models have been used to delimit threat zones on the basis of a previous multidisciplinary analysis that includes geomorphological analysis, interpretation of satellite images, analysis of historical events, determination of susceptibility maps due to mass movements and estimation of event volumes from geomorphological relationships and field studies; the obtained results are quite close to the reality of historical events [27] [81] [36].

One-dimensional models can be good tools for assessing the risk associated with mud and debris flow events, incorporating the behavior of the fluid. However, if the desired analysis is one in which the possible damage associated with the magnitude of the event is known with some certainty, this type of model is insufficient, and a complete modeling of factors is advisable. Instead, is it preferable to use conditioning factors and flow propagation with two-dimensional models that consider the behavior of the flow in the flood valley and the entire process of deposition.

Two-dimensional hydraulic models

Two-dimensional hydraulic models make it possible to determine the flooded area, as well as the depths and velocities of the flow at any point in the riverbed or alluvial valley, which is why these models are very useful for developing hazard maps for different return periods on the basis of depth and flow velocity on the basis of the definition of criteria or thresholds of affectation according to these variables [27].

Among the two-dimensional hydraulic models are those whose mathematical approximation is that of the hydrodynamic transport of a single phase; that is, all the flows behave as a single mass, and the model is based on the solution of the Navier–Stokes equations or shallow water equations integrated in the vertical direction and solved under a scheme of elements or finite volumes for a triangular or rectangular mesh that can be of variable size. Shallow water equations are described by means of the conservation of momentum and the conservation of mass equations in two dimensions.

Although many two-dimensional models exist that solve the Navier–Stokes equations, not all of them include rheological models that characterize the non-Newtonian behavior of a mud or debris flow. Therefore, when the model to be used for the simulation of flow propagation is chosen, the type of flow and the most appropriate rheology must be considered because there may be substantial differences in the results of the propagation model. This was reported by Páez, 2016 [39], who compared important differences in flow depths and velocities for the same event for the Bingham, quadratic and pseudo-Manning models. The rheology of the flow could make it convenient to compare the results of the modeling for different rheological models, since the uncertainty also generates uncertainty with regard to the possible damage that a population affected by a mud or debris flow could receive.

Furthermore, there are two-dimensional models of two phases, such as OpenLISEM [42], and three phases, such as AVAFLOW [82]. This type of mathematical approximation divides the flow into the liquid phase and one- or two-phase solids, for which independent mass and momentum conservation equations are derived for each phase [83]. According to Pudasaini & Mergili, 2019 [84], the fluid phase composed mainly of water and very fine particles of the type silt and clays or colloids is modeled with a viscoplastic model such as Herschel–Bulkley or Bingham through a constant yielding effort. Furthermore, if it is a two-phase flow, the solid phase can be represented by the Mohr‒Coulomb model, whereas if a third phase is added, there would be two solid phases, one consisting of fine solids and the other consisting of coarse solids. Fine solids include sands and fine gravels and are modeled with a viscoplastic Coulomb rheology dependent on pressure and shear stress because the elastic limit could depend on the pressure and friction of the material. While coarse solids include coarse gravels and boulders, which are modeled with plastic Mohr‒Coulomb laws regardless of the cutting speed, since coarse particles do not have viscous interactions, the friction force depends on the normal load and friction resulting from collisions between the particles [84].

In Table 5, some two-dimensional hydraulic models are presented, as well as the rheologies they include and the modeled processes, among which some models include the transport of bottom and suspended sediments; however, this type of modeling can only be used as a complement to the pseudo-Manning model, in which the flow is modeled as a water flow with a high sediment load and a higher Manning n than usual to represent friction losses due to the interaction and collision of the particles.

Modeling the propagation of a mud or debris flow has great challenges because its behavior has characteristics that are difficult to simulate, such as the generation of lateral erosion during the transit of the event, the change in the rheological properties of the flow as more sediment is incorporated and the generation of pulses that can stop and start moving again. For the simulation of this last condition, the FLATModel includes a “stop and go” algorithm, which defines whether a cell is at total rest according to the dynamics of its neighboring cells, considering that a cell is in motion if the following three conditions are met: 1. The cell was in motion in the time interval just before, 2. some of the neighboring cells were in motion in the time interval just before, and 3. the geometric slope of the flow is greater than the angle of internal friction [29]. This last criterion is very important since the inclination of the final deposit is similar to the angle of internal friction when the Voellmy model is used. However, the initial angle would be that of completely dry material; as a result, as the flow progresses, the value may vary. To correct this value, FLATModel increases the angle value if the flow velocity is less than a limit and varies linearly from the initial value for dry material [30].

Basal erosion modeling

In some circumstances in which the slope and the type of material favor it, the shear stress generated by the torrential flood on the bed is sufficient to incorporate material from the bed into the flow, from a few cubic meters to a volume of ten or more times the volume initially mobilized in the avenue [92]. Even if the carryover of the material is very high, such that the equilibrium or maximum concentration described by Takahashi, 1977 [22] (Equation 10) is greater than 0.4, this carryover of material from the bed can become the trigger for a torrential flood.

A first approximation of the depth of basal erosion was made by Takahashi, 1991 [93], in which it is calculated from the concentration of sediments in the flow 𝐶, the sediment concentration in the bed 𝐶_𝑏, the flow depth ℎ, the slope of the bed 𝜃, the angle of internal friction of the bed material 𝜙_𝑏 and the density difference ∆ between the sediment 𝜌 and the water 𝜌_𝑤.

Furthermore, FLATModel also incorporates basal erosion and does so from two approaches: the static approach considers the balance between the basal resistance forces 𝜏_𝑟𝑒𝑠, which can be represented by the Mohr‒Coulomb criterion for an infinite slope, and the friction force exerted on the bed 𝜏_𝑏. The dynamic approach consists of considering that the new material incorporated into the flow is accelerated to the average speed of the flow, whereby the incorporated mass depends on the availability of momentum or momentum. For the static approach, if these forces are not in equilibrium, that is, if the friction force is much greater than the resistance force, the model can calculate the basal erosion depth. ℎ_𝑒𝑏 is an expression that depends on the cohesion of the bed material 𝑐, the angle of internal friction of the bed material 𝜙_𝑏, the slope of the bed 𝜃, the flux density 𝜌, the depth of flow ℎ and the acceleration of gravity 𝑔 [30]. There are other approaches to simulate entrainment, including quasi mechanical models, empirical equations [94], three-dimensional models with SPH [95] and cellular automaton models [96].

Coupled models

As mentioned above, the modeling of mud and debris flows requires the use of different models that, when coupled, allow the magnitude of impact by an event to be determined. However, some software already includes several models coupled in one, which have been used to determine from the susceptibility of a basin to landslides and the propagation of flow with a single model on a medium scale or regional scale. Example models used include for the ASCHFLOW and Flow-R models.

ASCHFLOW is a two-dimensional, one-phase model that simulates material shedding, basal erosion, transport, and deposition of a mud and debris flow on a regional scale. This model simulates the transport of slipped material to secondary drains and to a main channel, as well as its deposition in a flood valley; all of this can be performed from a GIS interface. To determine the slipped areas, ASCHFLOW relies on the model of infinite slopes; it also includes a basal erosion rate calculated from empirical or semiempirical relationships according to the shear stress, slope and volume of the flow at each instant of time. As material is incorporated due to erosion of the channel, the model recalculates the depth of the flow. Finally, the transport and deposition of the flow is calculated from a Bingham or Voellmy rheological model solved by a finite difference scheme in a square mesh [97].

Flow-R is a spatially distributed empirical model in which the areas susceptible to landslides are determined from morphological criteria defined by the user, whereas the propagation of the flow is determined from friction information and algorithms that determine the direction of the flow will take [98]. This tool allows for determining susceptibility on a regional scale; however, its is not recommended to use it to carry out a detailed modeling of a channel and a particular event, since the propagation algorithm generates a wide range of possible event routes, , it does not consider the rheology of the flow.

This type of coupled model is very useful for determining the basins with the greatest threat within a region; with these models, not only is the intrinsic susceptibility of the basin taken into account, but a flow transit is also incorporated. However, these models were initial approximations; today, with a complete hydrological - geotechnical - hydraulic model, the detailed modeling of the events with a certain return period can be modeled in detail, and events that occurred in the past can be reconstructed. When this modeling is performed, large amounts of complexity of the relevant processes that occur during the generation, as well as the during the transit and deposition of a mud or debris flow, are incorporated. An example of this type of coupled model is OpenLISEM [42], which was initially developed as an erosion model; it evolved as a two-dimensional model that combines the model of infinite slopes with a hydrological model of events and a 2-phase flow transit model based on the equations of Pudasaini, 2012 [83].

Scope and limitations of mud and debris flow modeling

As previously described, the modeling of mud and debris flows involves interdisciplinary work to incorporate all the geological, geotechnical, hydrological and hydraulic aspects involved in the generation of a flow. Most investigations have focused on the simulation of flow propagation and the simulation of complex behavior of these flows, achieving great advances in the incorporation of rheological relationships. However, the majority of computational models include few rheology options, which This involves high uncertainty, considering that, depending on the selected rheological model, there may be substantial differences in the results obtained in terms of depth, speed and form of flow deposition. Likewise, hydraulic models are still very limited in their ability to describe the changes in the characteristics of a flow due to the incorporation of sediments throughout the transit of a flow, as well as the processes of lateral erosion.

Few studies consider conditioning factors other than landslides, such as volcanic eruptions, earthquakes and the failure of natural dams. Additionally, studies tend to simplify the modeling of landslides by incorporating the processes that exert the greatest influence and ignoring internal processes that can describe the slope's predisposition to detachment, as well as external factors such as water infiltration. In most cases, is simulated with static infiltration models that do not the antecedent conditions of rain, the evaporation of water from the soil and the transpiration of the vegetation that can be dynamic and that alter the conditions that pertain to the detachment of the slopes. Therefore, computational models alone are insufficient for performing a complete threat assessment; instead, they must be complemented with exhaustive studies in the field. There is no structured methodology for the coupling of hydrological–geotechnical–hydraulic models. This coupling would include the integration of the rain preceding the event, the detonation of the flow, the material that slipped from unstable areas toward the intermittent drains and main channels and the transit and deposition of this material in the flood valley. However, these processes are usually modeled as independent processes without a clear interlocking structure.

However, this type of analysis must advance from being purely investigative exercises to becoming true risk management tools, and guidelines regarding threat scenarios must be established to allow populations to incorporate prevention measures in their territories. Prevention measures would pertain to territorial planning and risk mitigation and would address the underlying factors that generate risks, such as deforestation. Finally, in cases where these measures are difficult to implement, mud and debris flow modeling should be the first step in generating early warning strategies that allow populations to live with risk, mitigate the impacts of mud and debris flows and increase the resilience of these populations.

Conclusions

In some countries, flows of mud and debris have generated great economic and life losses. Due to the effects of climate change, such as increases in the global temperature and changes in the duration and intensity of precipitation, in the future Alterations will be generated in the patterns of natural hazards, which will be increasingly frequent and of greater magnitude in the future. As a result, susceptibility to landslides will be increased among other conditioning factors of mud and debris flows associated with the climate, especially in high-slope basins with little vegetation or with high levels of interference from anthropic activities. Threat assessment supported by computational modeling is a tool that will allow populations at risk to implement prevention measures as well as early warning strategies for risk management and to minimize the impacts of these events.

However, modeling the flow of mud and debris is an arduous task that requires expertise in several disciplines, including geology, geomorphology, hydrology, hydraulics and informatics, because of the multiple conditioning factors and complex behavior associated with high concentrations of sediments. Hydrological modeling allows simulation of the infiltration, evaporation and runoff processes generated in the basin before and during events and is a fundamental input for landslide analysis; the latter can be carried out via statistical methods, heuristic methods, deterministic methods, or a combination of methods. Finally, once a hydrological model is integrated with the conditioning factors and the contribution of sediments to the flow, the flow must be transited through the main channel to the flood valley where the flow is deposited. There are multiple approaches for accomplishing this, such as the use of empirical relationships analyzed with two-dimensional multistage models and the use of regionally coupled models.

References

Notes

* Research article

Author notes

^a Corresponding author. E-mail: jessica.paez@unad.edu.co

Additional information

How to cite this article: J. P. Páez Pedraza, “Mud and Debris Flow Modeling as a Hazard Assessment Tool” Ing. Univ. vol. 29, 2025. https://doi.org/10.11144/Javeriana.iued29.mdfm

Descargar