2.1. Landau–Ginzburg granular rheology model

By treating granular material as a bi-phasic system undergoing a solid–fluid phase transition, a Landau–Ginzburg equation describing the evolution of the fluidisation parameter, $\lambda$, is formulated as

(2.1)\begin{equation} t_0 \frac{{\rm D} \lambda}{{\rm D} t}=l^2 \nabla^2 \lambda +r(\mu)\lambda -B \lambda^{1+b}. \end{equation}

Here, $t_0$ represents a relaxation time, $l$ denotes a microscopic length associated with non-local transport and $B>0$ and $b>0$ are rheological parameters. The function $r(\mu )$ controls the occurrence of granular phase transitions at a critical effective frictional coefficient, $\mu _c$, requiring that $r(\mu )>0$ for $\mu >\mu _c$ and $r(\mu )<0$ for $\mu <\mu _c$. The functional form of $r(\mu )$ will be determined by comparing the model solution with empirical data. Note that, in the partially fluidised theory (Aranson & Tsimring Reference Aranson and Tsimring2002; Aranson et al. Reference Aranson, Tsimring, Malloggi and Clément2008) and the granular fluidity theory (Kamrin & Henann Reference Kamrin and Henann2015; Mowlavi & Kamrin Reference Mowlavi and Kamrin2021), the source terms are assumed a priori to present the form of an analytic expansion with integer power exponents. However, Lee & Yang (Reference Lee and Yang2017) have demonstrated that the power exponent in the dissipative term, $1+b$, is crucial for accurately reproducing the observed rheology, so here, $b$ is treated as a free parameter to be determined. Note also that an additional bistability term may be introduced in (2.1) to account for hysteresis of flow onset and arrest (Aranson & Tsimring Reference Aranson and Tsimring2002; Lee & Yang Reference Lee and Yang2017; Mowlavi & Kamrin Reference Mowlavi and Kamrin2021), which is, however, beyond the scope of this study focusing on the flow-to-arrest process.

Three potential candidates for the physical definition of $\lambda$ are discussed as follows. In the partially fluidised theory (Aranson & Tsimring Reference Aranson and Tsimring2002; Aranson et al. Reference Aranson, Tsimring, Malloggi and Clément2008), $\lambda$ conceptually denotes the concentration of granular solid phase to fluid phase, requiring an additional empirical relation to relate it to the flow variables, thereby introducing further complexity. In the granular fluidity theory (Kamrin & Koval Reference Kamrin and Koval2012; Kamrin & Henann Reference Kamrin and Henann2015; Mowlavi & Kamrin Reference Mowlavi and Kamrin2021), $\lambda$ is defined as a granular fluidity parameter, the inverse of effective viscosity rescaled by pressure, $\dot {\gamma }/\mu$. This definition introduces dimensional aspects and makes the coefficient $B$ a flow variable (Kamrin & Henann Reference Kamrin and Henann2015), posing challenges for analytical analysis. For the sake of simplicity, here $\lambda$ is chosen as the inertial number, i.e. $\lambda = I$, which is dimensionless and defined directly by the flow variables (Lee & Yang Reference Lee and Yang2017). With this choice, in a steady homogeneous shear flow, (2.1) yields a local rheology relation given by

(2.2)\begin{equation} I_{loc}(\mu)=\left( \frac{r(\mu)}{B} \right)^{1/b}. \end{equation}

2.2. Application to inclined-plane flow

Consider a steady granular flow on a rough inclined plane, where the flow height $h$ is uniform and the slope angle $\theta$ is measured from the horizon. A Cartesian coordinate is set at the bottom somewhere upstream, where $x$ points down the inclined plane and $z$ is the direction normal to the plane. Assuming the material is incompressible with a constant solid volume fraction $\phi$, the shear stress $\tau$ and the pressure $P$ can be derived from the momentum balance equations at steady state, given by

(2.3)$$\begin{gather} \tau = \rho \phi g \sin\theta (h-z), \end{gather}$$

(2.4)$$\begin{gather}P = \rho \phi g \cos\theta (h-z), \end{gather}$$

which leads to the effective frictional coefficient

(2.5)\begin{equation} \mu = \tan \theta. \end{equation}

In this configuration, the Landau–Ginzburg equation (2.1) reduces to

