3.2.1 Notation

We denote the granular domain ${\it\Omega}_{F}(t)$ and its projection $I_{F}(t)$ on the reference plane, for a positive $t\in [0,T]$ , i.e.

The multilayer approach considers a vertical partition of the domain in $N\in \mathbb{N}^{\ast }$ layers with preset thicknesses $h_{{\it\alpha}}(t,x)$ for ${\it\alpha}=1,\ldots ,N$ , (see figure 1). Note that $\sum _{{\it\alpha}=1}^{N}h_{{\it\alpha}}=h$ . In practice, we set a vertical partition of the domain as follows: we introduce the positive coefficients $l_{{\it\alpha}}$ such that

Note that the thickness of each layer is automatically adapted to the movement of the free surface, since it depends on the total thickness of the mass flow. These layers are separated by $N+1$ interfaces ${\it\Gamma}_{{\it\alpha}+1/2}(t)$ , which are described by the equations $z=z_{{\it\alpha}+1/2}(t,x)$ for ${\it\alpha}=0,1,\ldots ,N$ , $x\in I_{F}(t)$ . We assume that these interfaces are smooth enough. Observe that the fixed bottom and the free surface are respectively $b=z_{1/2}$ and $b+h=z_{N+1/2}$ , corresponding to the interfaces at the bottom ${\it\Gamma}_{1/2}$ and at the free surface ${\it\Gamma}_{N+1/2}$ respectively. Note that $z_{{\it\alpha}+1/2}=b+\sum _{{\it\beta}=1}^{{\it\alpha}}h_{{\it\beta}}$ and $h_{{\it\alpha}}=z_{{\it\alpha}+1/2}-z_{{\it\alpha}-1/2}$ , for ${\it\alpha}=1,\ldots ,N$ . The subdomain between ${\it\Gamma}_{{\it\alpha}-1/2}$ and ${\it\Gamma}_{{\it\alpha}+1/2}$ is denoted by ${\it\Omega}_{{\it\alpha}}(t)$ , for a positive $t\in [0,T]$ ,

We need to introduce a specific notation about the approximations of the variables on the interfaces (see figure 1). For a function $f$ and for ${\it\alpha}=0,1,\ldots ,N$ , we set

Note that if the function $f$ is continuous,

In addition, for a given time $t$ , we denote

the space–time unit normal vector and the space unit normal vector to the interface ${\it\Gamma}_{{\it\alpha}+1/2}(t)$ outward to the layer ${\it\Omega}_{{\it\alpha}+1}(t)$ for ${\it\alpha}=0,\ldots ,N$ .

3.2.2 Weak solution with discontinuities

We are looking for a particular weak solution $(\boldsymbol{u},p,{\it\rho})$ of (3.11) (see Fernández-Nieto et al. Reference Fernández-Nieto, Koné and Chacón Rebollo2014). We remind the reader that this solution must meet the following conditions:

1. (i) $(\boldsymbol{u},p,{\it\rho})$ is a standard weak solution of (3.11) in each layer ${\it\Omega}_{{\it\alpha}}(t)$ ,

2. (ii) $(\boldsymbol{u},p,{\it\rho})$ satisfies the normal flux jump condition for mass and momentum at the interfaces ${\it\Gamma}_{{\it\alpha}+1/2}(t)$ , namely:

   (3.18) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}[({\it\rho};{\it\rho}\boldsymbol{u})]_{{\it\alpha}+1/2}~\boldsymbol{n}_{T,{\it\alpha}+1/2}=0,\\ \displaystyle \left[\left({\it\rho}\boldsymbol{u};{\it\rho}\boldsymbol{u}\otimes \boldsymbol{u}-\frac{1}{{\it\varepsilon}}{\bf\sigma}\right)\right]_{{\it\alpha}+1/2}~\boldsymbol{n}_{T,{\it\alpha}+1/2}=0,\end{array}\right\}\end{eqnarray}$$

where $[(a;b)]_{{\it\alpha}+1/2}$ denotes the jump of $(a;b)$ across the interface ${\it\Gamma}_{{\it\alpha}+1/2}(t)$ .

There are two main differences with the derivation in Fernández-Nieto et al. (Reference Fernández-Nieto, Koné and Chacón Rebollo2014): first, the ${\it\mu}(I)$ -rheology produces a non-constant viscosity coefficient, which implies an additional difficulty in order to develop the momentum balances at the interfaces. Second, the system to be solved is an asymptotic approximation of the Navier–Stokes equations, which helps to resolve the previous difficulty since we look for a first-order approximation in ${\it\varepsilon}$ .

A particular family of velocity functions is considered, by assuming that the thickness of each layer is small enough to make the horizontal velocities independent of the vertical variable $z$ (usual shallow domain hypothesis). As a consequence, thanks to the incompressibility condition in each layer, we obtain that vertical velocities are linear in $z$ and maybe discontinuous. We denote the velocity in each layer as

where $u_{{\it\alpha}}$ and $w_{{\it\alpha}}$ are the horizontal and vertical velocities, respectively, on layer ${\it\alpha}$ . Then,

for some smooth function $d_{{\it\alpha}}(t,x)$ .

Pressure

From the asymptotic analysis, we directly get a hydrostatic pressure that is obtained from the vertical momentum equation in (3.4):

By the continuity of the dynamic pressure, we can deduce at first order that

Mass conservation across the interfaces and normal velocity

The mass conservation jump condition gives us the definition of the normal mass flux at the interface ${\it\Gamma}_{{\it\alpha}+1/2}(t)$ , denoted by $G_{{\it\alpha}+1/2}:=G_{{\it\alpha}+1/2}^{+}=G_{{\it\alpha}+1/2}^{-}$ for:

