# An SPH study on viscoplastic surges overriding mobile beds: The many regimes of entrainment

Saoirse Robin Goodwin, Suzanne Lapillonne, Guillaume Piton, Guillaume Chambon

► To cite this version:

Saoirse Robin Goodwin, Suzanne Lapillonne, Guillaume Piton, Guillaume Chambon. An SPH study on viscoplastic surges overriding mobile beds: The many regimes of entrainment. Computers & Geosciences, 2023, 181, pp.105476. 10.1016/j.cageo.2023.105476. hal-04301607

HAL Id: hal-04301607

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

Submitted on 23 Nov 2023

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Research paper

# An SPH study on viscoplastic surges overriding mobile beds: The many regimes of entrainment

Saoirse Robin Goodwin<sup>a,b</sup>, Suzanne Lapillonne<sup>a</sup>, Guillaume Piton<sup>a</sup>, Guillaume Chambon<sup>a,*</sup>

## ARTICLE INFO

**Keywords:**

Landslides

Smoothed particle hydrodynamics

Viscoplastic

Entrainment

## ABSTRACT

Flow-type landslides entrain mobile bed material, but the processes involved are diverse and require systematic study. We perform direct numerical simulations using the open-source SPH package DualSPHysics with a regularized Herschel-Bulkley rheology. We then compare model output with physical test data, and hence investigate the effects of varying the bed yield stress $\tau_{y,b}$ and bed depth $h_b$, interpreting the results using a newly-identified set of dimensionless numbers. Results reveal diverse interaction regimes between surges and mobile beds, including “rigid bed”, “lubrication”, “shallow ploughing”, “surfing”, “plunging”, and “deep ploughing”. Shallow, borderline-stable beds “lubricate” the surge: once destabilized, these beds cause strong acceleration of the combined flow front. Deeper borderline-stable beds allow the surge material to “plunge” downward, massively displacing bed material upward and downstream. For stabler beds, “ploughing” and “surfing” are associated with intermediate and high values of $\tau_{y,b}$, respectively. In both cases, beds retard the surge, with mobile dams forming for “ploughing” regimes. Across all regimes identified, the influence of $\tau_{y,b}$ is non-monotonic, with intermediate values decelerating the combined flow fronts the most. Furthermore, the different interaction regimes exhibit unique velocity profiles. We develop phase diagrams based on three dimensionless numbers, demarcating these regimes.

## 1. Introduction

Numerical simulations of flow-type landslides are increasingly used, for instance for conducting risk assessments and for predicting the consequences of potential changes in flow properties under a changing climate. However, modelling flow-type landslides is faced with numerous challenges, including capturing the complex rheological behaviour of the materials, and accounting for entrainment mechanisms that result from interactions between the surges and the underlying terrain (“beds”). Fig. 1 illustrates two cases where clay-rich flow-type landslides entrained significant amounts of material in Alpine channels, reflected by the depth of the incisions.

Mass and momentum exchange resulting from entrainment strongly influence the flow volume and dynamics. Entrainment processes are inherently unsteady (Iverson, 1997, 2012; Iverson and Ouyang, 2014), vary greatly from the head to the body to the tail of surges (Berger et al., 2011; Song and Choi, 2021), and are influenced by the mechanical properties of the bed and the channel geometry. Furthermore, several distinct surge/bed interaction mechanisms have been catalogued for geophysical flows, including debris flows (Iverson and Ouyang, 2014), dry granular flows (Nikooei and Manzari, 2021), and

snow avalanches (Gauer and Issler, 2004; Li et al., 2022). Notwithstanding, entrainment processes can seldom be investigated directly in the field (with Berger et al., 2011; McCoy et al., 2012, constituting valuable exceptions), and the underlying physical mechanisms remain poorly understood. As such, entrainment modelling remains largely based on empirical parametrizations, the validity of which is often difficult to assess (see for instance Issler, 1998; McDougall and Hungr, 2005; Pirulli and Pastor, 2012; Iverson and Ouyang, 2014).

Alternatively, Direct Numerical Simulations (DNS) in 2D or 3D can provide valuable insight into the mass- and momentum-exchange processes that occur during surge-bed interaction. Indeed, Nikooei and Manzari (2020, 2021) used a 2D incompressible Smoothed Particle Hydrodynamics (SPH) model to study channelized granular surges overriding mobile beds, noting the existence of different entrainment regimes. Li et al. (2022) used a 2D MPM model to investigate interaction between snow avalanches and the underlying snow cover, distinguishing between “erosion” and “entrainment” processes. Results from Mangeney et al. (2007) using a 2D partial-fluidization model also show that bed depth affects entrainment dynamics, at least for dry granular flows.

<footnote><sup>a</sup> University Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, IGE, Grenoble, 38000, France</footnote>

<footnote><sup>b</sup> Institute of Applied Geosciences, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany</footnote>

<footnote>* Corresponding author.
E-mail addresses: srgoodwin@protonmail.com (S.R. Goodwin), guillaume.chambon@inrae.fr (G. Chambon).</footnote>
Nomenclature

<table><tr><td>α</td><td>Artificial viscosity factor</td></tr><tr><td>U&#773;<sub>b</sub></td><td>Average velocity in mobile bed</td></tr><tr><td>U&#773;<sub>s</sub></td><td>Average velocity in surge</td></tr><tr><td>Γ</td><td>Dimensionless number denoting ratio of bed and surge depths</td></tr><tr><td>γ</td><td>Polytropic constant</td></tr><tr><td>κ</td><td>Constant for CFL condition</td></tr><tr><td>μ<sub>m</sub></td><td>Viscosity of point m</td></tr><tr><td>μ<sub>eff</sub></td><td>Effective viscosity</td></tr><tr><td>μ<sub>reg</sub></td><td>Regularized (maximum) effective viscosity</td></tr><tr><td>ρ</td><td>Density</td></tr><tr><td>ρ<sub>0</sub></td><td>Neutral density</td></tr><tr><td>ρ<sub>max</sub></td><td>Maximum density possible in SPH calculations</td></tr><tr><td>ρ<sub>min</sub></td><td>Minimum density possible in SPH calculations</td></tr><tr><td>ΣU<sub>b</sub></td><td>Total velocity in mobile bed</td></tr><tr><td>ΣU<sub>s</sub></td><td>Total velocity in surge</td></tr><tr><td>τ<sub>0</sub></td><td>“Neutral” yield stress (for parametric study only)</td></tr><tr><td>τ<sub>y</sub></td><td>Yield stress</td></tr><tr><td>τ<sub>y,b</sub></td><td>Mobile bed yield stress</td></tr><tr><td>τ<sub>y,s</sub></td><td>Surge yield stress</td></tr><tr><td>τ<sub>zz</sub> - P</td><td>Bed-perpendicular total normal stresses</td></tr><tr><td>θ</td><td>Channel slope</td></tr><tr><td>V&#773;</td><td>Vector of partial derivatives</td></tr><tr><td>γ&#775;</td><td>Shear rate tensor</td></tr><tr><td>τ&#775;</td><td>Extra-stress tensor</td></tr><tr><td>g&#775;</td><td>Gravity-induced acceleration</td></tr><tr><td>U&#773;</td><td>Velocity</td></tr><tr><td>A<sub>CFL</sub></td><td>“Factor of safety” for CFL calculation</td></tr><tr><td>c<sub>0</sub></td><td>Speed of sound</td></tr><tr><td>D</td><td>Total derivative</td></tr><tr><td>dp</td><td>Point size</td></tr><tr><td>G<sub>d</sub></td><td>Dynamic yield number: dimensionless number denoting characteristic stress at the surge-bed interface</td></tr><tr><td>G<sub>s</sub></td><td>Static yield number: dimensionless number denoting inherent bed stability</td></tr><tr><td>h</td><td>Kernel radius (smoothing length)</td></tr><tr><td>h&#773;<sub>b</sub></td><td>Mobile bed depth</td></tr><tr><td>h&#773;<sub>s</sub></td><td>Surge depth</td></tr><tr><td>k</td><td>Consistency index</td></tr><tr><td>k<sub>0</sub></td><td>“Neutral” coefficient for speed of sound (for parametric study only)</td></tr><tr><td>k<sub>sound</sub></td><td>Coefficient for speed of sound</td></tr><tr><td>L&#773;<sub>b</sub></td><td>Mobile bed length</td></tr><tr><td>M</td><td>Mach number</td></tr><tr><td>n</td><td>Flow index</td></tr><tr><td>P</td><td>Pressure</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>U<sub>x</sub></td><td>Channel-parallel velocity</td></tr></table>

<table><tr><td>W</td><td>Kernel-sampled property</td></tr><tr><td>x</td><td>Channel-parallel direction</td></tr><tr><td>x<sub>front, flow</sub></td><td>Position of the flow front on a mobile bed</td></tr><tr><td>x<sub>front, surge</sub></td><td>Position of the surge front on a mobile bed</td></tr><tr><td>x<sub>front</sub></td><td>Position of the flow front on an immobile bed</td></tr><tr><td>z</td><td>Channel-normal direction</td></tr></table>

Physical experiments offer another alternative for investigating entrainment. Bates and Ancey (2017) presented experimental work on the interaction between channelized viscoplastic surges and mobile beds. Both surges and beds were made of the same viscoplastic fluid. These authors varied the slope $θ$, the length of the mobile bed $L_b$ and its depth $h_b$. Results revealed that both $h_b$ and $L_b$ greatly affect the entrainment mechanisms, and highlighted the fact that the flow front resulting from bed entrainment can be accelerated, relative to cases without mobile beds. Only shallow mobile beds were considered in this study, such that the surges always “scoured out” the bed material. The effects of larger $h_b$ were unexplored, as were the effects of varying the rheological parameters of the mobile bed (e.g. its yield stress).

