2.1. A brief recall of the model of Deleage & Richard (Reference Deleage and Richard2025)

Figure 1.

Definition sketch.

Figure 1. Long description

A diagram illustrating granular flow models. The diagram features a three-dimensional coordinate system with axes labeled z, x, and g. It shows a collection of circular particles arranged in a staggered pattern. Key components include height (h), distance (d), length (L), and an angle (θ). Arrows indicate the direction of flow and movement within the granular material. The diagram also includes a curved path labeled u, representing the flow direction. The relationships between the particles and the directional arrows illustrate the interactions and collective effects in granular flow.

The model of Deleage & Richard (Reference Deleage and Richard2025) is a three-equation model for a granular flow down an inclined plane of angle $\theta$ under the action of gravity $g$ (see figure 1). The grain diameter $d$ , the particle density $\rho _p$ and the solid volume fraction $\phi$ are assumed to be constant. The model is derived consistently from the incompressible $\mu (I)$ -rheology, defined by the three parameters $\mu _1, \mu _2$ and $I_0$ , such that the internal friction coefficient $\mu$ is given as a function of the inertial number $I$ as (Jop et al. Reference Jop, Forterre and Pouliquen2006)

(2.1) \begin{align} \mu (I) = \mu _1 + \frac {\mu _2-\mu _1}{1+\frac {I_0}{I}}, \quad \text{with} \quad I = \frac {|\boldsymbol{\dot {\gamma }}|d}{\sqrt {p/\rho _p}}. \end{align}

The scalars $|\boldsymbol{\dot {\gamma }}|$ and $p$ respectively stand for the local shear and pressure inside the flow.

The three variables of the model are the flow depth $h$ , the mean velocity $U$ and the enstrophy $\psi$ , which is related to the shear by

(2.2) \begin{align} \psi = \frac {\langle \left ( u- U \right )^2 \rangle }{h^2} \end{align}

where $u$ is the fluid velocity in the $Ox$ -direction (see figure 1) and $\langle X \rangle$ denotes the depth average of the quantity $X$ , defined by

(2.3) \begin{align} \langle X \rangle = \frac {1}{h}\int _0^h X\; \mathrm{d}z. \end{align}

Note that $\langle u \rangle =U$ . The evolution of the model variables is governed by the following system:

(2.4) \begin{align} \frac {\partial h}{\partial t}+ \frac {\partial hU}{\partial x}=0, \\[-28pt] \nonumber \end{align}

(2.5) \begin{align} \frac {\partial hU}{\partial t} \kern-1pt + \kern-1pt \frac {\partial }{\partial x} \!\left (\! hU^2+\beta h^3\psi +\alpha _p\frac {gh^2}{2}\cos \theta \! \right ) =g\cos \theta \! \left [h(\tan \theta -\mu _b)+\frac {10\lambda }{A^2B} \!\left (\! \psi -\frac {U^2}{4h^2} \!\right ) \!\right ]\!, \\[-28pt] \nonumber \end{align}

(2.6) \begin{align} \frac {\partial h\psi }{\partial t}+ \frac {\partial hU\psi }{\partial x}= \frac {2gU\cos \theta }{h^2}\left [\frac {\alpha _1}{h^{1/2}}\left (U-\frac {2}{5}Ah^{3/2}\right )+\alpha _2\left (\psi -\frac {U^2}{4h^2}\right )\right ]\!. \\[0pt] \nonumber \end{align}

Equation (2.4) is the depth-averaged mass equation, (2.5) is the depth-averaged horizontal momentum equation and (2.6) is the enstrophy equation. These three equations are written in conservative form. The expressions of the coefficients $A,B,C,\alpha _1$ and $\alpha _2$ are

(2.7) \begin{align} A=I_0\frac {\tan \theta -\mu _1}{\mu _2-\tan \theta }\frac {\sqrt {\phi g\cos \theta }}{d},\quad B=\frac {\mu _2-\mu _1}{(\mu _2-\tan \theta )(\tan \theta -\mu _1)}, \\[-28pt] \nonumber \end{align}

(2.8) \begin{align} C=\tan \theta \left (2\tan \theta -\frac {1}{B}\right )\!, \\[-28pt] \nonumber \end{align}