(2.6)\begin{equation} 0=l^2 \frac{{\rm d}^2 I}{{\rm d} z^2} +r(\theta)I -B I^{1+b}. \end{equation}

Multiplied by ${\rm d}I/{\rm d}z$, the (2.6) can be manipulated to yield

(2.7)\begin{equation} 0=\frac{{\rm d}}{{\rm d}z}\bigg[ \frac{l^2}{2} \left(\frac{{\rm d} I}{{\rm d}z}\right)^2 +\frac{r(\theta)}{2}I^2-\frac{B}{2+b} I^{2+b} \bigg], \end{equation}

or

(2.8)\begin{equation} \frac{l^2}{2}\left(\frac{{\rm d} I}{{\rm d} z}\right)^2 +\frac{r(\theta)}{2}I^2-\frac{B}{2+b} I^{2+b} =\textrm{const}. \end{equation}

To proceed with the analysis, the boundary conditions for the order parameter at the free surface and bottom are required. Since no fluidisation enters or leaves the free surface, the surface condition is typically assumed to be flux free (Aranson & Tsimring Reference Aranson and Tsimring2002; Kamrin & Henann Reference Kamrin and Henann2015; Lee & Yang Reference Lee and Yang2017)

(2.9)\begin{equation} \frac{{\rm d} I}{{\rm d} z} =0,\quad \textrm{at}\ z=h. \end{equation}

For irregular grains flowing down a rough plane considered here, the bottom can be fairly assumed to be highly dissipative, indicating a pure solid phase

(2.10)\begin{equation} I =0,\quad \textrm{at}\ z=0. \end{equation}

This bottom condition is supported by direct experimental observation of sand flows on a rough inclined plane in which a thin jammed layer forms at the bottom below the flowing layer (Aranson et al. Reference Aranson, Tsimring, Malloggi and Clément2008; Malloggi et al. Reference Malloggi, Andreotti and Clément2015). Now, with the surface condition, (2.8) reduces to

(2.11)\begin{equation} \frac{l^2}{2}\left(\frac{{\rm d} I}{{\rm d} z}\right)^2 +\frac{r(\theta)}{2}I^2-\frac{B}{2+b} I^{2+b} = \frac{r(\theta)}{2}I_s^2-\frac{B}{2+b} I_s^{2+b}, \end{equation}

where $I_s$ denotes the inertial number at the free surface. In terms of the following normalised variables:

(2.12a–c)\begin{equation} I^*=\frac{I}{I_s},\quad z^*=\frac{z}{h},\quad \epsilon = \frac{I_s}{I_{loc}(\theta)}, \end{equation}

where $I_{loc}$ is given by (2.2), (2.11) can be rewritten as

(2.13)\begin{equation} \frac{{\rm d} I^*}{{\rm d}z^*}= \frac{\rm \pi}{2} \kappa \sqrt{ 1-I^{*2} -\frac{2}{2+b}\epsilon ^{b} ( 1 - I^{*2+b})}, \end{equation}

where

(2.14)\begin{equation} \kappa = \frac{2\sqrt{B}I_{loc}(\theta)^{b/2} h}{{\rm \pi} l}. \end{equation}

With the bottom condition $I^*=0$ at $z^*=0$, (2.13) gives rise to an integral equation

(2.15)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \int^{I^*}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2} -\dfrac{2}{2+b}\epsilon ^{b} ( 1 - X^{2+b})}}. \end{equation}

Substitution of the surface condition $I^*=1$ at $z^*=1$ into (2.15) gives

(2.16)\begin{equation} \frac{\rm \pi}{2}\kappa = \int^{1}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2} -\dfrac{2}{2+b}\epsilon ^{b} ( 1 - X^{2+b})}}. \end{equation}

The above two equations suggest $I^*=I^*(z^*,\kappa,b)$ and $\epsilon =\epsilon (\kappa,b)$. With (2.14) and (2.16), the stopping height $h_{stop}$ can be determined by prescribing the flow-arrest condition $\epsilon =0$, given by

(2.17)\begin{equation} \kappa_{\epsilon=0}=\frac{2\sqrt{B}I_{loc}(\theta)^{b/2} h_{stop}}{{\rm \pi} l} =1. \end{equation}

By using (2.2) to replace $I_{loc}$ in (2.17), the stopping height is expressed in terms of $r(\theta )$ as