By integrating the incompressibility equation and using this mass conservation condition, we get the definition of the vertical velocity $w_{{\it\alpha}}$ for $\boldsymbol{u}_{{\it\alpha}}$ a solution of (3.11) (see Fernández-Nieto et al. Reference Fernández-Nieto, Koné and Chacón Rebollo2014 for details). Thus, if we consider that there is no mass transference with the bottom (i.e. $G_{1/2}=0$ ), we obtain

and for ${\it\alpha}=1,\ldots ,N$ and $z\in (z_{{\it\alpha}-1/2},z_{{\it\alpha}+1/2})$ ,

where

Momentum conservation across the interfaces

The momentum jump conditions become:

where, for ${\it\alpha}=1,\ldots ,N-1$ , the total stress tensor is

with $p_{{\it\alpha}+1/2}=p_{{\it\alpha}+1/2}^{+}=p_{{\it\alpha}+1/2}^{-}$ and ${\bf\tau}_{{\it\varepsilon},{\it\alpha}+1/2}^{\pm }$ approximations of $p_{{\it\alpha}}$ and ${\bf\tau}_{{\it\varepsilon},{\it\alpha}}$ at ${\it\Gamma}_{{\it\alpha}+1/2}$ . Following Fernández-Nieto et al. (Reference Fernández-Nieto, Koné and Chacón Rebollo2014), we finally get that the momentum jump condition gives

where $\widetilde{{\bf\tau}}_{{\it\varepsilon},{\it\alpha}+1/2}$ is an approximation of $({\it\eta}\unicode[STIX]{x1D63F}_{{\it\varepsilon}}(\boldsymbol{u}_{{\it\alpha}}))_{|{\it\Gamma}_{{\it\alpha}+1/2}}$ , defined by

In this equation,

and $({\mathscr{U}_{\mathscr{Z}}^{H}}_{,{\it\alpha}+1/2},{\mathscr{U}_{\mathscr{Z}}^{V}}_{,{\it\alpha}+1/2})$ is defined to approximate the derivatives in $z$ . We remind the reader that we use a mixed formulation because of the possible vertical discontinuous profile. Thus, the auxiliary unknown $\mathscr{U}_{\mathscr{Z}}$ satisfies

And to approximate $\mathscr{U}_{\mathscr{Z}}$ , we approximate $\boldsymbol{u}$ by $\widetilde{\boldsymbol{u}}$ , a $\mathbb{P}_{1}(z)$ interpolation such that $\widetilde{\boldsymbol{u}}_{|z=1/2(z_{{\it\alpha}-1/2}+z_{{\it\alpha}+1/2})}=\boldsymbol{u}_{{\it\alpha}}$ . Then ${\mathscr{U}_{\mathscr{Z}}}_{{\it\alpha}+1/2}=({\mathscr{U}_{\mathscr{Z}}^{H}}_{,{\it\alpha}+1/2},{\mathscr{U}_{\mathscr{Z}}^{V}}_{,{\it\alpha}+1/2})$ is an approximation of $\mathscr{U}_{\mathscr{Z}}(\widetilde{\boldsymbol{u}})$ at ${\it\Gamma}_{{\it\alpha}+1/2}$ .

Let us remark that the previous expression of the tensor $\widetilde{{\bf\tau}}_{{\it\varepsilon},{\it\alpha}+1/2}$ in (3.30) has the same structure as the original case in Fernández-Nieto et al. (Reference Fernández-Nieto, Koné and Chacón Rebollo2014), except for the viscosity that now is not constant because it is defined by the ${\it\mu}(I)$ -rheology in (2.11). Then, we must give an approximation of the viscosity at the interface up to first order in ${\it\varepsilon}$ , which we denoted by ${\it\eta}_{{\it\alpha}+1/2}$ . We have considered the following first-order approximation of $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ at $z=z_{{\it\alpha}+1/2}$ ,

Then, it reads

with

with $h_{{\it\alpha}+1/2}$ being the distance between the midpoints of layers ${\it\alpha}$ and ${\it\alpha}+1$ . Moreover,

and

Note that ${\it\eta}_{N+1/2}=0$ , because we suppose that the atmospheric pressure is zero.

4.1.1 Steady uniform Bagnold flow: analytical solution

Let us first compare the model with the analytical solution for a uniform flow over an inclined plane of slope ${\it\theta}$ and thickness $H>0$ , i.e. a Bagnold flow (see Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001, GDR MiDi 2004 or Lagrée et al. Reference Lagrée, Staron and Popinet2011). This solution is obtained by imposing zero pressure and zero shear stress at the free surface and a no-slip condition at the bottom (Coulomb friction can be easily changed by no-slip condition in the ${\it\mu}(I)$ -rheology multilayer model, see §§ 2.3 and 3.3).

Figure 2.

Sketch of the analytical solution.

By denoting $u$ and $w$ the downslope and normal velocities, $p$ the pressure and ${\it\tau}$ the shear stress and by taking the rheological parameters defined in § 2.2, the steady uniform Bagnold flow is described by

For the numerical simulation, as in the analytical solution, we consider a uniform flow with constant thickness $H=1$  m and velocity $u=w=0~\text{m}~\text{s}^{-1}$ at the initial time $t=0$  s. The boundary conditions at the free surface and at the bottom have been set as in (4.3). At the right and left boundaries, we use open boundary conditions.

Figure 3.

Comparison between the analytical solution (dashed and solid lines) and the simulations obtained using the MSM with the ${\it\mu}(I)$ -rheology (symbols). (a) Analytical and simulated downslope horizontal velocity $u$ , pressure $p$ and strain rate $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ . (b) Analytical and simulated shear stress and friction coefficient ${\it\mu}(I)$ . (c) Comparison between the simulated (symbols) and the exact (dashed line) horizontal velocity at the free surface as a function of the slope angle.

Figure 4.

