# Towards viscous debris flow simulation using DualSPHysics v5.2: internal behaviour of viscous flows and mixtures

Suzanne Lapillonne, Georgios Fourtakas, Vincent Richefeu, Guillaume Piton, Guillaume Chambon

► To cite this version:

Suzanne Lapillonne, Georgios Fourtakas, Vincent Richefeu, Guillaume Piton, Guillaume Chambon. Towards viscous debris flow simulation using DualSPHysics v5.2: internal behaviour of viscous flows and mixtures. Geoscientific Model Development, 2025, 18 (19), pp.7059-7075. 10.5194/gmd-18-7059-2025. hal-05316843

HAL Id: hal-05316843

https://hal.science/hal-05316843v1

Submitted on 15 Oct 2025

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
<figure><img src="image_2.png" /></figure>

# Towards viscous debris flow simulation using DualSPHysics v5.2: internal behaviour of viscous flows and mixtures

Suzanne Lapillonne¹, Georgios Fourtakas², Vincent Richefeu³, Guillaume Piton¹, and Guillaume Chambon¹

¹Univ. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, Grenoble, France

²School of Engineering, University of Manchester, Manchester, M13 9PL, UK

³Univ. Grenoble Alpes, 3SR, Gières, France

**Correspondence:** Suzanne Lapillonne (lapillonne.s@gmail.com)

Received: 19 July 2024 – Discussion started: 22 October 2024

Revised: 23 June 2025 – Accepted: 27 July 2025 – Published: 10 October 2025

**Abstract.** This paper investigates the accuracy of a solid-fluid model using the smoothed-particle hydrodynamics (SPH) software DualSPHysics v5.2 coupled with ProjectChrono for viscous debris flow modelling. It focuses on different validation steps of the method, both for pure fluid and a mixture of fluid and boulders, to build the reliability of the model in order to prepare for the simulation of a simplified debris flow. First, velocity profiles, the free-surface shape and the velocity of surges of a viscous fluid are validated against well-documented experimental data. This is, to the best of our knowledge, one of the few validations of the SPH approach for very viscous flows near the creeping threshold (Re ≈ 1). Second, the influence on the macroscopic viscosity of the introduction of granular elements in a viscous fluid is studied against a semiempirical formula. Finally, the method is applied to simplified 2D debris flow surges with field-like features. Surges modelled in this paper are composed of a viscous Newtonian fluid and polydisperse boulders. The flow of surges of different concentrations is studied, and Froude numbers of real field measurements are retrieved. Such complex models are shown to be relevant to the study of debris flow dynamics.

## 1 Introduction

Debris flows are rapid flows saturated with non-plastic debris in a steep channel (Hungr, 2005). These fast flows yield suddenly, behaving as so-called surges, creating a granular front, which has the potential to be very destructive, followed by a viscous matrix engulfing granular material. Debris flows

evolve naturally in steep, erosion-prone catchments and can be triggered by abundant runoff (Bel, 2017; Bernard et al., 2025), landslides (Iverson, 1997; Recking et al., 2013), or snow and ice melt (Recking et al., 2013). Once initiated, the flow propagates downstream, often recruiting material from the channel through entrainment (Simoni et al., 2020; Reid et al., 2016). In the European Alps, the material transported by the flow usually comes from the weathering of mountain hillslopes (Recking et al., 2013). This leads to the presence of viscous-dominated debris flow, with a high clay content in the interstitial fluid and granular materials of a wide verity of sizes, from sand to boulders (Coussot et al., 1998). Efforts of the scientific community in the debris flow field are growing in every direction, some focusing on triggering conditions and predictability (e.g. Bel et al., 2017; Berti et al., 2012) and others on impact forces and hazard mitigation infrastructures (e.g. Alibaba et al., 2015; Poudyal et al., 2019) or the internal dynamics within the flow (e.g. Laigle and Labbe, 2017; Leonardi et al., 2014). To better mitigate the risks associated with debris flows, there is a need to better understand the internal mechanics of debris flows and their interactions with infrastructure. Due to their destructive nature, it is not possible undertake measurements directly within the flow, and experimental small-scale facilities cannot fully replicate how debris flows propagate and maintain their dynamics, due to scale effects. This has led to the scientific community turning towards numerical methods as a relevant mode of action to describe the internal dynamics of debris flow motion.

Strategies to model debris flow phenomena are adapted to the scale of interest. First, event-scale models replicate the debris flow runout and spreading at a large scale, span-
ning over the entire reach. Such models are typically used as an input for the design of hazard maps on a debris fan. Event-scale models rely on depth-averaged methods. Depth-averaged models can either be single phase or multiphase. The former models consider the debris flow material as one homogenized phase. The rheology of the phase differs between models, with some using a Voellmy rheology (e.g. Christen et al., 2010; Pastor et al., 2009; Pirulli et al., 2005) and others using a Herschel-Bulkley or Bingham rheology (e.g. Laigle and Coussot, 1997; Gibson et al., 2021). Such models are typically user-friendly and, thus, popular for engineering applications. RAMMS (Christen et al., 2010) and HEC-RAS 2D (Gibson et al., 2021) are two of the common examples used in engineering for debris flow modelling. However, because these are single-phase models, there is no change in density or rheological parameters along the flow path or duration; thus, there is a lack of representation of the local physical processes that impact the macroscopic behaviour of the flow. To improve this aspect, multiphase depth-averaged models have started to bloom over the past decade. They typically depict the debris flow material with one phase representing the liquid contribution and one or two phases representing the solid contribution. They are then able to represent local effects within the phase interactions, such as dilatancy, pore-fluid pressure, entrainment, de-watering or washing over deposition. The most common models employing this method are D-CLAW (George and Iverson, 2014), r.avaflow 2.0 (Mergili and Pudasaini, 2014), Titan2D (Pastor et al., 2018) and the recently improved version of RAMMS (Meyrat et al., 2022). However, depth-averaged models are, by nature, only adapted to large-scale representation of the flow, as they intrinsically neglect the momentum changes in depth; thus, when dealing with the mechanical changes inside the debris flow material during the flow, they become limited. The large boulders that are transported by debris flows cannot be explicitly described in such models.

Conversely, surge-scale models or models with smaller spatial domains aim to represent the physics of debris flow dynamics at a mesoscopic and microscopic scale in order to better understand the internal motion of the material and its interaction with infrastructure. They rely on explicit models where the velocity gradients in the depth of the flow are solved (e.g. Albaba et al., 2015; Laigle and Labbe, 2017; Chambon et al., 2011; Mitsoulis et al., 1993; Roquet and Saramito, 2003; Kong et al., 2021; Leonardi et al., 2014). For such a strategy, models can be classified into three types: (i) fluid models, modelling the debris flow material as a continuum of fluid with adapted rheology; (ii) granular models, which model the debris flow material as an ensemble of grains; and (iii) coupled methods, which represent both phases.

First, fluid models represent debris flows as a continuous fluid with complex rheology (Chambon et al., 2011; Labbé, 2015; Mitsoulis et al., 1993; Rodriguez-Paz and Bonet, 2004; Roquet and Saramito, 2003). The rheology can be fitted to

observations to recover macroscopic flow features. Time-independent rheological constitutive models have proven to be effective to model debris flows when the volumetric fractions of silts and clay is sufficient and local discontinuities, segregation and concentration heterogeneity can be neglected (Ancey, 2003; Laigle and Labbe, 2017). As a first approximation, muddy debris flows, which are mechanically driven by viscous motion, can be studied at the macroscopic scale by considering a homogenous flow of a Newtonian viscous fluid (Han et al., 2014) or of a more complex material such as a Herschel-Bulkley fluid (Coussot and Meunier, 1996; Chambon et al., 2011). Such methods are often able to reproduce experimental results (Laigle and Labbe, 2017; Chambon et al., 2011) and retrieve characteristics of the flow that are close to field reality, with Froude numbers ranging from 0.5 to 3 (Marchi et al., 2002; Lapillonne et al., 2023; McArdell et al., 2023). These methods perform well for studies of interactions with obstacles or spreading onto a deposition area; however, they neglect or simplify the granular skeleton in the flow. While the granular phase may be represented rheologically, the design of such models lacks the input of solid-solid interactions of the granular content onto the macroscopic flow. Thus, they cannot explore the role of the granular skeleton in debris flow motion (Kaitna et al., 2011) or conditions related to barrier clogging (Piton et al., 2022).

Second, granular models solely represent the granular phase of the debris flow material (Albaba et al., 2015; Cecato et al., 2018), omitting the fluid input or strongly simplifying it by adapting the force interactions between the grains. They mainly focus on the estimation of impact forces of the boulders carried by the flow (Albaba et al., 2015; Cecato et al., 2018). These models are effective to study impact forces applied to an infrastructure and the general rearrangement of the granular skeleton in a granular flow, but the absence of resolved interstitial fluid makes their parallel with real debris flows difficult, especially when it comes to studying established flows. Due to their granular nature, they also tend to require very steep flowing conditions with highly supercritical Froude numbers (up to 7). Such steep sites exist, but they seldom have assets located nearby or mitigation infrastructure, as those are typically rather located near alluvial fans or valley bottom where the slopes are milder.

Finally, thanks to scientific and technological advances, coupled methods have grown in the past years. These few hybrid methods are promising for the following reasons: resulting forces are computed with all inputs, and the interplay between large grains and the interstitial fluid is represented. They are able to explore larger scientific questions, but they require thorough validation. Kong et al. (2021) modelled debris flows using a coupled solver with computational fluid dynamics (CFD) for the fluid phase and the discrete element method (DEM) for the granular phase. The model was validated against analytical predictions, field measurements and experimental results for the impact load on the barrier
(Li et al., 2021). They were able to investigate the motion of the debris flow material over a wide range of Froude numbers (0.6–8.7). However, owing to the numerical complexity of CFD–DEM solvers, the computational cost of such models is relatively high. Due to this technical inconvenience, they relied on initializing the velocity profiles of the flow as uniform close to the obstacle (1 m upstream). While the no-slip boundary condition led to a retrieval of a velocity gradient, the distance at which the uniform velocity profile is initialized has an unknown effect on the results. Thus, the model renders interesting results for barrier impact investigation but could not be applied to study debris flow motion on a larger timescale. Other examples of models studying the gravity-driven free-surface flow of a mixture of large solid particles and fluid, such as those outlined in Peng et al. (2021) and Robb et al. (2016), can also show the high potential of CFD–DEM to fully represent mechanics of such mixtures.