(2.18)\begin{equation} h_{stop}=\frac{{\rm \pi} l}{2\sqrt{r(\theta)}}. \end{equation}

As a result, the relation (2.14) can be rewritten as

(2.19)\begin{equation} \kappa = \frac{h}{h_{stop}(\theta)}, \end{equation}

which reveals that $\kappa$ is the ratio of the flow height to the stopping height.

2.3. Determination of $b$

The parameter $b$ is determined by comparing the mean velocity relation with the phenomenological observation. The mean velocity is evaluated by integrating the local velocity over the depth divided by the flow height, $\bar {u} = (1/h)\int ^{h}_0 u(z)\, {\rm d} z$. By exploiting integration by parts and assuming a no-slip condition at the bottom, $u=0$ at $z=0$,

(2.20)\begin{equation} \bar{u} = \left.\frac{1}{h} (z-h) u\right|^h_0 -\frac{1}{h} \int^{h}_0 (z-h) \frac{{\rm d}u}{{\rm d} z}\,{\rm d}z=\frac{1}{h} \int^{h}_0 (h-z) \frac{{\rm d}u}{{\rm d} z}\,{\rm d}z. \end{equation}

By replacing ${\rm d}u/{\rm d}z$ with $I$ using (1.1)

(2.21)\begin{equation} \bar{u} =\frac{\sqrt{\phi g \cos\theta }}{Dh} \int^{h}_0 I (h-z)^{3/2} \, {\rm d}z. \end{equation}

In terms of the normalised quantities (2.12a–c), the above equation can be rewritten as

(2.22)\begin{equation} Fr = \sqrt{\phi \cos\theta} \frac{h}{D} I_{loc}(\theta) \epsilon(\kappa,b) \int^{1}_0 I^*(z^*,\kappa,b) (1-z^*)^{3/2} \, {\rm d} z^*, \end{equation}

where $Fr\equiv \bar {u}/\sqrt {gh}$ is the Froude number. Note that the $\sqrt {\cos \theta }$ term can be treated as a constant since it varies slightly within the steady flow range (Jop et al. Reference Jop, Forterre and Pouliquen2005). With the help of (2.2), (2.18) and (2.19), (2.22) can be manipulated to be

(2.23)\begin{equation} Fr = \sqrt{\phi\cos\theta}\bigg(\frac{2\sqrt{B}D}{{\rm \pi} l}\bigg)^{-{2}/{b}} \left[\frac{D}{h_{stop}(\theta)}\right]^{{2}/{b}-1} \kappa \epsilon(\kappa,b) Q(\kappa,b), \end{equation}

where

(2.24)\begin{equation} Q=\int^{1}_0 I^*(z^*,\kappa,b) (1-z^*)^{3/2} \, {\rm d} z^*. \end{equation}

There has been substantial experimental and numerical evidence showing that the Froude number is one to one with the height ratio $h/h_{stop}$ over a wide range of $h/h_{stop}$ (Pouliquen Reference Pouliquen1999; Silbert et al. Reference Silbert, Landry and Grest2003; Deboeuf et al. Reference Deboeuf, Lajeunesse, Dauchot and Andreotti2006; Aranson et al. Reference Aranson, Tsimring, Malloggi and Clément2008; Malloggi et al. Reference Malloggi, Andreotti and Clément2015). To be consistent with this observation,

(2.25)\begin{equation} b=2, \end{equation}

has to be selected such that $Fr$ depends on $\kappa$ only

(2.26)\begin{equation} Fr =\alpha_0 \kappa \epsilon (\kappa ) Q(\kappa), \end{equation}

with

(2.27)\begin{equation} \alpha_0 = \frac{{\rm \pi} l}{2D} \sqrt{\frac{\phi\cos\theta}{B}}. \end{equation}

The above analysis significantly demonstrates that the choice of $b=2$, leading to the cubic term in the Landau–Ginzburg equation (2.1), is crucial to reproducing the observed mean flow characteristic. This explains why the common formulation with a quadratic dissipative term fails to yield the one-to-one $Fr-\kappa$ relation (Kamrin & Henann Reference Kamrin and Henann2015). Note that the previous study by Lee & Yang (Reference Lee and Yang2017) employed scaling analysis to derive the same result (2.25) to reproduce the Bagnold scaling in the dense inertial regime. In contrast, the mean velocity analysis adopted here is not limited to a specific flow regime, indicating that the cubic term in the Landau–Ginzburg equation is a general result valid across the quasi-static and the dense inertial regime.

