E3S Web of Conferences 415, 02018 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502018

Chute flows of dry granular media: numerical simulations by awell-posed multilayer model and comparisons with experiments
Luca Sarno1,2*, Yongqi Wang1, Yih-Chin Tai3, Maria Nicolina Papa2, and Paolo Villani2
1Chair of Fluid Dynamics, Dept. of Mechanical Engineering, Technical University of Darmstadt, 64287 Darmstadt, Germany
2Laboratory of Environmental and Maritime Hydraulics, Dept. of Civil Engineering, University of Salerno, 84084 Fisciano, Italy
3Dept. of Hydraulic and Ocean Engineering, National Cheng-Kung University, 70101 Tainan, Taiwan

Abstract. Debris flows and avalanches are dangerous natural phenomena, characterized by the gravitydriven motion of granular media immersed in a fluid. For an appropriate hazard assessment or disastermitigation by scenario investigation, it is crucial to capture the underlying dynamics of the granular solid
phase. For this purpose, a multilayer depth-averaged approach represents a promising and computationallyefficient tool over fully three-dimensional models. Here we use a mathematically well-posed multilayer
model, which implements the µ(I)-rheology and a dilatancy law depending on the inertial number, I, andcompare the numerical results of the model with laboratory experiments of steady uniform chute flows overan erodible bed. The well-posedness of the model for any value of I, which is essential to get convergent
numerical solutions, is achieved by considering an approximation of the in-plane stress gradients, directly
emerging from the µ(I)-rheology. The predicted velocity profiles show a very good agreement with theexperimental ones, measured by particle image velocimetry (PIV). The volume fraction profiles by themultilayer model are also in good qualitative agreement with those measured by using the stochastic-opticalmethod (SOM), while they tend to overestimate the volume fraction measurements in the more dilute upper
region, closer to the free surface.

1 Introduction
Debris flows and granular avalanches are dangerous
phenomena for humans and infrastructures inmountainous regions. The physical understanding of the
granular dynamics and the related mathematicalmodelling are essential for the risk assessment of such
events. In the context of continuum mechanics, depthaveraged models, equipped with suitable constitutivelaws, are widely employed for the description of
granular flows, especially for large field-scale
simulations (e.g., [1]), as they are computationally more
efficient than three-dimensional (3D) models. Yet,depth-averaged models miss significant dynamicalprocesses occurring along the transverse flow direction.
While single-layer models are suitable for describing
granular flows with a blunt velocity profile, they
become insufficient if the velocity profiles are highly
sheared and swiftly change in shape owing to arheological stratification (e.g., [2]). The multilayer
approach [3-5] represents a way to partially overcome
the aforementioned limitations without renouncing thecomputational advantages of the depth-averaged
models. In this work, we employ a mathematically wellposed multilayer model, introduced by [5], and compare
the modelled velocity and volume fraction profiles withthe experimental data of steady uniform chute granular
flows over an erodible bed. The model uses the wellestablished µ(I)-rheology [6] and also incorporates adilatancy law, depending on the inertial number I, totake into account the effects of the variable volume
fraction [7]. It is worth noting that the well-posedness ofthe model for any value of the inertial number, I, is also
practically useful to obtain convergent numericalsolutions [5]. For a suitable comparison with theexperiments in a narrow chute, the multilayer model ishere equipped with rate-independent Coulomb-friction
terms accounting for the sidewall resistances (e.g., [8]).
Fig. 1 shows a sketch of the multilayer approach, where
the spatial domain of the granular avalanche ispartitioned into a pre-set number, N, of layers thatdynamically evolve according to the chosen rheology.

Fig. 1. Sketch of the multilayer approach [5], where thegeneric layer is represented by the symbol, i , the layer
thickness by hi , and the layer interfaces by i 1/ 2 .

*Correspondingauthor: sarno@fdy.tu-darmstadt.de© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 02018 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502018

2 Multilayer model equations
numerical integrationand  I    min 1 h    h u  M

t

i

ix

iix ,ix ,i 1/ 2

diameter, is a diffusion-like coefficient accounting forthe effects of the in-plane stress gradients, and ensures
the well-posedness of the model for all I [5]. The shear
stresses, t xz*, in (3) are calculated by suitablei 1 2

discretizations of the µ(I)-rheology [5],  s 2  s  2
pI0 I  1 *T D.
D(1) M x ,i 1/ 2  ,(2)is the static friction coefficient,  2 is an upper bound for 

 , and I 0 is a constant [6].
The PDE system, (1)-(3), is numerically integrated bymeans of an operator splitting scheme, where itsadvection part is advanced by a finite volume scheme,
employing a lateralized Harten–Lax–van Leer (LHLL)
approximate Riemann solver [10]. High resolution in
space is achieved thanks to a MUSCL predictorcorrector treatment. More details on the numerical
scheme can be found in [5].

 h uh i i x ,i    gz i i  t hi i u x ,i   x  h x 2 





i ii



h
 i 1 h j j  N  g z  i i  x  j i 1 j h j  hi i  x   j 1  ij








2



g x hi i  g z hi i  x b h  h u N  x  i hi  j i 1 j h j  i i  x  i i x ,i   h 2
 i i 