(a) Computing time as a function of the number of layers in the MSM and (b) relative error between the computed and exact velocity for simulations over a slope of ${\it\theta}=24.64^{\circ }$ .

We choose the rheological parameters $I_{0}=0.279$ and ${\it\mu}_{s}=0.38\approx \tan (20.8^{\circ })$ , ${\it\mu}_{2}=0.62\approx \tan (31.8^{\circ })$ and the particle diameter $d_{s}=4$  cm with solid volume fraction ${\it\varphi}_{s}=0.62$ . The slope angle is taken as ${\it\theta}=0.43~\text{rad}\approx 24.64^{\circ }$ . Figure 3 shows the good agreement between the simulated and exact solutions for the profiles of the velocity, pressure, shear stress, ${\it\mu}(I)$ and $\Vert \unicode[STIX]{x1D63F}(u)\Vert$ . It also shows the downslope velocity at the free surface as a function of the slope angle. Note that for slopes smaller than $\arctan ({\it\mu}_{s})$ , the surface velocity is close to zero (it is not equal to zero, because of the regularization method), namely $u(z=H)\sim 10^{-5}$ , so that the mass is almost at rest.

These results are computed using 50 layers in the MSM. Note that we use a slope in the well-posed region described in Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) for the full ${\it\mu}(I)$ -rheology, which is smaller than the well-posed region for the depth-averaged ${\it\mu}(I)$ -rheology (well-posed for ${\it\delta}_{s}=20.8^{\circ }\leqslant {\it\theta}<{\it\delta}_{2}=31.8^{\circ }$ in the case of depth-averaged ${\it\mu}(I)$ -rheology). If we use a slope close enough to ${\it\delta}_{2}$ the system becomes unstable because of the ill-posedness of the full ${\it\mu}(I)$ -rheology in this region.

Figure 4 shows the computing time required to simulate 50 seconds (on a laptop with Intel $^{\unicode[STIX]{x00AE} }$ Core™ i7-4500U and 8 GB of RAM) and the relative error between the computed velocity and the exact solution using a different number of layers and 30 nodes in the $x$ -direction. In table 1 we show that second-order accuracy is reached in norms $L^{1}$ and $L^{2}$ , while in norm $L^{\infty }$ the order is over 1.5. This behaviour in norm $L^{\infty }$ is due to the boundary conditions at the free surface, where the analytical solution verifies that $\partial _{z}u=0$ .

Table 1.

Order of the error for the velocity of the Bagnold flow.

4.1.2 2D steady uniform flow in a channel: sidewalls effect on the velocity profile

The relevance of the lateral wall friction when a granular material flows in a channel has been proved by Jop et al. (Reference Jop, Forterre and Pouliquen2005) (see also Baker, Barker & Gray Reference Baker, Barker and Gray2016). Under the hypothesis of a steady uniform flow and neglecting the variations in the transverse direction, they propose to model this effect with a modified ${\it\mu}(I)$ -rheology by adding an extra term, which defines an effective friction term

where $W$ is the channel width, ${\it\mu}_{w}$ is the constant coefficient of friction with the side walls and $H$ the thickness of the flow. Note that this extra term (last term) increases when we get closer to the bottom from the free surface. In addition, for a channel slope ${\it\theta}$ with respect to the horizontal axis, the analytical expression for the velocity profile reads (see appendix B in Jop et al. Reference Jop, Forterre and Pouliquen2005):

where

Note that we can integrate the last equation and the discharge is obtained as a function of the slope ${\it\theta}$ . In these experiments, Jop et al. (Reference Jop, Forterre and Pouliquen2005) fixed the discharge and computed the velocity profile for different widths. We choose the same material and rheological properties and perform these experiments. The grain diameter is $d_{s}=0.53$  mm and the volume fraction is ${\it\varphi}_{s}=0.6$ . The rheological parameters are ${\it\mu}_{s}=\tan (20.9^{\circ })$ , ${\it\mu}_{2}=\tan (32.76^{\circ })$ and $I_{0}=0.279$ . The friction coefficient with the sidewalls is ${\it\mu}_{w}=\tan (10.4^{\circ })$ . We set the thickness of the flow $H=45d_{s}$ , the dimensionless discharge $Q^{\ast }=31.5$ and different widths of the channel are considered ( $W=19d_{s}$ , $57d_{s}$ , $283d_{s}$ ).

Figure 5 shows the exact agreement between the analytical solution and the simulations for the velocity profile. In this case we consider 50 layers to obtain the solution. In table 2 we present the errors and accuracy order of the method for $W=283d_{s}$ . We obtain similar results to those presented in the previous subsection.

Figure 5.

Comparison between the analytical solution (solid lines) and the simulations obtained using the MSM with the ${\it\mu}(I)$ -rheology (symbols) for the velocity profile for $Q^{\ast }=31.5$ , and different widths: $W=19d_{s}$ , $W=57d_{s}$ and $W=283d_{s}$ .

Table 2.

Order of the error for the velocity of the Bagnold flow with lateral wall effect. Case $W=283d_{s}$ .

4.1.3 Laboratory experiments: transient velocity profiles

We compare our model with the experiments and simulations presented in Jop et al. (Reference Jop, Forterre and Pouliquen2007) where granular material is flowing within a narrow channel of width $W=19d_{s}\approx 1$  cm. The grain diameter is $d_{s}=0.53$  mm, the volume fraction is ${\it\varphi}_{s}=0.6$ and the rheological parameters are ${\it\mu}_{s}=\tan (20.9^{\circ })$ , ${\it\mu}_{2}=\tan (32.76^{\circ })$ and $I_{0}=0.279$ . The friction with the wall is modelled as in the precedent test, following (4.4). In this case, the authors set the friction coefficient with the sidewalls as ${\it\mu}_{w}=\tan (13.1^{\circ })$ .

