Discrete modelling of debris flows for evaluating impacts on structures

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.


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 and (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 min) 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. 2010Doyle et al. , 2011Manville 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) and 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(Thouret et al. , 2014Ettinger 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 analysing 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 miniaturised 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 setup 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 illustrate 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.

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. "Physical and mechanical parameters of the model" section 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 Fig. 1 describes each of the steps of the simulation process.

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. 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. 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: where f is the density of the fluid phase, v i is the volume of the immersed particle, g is the gravitational acceleration and z is 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): where f is the density of the fluid phase (kg/m 3 ), d is the particle diameter (m), C d is the drag coefficient, V f is the velocity vector of the fluid and V b is the velocity of the immersed particle (m/s). The drag coefficient C d can be written as follows Zhao et al. (2014): where Re p is the particle Reynolds number, expressed as As proposed in Zhao et al. (2014), a correction factor n − +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: where is a term ranging between 3.4 and 3.7, expressed as

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 . 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 where Δ b is the relative bend-rotation increment and k r is the rolling resistance stiffness. The rolling resistance stiffness is calculated as follows: where k s is the contact stiffness in the tangential direction and R is the contact effective radius of the two elements in contact: The magnitude of the updated rolling resistance moment is then checked against the threshold limit M * given by where rr is the rolling resistance friction coefficient and F n is the normal linear contact force. The values of contact model parameters are discussed in "Physical and mechanical parameters of the model" section.

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

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

Physical and mechanical parameters of the model
As mentioned in "Contact laws" section, 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 and 55 m 3 /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 and (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 9500 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 × 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 whilst 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 3 , 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 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 3 . 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.

Model calibration
The simulation of the flow lasts 12 s, 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 where bulk is the apparent density of the flow, s is the particle density and 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, whilst 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 and  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; and (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  and 4), a strategy to calibrate the model was developed: we started by setting the rolling resistance parameters which influence 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 3 to meet an apparent bulk density of the flow around 1900 kg/ m 3 . We were then able to set the fluid density to a value of 1100 kg/m 3 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 3 and a volume discharge of about 40 m 3 /s (Fig. 3).

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 Sect. (1 m × 1 m). Fluid density The fluid density is used to calculate the drag force: the higher the fluid density is, the 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 Ref. case, 4 and 5 Friction between particle-particle The influence of this parameter is negligible on the overall results given the range of values considered Ref. case, 6 and 7 Friction between wall and particle The influence of this parameter is negligible on the overall results in the range of values considered here Ref. case, 8 and 9 Dynamic viscosity 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 Ref. case and 10 Rolling resistance 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) 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)

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 Eq. (4). Figure 6 shows the free surface (a) and the velocity field (b) in the part of the channel (15 m in front of the pillar but only 1.5 m in the back). We conducted the DEM simulation with fluid interactions over a duration of 12 s.

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.

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 h. The simulations were performed on computer with a processor Intel (R) Core (TM)i9-9900 CPU and 32Go RAM (DDR4, 2666 MHz).
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 analyse 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 stabilise at around 74 kPa. Figure 8 shows a longitudinal cross section of the model flow for different durations between 4 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 particle 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(Hübl et al. , 2017: where 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: 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 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 s to t = 12 s), (ii) the static fluid stress profile and (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. 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 and (d) final state 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. 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 Mu V DF 2 , with Mu = 1867 kg·m −3 . We note p = P MAX 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, p = 21 ; for h Mu = 3 m, p = 16 ; and for h Mu = 1.5 m, p = 7 . Froude number varies from 0.49 to 0.78 corresponding to the variation of flow height. The maximal normalised pressure, or maximal dynamic pressure coefficient 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).

Effect of flow heigh
A great drop in reduced pressure p can be observed for the range of Fr 0.49-0.78. The graph (Fig. 10) also indicates the maximal DF height against the pillar h Mu,2 (including accumulation), which explains why ratios z/h Mu are greater than 1. 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 DF height against the pillar h Mu,2 (including accumulation), which are very close in these three cases.

Effect of flow velocity
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 varies from 0.78 to 1.56 corresponding to the effect of flow velocity. The distribution of normalised dynamic pressure coefficient on the  Cui et al. (2015) and Wang et al. (2018): the reduced pressure coefficient 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).

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

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;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, whilst 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 andThouret 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: where is the debris flow density, h Mu,2 is the maximal debris flow height against the pillar (including accumulation), V DF is the flow velocity and c m is 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: where 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 = 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: where 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 pressurep 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 = 3m∕s, p max = 150kPa ; for V DF = 4.5 m/s,p max = 217kPa ; and for V DF = 6.0 m/ s, p max = 310kPa . These maximal pressures are in accordance with the values provided by numerical model (Table 5).

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, (15) p = c d V 2 cos 2 (16) p max = 5 mu .V DF 0.8 gh Mu,2 0.6 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-flow 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;Zhao et al. 2020). In most cases, modelling 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).

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 and size of the largest solid particles. We considered the large-size solid part included in the flow mixture as particles, whilst 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.