2.4. Exact solutions of the internal flow profiles

With $b=2$, (2.15) and (2.16) reduce to

(2.28)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \int^{I^*}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2}} \sqrt{1-\frac{1}{2}\epsilon ^{2} ( 1 + X^{2})}}, \end{equation}

and

(2.29)\begin{equation} \frac{\rm \pi}{2}\kappa = \int^{1}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2}}\sqrt{1 -\frac{1}{2}\epsilon ^{2} ( 1 + X^{2})}}, \end{equation}

respectively, which can be further manipulated to yield

(2.30)\begin{equation} \frac{\rm \pi}{2}\sqrt{1-\frac{1}{2}\epsilon^2}\kappa = K\left(\sqrt{\frac{\epsilon^2}{2-\epsilon^2}}\right), \end{equation}

and

(2.31)\begin{equation} z^* K\left(\sqrt{\frac{\epsilon^2}{2-\epsilon^2}}\right) = F\left(I^*;\sqrt{\frac{\epsilon^2}{2-\epsilon^2}}\right). \end{equation}

Here, $K(\sqrt {{\epsilon ^2}/({2-\epsilon ^2})})$ and $F(I^*;\sqrt {{\epsilon ^2}/({2-\epsilon ^2})})$ are the complete and incomplete elliptic integrals of the first kind, respectively (Abramowitz & Stegun Reference Abramowitz and Stegun1964). Equation (2.30) indicates that $\epsilon$ is a function of $\kappa$ only. On the other hand, (2.31) describes the cross-depth profiles of the normalised inertial number $I^*$, whose shapes depend on $\epsilon$ and hence $\kappa$.

The exact solutions (2.30) and (2.31) reveal that the internal flow rheology evolves with the height ratio $\kappa$. Figure 1 plots the profile of $\epsilon$ vs $\kappa$. It shows that, as $\kappa$ increases from unity, $I_s$ monotonically increases from zero and eventually saturates to $I_{loc}$. Figure 2(a) displays the cross-depth profiles of $I/I_s$ at various $\kappa$. For a very thick flow as $\kappa \gg 1$, the profile is nearly uniform except at the very bottom where a boundary layer with large gradients of $I$ forms to conform to the solid-phase bottom. According to the $\epsilon -\kappa$ profile shown in figure 1, the values of $I$ in the uniform region approximate $I_{loc}$. This indicates that very thick flows fall into the regime dominated by local rheology. On the other hand, for a thin flow with $\kappa \gtrsim 1$, the profiles of $I$ are non-uniform across the layer. This case corresponds to the regime dominated by non-local effects, where the solid-phase bottom strongly correlates with the rheology throughout the depth. The transition between the local and the non-local features of the two regimes can be further illustrated by plotting the cross-depth profiles of $I/I_{loc}$ at a particular slope angle $\theta$ in figure 2(b). This shows that as $\kappa$ is increased, the profiles gradually collapse onto the uniform distribution $I/I_{loc}=1$, while when $\kappa$ is decreased to unity, the profile gradually vanishes. The above results provide physical insights into how the flow height controls the flow-to-no-flow transition: as the flow becomes thinner, the non-local transport capability is progressively enhanced, causing the bottom solid phase to spread and diminish fluidisation over the entire domain.

Figure 1.

Profile of $\epsilon$ vs $\kappa$ given by (2.30).

Figure 2.

Depth profiles of (a) the inertial number $I$ normalised by the surface magnitude $I_s$, (b) the inertial number $I$ rescaled by the local rheology $I_{loc}$ at $\theta =32.2^\circ$, where the functional form and the parameters in $I_{loc}$ are obtained in § 3, and (c) downslope velocity $u$ normalised by the surface magnitude; (a–c) $\kappa =1.004,1.111,1.987,3.736,5.81,13.07$. Dashed line: the Bagnold profile given by (2.33).