For the simulations, we impose zero velocity at the initial time and the material flows because of the gravitational force. We use 50 layers in the MSM. Figure 6 shows the comparison between the results in Jop et al. (Reference Jop, Forterre and Pouliquen2007) and the simulation using the MSM ${\it\mu}(I)$ -rheology model for the transient velocity profiles and its final state in two configurations. As is concluded in Jop et al. (Reference Jop, Forterre and Pouliquen2007), the agreement for the transient velocity profiles is better for high slopes and less accurate for low inclinations. In the case of high inclination angles our model results are similar to their simulations (see figure 6 b), and are slightly better in the other cases (see figure 6 a).

Figure 6.

Comparison between the laboratory experiments (symbols) and simulations (dash-dotted lines) in Jop et al. (Reference Jop, Forterre and Pouliquen2007), and simulations using the MSM with the ${\it\mu}(I)$ -rheology (solid lines) for the transient velocity profiles for two different slopes ${\it\theta}$ , and times $t^{\ast }=t/\sqrt{d_{s}/g}$ : (a) ${\it\theta}=26.1^{\circ }$ , $t^{\ast }=1.2,15.1,166.1$ ; (b) ${\it\theta}=32.15^{\circ }$ , $t^{\ast }=2.3,7.5,24.2,77.8,175$ .

4.2.1 Improvement of the approximation of the ${\it\mu}(I)$ coefficient

The advantage of the multilayer models is that we obtain a variable profile of the downslope velocity, in contrast with the prescribed profile of the monolayer model. It makes it possible to obtain a better approximation of $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ (see (3.33)). As a consequence, this improves the approximation of the inertial number $I$ (see (2.7) and (3.37)), which is a key number in the variable friction coefficient ${\it\mu}(I)$ .

As the main advantage of the multilayer model is the improvement of the approximation of $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ , we present two approximations that can be made with the multilayer model. First, let us recall that a first-order approximation corresponds to the definition (3.33). This approximation considers only the leading-order term, i.e. $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert _{{\it\alpha}+1/2}\approx \Vert \partial _{z}u_{{\it\alpha}+1/2}\Vert =\Vert {\mathscr{U}_{\mathscr{Z}}^{H}}_{,{\it\alpha}+1/2}\Vert$ . Note that in dimensionless form, we have

We can improve the approximation of $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ at the interfaces $z=z_{{\it\alpha}+1/2}$ by considering the approximation taking into account second-order terms in the previous equation. For the numerical tests, we consider the following approximation $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ at the interfaces,

Note that this definition corresponds to an approximation of

at $z=z_{{\it\alpha}+1/2}$ . Nevertheless, in (4.7), the term $2\partial _{z}u\,\partial _{x}w$ is not taken into account although it is of the same order as $4(\partial _{x}u)^{2}$ . This is because when an approximation of this term is added, we obtain results that are very similar to those obtained when considering (4.8). Furthermore, adding this term implies an additional computational cost since pre-calculated vertical velocities are required. Note that (4.8) is a second-order correction while we have developed a first-order model that neglects other second-order terms. This correction, however, highlights the importance of second-order terms in granular collapses over erodible beds, even if it is a partial correction.

The effect of this second-order approximation on the results is discussed below in §§ 4.2.3 and 4.2.4. In particular, we observe an improvement of the results compared to the original first-order approximation of $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ .

Figure 7.

Sketch of the initial and final state of the granular collapse. A granular column with a thickness $h_{0}=14$  cm and a length $r_{0}=20$  cm is released on an inclined plane of slope ${\it\theta}$ . The plane is covered by an erodible bed of thickness $h_{i}$ made of the same material. When the flow stops, the maximum final thickness is $h_{f}$ and its final extent $r_{f}$ .

4.2.2 Experimental and test data

In the laboratory experiments performed in Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010), subspherical glass beads of diameter $d_{s}=0.7$  mm were used. The particle density ${\it\rho}_{s}=2500~\text{kg}~\text{m}^{-3}$ and volume fraction ${\it\varphi}_{s}=0.62$ were estimated, leading to an apparent flow density ${\it\rho}={\it\varphi}_{s}{\it\rho}_{s}=1550~\text{kg}~\text{m}^{-3}$ .

The variable $r_{f}$ denotes the runout distance, i.e. the length of the deposit measured from the position of the front of the released material at the initial time located at $x=0$  m, $t_{f}$ denotes the flow time from $t=0$  s to the time when the material stops and $h_{f}$ denotes the maximum final thickness of the deposit (see figure 7).

In order to use the ${\it\mu}(I)$ -rheology, the rheological parameters ( ${\it\mu}_{s}$ , ${\it\mu}_{2}$ and $I_{0}$ ) must be defined. We consider the data proposed in Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015). The minimum and maximum friction angles are ${\it\mu}_{s}=\tan (20.9^{\circ })$ and ${\it\mu}_{2}=\tan (32.76^{\circ })$ , according to the measurements made in the experiments presented in Pouliquen & Forterre (Reference Pouliquen and Forterre2002) and Jop et al. (Reference Jop, Forterre and Pouliquen2005). These parameters can be obtained by fitting the curve $h_{stop}({\it\theta})$ , where $h_{stop}$ is the thickness of the deposit lying on the slope when the supply is stopped after steady uniform flow (see Pouliquen (Reference Pouliquen1999b ) for more details). Nevertheless, in Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) the value of ${\it\mu}_{s}$ and ${\it\mu}_{2}$ are incremented, in order to consider the effect of lateral wall friction. Let us remember that lateral wall friction is modelled in Jop et al. (Reference Jop, Forterre and Pouliquen2005) as an additional friction term ${\it\mu}_{w}(h-z)/W$ , where ${\it\mu}_{w}=\tan (10.4^{\circ })$ is the coefficient of friction in the side walls. Moreover, the thickness of the flowing layer (see Mangeney et al. Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010) is approximately 0.05 m and the width of the channel $W=10$  cm. Therefore, the additional friction term is approximately 0.1. As a result, in Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) the authors propose to consider ${\it\mu}_{s}=\tan (25.5^{\circ })\approx \tan (20.9^{\circ })+0.1$ and ${\it\mu}_{2}=\tan (36.5^{\circ })\approx \tan (32.76^{\circ })+0.1$ . Moreover, we set $I_{0}=0.279$ (see Jop et al. Reference Jop, Forterre and Pouliquen2006).