Canelas et al. (2017) and Trujillo-Vela et al. (2020) employed the smoothed-particle hydrodynamics (SPH) framework with DEM. Canelas et al. (2017) studied the flow of granular material and water, which could bear some resemblance to granular debris flows. Their study contributed to understanding how this flow interacts with a filtering dam structure. However, they used a low-viscosity fluid to represent the interstitial fluid, which distinguishes their work from the alpine debris flows that this paper addresses. Trujillo-Vela et al. (2020) added a soil phase with a Drucker–Prager criterion to study the runout distance of a mixture of soil, water and a few (3) boulders. They were able to validate their pure fluid behaviours accurately against analytical derivations. However, their low boulder concentration does not allow for granular macroscopic effects to be retrieved, such as the crystallization of a granular skeleton or the viscosity change due to the introduction of boulders.

Finally, Leonardi et al. (2014) investigated a lattice Boltzmann solver coupled with DEM for debris flow modelling. They successfully validated the internal dynamics of the fluid during motion, the macroscopic viscosity response of the introduction of boulders in the fluid, and the macroscopic response of granular collapse (Leonardi et al., 2014) and flow fronts (Leonardi et al., 2015). This method is able to investigate the complexity of the debris flow surge but has a significantly high sophistication and computational cost.

In this work, we build upon these attempts towards a model for viscous debris flow modelling using a coupled solid–fluid method. The work presented in this paper aims to conduct a similar approach to Leonardi et al. (2014) with a model using SPH for the fluid phase, in order to obtain a fully validated model that has a lower sophistication, so that the computational time and repeatability become more sensible for engineering and exploratory uses. Debris flows considered in this paper fall under the viscous debris flow category (Takahashi, 2014), common in the French Alps. They are modelled with a simplified composition: a solid phase capturing the bigger sediment fraction and an inter-

stitial fluid that is typically highly viscous (Coussot and Munier, 1996; Bardou, 2002). A Lagrangian method was chosen to be more appropriate, as such methods deal straightforwardly with interfaces both at the free surface (especially for the steep free surface at the front of debris flow surges) and the granular material. In this paper, the SPH method was chosen because of its simplicity to describe the free surface and the interactions with solid elements. The DualSPHysics solver (Domínguez et al., 2021) is used because it is a widely available open-source software and has already been coupled with ProjectChrono (Tasora et al., 2016), a solid–solid solver based on the differential variational inequality (DVI) method that can model collisions of boulders. The solver was also chosen due to the following practical arguments: it is highly documented online, it is GPU- and CPU-accelerated, it is open source, it has a relatively low computational time (compared to traditional Eulerian methods and DEM) and it has a growing user community. This makes it a great candidate to be used more widely in studies of mass movements. However, the DualSPHysics solver is rarely used to model very viscous flows, and it has not been validated against experimental or field observations. This papers seeks to check the reliability of this tool by assessing it against validation data. Future users will then better know the capabilities and limitations of their solver.

In this paper, debris flow material is simplified as a viscous Newtonian fluid loaded with granular elements for which explicit collisions are computed. There are no macroscopic computations representing collision through rheology, and the fluid phase is purely Newtonian. First, the paper presents the methodology of the construction of the model. Second, we investigate the performance of the software regarding the representation of viscous surges, with respect to both the macroscopic and internal behaviour of the flow. We reproduce numerical experimental data from Freydier et al. (2017) and match velocity profiles and free-surface shapes to the experimental results. Such precise validation of SPH on laminar viscous fronts is rare, but it is necessary to move forward in the modelling of viscous geophysical and industrial flows. Third, the macroscopic viscosity of liquid–solid mixtures for different volumetric fractions is validated against the Krieger–Dougherty semiempirical equation (Krieger and Dougherty, 1959). Finally, we leverage the fact that the method is able to both represent slow laminar flow dominated by the viscous regime and yield accurate changes in the macroscopic viscosity when the addition of boulders is explored to investigate a boulder-laden flow surge. In the last section, we present exploratory results to model a surge-scale 2D debris flow (with respect to width) with an SPH-based hybrid model, investigating three different solid concentrations and their influence on the macroscopic behaviour of the surge.
## 2 Materials and methods

### 2.1 SPH modelling method

Smoothed-particle hydrodynamics is a computational fluid dynamics method based on the Lagrangian framework. In SPH, the continuous domain is discretized into numerical nodes (“particles”), which are points of known information. Typical properties of the continuum (e.g. velocity and density) are associated with each of these points. The SPH method relies on the resolution of the Navier–Stokes equation via interpolation onto these nodes. Particles interact with each other in a defined neighbourhood, named a smoothing kernel, by resolving the Navier–Stokes equation system.

#### 2.1.1 Governing equations

In SPH, the Navier–Stokes equations for a weakly compressible fluid are used. In their Lagrangian form, they can be written as the following system (Liu and Liu, 2003):

$$ \frac{D\rho}{Dt} = -\rho\nabla \cdot \underline{u}, \qquad (1) $$

$$ \frac{Du}{Dt} = \frac{1}{\rho}\nabla P + \frac{1}{\rho}\nabla \cdot \underline{\tau} + \underline{g}, \qquad (2) $$

where $\rho$ is the density, $\underline{u}$ is the velocity, $P$ is the pressure, $\underline{\tau}$ is the viscous stress tensor, $\underline{g}$ represents body forces and $\frac{D}{Dt}$ stands for the material derivative. Note that (as notation) we use ... for vectors and ... for tensors. To close this system, an equation of state is used:

$$ P = \frac{c_0^2 \rho_0}{\gamma} \left( \left( \frac{\rho}{\rho_0} \right)^{\gamma} - 1 \right), \qquad (3) $$

where $c_0$ is the numerical speed of sound, $\gamma = 7$ and $\rho_0$ is the reference density. This compressibility of SPH is purely artificial and is required because of the explicit resolution of the Navier–Stokes equation.

For any point $\underline{r}$ in an ensemble $\Omega$ and a function $f: \mathbb{R}^3 \to \mathbb{R}$ continuous over $\Omega$, the following can always be written:

$$ f(\underline{r}) = \int_{\Omega} \delta(\underline{r} - \underline{r}') f(\underline{r}') d\underline{r}', \qquad (4) $$

where $\delta$ is the Dirac function. In continuous SPH, the key idea is to substitute the value of $f(\underline{r})$ with $\langle f(\underline{r}) \rangle$ as follows (Liu and Liu, 2003):

$$ \langle f(\underline{r}) \rangle = \int_{\Omega} W_h(\underline{r} - \underline{r}') f(\underline{r}') d\underline{r}', \qquad (5) $$

where $W_h$ is weighted interpolation function of reaching length $h_W$, known as a smoothing kernel.

This integral procedure can then be discretized on numerical nodes. Applying this method to the governing equations

(Eqs. 1 and 2), the fluid is then discretized into so-called SPH particles. These particles move along with the flow, representing a Lagrangian computation node with an assigned mass. Equations (1) and (2) become, for a point $a$, a computational node in a Newtonian fluid the following system:

$$ \frac{\partial \rho_a}{\partial t} = \rho_a \sum_b ((\underline{u}_a - \underline{u}_b) \cdot \nabla W_{ab}) V_b + h_W c_0 \mathcal{D}_a, \qquad (6) $$

$$ \frac{\partial \underline{u}_a}{\partial t} = -\sum_b m_b \left( \frac{P_a + P_b}{\rho_b \rho_a} \right) \underline{\nabla}_a W_{ab} + \sum_b m_b \left( \frac{4v_0 \underline{r}_{ab} \cdot \underline{\nabla}_a W_{ab}}{(\rho_a + \rho_b)(\underline{r}_{ab}^2 + \eta^2)} \right) \underline{u}_a + \underline{g}, \qquad (7) $$

$$ \frac{\partial \underline{r}_a}{\partial t} = \underline{u}_a, \qquad (8) $$

where the subscripts $a$ and $b$ refer to the interpolating and neighbouring particles, respectively; $\underline{r}_i$ is the positional vector of particle $i$; $\underline{r}_{ij} = \underline{r}_j - \underline{r}_i$; $V_i$ is the volume occupied by the particle; $m_i$ is the mass assigned to the numerical particle; $W_{ij}$ is the smoothing kernel; $\eta = 0.01 h_W$, to avoid singularity; $\mathcal{D}_i$ is a density diffusion term (see Supplement); and $v_0$ is the kinematic viscosity of the fluid. Equation (7) is written for Newtonian fluids. By definition of how the solution is computed, there is an inherent “smoothing error” in the solution. The smoothing kernel used is the quintic Wendland (1995) kernel:

$$ W(r, h_W) = \alpha_D \left(1 - \frac{1}{2} \frac{r}{h_W}\right)^4 \left(2 \frac{r}{h_W} + 1\right) \qquad (9) $$
for $0 < \frac{r}{h_W} < 2$.

Here, $\alpha_D$ depends on the dimension of the model: $\alpha_{2D} = \frac{7}{4}\pi h_W^2$ and $\alpha_{3D} = \frac{21}{16}\pi h_W^3$. The value of the smoothing length $h_W$ is set through a coefficient relating $h_W$ and $d_p$: $h_W = C_h \sqrt{2} d_p$ (in 2D).

In this paper, the boundary conditioning is done using the mDBC (“modified Dynamic Boundary Condition”) method from English et al. (2021). Further information about the equations in the method used in this paper (time-stepping – Courant condition coefficient set to 0.2; boundary conditions; shifting method from Lind et al., 2012; and density diffusion methods) can be found in Sect. S1 in the Supplement.

### 2.2 Collision algorithm