In this work, we use numerical simulations to obtain further insight into the physical mechanisms and regimes involved in entrainment processes. We specifically investigate the influences of bed depth ($h_b$) and the mechanical properties of the bed on characteristics such as the flow reach, velocity profiles, and momentum transfer into the bed. We limit the study scope to a viscoplastic rheology, which can be viewed as a prototype constitutive law for capturing a transition between a static and a flowing behaviour at a particular stress. Furthermore, such viscoplastic rheologies have been shown appropriate for modelling clay-rich flows (see Coussot, 1994; Huang and García, 1998; Coussot et al., 1998; Pastor et al., 2015). We use a 2D Smoothed Particle Hydrodynamics (SPH) model from the open-source software DualSPHysics of Domínguez et al. (2021). We implement a regularized viscoplastic (Herschel-Bulkley) rheology, and compare our model output against experimental data from Bates and Ancey (2017). The model is then used to simulate gravity-driven surges overriding mobile beds across a large range of parameters. Finally, we rationalize the results by presenting phase diagrams based on three governing dimensionless numbers, distinguishing between various entrainment regimes for the parameter combinations that we investigated.

## 2. Numerical modelling

The SPH method (Monaghan, 1988) is well-suited to modelling free-surface flows (Liu and Liu, 2003). It can efficiently handle large deformations and multiple materials, as well as complex rheologies. Indeed, it has been used to model viscoplastic free-surface surges in different contexts (e.g. Laigle et al., 2007; Chambon et al., 2011; Jiang et al., 2011; Labbé and Laigle, 2013; Ghaïtanellis et al., 2018). We used the open-source GPU-accelerated SPH model DualSPHysics (Domínguez et al., 2021), which can model Newtonian (e.g. Reece et al., 2020; Bartzke et al., 2021) and non-Newtonian materials (e.g. Fourtakas et al., 2013; Fourtakas and Rogers, 2016; Zubeldia et al., 2018). We implemented a regularized viscoplastic Herschel-Bulkley constitutive law into the code for the GPU-accelerated non-Newtonian version of DualSPHysics (version: “DualSPHysics5-NNewtonian v1.005 (v5.0.164) (21-11-2020)”), originally presented in Zubeldia et al. (2018).

### 2.1. Governing equations

SPH solves the classical mass and momentum conservation equations given by:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \underline{U}) = 0 \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Evidence of entrainment by clay-rich flow-type landslides in the French Alps: (a) community of the Désert de Valjeffrey; (b) community of Pontamafrey Montpascal.</figcaption></figure>

$$ \frac{DU}{Dt} = \frac{1}{\rho} \left( -\nabla P + \nabla \cdot \underline{\tau} \right) + g \quad (2) $$

where $\rho$ is density; $t$ is time; $\nabla$ is the Nabla operator; $U$ is velocity; $D$ indicates the total derivative; $P$ is pressure; $\underline{\tau}$ is the extra-stress tensor; and $g$ is gravity-induced acceleration. In the weakly-compressible SPH method, pressure is computed using Tait's equation of state (see Domínguez et al., 2021):

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

where $c_0$ is a speed of sound (see Appendix A); $\rho_0$ is the “neutral” density; and $\gamma$ is a constant.

In this study, we used the quintic Wendland kernel W for SPH interpolation, with a smoothing length $h$ (see Wendland, 1995; Macià et al., 2011). We also adopted extra algorithms for stability: namely the “density diffusion” algorithm by Fourtakas and Rogers (2016), and the shifting algorithm by Skillen et al. (2013). Corresponding parameters are given in Appendix B.

## 2.2. Regularized Herschel-Bulkley rheology

The extra-stress tensor $\underline{\tau}$ (see Eq. (2)) is computed according to the viscoplastic Herschel-Bulkley constitutive law (e.g. Herschel and Bulkley, 1926; Ancey, 2007; Balmforth et al., 2014):

$$ \begin{cases} \dot{\underline{\tau}} = 0 & \text{if } |\underline{\tau}| < \tau_y \\ \underline{\tau} = \left( \frac{\tau_y}{|\dot{\underline{\tau}}|} + k|\dot{\underline{\tau}}|^{n-1} \right) \frac{\dot{\underline{\tau}}}{\underline{\tau}} & \text{if } |\underline{\tau}| \ge \tau_y \end{cases} \quad (4) $$

where $\tau_y$ is the yield stress; $k$ is the consistency; $n$ is the flow index; $\dot{\underline{\tau}}$ denotes the shear rate: $\dot{\underline{\tau}}_{ij} = \partial_j U_i + \partial_i U_j$. The tensor norms are defined as $|\dot{\underline{\tau}}| = [0.5 \sum_i \dot{\tau}_{ij} \dot{\tau}_{ij}]^{0.5}$ and $|\underline{\tau}| = [0.5 \sum_i \tau_{ij} \tau_{ij}]^{0.5}$. For $|\underline{\tau}| < \tau_y$, the extra-stress tensor $\underline{\tau}$ is indeterminate. In practice, different solutions can be employed to deal with this indeterminacy (Saramito and Wachs, 2016), including regularization approaches (e.g. Frigaard and Nouar, 2005; Putz et al., 2009; Saramito and Wachs, 2016). Here, we chose to regularize the constitutive law by capping the effective viscosity $\mu_{eff}$ at a value $\mu_{reg}$, such that:

$$ \begin{cases} \mu_{eff} = \mu_{reg} & \text{if } \mu_{eff} \ge \mu_{reg} \\ \mu_{eff} = \frac{\tau_y}{|\dot{\underline{\tau}}|} + k|\dot{\underline{\tau}}|^{n-1} & \text{if } \mu_{eff} < \mu_{reg} \end{cases} \quad (5) $$

The viscosity cap $\mu_{reg}$ was implemented as an adjustable parameter in the DualSPHysics code. Low values of $\mu_{reg}$ tend to artificially reduce the effects of the yield stress on flow behaviour, whereas high values

require extremely fine time-steps, and lead to tensile instabilities (see Appendix A for an investigation thereon). Setting $\mu_{reg}$ thus necessitates compromise.

## 2.3. Simulation setup and procedure

Fig. 2a schematically shows the simulation setup used in this study, based on the experimental flume of Bates and Ancey (2017). These authors filled a gated reservoir with Carbopol; the fluid volume per unit width was 0.05 m$^2$, and the density $\rho$ was 1000 kg/m$^3$. The mobile bed depth was 3 < $h_b$ < 9 mm; its length was 100 < $L_b$ < 300 mm; and the channel slope was 12 < $\theta$ < 24°. We recreated this setup numerically. In the simulations we set $L_b$ = 300 mm and $\theta$ = 20°, and considered bed depths in the range 3 < $h_b$ < 46 mm.

Fig. 3a replots rheometric data from Bates (2015). The measured yield stress $\tau_y$ appear to be strongly variable among the tests, lying between 58 and 130 Pa typically. As shown in Appendix C, such uncertainties on the rheological parameters can have significant effects on the simulation results. In the following, we set $\tau_{y,s}$ = 58 Pa for the surge material. In the experiments of Bates and Ancey (2017), the bed and surge yield stress, denoted as $\tau_{y,b}$ and $\tau_{y,s}$ respectively, were identical. Here, we extended the results by exploring values of bed yield stress $\tau_{y,b}$ ranging between 18 and 258 Pa. Additionally, we also considered a control case with a rigid bed, characterized by $\tau_{y,b}$ = 999 Pa. We set $k$ = 35 Pa s$^n$ and $n$ = 0.33 for both surges and mobile beds, following Bates and Ancey (2017).

Following an extensive sensitivity study (Appendix A), we set the viscosity cap at $\mu_{reg}$ = 10 kPa s. Fig. 3b shows the relation between the simple shear stress $|\underline{\tau}|$, the shear rate $\dot{\underline{\tau}}$ and the effective viscosity $\mu_{eff}$ of the surge in simple shear as predicted by the regularized Herschel-Bulkley law implemented in our SPH model.

In the experiments of Bates and Ancey (2017), dried Carbopol coated the boundaries, enforcing a no-slip condition. In our SPH simulations, three layers of fixed points form the outer boundaries of the channel. No-slip was enforced through eight extra layers of fluid particles (Fig. 2b), whose yield stress was set to 1 kPa, i.e. much higher than $\tau_{y,b}$ or $\tau_{y,s}$. Hence, these “boundary” particles moved negligibly.

Simulations comprised a “dam-break” segment (reused for every test) and a “downstream segment” (unique for each test). For the “dam-break segment”, the surge material initially occupied the holding reservoir, with a horizontal free surface (see Fig. 2a). A single “dam-break segment” simulation was run for 2.4 s of physical time, corresponding to the surge front reaching a position 83 mm upstream of the mobile bed (Fig. 2b). We then extracted the free surface of the surge at this time. For each “downstream segment” simulation, the surge was regenerated using the free surface from the “dam-break segment”. Following a process designed to recover the original velocity profile at the start of the mobile bed zone, a uniform velocity substantially higher than the original mean velocity was applied at the initiation of
<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Side-view schematic of the numerical setup in this study, based on the experimental setup from Bates and Ancey (2017). All dimensions are given in mm. The points x = 0 (the start of the mobile bed) and z = 0 (the surface of the channel) are both marked. (b) Snapshot of the initial frame of the “downstream segment”. The no-slip and channel boundaries are shown in black and fuchsia, respectively.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Rheology of the surge material: shear stress as a function of shear rate in simple shear. (a) Experimental data (replotted from Bates (2015)). (b) Regularized constitutive model used in this study. The right axis shows the effective viscosity $\mu_{eff}$. In (a), the yellow curve corresponds to the nominal rheological parameters given in Bates and Ancey (2017), specifically $\tau_y = 58$ Pa, $k = 35$ Pa s$^a$ and $n = 0.33$. Dashed and dot-dashed curves correspond to Herschel-Bulkley fits of the upper and lower bounds, respectively. In (b), the numerical constitutive law is shown for $\tau_y = 58$ Pa, $\mu_{cg} = 10$ kPa s, $k = 35$ Pa and $n = 0.33$. Dotted lines correspond to the unregularized Herschel-Bulkley behaviour.</figcaption></figure>