That might be a coarse way to introduce the wall friction effect, owing that the multilayer model is able to approximate the term ${\it\mu}_{w}(h-z)/W$ in each layer. However, introducing the friction term ${\it\mu}_{w}(h-z)/W$ does not give good results in this granular collapse test. Actually it seems that there is not enough friction at the initial times to well capture the S-shaped profile. Note that the numerical simulations with the proposed hydrostatic multilayer approach cannot be accurate at short times because the dominant effect is the non-hydrostatic pressure. At short times the flowing layer can be overestimated, and, consequently, the effect of the lateral wall friction is not properly taken into account. These comparisons only make sense at the latter stage of the flow. We plot the results for intermediate times in order to illustrate this discussion (see figure 15).

This experiment has been simulated for different slopes ${\it\theta}$ and thicknesses $h_{i}$ of the erodible bed: ${\it\theta}=16^{\circ }$ and $h_{i}=1.4,2.5,5$  mm, ${\it\theta}=19^{\circ }$ and $h_{i}=1.5,2.7,5.3$  mm, ${\it\theta}=22^{\circ }$ and $h_{i}=1.82,3.38,4.6$  mm, and ${\it\theta}=23.7^{\circ }$ and $h_{i}=1.5,2.5,5$  mm. Note that the model does not take into account the effect of removing the gate during the initial instants even though it has a non-negligible impact on the flow dynamics as shown in Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015). For instance, when the gate is taken into account, even with no friction along it, the flow is substantially slowed down; however, the deposit is almost unchanged. All the simulations are performed using 20 layers.

We compare hereafter (i) the constant and variable friction rheologies and (ii) the monolayer and multilayer approaches. In table 3, we summarize the notation and symbols used for the different models. MSM is the notation for multilayer models. Note that the monolayer model with a constant friction coefficient (denoted ${\it\mu}_{s}$ -monolayer model) corresponds to the Savage–Hutter classical model when $K_{act/pas}=1$ . This model is widely used in the literature, e.g. Gray, Tai & Noelle (Reference Gray, Tai and Noelle2003), Mangeney-Castelnau et al. (Reference Mangeney-Castelnau, Bouchut, Vilotte, Lajeunesse, Aubertin and Pirulli2005). Note also that the monolayer model with a variable friction coefficient (denoted ${\it\mu}(I)$ -monolayer model) corresponds to the one proposed by Pouliquen (Reference Pouliquen1999b ), which has been recently used by Mangeney-Castelnau et al. (Reference Mangeney-Castelnau, Vilotte, Bristeau, Perthame, Bouchut, Simeoni and Yerneni2003) and Mangeney et al. (Reference Mangeney, Bouchut, Thomas, Vilotte and Bristeau2007). This model also coincides with the one proposed in Gray & Edwards (Reference Gray and Edwards2014) by dropping second-order terms, that is, viscous terms.

Table 3.

Summary of notation of the different models and colours/symbols.

4.2.3 Deposit profiles

Let us compare the deposits simulated with the ${\it\mu}(I)$ -rheology and with a constant friction coefficient ${\it\mu}_{s}$ for different slopes ${\it\theta}$ and erodible bed thicknesses $h_{i}$ . Figure 8 shows that the deposit calculated with the variable friction coefficient ${\it\mu}(I)$ is closer to the experimental deposit than the one calculated with a constant friction coefficient ${\it\mu}_{s}$ . The runout distance with the constant coefficient ${\it\mu}_{s}$ is always too long except for ${\it\theta}=19^{\circ }$ and $h_{i}=5.3$  mm (see figure 8 d). To properly reproduce the runout distance with a constant friction coefficient, we need to increase its value. For example, with a slope ${\it\theta}=16^{\circ }$ and an erodible bed thickness $h_{i}=2.5$  mm (figure 8 a), we need to use the value ${\it\mu}_{s}=\tan (27.3^{\circ })$ to produce the runout observed in the laboratory experiments. That means an increment of 0.039 in the ${\it\mu}_{s}$ value. These results are consistent with the simulations of Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015), showing that the runout is strongly overestimated when the viscosity tends to zero (i.e. when ${\it\mu}$ tends to ${\it\mu}_{s}$ ).

Figure 8.

Deposit obtained in the experiments (solid-circle blue line), with the ${\it\mu}_{s}$ -MSM (dotted-circle red line) and with the ${\it\mu}(I)$ -MSM (solid-cross green line), for different slopes ${\it\theta}$ and erodible bed thicknesses $h_{i}$ .

Figure 9.

Deposit obtained in the experiments (solid-circle blue line), with the ${\it\mu}_{s}$ -MSM (dotted-circle red line), the ${\it\mu}(I)$ -MSM (solid-cross green line), the ${\it\mu}_{s}$ -monolayer model (dashed red line) and with the ${\it\mu}(I)$ -monolayer model (solid green line) for a slope ${\it\theta}=22^{\circ }$ and an erodible bed thickness $h_{i}=1.82$  mm. Circled sections indicate specific locations for data shown in figure 10.