(2.9) \begin{align} \alpha _1=-\frac {33C}{4AB(17C+14)}, \quad \alpha _2 = \frac {231(3C-4)}{4A^2B(17C+14)}. \\[0pt] \nonumber \end{align}

The parameters $\lambda ,\beta$ and $\alpha _p$ are given by

(2.10) \begin{align} \lambda = \frac {847(4-3C)}{20(157C+84)},\quad \beta = \frac {11}{16} + \frac {5\lambda }{44}, \quad \alpha _p=1-\frac {3C}{4}+\frac {\lambda C}{7}, \end{align}

and the basal friction $\mu _b$ is given by

(2.11) \begin{align} \mu _b = \mu _1 + \frac {\mu _2-\mu _1}{I_0\frac {2h^{3/2}\sqrt {\phi g\cos \theta }}{5dU}+1}. \end{align}

The coefficients $\beta$ and $\alpha _p$ (instead of $1$ ) and the relaxation term with the coefficient $\lambda$ in (2.5) result from first-order corrections that ensure the model’s consistency. They are obtained using the asymptotic expansions detailed in Deleage & Richard (Reference Deleage and Richard2025).

These equations form a hyperbolic system of conservation laws with relaxation source terms if $C\lt 4/3$ . The case $C\gt 4/3$ is discussed briefly below and in more detail in Deleage & Richard (Reference Deleage and Richard2025). The system is compatible with the energy balance equation

(2.12) \begin{align} \begin{aligned} &\frac {\partial he}{\partial t}+\frac {\partial }{\partial x}\left (hUe+\varPi U\right ) = gU\cos \theta \\ &\times \left [h(\tan \theta -\mu _b)+\frac {\beta \alpha _1}{h^{1/2}}\left (U-\frac {2}{5}Ah^{3/2}\right )+\left (\beta \alpha _2+\frac {10\lambda }{A^2B}\right )\left (\psi -\frac {U^2}{4h^2}\right )\right ] \!, \end{aligned} \end{align}

where the specific energy is

(2.13) \begin{align} e=\frac {U^2}{2}+\frac {\beta h^2\psi }{2}+\alpha _p\frac {gh}{2}\cos \theta \end{align}

and

(2.14) \begin{align} \varPi = \beta h^3\psi + \alpha _p\frac {gh^2}{2}\cos \theta . \end{align}

Since this model is consistent with the $\mu (I)$ -rheology, it takes into account steady uniform flows with a Bagnold velocity profile and gives the same instability threshold for granular roll-wave instability as the instability threshold of the $\mu (I)$ -rheology computed by Forterre (Reference Forterre2006) (Deleage & Richard Reference Deleage and Richard2025). However, the $\mu (I)$ -rheology only allows the existence of such steady uniform flows for an angle $\theta$ such that

(2.15) \begin{align} \mu _1\lt \tan \theta \lt \mu _2. \end{align}

Indeed, these flows can only occur when the weight is compensated by the friction, which gives the equilibrium condition $\mu = \tan \theta$ . For the standard form of the $\mu (I)$ curve (2.1), the friction coefficient $\mu$ is restricted to the interval $(\mu _1,\mu _2)$ . Hence, no steady uniform flow can be obtained when the condition (2.15) is not satisfied (apart from uniform flow at rest). As a consequence, the definitions (2.7)–(2.10) of the coefficients of the model are only valid when condition (2.15) holds true. For an inclination angle out of the range $(\mu _1,\mu _2)$ , the model (2.4)–(2.6) cannot be used. We refer to Deleage & Richard (Reference Deleage and Richard2025) for more details concerning these limitations of system (2.4)–(2.6). Of course, the limitations of system (2.4)–(2.6) result from the limitations of rheology $\mu (I)$ , both for quasi-static and diluted regimes. These problems could be resolved through approaches that include more physical effects. Partial regularisations have been proposed to this end (see, for example, Barker & Gray Reference Barker and Gray2017).

Another problem of system (2.4)–(2.6), which is also inherited from the $\mu (I)$ rheology, can occur even when the condition (2.15) is satisfied, and is linked with the value of the parameter $\mu _2$ . Forterre (Reference Forterre2006) showed that, for the $\mu (I)$ -rheology, the roll-wave instability occurs when the Froude number $Fr=|U|/\sqrt {gh_0\cos \theta }$ is greater than a critical Froude number $Fr_c$ given by