the downstream segment. After ~0.8 s, at the start of the mobile bed zone, the original velocity profile obtained at the end of the dam-break segment was effectively recovered (Fig. 4). Note that it would have been simpler to superimpose the original velocity distribution onto the regenerated surge, but the version of DualSPHysics that we used had no such functionality.

The full simulation plan, as well as the full simulation parameters, are given in Appendices B and D respectively. Each “downstream segment” simulation took approximately 24 h.

## 2.4. Dimensionless numbers

In the following, our results are described in terms of four dimensionless numbers. The Froude number $Fr$ for the incoming surge body is calculated from:

$$Fr = \frac{U_{x,s}}{\sqrt{gh_s \cos \theta}} \qquad (6)$$

where $U_{x,s}$ is the longitudinal surge velocity, $h_s$ is the surge body depth, and $\theta$ is the channel slope. In our simulation this Froude number was constant: $Fr \approx 0.3$. The “static yield” number $G_s$ is defined as:

$$G_s = \frac{\tau_{y,b}}{\rho g h_b \sin \theta} \qquad (7)$$

This number describes how far the bed is from yielding under its own weight; $G_s < 1$ implies an unstable bed. The “dynamic yield” number $G_d$ is defined as:

$$G_d = \frac{\tau_{y,b}}{|\tau_i|} \qquad (8)$$

where $|\tau_i|$ is the average stress norm at the surge/bed interface (calculated for a surge overriding a rigid bed). $G_d$ describes the capacity of the surge to mobilize and entrain bed material. Lastly, $\Gamma$ is defined as:

$$\Gamma = \frac{h_b}{h_s} \qquad (9)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Development of the initial velocity profile at the start of the “downstream segment”, for the control case: downstream velocity $U_x$ versus depth along the flow $z$. Several time-steps are shown. For reference, the velocity profile obtained at the end of the “dam-break segment” is also shown. Data was sampled from a region perpendicular to the channel bed, located between $x = -0.1$ and $x = 0.0$ m. Note that a plug zone has already formed at $t = 0.2$ s.</figcaption></figure>

This number describes the capacity for “plunging” of the surge material into the mobile bed.

## 3. Results

Hereafter, we define “surges” as the material coming from the holding reservoir; “mobile beds” as entrainable bed material; and “combined flows” as flows resulting from the interaction between the moving surges and bed material. The “surge front” is the surge's frontmost point; the “bed front” is the frontmost point of the bed that lies above $z = 0$ (Fig. 2); and the “combined flow front” is whichever of the surge or bed fronts is further downstream. The “boundary” refers to the eight layers of high-viscosity points used to enforce the no-slip condition. Lastly, recall that the “control case” corresponds to a surge overriding an immobile (i.e. rigid) bed.

### 3.1. Model comparison: experiments versus simulations

Fig. 5 compares the flow geometry in two experimental tests from Bates and Ancey (2017) with numerical results from our study. Recall that the rheological and geometric parameters in the physical and numerical tests were identical. However, the time when the numerical screenshots were taken was arbitrary since (i) the timestamps for the images of the physical tests are not available, and (ii) the flow velocities in the numerical simulations were faster than those in the physical tests. Possible reasons for this discrepancy include the role of sidewall shearing in the experiments and the uncertainties in the rheological measurements (see Fig. 3a and Appendix C).

Fig. 5a shows that the numerical model reproduces the geometry of the physical flow relatively well, with only minor quantitative differences. Both images show a point of maximum downward penetration of the surge at $x \approx 12$ mm; a sloped surge-bed interface ($12 < x < 55$ mm); and an inflection point at the free surface, 60 mm downstream, where the surge front “draws up” bed material. Fig. 5b shows another very good match between the numerical model and the physical data. This includes the position of the surge-bed interface, which corresponds to the “dimple” at the free surface upstream of the combined flow front. The amount of bed material in front of the surge is higher than in Fig. 5a. We can also note the zone involving a bed-parallel free surface ($60 < x < 110$ mm) behind the combined flow front.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Comparisons of the flow geometry from experiments of Bates and Ancey (2017) with SPH results: (a) $\theta = 20^{\circ}$; and $h_b = 6$ mm; (b) $\theta = 20^{\circ}$ and $h_b = 6$ mm. For both cases, $\tau_{y,s} = \tau_{y,b} = 58$ Pa. Note: the experimental images were treated to remove text and other annotations, add scales, reduce the fish-eye effect and sharpen the image, and adjust the image length.</figcaption></figure>

Fig. 6 quantitatively shows the combined flow front positions for mobile beds, with respect to control cases. Both physical and numerical data are shown. For the physical test ($G_s = 2.9$), the mobile bed accelerates the flow, relative to the control case. (Note: in the experimental data, the surge overriding the mobile bed is slightly ahead of the control case, even before reaching the bed, likely due to minor experimental error.) The SPH model reproduces the physically-observed behaviour well when setting $G_s = 0.9$. For $G_s = 2.4$ and 2.9, the mobile bed initially accelerates the flows, relative to the control case. However, these combined flow fronts eventually decelerate. These discrepancies are likely related to the differences concerning the flow velocities, per the discussion regarding Fig. 5. Note that sidewall shearing and rheological uncertainties affect both the control cases and entrainment cases, although potentially to different extents. Notwithstanding, Figs. 5 to 6 indicate the SPH model's ability to qualitatively reproduce experimentally-observed interaction mechanisms and flow dynamics, including the possibility for mobile beds to accelerate the flow.

### 3.2. Flow dynamics and regimes

Fig. 7 shows surges overriding mobile/immobile beds for values of $h_b$ and $\tau_{y,b}$. We identified several interaction regimes based on the results, exemplified by cases A to F. Cases A to F all show: (a) surge and bed material; (b) downstream velocity $U_x$; (c) bed-perpendicular total normal stresses ($\tau_{zz} - P$); and (d) the proximity to yielding $|\tau_z|/\tau_{y,s}$ or $|\tau_z|/\tau_{y,b}$. The times are: (i) $t = 2$ s (the start of the surge-bed interaction), (ii) $t = 5$ s (an intermediate point) and (iii) $t = 8$ s (where the control surge is near the end of bed zone). Note that the supplementary data includes videos of the results coloured by phase, velocity and yield stress.

#### 3.2.1. Control case – rigid bed

Fig. 7A shows the control case, where the bed is effectively rigid and immobile. The surge velocity increases slightly during downstream
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Comparing a test from Bates and Ancey (2017) (20-6Ca; solid blue line) and three SPH simulations (the broken lines, with different $\tau_{y,b}$). Shown are the combined flow front position on a mobile bed (y-axis) vs. an analogous case on an immobile bed (x-axis). For all cases, $h_b = 6$ mm and $\theta = 20^\circ$. Section 3 gives the definition of the combined flow front position. The solid black reference line assumes an immobile bed for both axes. The grey-shaded area shows the mobile bed region. The inset shows a zoomed-out version of the data.</figcaption></figure>

motion (Fig. 7A(b)); the surge exhibits a hydrostatic profile for the normal stress (Fig. 7A(c)); and yielding occurs in a zone located near the bed interface, with an unyielded plug zone thereover (Fig. 7A(d)). (Note: plugs are parts of the flow where the stress state lies below the yielding criterion and the material locally behaves like a rigid body. In free-surface flows, plugs typically form close to the free-surface, with a shear zone underneath.)

### 3.2.2. Surfing regime

For sufficiently high value of the bed yield stress ($\tau_{y,b}$), a “surfing” regime exists, regardless of the bed depth ($h_b$). Fig. 7B was run with $h_b = 18$ mm, and $\tau_{y,b} = 118$ Pa; corresponding to $G_s \approx 2.0$, $G_d \approx 1.3$, and $\Gamma = 0.6$. Fig. 7B(a)(i) and B(a)(ii) show that the surge front only draws up a little bed material, with minimal deformation. As in the control case, the normal stresses are largely hydrostatic, there is a shear layer located at the surge-bed interface, and the bed remains mostly unyielded. However, the combined flow velocity is clearly lower than in the control case. Note that as the surge decelerates (Fig. 7B(a)(iii)), it draws progressively more bed material upward. Furthermore, less bed material yields than earlier in the simulation.

### 3.2.3. Shallow beds: ploughing/lubrication regimes

Fig. 7C and D show “shallow ploughing” and “lubrication” regimes, respectively. In both cases, $h_b = 6$ mm, whilst $\tau_{y,b}$ are 58 and 18 Pa, respectively. These correspond to $G_s \approx 2.9$, $G_d \approx 0.6$; $\Gamma = 0.2$ (Fig. 7C) and $G_s \approx 0.9$, $G_d \approx 0.2$; $\Gamma = 0.2$ (Fig. 7D). In Fig. 7C, most bed material is displaced by the surge and forms a mobile dam. The surge pushes the dam forward, with the combined flow front comprising bed material. The dam retards the surge front (Fig. 7C(b)), relative to the control case. Furthermore, the normal stresses in the surge are hydrostatic (Fig. 7C(c)); and yielding occurs at the basal interface with the channel boundary (Fig. 7C(d)).

In Fig. 7D, the bed was initially borderline stable, and rapidly accelerates due to its own weight when the surge material arrives. Bed material accumulates directly ahead of the surge (Fig. 7D(a)), similarly to the mobile dam in Fig. 7C(a). However, the displaced bed material *accelerates* the surge relative to the control case. The bed acts as a lubricating layer between the surge and the boundary, hence the term “lubrication”. Notably, in this case the normal stresses are not hydrostatic (Fig. 7D(c)), especially at the surge-bed interface near the

surge front. Yielding concentrates at the basal flow/boundary interface. However, the size of the yielded zone is spatially heterogeneous, being larger in the displaced bed material.

### 3.2.4. Deep beds: ploughing/plunging regimes

Fig. 7E and F show “deep ploughing” and “plunging” regimes, respectively. In Fig. 7E, $h_b = 30$ mm and $\tau_{y,b} = 118$ Pa; therefore $G_s \approx 1.2$, $G_d \approx 1.3$, and $\Gamma = 1.0$). In Fig. 7F, $h_b = 30$ mm and $\tau_{y,b} = 98$ Pa; therefore $G_s \approx 1.0$, $G_d \approx 1.1$ and $\Gamma = 1.0$.