1 *

t xz
 t xz*

g



1g

2



i 1 2i 1 2

(5)

In (5) T* is the shear-stress tensor, D is the strain-ratetensor,
is the Frobenius norm, p is the pressure,  sg2(4)

projected along x, which can be set equal to zero for thecase of immiscible layers or calculated by suitableclosures. Finally, Eq. (3) is the momentum balance
equationofthelayeri.
In(3),  d g  s  2  g z I 0 0 , with d g being the grainand, thus, is usually neglected in the first-order depthaveraged models available in the literature [9]. Withreference to the generic layer i, the model partial
differential equations (PDE) are relax t i  u  x i   i  i , EOS ,.

1 a ' I
In (4) a’ is an empirical parameter,  min and  max are,respectively, the lower and upper bounds for  [5]. Eq.
(2) is the mass balance of the generic layer i, where
M x ,i 1/ 2 are the mass fluxes per unit length of interface

The model equations were obtained by depth-averaging
the mass and momentum balances in each layer,
followed by first-order asymptotic expansion for theparameter   H L  1 , with H and L, respectively,
being the typical depth and length scales of theavalanche. More details on the derivation of the modelequations and its key mathematical properties can befound in [5]. A distinctive feature of the proposed
multilayer model is that an approximation of the inplane stress gradients, directly emerging from the µ(I)rheology, is considered in the layer-averaged
momentum equations, to ensure the well-posedness for
any I. It is worth noting that this additional term issmaller than the leading order approximation, O    ,*i max  min(3)3 Comparisons with experimental dataand discussion

The predictions by the multilayer model are comparedwith the steady uniform velocity and volume fraction
profiles of chute flow experiments over an erodible bed,
carried out at the University of Salerno (Italy).

 ui 1/ 2,downwind M x,i 1/ 2  ui 1/ 2,downwind M x,i 1/ 2  N h  h u h
g z   j  i 1 j h j  i i  i i x ,i i i ,W2  hi i u x ,i hi i side

where g x  g sin  , g z  g cos  ,  g is the grain
density, W is the channel width and side is a friction
coefficient parameterizing the sidewall resistances. In(1)-(3) the unknowns are the layer-averaged volume
fractions,  i , the layer masses, hi i , and the layeraveraged momenta, hi i ux ,i , with u x ,i denoting the xcomponent velocity in the layer i. Eq. (1) describes thetransport of  i and its concurrent relaxation toward theequilibrium value,  i , EOS , according to a relaxation
time,  relax [5]. The dilatancy law, here used for defining

Fig. 2. The experimental set up.

 i , EOS , represents a slight modification of the formula

The laboratory setup, similar to that used in [11], is
composed of a 2-m-long Plexiglas chute with a width,

proposed by [7] for dense granular flows,2E3S Web of Conferences 415, 02018 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502018

W  80mm , and an inclination of 32.5° (Fig. 2). A 6cm-high weir was placed at the chute outlet to enforce
the formation of a static wedge of granular material
below the active surface flow. The granular material,
composed of polyoxymethylene (POM) spheroidal
beads with diameter, d g  3.3mm , was initially stored

simulations were performed with a number of layers,N  60 , and a spatial discretization, x  0.1m . The
time step was chosen equal to t  0.002s , fulfilling theCourant-Friedrichs-Lewy (CFL) condition for all theinvestigated flows. The following rheologicalparameters
wereemployed:s  tan  27.0 ,in the upper part of the chute and then flowed down
through an adjustable gate. Different flow rates and flow
regimes were investigated by adjusting the gate opening
between 6 cm and 14 cm. Despite the fact that the flume
inclination was fixed (32.5°), the inclination of the freesurface of the granular flow, as well as the one of theerodible bed, was found to self-adjust depending on themass flow rate. For each experiment, a transient of a few
seconds was followed by a steady state regime, in which
the fully developed surface flow was superimposed on a
static granular wedge. In an intermediate region of the
flume far enough from the inlet and the outlet, almost
uniform flow conditions were observed. The
measurements were taken at a cross section, 1 m
downstream from the gate, by means of two high-speed
cameras (AOS Technologies Corp.) capable ofrecording 1000 frames per second. The sidewall velocity
profiles were obtained using the particle image
velocimetry (PIV) open-source package PIVlab [12],and a 4-pass PIV analysis was carried out [13]. The
velocity field at the free surface, captured by the front
camera, allowed the estimation of a homogenization
coefficient,  , accounting for the deviations from apurely 2D flow (e.g., [11, 14]), and, consequently,permitted the estimation of the width-averaged velocity
profiles, which are directly compared with the computed
profiles by the model.
The volume fraction profiles near the sidewall weremeasured using the stochastic-optical method (SOM)[15]. This method, based on a binarization of thephotographs taken by the side camera, uses a stochastictransfer function, previously obtained by Monte Carlosimulations, and requires controlled illumination
conditions. In this application, a no-flicker LED lamp
(PhotoSonics MultiLED-LT) was placed at a distance of
32cm and with an incidence angle of 25° with respect tothe normal to the sidewall. To ensure the reliability ofmeasurements, the velocity and volume fraction profiles
were estimated by time-averaging the measurements
over a time interval of 5 seconds, and also by averaging
over 4 repetitions of the same experiment.
The simulations by the multilayer model wereconducted on a time interval of 25s, such that fully
developed flows occurred for all runs. The simulations
started from an initially static deposit with an arbitrary
thickness, h, which was made to flow on a bed
inclination,  , equal to the inclination of the erodible
bed in the experiment to be described. Namely, theinitial conditions areh   i 1 hi  x, 0   0.099 m,