Note that similar transitions in the inertial number profiles have been predicted by the numerical solutions of the $I$-based Landau–Ginzburg equation with a hysteretic bistability term reported in Lee & Yang (Reference Lee and Yang2017). This suggests that the hysteresis effect does not significantly influence the non-local flow behaviour. However, the hysteresis effect can result in sudden-jump transitions between flow and arrest (Pouliquen & Forterre Reference Pouliquen and Forterre2002; Lee & Yang Reference Lee and Yang2017; Edwards et al. Reference Edwards, Russell, Johnson and Gray2019; Mowlavi & Kamrin Reference Mowlavi and Kamrin2021), which cannot be captured by the non-hysteresis exact solutions (2.30) and (2.31) which predict the transition to be smooth.

The rheological features of the locality- and the non-locality-dominated regime is also reflected in the shapes of the velocity profiles. With the no-slip bottom, the velocity is evaluated by $u=\int ^{z}_0 ({\rm d}u/{\rm d}z) \, {\rm d}\tilde {z}$. By replacing ${\rm d}u/{\rm d}z$ with $I$ and using the normalised quantities (2.12a–c), the velocity is derived to be

(2.32)\begin{equation} u=\frac{\sqrt{\phi g h^3 \cos\theta }}{D} I_{loc}(\theta) \int^{z^*}_0 I^*(\tilde{z},\kappa) (1-\tilde{z})^{1/2} \, {\rm d} \tilde{z}. \end{equation}

Figure 2(c) displays the normalised velocity profiles given by (2.30), (2.31) and (2.32) for the values of $\kappa$ corresponding to those in figure 2(a). For the large values of $\kappa$ in the locality-dominated regime, the velocity profiles collapse onto the Bagnold profile (obtained by setting $I^*=1$ in (2.32)), i.e.

(2.33)\begin{equation} u_{Bag} =\frac{2}{3}\frac{\sqrt{\phi g h^3 \cos\theta }}{D} I_{loc}(\theta) [1-(1-z^*)^{3/2} ], \end{equation}

plotted as the dashed line in figure 2(c). In the non-locality-dominated regime with $\kappa$ close to unity, the profile shapes become noticeably concave near the bottom, characterising creep-flow behaviours. The predicted shape transition is qualitatively consistent with the simulation findings of inclined-plane flows of polyhedral irregular grains (Azéma et al. Reference Azéma, Descantes, Roquet, Roux and Chevoir2012). Similar theoretical predictions of the velocity shapes have been reported in the previous Landau–Ginzburg modelling studies (Kamrin & Henann Reference Kamrin and Henann2015; Lee & Yang Reference Lee and Yang2017).

Note that, in the locality-dominated regime with a vanishing boundary layer as $\epsilon \rightarrow 1$ and $I^*\rightarrow 1$, the Pouliquen flow rule is obtained from (2.26), given by

(2.34)\begin{equation} \frac{\bar{u}}{\sqrt{gh}} =\frac{2}{5}\alpha_0 \frac{h}{h_{stop}(\theta)}. \end{equation}

This demonstrates that the flow with a highly dissipative bottom considered here can be captured by the Pouliquen flow rule when it is sufficiently thick. However, this situation takes place when $\kappa \gg 10$, according to the fact that the boundary layer thickness of $I$ is predicted to be approximately one tenth of the flow height at $\kappa \approx 10$ as shown in figure 1(b). This result agrees with the observation that the Pouliquen flow rule does not emerge within the typical experimental ranges for sandy flows, $1<\kappa <20$ (GDR-Midi 2004; Deboeuf et al. Reference Deboeuf, Lajeunesse, Dauchot and Andreotti2006; Börzsönyi & Ecke Reference Börzsönyi and Ecke2007; Aranson et al. Reference Aranson, Tsimring, Malloggi and Clément2008; Malloggi et al. Reference Malloggi, Andreotti and Clément2015), which should fall in the non-locality-dominated regime.

2.5. Asymptotic solution for the flow rule

In the non-locality-dominated regime as $\kappa \rightarrow 1$, the asymptotic solution of the mean flow rule (2.26) is sought as follows. In this limit, since $\epsilon ^2 \rightarrow 0$ as $\kappa \rightarrow 1$ according to (2.30), (2.28) can be asymptotically expanded in terms of $\epsilon ^2$, resulting in

(2.35)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \int^{I^*}_0 \frac{{\rm d}X}{\sqrt{1-X^{2}}}+\frac{\epsilon^2}{4} \int^{I^*}_0 \frac{1+X^2}{\sqrt{ 1-X^{2}}}\, {\rm d} X +O(\epsilon^4). \end{equation}