(2.16) \begin{align} Fr_c^2=\frac {16}{25}(4-3C), \end{align}

and $C$ is given by (2.8). Deleage & Richard (Reference Deleage and Richard2025) showed that, if $\mu _2\gt \sqrt {2/3}$ , there are some values of $\theta$ such that (2.15) is satisfied and $4-3C\lt 0$ , which implies that every steady uniform flow is unstable, whatever the Froude number is (see (2.16)). In system (2.4)–(2.6), this unstable behaviour is materialised by the fact that the coefficients $\alpha _p$ and $\alpha _2$ change sign, leading to a model that can be ill posed. In particular, if $C\gt 4/3$ the hyperbolicity of the system is not guaranteed. In the next subsection, we propose a regularised version of system (2.4)–(2.6), which is defined and well posed for every inclination angle $\theta$ .

2.2. Regularised model

The need to regularise the system comes from the fact that some relaxation terms, which are essential for ensuring the consistency of the model, guarantee good properties near equilibrium but can lead to poor behaviour far from equilibrium. This is particularly the case when no equilibrium is possible (i.e. if the condition (2.15) is not satisfied), near a dry front or when flow stoppage occurs. The conservative part of the equations always remains well posed and therefore does not need to be regularised. The regularisation method therefore consists of switching off the relaxation terms that cause problems far from the model’s domain of validity and replacing them with terms that behave well under these conditions. The method consists of defining a regularisation parameter $f$ that allows an interpolated model to be written between the limit case of the consistent model, corresponding to $f=1$ , and the limit case of the model completely out of equilibrium, corresponding to $f=0$ . To preserve the consistency of the model near equilibrium up to $O(\varepsilon )$ , the coefficient $f$ must be of the form $f=1+O(\varepsilon ^2)$ .

The consistent model is derived using an asymptotic method, which is algorithmic in nature. Nothing of the sort is possible far from equilibrium or when no equilibrium is possible. The alternative model must therefore be written to ensure that it remains well posed and physically relevant. Its parameters are determined empirically by comparison with experimental or numerical results. In some cases, particularly for $\tan \theta \geq \mu _2$ , there are very few results available to calibrate the model. In such cases, simple arbitrary choices must be made, at least provisionally, ensuring that the model behaves reasonably well, at least in terms of numerical resolution, in all circumstances. These arbitrary choices can always be refined in the future if new experimental work requires it. It should also be mentioned that, in the case $\tan \theta \geq \mu _2$ , the assumption of incompressibility is not necessarily relevant due to the importance of dilatancy phenomena, which would imply revisiting the entire modelling process to include dilatancy.

Therefore system (2.4)–(2.6) is regularised with the coefficient $f$ defined by