In Fig. 7E(a), a mobile dam develops, similar to the “shallow ploughing” regime already seen, although the bed is not fully pushed upward and downstream. Indeed, much bed material remains below the surge, with the mobile dam still retarding the surge. Velocity is transferred into the mobile bed well ahead of the surge front (Fig. 7E(b)). The normal stresses are rather non-hydrostatic (Fig. 7E(c)), notably when the surge has just reached the mobile bed (Fig. 7E(c)(i)). However, it approaches a hydrostatic profile later on (Fig. 7E(c)(iii)). The yielded zone mostly appears along the surge-bed interface (Fig. 7E(c)). A secondary yielded zone exists in the bed below the surge, although the values of the shear are lower. This morphology of having two yielded zones separated by unyielded bed material — differentiates it from the “shallow ploughing” regime.

Fig. 7F(a)(ii) shows the surge strongly “drawing up” bed material. The surge material penetrates further downward than for “deep ploughing”, with a large slope-perpendicular velocity component. We thus denote this regime as “plunging”. Bed material which was initially borderline-stable is massively displaced downstream, partially due to its own weight. The bed velocity is much higher than for “deep ploughing” (Fig. 7F(b)). Furthermore, the normal stresses of the flow are highly non-hydrostatic (Fig. 7F(c)). The yielded zone at the surge-bed interface (Fig. 7F(d)) is qualitatively similar to the “deep ploughing” case (Fig. 7E(d)). However, there is also a massive yielded zone in the bed, corresponding to its rapid displacement.

## 3.3. Flow reach

Fig. 8 shows the evolution of (a) the *surge front* and (b) the *combined flow front*, both relative to the control case. In Fig. 8a(i), only the case $G_s = 0.9$ corresponds to the *surge material accelerating*, relative to the control case; this is the “lubrication” regime. For higher $G_s$, the *surge* is decelerated due to the “mobile-dam” seen in Fig. 7C; this is the “shallow ploughing” regime. For higher values of $G_s$, the surge front passes through the “surfing” regimes to re-approach the control case (i.e. the “rigid” regime). This non-monotonous behaviour – where the flow reach decreases and then increases as $G_s$ is increased continuously – is notable. Both highly-mobile beds (low $G_s$) or immobile beds ($G_s \to \infty$) cause the furthest flow reaches; in contrast, intermediate values of $G_s$ decelerate surges. In Fig. 8b(i), as expected, for the “lubrication” and “shallow ploughing” regimes, the *combined flow front* generally progresses further than the *surge front*. Notwithstanding, for higher $G_s$ (in the “surfing” regime), the positions of the *combined flow front* and the *surge front* start to converge, as expected.

Fig. 8a(ii) shows data for deeper beds, for which the surge front is not accelerated relative to the control case. In contrast, Fig. 8b(ii) shows that the combined flow front can be accelerated. For the “lubrication” regime ($G_s = 1.0$), the combined flow front exits the mobile bed zone faster than for the “lubrication” regime on the shallower bed (Fig. 8b(i)). In contrast, for the ‘[deep] ploughing’ regime (intermediate $G_s$), the combined flow travels less far than for ‘[shallow] ploughing’ (also intermediate $G_s$, but on Fig. 8a). The surfing regime appears for higher $G_s$, compared to the shallow case. As before, non-monotonous behaviour is observed as $G_s$ is increased.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Flow regime visualizations: (a) the surge and bed coloured by phase; (b) the downstream velocity component $U_x$; (c) the normal stresses perpendicular to the channel bed $\tau_{zz} - P$; (d) the concentration of yield stresses, scaled by the yield stress, i.e. $||f/\tau_{yx}||$ or $||f/\tau_{yb}||$. Times are: (i) $t = 2$ s; (ii) 5 s; and (iii) 8 s. The six cases are: (A) control (no entrainment; $G_s = \infty$, $G_d \approx 10.8$; $\Gamma = 0$); (B) surfing ($G_s \approx 2.0$; $G_d \approx 1.3$; $\Gamma = 0.6$); (C) shallow ploughing ($G_s \approx 2.9$, $G_d \approx 0.6$; $\Gamma = 0.2$); (D) lubrication ($G_s \approx 0.9$, $G_d \approx 0.2$; $\Gamma = 0.2$); (E) deep ploughing ($G_s \approx 1.2$, $G_d \approx 1.3$; $\Gamma = 1.0$); (F) plunging ($G_s \approx 1.0$, $G_d \approx 1.1$; $\Gamma = 1.0$). Note that the supplementary data includes videos of the results coloured by phase, velocity and yield stress.</figcaption></figure>

### 3.4. Velocity profiles

Fig. 9 shows the longitudinal velocity profiles within the flows, for various $\tau_{yx}$ and $h_b$. The profiles are extracted at $100 < x < 105$ mm (see Fig. 5a & the insets in Fig. 9), about a third of the way into the mobile bed, to avoid the “edge effects” at the start of the mobile bed. We divide the measuring zone in the longitudinal direction, and average the velocities in each slice at each time-step, to produce the profiles. In the following, we primarily discuss column (ii), corresponding to

the flow body, close to a steady state. Column (i) corresponds to the highly-unsteady flow front, and is included for completeness.

Fig. 9a corresponds to a low value of $h_b$. The flow profile for the lowest value of $G_s$ (the yellow line), corresponding to the “lubrication” regime, includes a shear zone located at the lower boundary, and a plug zone thereover. Conversely for the highest $G_s$ (the black line, the control case), the shear zone starts at the surge-bed interface, with a plug zone above. Intermediate cases correspond to “surfing” (dark red), with a shear zone near $z = 0$ mm; and “shallow ploughing” (red
<figure><img src="image_10.png" /><figcaption>Fig. 7. (continued).</figcaption></figure>

to gold), with shear zones progressively moving towards the channel boundary. At the combined flow front (Fig. 9a(i)), we see the same spectrum of behaviour, albeit without fully formed plug regions.

Fig. 9b(ii) shows the case of $h_b = 12$ mm. Again, the flow profiles for the lowest and highest values of $G_s$ (the yellow and black lines, respectively), correspond to the “lubrication” and “rigid bed” regimes. The profiles for $G_s = 1.4$ and $1.9$ belong to the “shallow ploughing” regime, whereas the behaviour of $G_s = 2.4$ to $4.4$, starts to resemble “deep ploughing”. Specifically, the velocity profile for these latter cases shows a transition with an inflection point at the interface between the surge and the bed material, with a fully-developed plug for the surge material. The shear rate in the bed material is lower than in the shear zone located in the surge material.

In Fig. 9c, $h_b = 30$ mm. The flow profile for the lowest value of $G_s$ exemplifies the “plunging” regime. Again, an inflection point is clearly visible deep into the mobile bed zone. Notably, the velocity is much faster than for the “deep ploughing” regime in Fig. 9b. Cases with $G_s = 2.2$ and $2.7$, corresponding to the “surfing” regime, show minimal bed deformation. The intermediate case $G_s = 1.6$ clearly shows two shear zones and two plug zones in the mobile bed and the surge, respectively; this is the “deep ploughing” regime.

Overall, these graphs show that the different regimes form a continuum. Furthermore, the trend of how the velocity profiles evolve with $G_s$ is consistent across the different bed thicknesses $h_b$. Finally, note that for a given value of $G_s$, the longitudinal velocity $U_x$ tends to decrease
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Comparing the surge and flow progression downstream against the control case. (a): the surge front position ($x_{\text{front, surge}}$) on a mobile bed; (b) the flow front position ($x_{\text{front, flow}}$) on a mobile bed, against the flow front position on an immobile bed ($x_{\text{front}}$), i.e. the control case. Bed depths are (i) $h_b = 6$ mm, and (ii) $h_b = 18$ mm. The increasing $G_s = \tau_{y,b}/\rho g h_b \sin\theta$ in each plot correspond to progressively less mobile beds. The insets show the bed depth and the monitoring region (i.e. the entire computational domain.).</figcaption></figure>

as $h_b$ increases. This is consistent with the observations regarding the flow reach in Section 3.3.

### 3.5. Momentum transfer

Fig. 10a shows the ratio of total bed and flow momentum, as a function of $G_d$. For the momentum ratio, unity indicates that all the momentum belongs to the mobile bed, whilst zero indicates that all the momentum is in the surge. Recall that lower values of $G_d$ imply more bed mobilization. We here chose this metric because it highlights the trends relating to the momentum more clearly. The time-step shown corresponds to quasi-steady-state behaviour. For $G_d > 0.5$ and low $\Gamma = h_b/h_s$ (i.e. shallow beds), almost all momentum is within the surge (Fig. 10a). For higher $\Gamma$, more momentum is in the bed for all $G_d$, implying that increasing $h_b$ for a given value of $G_d$ increases momentum transfer from surges to mobile beds.

Fig. 10b shows the ratio between the average momentum per SPH particle of the surge and mobile bed materials $\overline{U}_b/\overline{U}_s$ against $G_d$. Interestingly, the momentum per SPH point in the mobile bed is evidently somewhat insensitive to $\Gamma$. Whilst the bed thickness governs the total momentum in the bed (Fig. 10a), the per-point momentum in the bed increases only slightly as $\Gamma$ increases.

### 3.6. Mass exchange

Fig. 11 shows the percentage change in combined flow volume per second, i.e. the entrainment rate. To calculate this quantity, we defined a threshold for “moving” material, corresponding to $\overline{U} \ge 1$ mm/s. For regularized viscoplastic models, defining yielded/unyielded or moving/immobile material is somewhat arbitrary. Nevertheless, we checked that the threshold does not affect the general trends. We

average data over a 3 s duration, and quote the middle time-step on the figures. Comparing Fig. 11a and b shows that most entrainment occurs at lower $t$, when the surge front initially reaches the mobile bed. Either increasing $\Gamma$ (and hence $h_b$), or decreasing $G_d$, increases the entrainment rate. The maximum observed entrainment rate is around 4% of the total flow volume per second, occurring for the highest $\Gamma$ for a borderline-stable mobile bed. The entrainment rate thereafter (at $t = 7.5$ s) is between 0 and ~1%/s for all cases.

