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

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.  Albaba 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, spanning 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; Ceccato 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; Ceccato 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 interstitial fluid that is typically highly viscous (Coussot and Meunier, 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 be 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):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{D}\mathit{\rho }}{\mathrm{D}t}}=-\mathit{\rho }\mathrm{\nabla }\cdot \underset{\mathrm{¯}}{u},\text{(2)}& {\displaystyle \frac{\mathrm{D}\underset{\mathrm{¯}}{u}}{\mathrm{D}t}}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }P+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }\cdot \underset{\mathrm{¯}}{\underset{\mathrm{¯}}{\mathit{\tau }}}+\underset{\mathrm{¯}}{g},\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }P={\displaystyle \frac{c_{\mathrm{0}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\gamma }}}\left({\left({\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}\right)}^{\mathit{\gamma }}-\mathrm{1}\right),\end{array}$$

where c₀ is the numerical speed of sound, γ=7 and ρ₀ 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 $\underset{\mathrm{¯}}{r}$ in an ensemble Ω and a function $f:{\mathbb{R}}^{\mathrm{3}}\to \mathbb{R}$ continuous over Ω, the following can always be written:

$$\begin{array}{}\text{(4)}& f\left(\underset{\mathrm{¯}}{r}\right)=\underset{\mathrm{\Omega }}{\int }\mathit{\delta }(\underset{\mathrm{¯}}{r}-{\underset{\mathrm{¯}}{r}}^{\prime })f\left({\underset{\mathrm{¯}}{r}}^{\prime }\right)\mathrm{d}{\underset{\mathrm{¯}}{r}}^{\prime },\end{array}$$

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

$$\begin{array}{}\text{(5)}& \langle f\left(\underset{\mathrm{¯}}{r}\right)\rangle =\underset{\mathrm{\Omega }}{\int }{W}_h(\underset{\mathrm{¯}}{r}-{\underset{\mathrm{¯}}{r}}^{\prime })f\left({\underset{\mathrm{¯}}{r}}^{\prime }\right)\mathrm{d}{\underset{\mathrm{¯}}{r}}^{\prime },\end{array}$$

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:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial {\mathit{\rho }}_a}{\partial t}}={\mathit{\rho }}_a\sum _b\left(({\underset{\mathrm{¯}}{u}}_a-{\underset{\mathrm{¯}}{u}}_b)\cdot \underset{\mathrm{¯}}{\mathrm{\nabla }}{W}_{ab}\right)V_b+h_W{c}_{\mathrm{0}}{\mathcal{D}}_a,\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial {\underset{\mathrm{¯}}{u}}_a}{\partial t}}& =-\sum _b{m}_b\left({\displaystyle \frac{P_a+P_b}{{\mathit{\rho }}_{\mathrm{b}}{\mathit{\rho }}_a}}\right){\underset{\mathrm{¯}}{\mathrm{\nabla }}}_a{W}_{ab}\\ & +\sum _b{m}_b\left({\displaystyle \frac{\mathrm{4}{\mathit{\nu }}_{\mathrm{0}}{\underset{\mathrm{¯}}{r}}_{ab}\cdot {\underset{\mathrm{¯}}{\mathrm{\nabla }}}_a{W}_{ab}}{({\mathit{\rho }}_a+{\mathit{\rho }}_{\mathrm{b}})({\underset{\mathrm{¯}}{r}}_{ab}^{\mathrm{2}}+{\mathit{\eta }}^{\mathrm{2}})}}\right){\underset{\mathrm{¯}}{u}}_a+\underset{\mathrm{¯}}{g},\end{array}\text{(8)}& {\displaystyle \frac{\partial {\underset{\mathrm{¯}}{r}}_a}{\partial t}}={\underset{\mathrm{¯}}{u}}_a,\end{array}$$

