# Discrete modelling of debris flows for evaluating impacts on structures

Rime Chehade, Bastien Chevalier, Fabian Dedecker, Pierre Breul, Jean-Claude Thouret

► To cite this version:

Rime Chehade, Bastien Chevalier, Fabian Dedecker, Pierre Breul, Jean-Claude Thouret. Discrete modelling of debris flows for evaluating impacts on structures. Bulletin of Engineering Geology and the Environment, 2021, 80 (8), pp.6629. 10.1007/s10064-021-02278-3. hal-03464203

HAL Id: hal-03464203

https://uca.hal.science/hal-03464203

Submitted on 3 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discrete modelling of debris flows for evaluating impacts on structures

Rime Chehade (1), Bastien Chevalier (1), Fabian Dedecker (2), Pierre Breul (1), Jean-Claude Thouret (3)

(1) Université Clermont Auvergne, CNRS, Clermont Auvergne INP, Institut Pascal, 63000 Clermont-Ferrand, France (rime.chehade@uca.fr, bastien.chevalier@uca.fr, pierre.breul@uca.fr)

(2) Itasca Consultants SAS, 29 Avenue Joannes Masset, 69009 Lyon, France (f.dedecker@itasca.fr)

(3) Université Clermont Auvergne, CNRS, OPGC, IRD, Laboratoire Magmas et Volcans, 63000 Clermont-Ferrand, France (j-claude.thouret@uca.fr)

Corresponding author: Bastien Chevalier bastien.chevalier@uca.fr

Published in Bulletin of Engineering Geology and the Environment

Received: 5 November 2020 / Accepted: 5 May 2021

https://doi.org/10.1007/s10064-021-02278-3

Bulletin of Engineering Geology and the Environment, 80(8), 6629-6645
**Abstract**

Lahars (volcanic debris flows) are natural phenomena that can generate severe damage and wreak havoc in densely populated urban areas. The evaluation of the forces and pressures generated by these mass flows on constructions (e.g., buildings, bridges and other infrastructure) is crucial for civil protection, assessment of physical vulnerability and risk management. The current tools developed to model the spread of flows at large scale in densely populated urban areas remain inaccurate in the evaluation of mechanical efforts. Here, we developed a discrete numerical model for evaluating debris-flow (DF) impact forces at the local scale of one structure (pillar or column) like a building, a bridge and other infrastructure. In this model, the large-sized solid particles that damage infrastructures and edifices are explicitly modelled using Distinct Element Method (DEM). We considered the fluid and fine-grained solid particles not only in the frame of the pressure exerted on structures, but also through their effects on the movement of particles, i.e. buoyancy and drag. The fluid velocity field and the fluid free surface obtained from Computational Fluid Dynamics (CFD) calculation based on Navier Stokes equations are imported in the DEM simulation. This model is able to reproduce a range of magnitudes of DFs in terms of volumes, velocities and flow heights. Finally, the model provides insights on impact forces generated by particles on structures and on hydrostatic and/or dynamic pressure due to the combined effect of fluid and solid phases. The model provides results consistent with existing empirical models.

**Keywords:** lahar, debris flow, Discrete Element Method, fluid, impact forces.

**1. Introduction**

Lahars (volcanic debris flows) are fast moving mixtures of sediment and water, other than normal streamflow, originating from volcanoes (Doyle et al. 2010; Manville et al. 2013). Their initiation results from a combination of geomorphological factors with rainfall/runoff, and they occur mainly on composite volcanoes (Lavigne and Thouret 2000). Lahar is an Indonesian term covering two types of flows (Manville et al. 2013):

(1) hyperconcentrated flows (HCF thereafter), which contain a sediment concentration of 20 % to 50 % by volume, are more turbulent and the deposits of which are less heterometric and better sorted; (2) debris flows (DFs thereafter) which include a sediment concentration greater than 50% by volume, the deposits of which deposits are massive, heterometric and very poorly sorted. According to Smith (1986) and Manville (2013), lahars change from HCF to DFs at solid concentrations exceeding 50% by volume. This paper focuses on the impact of DFs on structures.
Lahars can vary greatly in volume, peak discharge and runout. They present a complex flow behaviour with intricate interactions between fluid and particles, and in-between particles, depending on both fluid and particles properties (Manville et al. 2013; Iverson 2014). Lahar velocity (4 to 10 m/s on average) are moderate to high on the slopes of volcanoes, particularly following the occurrence of intense (10 to 40 mm/h) and/or long-lasting (60-90 minutes) rainfall events (Lavigne and Thouret 2000). Solid particles show a wide range of sizes from silt and sand up to boulders several meters across (Lavigne and Thouret 2000; Scott et al. 2005; Doyle et al. 2010, 2011; Manville et al. 2013).

Being natural and highly energetic flows with a high destructive power, lahars as well as non-volcanic DFs present a great danger and are difficult to control: e.g. Pinatubo in 1991 (Newhall and Punongbayang 1996), Mayon in the Philippines in 2006 (Paguican et al. 2009), Arequipa in Peru (Thouret et al. 2013; Pallares et al. 2015), Chaitén in Chile in 2008 (Major et al. 2016). They usually flow through river channels but they can spread out onto the valleys banks and fans, covering large areas with thick and coarse material, causing serious damage to cities, e.g., Sarno in Italy (Zanchetta et al. 2004), Arequipa in Peru (Thouret et al. 2014), El Porvenir and Rolando Rodriguez in Nicaragua (Scott et al. 2005).

DFs cause mainly three types of damage: direct impact, erosion and deposition (Hu et al. 2011). DFs generate three major mechanical effects (Zanchetta et al. 2004): (1) hydrostatic pressure depending on fluid height and density, (2) dynamic pressure depending on fluid velocity and density, and (3) collision of particles acting as missiles and depending on their mass and velocity (Thouret et al. 2020). These effects may vary in space and time, since debris flows usually break down into three segments, which are head, body and tail, each with different characteristics (Vallance 2000; Cui et al. 2015; Vallance and Iverson 2015).

The estimation of interactions between DF and construction or infrastructure is an essential step for assessing and mapping hazards and damage in densely populated areas at risk (Thouret et al. 2013, 2014; Ettinger et al. 2015). The dynamic interaction, or impact, between flow and structures (building or bridge) plays a key role in the evaluation of loss. Scientists and engineers have long faced the challenge of understanding and analyzing the DF impact forces, which can destroy buildings and infrastructure (Zanchetta et al. 2004; Federico and Amoruso 2008; Hu et al. 2011; Bugnion et al. 2012). Measuring the impact of DF on structures aims to use the quantitative data obtained from *in situ* measurements or laboratory experiments (flow velocity and height, impact pressure, and collision) to help engineers design structures that can mitigate the damage caused by DFs (Zanchetta et al. 2004; Federico and Amoruso 2008; Hu et al. 2011; Bugnion et al. 2012).
The evaluation of the forces generated by DFs on structures is challenging due to (1) the diversity and complexity of debris flows (fluctuations in sediment concentration, pore pressure and behaviour), (2) the topography of the channels and fans along which they propagate, and (3) the diversity of the types of structures and their exposure.

As shown in Table 1, many experimental studies have focused on quantifying the mechanical stresses induced on structures: *in situ* measurements (river channel), experiments in outdoor flumes, small-scale and miniaturized experiments (Zanchetta et al. 2004; Iverson et al. 2010; Bugnion et al. 2012).

*In situ* measurements evaluate precisely the forces induced by the events, but the flow conditions are difficult to control. Small-scale experiments solve part of this challenge but scaling down makes it difficult to satisfy similarity laws. In the case of *in situ* experiments, repeated iteration of tests becomes expensive and difficult, particularly with high volumetric flow rates, which can destroy the experimental set up and instruments (Iverson et al. 1992; Major et al. 1995).

Numerical modelling can help overcome the drawbacks of experimental measurements by allowing advanced analyses of the parameters characterizing a DF. The literature provides several numerical approaches able to evaluate the impact of DF on barrier or obstacles, ranging from purely granular to granular-fluid to purely fluid.

The following references illustrates this range of models. As for purely granular models, Teufelsbauer et al. (2011) used Distinct Element Method (DEM) to simulate dry granular avalanches down a ramp by introducing the rotation resistance that provides a realistic description of the flow behaviour. Cheng et al. (2018) used DEM (PFC2D) to simulate dry granular flows. Albaba et al. (2019) used also a purely granular DEM model with modified contact laws between particles to take into account the effect of fluid between particles through adhesion and cohesion. As for fluid-granular models, Trujillo-Vela et al. (2020) performed a complete SPH-DEM coupling to simulate DFs over a distance of 300 m from a slope and modelling the runout. Leonardi et al. (2016) developed a DEM-LBM (Lattice Boltzmann Method) model to simulate the fluid-grain mixture, but this approach has a very high computation cost. Kattel et al. (2018) simulated two phase DFs based on the depth-averaged models to measure the impact on tetrahedral obstacles with different sizes, location and orientation with respect to the flow.