## 4. Discussion

### 4.1. Summary of regimes

We have identified several governing dimensionless numbers relating to surge/bed interaction, including $G_s$, $G_d$, and $\Gamma$. These depend on quantities including the bed depth $h_b$ and the bed yield stress $\tau_{y,b}$, which are the key variables explored in this study. Fig. 12 shows how the dimensionless numbers $\Gamma$, $G_s$ and $G_d$ relate to (a) the regimes we identified, as well as (b) the flow reach, (c) the total momentum in the bed and (d) the total volume of moving flow material. All data explored in this study form a 2D curved surface in this phase diagram. In fact, in principle, any two of the three dimensionless numbers can be used to identify a given surge-bed interaction regime, given that all of the points lie on a curved 2D surface. A single dimensionless number is insufficient to identify any particular regime. Notwithstanding, recall that the transitions between the regimes are smooth, and so the dimensionless numbers cannot absolutely pinpoint any particular interaction regime. The black zones indicate parameter combinations that could not be investigated. Specifically, the mobile bed could not be made shallower than 3 mm (i.e. ~ 7 point diameters) due to the finite resolution of the SPH. We used the red datapoints to interpolate
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Velocity profiles for (a) $h_b = 3$ mm, (b) $h_b = 12$ mm and (c) $h_b = 30$ mm, at time-steps (i) $t = 4$ and (ii) $t = 10$ s. The depth $z$ in the flow is shown against the local downstream velocity $U_s$. The increasing $G_s = \tau_{y,b}/\rho g h_b \sin \theta$ in each plot correspond to progressively less mobile beds. The light grey regions show the original depth of the mobile beds, whilst the dark grey regions show the outer boundaries. The insets show the depth of the bed, and the monitoring region (the shaded blue area).</figcaption></figure>

the surfaces in Fig. 12b to d, including some extra datapoints which lie beyond the limits of the axes (shown using solid red circles where possible).

Before describing each of the regimes in detail, a short discussion of the origin of the names adopted is warranted. Some of the names for the regimes observed were borrowed from previous works; others are newly coined. The “rigid” regime is simply where the yield stress of the bed is high enough that deformation within the bed remains negligibly small. The “surfing” regime is so called as it looks as though the surge surfs over the bed, whilst provoking limited deformation therewithin. The term “ploughing” has been borrowed from other papers (e.g. Sovilla et al., 2006), where a similar regime was observed for snow avalanches, with the surge pushing bed material downstream and upwards; “shallow ploughing” and “deep ploughing” correspond to such phenomena in shallow and deep beds, respectively. We use the terms “shallow” and “deep” in a qualitative way, as the transition between the regimes is continuous; for reference, however, “deep ploughing”

appreciably becomes more obvious when the bed is more than half of the surge depth. We adopted the term “lubrication” to highlight the role of the mobile bed in acting as an intermediary layer that reduced shear-induced dissipation (that would otherwise occur from direct interaction between the surge and the outer boundaries). Finally, the term “plunging” was adopted to highlight the strong slope-normal motion of the surge into the mobile bed zone.

1. The “rigid” regime (also referred to as “control regime” elsewhere in this manuscript), is observed for very high $G_s$ and $G_d$. No momentum is transferred from the surge into the bed, and there is no deceleration of the surge relating to bed deformation. The velocity profile includes a shear zone at the surge base, with a plug region above it, characteristic of free-surface viscoplastic surges.

2. The “surfing” regime is observed for intermediate $G_s$ and high $G_d$; surfing can occur for a wide range of $\Gamma$. The surge transfers
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Metrics related to bed momentum. (a): the total momentum in the bed $\rho \sum U_s$ normalized by the total momentum in the surge and the bed $\rho (\sum U_s + \sum U_b)$; (b): the average momentum of the bed $\rho \bar{U}_b$ normalized by that of the surge $\rho \bar{U}_s$, shown against $G_d = \tau_{y,b}/|\tau|$. Several values of $\Gamma = h_b/h_s$ are shown; increasing $\Gamma$ correspond to deeper beds. The inset shows that data is sampled from the entire computational domain.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Rate of change of flow volume; (a) $t = 2.5$ s; (b) $t = 7.5$ s. The rate of change of the moving combined flow material $\frac{V_t - V_l}{V_t}$, where “moving material” is defined as points where $|U| \ge 1$ mm/s is shown against $G_d = \tau_{y,b}/|\tau|$. Several $\Gamma = h_b/s$ are shown; increasing $\Gamma$ correspond to deeper beds.</figcaption></figure>

some momentum into the mobile bed, with the mobile bed correspondingly flowing slightly. The bed motion decelerates the flow, relative to the “rigid” regime. The velocity remains mostly in the surge, again with shear and plug regions.

3. The “lubrication” regime is characterized by low $\Gamma$ (i.e. low $h_b$), low $G_s$ and low $G_d$ (i.e. shallow, borderline-stable beds). All bed material is entrained by the surge. Once “kickstarted” by the arriving surges, the beds rapidly accelerate, partially due to their own weight, and partially through transfer of momentum from the surge. These weak mobile beds behave as lubricants, accelerating surges relative to the “rigid” regime. The combined flow has one shear and one plug zone, with the shear zone localized at the lower boundary.

4. The “shallow ploughing” regime is observed for intermediate $G_s$, low $G_d$ and low $\Gamma$. Surges push mobile bed material forward, forming a mobile dam that retards the progress of the combined flow front, relative to the “rigid” regime. Some momentum is transferred into the mobile beds, with most bed material entrained by the surge front. There is a single shear zone, which is characterized by an inflection point at the surge/bed interface.

5. The “deep ploughing” regime is seen for intermediate $G_s$ and low $G_d$, but high $\Gamma$. Deep mobile beds take more momentum from the surge than for the “shallow ploughing” regime. There are separate shear zones for the surge and bed, with a clear inflection point between the two. However, only the surge exhibits a plug region: the shearing from the surge causes the mobile bed material to deform across its entire depth.

6. Finally, the “plunging” regime is observed for low $G_s$, high $G_d$ and high $\Gamma$. Relatively weak mobile bed material enables the surge to plunge far into the bed, massively displacing the bed material upward and downstream. Large amounts of momentum are transferred into the bed from the surge. Furthermore, as the bed was originally borderline-stable, the bed gains extra momentum due to its own weight, after having been “kickstarted”. Two shear zones with plugs atop them exist for the surge and bed, respectively. There is a clear inflection point located deep in the bed zone, corresponding to the bed-surge interface, in the velocity profiles of the surge and the bed.

**Fig. 12b** superimposes a colour map of the combined flow reach at $t = 10.0$ s, when the flow has become well-developed. The flow reach is scaled such that the reach corresponding to the rigid bed falls in the exact center of the colorbar. However, two important items should be noted: (i) the point corresponding to the rigid bed modelled
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Phase diagrams showing how Γ, Gs and Gd relate. Recall that Γ = h<sub>b</sub>/h<sub>s</sub>, with higher Γ implying deeper beds; G<sub>d</sub> = τ<sub>y,b</sub>/|Γ|<sub>i</sub>, with higher G<sub>d</sub> implying higher bed yield stresses; and G<sub>s</sub> = τ<sub>y,b</sub>/ρgh<sub>b</sub> sinθ. The colorbars are scaled to the minimum and maximum values of each quantity. The red squares approximately demarcate the shallow/deep regimes. (a): A qualitative phase diagram indicating the following regimes: shallow/deep ploughing; lubrication; plunging; and surfing. The sole datapoint relating to the rigid regime is off the chart, but the dashed purple arrow indicates that the rigid regime corresponds to high G<sub>s</sub> and G<sub>d</sub>. Red circles indicate the datapoints used for interpolation in other parts of this Figure; as part of the interpolation process, some of the datapoints could not be included on the main plot. The solid red circles represent some of the datapoints that are “off the chart” as an indicative reference. (b): The flow reach (defined in Section 3). (c): Total momentum in the bed, normalized by the total momentum in the surge (ΣU<sub>b</sub>/ΣU<sub>s</sub>). (d): The total volume of moving flow material, incorporating a velocity threshold to ignore near-static material.</figcaption></figure>

in our simulations lies off the shown plot area, as the interpolation
process used to generate the plots in Fig. 12 cannot deal well with
edge cases; and (ii) certain other combinations of bed properties can
have the same reach as the rigid bed, hence the green points near the
right corner of Fig. 12(b). The “lubrication” and “plunging” regimes
involve the highest and lowest flow reaches, respectively. The “rigid”,
“surfing” and “ploughing” regions exhibit intermediate flow reaches.
The non-monotonous behaviour – with intermediate values of τ<sub>y,b</sub>
decelerating surges relative to immobile (rigid) beds, and low values
of τ<sub>y,b</sub> accelerating surges – must be highlighted. Fig. 12c shows the
total momentum in the bed, normalized by the momentum in the
surge. Where G<sub>s</sub> ≈ 1 (a borderline-stable bed), the bed gains the
most momentum relative to the surge, due to a combination of (i)
its self-weight and (ii) momentum transferred from the surge. For the
“rigid” and “surfing” regions, the bed momentum is far lower, as the
bed material is harder to mobilize by the surge; the ploughing regime
constitutes an intermediate case. Fig. 12d shows the combined volume
of *moving flow* material, normalized by the total volume of *surge and
mobile bed* material. Predictably, the “rigid-like” zone involves almost
no entrainment. The “plunging regime” sees the highest combined flow
volume. In the “lubrication regime”, even though all the bed material
is entrained, the volume is limited by h<sub>b</sub>, hence the intermediate values
for combined flow volume for this region.

Fig. 12 shows that the regimes form a continuous spectrum. It is
therefore difficult in some cases to clearly demarcate the transition
between regimes.

4.2. Comparison with other studies