where the subscripts a and b refer to the interpolating and neighbouring particles, respectively; ${\underset{\mathrm{¯}}{r}}_i$ is the positional vector of particle i; ${\underset{\mathrm{¯}}{r}}_{ij}={\underset{\mathrm{¯}}{r}}_j-{\underset{\mathrm{¯}}{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; η=0.01h_W, to avoid singularity; 𝒟_i is a density diffusion term (see Supplement); and ν₀ 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:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& W(r,h_W)={\mathit{\alpha }}_D{\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{r}{h_W}}\right)}^{\mathrm{4}}\left(\mathrm{2}{\displaystyle \frac{r}{h_W}}+\mathrm{1}\right)\\ & \text{for }\mathrm{0}<{\displaystyle \frac{r}{h_W}}<\mathrm{2}.\end{array}\end{array}$$

Here, α_D depends on the dimension of the model: ${\mathit{\alpha }}_{\text{2D}}=\frac{\mathrm{7}}{\mathrm{4}}\mathit{\pi }{h}_W^{\mathrm{2}}$ and ${\mathit{\alpha }}_{\text{3D}}=\frac{\mathrm{21}}{\mathrm{16}}\mathit{\pi }{h}_W^{\mathrm{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{\mathrm{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:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& M{\displaystyle \frac{\mathrm{d}\underset{\mathrm{¯}}{V}}{\mathrm{d}t}}=\sum _k{m}_k{\underset{\mathrm{¯}}{f}}_k,\\ & I{\displaystyle \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t}}=\sum _k{m}_k\left(\right({\underset{\mathrm{¯}}{r}}_k-{\underset{\mathrm{¯}}{R}}_{\mathrm{0}})\times {\underset{\mathrm{¯}}{f}}_k)\cdot \underset{\mathrm{¯}}{z},\end{array}\end{array}$$

where the subscript k denotes an SPH boundary particle with a force per unit mass f_k, $\underset{\mathrm{¯}}{V}$ is the linear velocity, M is the mass of the rigid body, R₀ is the position of the centre of mass, I is the moment of inertia and Ω is the rotational velocity at a position $\underset{\mathrm{¯}}{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 $\underset{\mathrm{¯}}{V}$ and Ω, respectively, back to SPH solver. The velocity of each particle boundary in SPH is updated as follows:

$$\begin{array}{}\text{(11)}& {\underset{\mathrm{¯}}{u}}_k=\underset{\mathrm{¯}}{V}+\mathrm{\Omega }\underset{\mathrm{¯}}{z}\times ({\underset{\mathrm{¯}}{r}}_k-{\underset{\mathrm{¯}}{R}}_{\mathrm{0}}),\end{array}$$

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≈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}_{\text{max}}=\mathrm{0}$, where $\underset{\mathrm{¯}}{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:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& p=\mathit{\rho }g(z_{\text{max}}-z)\cos \mathit{\theta },\\ & u_x={\displaystyle \frac{\mathit{\rho }g}{\mathrm{2}\mathit{\mu }}}\left(z_{\text{max}}^{\mathrm{2}}-(z-z_{\text{max}}{)}^{\mathrm{2}}\right)\sin \mathit{\theta },\\ & \stackrel{\mathrm{‾}}{\left|\right|\underset{\mathrm{¯}}{u}\left|\right|}={\displaystyle \frac{\mathit{\rho }g\sin \mathit{\theta }}{\mathrm{3}\mathit{\mu }}}{z}_{\text{max}}^{\mathrm{2}},\end{array}\end{array}$$

where p is the pressure, ρ is the density of the fluid, g is the gravitational force, μ is the viscosity, u is the velocity (with $\stackrel{\mathrm{‾}}{\mathrm{\cdots }}$ denoting an average, the suffix ⋯_x denoting the $\underset{\mathrm{¯}}{x}$ component of the velocity, $\left|\right|\mathrm{\cdots }\left|\right|$ being the magnitude operator and $\underset{\mathrm{¯}}{\mathrm{\cdots }}$ denoting a vector), z is the distance from the incline (in the normal direction), z_max is the free-surface elevation and θ 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.

Table 1 Characteristics of the experiments: density (ρ), viscosity (μ), slope angle (θ), conveyor belt velocity (u_b) and Reynolds number (Re).

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 Ancey, 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 ϕ).

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 ϕ_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 (ϕ_D<0.05) and was extended to any dimension D by Brady (1983) such that

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\mathit{\mu }}_{\mathrm{eq}}={\mathit{\mu }}_{\mathrm{0}}\left(\mathrm{1}+{\displaystyle \frac{D+\mathrm{2}}{\mathrm{2}}}{\mathit{\varphi }}_D\right),\end{array}$$

where μ_eq is the macroscopic viscosity of the mixture and μ₀ 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(ϕ²) or even O(ϕ³) 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, ϕ_m and p:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{{\mathit{\mu }}_{\mathrm{eq}}}{{\mathit{\mu }}_{\mathrm{0}}}}={\left(\mathrm{1}-{\displaystyle \frac{\mathit{\varphi }}{{\mathit{\varphi }}_{\mathrm{m}}}}\right)}^{{\displaystyle }-p{\mathit{\varphi }}_{\mathrm{m}}},\end{array}$$

where ϕ_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).

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 θ. 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 ID₁ and of L=0.4 m and h=0.03 m for experiment ID₂. 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.