Zhao et al. (2014) presented a DEM-CFD coupling system for the study of grain sedimentation in fluids. Li and Zhao (2018) considered a DEM-CFD coupling to simulate the impact of DFs on a flexible barrier: the coupling takes into consideration the effect of fluid on particles and the interaction between DF particles and the barrier, which is also modelled by DEM.

With respect to purely fluid models, Federico et al. (2008) used the Finite Element Analysis to analyse the impact of fluid against structures. Dai et al. (2017) carried out a fluid-structure coupling based on a purely fluid model
using Smoothed Particle Hydrodynamics SPH and measured the impact of DFs on a structure (check dams). Mead et al. (2017) used the numerical SPH model to predict the motion of particles within the fluid, for three types of flows (non-Newtonian, HCF and DFs): each of them play a different role in determining the hydrostatic and dynamic pressures applied to buildings.

The model proposed in this paper to estimate the impact stresses of DF on obstacle is a granular-fluid model based on DEM and using separate CFD calculation results. It is an intermediate between “purely” granular and “purely” fluid. This model is able to estimate the impact of a given DF in particular the contribution of the solid part at the scale of the structure (a pillar), without reproducing the complexity of the flow mechanisms on the long distance. Hence, the model is achieved within a reasonable calculation time. Fluid phase affects the particles of the DEM model through drag force and buoyancy. Fluid calculation is computed once and before the DEM simulation is conducted. The fluid forces are added continuously during the flowing process of the particles and the interaction is only considered in one way: action on the fluid on particles.

Input parameters are the DF characteristics and the obstacle shape, size and orientation. These impact stresses obtained represent input data for (civil) engineering models to assess the mechanical behaviour of constructions or structural elements.

Firstly, we introduce the numerical approach used to simulate DFs. Secondly, we present the calibration process of the model parameters to reproduce a flow of given characteristics. Thus, we measure the interaction of the flow upon impact with the pillar. Finally, we display the measurements of efforts and stress level obtained from the numerical model.

## 2. Methods

### 2.1. Particles and fluid modelling

The study focuses on modelling the particles impacts generated by water-rich flows on obstacles located across an artificial channel. We modelled the coarsest part of the solid fraction of the flow explicitly with Distinct Element Method (DEM): we referred to this fraction as particles. The DEM is useful to describe the mechanical behaviour of discrete materials, usually soil or rocks, as it considers the material as a collection of rigid particles (Cundall and Strack 1979). We use the Itasca software Particle Flow Code in 3 Dimensions or *PFC3D* (Itasca 2016).

As for the fine-grained solid particles, composed of silt and fine sand, we did not explicitly modelled them but integrated them in the fluid part. The fluid phase means the mixture of water with the fine particles that are
transported. The movement of these particles can be assimilated to a fluid flow. The model will separately consider
the movements of particles and the fluid and will model the effect of the fluid on the particles.

The fluid phase plays a major role in the movement of the flow of particles: the effect of the fluid phase on the
particles had to be considered. Furtney et al. (2013) discussed the methods to model fluid-particle systems with
DEM by using a coupling between fluid model and PFC3D, depending on the physics of the problem, for example,
whether or not the fluid-particle interactions at the local scale prevails in the performance of the system.
Subsection 3.2 describes the contact laws used in the model and the associated parameters. The model must be
able to reproduce a steady flow, with given characteristics, in particular velocity and discharge. The velocity
vectors of the fluid phase have to be determined to compute the interactions produced by the fluid phase on
particles. We obtained a velocity field of fluid phase by simulating the fluid flow with 3D Navier-Stokes equations
with free surface changes using Telemac3D. This is a finite element open source code for free-surface flow
developed by the LNHE (Laboratoire National pour l'Hydraulique et l'Environment) and EDF (Electricité de
France) and based on 3D hydrodynamic equations, to compute the velocity field and the free surface of the fluid
phase (Moulinec et al. 2011; Rameshwaran et al. 2013). Then, we used the obtained velocity vectors to compute
the fluid effect on particles. The constant values are the velocity field of the fluid and flow height. The fluid effect
on a particle (i.e., drag force and buoyancy force) is updated at each time step, depending on its position in the
channel and its velocity, on the one hand, and on the fluid velocity and flow height at that precise location, on the
other hand. The flow chart represented in Figure 1 describes each of the steps of the simulation process.

## 2.2. Fluid effects on particles

In our model, the mixture of fine particles and the fluid induces mechanical effects on the particles through
buoyancy and drag. All other interactions are neglected (Zhao et al. 2014, 2018). The interactions from the
particles to the fluid are not considered, and the model can therefore be assimilated to a one-way coupling: we
model the flow in a reduced distance range around the structure (i.e. a Reinforced Concrete RC pillar), and we
focus on the impact pressure on this structure. At this scale, we considered that the main particle-fluid interactions
are the drag force and the buoyancy: those interactions will influence greatly how particles are driven through the
channel and so how they will collide with the obstacle. Similar assumptions are made in other CFD-DEM coupling
models in the literature, e.g., (Zhao et al. 2018 and Li and Zhao 2018). Buoyancy or Archimedes force generates
an upward, dispersive force exerted on an immersed particle. The buoyancy force is written as follows:
$$F_b = -\rho_f v_i g z \qquad (\text{Eq.1})$$

where $\rho_f$ is the density of the fluid phase, $v_i$ the volume of the immersed particle, $g$ the gravitational acceleration, and $z$ the vertical vector unit pointing upward.

Drag is a force acting by a fluid to a solid resulting from the difference of velocity between them. The common form of the drag force $F_d$ exerted on a single spherical particle can be written as follows Pudasaini (2012) :

$$F_d = \frac{1}{2} C_d \rho_f \frac{\pi d^2}{4} \|V_f - V_b\| (V_f - V_b) \qquad (\text{Eq.2})$$

where $\rho_f$ is the density of the fluid phase (kg/m³), $d$ is the particle diameter (m), $C_d$ is the drag coefficient, $V_f$ is the velocity vector of the fluid and $V_b$ the velocity of the immersed particle (m/s). The drag coefficient $C_d$ can be written as follows Zhao et al. (2014) :

$$C_d = \frac{24}{Rep} (1 + 0.150Re_p^{0.681}) + \frac{0.407}{1 + \frac{8710}{Rep}} \qquad (\text{Eq.3})$$

where $Re_p$ is the particle Reynolds number, expressed as $Re_p = \frac{\|V_f - V_b\|d}{\eta_f}$ with $\eta_f$ is the fluid kinematic viscosity.

As proposed in Zhao et al. (2014), a correction factor $n^{-\xi+1}$ depending on the solid fraction in the flow ($n$) can be applied to the drag force, given in Eq.2, to consider the influence of the particle concentration and the contacts between particles. The new drag force $F_{dc}$ taking into account this correction is expressed as Zhao et al. (2014) follows:

$$F_{dc} = \frac{1}{2} C_d \rho_f \frac{\pi d^2}{4} \|V_f - V_b\| (V_f - V_b) n^{-\xi+1} \qquad (\text{Eq.4})$$

where $\xi$ is a term ranging between 3.4 and 3.7, expressed as:

$$\xi = 3.7 - 0.65 \exp \left[ \frac{-(1.5 - \log Re_p)^2}{2} \right] \qquad (\text{Eq.5})$$

## 2.3. Contact laws

The linear contact model is used for two types of particles interactions that we encounter in the model: contact between particles belonging to the solid phase, and contact between particles and wall, which represents the collision between particles and a rigid structure. This contact model uses contact stiffness in the normal and shear directions: $k_n$ and $k_s$, respectively. A Coulomb limit boundary applies to the shear force using a friction coefficient $\mu$. We also considered a viscous damping in the normal direction at the contact: the ratio of viscous damping to critical viscous damping is set to 0.4, based on work carried out by Albaba et al. (2015). This damping dissipates a part of energy during a shock at a contact.
The particles are spherical, so we added a rolling resistance in the contact model to account for the effect of the non-spherical particle shape. The rolling resistance mimics the behaviour of particles with angular shapes (Itasca 2016). With rolling resistance, internal moment of particles is linearly incremented with the particle rotation, up to a threshold value. The increment of internal moment is given by:

$$ \Delta M_r = -k_r \Delta \theta_b \qquad (\text{Eq.6}) $$

where $\Delta \theta_b$ is the relative bend-rotation increment and $k_r$ the rolling resistance stiffness. The rolling resistance stiffness is calculated as follows:

$$ k_r = k_s R^2 \qquad (\text{Eq.7}) $$

where $k_s$ is the contact stiffness in the tangential direction and $R$ the contact effective radius of the two elements in contact:

$$ \frac{1}{R} = \frac{1}{R_1} + \frac{1}{R_2} \qquad (\text{Eq.8}) $$

The magnitude of the updated rolling resistance moment is then checked against the threshold limit $M^*$ given by:

$$ M^* = \mu_{rr} R F_n \qquad (\text{Eq.9}) $$

where $\mu_{rr}$ is the rolling resistance friction coefficient and $F_n$ the normal linear contact force. The values of contact model parameters are discussed in section 3.2.

## 2.4. Model procedure

Firstly, we generated a representative elementary volume (REV) of particles. This REV exhibits the packing characteristics that particles are assumed to adopt in the flow. Secondly, the REV is positioned in a supply chamber, whose width is the one of the flow paths at the entrance of the channel (Fig.2). The particles located in the supply chamber all move toward the entrance of the channel with a fixed velocity in the direction of the flow path axis and equal to the fluid phase velocity. As soon as the particles are in the channel, the fluid effects applied to the particles according to the velocity vectors corresponding to the position of each of the particles in the channel. Particles can move along the channel according to the fluid forces applied on them and their interactions with other particles and walls. The REV periodically repeated in the supply chamber, so that a steady particles flow feeds the channel.

## 3. Model parameters

Table 2 gives the characteristics of the simulated debris flow. We fixed the flow velocity $V_{DF}$ to 3 m/s. This corresponds to the lower value of velocity range for such flows, with an estimated Froude number of 0.78 (Hübl
et al. 2009). However, the parameter set considered in this work represents the properties of a debris flow, but the procedures are flexible enough to allow the application of another parameter set, as we will see later.

**3.1. Geometry of the simulated channel**

We calibrated the numerical parameters by modelling a channel with simple geometry: the channel is straight with a rectangular cross section and a low-gradient (2%) slope in the direction of the flow path axis, which is similar to the slope measured in some sections of the Dahlia ravine channel in the city of Arequipa (Peru) (Mead et al. 2017). The width of the channel is 10 m. We have chosen a flow channel 25 m in length, so that the flow can settle in the channel, and at the same time short enough to limit the number of particles simultaneously present during the simulation run. The supply chamber is located at the top entrance of the channel, and once the particles exit the channel, we delete them from the model (Fig.2).

**3.2. Physical and mechanical parameters of the model**

As mentioned in subsection 2.3, a REV of spherical particles is prepared in advance of the simulation run. We chose the REV volume to supply the channel with particles and to obtain a discharge between 35 to 55 m³/s in a 10 m wide flow channel. The REV height is set to same height as the flow, which is 1.5 m.

The grain-size distribution of particles is uniformly distributed between $d_{min}$ and $d_{max}$. The threshold diameter $d_{min}$ separating particles from fine fraction corresponds to the minimum particle size explicitly modelled in the simulation. In natural flow channels, the size range of coarse particles is much wider, from silt and sand up to boulders (Zanchetta et al. 2004; Pallares et al. 2015). We fixed the particle size range based on the following conditions: (1) to represent part of the range observed in the field, (2) to simulate the particles which significantly contribute to the impact force on structures, (3) to set a minimal diameter large enough to constrain the total number of particles in the simulation and to reduce calculation time, (4) to set a maximal particle size small enough to obtain a representative number of particles between the channel bottom and the top flow surface. Based on these conditions, the minimal diameter was set to $d_{min}= 0.10$ and the maximal diameter $d_{max}= 0.40$ m.

Within the REV, particles are generated at an initial volume fraction of 0.50. The REV is considered representative if its length is large enough. A set of particles is prepared to obtain a representative number of particles of 9,500 spherical particles, for a REV length of 5 m.

Regarding the mechanical parameters of the contact model, the contact stiffness is fixed to $k_n= 10^7$ N/m in the normal direction and $k_s = 5 \times 10^6$ N/m in the tangential direction with $k_s/k_n = 0.5$ (Itasca 2016). The value of contact
stiffness has to be large enough to reduce the overlap between particles while maintaining numerical time step long enough.

Besides the contact stiffness, we studied the effect of other mechanical properties to reproduce the desired flow conditions. At particle-particle contacts, different friction coefficients from 0.05 to 0.4 were tested. The Coulomb friction coefficient between particles and walls varied between 0 and 1. The density of particles ranges between 2500 and 2700 kg/m³, a range that encompasses the average density of rock (Iverson et al. 1997). For rolling resistance in particle-particle contacts, different values of rolling resistance friction coefficient $\mu_{rr}$ between 0 and 0.6 were tested (Itasca 2016).

The density of fluid phase, composed of water and fine solid particles, ranged between 1000 and 2000 kg/m³. Values between 0.03 and 0.075 of the dynamic viscosity used to compute the drag forces (in Eq.4) were tested, based on results from Mead et al. (2017) and Cui et al. (2015). Table 3 compiles the results of this parametric study.

### 3.3. Model calibration

The simulation of the flow lasts 12 seconds, which is sufficient to fill completely the channel and to obtain a continuous and steady flow. At the end of the simulation, the channel contains about 52.655 particles. During the simulation, we recorded the height of the particle layer, the average particle velocity and the total discharge (both particle and fluid phase) across two cross sections arbitrarily chosen at 8.0 m and 15 m from the entrance of the ravine. We calculated the average bulk density of the flow in the channel at the end of the simulation. The solid volume fraction $v_s$ is defined by Zanchetta et al. (2004) as:

$$v_s = \frac{\rho_{\text{bulk}} - \rho_f}{\rho_s - \rho_f} \qquad \text{Eq.10}$$

where $\rho_{\text{bulk}}$ is the apparent density of the flow, $\rho_s$ is the particle density and $\rho_f$ the fluid density.

The model is calibrated to set all properties and obtain a flow presenting typical characteristics of DF, in terms of solid concentration, density and flow rate, while the model sensitivity to every parameter is studied. The typical DF characteristics are based on Iverson et al. (1997), Hu et al. (2011) and Mead et al. (2017).

Table 3 presents 12 runs as well as the results obtained in terms of flow characteristics and variability with respect to a reference case. In sum, the model parameters that strongly influence the results can be ranked as follows (Tables 3, 4): (1) the rolling resistance influences the discharge and solid volume fraction, as it helps dilatancy to appear in the granular contact network, (2) fluid and block density have a major influence on the bulk density, (3)
the dynamic viscosity has a limited effect on velocity, (4) the friction coefficient does not influence the results obtained in this study.