Most research on entrainment processes has hitherto focused on dry
or wet granular surges (Gauer and Issler, 2004; Berger et al., 2011;
Iverson et al., 2011; Iverson, 2015; Nikooei and Manzari, 2020; Song
and Choi, 2021; Li et al., 2022); the complexity of these materials

<figure><img src="image_16.png" /><figcaption>Fig. A.13. Numerical setup for the sensitivity study. All dimensions are given in mm. There is no mobile bed, but a no-slip condition is imposed.</figcaption></figure>

complicates the systematic study of entrainment mechanisms. In con-
trast, considering viscoplastic fluids affords some key simplifications;
for instance, the yield stress is constant and is unambiguously de-
fined; viscoplastic fluids are incompressible; and stress transmission
involves no heterogeneous dispersion of stresses through a granular
skeleton. However, studying entrainment mechanisms using just this
simplified rheology, and only varying τ<sub>y,b</sub> and h<sub>b</sub>, reveals a wide range
of entrainment regimes. This highlights the fundamental complexity of
the problem. Notwithstanding, comparisons with findings from studies
considering other materials are warranted.

Gauer and Issler (2004) distinguish between “erosion” for the break-
ing up of a bed (substrate), and “entrainment” for this material be-
coming part of the surge. Li et al. (2022) built on this work, providing
fascinating insights into the distinctions between the two behaviours
<figure><img src="image_17.png" /><figcaption>Fig. A.14. Effects of SPH point resolution on tensile instabilities. (a) $dp = 0.25$ mm; (b) $dp = 1.00$ mm. (i) $t = 0.4$ s; (ii) $t = 0.8$ s. For both cases, $c_0 = 100$ m/s and $\mu_{reg} = 1$ kPa s.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. A.15. Effects of varying the point size $dp$ on the velocity profile: normalized velocity versus normalized flow depth. Note: the surge is shallow and unsteady, so the plug region is not well-developed.</figcaption></figure>

for snow avalanches. A key factor for the occurrence of “erosion” is the breakup of cohesion in the granular bed (i.e. a change in the rheological behaviour of the bed), due to shearing action from the surge. In our simplified model using viscoplastic fluids, “erosion” could not be directly observed, so only entrainment was investigated.

Of the various entrainment regimes given in Gauer and Issler (2004), only “ploughing” was directly observed in our flows, although “abrasion”-like behaviour occurred for the “surfing” regime. For the dry granular column collapses modelled by Nikooei and Manzari (2020), a “ploughing” regime was also observed when the friction angle of the mobile bed was sufficiently low. Iverson et al. (2011) observed behaviour that we term a “lubrication” regime, where surges accelerate. Their surges comprised gravel, sand and mud, traversing shallow mobile beds. “lubrication” was observed for saturated beds, where the effective stress (and hence the yield stress) was lowest; surges were faster than on bare beds. We have shown that acceleration can occur for deeper beds as well, via “plunging”.

At the Illgraben catchment in Switzerland, Berger et al. (2011) observed that most entrainment occurs at the surge front of debris

flows; Song and Choi (2021) made similar findings using small-scale physical tests, as did (Nikooei and Manzari, 2021) using numerically-modelled dry granular surges. Snow avalanches exhibit similar behaviour (Sovilla and Bartelt, 2002; Gauer and Issler, 2004; Eglit and Demidov, 2005). Our results for viscoplastic flows are consistent in this regard, with most mass transfer from the bed to the flow occurring when the surge front arrives.

## 5. Concluding comments

In this paper, we used an SPH model with a Herschel-Bulkley rheology to investigate the interaction between gravity-driven surges and mobile beds, focusing on the entrainment mechanisms.

1. We compared the SPH model output with physical test data from Bates and Ancey (2017), finding a good match between the geometry of the flow in their physical and our numerical results.

2. We systematically varied the bed depth $h_b$ and the bed yield stress $\tau_{y,b}$. We identified several governing dimensionless numbers relating to surge/bed interaction, including $G_s$, $G_d$, and $\Gamma$, allowing us to rationalize several distinct interaction regimes, some of which have previously been observed by other authors using different materials. These regimes are: (i) “rigid” bed; (ii) “surfing”; (iii) “lubrication”; (iv) “shallow ploughing”; (v) “deep ploughing”; (vi) “plunging”

3. There exist regimes in which the combined flow can be accelerated. Most notably, this includes the “lubrication” regime, in which shallow beds with low yield stresses are massively destabilized by the surge.

4. The behaviour threading the six different regimes together is highly non-monotonous. High values of $G_s$ (i.e. rigid beds) give intermediate surge velocities. In comparison to this case, low values of $G_s$ accelerate the surge, whilst intermediate values of $G_s$ decelerate the surges. This non-monotonous behaviour was unexpected and underscores the complexity of the surge-bed interaction mechanisms.

The results from this study underscore the difficulty of capturing the mechanisms of interaction between surges and mobile beds in a constitutive framework. Notwithstanding, our results also provide a basis for further investigations on entrainment, focusing on specific interaction regimes. Indeed, it would be worth testing how the regimes observed in this study manifest for broader parameter combinations. A first step
<figure><img src="image_19.png" /><figcaption>**Fig. A.16.** Comparison of the flow at (a) *t* = 0.15 s, (b) *t* = 0.30 s, and (c) *t* = 0.60 s, for different values of the viscosity cap *μ*<sub>reg</sub>: (i) 1, (ii) 10, (iii) 100, (iv) 1000, and (v) 10,000 Pa s. For these simulations the point size *dp* is 0.4 mm, whilst the speed of sound *c*<sub>0</sub> is 100 m/s.</figcaption></figure>

would be to investigate the effects of varying the slope angle and the Froude number of the incoming surge. The effects of the bed length and the surge rheology also warrant further investigation. Furthermore, future work studying the applicability of the different entrainment regimes identified in this study for field applications, including those in the vicinity of structures, would be valuable. Such future studies would profit from using improved versions of the software package adopted in this study, specifically with regards to enforcing computational stability (to improve the “remeshing” method developed in our study), and the boundary conditions (to avoid the limitations of the time-steps required for our no-slip boundaries), which would enable full 3D modelling.

CRediT authorship contribution statement

**Saoirse Robin Goodwin:** Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – review & editing, Visualization, Project administration, Funding acquisition. **Suzanne Lapillonne:** Software, Writing – review & editing. **Guillaume Piton:** Writing – review & editing, Supervision, Funding acquisition. **Guillaume Chambon:** Conceptualization, Methodology, Formal analysis, Writing – review & editing, Supervision, Project administration, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

This work was funded by “MOPGA 2021–22” (MOPGA-976501H, ANR, French Republic).

<table><caption>**Table A.1**
Parameters varied across the sensitivity study.</caption><thead><tr><td>Speed of sound c<sub>0</sub> (m/s)</td><td>Viscosity cap μ<sub>reg</sub> (kPa s)</td><td>Point size dp (mm)</td></tr></thead><tbody><tr><td>1, 10, 100</td><td>0.01, 0.1, 1, 10, 100</td><td>0.25, 0.40, 0.50, 1.00</td></tr></tbody></table>

Computer code availability

We modified an open-source version of DualSPHysics, specifically “DualSPHysics5_NNewtonian v1.005 (v5.0.164) (21-11-2020)”. This code is available for download from the DualSPHysics project. The seven modified files from our project — which belong in the folder “DSPH_v5.0_NNewtonian/source/” can be found on GitHub at: https://github.com/srgoodwin/dualsph_hb/.

Appendix A. Sensitivity study

The sensitivity study was performed with a smaller volume of material than the main study, as certain parameter combinations were computationally expensive. An effectively “immobile” rigid bed, characterized by a very large yield stress (see Section 2.3) was considered (Fig. A.13). The computational domain was 0.6 m long, and the channel slope was 16°. Each simulation lasted 0.8 s, covering the initial dam-break and the transition to a “creeping” regime. Velocity profiles inside the surges were measured in a sampling region with a downstream length of 10 mm, placed 250 mm downstream from the gate; we sampled data 0.4 s after dam-break, corresponding to the surge body. The influences of three numerical parameters, namely the speed of sound *c*<sub>0</sub>, the viscosity cap *μ*<sub>reg</sub>, and the SPH point size *dp* (spatial resolution), were systematically investigated. Table A.1 summarizes the simulation plan. Note that these parameters were applied to both the surge and the bed material.

A.1. Spatial resolution

The SPH method is prone to tensile instabilities (Monaghan, 2000), notably for transient free-surface flows (Xu and Yu, 2018). Artificial shifting techniques can alleviate this issue (Xu and Yu, 2018), whilst the
<figure><img src="image_20.png" /><figcaption>Fig. A.17. Velocity profiles at t = 0.4 s, for a range of c₀ and μₜreg values. Parts (a), (b) and (c) show data for values of c₀ of 1, 10 and 100 m/s respectively. Part (d) compares the “converged” profiles obtained for c₀ ≥ 10 m/s and μₜreg ≥ 0.1 kPa s.</figcaption></figure>

spatial resolution also plays an important role. We studied the effects of SPH point size dp, finding that coarser resolutions correlate with larger tensile instabilities (Fig. A.14). However, we observed that the resolution only minimally affects the bulk flow dynamics, reflected by the velocity profiles (Fig. A.15). We thus set dp = 0.4 mm, balancing stability and computational speed.

## A.2. Speed of sound c₀ and viscosity cap μₜreg

The SPH timestep Δt is governed by the following Courant-Friedrichs-Lewy condition:

$$ \Delta t \le A_{CFL} \cdot \min \left[ \frac{h}{c_0}, \frac{\rho_m h^2}{\mu_m}, \sqrt{\frac{h}{a_m}}, \kappa \frac{\rho_0}{d \rho_m / dt} \right] \quad (A.1) $$

where A_CFL is a “factor of safety” (value: 0.2); ρ_m is the density of point m; μ_m is its effective viscosity; a_m is its acceleration; κ is a constant (value: 10⁻³); and dρ_m/dt is the particle density variation. Both the speed of sound c₀ and the viscosity cap μₜreg directly affect the value

of Δt; setting their values requires compromising between accuracy, stability and computational time.