Figure 1 Display of the overall flow for simulated flow with SPH and free-surface elevation of the experimental measurements for ID₁ with d_p=0.001 m and C_h=1.0, $x^{*}=x/h_{\text{max}}$, and $z^{*}=z/h_{\text{max}}$, where h_max is the theoretical maximal flow height in the uniform zone.

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{\mathrm{2}}{d}_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>\mathrm{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^{\ast }=h/h_{\text{max}}$ and $x^{\ast }=x/h_{\text{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^{\ast }=u/u_{\mathrm{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:

$$\begin{array}{}\text{(15)}& {\mathit{\zeta }}_{\mathrm{f}}={\displaystyle \frac{\left|u_{\mathrm{f}}\right|}{u_{\mathrm{b}}}}=\left|u_{\mathrm{f}}^{*}\right|,\end{array}$$

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

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

Figure 2b shows the relative error in the flow height in the uniform zone for which ${\partial }_{x}h=\mathrm{0}$, which is defined as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{\mathit{\zeta }}_h={\displaystyle \frac{|h_{\mathrm{u}}-h_{\text{max}}|}{h_{\text{max}}}}=|h_{\mathrm{u}}^{\ast }-\mathrm{1}|,\end{array}$$

where $h_{\mathrm{u}}^{*}$ is the normalized height in the uniform zone, taken in the uniform zone with $x^{*}>\mathrm{9}$.

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

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{\mathit{\zeta }}_{\text{fs}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_x}}\sum _{i=\mathrm{0}}^{N_x}{\left({\displaystyle \frac{h^{*}\left(x_i\right)-h_{\text{ref}}^{*}\left(x_i\right)}{h_{\text{ref}}^{*}\left(x_i\right)}}\right)}^{\mathrm{2}}},\end{array}$$

where h^*(x_i) is the normalized free-surface elevation of the simulated flow at position x_i, $h_{\text{ref}}^{*}\left(x\right)$ 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 ID₁ and ID₂ are tested (see Table 1). In Fig. 2, full symbols show the convergence study for experiment ID₁, whereas void symbols show the convergence study for experiment ID₂. 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 ID₁. The same features were monitored for experiment ID₂ 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_{\mathrm{f}}^{*}=\mathrm{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.

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_{\mathrm{f}}^{*}=\mathrm{0.76}$ and $x_{\mathrm{f}}^{*}=\mathrm{1.5}$ is close to the front tip, where ∂_xz is significant and $x_{\mathrm{f}}^{*}=\mathrm{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).

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

$$\begin{array}{}\text{(18)}& {\mathit{\zeta }}_{\mathrm{u}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_h}}\sum _{j=\mathrm{0}}^{N_h}{\left(u^{*}\left(h_j\right)-u_{\text{ref}}^{*}\left(h_i\right)\right)}^{\mathrm{2}}},\end{array}$$

where u^*(h_j) is the normalized longitudinal velocity of the simulated flow at elevation h_j in the flow, $u_{\text{ref}}^{*}\left(h_j\right)$ 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 ζ_u,j is plotted for each (d_p, C_h). Convergence is shown for each smoothing coefficient. The 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 agrement 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 DualSPHysics, 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 $\underset{\mathrm{¯}}{x}$ direction. The geometry is described in Fig. 4, where the boundaries are periodic along the $\underset{\mathrm{¯}}{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):

$$\begin{array}{}\text{(19)}& {\mathit{\mu }}_{\mathrm{eq}}={\displaystyle \frac{B^{\mathrm{2}}}{\mathrm{8}{u}_{x,\text{max}}}}\left|{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}x}}\right|.\end{array}$$