Figure 9 shows, for a slope ${\it\theta}=22^{\circ }$ and $h_{i}=1.82$  mm, the final deposit obtained using the constant or variable friction coefficients for multilayer and monolayer models. The difference between the multilayer and monolayer models is stronger when using the ${\it\mu}(I)$ -rheology. For instance, the multilayer approach changes the full deposit profiles for the ${\it\mu}(I)$ -rheology, while it only changes the front position for ${\it\mu}_{s}$ . The multilayer approach makes it possible to obtain a deposit shape which is very close to the experiments with the ${\it\mu}(I)$ -rheology. More generally, the shape of the deposit is closer to the observations with ${\it\mu}(I)$ -MSM than with ${\it\mu}_{s}$ -MSM.

4.2.4 Effect of the erodible bed

Figure 10 shows two zooms, one near the maximum thickness of the deposit (zone of circle (I) in figure 9), and one near the front (zone of circle (II) in figure 9), for ${\it\theta}=22^{\circ }$ and different values of $h_{i}$ . With the ${\it\mu}(I)$ -MSM, the runout distance $r_{f}$ increases as the thickness of the erodible bed $h_{i}$ increases (figure 10 a (II)) as observed in laboratory experiments. On the contrary, with the ${\it\mu}_{s}$ -MSM (figure 10 b (II)), the runout distance $r_{f}$ decreases with increasing $h_{i}$ . Note that in both cases the maximum final thickness $h_{f}$ decreases with increasing $h_{i}$ as it occurs in the experiments (figure 10 a (I),b (I)).

Figure 10.

Influence of the thickness of the erodible bed on the runout distance $r_{f}$ and on the maximum final thickness $h_{f}$ (inset graphs) with the ${\it\mu}(I)$ -MSM (a) and with the ${\it\mu}_{s}$ -MSM (b), for a slope ${\it\theta}=22^{\circ }$ , in zones marked with circles in figure 9.

Figure 11 shows that the decrease in runout distance with increasing $h_{i}$ for constant friction ${\it\mu}_{s}$ is observed for all slopes, e.g.  ${\it\theta}=0^{\circ },10^{\circ },16^{\circ },19^{\circ },22^{\circ }$ and $23.7^{\circ }$ . For the constant friction coefficient case, the ${\it\mu}_{s}$ -MSM and ${\it\mu}_{s}$ -monolayer models follow the same trend. Note that this non-physical decrease in runout distance with increasing $h_{i}$ has been demonstrated analytically in Faccanoni & Mangeney (Reference Faccanoni and Mangeney2013) for the monolayer model. Moreover, laboratory experiments show that when the thickness of the erodible bed increases, for slopes ${\it\theta}\geqslant {\it\theta}_{c}$ , where ${\it\theta}_{c}\in [12^{\circ },16^{\circ }]$ is a critical slope, the runout distance $r_{f}$ and the stopping time $t_{f}$ both increase while the maximum final thickness $h_{f}$ decreases. Note that there is no pattern concerning the runout when the thickness $h_{i}$ is increased for slopes ${\it\theta}<{\it\theta}_{c}$ ( ${\it\theta}=0^{\circ },10^{\circ }$ ) in the laboratory experiments.

Figure 11.

Influence of the thickness $h_{i}$ of the erodible bed on the final runout $r_{f}$ for slopes ${\it\theta}=0^{\circ },10^{\circ },16^{\circ },19^{\circ },22^{\circ }$ and $23.7^{\circ }$ observed in the experiments of Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010) (solid-circle blue line) and obtained with different simulations using the ${\it\mu}_{s}$ -MSM, with 20 layers (dotted-circle red line) and the ${\it\mu}_{s}$ -monolayer model, i.e. the Savage–Hutter model (dashed red line). There is no laboratory data for ${\it\theta}=23.7^{\circ }$ . Normalization is made using $h_{0}=14$  cm.

Figure 12 shows that the increase of runout distance observed in the experiments for increasing $h_{i}$ is qualitatively well reproduced with the ${\it\mu}(I)$ -MSM. With the ${\it\mu}(I)$ -MSM, the runout increase with $h_{i}$ is actually larger for higher slopes, as observed experimentally: at ${\it\theta}=16^{\circ }$ , the runout distance is almost unaffected by the thickness of the erodible bed while it increases by 26.9 % at ${\it\theta}=22^{\circ }$ when the thickness of the erodible bed increase from 1.82 to 4.6 mm. Note that in the ${\it\mu}(I)$ -MSM, the increase of the runout distance appears on slopes ${\it\theta}>16^{\circ }$ , higher than ${\it\theta}_{c}$ in the experiments. Actually, it appears starting with the slope ${\it\theta}=18^{\circ }$ . When using the ${\it\mu}(I)$ -monolayer model, the runout distance is higher than for the ${\it\mu}(I)$ -MSM whatever the slope and thickness of the erodible bed. Based on the values of the runout distance in these cases, it is hard to discriminate which of the monolayer or multilayer models is closer to the experiments. However, in the ${\it\mu}(I)$ -monolayer model, the runout distance at ${\it\theta}=16^{\circ }$ and $19^{\circ }$ decreases when $h_{i}$ increases, contrary to the experimental data. For ${\it\theta}=22^{\circ }$ and ${\it\theta}=23.7^{\circ }$ , the monolayer and multilayer ${\it\mu}(I)$ models reproduce qualitatively the increase in runout with $h_{i}$ . Note that for ${\it\theta}=0^{\circ },10^{\circ }$ ( ${\it\theta}<{\it\theta}_{c}$ ), the ${\it\mu}(I)$ models predict a very slight decrease in the runout distance.