The speed of sound c₀ should be chosen so the Mach number M = u/c₀ remains between 0.001 and 0.1 (Monaghan, 1999; Liu and Liu, 2003; Violeau and Leroy, 2014). The velocities in our simulations were ≤ 2 m/s, implying c₀ > 20 m/s. However, our simulations involve two characteristic velocities — during (i) initial dam-break, and (ii) the slower propagation phase. To reduce the computational cost for the propagation phase, we chose the lowest value of c₀ possible, namely c₀ = 25 m/s. Note that this value of c₀ maintains the density ρ within ±3 % of ρ₀ (Morris et al., 1997, see).

The viscosity cap μₜreg is set by considering: (i) the yield stress and viscosity of the material; (ii) the typical magnitude of stresses and shear rates expected; and (iii) the timescale of the flow. Fig. A.16 shows that small values of μₜreg cause the surge and bed to behave essentially like Newtonian fluids (Fig. A.16(i)). Conversely, viscoplastic behaviour, evidenced by an effectively immobile bed and the surge retaining some initial geometric features during deformation, is recovered for large
<figure><img src="image_21.png" /><figcaption>**Fig. A.A.18.** Comparison of the time histories of (a) the flow depth in the measuring section and (b) position of the flow front for varying consistencies *k* and yield stresses *τ* around the nominal values considered in the main study. Values were selected to cover the measurement uncertainties from Bates and Ancey (2017).</figcaption></figure>

enough values of $\mu_{reg}$. Qualitatively, convergence is observed for $\mu_{reg} \ge 100$ Pa s (Figs. A.16(iii) to (v)).

Fig. A.17 shows the velocity profiles inside the surges for different values of $c_0$ and $\mu_{reg}$. For $\mu_{reg} = 1$ Pa s, the velocity profile descends into the bed (see also Fig. A.16). For $\mu_{reg} = 10$ Pa s, the no-slip condition at the bed is still not obtained. By contrast, the velocity profiles for $\mu_{reg} = 0.1$ kPa s and $\mu_{reg} = 0.5$ kPa s exhibit both no-slip and a plug region, and are almost identical. More specifically, Fig. A.17(d) shows that the velocity profiles obtained for $\mu_{reg} \ge 0.1$ Pa s and $c_0 = 10$ or 100 m/s are all extremely similar. The slight increase in velocity as both $c_0$ and $\mu_{reg}$ are increased is due to the occurrence of tensile instabilities (Fig. A.14), which tend to be amplified for stiffer materials. In the main simulations, these tensile instabilities were alleviated using the “regeneration” method detailed in Section 2.3. Finally, we adopted a value $\mu_{reg} = 10$ kPa s for these main simulations.

## Appendix B. Simulation plan

See Table B.2.

## Appendix C. Effects of rheological parameters $\tau_c$ and $k$ on flow velocity & depth

Rheometrical measurements of the yield stress $\tau_y$ and the consistency $k$ can be affected by relatively large uncertainties (see Fig. 3a). For the same configuration as that considered in Appendix A, Fig. A.A.18 shows how variations of these parameters by ±20% affect (a) the flow

<table><caption>**Table B.2**
Simulations performed in our study.</caption><thead><tr><th>Test series</th><th>Bed depth h<sub>b</sub> (mm)</th><th>Bed yield stress τ<sub>y,b</sub> (Pa)</th></tr></thead><tbody><tr><td>Immobile bed</td><td>0</td><td>1000</td></tr><tr><td rowspan="10">Mobile bed</td><td>3</td><td>18, 28, 38, 48, 58, 68, 78, 88, 118, 178, 218</td></tr><tr><td>4</td><td>38, 48, 58, 68, 78, 88, 98</td></tr><tr><td>5</td><td>28, 38, 48, 58, 68, 78, 88</td></tr><tr><td>6</td><td>18, 38, 48, 58, 78, 98, 118, 158, 218</td></tr><tr><td>9</td><td>58, 78, 98, 118, 138, 158, 178, 198, 238</td></tr><tr><td>12</td><td>38, 58, 78, 98, 118, 138, 158, 178, 198, 218</td></tr><tr><td>18</td><td>58, 68, 78, 98, 118, 138, 158, 198, 218</td></tr><tr><td>24</td><td>78, 88, 98, 118, 138, 158</td></tr><tr><td>30</td><td>58, 78, 98, 118, 158, 188, 218, 268</td></tr><tr><td>42</td><td>138, 178, 218, 258</td></tr></tbody></table>

<table><caption>**Table D.3**
Numerical SPH settings. “Num” = numerical parameter; “phys” = physical parameter.</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Gravity (phys)</td><td>g</td><td>9.81</td><td>m/s</td></tr><tr><td>Const. for state equation (num)</td><td>γ</td><td>7</td><td>-</td></tr><tr><td>Coeff. for speed of sound (num)</td><td>k<sub>sound</sub></td><td>20</td><td>-</td></tr><tr><td>Speed of sound (num)</td><td>c<sub>0</sub></td><td>25</td><td>m/s</td></tr><tr><td>Nominal point spacing (num)</td><td>dp</td><td>0.0005</td><td>m</td></tr><tr><td>Kernel radius (num)</td><td>h</td><td>0.0014</td><td>m</td></tr><tr><td>Coeff. for CFL (num)</td><td>A<sub>CFL</sub></td><td>0.2</td><td>-</td></tr><tr><td>Artificial viscosity factor (num)</td><td>α</td><td>0</td><td>-</td></tr><tr><td>Artificial viscosity factor at bounds (num)</td><td>-</td><td>0</td><td>-</td></tr><tr><td>Minimum density possible (num)</td><td>ρ<sub>min</sub></td><td>500</td><td>kg/m<sup>3</sup></td></tr><tr><td>Maximum density possible (num)</td><td>ρ<sub>max</sub></td><td>3000</td><td>kg/m<sup>3</sup></td></tr></tbody></table>

<table><caption>**Table D.4**
Global settings for the program.</caption><thead><tr><th>Setting name</th><th>Setting</th><th>Notes</th></tr></thead><tbody><tr><td>Step algorithm</td><td>Symplectic</td><td>-</td></tr><tr><td>Kernel</td><td>Wendland</td><td>Default</td></tr><tr><td>Rheology treatment</td><td>Single &amp; multiphase</td><td>-</td></tr><tr><td>Velocity gradient type</td><td>SPH</td><td>-</td></tr><tr><td>Viscosity treatment</td><td>Constitutive equation</td><td>-</td></tr><tr><td>Density diffusion term</td><td>Fourtakas (full)</td><td>-</td></tr><tr><td>Density diffusion term value</td><td>0.1</td><td>Default</td></tr><tr><td>Shifting algorithm</td><td>Full</td><td>-</td></tr><tr><td>Shifting coefficient</td><td>-10</td><td>-</td></tr><tr><td>Shifting threshold for free surface</td><td>1.5</td><td>Default</td></tr><tr><td>Rigid algorithm</td><td>SPH</td><td>Default</td></tr><tr><td>Floating pause at start of simulation</td><td>0 s</td><td>Default</td></tr><tr><td>Coefficient to find miniumum timestep</td><td>0.05</td><td>Default</td></tr><tr><td>Relaxation timestep parameter for visc.</td><td>0.2</td><td>Default</td></tr><tr><td>Fixed timestep</td><td>Off</td><td>Default</td></tr><tr><td>Simulation time</td><td>6 s</td><td>-</td></tr><tr><td>Timestep output</td><td>0.01 s</td><td>-</td></tr></tbody></table>

depth measured in the measuring region, and (b) the position of the flow front relative to the top of the channel. As seen, increasing either $k$ or $\tau_y$ increases surge body thickness and reduces the surge propagation velocity. These results provide context for the discrepancies between physical and numerical measurements shown in Section 3.1.

## Appendix D. Numerical settings

Table D.3 summarizes the main numerical parameters adopted in this study. Additional DualSPHysics settings are indicated in Table D.4. We selected the input parameters following a wide-ranging sensitivity study (Appendix A).

## Appendix E. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.cageo.2023.105476.
## References

Ancey, C., 2007. Plasticity and geophysical flows: A review. J. Non-Newton. Fluid Mech. 142 (1), 4-35. http://dx.doi.org/10.1016/j.jnnfm.2006.05.005, Viscoplastic fluids: From theory to application.

Balmforth, N.J., Frigaard, I.A., Ovarlez, G., 2014. Yielding to stress: Recent developments in viscoplastic fluid mechanics. Annu. Rev. Fluid Mech. 46 (1), 121-146. http://dx.doi.org/10.1146/annurev-fluid-010313-141424.

Bartzke, G., Fourtakas, G., Canelas, R., Rogers, B.D., Huhn, K., 2021. Simulation of flow past a sphere on a rough bed using smoothed particle hydrodynamics (SPH). Comput. Part. Mech. 1-14. http://dx.doi.org/10.1007/s40571-021-00417-x.

