  In this work we present a multilayer shallow model to approximate the
Navier-Stokes equations with hydrostatic pressure and the $\mu(I)$-rheology.
The main advantages of this approximation are (i) the low cost associated with
the numerical treatment of the free surface of the modelled flows, (ii) exact
conservation of mass and (iii) the ability to compute 3D profiles of the
velocities in the directions along and normal to the slope. The derivation of
the model follows [14] and introduces a dimensional analysis based on the
shallow flow hypothesis. The proposed first order multilayer model fully
satisfies a dissipative energy equation. A comparison with an analytical
solution with a non-constant normal profile of the downslope velocity
demonstrates the accuracy of the numerical model. Finally, by comparing the
numerical results with experimental data, we show that the proposed multilayer
model with the $\mu(I)$-rheology reproduces qualitatively the effect of the
erodible bed on granular flow dynamics and deposits, such as the increase of
runout distance with increasing thickness of the erodible bed. We show that the
use of a constant friction coefficient in the multilayer model leads to the
opposite behaviour. This multilayer model captures the different normal
profiles of the downslope velocity during the different phases of the flow
(acceleration, stopping, etc.) including the presence of static and flowing
zones within the granular column.