In the DualSPHysics software, the CHRONO library (Anitescu and Tasora, 2008) is coupled with the SPH solver to deal with complex solid–fluid interactions (Martínez-Estévez et al., 2023). The rigid bodies are considered to be a subset of the SPH particle. The core concept of the CHRONO solver method is based on the differential variational inequality (DVI) approach (Pang and Stewart, 2008), which shares a mathematical framework with non-smooth contact
dynamics (Moreau, 1977; Jean and Moreau, 1992; Radjai and Richefeu, 2009). In a nutshell, the interaction between bodies are simultaneously computed for all contacts in a unified manner by considering the transitions between local and global/generalized frames. In addition to Newton's law relations involving impulses instead of forces, some conditions need to be considered to satisfy the non-penetration of the bodies (Signorini condition) and the Coulomb friction through inequalities. This model has been extensively used for different problems relating to rigid-body dynamics, especially granular material and complex mechanical systems such as rovers and wheeled and tracked vehicles. A review of validated test cases can be found in Tasora et al. (2016).

In the coupling, the forces arising on the fluid-driven object are computed and integrated in time as follows:

$$M \frac{dV}{dt} = \sum_k m_k \underline{f}_k, \\ I \frac{d\Omega}{dt} = \sum_k m_k ((\underline{r}_k - \underline{R}_0) \times \underline{f}_k) \cdot \underline{z}, \quad (10)$$

where the subscript *k* denotes an SPH boundary particle with a force per unit mass $f_k$, $\underline{V}$ is the linear velocity, *M* is the mass of the rigid body, $\underline{R}_0$ is the position of the centre of mass, *I* is the moment of inertia and $\Omega$ is the rotational velocity at a position $\underline{r}_k$.

The CHRONO solver performs sub-iterations within an SPH time step to adjust the rigid-body kinematics, and it gives an updated centre of mass with linear and angular velocities of $\underline{V}$ and $\Omega$, respectively, back to SPH solver. The velocity of each particle boundary in SPH is updated as follows:

$$\underline{u}_k = \underline{V} + \Omega \underline{z} \times (\underline{r}_k - \underline{R}_0), \quad (11)$$

and a system update for fluid and boundary particles is performed. For more information, the reader is directed to Martínez-Estévez et al. (2023).

## 2.3 Physical properties

The aim of this paper is to model viscous laminar debris flow surges with a Reynolds number close to the creeping flow limit ($Re \approx 0.1$). We consider 2D flows of surges over an inclined plate.

Surges are assumed to be composed of a non-uniform front followed by a uniform zone, where the free-surface elevation ($z_{max}$) becomes parallel to the bed: $\partial_x z_{max} = 0$, where $\underline{x}$ is the direction longitudinal to the flow.

Hunt (1994) proposed theoretical study of the velocity profiles in a viscous laminar front, applied to debris flow processes. Assuming a steady flow of a fully developed laminar flow of a Newtonian fluid down an incline in 2D, the Navier-Stokes equations yield in the uniform zone is as follows:

<table><caption>**Table 1.** Characteristics of the experiments: density ($\rho$), viscosity ($\mu$), slope angle ($\theta$), conveyor belt velocity ($u_b$) and Reynolds number ($Re$).</caption><thead><tr><th>Exp. ID</th><th>ρ<br>(kg m<sup>-3</sup>)</th><th>μ<br>(Pas)</th><th>θ<br>(°)</th><th>u<sub>b</sub><br>(mm s<sup>-1</sup>)</th><th>Re<br>(-)</th></tr></thead><tbody><tr><td>ID<sub>1</sub></td><td>1383</td><td>5.6</td><td>15.3</td><td>75</td><td>0.34</td></tr><tr><td>ID<sub>2</sub></td><td>1380</td><td>4.9</td><td>11.9</td><td>164</td><td>1.3</td></tr></tbody></table>

$$p = \rho g (z_{\max} - z) \cos \theta, \\ u_x = \frac{\rho g}{2\mu} (z_{\max}^2 - (z - z_{\max})^2) \sin \theta, \\ \overline{||u||} = \frac{\rho g \sin \theta}{3\mu} z_{\max}^2, \quad (12)$$

where *p* is the pressure, *ρ* is the density of the fluid, *g* is the gravitational force, *μ* is the viscosity, *u* is the velocity (with $\cdots$ denoting an average, the suffix $\cdots x$ denoting the *x* component of the velocity, $||\cdots||$ being the magnitude operator and $\cdots$ denoting a vector), *z* is the distance from the incline (in the normal direction), $z_{max}$ is the free-surface elevation and $\theta$ is the angle of the incline (taking into account the boundary condition $u(0) = 0$).

## 2.4 Selected experimental setup

Freydier et al. (2017) present experimental data precisely measuring viscous surges, i.e. flowing material close to the maximum depth rushing over a dry bed. Using advanced velocity measurement techniques comprising the image analysis of seeded, transparent fluids enlighten with a longitudinal, vertical laser sheet, these experiments captured both the macroscopic behaviour of a viscous flow front – monitoring properties of the flow, such as free-surface elevation and front velocity – and internal dynamics within the flow front – measuring velocity profiles. The experimental setup is composed of a conveyor belt tilted to a chosen angle, transparent side walls and a wall upstream of the conveyor belt forcing a flow front to steadily form. This experimental setup is thoroughly described in Freydier et al. (2017) and Chambon et al. (2009). Fluids used are transparent mixtures of glucose and water, allowing one to measure accurate velocity fields and free-surface profiles. Freydier (2017) showed that these flows have a viscosity independent of the strain and strain rate, varying with the concentration of glucose in the mixture. High-definition, high-velocity images at different locations of the flow and complete velocimetry within the surge provide a full description of both the free-surface shape and the velocity fields within the flow.

The characteristics of the experiments used in this work for numerical validation of the software are shown in Table 1.
## 2.5 Rheology of neutrally buoyant mixtures

One of the key questions of debris flow modelling is the accurate representation of the interactions between the fluid phase and the granular content. The presence of non-Brownian particles in the flow increases overall viscosity, thereby impacting the flow front velocity and flow height (Coussot and Ancy, 2022). Thus, accurately representing the rheological effects of the presence of solid particles in the mixture is essential to ensure the performance of DualSPHysics with respect to modelling debris flow events.

In this section, we study the macroscopic viscosity of a hard disc suspension in a 2D Poiseuille flow against known semiempirical equations. The macroscopic viscosity is studied as a function of the solid fraction (aerial fraction $\phi$).

The rheology of suspensions is a thoroughly investigated field. A review of the different approaches, challenges and rheological models can be found in Guazzelli and Pouliquen (2018). Theoretical approaches draw from the pioneering work by Einstein (1906), which relates the particle fraction $\phi_D$, where $D$ denotes the dimension of the problem, to the apparent viscosity linearly. This development of non-Brownian physics is only valid for very dilute regimes ($\phi_D < 0.05$) and was extended to any dimension $D$ by Brady (1983) such that

$$ \mu_{eq} = \mu_0 \left( 1 + \frac{D+2}{2} \phi_D \right), \quad (13) $$

where $\mu_{eq}$ is the macroscopic viscosity of the mixture and $\mu_0$ is the viscosity of the interstitial fluid.

At larger volumetric fractions, the effect of neighbouring particles becomes significant, and the coupled interactions are expected to induce a viscosity contribution of $O(\phi^2)$ or even $O(\phi^3)$ for concentrated suspension. However, the precise description of the rheology of such mixtures remains difficult (Guazzelli and Pouliquen, 2018). Numerous attempts at describing this problem for higher volumetric fractions have been made in 3D. They usually aim at recovering the Einstein equation at a low concentration while also describing the diverging viscosity at high volumetric fractions (Guazzelli and Pouliquen, 2018). Among them, the Krieger and Dougherty (1959) model is one of the most widely used semiempirical formulas to predict macroscopic viscosity from solid concentration. It describes the behaviour of the mixture as a power law with two fitting parameters, $\phi_m$ and $p$:

$$ \frac{\mu_{eq}}{\mu_0} = \left(1 - \frac{\phi}{\phi_m}\right)^{-p\phi_m}, \quad (14) $$

where $\phi_m$ is the maximal flowing fraction, with typical values of between 0.58 and 0.65 (Guazzelli and Pouliquen, 2018).

In order to recover the Einstein law for diluted suspension, the exponent $p$ is often considered to be equal to $p = 2.5$. This value has shown especially accurate results for diluted flows but has been shown to be overestimated for concentrated suspensions (Guazzelli and Pouliquen, 2018).

<figure><img src="image_3.png" /><figcaption>Figure 1. Display of the overall flow for simulated flow with SPH and free-surface elevation of the experimental measurements for ID1 with $d_p = 0.001$ m and $C_h = 1.0$, $x^* = x/h_{max}$, and $z^* = z/h_{max}$, where $h_{max}$ is the theoretical maximal flow height in the uniform zone.</figcaption></figure>

# 3 Results

## 3.1 Viscous Newtonian surges

In this section, the goal is to investigate the macroscopic and internal features of a surge of a viscous Newtonian fluid and compare our findings with the experiments of Freydier et al. (2017) presented in Sect. 2.4

### 3.1.1 Numerical setup

The numerical setup reproduces the experimental setup as a 2D flow front (Fig. 1), with the flow features shown in Table 1. A surface of simulated fluid is generated onto an inclined conveyor belt. The inclination of the plane is represented by an inclined gravity vector in order to have a horizontal flow, following the slope $\theta$. The fluid is contained by an upstream wall. The surface of the experiment has been set up so that the wall is further than the distance required to reach a uniform height, as shown in Freydier et al. (2017). The flow height is initialized at the theoretical height of the uniform zone (from Eq. 12). At the initial time step, the flow is set up to be a rectangle of length $L = 0.27$ m and height $h = 0.017$ m for experiment ID1 and of $L = 0.4$ m and $h = 0.03$ m for experiment ID2. The conveyor belt is simulated by a moving bed, moving in the upstream direction, and is implemented in a longitudinally periodic domain, as if infinite, to avoid any discontinuity with respect to the boundary. It initiates a backwards motion, using a ramp to reach the velocity $u_b$. This ramp was put in place in order to avoid instabilities caused by the sudden infinite acceleration of an instantaneous motion.