Bates, B.M., 2015. Fluid Dynamics of Basal Entrainment by Geophysical Gravity-Driven Flows (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, URL: https://doi.org/10.5075/epfl-thesis-6765.

Bates, B.M., Ancey, C., 2017. The dam-break problem for eroding viscoplastic fluids. J. Non-Newton. Fluid Mech. 243, 64-78. http://dx.doi.org/10.1016/j.jnnfm.2017.01.009.

Berger, C., McArdell, B.W., Schlunegger, F., 2011. Direct measurement of channel erosion by debris flows, illgraben, Switzerland. J. Geophys. Res. 116, F01002. http://dx.doi.org/10.1029/2010JF001722.

Chambon, G., Bouvarel, R., Laigle, D., Naaim, M., 2011. Numerical simulations of granular free-surface flows using smoothed particle hydrodynamics. J. Non-Newton. Fluid Mech. 166 (12), 698-712. http://dx.doi.org/10.1016/j.jnnfm.2011.03.007.

Cousso, P., 1994. Steady, laminar, flow of concentrated mud suspensions in open channel. J. Hydraul. Res. 34, 535-559. http://dx.doi.org/10.1080/00221686.1994.9728354.

Cousso, P., Laigle, D., Arattano, M., Deganutti, A., Marchi, L., 1998. Direct determination of rheological characteristics of debris flow. p. 124. http://dx.doi.org/10.1061/(ASCE)0733-9429(1998)124%3A8(865).

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., Gómez-Gesteira, M., 2021. DualSPHysics: from fluid dynamics to multiphysics problems. Comput. Part. Mech. http://dx.doi.org/10.1007/s40571-021-00404-2.

Eglit, M., Demidov, K., 2005. Mathematical modeling of snow entrainment in avalanche motion. Cold Reg. Sci. & Technol. 43 (1), 10-23. http://dx.doi.org/10.1016/j.coldregions.2005.03.005.

Fourtakas, G., Rogers, B.D., 2016. Modelling multi-phase liquid-sediment scour and resuspension induced by rapid flows using smoothed particle hydrodynamics (SPH) accelerated with a graphics processing unit (GPU). Adv. Water Resour. 92, 186-199. http://dx.doi.org/10.1016/j.advwatres.2016.04.009.

Fourtakas, G., Rogers, B.D., Laurence, D.R., 2013. Modelling sediment resuspension in industrial tanks using SPH. Houille Blanche 99 (2), 39-45. http://dx.doi.org/10.1051/lhb/2013014.

Frigaard, I.A., Nouar, C., 2005. On the usage of viscosity regularisation methods for visco-plastic fluid flow computation. J. Non-Newton. Fluid Mech. 127 (1), 1-26. doi:j.jnnfm.2005.01.003.

Gauer, P., Issler, D., 2004. Possible erosion mechanisms in snow avalanches. Ann. Glaciol. 38, 384-392. http://dx.doi.org/10.3189/172756404781815068.

Ghaïtanellis, A., Violeau, D., Ferrand, M., Abderrezzak, K.E.K., Leroy, A., Joly, A., 2018. An SPH elastic-viscoplastic model for granular flows and bed-load transport. Adv. Water Resour. 111, 156-173. http://dx.doi.org/10.1016/j.advwatres.2017.11.007.

Herschel, W.H., Bulkley, R., 1926. Konsistenzmessungen von Gummi-Benzollösungen. Kolloid Z. 39, 291-300. http://dx.doi.org/10.1007/BF01432034.

Huang, X., García, M.H., 1998. A Herschel-Bulkley model for mud flow down a slope. J. Fluid Mech. 374, 305-333. http://dx.doi.org/10.1017/S0022112098002845.

Issler, D., 1998. Modelling of snow entrainment and deposition in powder-snow avalanches. Ann. Glaciol. 26, 253-258.

Iverson, R.M., 1997. The physics of debris flows. Rev. Geophys. 35 (3), 245-296. http://dx.doi.org/10.1029/97RG00426.

Iverson, R.M., 2012. Elementary theory of bed-sediment entrainment by debris flows and avalanches. J. Geophys. Res. 117 (F03006), http://dx.doi.org/10.1029/2011JF002189.

Iverson, R.M., 2015. Scaling and design of landslide and debris-flow experiments. Geomorphology 244, 9-20. http://dx.doi.org/10.1016/j geomorph.2015.02.033.

Iverson, R.M., Ouyang, C., 2014. Entrainment of bed material by earth-surface mass flows: Review and reformulation of depth-integrated theory. Rev. Geophys. (AGU) 53, 27-58. http://dx.doi.org/10.1002/2013RG000447.

Iverson, R.M., Reid, M.E., Logan, M., LaHusen, R., Godt, J.W., Griswold, J.P., 2011. Positive feedback and momentum growth during debris-flow entrainment of wet bed sediment. Nat. Geosci. 4, 116-121. http://dx.doi.org/10.1038/NGE01040.

Jiang, T., Ouyang, J., Li, Q., Ren, J., Yang, B., 2011. A corrected smoothed particle hydrodynamics method for solving transient viscoelastic fluid flows. Appl. Math. Model. 35, 3833-3853. http://dx.doi.org/10.1016/j.apm.2011.02.014.

Labbé, M., Laigle, D., 2013. SPH modelling of the interaction between free-surface flows of a viscoplastic fluid and a structure. WIT Trans. Eng. Sci. 79, 399-408. http://dx.doi.org/10.2495/MPF130331.

Laigle, D., Lachampe, P., Naaim, M., 2007. SPH-based numerical investigation of mudflow and other complex fluid flow interactions with structures. Comput. Geosci. 11, 297-306. http://dx.doi.org/10.1007/s10596-007-9053-y.

Li, X., Sovilla, B., Ligneau, C., Jiang, C., Gaume, J., 2022. Different erosion and entrainment mechanisms in snow avalanches. J. Mech. Res. Commun. 124, 103914. http://dx.doi.org/10.1016/j.mechrescom.2022.103914.

Liu, G.-R., Liu, M.B., 2003. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific Publishing.

Macià, F., Antuono, M., Colagrossi, A., 2011. Benefits of using a wendland kernel for free-surface flows. In: Proceedings of the 6th International SPHERIC Workshop. pp. 30-37.

Mangeney, A., Tsimring, L.S., Vollson, D., Aranson, I.S., Bouchut, F., 2007. Avalanche mobility induced by the presence of an erodable bed and associated entrainment. Geophys. Res. Lett. 34, L22401. http://dx.doi.org/10.1029/2007GL031348.

McCoy, S.W., Kean, J.W., Coe, J.A., Tucker, G.E., Staley, D.M., Wasklewicz, T.A., 2012. Sediment entrainment by debris flows: In situ measurements from the headwaters of a steep catchment. J. Geophys. Res.: Earth Surf. 117, F03016. http://dx.doi.org/10.1029/2011JF002278.

McDougall, S., Hungr, O., 2005. Dynamic modelling of entrainment in rapid landslides. Can. Geotech. J. 42, 1437-1448. http://dx.doi.org/10.1139/t05-064.

Monaghan, J.J., 1988. An introduction to SPH. Comput. Phys. Comm. 48, 89-96.

Monaghan, J.J., 1999. Solitary waves on a cretan beach. J. Waterw. Port Coast. Ocean Eng. (ASCE) 125 (3), 145-154.

Monaghan, J.J., 2000. SPH without a tensile instability. J. Comput. Phys. 159, 290-311.

Morris, J.P., Fox, P.J., Zhu, Y., 1997. Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys. 136 (1), 214-226. http://dx.doi.org/10.1006/jcph.1997.5776.

Nikooei, M., Manzari, M.T., 2020. Studying effect of entrainment on dynamics of debris flows using numerical simulation. Comput. Geosci. 134 (10), 104337. http://dx.doi.org/10.1016/j.cageo.2019.104337.

Nikooei, M., Manzari, M.T., 2021. Investigating the effect of mixing layer rheology on granular flow over entrainable beds using SPH method. Comput. Geosci. 115, 104792. http://dx.doi.org/10.1016/j.cageo.2021.104792.

Pastor, M., Stickle, M.M., Dutto, P., Mira, P., Merodo, J.A.F., Blanc, T., Sancho, S., Benitez, A.S., 2015. A viscoplastic approach to the behaviour of fluidized geomaterials with application to fast landslides. Contin. Mech. Thermodyn. 27 (1-2), 21-47.

Pirulli, M., Pastor, M., 2012. Numerical study on the entrainment of bed material into rapid landslides. Géotechnique 62 (11), 959-972. http://dx.doi.org/10.1680/geot.10.P.074.

Putz, A., Frigaard, I.A., Martinez, D.M., 2009. On the lubrication paradox and the use of regularisation methods for lubrication flows. J. Non-Newton. Fluid Mech. 163 (1-3), 62-77. http://dx.doi.org/10.1016/j.jnnfm.2009.06.006.

Reece, G., Rogers, B.D., Lind, S., Fourtakas, G., 2020. New instability and mixing simulations using SPH and a novel mixing measure. J. Hydrodyn. 32, 684-698. http://dx.doi.org/10.1007/s42241-020-0045-x.

Saramito, P., Wachs, A., 2016. Progress in numerical simulation of yield stress fluid flows. Rheol. Acta 56 (3), 1-20. http://dx.doi.org/10.1007/s00397-016-0985-9.

Skillen, A., Lind, S., Stansby, P.K., Rogers, B.D., 2013. Incompressible smoothed particle hydrodynamics (SPH) with reduced temporal noise and generalised fickian smoothing applied to body-water slam and efficient wave-body interaction. Comput. Methods Appl. Mech. Engrg. 265, 163-173. http://dx.doi.org/10.1016/j.cma.2013.05.017.

Song, P., Choi, C.E., 2021. Revealing the importance of capillary and collisional stresses on soil bed erosion induced by debris flows. J. Geophys. Res.: Earth Surf. 1-16. http://dx.doi.org/10.1029/2020JF005930.

Sovilla, B., Bartelt, P., 2002. Observations and modelling of snow avalanche entrainment. Nat. Hazards Earth Syst. Sci. 2 (3/4), 169-179. http://dx.doi.org/10.5194/nhess-2-169-2002.

Sovilla, B., Burlando, P., Bartelt, P., 2006. Field experiments and numerical modeling of mass entrainment in snow avalanches. J. Geophys. Res. 111, F03007. http://dx.doi.org/10.1029/2005JF000391.

Violeau, D., Leroy, A., 2014. On the maximum time step in weakly compressible SPH. J. Comput. Phys. 256 (1), 388-415. http://dx.doi.org/10.1016/j.jcp.2013.09.001.

Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math. 4 (1), 389-396.

Xu, X., Yu, P., 2018. A technique to remove the tensile instability in weakly compressible SPH. Comput. Mech. 62, 963-990. http://dx.doi.org/10.1007/s00466-018-1542-4.

Zubeldia, E.H., Fourtakas, G., Rogers, B.D., Farias, M.M., 2018. Multi-phase SPH model for simulation of erosion and scouring by means of the shields and Drucker-Prager criteria. Adv. Water Resour. 117, 98-114. http://dx.doi.org/10.1016/j.advwatres.2018.04.011.