The maximal velocity along $\underset{\mathrm{¯}}{x}$ (u_x,max) is taken as the average of the velocities along $\underset{\mathrm{¯}}{x}$ of the fluid particles with positions 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.

Figure 4 Experimental setup for the loaded Poiseuille case; the pressure gradient is introduced through , and the colour bar represents velocity in the longitudinal direction.

Table 2 Summary of the designs of the numerical experiments; note that for concentrations above 0.4, $d_p/r=\mathrm{4.5}\mathrm{\%}$.

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 ϕ_m. The black lines represent the uncertainty due to the spatial averaging.

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

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:

$$\begin{array}{}\text{(20)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_s}}\sum _{i=\mathrm{0}}^{N_s}{\left({\displaystyle \frac{{\mathit{\mu }}_{\mathrm{eq}}}{{\mathit{\mu }}_{\mathrm{0}}}}-{\displaystyle \frac{{\mathit{\mu }}_{\text{est}}}{{\mathit{\mu }}_{\mathrm{0}}}}\right)}^{\mathrm{2}}},\end{array}$$

where N_s is the number of tested concentrations and μ_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 ϕ reaches ϕ_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.

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.

The flow is initialized with a constant flow depth h_surge=1.4 m; u_b is then chosen so that the Froude number of the initial conditions reaches Fr≈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

Table 3 Characteristics of the 2D debris flow experiment.

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 \mathrm{1000}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$, with a channel width of B=8 m. Using a simplistic rectangular approximation, this lead to a length of the flow front and body of a corresponding event $\frac{V}{h_{\text{max}}B}\approx \mathrm{90}\mathrm{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_truncature=20 cm and up to h_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 \mathrm{2600}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ seem to float in debris flows, the average density of which is ≈1800–2000 kg m⁻³. To mimic this effect that cannot be explicitly modelled due to the lack of representation of smaller grain sizes, the boulder relative density $\frac{{\mathit{\rho }}_{\mathrm{b}}}{\mathit{\rho }}$, where ρ_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.

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.

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.

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{array}{}\text{(21)}& \begin{array}{rl}& \mathit{Fr}={\displaystyle \frac{\overline{u}}{\sqrt{g{h}_{\text{max}}}}},\\ & {\mathit{\mu }}_{\mathrm{eq}}={\displaystyle \frac{\mathit{\rho }g\sin \mathit{\theta }}{\mathrm{3}\overline{u}}}{h}_{\text{max}}^{\mathrm{2}}.\end{array}\end{array}$$

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.

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 $\mathit{\mu }/{\mathit{\mu }}_{\mathrm{0}}=\mathrm{1}$. Velocity fields and a general view of the different simulations for (c) ϕ=0.53, (d) ϕ=0.29, (e) ϕ=0.13 and (f) ϕ=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.

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 effective 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 non-linear trend: as the concentration of grains increases, the viscosity of the mixture exhibits a more-than-proportional growth. Both h_max and $\stackrel{\mathrm{‾}}{u}$ are used to compute Fr and μ, 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 characteristics is of major interest for applications on mudflow 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 DualSPHysics 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.