A convergence study was performed on the results. For each case, five values of kernel coefficients $C_h$ were tested (tested values for $C_h$: 1, 1.3, 1.5, 1.8 and 2.2). Kernel coefficients are a measure of the ratio between the smoothing length and particle spacing ($h_W = C_h \sqrt{2d_p}$ in 2D). Convergence with respect to the particle spacing distance $d_p$ was investigated for each value of $C_h$ (tested values for $d_p$: 0.001, 0.005 and 0.00025 m) following (Quinlan et al., 2006). The initial particle spacing was chosen such that $h/d_p > 15$,
i.e. at least 15 SPH particles describe the flow in the uniform zone for any of the simulations.

The validation of the behaviour of the flow will be regarded as the accuracy to represent both macroscopic and internal features. The bulk flow behaviour is studied through (i) the velocity of the flow front, (ii) the free-surface elevation in the uniform zone and (iii) the free-surface shape. The internal kinematics of the flow are studied through (iv) the velocity profiles at three different locations in the flow front (Fig. 3a) that correspond to the positions studied in Chambon et al. (2019).

### 3.1.2 Data processing

The numerical model is run for a sufficient period of time, typically more than 30 s, ensuring that steady state has been reached (see Figs. S2 and S3 in the Supplement). With SPH, the velocity of the continuum can only be estimated through the velocity of the particles in the discretized flow. With the movement of particles in the flow intrinsic to the SPH method, instantaneous measurements can be affected by fluctuations in the positions of the particles within the sampling window. To avoid these instantaneous effects, the results are averaged over 5 s.

For all results, heights and lengths are normalized by $h_{max}$, the theoretical flow height in the uniform zone, so that $h^* = h/h_{max}$ and $x^* = x/h_{max}$. The position in $x$ is defined relatively to the tip of the front. This theoretical flow height is derived from Eq. (12). The velocities are normalized by the velocity of the conveyor belt $u_b$ so that $u^* = u/u_b$.

### 3.1.3 Bulk features of the flow

In Fig. 2, errors in the macroscopic behaviours are plotted. Figure 2a shows the error in the velocity of the flow front, which is defined as follows:

$$ \zeta_f = \frac{|u_f|}{u_b} = |u_f^*|, \qquad (15) $$

where $u_f$ is the velocity of the front provided by the numerical model averaged over 1 s, taken as the derivative of the location of the toe of the front. As a reminder, because the setup is on a conveyor belt going at $u_b, u_f \to 0$.

Figure 2b shows the relative error in the flow height in the uniform zone for which $\partial_x h = 0$, which is defined as follows:

$$ \zeta_h = \frac{|h_u - h_{max}|}{h_{max}} = |h_u^* - 1|, \qquad (16) $$

where $h_u^*$ is the normalized height in the uniform zone, taken in the uniform zone with $x^* > 9$.

Figure 2c shows error in the overall free-surface profile as a root-mean-square error in the free-surface elevation:

$$ \zeta_{fs} = \sqrt{\frac{1}{N_x} \sum_{i=0}^{N_x} \left( \frac{h^*(x_i) - h_{ref}^*(x_i)}{h_{ref}^*(x_i)} \right)^2}, \qquad (17) $$

where $h^*(x_i)$ is the normalized free-surface elevation of the simulated flow at position $x_i$, $h_{ref}^*(x)$ is its analogous normalized free-surface elevation provided by the experimental data and $N_x$ is the number of points along $x$ on the front shape measurements provided by Freydier (2017).

Both experiments ID1 and ID2 are tested (see Table 1). In Fig. 2, full symbols show the convergence study for experiment ID1, whereas void symbols show the convergence study for experiment ID2. Overall, the convergence is shown for both experiments and for the three features of interest. The numerical solution converges to less than 1% of the theoretical value for $u_f$ (see Fig. 2a), less than 5% for $h_{max}$ (see Fig. 2b) and less than 10% on the overall shape of the flow front (see Fig. 2c). The smoothing coefficient has a drastic influence on the prediction of the model: $C_h$ drives not only the numerical computation but also the threshold of the method to detect the free surface (Domínguez et al., 2021). In applications where the free surface is very variable, as in this study, the value of the smoothing coefficient has a high impact on the overall performance of the model. A systematic study of the influence of such a coefficient is a relevant method to rule out any error from this choice.

### 3.1.4 Internal dynamics of the flows

Comparison between velocity profiles given by the numerical model and the profiles given by Freydier et al. (2017) at three different locations are studied for experiment ID1. The same features were monitored for experiment ID2 with similar performance (see Fig. S1 in the Supplement). In Fig. 3a, the three locations are highlighted, with the exact value of the positions shown in Table 1.

In Fig. 3b–d, the velocity profiles for $C_h = 1.8$ are plotted for each location and resolution. The overall behaviour of the flow is very well reproduced for all resolutions. Profiles are well within the error bars of the experimental data. Discrepancies in the height of the flow for $x_f^* = 0.76$ are due to the experimental measurements: close to the free surface, access to the velocity of the flow becomes difficult (Freydier, 2017). Even near the front, the shape of the velocity profile accurately reproduces the dynamics of the surge.

A global criterion to characterize the convergence for each location $i$ is defined as the root-mean-square error in the velocity profile along $x$:

$$ \zeta_u = \sqrt{\frac{1}{N_h} \sum_{j=0}^{N_h} (u^*(h_j) - u_{ref}^*(h_j))^2}, \qquad (18) $$

where $u^*(h_j)$ is the normalized longitudinal velocity of the simulated flow at elevation $h_j$ in the flow, $u_{ref}^*(h_j)$ is its analogous normalized longitudinal velocity provided by the experimental data and $N_h$ is the number of points in the vertical direction (which depends on the position $i$).

In Fig. 3e–g, the error $\zeta_{u,j}$ is plotted for each ($d_p, C_h$). Convergence is shown for each smoothing coefficient. The
<figure><img src="image_4.png" /><figcaption>Figure 2. Convergence study on macroscopic features of the flow, $d_p$ (in m): (a) error in the flow front velocity, as defined in Eq. (15); (b) error in the free-surface elevation in the uniform zone, as defined in Eq. (16); and (c) root-mean-square error in the surface shape for different resolutions and smoothing coefficients, as defined in Eq. (17). Full symbols represent the results of experiment ID₁, while void symbols represent the results of experiment ID₂.</figcaption></figure>

convergence is similar regardless of the choice of a smoothing coefficient $C_h$, and overall the error becomes smaller than 5% when using the finest resolution. Considering the precision of the experimental measurements, this is very satisfactory.

To have a good agreement between both the accuracy of macroscopic and internal features of the flow and computational time, the smoothing coefficient $C_h = 1.8$ is chosen for use in the rest of the study.

## 3.2 Macroscopic viscosity of mixtures

Following the validation of the highly viscous fluid in Dual-SPHysics, there is a missing piece to develop a more accurate representation of a 2D debris flow model – namely, the macroscopic viscosity changes due to the presence of solid particles in the flow. This needs to be correctly restituted by the model so that we can explore macroscopic behaviours in the debris flow model.

To test for the behaviour of a mixture, a 2D loaded Poiseuille test is performed. A fluid is generated between two horizontal plates, separated by a distance $B$, to which a pressure gradient is applied along the longitudinal direction. The pressure gradient is forced using a modified gravity in the $x$ direction. The geometry is described in Fig. 4, where the boundaries are periodic along the $x$ direction. The test is done using different concentrations of monodisperse discs of radius $r$ randomly positioned in the channel. Each object has a no-slip condition with the fluid and is neutrally buoyant. The overall macroscopic viscosity is retrieved from the maximal velocity (Guyon et al., 2012):

$$ \mu_{eq} = \frac{B^2}{8\mu_{x,max}} \left| \frac{dp}{dx} \right| \qquad (19) $$

The maximal velocity along $x$ ($u_{x,max}$) is taken as the average of the velocities along $x$ of the fluid particles with posi-

tions at the centre of the flume ±1 cm between 0.8 and 1 s, once stationarity is reached (see Fig. S4). The experimental designs are summarized in Table 2.

In Fig. 5, the comparison between the numerical results of the macroscopic viscosity against solid concentration and the empirical model is plotted. Equation (14) is fitted for both $p$ and $\phi_m$. The black lines represent the uncertainty due to the spatial averaging.

The overall shape of the results does match the typical shape of this predictive model, and the results can be fitted with precision. The root-mean-square error of the fit is evaluated as follows:

$$ \text{RMSE} = \sqrt{\frac{1}{N_s} \sum_{i=0}^{N_s} \left( \frac{\mu_{eq}}{\mu_0} - \frac{\mu_{est}}{\mu_0} \right)^2} \qquad (20) $$

where $N_s$ is the number of tested concentrations and $\mu_{est}$ is the estimated equivalent viscosity using the parameters of the fit.

The macroscopic response to the concentration of solid particles is satisfactory: the power-law shape of the classical models is well reproduced. Specifically, when $\phi$ reaches $\phi_m$, there is a divergence in the viscosity. This divergence is a key element of the flow of the mixture of grains and fluids. The fit of Eq. (14) leads to a root-mean-square error RMSE = 5.07%. When solving a full-scale debris flow, we can expect reasonable behaviour of the overall viscosity change when incorporating the boulders.

## 3.3 Towards full-scale debris flow modelling

Section 3.1 and 3.2 show that the method yields accurate results for both mesoscopic and microscopic mechanics within a surge in 2D. This last section explores the flow of a 2D debris flow model at the surge scale with three different solid concentrations.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Study of the internal dynamics of the flow: (a) positions of the velocity profiles of interest for three different locations in the front ($x_f^* = 0.76$ and $x_f^* = 1.5$ is close to the front tip, where $\partial_x z$ is significant and $x_f^* = 9.7$ in the uniform zone); (b–d) velocity profiles of the numerical experiment plotted onto the velocity profile of the experimental data by Freydier et al. (2017) for $C_h = 1.8$ (convergence in $d_p$ is shown); (e–g) root-mean-square error, as defined in Eq. (18), on the velocity profiles at the three positions showing the convergence study of the model in both $d_p$ and $C_h$ (with $d_p$ given in metres).</figcaption></figure>