(2.17) \begin{align} f=\left \{\begin{aligned} &\exp \left (2-\frac {25U^2}{4A^2h^3}-\frac {4A^2h^3}{25U^2}\right ), & \mu _1\lt \tan \theta \lt \mu _2, \\ &0, & \text{otherwise}. \end{aligned} \right . \end{align}

It follows from (2.17) that $0\leq f\leq 1$ , and that $f=0$ in the following three cases:

1. (i) when (2.15) is not satisfied;

2. (ii) when $h=0$ , i.e. in a dry front;

3. (iii) when $U=0$ , i.e. when the flow stops.

For a steady uniform flow, the depth-averaged velocity and the depth are linked via the relation $U=2Ah^{3/2}/5$ (Pouliquen Reference Pouliquen1999b ), which gives $f=1$ . The coefficient $f$ may thus be seen as a quantification of the distance between the solution of system (2.4)–(2.6) and a steady uniform flow.

The coefficient $f$ can thus regularise the model as an interpolation between system (2.4)–(2.6) when $f=1$ and another system when $f=0$ . When doing so, the consistency between system (2.4)–(2.6) and the $\mu (I)$ -rheology must be preserved by the regularisation procedure. Indeed, let us introduce as usual the typical height of the flow $h_0$ , and a typical horizontal length $L_0$ . The shallow-water parameter $\varepsilon$ is defined by

(2.18) \begin{align} \varepsilon = \frac {h_0}{L_0}. \end{align}

Deleage & Richard (Reference Deleage and Richard2025) showed that system (2.4)–(2.6) is consistent at order one with the $\mu (I)$ -rheology, which means that an exact solution of the governing equations of the $\mu (I)$ -rheology yields an approximate solution of system (2.4)–(2.6), up to an $O(\varepsilon ^2)$ remainder. Let us now suppose that condition (2.15) is satisfied. We introduce the dimensionless variables $\tilde {h}$ and $\tilde {U}$ given by

(2.19) \begin{align} \tilde {h}=\frac {h}{h_0}, \quad \tilde {U}=\frac {5U}{2Ah_0^{3/2}}. \end{align}

It follows that the coefficient $f$ can be written

(2.20) \begin{align} f=\exp \left (2-\frac {\tilde {U}^2}{\tilde {h}^3}-\frac {\tilde {h}^3}{\tilde {U}^2}\right )\!. \end{align}

Then assuming, as in Deleage & Richard (Reference Deleage and Richard2025), that $\tilde {U}$ satisfies an asymptotic expansion, given by

(2.21) \begin{align} \tilde {U}=\tilde {U}^{(0)}+\varepsilon \tilde {U}^{(1)}+O(\varepsilon ^2)=\tilde {h}^{3/2}+\varepsilon \tilde {U}^{(1)}+O(\varepsilon ^2), \end{align}

we deduce that

(2.22) \begin{align} f=\exp \big [O(\varepsilon ^2)\big ]=1+O(\varepsilon ^2). \end{align}

Therefore the error terms are at most $O(\varepsilon ^2)$ , and the consistency of the regularised model is preserved.

It remains to specify the system to be used far from steady uniform flows, i.e. when $f=0$ . In order to do this, the (2.4)–(2.6) are treated separately.

Equation (2.4) is the equation of conservation of mass. As it is expected to be valid in every flow regime, this equation is unchanged in the regularised version of the model.

Equation (2.5) is the balance of momentum. The left-hand side of the equation is not modified in in the regularised version, but the expression of the right-hand side is changed. The basal friction $\mu _b$ is replaced by the non-monotonic friction law denoted here $\mu (h,U)$ that was proposed by Pouliquen & Forterre (Reference Pouliquen and Forterre2002) and given by

(2.23) \begin{align} \mu (h,U)=\left \{ \begin{aligned} &\mu _b, & \frac {|U|}{\sqrt {gh\cos \theta }}\gt \beta ',\\ &\left (\frac {|U|}{\beta '}\sqrt {gh\cos \theta }\right )^\gamma (\mu _1-\mu _3) + \mu _3 + \frac {\mu _2-\mu _1}{\frac {h}{L}+1}, & 0\lt \frac {|U|}{\sqrt {gh\cos \theta }}\lt \beta ',\\ &\min \left (\mu _3 + \frac {\mu _2-\mu _1}{\frac {h}{L}+1},\left |\tan \theta -\alpha _p\frac {\partial h}{\partial x}\right |\right ), & U=0. \end{aligned} \right . \end{align}

Here, the constant $\beta '$ is a parameter appearing in the scaling

(2.24) \begin{align} \frac {|U|}{\sqrt {gh\cos \theta }}=Fr=\beta '\frac {h}{h_{\text{stop}}(\theta )}, \end{align}