The ${\it\mu}(I)$ -MSM is the model corresponding to the multilayer approach with the ${\it\mu}(I)$ -rheology when $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ is approximated by the main term in (4.7), and ${\it\mu}(I)$ -C-MSM when $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ is approximated by the correction (4.8). Figure 12 shows that the correction of $\Vert \unicode[STIX]{x1D63F}(\boldsymbol{u})\Vert$ corresponding to ${\it\mu}(I)$ -C-MSM improves the simulation of both the runout extent and the influence of the erodible bed. They both increase the runout when $h_{i}$ increases, although the effect of erosion is still much smaller than in the experiments.

For example, for ${\it\theta}=19^{\circ }$ , when $h_{i}$ varies from 1.5 to 5.3 mm, the experimental runout increases by 15.1 %. On the contrary, for the monolayer model the runout decreases by 1.1 %. For the ${\it\mu}(I)$ -MSM the runout increases by 2.5 %, and for the ${\it\mu}(I)$ -C-MSM it increases by 6.8 %. For ${\it\theta}=22^{\circ }$ , when $h_{i}$ varies from 1.82 to 4.6 mm, the experimental runout increases by 26.9 %. It increases by 2.2 %, 4.4 % and 8.6 % for the monolayer model, the ${\it\mu}(I)$ -MSM and the ${\it\mu}(I)$ -C-MSM, respectively.

In figure 13, the final time (time at which the front stops) is plotted as a function of the thickness of the erodible bed for ${\it\theta}=16^{\circ }$ , ${\it\theta}=19^{\circ }$ and ${\it\theta}=22^{\circ }$ . Moreover, for ${\it\theta}=22^{\circ }$ , we also plot the experimental data. Experimental data show that the final time increases when the thickness of the erodible bed increases. We can see that this is true for all the values of ${\it\theta}$ for the multilayer method. However, we observe that it is only true for the highest value, ${\it\theta}=22^{\circ }$ , in the case of the monolayer model, whereas the final time decreases when the erodible bed increases for ${\it\theta}=16^{\circ }$ and ${\it\theta}=19^{\circ }$ .

Figure 12.

Influence of the thickness $h_{i}$ of the erodible bed on the final runout $r_{f}$ for slopes ${\it\theta}=16^{\circ },19^{\circ },22^{\circ }$ and $23.7^{\circ }$ observed in the experiments of Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010) (solid-circle blue line) and obtained with different simulations using the ${\it\mu}(I)$ -MSM, with 20 layers (solid-cross green line), with the ${\it\mu}(I)$ -monolayer model (solid green line) and with the ${\it\mu}(I)$ -C-MSM (dashed black line). There is no laboratory data for ${\it\theta}=23.7^{\circ }$ . Normalization is made using $h_{0}=14$  cm.

Figure 13.

Influence of the thickness $h_{i}$ of the erodible bed on the final time $t_{f}$ for slopes ${\it\theta}=16^{\circ },19^{\circ }$ and $22^{\circ }$ by using the ${\it\mu}(I)$ -MSM with 20 layers (solid-cross green line), the ${\it\mu}(I)$ -monolayer model (solid green line) and the values observed in the experiments of Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010) for slope ${\it\theta}=22^{\circ }$ (solid-circle blue line). Normalization is made using ${\it\tau}_{c}=\sqrt{h_{0}/(g\cos {\it\theta})}$ and $h_{0}=14$  cm.

4.2.5 Flow dynamics and velocity profiles

Figures 14 and 15 show the time evolution of the granular column thickness for a slope ${\it\theta}=22^{\circ }$ and an erodible bed of thickness $h_{i}=1.82$  mm for ${\it\mu}_{s}$ and ${\it\mu}(I)$ , respectively, for both the monolayer and multilayer models. As observed for the deposit, the difference between the thickness profiles simulated with the multilayer and the monolayer model is stronger for ${\it\mu}(I)$ than for ${\it\mu}_{s}$ . The ${\it\mu}(I)$ -MSM makes it possible to increase the maximum thickness of the flow and decrease the thickness of the front. This is an important result as the shape of the front may be an indicator of the flow rheology (Pouliquen Reference Pouliquen1999a ; Jessop et al. Reference Jessop, Kelfoun, Labazuy, Mangeney, Roche, Tillier, Trouillet and Thibault2012).

Figure 14.

Thickness of the granular mass at different times in the experiments (solid-circle blue line), with the ${\it\mu}_{s}$ -MSM with 20 layers (dotted-circle red line), with the ${\it\mu}_{s}$ -monolayer model (Savage–Hutter model, dashed red line) and with the flow/non-flow interface (dashed brown line), for the slope ${\it\theta}=22^{\circ }$ and erodible bed thickness $h_{i}=1.82$  mm.

When a constant coefficient ${\it\mu}_{s}$ is used, very similar profiles are obtained with the ${\it\mu}_{s}$ -MSM and ${\it\mu}_{s}$ -monolayer model (Savage–Hutter model). As a result, the multilayer approach does not significantly improve the results when a constant friction coefficient is used. Note that during the initial instants, the simulated mass spreads faster than in the experiments. This is partly due to the role of initial gate removal that is not taken into account here. However, this effect does not explain the strong difference between the simulation and experiments (see Ionescu et al. Reference Ionescu, Mangeney, Bouchut and Roche2015 for more details). The hydrostatic assumption may also be responsible for this overestimation of the spreading velocity (see e.g. Mangeney-Castelnau et al. Reference Mangeney-Castelnau, Bouchut, Vilotte, Lajeunesse, Aubertin and Pirulli2005). In these figures we also show the flow/no-flow interfaces computed with the multilayer models ${\it\mu}_{s}$ -MSM, ${\it\mu}(I)$ -MSM and ${\it\mu}(I)$ -C-MSM, by considering a threshold velocity of $0.001~\text{m}~\text{s}^{-1}$ to compute these interfaces. We see in figure 14 that the flow/no-flow interface has a step-like shape for ${\it\mu}_{s}$ -MSM while it has a smoother shape for ${\it\mu}(I)$ -MSM (see figure 15) in better qualitative agreement with laboratory experiments and numerical modelling solving the full Navier–Stokes equations (e.g. figures 9 and 18 of Ionescu et al. Reference Ionescu, Mangeney, Bouchut and Roche2015). The flow/no-flow interface in the uppermost part even seems to be improved when using ${\it\mu}(I)$ -C-MSM. For ${\it\mu}(I)$ -MSM and ${\it\mu}(I)$ -C-MSM, the lower layers stop before the upper layers as observed experimentally. Note that in figures 14 and 15 we do not see the flow/no-flow interface at the final time because the material has already stopped. As was shown in § 4.1.1, the wall friction effect is crucial to determine the position of the flow/no-flow interface. Because in this test we do not consider the exact definition of the wall friction term, as in Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015), the flow/no-flow interface should not be very well captured.