<table><caption>**Table 2.** Summary of the designs of the numerical experiments; note that for concentrations above 0.4, $d_p/r = 4.5\%$.</caption><thead><tr><th>ρ<br>(kg m<sup>-3</sup>)</th><th>||g|| = 1/ρ ||∇P||</th><th>μ<sub>0</sub><br>(m<sup>2</sup> s<sup>-1</sup>)</th><th>L<br>(m)</th><th>B<br>(m)</th><th>φ<br>(-)</th><th>r<sub>max</sub>/L<br>(-)</th><th>d<sub>p</sub>/r<br>(-)</th><th>C<sub>h</sub><br>(-)</th><th>Collision algorithm</th></tr></thead><tbody><tr><td>1000</td><td>0.8</td><td>0.467</td><td>1.0</td><td>0.6</td><td>[0.05, 0.5]</td><td>2%</td><td>8.5% or 4.5%</td><td>1.8</td><td>CHRONO</td></tr></tbody></table>

### 3.3.1 Geometry of the setup

Thanks to the validation of the fluid phase for complex free surfaces and the validation of the macroscopic viscosity of the mixtures, a preliminary study of numerical experiments on debris flows in 2D is carried out.

Here, the setup is highly simplified to investigate the influence of boulders on the macroscopic behaviour of the flow. The numerical flow is solved in 2D using a conveyor belt

channel which is similar to that in Sect. 3.1 but to which boulders are added. The slope of the channel is taken to be similar to the slopes of debris flow torrents, $S = 10.1°$.

The initial conditions are based on the field data in Lapillonne et al. (2023) and Piton et al. (2018). In the following analysis, the goal is not to exactly represent a given event; rather, we aim to use sensible orders of magnitude for the flow of interest.
<figure><img src="image_6.png" /></figure>

**Figure 4.** Experimental setup for the loaded Poiseuille case; the pressure gradient is introduced through $\mathcal{S}$, and the colour bar represents velocity in the longitudinal direction.

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

**Figure 5.** Evolution of the macroscopic viscosity with solid concentration, with the numerical data and their error bars.

The flow is initialized with a constant flow depth $h_{\text{surge}} = 1.4$ m; $u_b$ is then chosen so that the Froude number of the initial conditions reaches $Fr \approx 1$. The choice of this Froude number is motivated by two factors: (i) the flow must represent debris flow Froude numbers that are found in nature (see Lapillonne et al., 2023); (ii) the introduction of boulders in the flow is expected to slow down the flow, thereby lowering its Froude number, compared to a grain-free surge. This way, we ensure that the flow characteristics will remain in reasonable ranges compared to values of the Froude number measured in the field ($Fr < 3$, with a distribution centred around $Fr = 1$). The viscosity of the interstitial fluid is kept the same throughout all of the simulations. Its value is chosen from the analytical derivation in Eq. (12) so that the grain-free flow would exhibit the features explained above ($Fr = 1$). These values are summarized in Table 3

The length of the numerical setup $L_x$ represents the front of a debris flow surge length. To determine it, the analysis was based on the volumes of the events on the upper station of the Réal torrent. The median volume of an event is $V \approx 1000 \text{ m}^3 \text{ s}^{-1}$, with a channel width of $B = 8$ m. Using a simplistic rectangular approximation, this leads to a length of the flow front and body of a corresponding event $\frac{V}{h_{\text{max}}B} \approx 90$ m. As the setup aims at only representing the front of the flow, the actual length of the numerical setup is $L_x = 33$ m.

Each boulder is represented as a disc, the size of which is randomly sampled on a grain size distribution measured both on the front and the body of an actual debris flow (see event description in Piton et al., 2018). Only grain sizes bigger than $d_{\text{truncature}} = 20$ cm and up to $h_{\text{max}}$ are sampled. Indeed, computational time is driven partly by the number of individual grains represented, the number of collisions and the size of the pores between the grains. Thus, modelling all of the singular grains in the debris flow material is not computationally reasonable. This computational barrier creates a conundrum because those small grains not only influence the overall viscosity of the interstitial fluid in which we encapsulate them, as demonstrated in the previous section, but they also affect the large boulders through granular dynamics, whereby small particles are driven downward and large particles migrate upward through kinetic sieving and squeeze expulsion (Gray, 2018). This phenomena is key to understanding how huge boulders with a density of $\approx 2600 \text{ kg m}^{-3}$ seem to float in debris flows, the average density of which is $\approx 1800–2000 \text{ kg m}^{-3}$. To mimic this effect that cannot be explicitly modelled due to the lack of representation of smaller grain sizes, the boulder relative density $\frac{\rho_b}{\rho}$, where $\rho_b$ is the boulder density, was arbitrarily set to 0.9. This crude assumption certainly deserves further investigation in future work. To place each boulder, their positions are drawn iteratively and tested to avoid any intersection with previously placed boulders. In the $z$ direction, their position is randomly drawn on a uniform probability law. In the $x$ direction, they are placed following a truncated lognormal law centred close to the front, as schematized in Fig. 6b. This means that most boulders will be generated near the front. This approach is a makeshift technique to preferentially place boulders near the front, as observed in the field (Johnson et al., 2012). The flow then reorganizes their position during the simulation.

The data of interest are measured over 55 s after the first 20 s with a sampling period of 0.5 s. The first 20 s of the flow are discarded because the flow has not yet reached stationarity (see Fig. S5 in the Supplement). The timescale chosen to do the observation is of the order of magnitude of the time travel between sensors in actual field monitoring stations. The data are averaged over that duration to have a representative measurement of the flow rather than an instantaneous one. Three experiments with an increasing number of boulders are carried out and summarized in Table 3.
<table><caption>**Table 3.** Characteristics of the 2D debris flow experiment.</caption><thead><tr><th>Number of boulders (-)</th><th>φ (-)</th><th>Slope (°)</th><th>μ<sub>0</sub> (m<sup>2</sup> s<sup>-1</sup>)</th><th>ρ (kg m<sup>-3</sup>)</th><th>Initial flow height (m)</th><th>u<sub>b</sub> (m s<sup>-1</sup>)</th><th>d<sub>p</sub></th><th>C<sub>h</sub></th></tr></thead><tbody><tr><td>[0, 50, 100, 154]</td><td>[0,0.13,0.29,0.53]</td><td>10.1</td><td>0.280</td><td>2000</td><td>1.4</td><td>3.5</td><td>0.01</td><td>1.8</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Debris flow geometry generation: (a) grain size distribution measured on the debris flow event described in Piton et al. (2018); (b) example of the distribution of the boulders along *x* and *z*, with the lognormal law of distribution for *x* illustrated in green (right axis) and *z* being random with respect to flow depth.</figcaption></figure>

The value of $h_{max}$ is taken as the maximal value of the flow height in the surge at each time step. The free-surface elevation is instantaneously very dependent on the local presence of the boulders: boulders clustering and piling up will temporarily and locally increase $h_{max}$. The statistical variability in $h_{max}$ was recorded to illustrate how the presence of boulders might introduce a kind of noise in this regularly monitored parameter. The upstream first 5 m section is ignored to evaluate maximal flow height in order to avoid eventual effects of the upstream wall. The Froude number and the equivalent viscosity are computed as follows:

$$ 
\begin{aligned}
Fr &= \frac{\bar{u}}{\sqrt{g h_{max}}}, \\
\mu_{eq} &= \frac{\rho g \sin\theta}{3 \bar{u}} h_{max}^2.
\end{aligned} 
\qquad (21) $$

### 3.3.2 Effect of the introduction of boulders on the macroscopic behaviour of the flow

In Fig. 7c–f, the overall view of the surges is shown. With increasing concentration, the shape of the surge shows higher flow depth, here being a proxy for increasing apparent viscosity. Overall, the surges show a recurrent circulation of the grains throughout the simulation. The boulders stay within the flow, and the simulations successfully behave as debris flow surges.

In Fig. 7a, the Froude numbers of the flow are shown against the boulder concentration in box plot form. The whiskers represent the quartiles of the sampling of $h_{max}$ over time. The Froude numbers decrease with increasing boulder concentration, as the obstruction to the flow increases the ef-

fective viscosity of the mixture. In this figure, the whiskers of the box plot represent the quartiles due to the time-averaging of the data. Throughout the simulation, the Froude number varies marginally, as the interquartile range of the box plot is narrow and stays > 0.5, i.e. at orders of magnitude consistent with the field data of Lapillonne et al. (2023).

In Fig. 7b, the viscosities are shown to increase with increasing boulder concentration. The relationship between viscosity and boulder concentration shows a positive nonlinear trend: as the concentration of grains increases, the viscosity of the mixture exhibits a more-than-proportional growth. Both $h_{max}$ and $\bar{u}$ are used to compute $Fr$ and $\mu$, but the latter shows stronger instantaneous variations (larger box plot) than the former, due to the difference in non linearities regarding $h_{max}$ in Eq. (21). These results demonstrate that complex coupled methods are a relevant tool to explore debris flow motion and internal dynamics, as well as the complementary contribution of the interstitial fluid and granular content to the overall behaviour of the surges.

## 4 Discussion

This section mirrors the organization of Sect. 3 so that each subsection discusses the corresponding results from the respective subsections of Sect. 3.

### 4.1 Viscous Newtonian surges

The accuracy of the results at different points in the flow are promising in terms of possible applications, especially close to the toe of the front. Precise knowledge about flow front
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The 2D debris flow model: box plot of (a) the Froude numbers and (b) the ratio of equivalent viscosity over fluid viscosity, where whiskers represent the upper and lower quartiles, the orange bar represents the mean value, and the vertical grey line in panel (b) highlights $\mu/\mu_0 = 1$. Velocity fields and a general view of the different simulations for (c) $\phi = 0.53$, (d) $\phi = 0.29$, (e) $\phi = 0.13$ and (f) $\phi = 0$. The white arrow represents the movement of the conveyor belt. Views have been tilted back to the horizontal plane in order to ease the reading of the comparison; the tilted gravity vector is shown.</figcaption></figure>

characteristics is of major interest for applications on mud-flow behaviour as well as debris flow modelling.

Overall, the method is thoroughly validated for viscous surges. Both the macroscopic behaviour and internal dynamics are recovered to a satisfactory accuracy. The errors in the simulations go well into the uncertainty brackets of the experimental measurements, meaning that the validation cannot be further optimized. This is a first step towards the modelling of boulder-laden debris flows with the SPH method.