2  tan  32.8 , I 0  0.55 , side  tan 10.4 , while
the dilatancy law, (4), was used with the following
parameters:  min  0.2 ,  max  0.6 and a '  0.5 .
With reference to the intermediate run, characterized byan inflow gate opening of 8cm and representative of theentire data set, Fig. 3 illustrates the velocity and volume
fraction profiles computed by the model and theexperimental ones. As it is clear from Fig. 3(a), theagreement between the numerical and experimental
velocity profiles is excellent. Minor discrepancies areobserved in the lower creep region ( h  z  0.03m ),where the model slightly underestimate theexperimental velocities.

Fig. 3. Comparisons between the steady uniform profiles
predicted by the multilayer model (at t=25s) and theexperimental ones (experiment with inflow gate opening
8cm): (a) x-component velocity, u x (b) volume fraction,  .

This is probably due to the chosen rheologicalassumptions (µ(I)-rheology and Coulomb-friction forthe sidewall resistances), which predict zero velocity
from a fixed depth, while the actual behaviour of thelower creep flow consists of an exponential tail [16].
Focusing on the volume fraction profiles (Fig. 3(b)), it
can be noted that the model with the dilatancy law, (4),qualitatively reproduces the general trend of theexperimental  profile. Yet, it is unable to correctly

Nu i  x, 0   0,  i  x, 0   0.6.

(6)

A spatial domain, x   0.5, 0.5 m , is considered witha no-slip bed kinematic condition and non-reflecting
boundaries at x  0.5m and x  0.5m . All the3E3S Web of Conferences 415, 02018 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502018

describe the upper, more dilute flow region,corresponding to a thickness of a few grain diameters ( 0.01m ). This disagreement suggests that the chosendilatancy law, (4), is probably no longer valid in thedilute regime.
The reported study on chute flows over an erodible bed
confirms the good capabilities of the multilayer
approach to reproduce the velocity and volume fraction
fields of dry granular avalanches. Especiallyconsidering its relatively low computational costs [5],the multilayer approach represents a valuable tool toregularize and approximate more computationallyexpensive 3D models based on similar rheologicalassumptions. Moreover, the versatility of the multilayer
framework would allow the future implementation of
more sophisticated dilatancy laws without much effort.
This work is part of the project, "StratifiedGRANULAR",
which received funding from the European Union’s Horizon
2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No. 797890.

References
1.

R. M. Iverson, Rev. Geophys., 35(3), 245-296(1997)
2. A. Armanini, H. Capart, L. Fraccarollo, M.
Larcher, J. Fluid. Mech., 532, 269–319 (2005)
3. E. Audusse, Discrete Cont. Dyn. B, 5(2), 189-214(2005)
4. E. D. Fernández-Nieto, J. Garres-Díaz, A.
Mangeney, G. Narbona-Reina, J. Fluid Mech.,
798, 643-681 (2016)
5. L. Sarno, Y. C. Tai, Y. Wang, M. Oberlack, Phys.
Fluids, 33(10), 103319 (2021)
6. P. Jop, Y. Forterre, O. Pouliquen, Nature,441(7094), 727-730 (2006)
7. F. da Cruz, S. Emam, M. Prochnow, J. N. Roux, F.
Chevoir, Phys. Rev. E, 72(2), 021309 (2005)
8. A. W. Roberts, ASME Trans. J. Eng. Ind., 91,373-381 (1969)
9. J. M. N. T. Gray, A. N. Edwards, J. Fluid Mech.,755, 503 (2014)
10. L. Fraccarollo, H. Capart, Y. Zech, Int. J. Numer.
Meth. Fl., 41(9), 951-976 (2003)
11. L. Sarno, L. Carleo, M. N. Papa, P. Villani, Rock
Mech. Rock Eng., 51(1), 203-225 (2018)
12. W. Thielicke, E. J. Stamhuis, J. Open Res Soft.,
2(1), e30 (2014)
13. L. Sarno, A. Carravetta, Y. C. Tai, R. Martino, M.
N. Papa, C. Y. Kuo, Adv. Powder Technol., 29(12),3107-3123 (2018)
14. L. T. Sheng, C. Y. Kuo, Y. C. Tai, S. S. Hsiau,
Exp. Fluids, 51(5), 1329-1342 (2011)
15. L. Sarno, M. N. Papa, P. Villani, Y. C. Tai, Granul.
Matter, 18(4), 1-12 (2016)
16. T. S. Komatsu, S. Inagaki, N. Nakagawa, S.
Nasuno, Phys. Rev. Lett., 86(9), 1757 (2001)

4