Figure 15.

Thickness of the granular mass at different times in the experiments (solid-circle blue line), with the ${\it\mu}(I)$ -MSM using 20 layers (solid-cross green line), with the ${\it\mu}(I)$ -monolayer model (solid green line), with the ${\it\mu}(I)$ -C-MSM (dashed black line), and the flow/non-flow interface with both multilayer models (dashed brown/solid brown lines) for the slope ${\it\theta}=22^{\circ }$ and erodible bed thickness $h_{i}=1.82$  mm.

Figure 16 shows that the second-order correction in ${\it\mu}(I)$ -C-MSM leads to simulated deposits that are generally closer to the experimental observations than those calculated with ${\it\mu}(I)$ -MSM. In particular, the deposits at ${\it\theta}=19^{\circ }$ and ${\it\theta}=22^{\circ }$ with $h_{i}=4.6$  mm are very well reproduced (figure 16 b,c,d,f). However, in some cases, ${\it\mu}(I)$ -MSM gives better results than ${\it\mu}(I)$ -C-MSM, for example for ${\it\theta}=22^{\circ }$ with $h_{i}=1.82$  mm. This is true for the overall dynamics as illustrated in figure 15 that shows the time change of the granular column thickness. We can see that with ${\it\mu}(I)$ -C-MSM, the avalanche is faster and the runout is overestimated and very similar to the runout obtained with the ${\it\mu}(I)$ monolayer model. As other second-order terms than those included in the ${\it\mu}(I)$ -C-MSM model are neglected, it is not easy to draw a firm conclusion on the improvement of results when using second-order terms.

Figure 16.

Simulated deposits at different slopes ${\it\theta}$ and erodible bed thicknesses $h_{i}$ with the ${\it\mu}(I)$ -C-MSM (dashed black line) and the ${\it\mu}(I)$ -MSM (solid-cross green line). The deposits observed in the experiments are represented by solid-circle blue lines.

The multilayer approach makes it possible to obtain a normal profile of the downslope velocity. Figures 17 and 18 show the normal profiles of the downslope velocity obtained at different times until the mass stops, for two different configurations of slopes and erodible beds. In order to obtain a more accurate profile, 40 layers are used in the ${\it\mu}(I)$ -MSM. The different kind of profiles observed in figures 17 and 18 are in good qualitative agreement with typical velocity profiles of granular flows GDR MiDi (2004) (see also Lusso et al. Reference Lusso, Bouchut, Ern and Mangeney2015a and Lusso et al. Reference Lusso, Ern, Bouchut, Mangeney, Farin and Roche2015b ), from Bagnold-like to S-shaped profiles. This shows the ability of the model to recover velocity profiles in a wide range of regimes.

Finally, let us compare the averaged velocity obtained with the ${\it\mu}(I)$ -monolayer model to the average of the velocities over all the layers in the ${\it\mu}(I)$ -MSM. In figure 17, for the green profile (respectively red and magenta profiles), the velocity in the monolayer model is $1.01~\text{m}~\text{s}^{-1}$ (respectively 0.02 and $0.14~\text{m}~\text{s}^{-1}$ ) and it is $0.95~\text{m}~\text{s}^{-1}$ in the multilayer model (respectively 0.03 and $0.05~\text{m}~\text{s}^{-1}$ ). Note that we obtain similar values for the first and second profiles. For the third profile, the averaged velocities strongly differ. Actually, at this position and time, the velocity profile corresponds to the stopping phase for the multilayer model but not for the monolayer model. As a result, the velocity obtained in the multilayer model is lower than that obtained in the monolayer model. Figure 19 shows the normal profile of normal velocity for the same configuration as figure 17. Note that the normal velocities are always negative and that their absolute values are greater in the upper layers.

Figure 17.

Normal profiles of the downslope velocity obtained with the ${\it\mu}(I)$ -MSM (40 layers) for ${\it\theta}=22^{\circ }$ and $h_{i}=1.82$  mm during granular collapse at different positions ( $x=0.095,0.495,0.995$  m). For these positions, we represent the velocity profiles for different times, taken every 0.15 s (blue lines). The first selected profile (green) shows a profile at the beginning of the flow and the second (red) and third (magenta) profiles were measured during the stopping stage. The final deposit is represented by the solid brown line.

Figure 18.

Normal profiles of the downslope velocity obtained with the ${\it\mu}(I)$ -MSM (40 layers) for ${\it\theta}=0^{\circ }$ and $h_{i}=1.5$  mm during granular collapse at different positions ( $x=0.045,0.245$  m) and times taken every 0.05 s.

Figure 19.

Normal profiles of the normal velocity obtained with the ${\it\mu}(I)$ -MSM (40 layers) for ${\it\theta}=22^{\circ }$ and for $h_{i}=1.82$  mm during granular collapse at different positions ( $x=0.095,0.495,0.995$  m) and times taken every 0.2 s.