and characterising the transition between the static and dynamic regimes. In this expression, $h_{\text{stop}}(\theta )$ is the depth of the static layer left by a steady uniform flow (see Pouliquen (Reference Pouliquen1999b )). An advantage of this non-monotonic law is that it captures frictional hysteresis. A drawback is that it involves, in addition to the three parameters $\mu _1,\mu _2,I_0$ of the rheology, four other parameters ( $\mu _3,\beta ', L,\gamma$ ) which are not always determined by experimental measurements. In order to be able to compare the performances of our model with laboratory data even when these parameters are not specified, we follow Poulain et al. (Reference Poulain2023) and set in this case $\beta '=0.136$ , $L=1.3d$ , $\gamma =10^{-3}$ and $\mu _3=\tan \theta _3$ , where $\theta _3=\theta _1+1.2^\circ$ and $\tan \theta _1=\mu _1$ .

Another non-monotonic friction law, which accounts for frictional hysteresis, was proposed by Edwards et al. (Reference Edwards, Viroulet, Kokelaar and Gray2017, Reference Edwards, Russell, Johnson and Gray2019), who introduced a parameter $h^*$ representing the minimum observable steady uniform flow thickness. From their experimental observations, $h^*\gt h_{\text{stop}}$ ( $h^* \approx 1.33 h_{\text{stop}}$ for glass beads, see Russell et al. Reference Russell, Johnson, Edwards, Viroulet, Rocha and Gray2019). They thus introduced a new parameter $\beta _* = 1.33\beta '$ and modified the friction law of Pouliquen & Forterre (Reference Pouliquen and Forterre2002) so that the transition between the two frictional behaviours occurs at $Fr=\beta _*$ for a steady uniform flow. As a consequence, their friction law, that we will denote here $\mu _*(h,U)$ , reads

(2.25) \begin{align} \mu _*(h,U)=\left \{ \begin{aligned} &\mu _b, & \frac {|U|}{\sqrt {gh\cos \theta }}\gt \beta _*,\\ &\left (\frac {|U|}{\beta _* \sqrt {gh\cos \theta }}\right )^\gamma \left (\mu _1 + \frac {\mu _2-\mu _1}{\frac {h\beta }{L\beta _*} + 1}-\mu _3 - \frac {\mu _2-\mu _1}{\frac {h}{L}+1}\right ) \\ & + \mu _3 + \frac {\mu _2-\mu _1}{\frac {h}{L}+1}, & 0\lt \frac {|U|}{\sqrt {gh\cos \theta }}\lt \beta _*,\\ &\min \left (\mu _3 + \frac {\mu _2-\mu _1}{\frac {h}{L}+1},\left |\tan \theta -\alpha _p\frac {\partial h}{\partial x}\right |\right ), & U=0. \end{aligned} \right . \end{align}

Another difference with the friction law of Pouliquen & Forterre (Reference Pouliquen and Forterre2002) is that Edwards et al. (Reference Edwards, Viroulet, Kokelaar and Gray2017) use the value $\gamma =1$ for the interpolation coefficient, while Pouliquen & Forterre (Reference Pouliquen and Forterre2002) use the value $\gamma = 10^{-3}$ . The law of Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019) was used by Russell et al. (Reference Russell, Johnson, Edwards, Viroulet, Rocha and Gray2019) to model retrogressive failures. We choose here to keep the friction law (2.23) that was proposed by Pouliquen & Forterre (Reference Pouliquen and Forterre2002). The difference between the friction law $\mu _b$ of Pouliquen (Reference Pouliquen1999b ), the law $\mu (h,U)$ of Pouliquen & Forterre (Reference Pouliquen and Forterre2002) and the law $\mu _*(h,U)$ of Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019) is studied in § 6.2. The law of Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019) is used in § 7 to model a retrogressive failure.

The relaxation term with a factor $\lambda$ in the momentum equation ensures consistency in order to obtain the right instability threshold near a steady uniform flow. When condition (2.15) is not true, there is no purpose in keeping this term and, moreover, this term can lead to an ill-posed model in some situations far from equilibrium. This term is thus multiplied by the coefficient $f$ and set equal to zero when $\tan \theta$ is outside the range $(\mu _1,\mu _2)$ .

Equation (2.6) is the equation of enstrophy. The left-hand side does not need to be regularised since it is the conservative part of the equation. However, the right-hand side is made of relaxation terms which ensure consistency near stationary and uniform flows but which could lead to bad behaviour far from equilibrium. As in the momentum equation, these terms are multiplied by the factor $f$ in order to set them equal to zero when $f=0$ . Recall that, at equilibrium, the enstrophy $\psi$ is equal to $U^2/(4h^2)$ , corresponding to a Bagnold velocity profile, and to a shape factor $\alpha =5/4$ (the shape factor can be defined with the help of the enstrophy by $\alpha =1+h^2\psi /U^2$ ). Far from equilibrium, few data on the velocity profile or on the shape factor of a granular flow are available. However, Lagrée et al. (Reference Lagrée, Saingier, Deboeuf, Staron and Popinet2017) showed that, in a dry front, the shape factor is equal to $5/4$ , except close to the front where it goes to 1 (see § 3). Therefore, the enstrophy still relaxes towards $U^2/(4h^2)$ (i.e. a Bagnold profile) even far for equilibrium, except close to a dry front, where the enstrophy is required to relax towards 0, which corresponds to a shape factor of 1. A way to do this is to impose the following enstrophy equation when $f=0$ :

