1. Introduction
Understanding the mechanical behaviour of granular media is essential in various pharmaceutical and manufacturing processes as well as in geophysical events like landslides, debris flows and snow avalanches (Mangeney et al. Reference Mangeney, Bouchut, Thomas, Vilotte and Bristeau2007; Persson, Alderborn & Frenning Reference Persson, Alderborn and Frenning2011; Lucas, Mangeney & Ampuero Reference Lucas, Mangeney and Ampuero2014; Steinkogler et al. Reference Steinkogler, Gaume, Löwe, Sovilla and Lehning2015; Delannay et al. Reference Delannay, Valance, Mangeney, Roche and Richard2017). As the number of grains in such flows are usually too high to be modelled in a discrete manner, accurate continuum models are needed to efficiently simulate these large-scale natural phenomena and contribute to improving, e.g. hazard assessment and mitigation in mountainous regions. Significant progress has been made, in particular over the last decades, with improved constitutive and rheological models as well as new numerical schemes. Despite this progress, it remains challenging to achieve general models that can describe the diverse behaviour of granular systems. The diversity stems from the varying apparent friction and solid fraction, as well as the effect of possible inter-granular cohesive forces, or even elastic forces in the slow flow regime (Campbell Reference Campbell2002, Reference Campbell2005; Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006; Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006; Pirulli & Mangeney Reference Pirulli and Mangeney2008; Rao & Nott Reference Rao and Nott2008; Rognon et al. Reference Rognon, Roux, Naaïm and Chevoir2008b; Moretti et al. Reference Moretti, Mangeney, Capdeville, Stutzmann, Huggel, Schneider and Bouchut2012; Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013; Radjai, Roux & Daouadji Reference Radjai, Roux and Daouadji2017; Moretti et al. Reference Moretti, Mangeney, Walter, Capdeville, Bodin, Stutzmann and Le Friant2020; Vo et al. Reference Vo, Nezamabadi, Mutabaruka, Delenne and Radjai2020; Mandal, Nicolas & Pouliquen Reference Mandal, Nicolas and Pouliquen2020; Gans et al. Reference Gans, Abramian, Lagrée, Gong, Sauret, Pouliquen and Nicolas2023; Poulain et al. Reference Poulain2023).
A major contribution to the understanding of dense granular flows came with the collection of experimental data and discrete simulations of granular systems under simple shear, in particular by GDR MiDi (2004). In various configurations, they showed that the shear stress is proportional to the normal stress with a factor of proportionality, a friction 
$\mu = \mu (I)$, depending on a dimensionless quantity called the inertial number 
$I$. Small inertial numbers describe quasi-static flows where the macroscopic deformation is slow compared with the microscopic rearrangements, while large inertial numbers refer to rapid flows. The 
$I$-dependence of the friction was further generalized to tensorial form by Jop et al. (Reference Jop, Forterre and Pouliquen2006) by considering an analogous relation between the shear stress tensor and the isotropic pressure. Now known as the 
$\mu (I)$-rheology, it has become a well-established continuum theory for describing dense granular flows. Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) showed that the system of the incompressible 
$\mu (I)$-Navier–Stokes equations is well-posed unless the inertial number becomes too high or too small, and Barker & Gray (Reference Barker and Gray2017) later suggested partial regularization by altering the functional form of the 
$\mu (I)$-curve. Such ill-posedness was also observed numerically by Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015), Barker & Gray (Reference Barker and Gray2017) and Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017).
With regard to numerical implementations, Lagrée, Staron & Popinet (Reference Lagrée, Staron and Popinet2011) was the first to incorporate the 
$\mu (I)$-rheology equations in a two-dimensional finite volume Navier–Stokes model, and they applied it to simulate granular collapse as a continuum. Formally, the 
$\mu (I)$-rheology can be interpreted as a viscoplastic constitutive law involving a Drucker–Prager yield criterion being dependent on pressure and strain rate (Ionescu et al. Reference Ionescu, Mangeney, Bouchut and Roche2015; Martin et al. Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017). In such a viscoplastic model, Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) simulated granular collapse using an arbitrary Lagrangian–Eulerian finite element formulation of the Navier–Stokes equations. Relying on a two-dimensional smoothed particle hydrodynamics (SPH) model and a three-dimensional particle finite element method (PFEM), Chambon et al. (Reference Chambon, Bouvarel, Laigle and Naaim2011) and Franci & Cremonesi (Reference Franci and Cremonesi2019), respectively, demonstrated the implementation of the 
$\mu (I)$-rheology in such particle-based schemes, thus enabling the trivial handling of the free surface which can pose difficulties in mesh-based schemes. As a purely mesh-free scheme, SPH requires a computationally expensive neighbour search for each time step, and handling boundary conditions can be challenging (Vacondio et al. Reference Vacondio, Altomare, De Leffe, Hu, Le Touzé, Lind, Marongiu, Marrone, Rogers and Souto-Iglesias2021). Similarly, PFEM can also become expensive due to the need of creating a new mesh every time the finite element mesh becomes too distorted, although efficient remeshing algorithms may improve this (Cremonesi et al. Reference Cremonesi, Franci, Idelsohn and Oñate2020). Moreover, classical fluid solvers have issues in capturing material disconnections and with representing fully static zones, which become important features in the fluid-to-solid transition. It has been shown that regularization methods with a sufficiently small regularization parameter can reproduce the static–flowing transition (Frigaard & Nouar Reference Frigaard and Nouar2005; Lusso et al. Reference Lusso, Ern, Bouchut, Mangeney, Farin and Roche2017) as well as augmented Lagrangian methods (Martin et al. Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017). In a different approach, Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015) considered an 
$I$-dependent elastoplastic model with a Drucker–Prager yield in a two-dimensional material point method (MPM).
The ill-posedness of the 
$\mu (I)$-rheology for very small and large values of inertial numbers is an indication that important physics are missing in the model. In particular, for small inertial numbers, i.e. slow flows, concepts from solid mechanics like elasticity and shear banding become important (Barker et al. Reference Barker, Schaeffer, Bohorquez and Gray2015). While the 
$\mu (I)$-rheology has been successful in modelling granular liquids, critical state soil mechanics (CSSM) (Schofield & Wroth Reference Schofield and Wroth1968) is well established in the geo-mechanics community for describing the, potentially large, deformation of granular media as elastoplastic solids. Initially conceived to model soils and clay, this theory is based on the existence of a critical state, i.e. as a frictional fluid state where shear deformations can continue without any further change of volume, pressure or shear stress. All such critical states must reside along a so-called critical state line (CSL) in the space given by the pressure and shear stress. Among CSSM models, the modified cam-clay (MCC) model first formulated by Roscoe & Burland (Reference Roscoe and Burland1968) is arguably one of the most used in numerical implementations, giving an ellipsoidal yield surface in principal stress space. The CSSM models have been successfully applied to modelling the solid to fluid transition in snow avalanches (Gaume et al. Reference Gaume, Gast, Teran, van Herwijnen and Jiang2018, Reference Gaume, van Herwijnen, Gast, Teran and Jiang2019; Trottet et al. Reference Trottet, Simenhois, Bobillier, Bergfeld, van Herwijnen, Jiang and Gaume2022). They have also been used to simulate the avalanche dynamics over complex topography (Li et al. Reference Li, Sovilla, Jiang and Gaume2021, Reference Li, Sovilla, Gray and Gaume2024; Cicoira et al. Reference Cicoira, Blatny, Li, Trottet and Gaume2022). However, these CSSM models are based on a rate-independent theory, which can therefore not capture the 
$I$-dependent behaviour of granular systems in their liquid state.
The success of CSSM in capturing snow avalanches is, arguably, due to its ability to consider cohesion as well as plastic compaction and dilatancy. Cohesion refers to the presence of an attractive force between the grains, with the attractive force being either persistent or not. In fact, cohesion is key in understanding snow avalanches. The extensive snow experiments of Rognon et al. (Reference Rognon, Chevoir, Bellot, Ousset, Naaïm and Coussot2008a) showed that flowing dry snow behaves as dense granular flows with varying grain sizes, suggesting that this variation is a result of aggregates forming and size segregation within the flow, thus implying an important role of the cohesive forces in the snow. The cohesion has varying physical origins. In snow, it is a result of sintering of ice crystals or due to liquid water creating bonds by capillarity between snow particles (Szabo & Schneebeli Reference Szabo and Schneebeli2007; Steinkogler et al. Reference Steinkogler, Gaume, Löwe, Sovilla and Lehning2015; Ligneau, Sovilla & Gaume Reference Ligneau, Sovilla and Gaume2022). In fine powders, it may also be of electrostatic or magnetic origin (Castellanos Reference Castellanos2005).
In an effort to describe the rheology of cohesive granular flows, Rognon et al. (Reference Rognon, Roux, Wolf, Naaïm and Chevoir2006, Reference Rognon, Roux, Naaïm and Chevoir2008b) suggested on the basis of two-dimensional discrete simulations to make the coefficients in the linear 
$\mu (I)$-relation proposed by da Cruz et al. (Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005) dependent on a dimensionless number giving the ratio between the maximum inter-granular attractive force to the average normal force. Similar conclusions were made by Khamseh, Roux & Chevoir (Reference Khamseh, Roux and Chevoir2015) and Badetti et al. (Reference Badetti, Fall, Hautemayou, Chevoir, Aimedieu, Rodts and Roux2018) based on simulations and experiments of wet granular materials, showing that the apparent friction should depend both on the inertial number and the dimensionless cohesion number. Based on discrete simulations in both two dimensions and three dimensions, Berger et al. (Reference Berger, Azéma, Douce and Radjai2015) and Vo et al. (Reference Vo, Nezamabadi, Mutabaruka, Delenne and Radjai2020), respectively, proposed instead a new generalized inertial number dependent on the same dimensionless cohesion number.
Pertaining to the plastic compressibility of granular systems, in CSSM this is realized in a simple manner through a hardening law of the yield criterion. Specifically, in CSSM it is typically assumed that the inverse solid volume fraction decreases with the logarithm of pressure. While compressibility is not considered in the 
$\mu (I)$-rheology, a relation between the solid volume fraction 
$\phi$ and the inertial number was shown by GDR MiDi (2004), da Cruz et al. (Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005) and Pouliquen et al. (Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006). This prompted the compressible 
$\mu (I), \phi (I)$-rheology, with a monotonically decreasing, typically linear, function 
$\phi (I)$ for the solid fraction. Unfortunately, this 
$\mu (I), \phi (I)$-rheology is mathematically ill-posed in time-dependent problems (Heyman et al. Reference Heyman, Delannay, Tabuteau and Valance2017; Schaeffer et al. Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019).
The critical behaviour of granular flow was recently demonstrated in an extensive set of three-dimensional discrete simulations of granular systems by Wang, Li & Liu (Reference Wang, Li and Liu2022), showing that, under shear loading, a critical state is reached at sufficient shear for all shear rates, and that the exact critical state depends on the rate in a way that can be described by the inertial number. Previous depth-averaged models have introduced the critical state concept into the 
$\mu (I)$-rheology, with the first attempt made by Pailha & Pouliquen (Reference Pailha and Pouliquen2009) focusing on granular materials immersed in a secondary fluid. This study was augmented by Bouchut et al. (Reference Bouchut, Fernández-Nieto, Mangeney and Narbona-Reina2016) and Garres-Díaz et al. (Reference Garres-Díaz, Bouchut, Fernández-Nieto, Mangeney and Narbona-Reina2020) in a multi-layer model accounting for in-depth heterogeneities and later presented for the fully dry case by Bouchut et al. (Reference Bouchut, Fernández-Nieto, Koné, Mangeney and Narbona-Reina2021) in order to provide a compressible 
$\mu (I)$-rheology model, all studies successfully compared with laboratory experiments. While all of these previous models relied on depth-averaged equations in one or more layers based on shallow water assumptions, they also relied on the critical state theory of Roux & Radjai (Reference Roux and Radjai1998). In this viscoplastic approach, the dynamics is governed by a frictional parameter which can evolve depending on the deviation of the solid fraction from a prescribed critical solid fraction, thus different from the elastoplastic theory of CSSM. In an alternative approach, Rauter (Reference Rauter2021) supplements the otherwise linearly decreasing 
$\phi (I)$ in a Navier–Stokes-type model with the critical state model of Johnson & Jackson (Reference Johnson and Jackson1987) and Johnson, Nott & Jackson (Reference Johnson, Nott and Jackson1990), giving an expression for the pressure as a function of the shear rate and solid fraction.
The idea of combining CSSM and the 
$\mu (I)$-rheology was conceived by Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) and Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019), albeit in a purely viscoplastic form, resulting in equations they termed compressible 
$I$-dependent rheology. While not based on the elastoplastic framework of CSSM, compressibility was considered through the evolution of a (closed) yield surface and a carefully chosen plastic flow rule. Their motivation came from the specific desire of making the rheology well-posed for all inertial numbers through the introduction of compressibility. While succeeding in proving the well-posedness of their model subject only to reasonable restrictions, this conceptual and analytical study did not offer a complete numerical realization beyond simple one-dimensional demonstrations, nor an experimental comparison or validation.
In an effort to capture the behaviour of granular media through both their solid and fluid states in a continuum model, it seems natural to seek a model that combines CSSM, which readily incorporates cohesion and compressibility, with the 
$\mu (I)$-rheology, allowing the inertial number dependence in flows. Now, guided by advances in particle-based numerical schemes and concepts from computational elastoplasticity, we are able to formulate and implement an elasto-viscoplastic MCC model that follows the 
$\mu (I)$-rheology under flow conditions. This allows the modelling of cohesive and compressible granular materials as elastoplastic solids which, in the fluid transition, recovers the expected 
$I$-dependent behaviour. In short, this is accomplished by making the CSL dependent on the inertial number and by relying on MPM which, due to its hybrid Eulerian–Lagrangian nature, has the ability to handle large deformations as the material transitions from solid-like to fluid-like, including material separation and reconsolidation. While MPM also offers the advantage of not requiring any special numerical treatment of the free surface, it comes with increasing computational costs (compared with classical finite element methods) as well as the need for a careful treatment of numerical dissipation. The type of cohesion treated in the model is one that does not cease to exist, e.g. in order to consider the situation of permanent attractive inter-granular force. However, as will be discussed, other types of cohesion could be considered in the proposed general framework that would make the model more suitable for landslides and debris flows. As an elastoplastic model implemented in MPM, the proposed model is inspired by the work of Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015), although differing both in the elastic and plastic equations. In addition to considering cohesion and compressibility through the critical state, the proposed model relies on advancements from the MPM community which improve stability and reduce numerical dissipation. In this paper, the model is demonstrated for collapse and flow problems where the magnitude of the elastic deformations is negligible compared with the (rate-dependent) plastic deformations, and the general solid behaviour and elastic nature of granular media will not be treated in detail here.
2. Theory
The granular material is here considered as a continuum governed by the conservation laws for mass and momentum. Conservation of mass is given by
\begin{equation} \frac{{\rm D} \rho}{{\rm D} t} + \rho(\boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{v}) = 0, \end{equation}and conservation of momentum under gravity states that
\begin{equation} \rho \frac{{\rm D} \boldsymbol{v}}{{\rm D} t} - \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{\sigma}_0 - \rho \boldsymbol{g} = 0 , \end{equation}
where 
$\rho = \rho (\boldsymbol {x},t)$ is the mass density, 
$\boldsymbol {v} = \boldsymbol {v}(\boldsymbol {x},t)$ is the velocity, 
$\boldsymbol {\sigma }_0 = \boldsymbol {\sigma }_0(\boldsymbol {x},t)$ is the symmetric Cauchy stress tensor, 
$\boldsymbol {g}$ is the constant gravitational acceleration and 
${\rm D}/{\rm D}t$ denotes the material time derivative. Without further assumptions, the system of equations given by these conservation laws is not closed. In this section, constitutive laws relating the deformation to the stress are outlined in order to close the system. In particular, the granular material will be treated as an elastoplastic continuum where the plastic, i.e. irreversible, deformations are highly rate dependent.
2.1. Definitions of deformation
With the deformation of a body over a time 
$t$ described by a one-to-one mapping from a material (initial) coordinate 
$\boldsymbol {X}$ to a spatial (current) coordinate 
$\boldsymbol {x}$, the deformation gradient 
$\boldsymbol {F}$ is defined with 
$F_{ij} = \partial x_i / \partial X_j$ and 
$\det (\boldsymbol {F}) > 0$. The singular values of this second-order invertible tensor are the principal stretches 
$\lambda _\alpha > 0$, 
$\alpha =1,\ldots,\dim$, where 
$\dim =1,2,3$ denotes the number of spatial dimensions. From the deformation gradient, the Hencky strain tensor is defined as
\begin{equation} \boldsymbol{\varepsilon}= \frac{1}{2} \ln \boldsymbol{F} \boldsymbol{F}^T = \sum_{\alpha=1}^{\dim} \ln \lambda_\alpha\ \boldsymbol{n}_\alpha \otimes \boldsymbol{n}_\alpha, \end{equation}
where the unit vectors 
$\boldsymbol {n}_\alpha$ define the spatial principal directions. Moreover, the velocity gradient tensor can be written as
\begin{equation} \boldsymbol{L} = \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{v} = \dot{\boldsymbol{F}} \boldsymbol{F}^{{-}1}, \end{equation}
where 
$\dot {\boldsymbol {F}}$ denotes the material time derivative of the deformation gradient. The trace of the Hencky strain, 
$\varepsilon _V = {\rm tr}\,\boldsymbol {\varepsilon } = \ln ( \det (\boldsymbol {F}) )$, is a measure of the volumetric deformation, with material time derivative 
$\dot {{\varepsilon }}_V = {\rm tr}(\boldsymbol {L})$.
2.2. Elasticity and yield surface
To account for both elastic, reversible, and plastic, irreversible, deformations, it is common to consider the multiplicative decomposition of the deformation gradient
\begin{equation} \boldsymbol{F} = \boldsymbol{F}^E \boldsymbol{F}^P, \end{equation}
where 
$\boldsymbol {F}^E$ refers to the deformation arising from the elastic forces, while 
$\boldsymbol {F}^P$ is the permanent, plastic, component. We seek a law that relates the stresses (due to elastic forces) to the elastic component of the Hencky strain
\begin{equation} \boldsymbol{\varepsilon}^E = \tfrac{1}{2} \ln ( \boldsymbol{F}^E (\boldsymbol{F}^E)^T ) {,} \end{equation}
through a so-called strain energy density function 
$\varPsi$, i.e. a hyperelastic law. In this work, a frame-indifferent isotropic hyperelastic model known as the Hencky model is adopted, where
\begin{equation} \varPsi(\boldsymbol{\varepsilon}^E) = \tfrac{1}{2} \lambda ({\rm tr}\,\boldsymbol{\varepsilon}^E)^2 + G \,{\rm tr}((\boldsymbol{\varepsilon}^E)^2), \end{equation}
which results in the Kirchhoff stress 
$\boldsymbol {\sigma } = \det (\boldsymbol {F}) \boldsymbol {\sigma }_0$ being related to the Hencky strain through
\begin{equation} \boldsymbol{\sigma} = \frac{\partial \varPsi(\boldsymbol{\varepsilon}^E)}{\partial \boldsymbol{\varepsilon}^E} = \lambda ({\rm tr}\,\boldsymbol{\varepsilon}^E) \boldsymbol{I} + 2 G \boldsymbol{\varepsilon}^E , \end{equation}
where 
$\boldsymbol {I}$ denotes the second-order identity tensor, and 
$\lambda$ and 
$G$ are the two Lamé parameters which may be related to Young's modulus 
$E$ and Poisson's ratio 
$\nu$. Equation (2.8) may also be written
\begin{equation} \boldsymbol{\sigma} = \mathcal{C} : \boldsymbol{\varepsilon}^E, \end{equation}
with the fourth-order isotropic stiffness matrix 
$\mathcal {C}$ well known from linear elasticity with elements 
$\mathcal {C}_{ijkl} = \lambda \delta _{ij}\delta _{kl} + G(\delta _{ik}\delta _{jl} + \delta _{il}\delta _{jk})$. The properties of this model are elegantly demonstrated in Xiao & Chen (Reference Xiao and Chen2002), showing that, with this model, the Kirchhoff stress and the Hencky strain constitute a work–conjugate pair.
It is instructive to relate Hencky's elasticity model to hypoelastic models where a widely used relation between an objective stress rate 
$\mathring {\boldsymbol {\sigma }}$ and the elastic rate-of-deformation 
$\boldsymbol {D}^E = \frac {1}{2}(\boldsymbol {L}^E+(\boldsymbol {L}^E)^T)$ can be written as
\begin{equation} \mathring{\boldsymbol{\sigma}} = \mathcal{C} : \boldsymbol{D}^E {.} \end{equation}
Choosing the Jaumann stress rate, this is similar to the rate equation used in the elastoplastic 
$\mu (I)$-rheology model of Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015) except they considered the Cauchy stress. Equation (2.10) should be exactly integrable to give an elastic relation. If this so-called self-consistency criterion is to be fulfilled, then it can be shown that the logarithmic rate of the Kirchhoff stress is the only choice among all corotational rates (Xiao, Bruhns & Meyers Reference Xiao, Bruhns and Meyers1999). To go even further, the logarithmic rate is argued by Xiao, Bruhns & Meyers (Reference Xiao, Bruhns and Meyers2000) as the only choice among corotational and non-corotational rates in the context of elastoplasticity. Interestingly, the integration of (2.10) with the logarithmic rate assuming a stress-free initial state yields exactly Hencky's elasticity model (Bruhns, Xiao & Meyers Reference Bruhns, Xiao and Meyers1999; Xiao & Chen Reference Xiao and Chen2002).
The onset of plastic deformations is characterized by a yield function 
$y = y(\boldsymbol {\sigma })$. Conventionally, 
$y \leqslant 0$ defines admissible states, and the yield surface 
$y=0$ defines the limit where deformations are no longer elastic. Thus, plastic deformations only occur on the yield surface. Here, the yield function will be expressed in terms of the stress invariants
\begin{equation} p ={-}\frac{1}{\dim} \,{\rm tr}\,{\boldsymbol{\sigma}} \quad\textrm{and}\quad q = \frac{1}{\sqrt{2} } || \boldsymbol{\tau} || , \end{equation}
where 
$p$ is termed the isotropic pressure and 
$q$ the equivalent shear stress. The latter was defined from the deviatoric stress tensor
\begin{equation} \boldsymbol{\tau} = \operatorname{dev}(\boldsymbol{\sigma}) = \boldsymbol{\sigma} - \frac{1}{\dim}\,{\rm tr}(\boldsymbol{\sigma})\boldsymbol{I}, \end{equation}
using the notation 
$|| \boldsymbol {\tau } || = \sqrt {\boldsymbol {\tau } : \boldsymbol {\tau }}$.
While yielding of many granular media, including sand and systems of glass beads, follows a Drucker–Prager criterion, other granular media, e.g. snow, are best described by a yield function providing a closed yield surface in stress space (Reiweger, Gaume & Schweizer Reference Reiweger, Gaume and Schweizer2015; Ritter, Gaume & Löwe Reference Ritter, Gaume and Löwe2020; Blatny et al. Reference Blatny, Löwe, Wang and Gaume2021). To this end, we consider the yield criterion given by the (cohesive) MCC model
\begin{equation} y(\boldsymbol{\sigma}) = y(p,q) = q^2 - \mu^2 (p+\beta p_c) (p_c-p) , \end{equation}
where 
$p_c \geqslant 0$ represents an isotropic compressive strength, 
$\beta \geqslant 0$ is a dimensionless measure of cohesion and 
$\mu >0$ defines the slope of the CSL along which all critical states reside. Alternatively, one may say that all stress states along the CSL are on a Drucker–Prager yield surface 
$q=\mu (p+\beta p_c)$. The MCC yield surface is illustrated in figure 1.
The MCC yield surface 
$y(p,q)=0$ from (2.13) in the space of the stress invariants 
$p$ and 
$q$ defined in (2.11a,b). Elastic states satisfy 
$y<0$. The dashed line represents the CSL. The dash-dotted line marks the value of 
$p$ at the apex of the ellipse, separating states undergoing plastic compaction and dilation as a result of the chosen associative plastic flow rule.
In this model, 
$-\beta p_c$ represents the isotropic tensile strength. The case 
$\beta = 0$ refers to the cohesionless case, and in this case the material cannot sustain tensile stress states. For 
$\beta > 0$, the system will always have a non-zero tensile strength (unless 
$p_c = 0$, i.e. the stress-free state). As such, cohesion will always be present, which may not be the case in, e.g. landslides. For granular systems in which the cohesion is better described by breakable bonds, 
$\beta$ may easily be set to zero after yield, however, such a situation will not be considered here.
In order to complete the description of the model, a plastic flow rule will be outlined in § 2.3, consequently determining the onset and magnitude of plastic dilation and compaction. In § 2.4, a hardening law is introduced that governs the evolution of 
$p_c$ and thus controls the, possibly extensive, variations of solid fraction. Finally, § 2.5 introduces dependence on the inertial number.
2.3. Plastic flow rule
Under the assumption that the rate of work per unit undeformed volume can be additively decomposed into an elastic and plastic part, it follows that
\begin{equation} \boldsymbol{\sigma} : \boldsymbol{L} = \boldsymbol{\sigma} : \boldsymbol{L}^E + \boldsymbol{\sigma} : \boldsymbol{L}^P, \end{equation}
where an elastic and a plastic part of the velocity gradient were defined. Following conventions, the plastic part 
$\boldsymbol {L}^P$ is chosen as the symmetric tensor
\begin{equation} \boldsymbol{L}^P = \dot{\gamma} \frac{\partial g / \partial \boldsymbol{\sigma} }{||\partial g / \partial \boldsymbol{\sigma}||} ,\end{equation}
where the equivalent plastic strain 
$\dot {\gamma } \geqslant 0$ is defined as 
$\dot {\gamma } = || \boldsymbol {L}^P ||$ and 
$g$ is the so-called plastic potential. In MCC models, the typical choice is an associative plastic flow rule, i.e. 
$g=y$ (the plastic potential is associated with the yield function). In other words, the plastic rate of deformation coincides with the gradient of the yield surface. An associative flow rule is a result of the principle of maximum plastic dissipation, in particular it can be shown that maximizing 
$\boldsymbol {\sigma } : \boldsymbol {L}^P$ subject to 
$y \leqslant 0$ leads to (2.15) where 
$\dot {\gamma }$ can be interpreted as a Lagrange multiplier (Simo Reference Simo1992; Bonet & Wood Reference Bonet and Wood2008).
With 
$g = y(\boldsymbol {\sigma }) = y(p,q)$, (2.15) can be written
\begin{equation} \boldsymbol{L}^P = \frac{1}{\sqrt{2}} \dot{\gamma} \frac{\partial y / \partial q }{||\partial y / \partial \boldsymbol{\sigma}||} \frac{ \boldsymbol{\tau}}{|| \boldsymbol{\tau} ||} - \frac{1}{\dim} \dot{\gamma} \frac{\partial y / \partial p}{||\partial y / \partial \boldsymbol{\sigma}||} \boldsymbol{I} {,} \end{equation}promoting the definition of the equivalent plastic shear strain rate as
\begin{equation} \dot{\gamma}_{_S} = \dot{\gamma} \frac{\partial y / \partial q }{||\partial y / \partial \boldsymbol{\sigma}||} = \sqrt{2} || \operatorname{dev}(\boldsymbol{L}^P) || , \end{equation}and the equivalent plastic volumetric strain rate as
\begin{equation} \dot{\gamma}_{_V} ={-} \dot{\gamma} \frac{\partial y / \partial p}{||\partial y / \partial \boldsymbol{\sigma}||} = {\rm tr}(\boldsymbol{L}^P) , \end{equation}
such that the plastic volumetric Hencky strain is 
$\varepsilon _V^P = \ln ( \det (\boldsymbol {F}^P) ) = \int \dot{\gamma}_{_V} \,{\rm d}t$. The associative plastic flow rule implies a plastic volume increase for stress states with 
$p < p_{cs}$ and a volume decrease for 
$p > p_{cs}$, where 
$p_{cs} = \frac {1}{2}p_c(1-\beta )$ denotes the pressure at the apex of the yield surface (marked by the dash-dotted vertical line in figure 1). Note that the CSL always goes through the apex, so there should indeed be no volume change for stress states at critical state.
At this point, it can be instructive to relate the current model to the critical state model of Roux & Radjai (Reference Roux and Radjai1998) which has been used in the depth-averaged (single- or multi-layer) models of Pailha & Pouliquen (Reference Pailha and Pouliquen2009), Bouchut et al. (Reference Bouchut, Fernández-Nieto, Mangeney and Narbona-Reina2016, Reference Bouchut, Fernández-Nieto, Koné, Mangeney and Narbona-Reina2021) and Garres-Díaz et al. (Reference Garres-Díaz, Bouchut, Fernández-Nieto, Mangeney and Narbona-Reina2020). Assuming for simplicity the cohesionless case, the stress states in that model are forced to be on the (Drucker–Prager) surface 
$q=\tan (\delta + \psi ) p$, where 
$\tan (\delta )=\mu$ is the critical friction coefficient and 
$\psi$ is a so-called dilatancy angle which satisfies 
$\psi > 0$ under dilation and 
$\psi < 0$ under compaction. If one imagines 
$q=\tan (\delta ) p$ intersecting an MCC-ellipsoid at its apex where 
$p = p_{cs}$, it follows that 
$q=\tan (\delta + \psi ) p$ must intersect the same ellipsoid at 
$p < p_{cs}$ under dilation and at 
$p > p_{cs}$ under compaction.
The choice of plastic flow rule (or plastic potential) is fundamental to the possibility and extent to which the material can undergo volumetric deformations, and accordingly to how critical state is reached. In Drucker–Prager/Mohr–Coulomb models, an associative flow rule would imply every plastic deformation being dilating. The use of non-associative flow rules in such models prevents this. In addition, the singularities in the Drucker–Prager/Mohr–Coulomb yield surfaces require special treatment which often may induce unwanted dilation, although they can be corrected with various strategies (Dunatunga & Kamrin Reference Dunatunga and Kamrin2015; Pradhana Reference Pradhana2017).
2.4. Hardening law
In soil mechanics, the isotropic compressive strength 
$p_c$ is typically increased with plastic compaction (i.e. hardening – yield surface expands) and decreased with plastic dilation (i.e. softening – yield surface shrinks). As such, 
$p_c = p_c(\varepsilon _V^P)$ which entails that 
$p_c$ will not change at critical state. Here, a similar relation will be assumed, specifically
\begin{equation} p_c(\varepsilon_V^P) = p_c^0 {e}^{-\xi \varepsilon_V^P}, \end{equation}
where 
$p_c^0 \geqslant 0$ is an initial compressive strength before any plastic compaction/dilation has occurred, and 
$\xi \geqslant 0$ is a dimensionless parameter.
This hardening law can be motivated from the observed behaviour of soils (Wood Reference Wood1991), cohesive and cohesionless powders (Walker Reference Walker1923; Poquillon et al. Reference Poquillon, Lemaitre, Baco-Carles, Tailhades and Lacaze2002; Castellanos Reference Castellanos2005) and porous brittle media such as snow (Blatny, Löwe & Gaume Reference Blatny, Löwe and Gaume2023), where the inverse solid volume fraction behaves as
\begin{equation} \frac{1}{\phi} = \frac{1}{\phi_0} - \varLambda \ln(p_c/p_c^0), \end{equation}
where 
$\varLambda >0$ is called the plasticity index and 
$\phi _0$ is the initial solid fraction. In (2.20), the compressive strength is 
$p_c = p_c^0$ when the solid fraction is 
$\phi =\phi _0$, with 
$p_c^0$ and 
$\phi _0$ being two independent parameters. In reality, 
$p_c^0$ and 
$\phi _0$ may be related, with typically an increasing 
$p_c^0$ for an increasing 
$\phi _0$. However, in this study, systems with varying 
$\phi _0$ will not be investigated.
Equation (2.20) can be rearranged as
\begin{equation} \frac{\phi_0}{\phi} - 1 ={-} \phi_0 \varLambda \ln(p_c/p_c^0). \end{equation}Here, the left-hand side may be approximated by the volumetric plastic Hencky strain
\begin{equation} \varepsilon_V^P = \det(\boldsymbol{F}^P) \stackrel{(i)}{\approx} \det(\boldsymbol{F}) = \ln(\phi_0/\phi) \stackrel{(ii)}{\approx} \frac{\phi_0}{\phi} - 1, \end{equation}under the assumption that (i) the elastic volumetric deformation is negligible compared with the plastic volumetric deformation and that (ii) the volumetric deformation is small. This allows (2.21) to be written as
\begin{equation} p_c = p_c^0 {e}^{-\varepsilon_V^P / (\varLambda \phi_0) }, \end{equation}
which is exactly (2.19) with 
$\xi = (\phi _0 \varLambda )^{-1}$. Note that the volumetric plastic Hencky strain 
$\varepsilon _V^P$ is zero before any deformation occurred (at which point 
$\phi =\phi _0$ and 
$p_c = p_c^0$).
2.5. Inertial number dependence
In the 
$\mu (I)$-rheology, the shear stress is made proportional to the isotropic pressure with a factor of proportionality 
$\mu$ depending on a dimensionless number known as the inertial number 
$I$. Assuming the deformation is dominantly plastic, the inertial number can be written
\begin{equation} I=\frac{\dot{\gamma}_{_S} d}{\sqrt{p /\rho_*}}, \end{equation}
where 
$d$ denotes the diameter of the grains and 
$\rho _*$ is the constant intrinsic grain density.
In the MCC model outlined in §§ 2.2–2.4, it is only at critical state that 
$\mu$ represents the factor of proportionality between the equivalent shear stress and the pressure. As most stress states in an unobstructed smooth flow would gather close to the CSL, the 
$\mu (I)$-rheology would be approximately recovered in such flows by letting 
$\mu$ depend on 
$I$ in the same way as in the 
$\mu (I)$-rheology. A typical expression for 
$\mu (I)$ is such that it is bound by a lower value 
$\mu _1$ and approaching an upper value 
$\mu _2$ asymptotically with increasing inertial number, in particular
\begin{equation} \mu(I) = \mu_1 + \frac{\mu_2-\mu_1}{\dfrac{I_0}{I}+1}, \end{equation}
from Jop, Forterre & Pouliquen (Reference Jop, Forterre and Pouliquen2005), where 
$I_0 > 0$ is a dimensionless constant.
Evidently, (2.24) is not well defined for negative pressures, which is possible in the cohesive case, i.e. when 
$\beta >0$. As a remedy, a new pressure 
$\bar {p} = p + \beta p_c$ is introduced, where the cohesive stress is added to the pressure, leading to the definition of a cohesion-dependent inertial number
\begin{equation} I=\frac{\dot{\gamma}_{_S} d}{\sqrt{\bar{p} /\rho_*}} {.} \end{equation}
With 
$\beta p_c>0$ here representing the magnitude of the isotropic cohesive stress (tensile strength), this cohesive inertial number was proposed by Berger et al. (Reference Berger, Azéma, Douce and Radjai2015) and Vo et al. (Reference Vo, Nezamabadi, Mutabaruka, Delenne and Radjai2020). They demonstrated, through discrete simulations, that this inertial number improved the scaling behaviour for granular flows where cohesion was added to the system through an (unbreakable) attractive inter-granular force. This is in contrast to other studies (Rognon et al. Reference Rognon, Roux, Wolf, Naaïm and Chevoir2006, Reference Rognon, Roux, Naaïm and Chevoir2008b; Khamseh et al. Reference Khamseh, Roux and Chevoir2015; Badetti et al. Reference Badetti, Fall, Hautemayou, Chevoir, Aimedieu, Rodts and Roux2018) which have proposed making the apparent friction dependent on a separate dimensionless cohesion number in addition to the (non-cohesive) inertial number. Note that the non-cohesive inertial number, (2.24), is related to its cohesive counterpart, (2.26), through a multiplication of 
$\sqrt {1+\beta p_c/p}$. With this inertial number, (2.25) becomes
\begin{equation} \mu(I) = \mu_1 + \frac{\mu_2-\mu_1}{\omega\dfrac{\sqrt{\bar{p}}}{\dot{\gamma}_{_S}}+1}, \end{equation}
where, for convenience, the constant 
$\omega = {I_0}/({d\sqrt {\rho _*}})$ was introduced. The evolution of 
$\mu$ with 
$I$ when changing 
$\omega$ is sketched in figure 2.
Evolution of 
$\mu (I)$ from (2.27) as the parameter 
$\omega$ is increased.
In a nutshell, the proposed model consists of the elastic law given by (2.8), the yield function given by (2.13) with 
$\mu$ being the function defined in (2.27), the plastic flow rule given by (2.15) and the hardening law given by (2.19). A summary of the model parameters is presented in table 1, where it should be noted that 
$I_0$, 
$d$ and 
$\rho _*$ appear in the proposed model through 
$\omega$ only.
Model parameters.
Formally, this elasto-viscoplastic model can be written as a Perzyna overstress model (Perzyna Reference Perzyna1963, Reference Perzyna1966), in particular
\begin{equation} {\dot{\gamma}_{_S}} = \begin{cases} \dfrac{1}{\eta} \left( q-q_y(p) \right), & \textrm{if } q > q_y(p) \\ 0, & \textrm{otherwise}, \end{cases} \end{equation}
with a 
$p$-dependent yield shear stress
\begin{equation} q_y(p) = \mu_1 \sqrt{ (p_c-p) \bar{p}} \geqslant 0, \end{equation}
and an 
$I$-dependent viscosity
\begin{equation} \eta = \frac{ (\mu(I)-\mu_1) \sqrt{ (p_c-p) \bar{p}} }{ \dot{\gamma}_{_S} } = \frac{(\mu_2-\mu_1) \sqrt{ (p_c-p) \bar{p}} }{\dot{\gamma}_{_S} + \omega \sqrt{\bar{p}}} \geqslant 0 {.} \end{equation}
The Perzyna yield surface given by 
$q_y(p)$ differs from the previous concept of a yield surface as stress states outside the Perzyna yield surface are admissible.
2.6. An operator splitting scheme
Solving for the stress (and deformation) can be accomplished through an elastic predictor–plastic corrector scheme (Simo Reference Simo1992; de Souza Neto, Perić & Owen Reference de Souza Neto, Perić and Owen2008). While the details of this operator splitting scheme are presented in Appendix A, the concept is simple: assume first the deformation is purely elastic, and if the yield condition is violated, correct for any plastic deformation through the plastic flow rule. In particular, a trial stress state 
$(p^{trial}, q^{trial})$ obtained elastically from the previous stress state 
$(p^n,q^n)$ at time 
$t^n$ gives the new stress state 
$(p^{n+1}, q^{n+1})$ after a plastic correction
$$\begin{gather} p^{n+1} = p^{trial} + K \dot{\gamma}_{_V} \Delta t \end{gather}$$
$$\begin{gather}q^{n+1} = q^{trial} - G \dot{\gamma}_{_S} \Delta t, \end{gather}$$
where 
$\Delta t = t^{n+1} - t^n$. This operation is known as a return mapping.
Solving the return mapping must be done iteratively, the convergence of which is difficult to guarantee. As an approximation, since the hardening law 
$p_c = p_c(\varepsilon _V^P)$ is only dependent on the update of 
$p$, the return mapping can be performed in the view of Perzyna, assuming 
$\mu = \mu _1$ is constant (cf. (2.29)), producing the desired 
$p^{n+1}$. With 
$p^{n+1}$ now known, 
$q^{n+1}$ can be found by inserting (2.32) into (2.13) and using (2.27), yielding a quadratic equation for 
$\dot{\gamma}_{_S}$
\begin{equation} a \dot{\gamma}^2_{_S} + b \dot{\gamma}_{_S} + c = 0, \end{equation}where
\begin{equation} \left.\begin{gathered} a = G \Delta t \\ b = (\mu_2-\mu_1) \sqrt{ (p_c-p^{n+1}) \bar{p}^{n+1} } + G \Delta t \omega \sqrt{\bar{p}^{n+1}} - (q^{trial} - q_y(p^{n+1})) \\ c ={-}\omega\sqrt{\bar{p}^{n+1}} (q^{trial} - q_y(p^{n+1})) \end{gathered}\right\}. \end{equation}
Since 
$a>0$ and 
$c<0$, exactly one positive solution for 
$\dot{\gamma}_{_S}$ is guaranteed, and the new stress state 
$(p^{n+1}, q^{n+1})$ is found. Notably, there are no issues with zero pressure or strain rate. The return mapping is illustrated in figure 3.
The return mapping in the view of Perzyna's overstress approach. A trial state 
$(p^{trial}, q^{trial})$ is projected to the stress state 
$(p^{n+1},q^{n+1})$ through an intermediate step 
$(p^{n+1}, q_y(p^{n+1}))$. Note that, since here 
$p^{trial} > p_{cs}$, volumetric compaction occurred and the yield surface therefore expanded (hardened) from the initial surface coloured grey to the final surface coloured black.
3. Numerical solver
Here, MPM is adopted as a numerical scheme for approximating solutions to (the weak form of) (2.2) subject to the elasto-viscoplastic constitutive model. In short, MPM consists of tracking moving particles (or ‘material points’) which store associated material properties such as mass, velocity and deformation gradient, while also adopting a background Eulerian grid facilitating the momentum computations. This is sketched in figure 4. Although it is commonly stated that the original formulation of MPM dates back to Sulsky, Chen & Schreyer (Reference Sulsky, Chen and Schreyer1994), MPM can be considered an extension of particle-in-cell (PIC) methods (Harlow Reference Harlow1964), notably the fluid-implicit particle (FLIP) method (Brackbill, Kothe & Ruppel Reference Brackbill, Kothe and Ruppel1988). It was first applied to the flow of granular media by Wieckowski, Youn & Yeon (Reference Wieckowski, Youn and Yeon1999) in the context of silo discharge, assuming simply a cohesionless constant Drucker–Prager yield criterion.
Illustrating sketch of the MPM discretization (grid and particles) of a flow on a surface subject to gravity. In this study, the domain is initialized with four to eight MPM particles per grid cell.
The material domain is discretized into 
$N_p$ particles with constant and equal mass 
$m_p$, initially equal volume 
$V_p$ and initially zero velocity 
$\boldsymbol {v}_p$. With a fixed Eulerian grid with uniform grid spacing 
$\Delta x$, we consider 
$C^1$-continuous interpolation functions 
$N_{i}(\boldsymbol {x}_p) \geqslant 0$ between a grid node 
$i$ and a particle at a coordinate 
$\boldsymbol {x}_p$, in particular
\begin{equation} N_i(\boldsymbol{x}_p) = \sum_{\alpha=0}^{\dim} N \left( \frac{x_{p \alpha}-x_{i \alpha}}{\Delta x} \right) , \end{equation}
where 
$\alpha$ denotes the component and with 
$N( {\cdot } )$ given by the quadratic B-spline
\begin{equation} N(x) =\begin{cases} \frac{3}{4} - |x|^2, & \textrm{if } |x| < \frac{1}{2} \\[6pt] \frac{1}{2}(\frac{3}{2} - |x|)^2, & \textrm{if } \frac{1}{2} \leqslant |x| < \frac{3}{2} \\[6pt] 0, & \textrm{otherwise.} \end{cases} \end{equation}
Such B-spline MPM circumvents the issue of cell crossing instabilities (Steffen, Kirby & Berzins Reference Steffen, Kirby and Berzins2008) associated with the original MPM which relied on linear hat functions (e.g. used in the elastoplastic 
$\mu (I)$-rheology model of Dunatunga & Kamrin Reference Dunatunga and Kamrin2015). Not considering the boundary term, the temporal and spatial discretization of the weak form of the momentum conservation equation can be written on the Eulerian grid nodes 
$i$ as
\begin{equation} \frac{ m^n_{i} \boldsymbol{v}^{n+1}_{i} - m^n_{i} \boldsymbol{v}^{n}_{i} }{\Delta t} ={-} \sum_p V^0_{p} \boldsymbol{\sigma}(\boldsymbol{\varepsilon}^{E,n}) \boldsymbol{\nabla} N_{i}(\boldsymbol{x}^n_p) + m_i^n \boldsymbol{g} ,\end{equation}
where 
$V^0_p$ is the initial particle volume and the nodal mass
\begin{equation} m_i^n = \sum_p m_p N_i(\boldsymbol{x}_p^n)\end{equation}
is obtained under the common assumption of mass lumping. The MPM discretization from the weak form can be thought of as the equivalent (updated Lagrangian) finite element discretization where now the quadrature points are the material points which can move in space. Note that conservation of mass is automatically fulfilled as the total mass of the system 
$\sum _{p} m_p$ is the sum of all particle masses, which individually will remain constant in time.
The detailed steps of the numerical scheme are outlined in Algorithm 1. At each time step of the scheme, particle momentum is transferred to the grid on which the momentum update, (3.3), is performed. Boundary conditions are applied on the grid velocities, before the new grid velocities are transferred back to the particles. As such, boundary conditions are treated in a similar way as in finite element methods, avoiding expensive neighbour particle searches needed in, e.g. SPH. Before commencing the next time step, the particles’ stress states are checked against a yield criterion, and a return mapping is performed if necessary. Adaptive time stepping in this explicit scheme is accomplished by respecting the Courant–Friedrichs–Lewy (CFL) condition and the elastic wave speed condition (i.e. the elastic wave should not travel more than one grid cell in one time step). Accordingly, we use the constants 
$C_{CFL} = C_{el} = 0.5$, where the elastic wave constant 
$C_{el}$ is defined analogously to the CFL constant 
$C_{CFL}$.
Symplectic B-spline AFLIP MPM
It should be emphasized that there are several methods for transferring the momentum/velocities between the grid nodes and particles, and the chosen method can greatly influence the numerical dissipation of the scheme. With the time-dependent constitutive model proposed here, it is crucial to minimize this dissipation in order to capture the dynamics of the deformation correctly. Earlier versions of MPM used PIC- or FLIP-style transfers, i.e.
\begin{equation} m^n_{i} \boldsymbol{v}^n_i = \sum_{p} m_p \boldsymbol{v}^n_p N_{i}(\boldsymbol{x}^n_p) \end{equation}for the particle-to-grid momentum transfer and either (PIC)
\begin{equation} \boldsymbol{v}_p^{n+1} = \sum_i \boldsymbol{v}_i^{n+1} N_{i}(\boldsymbol{x}_p^n), \end{equation}or (FLIP)
\begin{equation} \boldsymbol{v}_p^{n+1} = \boldsymbol{v}_p^n + \sum_{i} (\boldsymbol{v}_i^{n+1}-\boldsymbol{v}_i^n) N_i(\boldsymbol{x}_p^n) ,\end{equation}
for the grid-to-particle transfer. The PIC scheme is highly dissipative and does not preserve angular momentum in grid to particle, which could lead to rotational artefacts. While FLIP (as used, e.g. in the elastoplastic 
$\mu (I)$-rheology model of Dunatunga & Kamrin Reference Dunatunga and Kamrin2015) partly improves the conservation of angular momentum, it also reduces dissipation by only transferring the velocity increments, however, at the cost of more pronounced ‘ringing instability’ (Love & Sulsky Reference Love and Sulsky2006). Later, Jiang et al. (Reference Jiang, Schroeder, Selle, Teran and Stomakhin2015) proposed affine PIC (APIC) as an improvement of PIC, conserving angular momentum in both transfers, consequently resolving rotations correctly. The APIC particle-to-grid transfer is given by
\begin{equation} m^n_{i} \boldsymbol{v}^n_i = \sum_p m_p ( \boldsymbol{v}_p^n + \boldsymbol{C}_p^n (\boldsymbol{x}_i^n - \boldsymbol{x}_p^n) ) N_{i}(\boldsymbol{x}^n_p), \end{equation}
with a ‘velocity derivative’ tensor 
$\boldsymbol {C}_p = \boldsymbol {B}_p \boldsymbol {D}_p^{-1}$ where 
$\boldsymbol {D}_p$ is analogous to an inertia tensor and can be shown to be a constant diagonal tensor given by 
$\frac {1}{4} (\Delta x)^2 \boldsymbol {I}$ in the case of a uniform background grid with quadratic B-splines (Jiang et al. Reference Jiang, Schroeder, Selle, Teran and Stomakhin2015; Jiang, Schroeder & Teran Reference Jiang, Schroeder and Teran2017), and 
$\boldsymbol {B}_p$ is given by
\begin{equation} \boldsymbol{B}_p^{n+1} = \sum_i \boldsymbol{v}^{n+1}_i (\boldsymbol{x}^n_i - \boldsymbol{x}^n_p)^T N_{i}(\boldsymbol{x}_p^n) . \end{equation}
In an effort to reduce dissipation further, APIC was recently generalized to AFLIP by Fei et al. (Reference Fei, Guo, Wu, Huang and Gao2021) by using the grid-to-particle scheme of FLIP, i.e. (3.7), instead of that of PIC, i.e. (3.6). Because of the low-dissipation property of AFLIP, this scheme is adopted in this work. Like in Fei et al. (Reference Fei, Guo, Wu, Huang and Gao2021), a parameter 
$\alpha _f$ is introduced which represents a PIC/FLIP ratio, with 
$\alpha _f=0$ meaning pure APIC and 
$\alpha _f=1$ meaning pure AFLIP. In particular, 
$\alpha _f=0.95$ is used to introduce as little numerical dissipation as possible. More details and examples of the low-dissipation property of AFLIP are presented in the paper of Fei et al. (Reference Fei, Guo, Wu, Huang and Gao2021), and the reader is also referred to the excellent and recent reviews by de Vaucorbeil et al. (Reference de Vaucorbeil, Nguyen, Sinaie and Wu2020), Sołowski et al. (Reference Sołowski, Berzins, Coombs, Guilkey, Möller, Tran, Adibaskoro, Seyedan, Tielen and Soga2021) and Nguyen, de Vaucorbeil & Bordas (Reference Nguyen, de Vaucorbeil and Bordas2023) for further details on MPM.
4. Results
In this section, various granular flow and collapse simulations using the critical state 
$\mu (I)$-rheology model are presented. In § 4.1, simulation results from flow on inclined planes are compared with approximate analytical solutions of the velocity profile. Next, comparisons with previous experiments of granular collapse problems are made in § 4.2, and cohesive media are considered in § 4.3. Finally, the ability of the model to capture the two steady-state solutions of granular flow over a smooth two-dimensional bump is demonstrated in § 4.4. In all simulations, the Young's modulus 
$E$, Poisson's ratio 
$\nu$ and hardening factor 
$\xi$ are kept constant with values shown in table 2. Furthermore, the initial compressive strength 
$p_c^0$ is taken to be 100 Pa and the parameter 
$\beta$ is set to zero, except in § 4.3 where cohesive simulations are presented. Appendix B presents further reasoning behind these constitutive parameter choices.
Elastic and plastic model parameters used throughout this study, except in Appendix B where a sensitivity analysis of all parameters is presented.
4.1. Analytical solutions
Following the analytical treatment of GDR MiDi (2004), a cohesionless two-dimensional flow of height 
$h$ as sketched in figure 5 is considered. Under the assumption 
$\sigma _{xx} = \sigma _{zz}$ justified in da Cruz et al. (Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005) and Andreotti et al. (Reference Andreotti, Forterre and Pouliquen2013), 
$q = \rho g (h-z) \sin \theta$ and 
$p = \rho g (h-z) \cos \theta$. Assuming all stress states lie on the CSL results in 
$\mu (I) = q/p = \tan \theta$. Combined with (2.26) and (2.27), using 
$\beta =0$ and assuming all deformation is plastic, one obtains
\begin{equation} \frac{{\rm d} v_x(z)}{{\rm d}z} = \dot{\gamma}_{_S} = \omega \frac{\tan \theta - \mu_1}{\mu_2 - \tan \theta} \sqrt{\rho g (h-z)\cos \theta} {.} \end{equation}
Integrating over 
$z$, assuming a no-slip bottom surface, yields
\begin{equation} v_x(z) = v_{max} ( 1 - (1-z/h)^{3/2} ), \end{equation}where
\begin{equation} v_{max} = v_x(h) = \frac{2}{3} \omega \frac{\tan \theta - \mu_1}{\mu_2 - \tan \theta} \sqrt{\rho g h^3 \cos \theta}\end{equation}is the velocity at the surface of the flow.
Two-dimensional flow on a no-slip surface showing a Bagnold velocity profile. While (a) is a sketch of the set-up and flow, (b) shows the time evolution of the flow profile at a fixed 
$x$-position 
$2.5$ m downstream for a simulation with slope angle 
$\theta = 24^\circ$. Here, the black dashed line is the analytical solution from (4.2) where no fitting to the measured data has been made.
To see how the proposed model compares with the approximate analytical solution given by (4.2) and (4.3), two-dimensional simulations are performed assuming cohesionless spherical glass bead parameters given by Jop et al. (Reference Jop, Forterre and Pouliquen2005) based on the experiments by Forterre & Pouliquen (Reference Forterre and Pouliquen2003), i.e. 
$I_0=0.279$, 
$\mu _1=\tan (20.9^\circ )$ and 
$\mu _2=\tan (32.8^\circ )$, with initial density 
$\rho ^0=1550\ \textrm {kg}\ \textrm {m}^{-3}$ and intrinsic density 
$\rho _*=2500\ \textrm {kg}\ \textrm {m}^{-3}$ as summarized in table 3. In these simulations, a mass of 280 kg (per unit width) is being slowly fed onto a no-slip inclined surface and through a gate restricting the maximum height of the flow. The bead diameter is 
$d=2.5$ mm here, and the grid size 
$\Delta x$ is taken to be 1.2 mm. As expected, for slope angles 
$\theta$ smaller than 
$\tan ^{-1}(\mu _1) = 20.9^\circ$, the flow eventually stops, while for 
$\theta$ between 
$\tan ^{-1}(\mu _1)$ and 
$\tan ^{-1}(\mu _2)$, i.e. between 
$20.9^\circ$ and 
$32.8^\circ$, a steady state forms after some time. For the particular case of 
$\theta = 24^\circ$ at a point 2.5 m downstream, the plot in figure 5 shows how the velocity profile converges towards the Bagnold velocity profile given by (4.2) after a few seconds. The surface velocity measured for all steady-state flow cases is shown in figure 6(a), showing the agreement with the prediction by (4.3). As the flow transitions to a steady state, figures 6(b) and 6(c) show that 
$\mu$ approaches the expected constant value. Since the analytical solution assumes all stress states lie on the CSL, a slight deviation from the analytical solution is expected. The minor observed deviation can also be attributed to the non-idealized set-up where the particles are fed through a gate without a fixed flow height, resulting in a not perfectly smooth flow and constant flow height, and the plotted quantities are determined through an averaging scheme over a predefined region.
Parameters used in § 4.1.
Simulations on varying slope angles 
$\theta$. In (a), the maximum velocity 
$v_{max}$ is measured along the surface of the flow and the analytic curve is (4.3) with 
$h = 10-0.3\theta$ cm fitted from the measured data. The error bars represent one standard deviation as the surface velocity was found by averaging over the spatial extent of the flow, neglecting only the area close to the gate and the area close to the flow front. In (b,c), the evolution of the average 
$\mu$, denoted 
$\langle \mu \rangle$, obtained from a slice of thickness 0.02 m in each simulation is shown. The latter is plotted against a scaled time, where 
$t_f$ denotes the time at which all mass has passed through the gate. The two dashed horizontal lines marks the lower and upper bounds 
$\mu _1$ and 
$\mu _2$, respectively.
In the steady-state flow simulations presented here, the solid fraction 
$\phi = \rho /\rho _*$ settles at an approximately temporally and spatially constant value for each slope angle. The hardening factor 
$\xi$ affects the compressibility of the system, here it is chosen large enough to only induce small changes in the solid fraction of the system. In particular, the solid fraction, initially being 0.62, varies between 0.60 and 0.63. The lowest fractions are measured for the highest slope angles where the flow is more rapid, qualitatively similar to the 
$\phi (I)$-relation suggested by GDR MiDi (2004), da Cruz et al. (Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005) and Pouliquen et al. (Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006). This is shown in figure 7, which was obtained by averaging the solid fraction in a slice in the middle of the flow. Initially, at 
$t/t_f<0.1$, the erratic behaviour of the measured average solid fraction is also attributed to the unstable nature of the first particles coming through the gate combined with the lack of particles needed for a stable average.
Evolution of the solid volume fraction 
$\phi$ with time as a steady-state flow is reached on different inclinations. As in figure 6, 
$t_f$ denotes the time at which all mass has passed through the gate.
Studying here the continuous (highly dynamic) flow where the deformations are dominantly plastic, the results presented are largely independent of the elastic moduli. Thus, whether 
$E$ is 1 MPa or 1 GPa, the results are indistinguishable. The same can be said about the exact onset of the plastic deformation as long as it occurs at sufficiently small stress values, in particular it does not matter in the presented results whether 
$p_c^0$ is 10 Pa or 100 Pa. This is demonstrated in Appendix B.
4.2. Granular collapse
Here, the problem of cohesionless granular collapse on horizontal and tilted surfaces is addressed. For the horizontal case, the experiments of Xu et al. (Reference Xu, Jin, Tai and Lu2017) are considered, where glass beads initially confined between two gates separated 8 cm apart were released by opening the gates such that the beads could flow between two parallel transparent glass sheets spaced 2 cm apart. The initial height of the system was 5 cm. Three-dimensional simulations are performed with the spherical glass bead parameters of § 4.1, using the known bead diameter 
$d = 2$ mm and intrinsic density 
$\rho _* = 2500\ \textrm {kg}\ \textrm {m}^{-3}$ reported from the experiments, as well as 
$\Delta x = 0.7$ mm (
$N_p = 1.7$M). While the initial bulk density 
$\rho ^0$ was measured to vary between 1450 and 
$1625\ \textrm {kg}\ \textrm {m}^{-3}$, in the simulations it is for simplicity assumed that 
$\rho ^0 = 1550\ \textrm {kg}\ \textrm {m}^{-3}$, close to the average experimental value. The parameters are summarized in table 4. The walls are assumed frictionless and the ground no-slip, inferring the gate lift velocity 
$0.45\ \textrm {m}\ \textrm {s}^{-1}$ from the experimental photos. Figure 8 shows the results alongside the reported experimental observations, indicating an overall better agreement between the model and the experiments at time 
$t=0.42$ s (close to the rest state) than at the intermediate time. At the final time, note that the experimental data indicate that some granules spread out rather far from the main part of the system, resulting in long tails which were not captured by the simulations. This is believed to be due to the no-slip boundary condition at the bottom surface which in the experiments was a steel plate. Nevertheless, neglecting the tails, the repose angle of the main part of the collapsed system was reproduced remarkably well. The elastoplastic nature of the model naturally leads to the simulated granular collapse exhibiting localized bands of increased plastic (shear) strain. Shear bands are observed in granular flows (e.g. Mueth et al. Reference Mueth, Debregeas, Karczmar, Eng, Nagel and Jaeger2000; Fenistein & van Hecke Reference Fenistein and van Hecke2003; Jop et al. Reference Jop, Forterre and Pouliquen2006; Mandal, Nicolas & Pouliquen Reference Mandal, Nicolas and Pouliquen2021) and will be further discussed at the end of this section as well as in Appendices B and C. Regarding computational time, reaching the last state shown in figure 8 took almost 27 h on 8 cores (Intel Xeon Gold 6136 CPU 3.00 GHz, 187 GB RAM, OpenMP parallelization). In general, the computational time is, in addition to the degree of discretization, dependent on the Young's modulus and density as these parameters control the adaptive time stepping of the explicit scheme.
Granular collapse on a horizontal plane. (a) A middle slice from the simulations is shown and compared with the experiments from Xu et al. (Reference Xu, Jin, Tai and Lu2017) at different times. Note that the plotted slice is arbitrary as the front remains more or less uniform along the width as can be seen in the corresponding three-dimensional visualizations in (b). Here, the material is coloured according to the horizontal speed.
Parameters used to simulate the experiment by Xu et al. (Reference Xu, Jin, Tai and Lu2017).
Having studied granular collapse on a horizontal surface, we now consider collapse on tilted surfaces. In the dam break experiments of Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010), cohesionless glass beads with diameter 
$d=0.7 \pm 0.1$ mm and intrinsic density 
$\rho _* = 2500\ \textrm {kg}\ \textrm {m}^{-3}$ were released on inclined rough surfaces. The system had an initial density 
$\rho ^0 = 1550\ \textrm {kg}\ \textrm {m}^{-3}$, initial height 
$h_0=14$ cm and initial length 20 cm along the inclined flow direction, enclosed by Plexiglas sidewalls 10 cm apart. These experiments are here simulated in three dimensions (using 
$\Delta x = 2.4$ mm, 
$N_p = 1.6$M) with a frictionless back wall preventing movement in the upflow direction, assuming a no-slip bottom surface as well as frictionless sidewalls. With the authors reporting the measured repose and avalanche angle as 
$23.5 \pm 0.5^\circ$ and 
$25.5 \pm 0.5^\circ$, respectively, 
$\mu _1$ is thus chosen here as the inverse tangent of the mean of these two angles. The remaining parameters, i.e. 
$I_0$ and 
$\mu _2$, are chosen as the typical spherical glass bead parameters as used for the collapse on the horizontal plane and in § 4.1. The parameters are summarized in table 5. With these assumptions, the final equilibrium state from the experiments is reproduced in the simulations. This can be seen in figure 9 for three different inclination angles 
$\theta = 10^\circ$, 
$16^\circ$ and 
$22^\circ$. The three-dimensional simulations required approximately 18 computational hours per real second on 8 cores on the same computer architecture as given above. While the agreement in figure 9 appears excellent, it should be emphasized that the simulations were conducted under the assumption of frictionless sidewalls. The sidewalls may influence the final deposit as demonstrated by Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017). With the rough bed surface being covered by the same beads glued to the surface, the no-slip assumption at the base is also questionable. In particular, Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) argue that there is slip near the front.
Parameters used to simulate the experiments by Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010). These parameters were also used in § 4.4 to simulate flow over a smooth obstacle.
Final rest state of the flow of a dam break on planes with three different inclinations. Simulations are compared with experiments from Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010).
Comparisons with the dam break experiments at intermediate times are shown in figure 10, which also highlights the static–fluid transition. This figure demonstrates that the simulations give very similar results compared with the 
$\mu (I)$-rheology model of Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017) without wall friction. Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017) also conducted simulations with wall friction, naturally resulting in a better match to the experimental data.
Dam break on two different inclinations, 
$\theta =10^\circ$ (a,c) and 
$\theta =22^\circ$ (b,d), at two intermediate times 
$t$ prior to reaching the final equilibrium position. The simulations are coloured in grey scale to highlight the static–flowing transition: dark grey is static and light grey is flowing. Experimental data from Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010) as well as previous simulations by Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017) with and without wall friction are included for comparison.
The three inclination angles in the dam break simulations are all below 
$\tan ^{-1}(\mu _1) = 24.5^\circ$, and thus the flow will come to rest, i.e. the front position converges to a fixed value. Figure 11(a) shows the overall capability of the model in capturing the evolution of the downslope front position 
$x_f$ in the dam break simulations. Comparing the front velocities, figure 11(b), highlights the not perfectly smooth front evolution of the simulated collapse on the largest inclination. On the two smaller inclinations, while still not achieving a perfectly smooth evolution, the front velocity remains very close to the experimental measurements. Figure 11(c) shows the evolution of the solid volume fraction 
$\phi$ with time, indicating a smaller final density with increasing inclination angle. After an initial increase in solid fraction (for 
$t<0.1$ s) due to consolidation under gravity, the solid fraction decreases gradually during the collapse. In the experiments, the system had already settled under gravity and the experimental data is therefore not expected to feature the initial jump in solid fraction observed in the simulations. Although the experimental data are scattered, especially for the largest inclination, the data seem to indicate an overall decrease in solid fraction during the collapse. This is qualitatively the same behaviour shown in the simulations after the initial consolidation.
(a) Front position 
$x_f$, (b) velocity 
$v_f$ and (c) solid fraction 
$\phi$ of the flow following the dam break on inclined planes. In (c), note that this is the solid volume fraction of the complete system and not a local measurement. Experimental data from Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010) and Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017).
A three-dimensional rendering of the 
$\theta =22^\circ$ dam break is shown in figure 12(a) visualizing the downslope velocity. While the final equilibrium state features a smooth surface, a non-smooth surface can be observed at intermediate times, even in the cohesionless case. This was not reported in the experiments. Visualizing instead the equivalent plastic shear strain rate in figure 13, it is immediately clear that the surface bumps are correlated with the occurrence of shear bands, i.e. thin zones of localized plastic shear strain. This is related to the initial state of the granular system and the consolidation that immediately follows due to gravity, thus producing states with large compressive strengths 
$p_c$ due to the imposed hardening law. In Appendix B, it is demonstrated that adjusting the plastic parameters, especially 
$\xi$, the gravity-induced initial compaction can be significantly reduced, and consequently it is possible to partly smooth out the surface bumps observed in these granular collapse simulations. Strain localization is expected in elastoplastic modelling featuring strain softening, and therefore the proposed model will fail in accurately capturing localization patterns without a rigorous regularization technique. Nevertheless, this is shown in Appendix C to not influence the general dynamics of the collapse and flow.
Dam break on slope angle 
$\theta =22^\circ$ of (a) cohesionless and (b) cohesive glass beads. The full time evolution of these two simulations is shown in supplementary movies 1 and 2 available at https://doi.org/10.1017/jfm.2024.643, respectively.
Cohesionless dam break on slope angle 
$\theta =22^\circ$ coloured according to the equivalent plastic shear strain rate 
$\dot{\gamma}_{_S}$.
4.3. The case of cohesion
While the above experiments and simulations were conducted with cohesionless glass beads, it is instructive to consider the case of cohesive beads. To this end, the dam break case on the slope angle 
$\theta =22^\circ$ is studied with 
$\beta >0$. In figure 12, snapshots of the dam break at different times are presented for the case 
$\beta = 0$ and 
$\beta = 0.3$, allowing us to visually observe the difference in the evolution of the flow when cohesion is introduced. In particular, the extension of the flow is reduced and lumps of material adhere to each other at the surface of the flow. The two simulations presented in figure 12 are also shown in supplementary movies 1 and 2, respectively, at 
$\frac {1}{3}$ speed. Figure 14 shows how the evolution of the front position decreases with increasing 
$\beta$.
Front position 
$x_f$ of flow on inclined plane 
$\theta =22^\circ$ with increasing cohesion 
$\beta = 0$, 0.1, 0.2, 0.3.
The recent publication of Gans et al. (Reference Gans, Abramian, Lagrée, Gong, Sauret, Pouliquen and Nicolas2023) presented detailed experiments of granular collapse with cohesive glass beads produced using a novel technique by Gans, Pouliquen & Nicolas (Reference Gans, Pouliquen and Nicolas2020) of coating the beads with a thin layer of polyborosiloxane. Since the coating provides a stable cohesion in the form of a constant-in-time inter-granular adhesive force, these experiments can be considered suitable for assessing the proposed model, which, in its current form, assumes a persistent (unbreakable) cohesion. The experiments were conducted in a 15.4 cm wide channel having a rough bed constructed by gluing identical grains as the flowing grains, the granular system having an equal initial height 
$h_0$ and length of 8.9 cm. Two of their experiments, one with slightly more cohesive beads than the other, are presented in figure 15, where we call the least cohesive case Exp. A and the most cohesive case Exp. B. In the case of Exp. A, the surface of the system is plotted at two different times. The figure also includes results from two-dimensional simulations where the spherical glass bead parameters as presented in § 4.1 were utilized besides from those reported from the experimental measurements, in particular 
$d=0.8$ mm, 
$\rho _* = 2600\ \textrm {kg}\ \textrm {m}^{-3}$ and 
$\rho ^0 = 1508\ \textrm {kg}\ \textrm {m}^{-3}$. These parameters are summarized in table 6. As in the previous section, a no-slip ground was assumed. In the presented simulations the cohesion through 
$\beta p_c^0$ is tuned to match the experiments, here using 
$p_c^0=500$ Pa. As before, it is clear that runout is decreased with increasing cohesion. In particular, using 
$\beta =0.2$ and 
$\beta =0.3$ can be seen to give a front position evolution that aligns with the experiments, albeit not perfectly. As in the previous section, it should be mentioned that the assumption of a no-slip base and frictionless walls may not be fully accurate. Moreover, the previous comments made about strain localization in the cohesionless case also applies here, and the reader is again referred to Appendices B and C. More importantly, since the cohesion in the proposed model is tuned to match the experiments, it is pertinent to emphasize that the presented simulations do not constitute a rigorous validation of the model.
Cohesive granular collapse simulations compared with the experiments from Gans et al. (Reference Gans, Abramian, Lagrée, Gong, Sauret, Pouliquen and Nicolas2023) on a horizontal rough surface. (a,b) Snapshots of the simulation of Exp. A using 
$\beta =0.2$ at two different times where the markers indicate the surface measured in the experiments. (c) Evolution of the front position 
$x_f$ with time in both experiments are shown together with simulations with 
$\beta =0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3$.
Parameters used to simulate the experiments by Gans et al. (Reference Gans, Abramian, Lagrée, Gong, Sauret, Pouliquen and Nicolas2023) in figure 15.
Having investigated the effect of cohesion on granular collapse, one can also consider more prolonged flows. From discrete simulations, it has been demonstrated that cohesive granular flows can reach a steady state where the introduction of cohesive forces between the grains leads to the occurrence of a plug flow near the free surface (Rognon et al. Reference Rognon, Roux, Naaïm and Chevoir2008b; Mandal et al. Reference Mandal, Nicolas and Pouliquen2020). It appears possible for the proposed critical state 
$\mu (I)$-rheology model to qualitatively capture this behaviour, as presented in figure 16. Here, on an inclined slope of 
$\theta = 32^\circ$ with a no-slip base, steady-state flows are reached for various values of 
$\beta$. With increasing 
$\beta$ the occurrence of a plug flow is observed, and above a critical 
$\beta$ the flow stops. In these simulations, the same set-up and parameters as in § 4.1 were used (including 
$p_c^0 = 100$ Pa), except with a twice as large gate height and 
$\mu _2 = \tan (40^\circ ) > \tan (\theta )$ to ensure steady-state flow. These parameters are listed in table 7. A comparison with the results of the discrete simulations of Mandal et al. (Reference Mandal, Nicolas and Pouliquen2020) reveals clear differences, both in terms of the shape of the velocity profile under the plug and in terms of the relative change in velocity magnitude as the cohesion is increased. While the qualitative results from the cohesive simulations presented in this section are promising and in line with expectations, a robust validation of the cohesive model and its implementation remains to be addressed.
Steady-state velocity profiles for increasing cohesion parameter 
$\beta$ on a no-slip surface tilted at 
$\theta = 32^\circ$. The velocity profiles are measured at the same position 
$x=0.8$ m from the gate which had a height of 
$40d$.
Parameters used to simulate cohesive steady-state flow (figure 16).
4.4. Multiple solutions for flow over a smooth bump
In small-scale experiments of granular flow over a smooth obstacle, Viroulet et al. (Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017) showed that two very different steady-state regimes can occur. In particular, either a ‘jet’ flow detaching from the apex of the bump or a ‘shock’ forming upstream of the bump is observed, the latter being the case if a sufficient amount of erodible particles are placed in front of the bump. This initially static mass reduces the energy of the impacting flow and accordingly also the flow velocity, generating an increased flow thickness before the bump. Curious as to whether the proposed critical state 
$\mu (I)$-rheology model can capture both steady-state solutions, simulations are performed assuming beads with the same diameter and on the same terrain/bump as used in the experiments. Along the flow direction 
$x$, the elevation of the smooth symmetric bump is given by (in metres)
\begin{equation} b(x) = 0.0475 \,{\rm sech}\left(\frac{x-0.43}{0.04} \right), \end{equation}
as sketched in figure 17. At 
$x=0$, i.e. 43 cm from the centre of the bump, the beads flow through a gate of height 1.5 cm. With the bump being two-dimensional and the flow being confined in a narrow chute, the problem is here simulated in two spatial dimensions. Cohesionless spherical glass beads with identical parameters as in the dam break simulations of § 4.2 are assumed, i.e. the parameters presented in table 5. In the experiments of Viroulet et al., the beads show slippage along the base, prompting the introduction of a base friction in the simulation, for simplicity as a Coulomb friction law approximated on the velocities. In particular, introducing a base friction parameter 
$\mu _b$, the velocity 
$\boldsymbol {v}$ relative to the boundary is obtained as
\begin{equation} \boldsymbol{v} = \begin{cases} \boldsymbol{v}^*_{T} - \mu_b ||\boldsymbol{v}^*_{N}||_{L_2} \dfrac{\boldsymbol{v}^*_{T}}{||\boldsymbol{v}^*_{T}||_{L_2}} & \textrm{if}\; ||\boldsymbol{v}^*_{T}||_{L_2} > \mu_b ||\boldsymbol{v}^*_{N}||_{L_2} {,} \\ \boldsymbol{0} & \text{otherwise}, \end{cases} \end{equation}Sketch of terrain in experiments of Viroulet et al. (Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017).
where 
$\boldsymbol {v}^*_{T}$ and 
$\boldsymbol {v}^*_{N}$ refer to the tangential and normal velocity components, respectively, of 
$\boldsymbol {v}^*$, i.e. the relative velocity before the boundary condition is considered. Equation (4.5) is only used if 
$\boldsymbol {v}^*_N$ suggests penetration of the boundary, otherwise 
$\boldsymbol {v} = \boldsymbol {v}^*$. Such a boundary condition (4.5), which has been extensively used in MPM flow modelling before (Gaume et al. Reference Gaume, Gast, Teran, van Herwijnen and Jiang2018; Li et al. Reference Li, Sovilla, Jiang and Gaume2020; Cicoira et al. Reference Cicoira, Blatny, Li, Trottet and Gaume2022) and which is further validated in Appendix E, is motivated from having to apply boundary conditions on the grid where the velocities are readily available, unlike the stresses. It is reported from the experiments that, with no initial mass in front of the bump, the beads detach from the apex of the bump for slope angles around approximately 
$35^\circ$, which justifies roughly 
$\mu _b=0.5$. As in the experiments, the problem is here simulated with and without an initial mass (here 0.13 kg) in front of the bump. Figures 18 and 19 show that the simulations at a slope angle 
$\theta =39^\circ$ capture both of the two different steady-state regimes. In the case with initial mass in front of the bump, figure 18 shows the ‘shock’ forming upstream, demonstrating an overall good agreement with the experimental observations. The complete time evolution in this simulation can be seen in supplementary movie 3. Without initial mass in front of the bump, figure 19 shows the detaching ‘jet’ steady-state flow. In this case, the experimental images indicate a significant vertical spreading of the particles, a situation which is not observed in the simulations. This is believed to be due to the sidewalls in the experiments, where the wall friction reduces the speed of the particles closer to the walls. Consequently, these particles will experience a shorter trajectory over the bump compared with the particles in the middle of the flow. As the experimental images are taken from the side through the transparent walls, the different trajectories are not visually distinguishable. Three-dimensional simulations with sidewall friction would allow for a more proper comparison with experiments, and a more extensive parameter tuning could improve the agreement with experimental trajectories. Using a grid size of 
$\Delta x = 0.6$ mm and 
$N_p = 0.65$M, the two-dimensional simulations presented here required approximately 8 computational hours per second on 8 cores on the same computer architecture as presented in § 4.2. Supplementary movie 4 shows the complete time evolution in the jet flow simulation. In this movie, it can be noted how the flow slightly pulses in response to small waves that are propagating down from the inflow. This is here due to the not perfectly smooth feeding of particles through the gate, but could also be a physical result of the flow breaking down into roll waves on sufficiently long chutes (Barker & Gray Reference Barker and Gray2017; Viroulet et al. Reference Viroulet, Baker, Rocha, Johnson, Kokelaar and Gray2018). Nevertheless, it is not observed in the experiments.
Flow over bump with initial mass and 
$\theta =39^\circ$. Experiments of Viroulet et al. (Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017) (left) and simulations (right) at increasing time from top to bottom (
$t = 0, 0.4, 0.7, 1, 1.5, 4$ s). The full time evolution is shown in supplementary movie 3. Note the colour bar only applies to the simulations, not the experiments.
Flow over bump without initial mass and 
$\theta =39^\circ$. Experiments of Viroulet et al. (Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017) (left) and simulations (right) at increasing time from top to bottom (
$t = 0.3, 0.6, 0.9, 4$ s). The full time evolution is shown in supplementary movie 4. Note the colour bar only applies to the simulations, not the experiments.
5. Discussion and conclusions
This paper has presented a continuum model based on the combination of the 
$\mu (I)$-rheology with the critical state concept from geo-mechanics in a way that captures both solid and fluid properties of granular media in a single framework. Incorporating the elastoplastic CSSM theory, the proposed model is different from the previously proposed critical state 
$\mu (I)$-rheology models for dry granular flows (Barker et al. Reference Barker, Schaeffer, Shearer and Gray2017; Schaeffer et al. Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019; Bouchut et al. Reference Bouchut, Fernández-Nieto, Koné, Mangeney and Narbona-Reina2021). While these previous models are either implemented in depth-averaged schemes (in one or multiple layers) or realized only for one-dimensional problems, the elasto-viscoplastic model of this paper is implemented in a B-spline MPM framework allowing full three-dimensional simulations. Accounting for cohesion and compressibility through CSSM, the proposed model is more general than the two-dimensional MPM-implemented 
$\mu (I)$-rheology of Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015). The proposed model approximately reproduces the expected flow profiles on inclined slopes, and with the introduction of a smooth bump along the slope, the model is able to qualitatively provide the two steady-state regimes observed in the experiments of Viroulet et al. (Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017). While for cohesionless glass beads a Drucker–Prager yield may be more appropriate, for other granular materials, in particular snow, closed yield surfaces (such as the MCC yield surface) are more appropriate (Reiweger et al. Reference Reiweger, Gaume and Schweizer2015; Ritter et al. Reference Ritter, Gaume and Löwe2020; Blatny et al. Reference Blatny, Löwe, Wang and Gaume2021). Snow also undergoes significant volumetric deformations under typical loading conditions compared with cohesionless glass beads. The MCC model already facilitates an uncomplicated treatment of such compaction and dilation.
The simulations of granular collapse of spherical glass beads displayed shear bands affecting the smoothness of the surface at intermediate times before reaching the final equilibrium position. This was attributed to the initial state of the granular system and the immediate consolidation occurring as a result of the onset of gravity, thus leading to a sudden increase in 
$p_c$. As the gravity-induced compaction can be suppressed through adjustments of the plastic parameters of the hardening law, Appendix B demonstrated that the surface during the collapse can be smoothed. While such shear band-induced surface patterns were not reported in the spherical glass bead collapse experiments, they are generally observed in experiments involving dry granular media (Gudehus & Nübel Reference Gudehus and Nübel2004; Torres-Serra, Romero & Rodríguez-Ferran Reference Torres-Serra, Romero and Rodríguez-Ferran2020). Nevertheless, neglecting the surface effects, §§ 4.2 and 4.3 showed that the proposed CSSM-based 
$\mu (I)$-rheology model is able to capture the observed dynamics and final deposit in granular collapse experiments of spherical glass beads on both flat and inclined surfaces reasonably well.
While Lacaze & Kerswell (Reference Lacaze and Kerswell2009) showed, using discrete simulations, that the 
$\mu (I)$-rheology holds across the whole flow of such problems, others have presented numerical experiments that seem to indicate that the 
$\mu (I)$-rheology is not needed to capture the specific problem of granular collapse. Specifically, the recent paper of Rousseau et al. (Reference Rousseau, Métivet, Rousseau, Daviet and Bertails-Descoubes2022) claims that only a single friction parameter is needed to describe granular collapse and further argues that the 
$\mu (I)$-rheology is an over-complication for modelling such problems. This is analogous to the previous claims of Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) and Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017) who, implementing the 
$\mu (I)$-rheology in a Drucker–Prager viscoplastic approach with variable viscosity, showed that the results did not differ significantly if the viscosity was assumed constant. Similarly, Crosta, Imposimato & Roddeman (Reference Crosta, Imposimato and Roddeman2009) showed that granular collapse can be accurately modelled with a constant friction angle in a Mohr–Coulomb yield rule. Along these lines, it seems reasonable to assume that, with appropriately chosen constitutive parameters, a rate-independent CSSM (i.e. with a constant, 
$I$-independent, slope of the CSL) may also be adequate in modelling the granular collapse problem. This is shown and discussed in Appendix D.
Moreover, it should be mentioned that the wall friction may influence the results. In the simulations of granular collapse presented here (on various inclinations), all sidewalls were assumed frictionless. While investigating the influence of wall friction in granular flow was not within the scope of this study, such simulations could be accomplished by employing the frictional boundary conditions outlined in § 4.4. The effect of sidewalls during granular collapse has been analysed, e.g. numerically by Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017), showing that introducing friction may help in better capturing the experimental results. Along the same lines, slippage of the material along the base may also influence the results. The base was assumed no-slip in all granular collapse simulations presented here, and it was highlighted as a potential reason for the short runout compared with the experimental observations by Xu et al. (Reference Xu, Jin, Tai and Lu2017) of granular collapse on the horizontal surface. In particular, Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) reported a 10 % increase in runout going from a no-slip base to a basal friction of 0.48. The reader is referred to these references for further insight into the influence of wall and base friction on granular collapse. The flow over bump problem presented in § 4.4 was the only set-up in this work where a non-zero finite base friction was employed. Following the MPM community, this was accomplished by applying boundary conditions on the grid velocities, proving an approximation to the Coulomb friction law which was further validated in Appendix E. While one could foresee imposing boundary conditions on the stresses, this would come at an increasing computational cost as the stresses are not readily available on the grid.
The cohesion assumed in the proposed model is a form of cohesion that is always present (for 
$\beta$, 
$p_c > 0$). In order to apply the model to flows where this may not be an appropriate description of cohesion, the model could be altered to remove cohesion, i.e. setting 
$\beta$ to zero, after yield or another pre-defined criterion. With increasing cohesion, it was found that steady-state flow can still be observed, with increasing 
$\beta$ leading to plug flow near the surface, qualitatively similar to the observations in the discrete simulations of Rognon et al. (Reference Rognon, Roux, Naaïm and Chevoir2008b) and Mandal et al. (Reference Mandal, Nicolas and Pouliquen2020) of cohesive granular flows. An exact quantitative comparison between the proposed model and previous discrete simulations remains challenging as this entails relating 
$\beta$ to the prescribed inter-particle force parameter used in the discrete simulations. In this study, the cohesion parameter 
$\beta$ imposes a translation of the CSL in principal stress space. This is analogous to the usual treatment of cohesion in viscoplastic fluids, such as Bingham fluids, where cohesion can be accounted for by altering the yield stress. For example, Abramian, Staron & Lagrée (Reference Abramian, Staron and Lagrée2020) and Gans et al. (Reference Gans, Abramian, Lagrée, Gong, Sauret, Pouliquen and Nicolas2023) considered cohesion similarly by altering a simple 
$I$-dependent Mohr–Coulomb yield criterion in their volume-of-fluid incompressible Navier–Stokes scheme. Comparing with discrete simulations, Abramian et al. (Reference Abramian, Staron and Lagrée2020) concluded that such approach can accurately capture the effective cohesion induced by cohesive discrete particles. In fact, they proposed explicit relations between their (macroscopic) Mohr–Coulomb cohesion parameter and the (microscopic) inter-particle cohesive force parameter. Such a study could also be performed in the context of the model proposed here, in order to establish a link between 
$\beta$ (and 
$p_c^0$) and the inter-particle attractive force. On the contrary, the numerical studies by Mandal et al. (Reference Mandal, Nicolas and Pouliquen2020, Reference Mandal, Nicolas and Pouliquen2021) suggest that, while the inter-particle attractive force controls the initiation of the flow, the cohesive flow dynamics is better described by an effective attractive force that also takes into account the stiffness and inelasticity of the grains.
In the critical state 
$\mu (I)$-rheology model presented here, the granular material can undergo both elastic and plastic compaction and dilation. Unless the initial compressive strength is chosen to be very large, the elastic volumetric change is negligible and all volumetric deformation can be considered plastic. In § 2.4, this assumption was used to derive a hardening law that recovers the logarithmic relationship between the inverse solid volume fraction and the isotropic compressive strength, which is a relationship commonly used for the (slow) deformation of soils and powders (Walker Reference Walker1923; Poquillon et al. Reference Poquillon, Lemaitre, Baco-Carles, Tailhades and Lacaze2002; Castellanos Reference Castellanos2005). However, as mentioned in the Introduction, for granular media undergoing flow, a 
$\phi (I)$-relation has been reported, relating the solid volume fraction to the inertial number. In particular, above a critical inertial number, the solid fraction typically decreases linearly with this number. A hardening law that recovers the 
$\phi (I)$-relation may be possible, although this was not explored in this work, where we opted for a 
$\phi (p)$-relation. However, as briefly mentioned in §§ 4.1 and 4.2, this did induce a decreasing solid volume fraction with flow speed, qualitatively similar to the 
$\phi (I)$-relation. Formulating an accurate unified hardening law for both solid and fluid granular media, i.e. a law reproducing 
$\phi (I)$ in the high inertial number limit and 
$\phi (p)$ in the low inertial number limit, would require extensive further work.
The critical state 
$\mu (I)$-rheology assumed in this paper is a local theory in the sense it does not explicitly account for the observed situation where stress or velocity fluctuations caused by rearrangement at one location trigger rearrangements elsewhere in the system (Gaume, Chambon & Naaim Reference Gaume, Chambon and Naaim2011, Reference Gaume, Chambon and Naaim2020). This is what led to the non-local theories of Pouliquen & Forterre (Reference Pouliquen and Forterre2009) and Kamrin & Koval (Reference Kamrin and Koval2012). In the future, the model presented in this paper could be extended to introduce non-local effects in order to reproduce, e.g. the previous experimental observations of Pouliquen (Reference Pouliquen1999), Pouliquen & Forterre (Reference Pouliquen and Forterre2002) and GDR MiDi (2004) indicating that there is a minimum thickness required to observe a steady-state flow at a given inclination angle. This minimal thickness decreases with inclination angle. In the local 
$\mu (I)$-rheology this cannot be reproduced without possibly considering a threshold on the flow speed. Very recently, Dunatunga & Kamrin (Reference Dunatunga and Kamrin2022) and Haeri & Skonieczny (Reference Haeri and Skonieczny2022) included non-local effects in their elasto-viscoplastic Drucker–Prager 
$\mu (I)$-rheology model.
The flow and collapse simulations highlighted in this paper provide a small glimpse at what is feasible with the presented theoretical and computational framework. In the future, it may be applied to various granular dynamics problems. For example, the numerical implementation in MPM makes this framework particularly suited to study erosion (Li et al. Reference Li, Sovilla, Ligneau, Jiang and Gaume2022, Reference Li, Sovilla, Gray and Gaume2024; Vicari et al. Reference Vicari, Tran, Nordal and Thakur2022). While this may be accomplished well with other methods (Lusso et al. Reference Lusso, Ern, Bouchut, Mangeney, Farin and Roche2017; Fernández-Nieto et al. Reference Fernández-Nieto, Garres-Díaz, Mangeney and Narbona-Reina2018), in MPM the terrain can simply be filled with the same or different material as that of the incoming flow, with no extra coupling or changes to Algorithm 1 needed. Other phenomena that could be studied include the formation of (free-surface) roll waves which may develop in granular flow on rough inclined planes (Forterre & Pouliquen Reference Forterre and Pouliquen2003; Barker & Gray Reference Barker and Gray2017; Viroulet et al. Reference Viroulet, Baker, Rocha, Johnson, Kokelaar and Gray2018). Implementing the model in other hybrid/mesh-free numerical solvers such as SPH or PFEM should also allow such studies. In addition, future studies could investigate the capability of the model, or extensions of the model, to capture levee formation, self-channelization (Félix & Thomas Reference Félix and Thomas2004; Mangeney et al. Reference Mangeney, Bouchut, Thomas, Vilotte and Bristeau2007; Rocha, Johnson & Gray Reference Rocha, Johnson and Gray2019; Edwards et al. Reference Edwards, Rocha, Kokelaar, Johnson and Gray2023) and finger formation (Woodhouse et al. Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012; Baker, Johnson & Gray Reference Baker, Johnson and Gray2016) observed, e.g. in debris flows, snow avalanches and pyroclastic flows (Jessop et al. Reference Jessop, Kelfoun, Labazuy, Mangeney, Roche, Tillier, Trouillet and Thibault2012). One could also foresee coupling the model with particle size segregation in order to capture the segregation phenomena that occur in poly-disperse granular assemblies (Rognon et al. Reference Rognon, Roux, Naaïm and Chevoir2007; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021; Maguire et al. Reference Maguire, Barker, Rauter, Johnson and Gray2024). While an increased interest in understanding the effect of cohesion in granular flows has been witnessed (Rognon et al. Reference Rognon, Roux, Wolf, Naaïm and Chevoir2006, Reference Rognon, Roux, Naaïm and Chevoir2008b; Berger et al. Reference Berger, Azéma, Douce and Radjai2015; Khamseh et al. Reference Khamseh, Roux and Chevoir2015; Badetti et al. Reference Badetti, Fall, Hautemayou, Chevoir, Aimedieu, Rodts and Roux2018; Abramian et al. Reference Abramian, Staron and Lagrée2020; Vo et al. Reference Vo, Nezamabadi, Mutabaruka, Delenne and Radjai2020; Mandal et al. Reference Mandal, Nicolas and Pouliquen2020, Reference Mandal, Nicolas and Pouliquen2021; Macaulay & Rognon Reference Macaulay and Rognon2021; Dong et al. Reference Dong, Wang, Marks, Chen and Gan2023; Gans et al. Reference Gans, Abramian, Lagrée, Gong, Sauret, Pouliquen and Nicolas2023), this topic has also recently become increasingly relevant in the snow and avalanche community (Rognon et al. Reference Rognon, Chevoir, Bellot, Ousset, Naaïm and Coussot2008a; Bartelt et al. Reference Bartelt, Valero, Feistl, Christen, Bühler and Buser2015; Li et al. Reference Li, Sovilla, Jiang and Gaume2020). This is particularly motivated by the notion that climate change induces an increase in the proportion of wet cohesive snow avalanches compared with dry ones (Ballesteros-Cánovas et al. Reference Ballesteros-Cánovas, Trappmann, Madrigal-González, Eckert and Stoffel2018; Součková et al. Reference Součková, Juras, Dytrt, Moravec, Blöcher and Hanel2022). Hence, we expect the proposed approach to have important impacts in this field, in particular for the evaluation of avalanche impact pressures, runout distances and more generally risk management and mitigation (Ortner et al. Reference Ortner, Bründl, Kropf, Röösli, Bühler and Bresch2023).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.643.
Acknowledgements
L.B. wishes to thank Dr H. Rousseau and B. Trottet for stimulating discussions, all members of the UCLA MultiPLES Lab who advised on numerical implementations and Professor J.-F. Molinari for advice on the flow-over-bump problem and for hosting the first author in his lab.
Funding
This work was supported by the Swiss National Science Foundation (grant number PCEFP2_181227). J.M.N.T.G. was supported by a Royal Society Wolfson Research Merit Award (WM150058) and an EPSRC Established Career Fellowship (EP/M022447/1). Additional support was provided by NERC grants NE/X00029X/1 and NE/X013936/1.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
All research data supporting this publication are directly available within this publication.
Appendix A. Operator splitting: elastic predictor–plastic corrector
This appendix outlines the elastic predictor–plastic corrector scheme following Simo (Reference Simo1992).
The symmetric positive–definite second-order right and left Cauchy–Green strain tensors are given by
\begin{equation} \boldsymbol{C} = \boldsymbol{F}^T \boldsymbol{F} = \sum_{\alpha=1}^{\dim} \lambda_\alpha^2 \boldsymbol{N}_\alpha \otimes \boldsymbol{N}_\alpha\end{equation}and
\begin{equation} \boldsymbol{b} = \boldsymbol{F}\boldsymbol{F}^T = \sum_{\alpha=1}^{\dim} \lambda_\alpha^2 \boldsymbol{n}_\alpha \otimes \boldsymbol{n}_\alpha {,} \end{equation}
respectively. Here, the unit vectors 
$\boldsymbol {n}_\alpha$ and 
$\boldsymbol {N}_\alpha$ define the principal spatial and material directions, respectively. Furthermore, the plastic right Cauchy–Green tensor is defined as 
$\boldsymbol {C}^P = (\boldsymbol {F}^P)^T \boldsymbol {F}^P$, and the elastic left Cauchy–Green tensor as 
$\boldsymbol {b}^E = \boldsymbol {F}^E (\boldsymbol {F}^E)^T$. With these definitions, the elastic left Cauchy–Green strain can be written as
\begin{equation} \boldsymbol{b}^E = \boldsymbol{F}^E (\boldsymbol{F}^E)^T = \boldsymbol{F} (\boldsymbol{C}^P)^{{-}1} \boldsymbol{F}^T {.} \end{equation}
Thus, the change in 
$\boldsymbol {b}^E$ over time can be decomposed into two terms
\begin{equation} \frac{{\rm D}\boldsymbol{b}^E}{{\rm D}t} = \left.\frac{{\rm D} \boldsymbol{b}^E}{{\rm D}t} \right|_{\boldsymbol{C}^P = \textrm{const.}} + \left.\frac{{\rm D} \boldsymbol{b}^E}{{\rm D}t} \right|_{\boldsymbol{F} = \textrm{const.}}, \end{equation}
where the last term quantifies the change in 
$\boldsymbol {b}^E$ independent of the total deformation. As shown in Bonet & Wood (Reference Bonet and Wood2008), these terms can be related to the velocity gradient, in particular
\begin{equation} \boldsymbol{L}^E = \frac{1}{2} \frac{{\rm D} \boldsymbol{b}^E}{{\rm D}t} (\boldsymbol{b}^E)^{{-}1} {,}\quad \boldsymbol{L}^P ={-}\frac{1}{2} \left.\frac{{\rm D} \boldsymbol{b}^E}{{\rm D}t} \right|_{\boldsymbol{F} = \text{const.}} (\boldsymbol{b}^E)^{{-}1}, \end{equation}and as such we have
\begin{equation} \left.\frac{{\rm D} \boldsymbol{b}^E}{{\rm D}t} \right|_{\boldsymbol{F} = \textrm{const.}} ={-}2 \boldsymbol{L}^P\boldsymbol{b}^E \stackrel{(2.15)}={-} 2 \dot{\gamma} \boldsymbol{R} \boldsymbol{b}^E ,\end{equation}
where 
$\boldsymbol {R} = {\partial y / \partial \boldsymbol {\sigma } }/{||\partial y / \partial \boldsymbol {\sigma }||}$ was introduced.
We introduce discrete time 
$t^n$, 
$n \in \mathbb {N}$ and consider the integration of (A4) from 
$\boldsymbol {b}^{E,n}$ at time 
$t^n$ to 
$\boldsymbol {b}^{E,n+1}$ at time 
$t^{n+1}$ in a step of 
$\Delta t$ through an operator splitting scheme (as presented in § 3.2 in Simo Reference Simo1992). First, neglect the last term on the right-hand side of (A4), resulting in a trial state 
$\boldsymbol {b}^{E, {trial}}$ due to elasticity only. Then, neglect the first term on the right-hand side, integrating only the second term, i.e. (A6), from 
$\boldsymbol {b}^{E, {trial}}$ to 
$\boldsymbol {b}^{E, n+1}$ with an exponential integrator
\begin{equation} \boldsymbol{b}^{E,n+1} = \exp({ - 2\dot{\gamma} \boldsymbol{R} \Delta t }) \boldsymbol{b}^{E, {trial}} ,\end{equation}
which is (3.12b) in Simo (Reference Simo1992). Noting that the elastic Hencky strain 
$\boldsymbol {\varepsilon }^E$ is directly related to 
$\boldsymbol {b}^E$ through 
$\boldsymbol {\varepsilon }^E = \frac {1}{2} \ln \boldsymbol {b}^E$, we can write the above in terms of the Hencky strain, in particular
\begin{equation} \boldsymbol{\varepsilon}^{E, n+1} = \boldsymbol{\varepsilon}^{E, {trial}} - \dot{\gamma} \boldsymbol{R} \Delta t {.} \end{equation}Splitting the Hencky strain into a symmetric and deviatoric part, this equation can be expanded to
\begin{align} \left[ \operatorname{dev}{\boldsymbol{\varepsilon}^{E, {trial}}} + \frac{{\rm tr}{(\boldsymbol{\varepsilon}^{E, {trial}})}}{\dim} \boldsymbol{I} \right] - \left[ \operatorname{dev}{\boldsymbol{\varepsilon}^{E, n+1}} + \frac{{\rm tr}{(\boldsymbol{\varepsilon}^{E, n+1})}}{\dim} \boldsymbol{I} \right] - \Delta t \dot{\gamma} \frac{ \dfrac{\partial y}{\partial p} \dfrac{\partial p}{\partial \boldsymbol{\sigma}} + \dfrac{\partial y}{\partial q} \dfrac{\partial q}{\partial \boldsymbol{\sigma}} }{ ||\partial y / \partial \boldsymbol{\sigma}|| } = 0 ,\end{align}
where for 
$\boldsymbol {R}$ it was assumed that 
$y = y(p,q)$ and using the chain rule. From the definition of 
$p$ and 
$q$, it follows that
\begin{equation} \frac{\partial p}{\partial \boldsymbol{\sigma}} ={-}\frac{1}{\dim} \boldsymbol{I}\end{equation}and
\begin{equation} \frac{\partial q}{\partial \boldsymbol{\sigma}} = \frac{1}{\sqrt{2}} \frac{\boldsymbol{\tau}}{||\boldsymbol{\tau} ||} = \frac{1}{\sqrt{2}} \frac{\operatorname{dev} \boldsymbol{\varepsilon}^E}{|| \operatorname{dev} \boldsymbol{\varepsilon}^E ||} = \frac{1}{\sqrt{2}} \frac{\operatorname{dev} \boldsymbol{\varepsilon}^{E, {trial}}}{|| \operatorname{dev} \boldsymbol{\varepsilon}^{E, {trial}} ||}, \end{equation}
where the second last equality follows from Hencky's elasticity model where 
$\boldsymbol {\tau } = 2G \operatorname {dev} (\boldsymbol {\varepsilon }^E)$, and the last equality can be seen by applying the deviatoric operator on (A8). Now, (A9) can be written as
\begin{align} &\operatorname{dev}{\boldsymbol{\varepsilon}^{E, \text{trial}}} - \operatorname{dev}{\boldsymbol{\varepsilon}^{E, n+1}} + \frac{1}{\dim}{\rm tr}{(\boldsymbol{\varepsilon}^{E, {trial}})} \boldsymbol{I} - \frac{1}{\dim}{\rm tr}{(\boldsymbol{\varepsilon}^{E, n+1})} \boldsymbol{I}\nonumber\\ &\quad - \frac{\Delta t \dot{\gamma}}{\sqrt{2}} \frac{\partial y/\partial q}{||\partial y / \partial \boldsymbol{\sigma}||} \frac{\operatorname{dev} \boldsymbol{\varepsilon}^{E, {trial}}}{|| \operatorname{dev} \boldsymbol{\varepsilon}^{E, {trial}} ||} + \frac{\Delta t \dot{\gamma}}{\dim} \frac{\partial y/\partial p}{||\partial y / \partial \boldsymbol{\sigma}||} \boldsymbol{I} = 0. \end{align}This can be split into two as
\begin{equation} \left.\begin{gathered} \left[{\rm tr}{(\boldsymbol{\varepsilon}^{E, {trial}})} - {\rm tr}{(\boldsymbol{\varepsilon}^{E, n+1})} + \Delta t \dot{\gamma} \frac{\partial y / \partial p}{||\partial y / \partial \boldsymbol{\sigma}||} \right]\boldsymbol{I} = 0 \\ \left[ ||\operatorname{dev}{\boldsymbol{\varepsilon}^{E, {trial}}}|| - ||\operatorname{dev}{\boldsymbol{\varepsilon}^{E, n+1}}|| - \frac{1}{\sqrt{2}} \Delta t \dot{\gamma} \frac{\partial y /\partial q}{||\partial y / \partial \boldsymbol{\sigma}||} \right] \frac{\operatorname{dev} \boldsymbol{\varepsilon}^{E, {trial}}}{|| \operatorname{dev} \boldsymbol{\varepsilon}^{E, {trial}} ||} = 0 \end{gathered}\right\}, \end{equation}which results in only two scalar equations
\begin{equation} \left.\begin{gathered} {\rm tr}{(\boldsymbol{\varepsilon}^{E, {trial}})} - {\rm tr}{(\boldsymbol{\varepsilon}^{E, n+1})} + \Delta t \dot{\gamma} \frac{\partial y / \partial p}{||\partial y / \partial \boldsymbol{\sigma}||} = 0\\ ||\operatorname{dev}{\boldsymbol{\varepsilon}^{E, {trial}}}|| - ||\operatorname{dev}{\boldsymbol{\varepsilon}^{E, n+1}}|| - \frac{1}{\sqrt{2}} \Delta t \dot{\gamma} \frac{\partial y /\partial q}{||\partial y / \partial \boldsymbol{\sigma}||} = 0 \end{gathered}\right\}, \end{equation}
which can be expressed in terms of 
$p = -K\,\textrm {tr}{(\boldsymbol {\varepsilon }^E)}$ and 
$q = \sqrt {2}G||\operatorname {dev}{(\boldsymbol {\varepsilon }^E)} ||$ as
\begin{equation} \left.\begin{gathered} p^{n+1} = p^{trial} - K \Delta t \dot{\gamma} \frac{\partial y /\partial p}{||\partial y / \partial \boldsymbol{\sigma}||} \\ q^{n+1} = q^{trial} - G \Delta t \dot{\gamma} \frac{\partial y /\partial q}{||\partial y / \partial \boldsymbol{\sigma}||} \end{gathered}\right\}, \end{equation}
where 
$K = \lambda + 2G/\dim$ is the bulk modulus. Alternatively, using the definitions in (2.17) and (2.18), this can be expressed as
\begin{equation} \left.\begin{gathered} p^{n+1} = p^{trial} + K \dot{\gamma}_{_V} \Delta t \\ q^{n+1} = q^{trial} - G \dot{\gamma}_{_S} \Delta t \end{gathered}\right\}. \end{equation}Together with the requirement that the new stress state resides on the (updated) yield surface, i.e.
\begin{equation} y(p^{n+1},q^{n+1})=0 {,} \end{equation}
where (A16) represents the projection of trial stress state 
$(p^{trial}, q^{trial})$ to the closest stress state 
$(p^{n+1},q^{n+1})$ on the yield surface. This projection operation is often referred to as a return mapping and is in the general case solved with iterative schemes. In this work, (A16) and (A17) were solved with the Newton–Raphson method. Recall that, from the proposed hardening law of § 2.4, the compressive strength 
$p_c$ given by (2.19) depends on plastic volumetric deformation only. The plastic volumetric strain 
$\varepsilon _V^P$ increases or decreases only as a result of change in pressure, in particular the plastic volume change is given by
\begin{equation} \Delta \varepsilon_V^{P, n+1} = \frac{p^{n+1}-p^{trial}}{K} {.} \end{equation}
Thus, 
$p_c$ is independent of 
$q$. Further assuming 
$\mu = \mu _1$ is constant (recall discussion in § 2.6), it follows that 
$\partial y/\partial q$ is independent of 
$p$ and that 
$\partial y/\partial p$ is independent of 
$q$.
Appendix B. Influence of elastic and plastic parameters
Here, a condensed sensitivity study on the elastic and plastic parameters is presented, focusing both on granular column collapse and steady-state flow.
In the case of granular collapse, the dam break set-up from § 4.2 on slope angle 
$\theta =22^\circ$ with the same initial conditions and parameters as in table 5 is considered, albeit in two spatial dimensions to reduce computational time (using 
$N_p=17$k particles, the two-dimensional simulations take only 1.3 s instead of 18 h per second). The effect of changing the Young's modulus 
$E$, Poisson's ratio 
$\nu$, initial compressive strength 
$p_c^0$ and hardening factor 
$\xi$ on the front position evolution is demonstrated in figure 20. Varying the Young's modulus between 
$1$ MPa and 
$1$ GPa, it can be observed in figure 20(a) that it does not significantly change the collapse evolution. The upper value 
$1$ GPa was used by Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015) in their elastoplastic granular flow simulations. Moreover, increasing the Poisson's ratio means allowing the material to expand more, and accordingly a very small increase in front position may be observed with increasing Poisson‘s ratio. While a Poisson's ratio 
$\nu =0.3$ was used throughout this paper, for spherical glass beads 
$\nu =0.2$ may be more appropriate (Muqtadir, Al-Dughaimi & Dvorkin Reference Muqtadir, Al-Dughaimi and Dvorkin2020). Nevertheless, as shown in figure 20(b), 
$\nu =0.2$ and 
$\nu =0.3$ gave approximately indistinguishable results for the granular collapse. With a small hardening factor 
$\xi$, the system may evidence a significant initial compaction under gravity. As shown in figure 20(d), as long as it is large enough (
$\xi >10$), the collapse displays an almost identical dynamics of the front position. Figure 21 demonstrates further how the hardening factor influences the ability of the material to compact, clearly highlighting the reduction of the initial gravity-induced compaction with increasing 
$\xi$. Figure 22 indicates that the volumetric deformations are dominantly plastic. Using a very small initial compressive strength 
$p_c^0$, leaving 
$\xi$ unchanged, may also lead to strong compaction under gravity. In an appropriate range of values, changing 
$p_c^0$ does not affect the collapse significantly, as can be seen from figure 20(c). On the other hand, with too large 
$p_c^0$, the material may behave more as a fracturing solid. This is illustrated in figure 23, which shows that when using 
$p_c^0=10$ kPa one observes one large block of material fracturing before sliding off. At 
$t=0.2$ s, increasing 
$p_c^0$ results in a behaviour that may appear reminiscent of the behaviour for increasing 
$\beta$. However, while increasing 
$\beta$ results in a monotonically decreasing front position, increasing 
$p_c^0$ will not alter the runout dynamics unless it is so large (
$\,p_c^0=10$ kPa) as to induce solid-like fracturing. Visualizing the plastic shear strain rate, figure 24 suggests that the surface effects discussed in § 4.2 can be controlled through the plastic parameters. As can be seen in this figure, by tuning 
$p_c^0$ and 
$\xi$ in the cohesionless case one is able to suppress the surface bumps substantially. Nevertheless, as will be discussed in Appendix C, the surface effects are discretization dependent as they originate from strain localization.
Front position evolution in two-dimensional dam break on 
$\theta =22^\circ$ with various elastic and plastic parameters. Default parameters: Young's modulus 
$E=1$ MPa, Poisson's ratio 
$\nu =0.3$, initial compressive strength 
$p_c^0=100$ Pa and hardening factor 
$\xi =50$. Experimental data of Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010).
Solid fraction evolution in two-dimensional cohesionless dam break on 
$\theta =22^\circ$ with various hardening parameters 
$\xi$.
Elastic and plastic volume change in the cohesionless dam break on 
$\theta =22^\circ$ from § 4.2.
Cohesionless dam break on 
$\theta =22^\circ$ at times 
$t=0.2$ and 
$t=0.6$ with different values of initial compressive strength 
$p_c^0$. To highlight the shear bands, they are here visualized in terms of the equivalent plastic shear strain rate.
Cohesionless dam break on 
$\theta =22^\circ$ at time 
$t=0.2$ with different plastic parameters resulting in different surface features: (a) 
$p_c^0 = 0.1$ kPa, 
$\xi =50$ and (b) 
$p_c^0 = 0.2$ kPa, 
$\xi =70$. Here visualized in terms of the equivalent plastic shear strain rate.
In order to investigate the effect of 
$p_c^0$ further, a simple direct shear test is conducted in the absence of gravity. It is expected from critical state theory that the transition to critical state (where no further volume change occurs) will depend on whether the initial state is loose or dense (Wood Reference Wood1991). Considering loose and dense systems through small and large initial compressive strengths, respectively, a differing dynamics is measured, as shown in figure 25. Dense systems compact slightly before dilating while loose systems only compact, in both cases reaching a state with no further volumetric deformations. This observation is consistent with expectations, with similar results also obtained through previous discrete element simulations (Bernhardt, Biscontin & O'Sullivan Reference Bernhardt, Biscontin and O'Sullivan2016; Wang et al. Reference Wang, Li and Liu2022).
Evolution of the volumetric (Hencky) strain 
$\varepsilon _V$ during a direct shear test. The set-up is sketched in the inset figure. Here using 
$E=1$ MPa, 
$\nu =0.3$, 
$\xi =100$, 
$\beta =0$ and the parameters given in table 3.
In the case of steady-state flow, the majority of the stress states would align with the CSL, and the exact onset of plastic deformation is not crucial as long as it is small enough to ensure all deformations are dominantly plastic. In that case, the elastic regime is generally not important. Figure 26 shows the (in)sensitivity of elastic and plastic parameters in capturing the expected Bagnold steady-state velocity profile. As such, in this situation it can be claimed that the elastic model is not intended as a constitutive law supposed to give a quantitative description of the material, rather simply providing a convenient means to deal with the ill-posedness of the 
$\mu (I)$-rheology for low inertial numbers.
Steady-state flow profile reached over time 
$t$ on inclination 
$\theta =28^\circ$. Here, various elastic and plastic parameters are changed, where in all cases the initial condition is a linear velocity profile. The black dashed curve is the theoretical prediction from the 
$\mu (I)$-rheology, using the parameters given in table 7.
Appendix C. Strain localization and discretization dependency
When exposed to external stress or during the transition to flow, granular media typically display shear bands where the deformation localizes in thin zones that can have a width of a few grains (Mueth et al. Reference Mueth, Debregeas, Karczmar, Eng, Nagel and Jaeger2000; Fenistein & van Hecke Reference Fenistein and van Hecke2003). Strain localization and its dependency on discretization is a general feature of continuum elastoplastic models featuring strain softening. As a remedy, rigorous regularizing techniques introduce a characteristic length scale in the governing equations, e.g. through higher-order gradient methods (Triantafyllidis & Aifantis Reference Triantafyllidis and Aifantis1986; Aifantis Reference Aifantis1987, Reference Aifantis1992) or with the introduction of non-local variables through spatially weighted averaging (Jin & Arson Reference Jin and Arson2018; Monforte et al. Reference Monforte, Ciantia, Carbonell, Arroyo and Gens2019). In figure 27, the granular collapse is visualized at different grid resolutions. With increasing resolution, it is observed that the plastic strain becomes increasingly localized in these bands, as expected in such an elastoplastic model, with the thickness of the bands decreasing and the magnitude of the plastic strain in these bands increasing. As such, a proper experimental validation of the simulated bands would not be feasible. Nevertheless, the dynamics of the collapse remains independent of discretization, as shown in figure 28, demonstrating convergence of the front position evolution under grid refinement. Under flow on inclined planes with inclination 
$\theta > \tan ^{-1}(\mu _1)$, where all stress states are fully plastic and where there is no static–fluid interface, strain localization is not observed. In this case, figure 29 shows the discretization-independent results for the velocity profile.
Cohesionless dam break (corresponding to figure 24b) on 
$\theta = 22^\circ$ at time 
$t=0.2$ under grid refinement. Here, 
$\Delta x_0 = h_0/45$ using initially six particles per grid cell in all cases.
Evolution of front position in the cohesionless dam break on 
$\theta = 22^\circ$ under grid refinement, in all cases using initially six particles per grid cell. Here, (a,b) correspond to the dam break presented in figures 24(a) and 24(b), respectively.
Cohesionless flow on inclination 
$\theta =28^\circ$ under grid refinement, using parameters as in table 7. The velocity profile was initially imposed as linear (with a surface velocity of 
$4\ \textrm {m}\ \textrm {s}^{-1}$) and is here shown at two different later times, the right plot being the steady-state profile.
Appendix D. Results from rate-independent CSSM
Previous studies have shown the ability of 
$I$-independent models with constant friction or constant viscosity in capturing granular column collapse (Crosta et al. Reference Crosta, Imposimato and Roddeman2009; Ionescu et al. Reference Ionescu, Mangeney, Bouchut and Roche2015; Martin et al. Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017; Rousseau et al. Reference Rousseau, Métivet, Rousseau, Daviet and Bertails-Descoubes2022). It is therefore instructive to consider to which extent the same can be claimed for rate-independent CSSM, in particular the proposed model with a constant, 
$I$-independent, 
$\mu$. Indeed, figure 30 indicates that with the particular choice 
$\mu =\tan (26^\circ )$, the granular collapse front position is well captured.
Granular collapse on slope angle 
$\theta =22^\circ$ using rate-independent MCC with constant slope 
$\mu$ of CSL, otherwise using the parameters as in § 4.2 treating the rate-dependent case. Experimental data of Mangeney et al. (Reference Mangeney, Roche, Hungr, Mangold, Faccanoni and Lucas2010).
However, in capturing steady-state flow, the 
$I$-dependence is crucial. Figure 31 shows the inability of rate-independent MCC in achieving steady-state flow, either accelerating flow in the case 
$\theta >\tan ^{-1}(\mu )$ or stopping flows in the case 
$\theta <\tan ^{-1}(\mu )$.
Flow on slope angle 
$\theta =28^\circ$ using rate-independent MCC with constant slope of CSL, with (a) 
$\mu = \tan (35^\circ )$ and with (b) 
$\mu = \tan (21^\circ )$, showing stopping flows and accelerating flows, respectively. The other parameters are chosen as that of § 4.1 where the rate-dependent model produced steady-state flows. In both (a,b), the initial condition was here a linear velocity profile.
Appendix E. Basal friction
In order to validate the boundary conditions used in § 4.4 where a non-zero basal friction was introduced, simulations of a stiff elastic box sliding on inclined planes are performed. Figure 32 demonstrates that the proposed approach can accurately capture the analytical solutions from Coulomb's friction law.
Downslope position 
$x$ of a stiff elastic box sliding down planes with inclination 
$\theta =14^\circ, 16^\circ, 18^\circ,\ldots, 28^\circ, 30^\circ$, in all cases using a basal friction coefficient 
$\mu _b = \tan (15^\circ )$. As such, no flow is expected for inclinations 
$\theta$ below 
$15^\circ$. The initial position is 
$x=0$ with zero velocity. The dots represent the analytical solution 
$x = \textrm {max}(0, \frac {1}{2}g (\sin \theta - \mu _b \cos \theta ) t^2)$ and the solid lines represent the simulations.
Open access