Based on the parametric study carried out (runs 1 to #10) and the parameters influence described above (Tables 3, 4), a strategy to calibrate the model was developed: we started by setting the rolling resistance parameters which influences the characteristics of the obtained flows, particularly the volume fraction and discharge. Taking into account that blocks are modelled as spherical particles, the rolling resistance allows angularity of particles to be mimicked by limiting the rolling of spherical particles. This parameter was set to 0.2 to adjust the solid concentration and the discharge. In a second step, we set the particle density and the fluid density (a mixture of fluid and fine particles). Both parameters are directly related to the apparent bulk density and solid fraction. Based on the run #2, which yields the lowest bulk density, we set a particle density to 2500 kg/m³ to meet an apparent bulk density of the flow around 1900 kg/m³. We were then able to set the fluid density to a value of 1100 kg/m³ to reach the targeted flow properties: the apparent bulk density and, to a lesser extent, the average velocity. The calibration of the densities leads to a solid fraction around 55%, which corresponds to the considered flow. The dynamic viscosity was fixed to 0.048 Pa. s (e.g., Mead et al. 2017) to finally adjust the flow velocity. Lastly, we have considered a particle-particle friction coefficient of 0.4, neither too weak nor too strong, which does not influence our results. The flow modelled with the final set of parameters (run #11) shows the following characteristics: a solid fraction of 54.7 % in volume, a bulk density of 1867 kg/m³, and a volume discharge of about 40 m³/s (Fig.3).

# 4. Numerical set up for measuring impact forces

## 4.1. Geometry of the simulated channel

After calibrating the flow parameters in a channel without obstacle, we used the model to evaluate the forces induced by this particular flow on a pillar located in the centre of the channel at a distance of 15 m from the entrance (Fig. 4). The obstacle is a pillar with a square shaped horizontal cross section (1 m × 1 m). *FLAC3D* is used to model the pillar based on the finite difference method. The pillar is meshed with brick-shaped zones with 8 vertices. Mohr-coulomb elasto-plasticity model is assigned to all zones, with some fixed properties (density, elastic modulus, Poisson's ratio) relevant for a RC pillar. Proske et al. (2011) shows that the accumulation height of DFs upstream of the obstacle is higher with an obstacle located at the centre of the channel (Fig. 5). If $h_{Mu,2}$ corresponds to this measurement, we can write $h_{Mu,2} = h_{ac} + h_{Mu}$, with $h_{Mu}$, the height of DFs at the entrance
of the channel and $h_{ac}$, the accumulated height of DFs. We used $h_{Mu,2}$ to calculate the impact pressure of DFs on the obstacle.

## 4.2. Fluid phase calculation

In the case of the channel including a pillar, this fluid flow is not uniform anymore; thus, we have to compute the velocity field and free surface of the fluid phase velocity along the channel. We fix the flow rate and the fluid height at a 15 m distance ahead from the pillar. We fixed a constant flow rate $Q$ at the entrance of the channel and a constant flow depth $H$ at the exit of the channel. The density of the fluid is set to the same value as in the DEM simulation, i.e. 1100 kg/m$^3$. As a result, we obtained the velocity of the fluid and the free surface at each point of the triangular mesh used for the fluid calculation and particularly all around the pillar. Finally, we export these data into PFC3D and used them to compute drag forces as described in Equation 4. Figure 6 shows the free surface (a) and the velocity field (b) in the part of the channel (12m in front of the pillar but only 1.5m in the back). We conducted the DEM simulation with fluid interactions over a duration of 12 s.

## 4.3. Parametric study

We conducted a parametric study of the impact force induced by DFs on pillar by changing the flow parameters $h_{Mu}$ and $V$, which influence the Froude number (Table 5). Table 5 also gives the time-averaged values (between 9 and 12 s) of the maximal pressure and the mean pressure on the pillar, the standard deviations for each and the maximal DFs height $h_{Mu,2}$ (at the contact with the obstacle) in each case.

### 4.3.1. Reference case

Our reference flow includes a flow density of 1867 kg/m$^3$, a velocity $V_{DF}$ of 3 m/s, and flow height of 1.5 m corresponding to a Froude number $Fr = 0.78$. We previously computed the velocity field and free surface positions with the fluid code. Then we modelled the DF in the channel with PFC3D. The computing time for the reference case lasted about 8 hours. The simulations were performed on computer with a processor Intel (R) Core (TM)i9-9900 CPU and 32Go RAM (DDR4, 2666Mhz). During the simulation, we record the distribution of the normal stresses of particles on the pillar. Figure 7 shows the distribution of the pressure applied by particles on the obstacle, which is calculated directly in the DEM model from the spatial distribution of the contact forces between the particles and the pillar, as a function of time. The red curve is a moving average window with a period $k = 100$. The red curve is used to analyze the data points of the black curve, by averaging different subsets of the full dataset. It takes 5 s for the first particle to impact the pillar. From $t = 5$ s to $t = 9$ s, the stress progressively
increased. After t = 9 s, the impact stress curve, based on average moving window (k=100), seems to stabilize at around 74 kPa. Figure 8 shows a longitudinal cross section of the model flow for different durations between 4 s and 12 s after the simulation starts. At t= 4 s, the impact did not occur yet. At t= 6 s, the impact with the pillar occurred and the accumulation of particles ahead the pillar starts. At t= 6 s, we observe the progressive rise of the blocks upstream the pillar. At t= 9 s, the accumulation of particle upstream the pillar reaches a maximum corresponding to the onset of a plateau-shaped impact stress on the pillar, and a steady regime starts. At the end of the simulation (t = 12 s), the regime is steady compared to state at t= 9 s.

To evaluate the total impact pressure induced by debris flows on the pillar, we account for the contribution of the fluid phase that includes the water and fine-grained particles mixture. The pressure due to the fluid itself is hydrostatic and hydrodynamic (Hübl et al. 2017). To calculate the hydrostatic pressure $p_{static}$ of the fluid, we used the following formula (Hübl et al. 2009, 2017):

$$p_{\text{static}} = \rho_f \cdot g \cdot h_{Mu} \qquad \text{Eq.11}$$

where $\rho_f$ is the fluid density, g is the gravity and $h_{Mu}$ is the debris flow height.

The hydrodynamic component $p_{dynamic}$ due to the fluid pressure is calculated as follows:

$$p_{\text{dynamic}} = \rho_f V_0^2 \qquad \text{Eq.12}$$

where $V_0$ is the velocity of the flow at the vicinity of the upstream face of the pillar. This velocity value is close to 0. So we discarded the dynamic component of the fluid pressure.

The total impact pressure due to the fluid phase of the flow $p_{fluid}$ is expressed by:

$$p_{\text{fluid}} = p_{\text{static}} + p_{\text{dynamic}} \qquad \text{Eq.13}$$

Figure 9 shows three results for the reference case:

(i) the evolution of the average stress applied by the particles on the pillar (from t= 9 sec to t= 12 sec);

(ii) the static fluid stress profile;

(iii) the total pressure as the sum of the two previous curves, all of them as a function of the distance to the free surface.

Figure 9 also shows the maximum stage of the particles (accumulation height of particles) and the maximum stage of the fluid phase upstream of the pillar. In this case, the total pressure of DF increases towards the bottom of the pillar. In the lower part of the pillar (below 1.4 m), the pressure increases less rapidly with the depth in the flow.

In addition, there are strong timewise fluctuations of the stresses on the pillar: we link this variability of the pressure to the variability in time of the shocks between particles and pillar due to particle size and particle velocity.
### 4.3.2. Effect of flow height

Figure 10 shows the distribution of normalised total stress (pressure of particles + fluid pressure) as a function of the normalised height from the bottom of the channel $z/h_{Mu}$ for the three cases of flow heights. We calculated an average stress over a time period between 9 and 12 s. Each curve obtained represents four different calculations of the same flow conditions each using different randomly generated particle REVs. The total stress is normalised against $\rho_{Mu} V_{DF}^2$, with $\rho_{Mu} = 1867 \text{ kg.m}^{-3}$. We note $\hat{p} = \frac{\rho_{MAX}}{\rho_{Mu}V_{DF}^2}$ thereafter.

The total stress increases as the flow height increases and also closer to the bottom of the pillar. Maximal stress factor is observed in the lower 10% of normalised flow height on the pillar: for $h_{Mu} = 4$ m, $\hat{p} = 21$, for $h_{Mu} = 3$ m, $\hat{p} = 16$ and for $h_{Mu} = 1.5$ m, $\hat{p} = 7$. Froude number vary from 0.49 to 0.78 corresponding to the variation of flow height. The maximal normalised pressure, or maximal dynamic pressure coefficient $\hat{p}$ increases as Froude number decreases in the range of Froude numbers considered here, which is consistent with Hübl et al. (2009), Scheidl et al. (2013), Cui et al. (2015) and Wang et al. (2018): a negative correlation is obtained for the range of Fr 0.49-0.78 corresponding to the flow height variation in our model (Fig. 12).

A great drop in reduced pressure $\hat{p}$ can be observed for the range of Fr 0.49-0.78. The graph (Fig. 10) also indicates the maximal DFs height against the pillar $h_{Mu,2}$ (including accumulation), which explains why ratios $z/h_{Mu}$ are greater than 1.

### 4.3.3. Effect of flow velocity

Figure 11 shows the distribution of the normalised total stress (particles + fluid) on the pillar for the different flow velocities $V$ and for a flow height of 1.5 m at the entrance of the channel, as a function of the normalised flow height on the pillar. We calculated again an average stress over a period between 9 and 12 s, and each curve represents four different calculations of the same flow conditions each using different, randomly generated particle REVs. Figure 11 also shows the maximal DFs height against the pillar $h_{Mu,2}$ (including accumulation), which are very close in these three cases.

The higher the fluid velocity is, the greater the turbulence around the pillar will be: the particles are more agitated in the upper part of the flow. Figure 11 shows the normalised pressure versus normalised flow height on the pillar for different flow velocity $V_{DF}$ and for the same flow height $h_{Mu}$: the normalised stress curves are similar. Froude number vary from 0.78 to 1.56 corresponding to the effect of flow velocity. The distribution of normalised dynamic pressure coefficient on the pillar changes slightly as the Froude number increases, however, the maximal normalised pressure coefficients are very close in the three cases: for $V_{DF} = 3$ m/s, $\hat{p} = 7$, for $V_{DF} = 4.5$ m/s, $\hat{p} = 6.8$ and for $V_{DF} = 6$ m/s, $\hat{p} = 5.7$. These compare fairly well with the results obtained by Hübl et al. (2009),
Scheidl et al. (2013), Cui et al. (2015) and Wang et al. (2018): the reduced pressure coefficient $\hat{p}$ does not evolve largely whereas $Fr$ varies twofold (from 0.78 to 1.56), and a negative correlation is also obtained meanwhile (Fig.12).

# 5. Discussion

As mentioned in the introduction, the total impact pressure generated by DFs on a pillar entails two forces: impact of large solid particles and fluid pressure through its hydrostatic and dynamic components (Hübl et al. 2017).

## 5.1. Comparison with empirical models

The literature does not provide data from laboratory or in situ (channel or flume) experiments that would help us to validate the measured impact forces on obstacles based on our scenarios. Many authors (e.g., Liu et al. 2018; He et al. 2018; Zhao et al. 2020) chose to compare the results of their models by considering the rotation and deformation of the obstacle due to flow impact, after calibrating the flow, but without considering the presence of the blocks. Our approach is different, since the flow has been calibrated and considered the fluid and the presence of blocks at the same time, while our objective was to validate the impact forces measured via numerical simulations. Therefore, we have compared the obtained impact pressures with existing empirical models, which stem from experimental observations.

A range of methods exist to predict the impact pressure of debris flows on pillars (Zanchetta et al. 2004; Bugnion et al. 2012; Hu et al. 2011, to name a few). Most in situ experiments published to date have measured the density, velocity, and the height of natural debris flows but impact pressures are often considered proportional to the hydrodynamic pressure of a fluid or to the sum of hydrostatic and hydrodynamic pressures (references in Hu et al. 2011 and Thouret et al. 2020). The proportionality coefficients of the different methods encompass different characteristics of debris flows: e.g. hydraulic characteristics such as viscosity and turbulence.

For example, Hu et al. (2011) have proposed a method based on a field measurement dataset acquired in Jiangjia Ravine (China) in which the impact mean pressure ($p$) on buildings was estimated using hydraulic models and the following formula:

$$p = c_m (0.5 \rho g h_{Mu,2} + \rho V_{DF}^2) \qquad \text{Eq.14}$$

where $\rho$ is the debris flow density, $h_{Mu,2}$ the maximal debris flow height against the pillar (including accumulation), $V_{DF}$ the flow velocity and $c_m$ a dimensionless empirical coefficient including both hydrostatic and hydrodynamic pressures. Hu et al. (2011) found that the coefficient $c_m$ ranged between 0.28 and 9.88 based on
the Jiangjia Ravine dataset in China. Such a large ranging coefficient $c_m$ is due to the sensitivity of the particles' effect for such impact. To obtain a total impact pressure as high as 74 kPa, corresponding to our case study, $c_m$ should be equal to 1.91 which is in the range of values defined by the authors

Zhang et al. (1993) computed the total impact pressure, based on field measurements carried out in the Jiangjia River in China. The following expression was obtained:

$$p = c_d \rho V^2 \cos^2 \alpha \qquad \text{Eq.15}$$

where $\alpha$ is the angle between the flow direction and the direction of the normal to the impacting plane, and $c_d$ is an empirical factor usually between 3 and 5. In our case study, we should consider a factor $c_d$ of 4.4 to obtain a 74 kPa total impact pressure with $\alpha = 0$. This value is included in the range given by Zhang et al. (1993) based on in situ measurements.

Thus, based on both abovementioned examples, the results of total impact pressure using our model are in the range of the results derived from empirical model predictions.

In addition, we compared the maximal pressure on the pillar $p_{max}$ with the prediction of analytical model given by Hübl et al. (2009). It is a representation of the following hydro-dynamic formula:

$$p_{max} = 5\rho_{mu} \cdot V_{DF}^{0.8} (gh_{Mu,2})^{0.6} \qquad \text{Eq.16}$$

where $\rho_{mu}$ is the debris flow density, $V_{DF}$ is the velocity of debris flow, $g$ is the gravity and $h_{Mu,2}$ is the maximal debris flow height (including accumulation) whereas $h_{Mu}$ is the flow height before impact or without obstacle.

In the reference case, we obtained a maximal total pressure of the DF on the pillar of $p_{max}$ (numerical) = 119 kPa (Fig.9). The model based on hydro-dynamic formula yields $p_{max}$ = 150 kPa, with, $h_{Mu,2}$ = 2.4 m and $V_{DF}$ = 3 m/s.

The maximal pressure $p_{max}$ obtained from Eq.16 as a function of corresponding $h_{Mu,2}$ (Table 5) is calculated for range of three different velocities: for $V_{DF}$ = 3 m/s, $p_{max}$ = 150 kPa, for $V_{DF}$ = 4.5 m/s, $p_{max}$ = 217 kPa and for $V_{DF}$ = 6.0 m/s, $p_{max}$ = 310 kPa. These maximal pressures are in accordance with the values provided by numerical model (Table 5).

## 5.2. Limitations and outlooks

Our numerical model simulates DFs with a high concentration in particles whose characteristics can be adjusted according to the expected flow: velocity, discharge, bulk density and flow height. This model enables to assess the stresses (particles and fluid) resulting from flow interactions with pillars. In this way, the efforts generated by
debris flows can be recorded, and then used to quantify the vulnerability of structures and help engineers design structures to resist the effects of water-rich mass flows.

The CFD approach of Telemac3D used to compute the fluid velocity field and free surface implies a Newtonian fluid rheology. DFs are complex, non-newtonian, viscoplastic, two-phase flows, but in our case study, we model the fluid phase only to obtain a velocity field necessary to drive particles along the channel and to approximate the free surface.

The evaluation of a debris-flows impact pressure results from the analysis of the impacts exerted by both the fluid part and the solid part, separately. We still need to identify the effect of other parameters on the model: e.g., particle size and shape, dynamic viscosity, shape of particles, and orientation of pillars (Mead et al. 2017; Kattel et al. 2018). In this study, we considered a maximum diameter of 0.4 m. In future work, the sensitivity to the particle size distribution will be studied as well as the particle shape. These parameters can be easily modified in the numerical model. In addition, our model does not consider the deformability of structural pillars yet.

Many codes and methods have been presented in the literature to simulate the impact of a DF on a structure (Sosio et al. 2007; Toyos et al. 2008; He et al. 2018; Zhao et al. 2020). In most cases, modeling DF processes requires many assumptions and simplifications, making the application of the model more or less deviate from reality (Sosio et al. 2007; Toyos et al. 2008).

## 6. Conclusion

We developed a numerical model for debris flows, based on an explicit modelling of particles using the Distinct Element Method and considering the fluid effects. This model is a simplified description of the debris-flow phenomena that focuses on their overall mechanical effects on buildings or other obstacles along their path. The model is not meant to reproduce the initiation or the evolution of DFs, but focuses on the determination of forces exerted on structures at the scale of the construction element. This study has focused on reproducing the movement of coarse solid material associated with debris flows, based on simple and measured flow characteristics: velocity, discharge, size of the largest solid particles. We considered the large-size solid part included in the flow mixture as particles, while we integrated the fine-grained solid fraction in the fluid phase. The fluid phase generates two effects on particles, namely buoyancy and drag forces. We calibrated the model to obtain given flows characteristics, i.e. the flow velocity, discharge, and density.

To summarise, we highlight the following results:
1. The numerical model considers two flow phases: particles modelled by DEM and the fluid phase to estimate the velocity field and the height of the fluid.

2. The model is able to reproduce a debris flow with specific input parameters: bulk density, solid fraction average velocity, and height.

3. The model is able to separate on one hand the impact pressure induced by particles on pillars and, on the other hand, the impact pressure from the fluid phase. The model is also able to measure the variation of impacts as a function of time and to obtain the distribution of impact forces on the structure.

4. The influence of DF height and velocity of impact pressure obtained with the model compare fairly well with the results derived from prediction models.

5. The impact force distribution on obstacles and structures can be used as input data in civil engineering code for designing purposes or assessing vulnerability of constructions.

In further research, we have to use the model with a larger range of debris flows (discharge, density, particle size distribution, velocity, and flow height).

The advantages of the proposed method are typical of simplified approaches to model complex natural phenomena. We did not implement a two-way coupling model (Trujillo-Vela et al. 2020) because our key objective was to introduce a practical solution to evaluate the impact pressure of DFs, focusing on the effects of the blocks. Using the proposed model, we attempt to estimate block pressure in addition to the fluid contribution to obtain the total impact pressure. Therefore, the proposed numerical model complements rather than competes with more complex physically based approaches.

## Acknowledgements

The authors thank the Municipality of Arequipa, the Arequipa Civil Defence as well as the department of geology of the University Nacional San Agustin in Arequipa. This study benefits from the financial support by the “Pack ambition” research project 2018 of the Region AURA Auvergne Rhône-Alpes and from the Clervolc LabEx.

## References

- Albaba A, Lambert S, Nicot F, Chareyre B (2015) Relation between microstructure and loading applied by granular flow to a rigid wall using DEM modeling. Granular Matter 17:603-616.

https://doi.org/10.1007/s10035-015-0579-8
Albaba A, Schwarz M, Wendeler C, Loup B and Dorren L (2019) Numerical modeling using an elastoplastic-adhesive discrete element code for simulating hillslope debris flows and calibration against field experiments. *Natural Hazards and Earth System Sciences* 19(11):2339–2358.
https://doi.org/10.5194/nhess-19-2339-2019

Bugnion L, McArdell BW, Bartelt P, Wendeler C (2012) Measurements of hillslope debris flow impact pressure on pillars. *Landslides* 9(2):179-187. https://doi.org/10.1007/s10346-011-0294-4

Cheng YM, Li N, Fung IWH (2018) Laboratory and Field Test and Distinct Element Analysis of Debris Flow. *Natural Hazards and Earth System Sciences Discussions* 1–32. https://doi.org/10.5194/nhess-2018-36

Cui P, Zeng C, Lei Y (2015) Experimental analysis on the impact force of viscous debris flow. *Earth Surface Processes and Landforms* 40(12):1644-1655. https://doi.org/10.1002/esp.3744

Cundall PA, Strack ODL (1979) A Discrete Numerical Model for Granular Assemblies. *Géotechnique* 29(1):47-65. https://doi.org/10.1680/geot.1979.29.1.47

Dai Z, Huang Y, Cheng H, Xu Q (2017) SPH Model for Fluid-Structure Interaction and Its Application to Debris Flow Impact Estimation. *Landslides* 14(3):917–928. https://doi.org/10.1007/s10346-016-0777-4.

Doyle EE, Cronin SJ, Cole SE, Thouret JC (2010) The coalescence and organization of lahars at Semeru volcano, Indonesia. *Bulletin of Volcanology* (72):961-970. https://doi.org/10.1007/s00445-010-0381-8

Doyle EE, Cronin SJ, Thouret JC (2011) Defining conditions for bulking and debulking lahars. *Geological Society of America Bulletin* 123(7-8):1234-1247. https://doi.org/10.1130/B30227.1

Ettinger S, Mounaud L, Magill C, Yao-Lafourcade AF, Thouret JC, Manville V, Negulescu C, Zuccaro G, De Gregorio D, Nardone S, Uchuchoque JAL, Arguedas A, Macedo L, Manrique Llerena, N (2015) Building vulnerability to hydrogeomorphic hazards: Estimating damage probability from qualitative vulnerability assessment using logistic regression. *Journal of hydrology* 541:563-581.
https://doi.org/10.1016/j.jhydrol.2015.04.017

Federico F, Amoruso A (2008) Simulation of mechanical effects due to the impact of fluid-like debris flows on structures. *Italian Journal of Engineering Geology and Environment* 8(1):5-24.
https://doi.org/10.4408/IJEGE.2008-01.O-01

Furtney J, Zhang F, Han Y (2013) Review of methods and applications for incorporating fluid flow in the Discrete Element Method. In: Proceedings of the 5th International FLAC/DEM Symposium. Itasca International Inc. (Minneapolis, MN), Hangzhou, P.R. China, 10 p.

He, X., Liang, D., Wu, W., Cai, G., Zhao, C., & Wang, S. (2018) Study of the interaction between dry granular
flows and rigid barriers with an SPH model. International Journal for Numerical and Analytical Methods in
Geomechanics, 42(11):1217-1234. https://doi.org/10.1002/nag.2782

Hu K, Wei F, Li Y (2011) Real-time measurement and preliminary analysis of debris-flow impact force at Jiangjia Ravine, China. Earth Surface Process and Landforms 36:1268-1278.
https://doi.org/10.1002/esp.2155

Hübl J, Suda J, Proske D, Kaitna R, Scheidl C (2009) Debris flow impact estimation. In: Proceedings of the 11th International Symposium on Water Management and Hydraulic Engineering, Ohrid (North Macedonia) (1-5 September), 137-148.

Hübl, J, Suda J, Nagl G, Rudolf-Miklau F (2017) Standardized stress model for design of torrential barriers under impact by Debris Flow (according to Austrian Standard Regulation 24801). International Journal of Erosion Control Engineering 10(1):47-55. https://doi.org/10.13101/ijece.10.47

Itasca, 2016. Particle Flow Code 5.0 Documentation.

Iverson RM, Costa JE, LaHusen RG (1992) Debris-flow flume at H.J. Andrews Experimental Forest, Oregon. U.S. Geological Survey, Open-file Report 92-483. https://doi.org/10.3133/ofr92483

Iverson RM (1997) The physics of debris flows. Reviews of Geophysics 35(3):245-296.
https://doi.org/10.1029/97RG00426

Iverson RM, Logan M, LaHusen RG, Berti M (2010) The perfect debris flow? Aggregated results from 28 large-scale experiments. Journal of Geophysical Research 115(F03005):29.
https://doi.org/10.1029/2009JF001514

Iverson RM, George DL (2014) A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis. Proceedings of the Royal Society A 470: 20130819.
https://doi.org/rspa.2013.0819

Kattel P, Kafle J, Fischer JT, Mergili M, Tuladhar BM, Pudasaini SP (2018) Interaction of two-phase debris flow with pillars. Engineering Geology 242:197-217. https://doi.org/10.1016/j.enggeo.2018.05.023

Lavigne F, Thouret JC (2000) Les lahars : dépôts, origines et dynamiques. Bulletin de la Société Géologique de France 1741(5):545-557. https://doi.org/10.2113/171.5.545

Leonardi A, Wittel FK, Mendoza M, Vetter R, Herrmann HJ (2016) Particle-Fluid-Structure Interaction for Debris Flow Impact on Flexible Barriers. Computer-Aided Civil and Infrastructure Engineering 31(5):323-333. https://doi.org/10.1111/mice.12165
Li X, Zhao J (2018) A Unified CFD-DEM approach for modeling of debris flow impacts on flexible barriers.
International Journal for Numerical and Analytical Methods in Geomechanics 42:1643-1670.
https://doi.org/10.1002/nag.2806

Liu, C., Sun, Q., Zhou, G.G.D., 2018. Coupling of material point method and discrete element method for granular flows impacting simulations. Int. J. Numer. Methods Eng. 115, 172–188. https://doi.org/10.1002/nme.5800

Major JJ (1995) Experimental studies of deposition at a debris-flow flume. U.S. Geolocial Survey, Fact sheet
028-94. https://doi.org/10.3133/fs02894

Major JJ, Bertin, D, Pierson TC, Amigo A, Iroume A, Ulloa H, Castro J (2016) Extraordinary sediment delivery and rapid geomorphic response following the 2008–2009 eruption of Chaiten Volcano, Chile. Water Resources Research 52:5075-5094. https://doi.org/10.1002/2015WR018250

Manville V, Major JJ, Fagents SA (2013) Modeling lahar behavior and hazards. In: Fagents SA, Gregg TKP, Lopes RMC (Eds) Modeling volcanic processes: the physics and mathematics of volcanism, chapter 14.
Cambridge University Press, Cambridge, pp. 300-330. https://doi.org/10.1017/CBO9781139021562.014

Mead SR, Magill C, Lemiale V, Thouret JC, Prakash M (2017) Examining the impact of lahars on buildings using numerical modeling. Natural Hazards and Earth System Sciences 17:703-719.
http://dx.doi.org/10.5194/nhess-17-703-2017

Moulinec C, Denis C, Pham C et al (2011) TELEMAC: An Efficient Hydrodynamics Suite for Massively Parallel Architectures Q. Computers & Fluids 51: 30–34.

Newhall CG and Punongbayan RS (1996) Fire and Mud: Eruptions and Lahars of Mount Pinatubo, Philippines, University of Washington Press, Seattle.

Paguican EMR, Lagmay AMF, Rodolfo KS, Rodolfo RS, Tengonciang AMP, Lapus MR, Baliatan EG, Obille Jr EC (2009) Extreme rainfall-induced lahars and dike breaching, 30 November 2006, Mayon Volcano, Philippines. Bulletin of Volcanology 71(8):845-857. http://dx.doi.org/10.1007/s00445-009-0268-8

Pallares C, Fabre D, Thouret JC, Bacconnet C, Charca-Chura JA, Martelli K, Talon A, Yanqui-Murillo C (2015) Geological and geotechnical characteristics of recent lahar deposits from El Misti volcano in the city area of Arequipa, South Peru. Geotechnical and Geological Engineering 33(3):641-660.
https://doi.org/10.1007/s10706-015-9848-x

Proske D, Suda J, Hübl J (2011) Debris Flow Impact Estimation for Breakers. Georisk 5(2): 143–155.

Pudasaini SP (2012) A general two-phase debris flow model. Journal of Geophysical Research: Earth Surface 117. https://doi.org/10.1029/2011JF002186
Rameshwaran P, Naden P et al (2013) Inter-Comparison and Validation of Computational Fluid Dynamics Codes in two-stage meandering channel flows. Applied Mathematical Modelling 37(20–21):8652–8672.
https://doi.org/10.1016/j.apm.2013.07.016

Scott KM, Vallance JW, Kerle N, Macias JL, Strauch W, Devoli G (2005) Catastrophic precipitation-triggered lahar at Casita volcano, Nicaragua: occurrence, bulking and transformation. Earth Surface Processes and Landforms 30:59-79. https://doi.org/10.1002/esp.1127

Scheidl, C., Chiari, M., Kaitna, R., Müllegger, M., Krawtschuk, A., Zimmermann, T., Proske, D., 2013. Analysing Debris-Flow Impact Models, Based on a Small Scale Modelling Approach. Surv. Geophys. 34, 121–140.
https://doi.org/10.1007/s10712-012-9199-6

Smith, GA (1986) Coarse-Grained Nonmarine Volcaniclastic Sediment: Terminology and Depositional Process. Geological Society of America Bulletin 97(1):1–10.

Sosio, R., Crosta, G. B., & Frattini, P. (2007) Field observations, rheological testing and numerical modelling of a debris-flow event. Earth Surface Processes and Landforms, 32(2):290–306.

Teufelsbauer H, Wang Y, Pudasaini SP, Borja RI, Wu W (2011) DEM simulation of impact force exerted by granular flow on rigid structures. Acta Geotechnica 6:119-133
https://doi.org/10.1007/s11440-011-0140-9

Thouret JC, Enjolras G, Martelli K, Santoni O, Luque JA, Nagata M, Argedas A, Macedo L (2013) Combining criteria for delineating lahar- and flash flood-prone hazard and risk zones for the city of Arequipa, Peru. Natural Hazards and Earth System Sciences 13(2):339-360. https://doi.org/10.5194/nhess-13-339-2013

Thouret JC, Ettinger S, Guitton M, Santoni O, Magill C, Martelli K, Arguedas A (2014) Assessing physical vulnerability in large cities exposed to flash floods and debris flows: the case of Arequipa (Peru). Natural Hazards 73 (3):1771-1815. https://doi.org/10.1007/s11069-014-1172-x

Thouret JC, Antoine S, Magill C, Ollier C (2020) Lahars and Debris Flows: Characteristics and Impacts. Earth-Science Reviews 201. https://doi.org/10.1016/j.earscirev.2019.103003

Tiberghien D, Laigle D, Naaim M, Thibert E, Ousset F (2007). Experimental investigations of interaction between mudflow and an obstacle. In: Proceedings of the Fourth International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment (DFHM-4), Chengdu (China) (September 10-13), 281-292.

Toyos G., Gunasekera R., Zanchetta G., Oppenheimer C., Sulpizio R., Favalli M., Pareschi M. T. (2008). GIS-assisted modelling for debris flow hazard assessment based on the events of May 1998 in the area of Sarno,
Southern Italy: II. Velocity and dynamic pressure. Earth Surface Processes and Landforms: The Journal of the British Geomorphological Research Group, 33(11):1693-1708.

Trujillo-Vela, M. G., Galindo-Torres, S. A., Zhang, X., Ramos-Cañón, A. M., & Escobar-Vargas, J. A. (2020). Smooth particle hydrodynamics and discrete element method coupling scheme for the simulation of debris flows. Computers and Geotechnics, 125: 103669.

Vallance JW (2000) Lahars. In: Sigurdsson H, Houghton BF, McNutt SR, Rymer H, Stix J (eds) Encyclopedia of volcanoes. Academic Press, San Diego, pp 601–616

Vallance JW, Iverson R et al (2015) Lahars and their deposits. In: Sigurdsson H (ed) The Encyclopedia of Volcanoes, 2nd edn. Academic Press, San Diego, pp 649-664.

Wang, D., Chen, Z., He, S., Liu, Y., & Tang, H. (2018) Measuring and estimating the impact pressure of debris flows on bridge piers based on large-scale laboratory experiments. Landslides, 15(7):1331–1345.
https://doi.org/10.1007/s10346-018-0944-x

Zanchetta G, Sulpizio R, Pareschi MT, Leoni FM, Santacroce R (2004) Characteristics of May 5-6, 1998 volcaniclastic debris flows in the Sarno area (Campania, southern Italy): relationships to structural damage and hazard zonation. Journal of Volcanology and Geothermal Research 133:377-393.
https://doi.org/10.1016/S0377-0273(03)00409-8

Zhang S (1993) A comprehensive approach to the observation and prevention of debris flows in China. Natural Hazards 7:1–23. https://doi.org/10.1007/BF00595676

Zhao T, Houlsby GT, Utili S (2014) Investigation of Granular Batch Sedimentation via DEM-CFD Coupling. Granular Matter 16(6):921–932.

Zhao, L., He, J.W., Yu, Z.X., Liu, Y.P., Zhou, Z.H., Chan, S.L., 2020. Coupled numerical simulation of a flexible barrier impacted by debris flow with boulders in front. Landslides 17, 2723–2736.
https://doi.org/10.1007/s10346-020-01463-x
## TABLE CAPTIONS

**Table 1.** A selection of experimental studies evaluating forces induced by debris flows (DFs) or lahars (LHs).

**Table 2.** Characteristics of the debris flow chosen for the simulated flow runs (this study).

**Table 3.** Results of the model calibration.

**Table 4.** Analysis of model calibration results.

**Table 5.** Details of the DF cases considered in the parametric study.

## FIGURE CAPTIONS

**Fig. 1.** Flow chart displaying each of the stages through the simulation process.

**Fig. 2.** Sketch showing the numerical channel. At the top of the channel, a supply chamber is implemented to generate the particles.

**Fig. 3.** Results obtained for the final case study: (a) flow velocity and (b) discharge versus simulation time. The figure (3.a) shows the variation of the average particle velocity in two sections of the ravine as a function of time: the curve stabilizes at 2.8 m/s.

**Fig. 4.** Top view of the channel bed showing the modified geometry with the presence of a pillar acting as an obstacle.

**Fig. 5.** Process model and impact model of DFs on an obstacle from Proske et al. (2011).

**Fig. 6.** Top view of the results of the fluid phase modelling in a channel including a pillar after a simulation duration of 12 s: (a) flow height (m), (b) velocity field of fluid phase (m/s).

**Fig. 7.** Stress applied by particles on a pillar versus simulation time. The black curve represents raw data and the red curve shows a moving average window (filter) with k= 100.

**Fig. 8.** The process of the particles flowage upstream the pillar: (a) before interaction, (b) start of progressive rise of the particles upstream of the pillar, (c) continued gradual rise of the particles upstream of the pillar, (d) final state.

**Fig. 9.** Mean stress applied by particles (grey curve), fluid phase (blue curve) and total pressure (black dashed curve) on the pillar (reference case).

**Fig. 10.** Normalised stress induced by the DF (y-axis) versus normalised flow height on the pillar (x-axis), averaged between t= 9 s and t= 12 s, for 3 different flow heights $h_{Mu}$.

**Fig. 11.** Normalised stress induced by the DF (y-axis) versus normalised flow height on the pillar (x-axis), averaged between t= 9 s and t= 12 s for three different flow velocities V and for $h_{Mu}$ = 1.5 m at the entrance of channel.
**Fig. 12.** Maximum dimensionless pressure versus Froude number, compared with 4 empirical models (the standard deviation for maximal numerical pressure is shown).
<table><caption>**Table 1.** A selection of experimental studies evaluating forces induced by debris flows (DFs) or lahars (LHs)</caption><thead><tr><td>Reference</td><td>Study area</td><td>Method</td><td>Type of flow and obstacle</td><td>Tool</td><td>Measured parameters</td></tr></thead><tbody><tr><td>Zanchetta et al. (2004)</td><td>May 5-6, 1998 volcanic debris flows, Sarno area, (Italy)</td><td>Analytical and field observations</td><td>DF / LH on rigid structures (buildings, bridges)</td><td>Field observations, measurements from the study area and empirical equations</td><td>Hydrostatic pressure and dynamic overpressure, velocity, etc.</td></tr><tr><td>Tiberghien et al. (2007)</td><td>Laboratory flume, LTHE Grenoble</td><td>Small-scale experiment</td><td>Mudflow on fixed barrier</td><td>Pressure sensor installed inside the pillar body</td><td>Impact pressure total</td></tr><tr><td>Iverson et al. (2010)</td><td>USGS debris-flow flume at the H.J. Andrews Experimental Forest near Blue River, Oregon, USA</td><td>Large-scale experiments</td><td>DF in different cross sections of the flume</td><td>Electronic sensors along the flume centreline</td><td>Flow thickness, basal normal stresses and basal pore fluid pressures</td></tr><tr><td>Hu et al. (2011)</td><td>Jiangjia Ravine, China, August 25, 2004</td><td>Experimental and field measurements</td><td>DF on steel pile foundation</td><td>Three sensors at different flow depths</td><td>Grain impact loading and fluid pressure</td></tr><tr><td>Bugnion et al. (2012)</td><td>Channel constructed on the side of a rock quarry, near Veltheim, Switzerland</td><td>Experimental, in situ channel</td><td>DF on small obstacles</td><td>Two impact sensors installed in situ</td><td>Pressure (total)</td></tr></tbody></table>
<table><caption>**Table 2.** Characteristics of natural debris flows chosen for the simulated flow (this study).</caption><tr><td>Solid fraction</td><td>50 to 55 vol%</td></tr><tr><td>Density</td><td>1800 to 1900 kg/m³</td></tr><tr><td>Discharge</td><td>35 to 55 m³/s</td></tr><tr><td>Velocity</td><td>3 m/s</td></tr><tr><td>Flow height</td><td>1.5 m</td></tr><tr><td>Froude number</td><td>0.78</td></tr></table>
<table><caption>Table 3. Results of the model calibration</caption><thead><tr><th rowspan="3">run #</th><th colspan="6">Model parameters</th><th colspan="6">Simulation results in terms of flow characteristics</th></tr><tr><th>Particle density</th><th>Fluid density</th><th colspan="2">Friction coefficient</th><th>Dynamic viscosity</th><th>Rolling resistance</th><th>Average velocity</th><th>Average bulk density</th><th>Discharge</th><th>Solid concentration</th><th>Average velocity</th><th>Average bulk density</th><th>Discharge</th><th>Solid concentration</th></tr><tr><th>kg/m³</th><th>kg/m³</th><th>particle - particle</th><th>particle - wall</th><th>Pa.s</th><th>-</th><th>m/s</th><th>kg/m³</th><th>m³/s</th><th>vol%</th><th colspan="4">% of variation to the reference case study</th></tr></thead><tbody><tr><td>ref.</td><td>2700</td><td>1500</td><td>0.20</td><td>0</td><td>0.048</td><td>0</td><td>2.80</td><td>2139.6</td><td>45</td><td>53.3</td><td></td><td></td><td></td><td></td></tr><tr><td>1</td><td>2500</td><td>1500</td><td>0.20</td><td>0</td><td>0.048</td><td>0</td><td>2.80</td><td>2029.0</td><td>45</td><td>52.9</td><td>0.0 %</td><td>-5.2 %</td><td>0.0 %</td><td>-0.8 %</td></tr><tr><td>2</td><td>2700</td><td>1000</td><td>0.20</td><td>0</td><td>0.048</td><td>0</td><td>2.55</td><td>1921.4</td><td>43</td><td>54.2</td><td>-8.9 %</td><td>-10.2 %</td><td>-4.4 %</td><td>+1.7 %</td></tr><tr><td>3</td><td>2700</td><td>2000</td><td>0.20</td><td>0</td><td>0.048</td><td>0</td><td>2.90</td><td>2361.9</td><td>47</td><td>51.7</td><td>+3.6 %</td><td>+10.4 %</td><td>+4.4 %</td><td>-3.0 %</td></tr><tr><td>4</td><td>2700</td><td>1500</td><td>0.05</td><td>0</td><td>0.048</td><td>0</td><td>2.80</td><td>2140.8</td><td>45</td><td>53.4</td><td>0.0 %</td><td>+0.1 %</td><td>0.0 %</td><td>+0.2 %</td></tr><tr><td>5</td><td>2700</td><td>1500</td><td>0.40</td><td>0</td><td>0.048</td><td>0</td><td>2.80</td><td>2132.4</td><td>45</td><td>52.7</td><td>0.0 %</td><td>-0.3 %</td><td>0.0 %</td><td>-1.1 %</td></tr><tr><td>6</td><td>2700</td><td>1500</td><td>0.20</td><td>0.5</td><td>0.048</td><td>0</td><td>2.80</td><td>2122.8</td><td>45</td><td>51.9</td><td>0.0 %</td><td>-0.8 %</td><td>0.0 %</td><td>-2.6 %</td></tr><tr><td>7</td><td>2700</td><td>1500</td><td>0.20</td><td>1.0</td><td>0.048</td><td>0</td><td>2.80</td><td>2131.2</td><td>45</td><td>52.6</td><td>0.0 %</td><td>-0.4 %</td><td>0.0 %</td><td>-1.3 %</td></tr><tr><td>8</td><td>2700</td><td>1500</td><td>0.20</td><td>0</td><td>0.030</td><td>0</td><td>2.60</td><td>2142.0</td><td>46</td><td>53.5</td><td>-7.1 %</td><td>+0.1 %</td><td>+2.2 %</td><td>+0.4 %</td></tr><tr><td>9</td><td>2700</td><td>1500</td><td>0.20</td><td>0</td><td>0.075</td><td>0</td><td>2.90</td><td>2138.4</td><td>46</td><td>53.2</td><td>+3.6 %</td><td>-0.1 %</td><td>+2.2 %</td><td>-0.2 %</td></tr><tr><td>10</td><td>2700</td><td>1500</td><td>0.20</td><td>0</td><td>0.048</td><td>0.6</td><td>2.60</td><td>2032.8</td><td>55</td><td>44.4</td><td>-7.1 %</td><td>-5.0 %</td><td>+22.2 %</td><td>-16.7 %</td></tr><tr><td colspan="16">Final model</td></tr><tr><td>11</td><td>2500</td><td>1100</td><td>0.4</td><td>0.0</td><td>0.048</td><td>0.2</td><td>2.80</td><td>1867</td><td>40</td><td>54.7</td><td>0.0 %</td><td>-12.7 %</td><td>-11.1 %</td><td>+2.6 %</td></tr></tbody></table>
<table><caption>**Table 4** Analysis of model calibration results</caption><thead><tr><td>run #</td><td>Parameter</td><td>Observations</td></tr></thead><tbody><tr><td>reference case and 1</td><td>Density of particles</td><td>The average mean flow velocity is similar. The bulk density increases inasmuch the density of the particles is higher. The discharge does not vary.</td></tr><tr><td>ref. case, 2 and 3</td><td>Fluid density</td><td>The fluid density is used to calculate the drag force: higher the fluid density is, greater the force applied to the particle is, hence the flow velocity increases. So, the increase in flow velocity with higher fluid density leads to the increase in discharge. The bulk density increases with higher fluid density.</td></tr><tr><td>ref. case, 4 and 5</td><td>Friction between particle-particle</td><td>The influence of this parameter is negligible on the overall results given the range of values considered.</td></tr><tr><td>ref. case, 6 and 7</td><td>Friction between wall and particle</td><td>The influence of this parameter is negligible on the overall results in the range of values considered here.</td></tr><tr><td>ref. case, 8 and 9</td><td>Dynamic viscosity</td><td>The average flow velocity slightly increases with the increasing dynamic viscosity due to the increase in the drag force applied to the particle. The discharge does not vary.</td></tr><tr><td>ref. case and 10</td><td>Rolling resistance</td><td>Rolling resistance parameter mimics effects of angular shape on spherical particles: particles tend to roll less and eventually become blocked, a fact which explains the decrease in flow velocity. On the other hand, the surface area increases (the model shows a greater flow height, which explains the increase in flow rate).</td></tr></tbody></table>
<table><caption>**Table 5.** Details of the DF cases considered in the parametric study.</caption><thead><tr><th>Case</th><th>Flow height (m)</th><th>Flow velocity (m/s)</th><th>Froude number</th><th>Time averaged max. pressure / standard deviation (kPa)</th><th>Time averaged mean pressure / standard deviation (kPa)</th><th>h<sub>Mu,2</sub> (m)</th></tr></thead><tbody><tr><td>Reference</td><td>1.5</td><td>3</td><td>0.78</td><td>119/34</td><td>74/16</td><td>2.4</td></tr><tr><td>Case No. 1</td><td>3</td><td>3</td><td>0.55</td><td>258/59</td><td>145/22</td><td>3.8</td></tr><tr><td>Case No. 2</td><td>4</td><td>3</td><td>0.49</td><td>349/73</td><td>187/26</td><td>4.8</td></tr><tr><td>Case No. 3</td><td>1.5</td><td>4.5</td><td>1.18</td><td>257/69</td><td>126/28</td><td>2.6</td></tr><tr><td>Case No. 4</td><td>1.5</td><td>6</td><td>1.56</td><td>377/83</td><td>162/26</td><td>3.2</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12</figcaption></figure>