Within the scope of debris flow modelling, studying the internal dynamics of these complex flow fronts will be a significant advantage to understand their global dynamics, including their flowing and stopping conditions. To do so, the validation of the internal dynamics of such a laminar viscous flow front is a necessary step to approach a reasonable rendition of the internal physics. Validating some of the driving processes of the internal dynamics and simple indicators of the macroscopic behaviour of the flow is a step forwards in more realistically modelling actual debris flows. Key elements like front velocity, flow height and front shape were correctly represented in this case study. These parameters can be measured in the field (Hürlimann et al., 2019; Lapillonne et al., 2023; Schöffl et al., 2023; Aaron et al., 2023). The most advanced monitoring stations are even equipped to measure surface velocities (Theule et al., 2017; Schöffl et al., 2023), basal forces (Nagl et al., 2020) or the 3D shape of the debris surges (Aaron et al., 2023). In further investigation, sophisticated methods of monitoring and of numerical modelling will be used in conjunction to improve our knowledge of debris flow dynamics.

Within this work, we validated this numerical approach for such applications. Knowledge on experimental data and

data uncertainty for such mudflows and debris flows in the field leads us to consider the errors that we highlighted to be reasonable and acceptable: no model will ever be a perfect representation of such complex processes.

## 4.2 Macroscopic viscosity of mixtures

While the overall shape of the rheological response is correctly retrieved in the Poiseuille test, the actual values of the coefficient *p* in the different fitting of empirical equations do not correspond to the usual values found in the literature. Generally, the classical fitting parameters of Krieger-Dougherty are known to not always perfectly fit the experimental observations but, rather, to broadly represent the behaviour of mixtures (Guazzelli and Pouliquen, 2018). However, we acknowledge that some numerical studies, such as Kromkamp et al. (2006), are able to retrieve accurate 2D descriptors of the flow.

The lubrication forces in SPH are highly dependent on resolution. When two solid particles interact through the fluid at close range, the lubrication forces should have a great influence on their movement. However, in SPH, these lubrication forces diverge with a decreasing number of particles. If the solid particles become too close to one another, the resolution becomes insufficient at a certain scale. To avoid such problems, some numerical methods implement an artificial additional lubrication force (Bian and Ellero, 2014; Chèvremont et al., 2020). However, in the scope of debris flow modelling, where collisions drive the macroscopic flow, this lack of precision is taken as a model uncertainty.

The preference for a Poiseuille setup over a classical Couette flow is not common because of the risk of drifting of
the solid particle and irregular solid concentration. Irregular solid distribution due to migrations of the solid elements in a Poiseuille flume is a well-known phenomenon that would lead to an incorrect estimation of the actual aerial fraction (Guazzelli and Pouliquen, 2018). In the case presented in this paper, the timescale of this rearrangement is very long compared to the timescale of stabilization of a velocity profile. In that sense, making the assumption that no lateral migration occurs is sensible. Empirically, no migration is observed during the measurement period.

Overall, for the application of this paper, we consider that the divergence of the viscosity at increasing concentration, as well as the overall shape of the relationship between viscosity and concentration, is sufficient.

## 4.3 Debris flow model

Previous pioneering work has modelled mixtures of interstitial fluid and grains to approach debris flow behaviour, although always relying on simplifying hypotheses: Canelas (2015) employed a small-scale model of boulder-laden flows with the fluid as water; Leonardi et al. (2014) modelled a column collapse experiment involving mud and gravel; and Leonardi et al. (2015) modelled the impact of mixed grains and Bingham interstitial fluid against a flexible barrier, focusing on the impact dynamics. Our present work presents a proof of concept for a complex numerical model that encompasses multiple scales while also being accessible for specific applications.

As explained above, the main features of the viscosity changes due to the presence of grains are reproduced by the model, i.e. its non-linearity and divergence at high concentrations (Fig. 5). For field applications, the complicated mixture of water, sand, gravel and eventual clay is assumed, in a first, crude approximation, to be represented by our viscous flow. Its viscosity is, thus, seen as a calibration criterion that encapsulates the more complex behaviour of the mixture, the precise composition of which is variable and unknown. The capacity of the tool to model grains is only used here to capture the influence of large boulders on the overall flow behaviour (Fig. 7). Eventually, the capacity to model larger surges, with more grains of increasingly finer diameter, is limited by the computational power.

The lack of representation of cohesive rheology via a non-Newtonian constitutive model at the SPH resolution leads to an underestimation of the plug. In Fig. 7c–f, the effect of boulders on the velocity profile can be seen, with flows with a higher boulder concentration (Fig. 7c–e) having a larger section with an almost constant velocity. This effectively leads to a pseudo-plug that is driven by the presence of granular matter. However, the viscous matrix should also drive a plug zone. This plug influences the shear profile in the flow and affects the relative motion of the boulders in the flowing material, impacting both boulder mobility and spatial distribution. This assumption, combined with the 2D assumption, are

acknowledged as intrinsic limitations of the model. Nonetheless, the objective of this paper is to capture macroscopic features comparable to those observed in the field for which the model remains satisfactory.

The 2D experiments that we present are another proof of concept showing that coupled solid-fluid numerical models are able to qualitatively reproduce the macroscopic flow features of a field-scale debris flow. However, the 2D assumption for debris flows is unrealistic. The 3D motion of boulders likely drives a lot of the efforts in the granular skeleton. Such 2D models might be somewhat comparable to a very confined modelling of the flow, i.e. the boulders not being able to move laterally. This confinement probably leads to an overestimation of the slowing down of the flow and also likely tends to under-deposit grains, especially as it cannot represent the lateral deposition and formation of levees. Moreover, there is a diminution of lubrication effects because of the 2D simplification, as the lateral rearrangement of grains are not replicated. However, we believe that the extension of this model to a 3D case is perfectly feasible: the validation of each phase would be exactly the same. Nevertheless, it is significantly more resource-intensive to extend such models to 3D, and this was considered outside the scope of this study.

A more complex shape of the boulders, with elements that cannot roll, would be another step toward a more realistic approach. This could be generated as an ensemble of clumped, spherical elements. The interactions of these granular elements with the flat bottom would likely have an even stronger non-linear impact on the macroscopic flow. However, this first approach with simplified elements (and thus more accessible geometry) already shows the dramatic influence of boulders on the flow.

This study would benefit from a complete statistical analysis with randomly generated initial conditions and repeated occurrences. Indeed, the current result could be influenced by the initial grain distribution (size and location) sampled for the case. A statistical analysis running this experiment with an appropriate number of initial setups for the three concentrations would allow us to study more precise patterns, such as the recirculation of boulders into the flow. However, this does not influence the validity of the model, even with only one simplistic case.

The lubrication effects are under-represented, not only due to resolution and dimensional effects but also because of the nature of the interstitial fluid used in the flow. Indeed, the artificial interstitial fluid aggregates the watery composition of the overall mixture and the grain sizes ranging from clay to cobble ($d_{truncation}$ = 20 cm). This interstitial fluid should have a strong effect on the overall macroscopic rheology. The inclusion of clay and other granular material argues for the representation of the interstitial fluid as non-Newtonian with a yield criterion (Ancey, 2007). However, this simplified model already encompasses parts of the overall behaviour by incorporating granular elements and shows satisfactory features for debris flow modelling applications.
# 5 Concluding remarks

This paper shows that a solid-fluid model using Dual-SPHysics is reliable for debris flow modelling. The work addressed different steps that feed towards an accurate modelling of such complex flows:

- Viscous surges are validated against precise experimental data. Remarkable agreement between the simulations and the experimental data is found, with errors in the macroscopic features of < 5% after convergence. The velocity profiles at three different positions in the flow are also reproduced with errors of < 5% after convergence. It is, to the best of our knowledge, one of the few validations of the SPH approach for very viscous flows near the creeping threshold (Re ≈ 0.1).

- Viscous flows loaded with small grains are investigated against the Krieger-Dougherty relationship. The overall shape of the non-linear relationship between the solid concentration and macroscopic viscosity is fitted, with a root-mean-square error of 5.07%. The divergence between the expected values of the fitted and reference coefficients driving the relationship is assumed to be due to (i) the scarcity of available coefficients for 2D numerical models and (ii) a divergence in the representation of the lubrication forces at small scales within the method. The method is, however, still considered valid in the scope of debris flow investigation.

- Two-dimensional simplified debris flow surges with polydisperse boulders are modelled at various concentrations. Field-like behaviours were retrieved with Froude numbers in the range of the field measurements. Macroscopic viscosity is found to evolve non-linearly with the introduction of boulders. The importance of using complex models is highlighted, and the model is considered promising for debris flow physics exploration. There is no argument against the extension of this model to a 3D case, which makes its potentially powerful.

To conclude, the different elements of the model are considered validated for the use of Newtonian debris flow surges. This study specifically highlights that the presence of boulders in the flow brings a non-linearity to the viscosity that should be considered when using fully continuous fluid-like models. This model, however, is only a first stepping stone towards more complex, complete models. Specifically, encompassing the effects of non-spherical elements (or of a non-Newtonian interstitial fluid) would be of great interest to better explore debris flow physics. Nonetheless, this simplified model is able to investigate questions of interest to the debris flow community, such as the study of the dead zone when impacting a barrier or the co-interaction of the granular and flow fronts.