After evaluating the integrals, the equation becomes

(2.36)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \sin^{-1}(I^*)+\frac{\epsilon^2}{8}[3 \sin^{-1}(I^*)-I^*\sqrt{1-I^{*2}}] +O(\epsilon^4). \end{equation}

By inserting the free-surface condition $I^*=1$ at $z^*=1$ into (2.36), an asymptotic relation for $\kappa (\epsilon )$ is then obtained, given by

(2.37)\begin{equation} \kappa -1 = \tfrac{3}{8} \epsilon^2 +O(\epsilon^4), \end{equation}

which indicates an asymptotic relation for $\epsilon$ in terms of $\kappa -1$

(2.38)\begin{equation} \epsilon^2 = \tfrac{8}{3} (\kappa -1) + O[(\kappa-1)^2]. \end{equation}

To evaluate the integral relation (2.24) for $Q$, the explicit solution of $I^*$ is derived as follows. First, (2.36) is rewritten as

(2.39)\begin{equation} I^*=\sin\left( \frac{4 {\rm \pi}}{8+3\epsilon^2} \kappa z^* + \frac{\epsilon^2}{8+3\epsilon^2} I^*\sqrt{1-I^{*2}}+O(\epsilon^4) \right). \end{equation}

By replacing $\epsilon ^2$ with $\kappa -1$ given by the relation (2.38), the above equation reduces to

(2.40)\begin{equation} I^*=\sin\left( \frac{\rm \pi}{2} z^* + \frac{1}{3} \frac{\kappa-1}{\kappa }I^*\sqrt{1-I^{*2}}+ O[(\kappa-1)^2]\right), \end{equation}

which can be expanded in terms of $\kappa -1$ to yield

(2.41)\begin{equation} I^*= \sin\left( \frac{\rm \pi}{2} z^*\right) + \frac{1}{3} (\kappa-1) \cos\left(\frac{\rm \pi}{2}z^*\right) I^*\sqrt{1-I^{*2}}+ O[(\kappa-1)^2]. \end{equation}

According to the leading-order behaviour $I^*\sim \sin ({\rm \pi} z^*/2)$, (2.41) can then be asymptotically reduced to

(2.42)\begin{equation} I^*=\sin\left( \frac{\rm \pi}{2} z^*\right) + \frac{1}{3} (\kappa-1) \sin\left( \frac{\rm \pi}{2} z^*\right) \cos^2\left(\frac{\rm \pi}{2}z^*\right)+ O[(\kappa-1)^2]. \end{equation}

Substitution of (2.42) into (2.24) gives

(2.43)\begin{align} Q &= \int^{1}_0 \sin\left( \frac{\rm \pi}{2} z^*\right) (1-z^*)^{3/2} \, {\rm d} z^* \nonumber\\ & \quad + \frac{1}{3} (\kappa-1) \int^{1}_0 \sin\left( \frac{\rm \pi}{2} z^*\right) \cos^2\left(\frac{\rm \pi}{2}z^*\right) (1-z^*)^{3/2} \, {\rm d} z^*+O[(\kappa-1)^2]. \end{align}

After evaluating the integrals, $Q$ is derived to be

(2.44)\begin{equation} Q \approx 0.1625+0.03227 (\kappa -1) + O[(\kappa-1)^2]. \end{equation}

Substitution of (2.38) and (2.44) into (2.26) yields an asymptotic relation for $Fr$ in terms of $\kappa -1$, given by

(2.45)\begin{equation} Fr = \alpha_0 C_0 (\kappa -1)^{1/2} [1+ C_1 (\kappa -1) +O[(\kappa-1)^2] ], \end{equation}

where $C_0=0.265$ and $C_1=0.407$. Finally, the analytical solution of the flow rule (2.45) to the $(\kappa -1)^{3/2}$-order can be explicitly expressed in terms of $u/\sqrt {gh}$ and $h/h_{stop}$ as

(2.46)\begin{equation} \frac{\bar{u}}{\sqrt{gh}} = \alpha_0 C_0 \left(\frac{h}{h_{stop}(\theta)} -1\right)^{1/2} \left[1+ C_1 \left(\frac{h}{h_{stop}(\theta)} -1\right) \right]. \end{equation}

Note that $\alpha _0$ is to be determined by fitting of the solution (2.46) to empirical data.