(2.26) \begin{align} \frac {\partial h\psi }{\partial t}+\frac {\partial hU\psi }{\partial x} = -\frac {d}{\tau }\left (\psi -e^{-\frac {5d}{h}}\frac {U^2}{4h^2}\right )\!, \end{align}

where $\tau$ is a relaxation time. Indeed, far from a dry front, $h\gg d$ , hence $e^{-5d/h}\approx 1$ and the enstrophy relaxes towards $U^2/(4h^2)$ (a Bagnold profile). Close to a dry front, $h\to 0$ , hence $e^{-5d/h}$ goes to zero and the enstrophy relaxes towards 0 (a shape factor equal to 1). The factor 5 in the exponential is determined empirically by comparison with the results of Saingier et al. (Reference Saingier, Deboeuf and Lagrée2016). The sensitivity to this factor is studied in § 6.1.

The relaxation time $\tau$ cannot be the same as in (2.6), since it depends on the coefficient $\alpha _2$ which is not defined when $\tan \theta$ is not between $\mu _1$ and $\mu _2$ . However, the equilibrium expression $U=2Ah^{3/2}/5$ can be used to write the relaxation time of (2.6) as

(2.27) \begin{align} \!-\frac {1}{\tau }=\frac {2gU\alpha _2\cos \theta }{d h^2}\approx \frac {4gA\alpha _2\cos \theta }{5d h^{1/2}}=\frac {4}{5I_0}\sqrt {\frac {g\cos \theta }{\phi h}}\frac {231(3C-4)}{4(17C+14)}\frac {(\mu _2-\tan \theta )^2}{\mu _2-\mu _1}. \end{align}

This expression could be kept to define $\tau$ , but it vanishes when $\tan \theta =\mu _2$ . In order to avoid this, the same expression is used, but the factor $(\mu _2-\tan \theta )^2$ is arbitrarily replaced by $(\mu _2-\mu _1)^2$ . As the behaviour of the system is not very sensitive to the precise value of $\tau$ and given the scarcity of experimental results for $\tan \theta \geq \mu _2$ , this choice gives a plausible order of magnitude and has no other consequence than to guarantee the relaxation of the enstrophy for $\tan \theta \geq \mu _2$ .

To obtain reasonable values of $C$ when the slope does not satisfy the inequalities (2.15), the definition of $C$ must be expanded to satisfactorily cover all possible cases. For $\tan \theta = \mu _1$ , $C=2\mu _1^2$ and for $\tan \theta = \mu _2$ , $C=2\mu _2^2$ . In the range of slopes $\mu _1 \leq \tan \theta \leq \mu _2$ , where uniform steady flows are possible, $C$ satisfies $2 \mu _1^2 \leq C \leq 2 \mu _2^2$ . The simplest and most straightforward assumption is to take the value of $C$ for $\tan \theta = \mu _1$ if $\tan \theta \lt \mu _1$ and the value of $C$ for $\tan \theta = \mu _2$ if $\tan \theta \gt \mu _2$ .

As recalled above and explained in Deleage & Richard (Reference Deleage and Richard2025), the system can become ill posed if $C \geq 4/3$ (the hyperbolicity is not guaranteed anymore), which can happen if $\mu _2 \geq \sqrt {2/3}$ . To avoid this problem, the value of $C$ is capped at the arbitrary value of 1 if $\mu _2 \geq \sqrt {2/3}$ (see the discussion below). If $\mu _2 \lt \sqrt {2/3}$ , this problem cannot arise and, if $\tan \theta \gt \mu _1$ , $C$ has its consistent expression (2.8) if $\tan \theta \lt \mu _2$ and the value $2 \mu _2^2$ if $\tan \theta \gt \mu _2$ . In the latter case, the consistent expression (2.8) gives a larger value than $2 \mu _2^2$ since $B \leq 0$ in this case.

Therefore the complete definition of $C$ is