*Code and data availability.* The current version of the model is available from the project website (https://dual.sphysics.org/, last access: 20 January 2024; GNU Lesser General Public License), and the code and data are available from https://doi.org/10.57745/VCTHIS (Lapillonne, 2024).

*Supplement.* The supplement related to this article is available online at https://doi.org/10.5194/gmd-18-7059-2025-supplement.

*Author contributions.* Conceptualization: SL, GP and VR; data curation: GC; formal analysis: SL and GF; funding acquisition: VR and GP; investigation: SL, GP and GF; methodology: all authors; project administration: GP; software: SL and GF; supervision: GP, VR and GC; visualization: SL; writing – original draft: SL; writing – review and editing: all authors.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Financial support.* The research of Suzanne Lapillonne, Vincent Richefeu and Guillaume Piton has been supported by the LabEx Tec21 Investissements d'avenir (grant no. ANR-11-LABX-0030).

*Review statement.* This paper was edited by Lele Shu and reviewed by two anonymous referees.

# References

Aaron, J., Spielmann, R., McArdell, B. W., and Graf, C.: High frequency 3D LiDAR measurements of a debris flow: a novel method to investigate the dynamics of full scale events in the field, Geophys. Res. Lett., 50, e2022GL102373, https://doi.org/10.1029/2022GL102373, 2023.

Albaba, A., Lambert, S., Nicot, F., and Chareyre, B.: Modeling the impact of granular flow against an obstacle, in: Recent Advances in Modeling Landslides and Debris Flows, edited by: Wu, W., Springer International Publishing, 95–105, https://doi.org/10.1007/978-3-319-11053-0_9, 2015.

Ancey, C.: Role of particle network in concentrated mud suspensions, in: Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, Millpress Science Publishers, 257–268, 2003.
Ancey, C.: Plasticity and geophysical flows: a review, Jour-
nal of Non-Newtonian Fluid Mechanics, 142, 4–35,
https://doi.org/10.1016/j.jnnfm.2006.05.005, 2007.

Anitescu, M. and Tasora, A.: An iterative approach for cone com-
plementarity problems for nonsmooth dynamics, Computational
Optimization and Applications, 47, 207–235, 2008.

Bardou, E.: Méthodologie de diagnostic des laves torren-
tielles sur un bassin versant alpin, PhD thesis, EPFL,
https://doi.org/10.5075/epfl-thesis-2479, 2002.

Bel, C.: Analysis of debris-flow occurrence in active catchments of the French Alps using monitoring stations, PhD thesis, Université Grenoble Alpes, IRSTEA-ETNA, 2017.

Bel, C., Liébault, F., Navratil, O., Eckert, N., Bellot, H., Fontaine, F., and Laigle, D.: Rainfall control of debris-flow triggering in the Réal Torrent, Southern French Prealps, Geomorphology, 291, 17–32, 2017.

Bernard, M., Barbini, M., Berti, M., Boreggio, M., Simoni, A., and Gregoretti, C.: Rainfall-runoff modeling in rocky headwater catchments for the prediction of debris flow occurrence, Water Resour. Res., 61, e2023WR036887, https://doi.org/10.1029/2023WR036887 2025.

Berti, M., Martina, M., Franceschini, S., Pignone, S., Simoni, A., and Pizziolo, M.: Probabilistic rainfall thresholds for landslide occurrence using a Bayesian approach, J. Geophys. Res.-Earth, 117, https://doi.org/10.1029/2012JF002367, 2012.

Bian, X. and Ellero, M.: A splitting integration scheme for the SPH simulation of concentrated particle suspensions, Comput. Phys. Commun., 185, 53–62, 2014.

Brady, J. F.: The Einstein viscosity correction in n dimensions, Int. J. Multiphas. Flow, 10, 113–114, https://doi.org/10.1016/0301- 9322(83)90064-2, 1983.

Canelas, R.: Numerical modeling of fully coupled solid-fluid flows, PhD thesis, University of Lisboa, 2015.

Canelas, R. B., Domínguez, J. M., Crespo, A. J. C., Gómez-
Gesteira, M., and Ferreira, R. M. L.: Resolved simulation of
a granular-fluid flow with a coupled SPH-DCDEM model, J.
Hydrol., 143, 06017012, https://doi.org/10.1061/(asce)hy.1943-
7900.0001331, 2017.

Ceccato, F., Redaelli, I., di Prisco, C., and Simonini, P.: Im-
pact forces of granular flows on rigid structures: compar-
ison between discontinuous (DEM) and continuous (MPM)
numerical approaches, Comput. Geotech., 103, 201–217,
https://doi.org/10.1016/j.compgeo.2018.07.014, 2018.

Chambon, G., Ghemmour, A., and Laigle, D.: Gravity-driven surges of a viscoplastic fluid: an experimental study, Journal of Non-Newtonian Fluid Mechanics, 158, 54–62, https://doi.org/10.1016/j.jnnfm.2008.08.006, 2009.

Chambon, G., Bouvarel, R., Laigle, D., and Naaim, M.: Numerical simulations of granular free-surface flows using smoothed particle hydrodynamics, Journal of Non-Newtonian Fluid Mechanics, 166, 698–712, https://doi.org/10.1016/j.jnnfm.2011.03.007, 2011.

Chambon, G., Freydier, P., Naaim, M., and Vila, J.-P.: Asymp-
totic expansion of the velocity field within the front of vis-
coplastic surges: comparison with experiments, J. Fluid Mech.,
https://doi.org/10.1017/jfm.2019.943, 2019.

Chèvremont, W., Bodiguel, H., and Chareyre, B.: Lubricated con-
tact model for numerical simulations of suspensions, Powder
Technology, 372, 600–610, 2020.

Christen, M., Kowalski, J., and Bartelt, P.: RAMMS: Numerical simulation of dense snow avalanches in three-dimensional terrain, Cold Reg. Sci. Technol., 63, 1–14, 2010.

Coussot, P. and Ancey, C.: Rhéophysique des pâtes et des suspen-
sions, in: Rhéophysique des pâtes et des suspensions, EDP Sci-
ences, 2022.

Coussot, P. and Meunier, M.: Recognition, classification and me-
chanical description of debris flows, Earth-Sci. Rev., 40, 209–
227, 1996.

Coussot, P., Laigle, D., Arattano, M., Deganutti, A., and Marchi, L.: Direct determination of rheological characteristics of debris flow, J. Hydrol., 124, 865–868, 1998.

Domínguez, J. M., Fourtakas, G., Altomare, C., Canelas, R. B., Tafuni, A., García-Feal, O., Martínez-Estévez, I., Mokos, A., Vacondio, R., Crespo, A. J. C., Rogers, B. D., Stansby, P. K., and Gómez-Gesteira, M.: DualSPHysics: from fluid dynamics to multiphysics problems, Computational Particle Mechanics, https://doi.org/10.1007/s40571-021-00404-2, 2021.

Einstein, A.: Zur Theorie der Brownschen Bewegung, Annalen der
Physik, 324, 371–381, 1906.

English, A., Domínguez, J., Vacondio, R., Crespo, A., Stansby, P., Lind, S., Chiapponi, L., and Gómez-Gesteira, M.: Modified dynamic boundary conditions (mDBC) for general-purpose smoothed particle hydrodynamics (SPH): application to tank sloshing, dam break and fish pass problems, Computational Particle Mechanics, https://doi.org/10.1007/s40571-021-00403-3, 2021.

Freydier, P.: Dynamique interne au front d'écoulements à surface libre. Application aux laves torrentielles, PhD thesis, Université Grenoble Alpes, 2017.

Freydier, P., Chambon, G., and Naaim, M.: Experimental character-
ization of velocity fields within the front of viscoplastic surges
down an incline, Journal of Non-Newtonian Fluid Mechanics,
240, 56–69, https://doi.org/10.1016/j.jnnfm.2017.01.002, 2017.

George, D. L. and Iverson, R. M.: A depth-averaged debris-flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests, P. Roy. Soc. A-Math. Phy., 470, 20130820, https://doi.org/10.1098/rspa.2013.0820, 2014.

Gibson, S., Floyd, I., Sánchez, A., and Heath, R.: Comparing single-phase, non-Newtonian approaches with experimental results: validating flume-scale mud and debris flow in HEC-RAS, Earth Surf. Proc. Land., 46, 540–553, 2021.

Gray, J. M. N. T.: Particle segregation in dense gran-
ular flows, Annu. Rev. Fluid Mech., 50, 407–433,
https://doi.org/10.1146/annurev-fluid-122316-045201, 2018.

Guazzelli, E. and Pouliquen, O.: Rheology of dense granular suspensions, J. Fluid Mech., 852, https://doi.org/10.1017/jfm.2018.548, 2018.

Guyon, E., Hulin, J.-P., and Petit, L.: Hydrodynamique physique 3e édition, EDP Sciences, ISBN 9782759805617, 2012.

Han, Z., Chen, G., Li, Y., Xu, L., Zheng, L., and Zhang, Y.: A new approach for analyzing the velocity distribution of debris flows at typical cross-sections, Nat. Hazards, 74, 2053–2070, 2014.

Hungr, O.: Classification and terminology, in: Debris-Flow Hazards and Related Phenomena, edited by: Jakob, M. and Hungr, O., Springer Praxis Books, Springer, https://doi.org/10.1007/3-540-27129-5_16, 9–23, 2005.

Hunt, B.: Newtonian fluid mechanics treatment of debris flows and avalanches, J. Hydrol., 120, 1350–1363, 1994.
Hürlimann, M., Coviello, V., Bel, C., Guo, X., Berti, M., Graf, C., Hübl, J., Miyata, S., Smith, J. B., and Yin, H.-Y.: Debris-flow monitoring and warning: review and examples, Earth-Sci. Rev., 199, 102981, https://doi.org/10.1016/j.earscirev.2019.102981, 2019.

Iverson, R. M.: The physics of debris flows, Rev. Geophys., 35, 245–296, https://doi.org/10.1029/97RG00426, 1997.

Jean, M. and Moreau, J. J.: Unilaterality and dry friction in the dynamics of rigid body collections, in: 1st Contact Mechanics International Symposium, 31–48, https://doi.org/10.1023/A:1008289716620, 1992.

Johnson, C., Kokelaar, B., Iverson, R. M., Logan, M., LaHusen, R., and Gray, J.: Grain-size segregation and levee formation in geophysical mass flows, J. Geophys. Res.-Earth, 117, https://doi.org/10.1029/2011JF002185, 2012.

Kaitna, R., Hsu, L., Rickenmann, D., and Dietrich, W.: On the development of an unsaturated front of debris flow, Italian Journal of Engineering Geology and Environment, https://doi.org/10.4408/IJEGE.2011-03 351–358, 2011.

Kong, Y., Li, X., and Zhao, J.: Quantifying the transition of impact mechanisms of geophysical flows against flexible barrier, Eng. Geol., 289, 106188, https://doi.org/10.1016/j.enggeo.2021.106188, 2021.

Krieger, I. M. and Dougherty, T. J.: A mechanism for non newtonian flow in suspensions of rigid spheres, Transactions of the Society of Rheology, 3, 137–152, https://doi.org/10.1122/1.548848, 1959.

Kromkamp, J., van den Ende, D., Kandhai, D., van der Sman, R., and Boom, R.: Lattice Boltzmann simulation of 2D and 3D non-Brownian suspensions in Couette flow, Chem. Eng. Sci., 61, 858–873, https://doi.org/10.1016/j.ces.2005.08.011, 2006.

Labbé, M.: Modélisation numérique de l'interaction d'un écoulement de fluide viscoplastique avec un obstacle rigide par la méthode SPH: Application aux laves torrentielles, PhD thesis, Université Grenoble Alpes, 2015.

Laigle, D. and Coussot, P.: Numerical modeling of mudflows, J. Hydrol., 123, 617–623, https://doi.org/10.1061/(ASCE)0733-9429(1997)123:7(617) 997.

Laigle, D. and Labbe, M.: SPH-based numerical study of the impact of mudflows on obstacles, International Journal of Erosion Control Engineering, 10, 12, https://doi.org/10.13101/ijece.10.56, 2017.

Lapillonne, S.: Source code and parameters for viscous surges and mixtures with DualSPHysics v5.2, Zenodo [code, data set], https://doi.org/10.57745/VCTHIS, 2024.

Lapillonne, S., Fontaine, F., Liebault, F., Richefeu, V., and Piton, G.: Debris-flow surges of a very active alpine torrent: a field database, Nat. Hazards Earth Syst. Sci., 23, 1241–1256, https://doi.org/10.5194/nhess-23-1241-2023, 2023.

Leonardi, A., Wittel, F. K., Mendoza, M., and Herrmann, H. J.: Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions, Computational Particle Mechanics, 1, 3–13, https://doi.org/10.1007/s40571-014-0001-z, 2014.

Leonardi, A., Cabrera, M., Wittel, F. K., Kaitna, R., Mendoza, M., Wu, W., and Herrmann, H. J.: Granular front formation in free-surface flow of concentrated suspensions, Phys. Rev. E, 92, 052204, https://doi.org/10.1103/PhysRevE.92.052204, 2015.

Li, X., Zhao, J., and Soga, K.: A new physically based impact model for debris flow, Géotechnique, 71, 674–685, 2021.

Lind, S. J., Xu, R., Stansby, P. K., and Rogers, B. D.: Incompressible smoothed particle hydrodynamics for free-surface flows: a generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves, J. Comput. Phys., 231, 1499–1523, https://doi.org/10.1016/j.jcp.2011.10.027, 2012.

Liu, G.-R. and Liu, M. B.: Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific, ISBN 978-981-256-440-5, 2003.

Marchi, L., Arattano, M., and Deganutti, A. M.: Ten years of debris-flow monitoring in the Moscardo Torrent (Italian Alps), Geomorphology, 46, 1–17, 2002.

Martínez-Estévez, I., Domínguez, J., Tagliafierro, B., Canelas, R., García-Feal, O., Crespo, A., and Gómez-Gesteira, M.: Coupling of an SPH-based solver with a multiphysics library, Comput. Phys. Commun., 283, 108581, https://doi.org/10.1016/j.cpc.2022.108581, 2023.

McArdell, B. W., Hirschberg, J., Graf, C., Boss, S., and Badoux, A.: Illgraben debris-flow characteristics 2019–2022, Envidat [data set], https://doi.org/10.16904/envidat.378, 2023.

Mergili, M. and Pudasaini, S.: r.avaflow – The mass flow simulation tool, https://www.avaflow.org (last access: 10 July 2024), 2014.

Meyrat, G., McArdell, B., Ivanova, K., Müller, C., and Bartelt, P.: A dilatant, two-layer debris flow model validated by flow density measurements at the Swiss illgraben test site, Landslides, 19, 265–276, 2022.

Mitsoulis, E., Abdali, S. S., and Markatos, N. C.: Flow simulation of herschel-bulkley fluids through extrusion dies, The Canadian Journal of Chemical Engineering, 71, 147–160, https://doi.org/10.1002/cjce.5450710120, 1993.

Moreau, J. J.: Evolution problem associated with a moving convex set in a Hilbert space, Journal of Differential Equations, 26, 347–374, 1977.

Nagl, G., Hübl, J., and Kaitna, R.: Velocity profiles and basal stresses in natural debris flows, Earth Surf. Proc. Land., 45, 1764–1776, 2020.

Pang, J.-S. and Stewart, D. E.: Differential variational inequalities, Math. Program., 113, 345–424, 2008.

Pastor, M., Haddad, B., Sorbino, G., Cuomo, S., and Drempetic, V.: A depth-integrated, coupled SPH model for flow-like landslides and related phenomena, International Journal for Numerical and Analytical Methods in Geomechanics, 33, 143–172, 2009.

Pastor, M., Yague, A., Stickle, M., Manzanal, D., and Mira, P.: A two-phase SPH model for debris flow propagation, International Journal for Numerical and Analytical Methods in Geomechanics, 42, 418–448, 2018.

Peng, C., Zhan, L., Wu, W., and Zhang, B.: A fully resolved SPH-DEM method for heterogeneous suspensions with arbitrary particle shape, Powder Technology, 387, 509–526, 2021.

Pirulli, M.: Numerical modelling of landslide runout. A continuum mechanics approach, PhD thesis, Politecnico di Torino, https://iris.polito.it/handle/11583/2499763 (last access: 9 October 2025), 2005.

Piton, G., Fontaine, F., Bellot, H., Liébault, F., Bel, C., Recking, A., and Hugerot, T.: Direct field observations of massive bedload and debris flow depositions in open check dams, in: E3S Web of Conferences (Proc. of River Flow 2018 – Ninth International Conference on Fluvial Hydraulics), edited by: Paquier, A. and Riviere, N., Vol. 40, EDP Sciences, https://doi.org/10.1051/e3sconf/20184003003, 1–8, 2018.
Piton, G., Goodwin, S. R., Mark, E., and Strouth, A.: Debris flows, boulders and constrictions: a simple framework for modeling jamming, and its consequences on outflow, *J. Geophys. Res.-Earth*, 127, e2021JF006447, https://doi.org/10.1029/2021JF006447, 2022.

Poudyal, S., Choi, C., Song, D., Zhou, G. G., Yune, C.-Y., Cui, Y., Leonardi, A., Busslinger, M., Wendeler, C., Piton, G., Moase, E., and Strouth, A.: Review of the mechanisms of debris-flow impact against barriers, in: International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, 1027–1033, ISBN 9780578510828, 2019. Poudyal, S., Choi, C., Song, D., Zhou, G. G., Yune, C.-Y., Cui, Y., Leonardi, A., Busslinger, M., Wendeler, C., Piton, G., Moase, E., and Strouth, A.: Review of the mechanisms of debris-flow impact against barriers, in: International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, 2019, 1027–1033, 2019.

Quinlan, N. J., Basa, M., and Lastiwka, M.: Truncation error in mesh-free particle methods, *Int. J. Numer. Meth. Eng.*, 66, 2064–2085, 2006.

Radjai, F. and Richefeu, V.: Contact dynamics as a nonsmooth discrete element method, *Mechanics of Materials*, 41, 715–728, 2009.

Recking, A., Richard, D., and Degoutte, G.: *Torrents et Rivières de Montagne: dynamique et aménagement*, Quae Editions, ISBN 978-2-7592-1999-5, 2013.

Reid, M. E., Coe, J. A., and Brien, D. L.: Forecasting inundation from debris flows that grow volumetrically during travel, with application to the Oregon Coast Range, USA, *Geomorphology*, 273, 396–411, 2016.

Robb, D. M., Gaskin, S. J., and Marongiu, J.-C.: SPH-DEM model for free-surface flows containing solids applied to river ice jams, *J. Hydrol.*, 54, 27–40, 2016.

Rodriguez-Paz, M. X. and Bonet, J.: A corrected smooth particle hydrodynamics method for the simulation of debris flows, *Numerical Methods for Partial Differential Equations*, 20, 140–163, https://doi.org/10.1002/num.10083, 2004.

Roquet, N. and Saramito, P.: An adaptive finite element method for Bingham fluid flows around a cylinder, *Comput. Method. Appl. M.*, 192, 3317–3341, https://doi.org/10.1016/S0045-7825(03)00262-7, 2003.

Schöffl, T., Nagl, G., Koschuch, R., Schreiber, H., Hübl, J., and Kaitna, R.: A perspective of surge dynamics in natural debris flows through pulse-Doppler radar observations, *J. Geophys. Res.-Earth*, 128, e2023JF007171, https://doi.org/10.1029/2023JF007171, 2023.

Simoni, A., Bernard, M., Berti, M., Boreggio, M., Lanzoni, S., Stancanelli, L. M., and Gregoretti, C.: Runoff-generated debris flows: observation of initiation conditions and erosion–deposition dynamics along the channel at Cancia (eastern Italian Alps), *Earth Surf. Proc. Land.*, 45, 3556–3571, 2020.

Takahashi, T.: Debris Flow Mechanics, Prediction and Countermeasures, 2nd edn., CRC Press, ISBN 9780429224355, 2014.

Tasora, A., Serban, R., Mazhar, H., Pazouki, A., Melanz, D., Fleischmann, J., Taylor, M., Sugiyama, H., and Negrut, D.: Chrono: an open source multi-physics dynamics engine, in: High Performance Computing in Science and Engineering – Lecture Notes in Computer Science, edited by: Kozubek, T., Springer, 19–49, https://doi.org/10.1007/978-3-319-40361-8_2, 2016.

Theule, J. I., Crema, S., Marchi, L., Cavalli, M., and Comiti, F.: Exploiting LSPIV to assess debris-flow velocities in the field, *Nat. Hazards Earth Syst. Sci.*, 18, 1–13, https://doi.org/10.5194/nhess-18-1-2018, 2018.

Trujillo-Vela, M. G., Galindo-Torres, S. A., Zhang, X., Ramos-Cañón, A. M., and Escobar-Vargas, J. A.: Smooth particle hydrodynamics and discrete element method coupling scheme for the simulation of debris flows, *Comput. Geotech.*, 125, 103669, https://doi.org/10.1016/j.compgeo.2020.103669, 2020.

Wendland, H.: Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, *Advances in Computational Mathematics*, 4, 389–396, 1995.