(2.28) \begin{align} C=\left \{ \begin{aligned} & 2\mu _1^2 & \text{if }\tan \theta \leq \mu _1, \\[5pt]& \min \left [2\mu _2^2,\tan \theta \left (2\tan \theta -\frac {1}{B}\right )\right ] & \text{if } \mu _1\lt \tan \theta \quad \text{and} \quad \mu _2\lt \sqrt {\frac {2}{3}}, \\[5pt]&\min \left [1,\tan \theta \left (2\tan \theta -\frac {1}{B}\right )\right ] & \text{if } \mu _1\lt \tan \theta \quad \text{and} \quad \mu _2\geq \sqrt {\frac {2}{3}}. \end{aligned} \right . \end{align}

Comparisons with the results of Lagrée et al. (Reference Lagrée, Saingier, Deboeuf, Staron and Popinet2017) show that the transition from a shape factor $5/4$ to a shape factor $1$ close to a dry front is not satisfactorily predicted if the regularisation coefficient $f$ is used as such in the enstrophy equation. It can be determined empirically that the coefficient $f$ must be replaced in the enstrophy equation by $f^\delta$ , where the exponent $\delta$ is equal to $\delta = 1/20$ , to obtain a good agreement (see § 3). On the other hand, it is possible to keep the coefficient $f$ with no exponent in the momentum equation. The sensitivity of our results to the value of $\delta$ is studied in § 6.1.

The regularised system can finally be written

(2.29) \begin{align} \frac {\partial h}{\partial t}+ \frac {\partial hU}{\partial x}=0, \\[-28pt] \nonumber \end{align}

(2.30) \begin{align} &\frac {\partial hU}{\partial t}+ \frac {\partial }{\partial x}\left (hU^2+\beta h^3\psi +\alpha _p\frac {gh^2}{2}\cos \theta \right )\nonumber \\& \quad =g\cos \theta \left [h[\tan \theta - \mu (h,U)]+\frac {10\lambda f}{A^2B}\left (\psi -\frac {U^2}{4h^2}\right )\right ], \\[-28pt] \nonumber \end{align}

(2.31) \begin{align} & \frac {\partial h\psi }{\partial t}+ \frac {\partial hU\psi }{\partial x}= \frac {2f^\delta gU\cos \theta }{h^2}\left [\frac {\alpha _1}{h^{1/2}}\left (U-\frac {2}{5}Ah^{3/2}\right )+\alpha _2\left (\psi -\frac {U^2}{4h^2}\right )\right ]\nonumber \\& \phantom {\frac {\partial h\psi }{\partial t}+ \frac {\partial hU\psi }{\partial x}=} +(1-f^\delta )\frac {\alpha _3}{h^{1/2}}\left (\psi -e^{-\frac {5d}{h}}\frac {U^2}{4h^2}\right ). \\[0pt] \nonumber \end{align}

The coefficients $A,B,\alpha _1,\alpha _2,\lambda ,\alpha _p,\beta$ are defined as before by (2.7), (2.9), (2.10), with the regularised definition of $C$ given in (2.28). The regularisation coefficient $f$ is defined in (2.17). By construction, any term which is multiplicated by $f$ is set equal to zero when $f=0$ . The non-monotonic friction law is defined in (2.23), and the coefficient $\alpha _3$ is equal to

(2.32) \begin{align} \alpha _3=\frac {4d}{5I_0}\sqrt {\frac {g\cos \theta }{\phi }}\frac {231(3C-4)}{4(17C+14)}(\mu _2-\mu _1), \end{align}

with $C$ defined in (2.28).

The regularised model is still a hyperbolic system of conservation laws with source terms. Since the case $C \geq 4/3$ is avoided by the regularisation, the system is always hyperbolic. It admits the energy equation

(2.33) \begin{align} \begin{aligned} &\frac {\partial he}{\partial t}+\frac {\partial }{\partial x}\left (hUe+\varPi U\right ) = h^{3/2}(1-f^\delta )\alpha _3\left (\psi -e^{-\tfrac {5d}{h}}\frac {U^2}{4h^2}\right ) + gU\cos \theta \\ &\times \left [h[\tan \theta -\mu (h,U)]+\frac {\beta f^\delta \alpha _1}{h^{1/2}}\left (U-\frac {2}{5}Ah^{3/2}\right )+\left (\beta f^\delta \alpha _2+\frac {10\lambda f}{A^2B}\right )\left (\psi -\frac {U^2}{4h^2}\right )\right ], \end{aligned} \end{align}

and is well posed for every inclination angle $\theta$ . When $\mu _1 \lt \tan \theta \lt \mu _2$ , the fact that $f=1+O(\varepsilon ^2)$ (see (2.22)) implies that the regularised model is fully consistent with the $\mu (I)$ -rheology up to order one (except in the special cases discussed hereafter). Furthermore, one can check that the linearised version of (2.29)–(2.31) around a stationary uniform flow is exactly the same as for the non-regularised system (2.4)–(2.6). As a consequence, the theoretical predictions about the neutral stability curve, the phase velocity and the spatial growth rate are the same and the study of Deleage & Richard (Reference Deleage and Richard2025) can be used for the regularised model.

When $\mu _1\lt \tan \theta \lt \mu _2$ , the consistency of the regularised system (2.29)–(2.31) with the $\mu (I)$ -rheology is lost in the two following cases:

1. (i) When the Froude number is below the threshold $\beta '$ , the friction law of the regularised model changes and transitions towards a static coefficient. This choice was used by Pouliquen & Forterre (Reference Pouliquen and Forterre2002) and enables us to capture the deposit of height $h_{\text{stop}}(\theta )$ which is left by a steady uniform flow. Note that the $\mu (I)$ -rheology predicts that steady uniform flows can have an arbitrarily small height. This is in contradiction with experimental observations, as stated by Pouliquen & Forterre (Reference Pouliquen and Forterre2002): ‘no steady uniform flow is observed with a lower Froude number, i.e. with a thickness $h$ less than $h_{\text{stop}}$ ’. The loss of consistency with the rheology is thus motivated by the fact that the predictions of the rheology are not in agreement with the observations in this regime.

2. (ii) In the special case where $\mu _2\gt \sqrt {2/3}$ and $\tan \theta$ is close to $\mu _2$ , the regularised definition of $C$ (see (2.28)) also implies that consistency with the rheology is lost. In this case, the rheology predicts that all steady uniform flows are unstable, and the non-regularised system (2.4)–(2.6) can become ill posed. Little knowledge of the behaviour of the flow in this regime is available (see, however, the work of Börzsönyi & Ecke 2006). In the absence of more information, we choose to set $C=1$ in this case (see (2.28)) in order to keep a well-posed system, thus losing consistency with the rheology. Note that the slopes at which this problem can occur are quite restricted, since $\theta$ has to be smaller than $\theta _2=\arctan (\mu _2)$ and greater than $\arctan (\sqrt {2/3})\approx 39^\circ$ . The exact lower bound for $\theta$ is slightly greater than $39^\circ$ and depends on $\mu _1$ and $\mu _2$ , but remains close to $40^\circ$ for reasonable values of $\mu _1$ and $\mu _2$ . Note that Börzsönyi & Ecke (Reference Börzsönyi and Ecke2006), with a value $\mu _2 = 1 \gt \sqrt {2/3}$ , found stationary uniform flows only under an inclination angle of $40.4^\circ$ . Few experimental data in this regime are available, and we did not compare the performance of our regularised model in this case.

As explained before, the behaviour of the regularised model close to steady uniform flows is the same as the one of the non-regularised model, which was studied by Deleage & Richard (Reference Deleage and Richard2025). In the next sections, numerical validations of the regularised model are presented in three regimes far from steady uniform flows: dry fronts, flow stopping and accelerated flows. The performances of system (2.4)–(2.6) are compared with both experimental data and the classical two-equation model (see for example Pouliquen & Forterre Reference Pouliquen and Forterre2002; Gray et al. Reference Gray, Tai and Noelle2003), given by

(2.34) \begin{align} \left \{ \begin{aligned} & \frac {\partial h}{\partial t} + \frac {\partial hU}{\partial x} = 0, \\ &\frac {\partial hU}{\partial t} + \frac {\partial }{\partial x} \left (\alpha h U^2 + \frac {gh^2}{2}\cos \theta \right )=gh\cos \theta [\tan \theta -\mu (h,U)], \end{aligned} \right . \end{align}

where $\mu (h,U)$ is defined by (2.23) and the shape factor $\alpha$ is equal to 1, except in § 3 where we also use the value $\alpha =5/4$ .