# Surface curvature and secondary vortices in steady
dense shallow granular flows

Cyril Gadal, Chris G Johnson, J.M.N.T. Gray

► To cite this version:

Cyril Gadal, Chris G Johnson, J.M.N.T. Gray. Surface curvature and secondary vortices in steady dense shallow granular flows. Journal of Fluid Mechanics, 2026, 1027, pp.A34. 10.1017/jfm.2025.11043. hal-05472077

HAL Id: hal-05472077

https://hal.science/hal-05472077v1

Submitted on 22 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
<figure><img src="image_2.png" /></figure>

# Surface curvature and secondary vortices in steady dense shallow granular flows

Cyril Gadal¹, Chris G. Johnson¹ and J.M.N.T. Gray¹

¹Department of Mathematics and Manchester Centre for Nonlinear Dynamics, University of Manchester, Oxford Road, Manchester M13 9PL, UK

**Corresponding author:** Cyril Gadal, cyril.gadal@univ-rennes.fr

(Received 12 June 2025; revised 23 October 2025; accepted 8 December 2025)

Dense granular flows exhibit both surface deformation and secondary flows due to the presence of normal stress differences. Yet, a complete mathematical modelling of these two features is still lacking. This paper focuses on a steady shallow dense flow down an inclined channel of arbitrary cross-section, for which asymptotic solutions are derived by using an expansion based on the flow's spanwise shallowness combined with a second-order granular rheology. The leading-order flow is uniaxial with a constant inertial number fixed by the inclination angle. The streamwise velocity then corresponds to a lateral juxtaposition of Bagnold profiles scaled by the varying flow depth. The correction at first order introduces two counter-rotating vortices in the plane perpendicular to the main flow direction (with downwelling in the centre), and an upward curve of the free surface. These solutions are compared with discrete element method simulations, which they match quantitatively. This result is then used together with laboratory experiments to infer measurements of the second-normal stress difference in dense dry granular flow.

**Key words:** granular media, non-Newtonian flows

## 1. Introduction

When set in motion, granular materials exhibit a range of intriguing behaviours that arise from complex interactions at the grain scale. One peculiar phenomenon is the upward curvature of the free surface of dense granular flows in certain configurations, such as open channels with frictional walls (Félix & Thomas 2004; Deboeuf et al. 2006; Takagi, McElwaine & Huppert 2011; McElwaine, Takagi & Huppert 2012), non-rectangular channel geometries (see experiments in figure 1 and § 6) and Couette cells (Dsouza & Nott 2021; Cabrera & Polanía 2020). Grain-scale numerical simulations indicate that
this surface deformation is accompanied by the existence of secondary counter-rotating vortices in the plane perpendicular to the main flow direction (Krishnaraj & Nott 2016; Dsouza & Nott 2021; Kim & Kamrin 2023). These secondary flows have also been observed recently in experiments involving slow bulldozing of granular material using a conveyor belt (Einav et al. 2024). Interestingly, both surface curvature and vortices are enhanced in flows of elongated particles, supporting the idea that these effects stem from mechanisms occurring at the grain scale (Wortel et al. 2015; Fischer et al. 2016; Mohammadi et al. 2022).

In dense granular flows, the surface curvature and secondary flows are generally attributed to a non-zero, negative, second-normal stress difference (between the gradient and the vorticity direction) resulting from an anisotropic distribution of the contacts at the grain scale (Nagy et al. 2020; Srivastava et al. 2021). This is more generally related to the broken codirectionality and coaxiality between the stress and strain rate tensors, as observed in grain-scale simulations (e.g. Campbell 1989; Silbert et al. 2001; Weinhart et al. 2013; Srivastava et al. 2021). Note, however, that both surface deformation and secondary vortices have also been reported in rapid chute flows. In that regime, they are usually attributed to fluctuations of the granular temperature, which also cause significant spatial variations in solids volume fraction (Forterre & Pouliquen 2001; Börzsönyi et al. 2009; Brodu et al. 2015).

Dense inertial granular flows are commonly described using the so-called $\mu(I)$ granular local rheology (GDR MiDi 2004; Jop, Forterre & Pouliquen 2006). In this framework, the effective viscosity of the granular fluid depends on the local pressure due to the friction between the grains, but also on a dimensionless quantity, the inertial number $I$, which compares the time scale of macroscopic deformation with that of grain-scale rearrangements (GDR MiDi 2004). This constitutive law has been applied to a variety of flow configurations, including inclined channels, silo discharges and column collapses (e.g. Jop et al. 2006; Lagrée et al. 2011; Staron, Lagrée & Popinet 2014; Martin et al. 2017). However, the original $\mu(I)$ rheology has several limitations. First, it is mathematically ill-posed at both low and high values of $I$ when combined with the canonical shape of the $\mu(I)$ curve (Barker et al. 2015). Partial regularisation, using a specific shape of the $\mu(I)$ curve, is then required to obtain grid-resolution convergent simulations below a critical inertial number (Barker et al. 2017, 2021; Maguire et al. 2024; Lloyd et al. 2025). Second, the original form of this rheology cannot account for certain phenomena observed in experiments or grain-scale simulations. For instance, some flows can dilate and contract, motivating the development of compressible extensions in which the volume fraction $\Phi$ also depends on the inertial number, in addition to the friction coefficient $\mu$ (GDR MiDi 2004; Pouliquen et al. 2006; Forterre & Pouliquen 2008). However, such compressible $\mu(I)$ models are generally ill-posed unless formulated under specific frameworks (Barker et al. 2017; Heyman et al. 2017; Schaeffer et al. 2019). Moreover, some flows exhibit non-local effects responsible for the slow creeping motions observed in heap flows far from the free surface (Komatsu et al. 2001) or coexisting adjacent static and flowing regions (e.g. Edwards & Gray 2015; Rocha, Johnson & Gray 2019), which lead to corresponding extensions of the $\mu(I)$ rheology (Pouliquen & Forterre 2009; Bouzid et al. 2015; Kamrin 2019).

To account for broken codirectionality, broken coaxiality and the existence of normal stress differences, the $\mu(I)$ rheology can also be extended to a second-order formulation, similar to the general form of a second-order fluid, also known as the Rivlin–Ericksen fluid of second grade (Rivlin & Ericksen 1997; McElwaine et al. 2012; Srivastava et al. 2021). Recently, this second-order rheology has also been extended to incorporate non-locality (Kim & Kamrin 2023). Using the local version of this second-order rheology,
<figure><img src="image_3.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 1. Dense flow of glass beads of diameter $d \in [125, 165] \mu m$ down an inclined channel with a parabolic cross-section in a laboratory experiment. (a–d) Snapshots of an experiment, featuring the empty channel, the propagating flow front and the steady uniform flow (see also movie 1 of the Supplementary material). The channel is coated with a rough brown sandpaper, and the laser line shows the topography. (e) Topography laser measurements for an inclination of 26° and a flux of 51 g s⁻¹. The cross-section of the empty channel, corresponding to (a), is shown by the brown line. The surface of the steady-state granular flow, corresponding to (d), is shown by the red line and is curved upward compared with the flat black dashed line. See § 6 for additional details on the experiments.</figcaption></figure>
C. Gadal, C.G. Johnson and J.M.N.T. Gray

McElwaine et al. (2012) linked the upward surface curvature to the cross-stream gradient of the streamwise velocity, induced by friction on static levees present on the sides. However, the lack of a closure describing the flow behaviour at the interface between flowing and static regions prevented a complete resolution of the problem. Recently, Kim & Kamrin (2023) numerically solved the continuum mass and momentum equations using this second-order rheology and successfully reproduced the secondary flows observed in fully discrete, grain-scale simulations, in which the momentum balance is solved for each particle using the discrete element method (DEM). However, rather than solving the complete free-boundary problem, they imposed spatial boundary conditions (including the free surface shape) directly extracted from the DEM simulations. Hence, a complete model of the entire flow structure in open channels remains elusive.

Surface deformation and secondary flows are also observed in other types of non-Newtonian fluids, such as viscoelastic fluids or dense suspensions (e.g. Maklad & Poole 2021). Dense suspensions, in particular, share several features with dense granular flows due to the frictional interactions between the particles forming the suspension (Guazzelli & Pouliquen 2018). In their case, the second-normal stress difference has been measured using standard rheometer geometries (Singh & Nott 2003; Dbouk, Lobry & Lemaire 2013), but also using an alternative approach based on the measurement of the surface deformation in Weissenberg (rotating-rod) and tilted-trough geometries, in which confinement effects can be reduced and sensitivity improved (Boyer, Pouliquen & Guazzelli 2011; Couturier et al. 2011; Dai et al. 2013). However, these approaches rely on a known relation between the normal stress difference and the surface deformation, which is currently lacking in the case of dense granular flows.

This paper starts by extending the approach of McElwaine et al. (2012) by considering a steady granular flow inside a tilted channel of general (but smooth) cross-section, and derives asymptotic solutions for the flow structure and surface deformation. These solutions are then compared with the results of DEM simulations, which they match quantitatively. Finally, laboratory experiments are performed, in which simple measurements of the surface deformation and the channel shape are used to infer measurements of the second-normal stress differences using the results of the mathematical model.

## 2. Governing equations and second-order granular rheology

### 2.1. System and governing equations

The system considered here is sketched in figure 2. It consists of a dense flow of particles, of diameter $d$ and intrinsic density $\rho_p$, down an inclined channel of arbitrary cross-section. Let $Oxyz$ be a Cartesian coordinate system with the $x$-axis oriented down the slope, at an angle $\theta$ to the horizontal, and the $y$-axis and $z$-axis oriented in the cross-slope and upward normal directions, respectively. The flow has velocity components in the three directions, i.e. $\mathbf{u} = (u, v, w)$ and a bulk density $\rho = \Phi\rho_p$, where $\Phi$ is the solids volume fraction.

The dynamics is then governed by the mass conservation equation,

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (2.1) $$

and the momentum balance equation,

$$ \rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = \nabla \cdot \boldsymbol{\sigma} + \rho g, \qquad (2.2) $$

where $\nabla$ is the gradient operator, ‘.’ the dot product, $\boldsymbol{\sigma}$ the Cauchy stress tensor and $g = g(\sin \theta, 0, -\cos \theta)$ the gravitational acceleration.
<figure><img src="image_8.png" /><figcaption>Figure 2. Conceptual sketch corresponding to the laboratory experiment shown in figure 1 and movie 1 of the Supplementary material, showing the base shape, surface curvature and flow structure. The velocity $\mathbf{u} = (u, v, w)$ has a main component in the downslope direction $x$, and secondary flows $(v, w)$ in the cross-slope and normal directions. Note that the aspect ratio has been increased compared with the shallow flows considered here for the sake of clarity.</figcaption></figure>

## 2.2. The second-order granular rheology

To the authors' knowledge, the only granular constitutive model to include the non-coaxiality and non-codirectionality between principal directions of stress and strain rate tensors, and therefore also accounting for both the first- and second-normal stress differences, is constructed as a series of Rivlin-Ericksen tensors (Rivlin & Ericksen 1997; McElwaine et al. 2012; Wu et al. 2017; Srivastava et al. 2021; Kim & Kamrin 2023). Note that this rheology neglects any elastic effects and assumes that non-hydrostatic stresses emerge only from the flow and are then zero when the material is at rest (Srivastava et al. 2021). The strain rate and spin rate tensors, $\mathbf{D}$ and $\mathbf{W}$, are defined from the velocity gradient $\nabla \mathbf{u}$ as

$$ \mathbf{D} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T), \quad \mathbf{W} = \frac{1}{2}(\nabla \mathbf{u} - \nabla \mathbf{u}^T). \qquad (2.3) $$

The second-order Cauchy stress tensor is then

$$ \boldsymbol{\sigma} = p \left( -\boldsymbol{I} + \frac{\mu_1}{\|\boldsymbol{D}\|} \left[ \boldsymbol{D} - \frac{1}{3} \mathrm{tr}(\boldsymbol{D}) \boldsymbol{I} \right] - \frac{\mu_2}{\|\boldsymbol{D}\|^2} \left[ \boldsymbol{D}^2 - \frac{1}{3} \mathrm{tr}(\boldsymbol{D}^2) \boldsymbol{I} \right] - \frac{\mu_3}{\|\boldsymbol{D}\|^2} \mathring{\boldsymbol{D}} \right), \quad (2.4) $$

where $p = -\mathrm{tr}(\boldsymbol{\sigma})/3$ is the pressure, since all but the unit tensor $\boldsymbol{I}$ are deviatoric, $\|\boldsymbol{D}\|$ is the second invariant of the strain-rate tensor,

$$ \|\boldsymbol{D}\| = \sqrt{\frac{1}{2} \mathrm{tr}(\boldsymbol{D}^2)}, \qquad (2.5) $$

and $\mathring{\boldsymbol{D}}$ is the Jaumann derivative,

$$ \mathring{\boldsymbol{D}} = \mathbf{D} - \boldsymbol{W}\boldsymbol{D} + \boldsymbol{D}\boldsymbol{W}, \qquad (2.6) $$
where $\dot{D} = \partial D/\partial t + (u \cdot \nabla)D$ is the material derivative. In principle, the rheological coefficients $\mu_1$, $\mu_2$ and $\mu_3$ may depend on any dimensionless parameters of the system, such as the volume fraction, the dimensionless granular temperature (Kim & Kamrin 2020, 2023), or the inertial number,

$$I = \frac{\dot{\gamma}d}{\sqrt{p/\rho_p}}, \qquad (2.7)$$

corresponding to a ratio of time scales at the particle scale, $d/\sqrt{p/\rho_p}$, and at the flow scale, $1/\dot{\gamma}$. Here, the definition of the shear rate is

$$\dot{\gamma} = 2 \|D\|. \qquad (2.8)$$

In practice, the rheological coefficients have been shown using DEM simulations to predominantly depend on the inertial number (GDR MiDi 2004; Srivastava et al. 2021). Note that this is not the case for flows very close to rest (typically $I < 10^{-2}$), where these coefficients also depend on the granular temperature/fluidity as non-locality needs to be taken into account (Kim & Kamrin 2020, 2023). However, the inertial number is assumed to be large enough for these effects to be negligible.

Although various functional forms can be used to represent $\mu_1(I)$, $\mu_2(I)$ and $\mu_3(I)$, the following derivations can be made without prescribing any specific form. This is because the inertial number is constant at leading order (see § 4.1) and the first-order cross-stream and normal velocities are not coupled to the first-order corrections to the inertial number (see § 4.2). In the classical $\mu(I)$ rheology (Jop et al. 2006), which is recovered by neglecting the second-order terms and compressibility ($\text{tr}(D) = 0$), the Cauchy stress,

$$\sigma = p \left[ -I + \mu_1(I) \frac{D}{\|D\|} \right], \qquad (2.9)$$

and the form of $\mu_1(I)$ determines the region of $I$ over which the rheology is well-posed (Barker et al. 2015, 2017). The form of $\mu_2(I)$ and $\mu_3(I)$ will influence the well-posedness of the second-order rheology (3.3), but whether the resulting equations are well-posed or ill-posed is as yet unknown. Analysis of well-posedness would be essential for studying time-dependent solutions of (3.3), but this lies beyond the scope of this study, which focuses on steady-state solutions, where ill-posedness is not an issue.

The loss of coaxiality between stress and strain-rate included in the constitutive model (2.4) corresponds to the existence of normal stress differences: the diagonal components of the stress are no longer equal (in contrast to the first-order $\mu(I)$ rheology (2.9)). These effects are usually described by the canonical first- and second-normal stress differences, which are defined relative to the local coordinate system formed by the flow: flow direction, the velocity-gradient direction and the vorticity direction. The first normal stress difference is the stress difference between the flow and gradient directions, while the second is the difference between the gradient and vorticity directions. The local flow coordinate system generally depends on space and time; only under additional kinematic assumptions (see § 4) can they be identified with the Cartesian coordinate system Oxyz.

### 2.3. Boundary conditions

The kinematic and dynamic boundary conditions at the flow surface $s$ are

$$\frac{\partial s}{\partial t} + n \cdot u = 0, \quad z = s, \qquad (2.10)$$
$$ \sigma \cdot n = 0, \quad z = s, \qquad (2.11) $$

where **n** is a unit vector normal to the flow surface. Additionally, there is no slip at the base **b**, such that

$$ \mathbf{u} = 0, \quad z = b. \qquad (2.12) $$

### 3. Assumptions and reduced problem

#### 3.1. Governing equations

Motivated by the experiments shown in figure 1 and movie 1 of the Supplementary material, the system is assumed to be steady in time and invariant along the x-axis, such that the flow velocity only depends on y and z while the flow height h and channel cross-section b only depend on y. The solids volume fraction $\Phi$ and the bulk flow density ($\rho = \Phi\rho_p$) are also assumed to be constant and uniform (GDR MiDi 2004). Under those assumptions, the mass and momentum equations reduce to

$$ \nabla \cdot \mathbf{u} = 0, \qquad (3.1) $$

$$ \rho(\mathbf{u} \cdot \nabla)\mathbf{u} = \nabla \cdot \boldsymbol{\sigma} + \rho g. \qquad (3.2) $$

The Cauchy stress tensor also simplifies to

$$ \boldsymbol{\sigma} = p \left( -\mathbf{I} + \frac{\mu_1}{\|\boldsymbol{D}\|} \boldsymbol{D} - \frac{\mu_2}{\|\boldsymbol{D}\|^2} \left[ \boldsymbol{D}^2 - \frac{1}{3} \mathrm{tr}(\boldsymbol{D}^2) \mathbf{I} \right] - \frac{\mu_3}{\|\boldsymbol{D}\|^2} \left[ (\mathbf{u} \cdot \nabla)\boldsymbol{D} - \mathbf{W}\boldsymbol{D} + \mathbf{D}\mathbf{W} \right] \right), \qquad (3.3) $$

as mass conservation implies that $\mathrm{tr}(\boldsymbol{D}) = 0$.

#### 3.2. Boundary conditions

Under those assumptions, the boundary conditions at the flow surface (2.10) and (2.11) reduce to

$$ v s_y = w, \quad z = s, \qquad (3.4) $$

$$ \sigma^{xy} s_y = \sigma^{xz}, \quad z = s, \qquad (3.5) $$

$$ \sigma^{yy} s_y = \sigma^{yz}, \quad z = s, \qquad (3.6) $$

$$ \sigma^{zy} s_y = \sigma^{zz}, \quad z = s. \qquad (3.7) $$

The subscript notation is used for derivatives in this paper, such that $\partial X/\partial y = X_y$ for the y coordinate, and similarly for the z coordinate. Additionally, the flow depth is zero at the lateral edges,

$$ s = b, \quad y = \pm \frac{W}{2}, \qquad (3.8) $$

where W is the cross-stream horizontal extent of the flow. Finally, the conservation of the mass flux,

$$ Q = \rho \int_{-W/2}^{W/2} \int_{b(y)}^{s(y)} u(y, z) \, dz \, dy, \qquad (3.9) $$

will be used to make comparisons between the mathematical model, the experiments and the numerical simulations.

#### 3.3. Scalings and dimensionless equations

Neglecting the normal stress differences, properties of downslope steady-state granular flows result from a balance between friction and gravity (e.g. GDR MiDi 2004),
such that

$$ \mu_1 \mathcal{P} \sim \rho g \mathcal{H} \sin \theta, \qquad (3.10) $$

where $\mathcal{H}$ and $\mathcal{P}$ are the characteristic height and pressure scales. As the pressure is hydrostatic,

$$ \mathcal{P} \sim \rho g \mathcal{H} \cos \theta, \qquad (3.11) $$

this leads to $\mu_1(\mathcal{I}) \sim \tan \theta$, or equivalently to a characteristic inertial number,

$$ \mathcal{I} = \mu_1^{-1} (\tan \theta). \qquad (3.12) $$

The inertial number is then controlled by the slope, but also the friction law of the considered granular media.

As the shear rate is dominated by the vertical velocity gradient $\mathcal{U}/\mathcal{H}$, (2.7) implies a Bagnold scaling for the velocity

$$ \mathcal{U} \sim \mathcal{I} \frac{\mathcal{H}^{3/2}}{d} \sqrt{\Phi g \cos \theta}. \qquad (3.13) $$

The characteristic vertical and horizontal length scales, $\mathcal{H}$ and $\mathcal{W}$, can be determined by rewriting dimensionally the edge condition (3.8) and the flux conservation (3.9) as

$$ \mathcal{H} = b \left( \frac{W}{2} \right), \qquad (3.14) $$

$$ Q \sim \rho \mathcal{U} \mathcal{H} \mathcal{W}, \qquad (3.15) $$

which form, together with (3.13), a system of three equations with three unknowns, $\mathcal{U}$, $\mathcal{H}$ and $\mathcal{W}$.

The characteristic scales defined by (3.11), (3.12), (3.13), (3.14) and (3.15) can then be computed by prescribing values for $d$, $\rho_p$, $\Phi$, $g$, $\theta$, $Q$, the base shape $b(y)$ and a friction law $\mu_1(I)$. Note that, instead of a friction law, it is also possible to specify a value of $\mathcal{I}$ and compute $\mu_1$ from $\tan \theta$.

These scales suggest introducing dimensionless variables, indicated by the tilde, as

$$ (z, s, y, W) = \mathcal{H} \left( \tilde{z}, \tilde{s}, \frac{\tilde{y}}{\varepsilon}, \frac{\tilde{W}}{\varepsilon} \right), \qquad (3.16) $$

$$ \mathbf{u} = \mathcal{U} \tilde{\mathbf{u}}, \qquad (3.17) $$

$$ (p, \sigma) = \mathcal{P}(\tilde{p}, \tilde{\sigma}), \qquad (3.18) $$

$$ I = \mathcal{I} \tilde{I}, \qquad (3.19) $$

where $\varepsilon = \mathcal{H}/\mathcal{W}$ accounts for the spanwise shallowness of the flow. Note that the inertial number $I$, although already dimensionless, is still scaled so that all leading-order quantities are of order $O(1)$. The dimensionless momentum and mass balance equations then become

$$ \varepsilon \tilde{v}_{\tilde{y}} + \tilde{w}_{\tilde{z}} = 0, \qquad (3.20) $$

$$ Fr^2(\varepsilon \tilde{v} \tilde{u}_{\tilde{y}} + \tilde{w} \tilde{u}_{\tilde{z}}) = \tilde{\sigma}_{\tilde{y}}^{xy} \varepsilon + \tilde{\sigma}_{\tilde{z}}^{xz} + \tan(\theta), \qquad (3.21) $$

$$ Fr^2(\varepsilon \tilde{v} \tilde{v}_{\tilde{y}} + \tilde{w} \tilde{v}_{\tilde{z}}) = \tilde{\sigma}_{\tilde{y}}^{yy} \varepsilon + \tilde{\sigma}_{\tilde{z}}^{yz}, \qquad (3.22) $$

$$ Fr^2(\varepsilon \tilde{v} \tilde{w}_{\tilde{y}} + \tilde{w} \tilde{w}_{\tilde{z}}) = \tilde{\sigma}_{\tilde{y}}^{zy} \varepsilon + \tilde{\sigma}_{\tilde{z}}^{zz} - 1, \qquad (3.23) $$
where $Fr = U/\sqrt{g\mathcal{H}\cos\theta} = \mathcal{I}\mathcal{W}\sqrt{\Phi}\epsilon/d$ is the Froude number. The dimensionless velocity gradient is

$$ \tilde{\nabla}\tilde{u} = \begin{bmatrix} 0 & \epsilon\tilde{u}_{\tilde{y}} & \tilde{u}_{\tilde{z}} \\ 0 & \epsilon\tilde{v}_{\tilde{y}} & \tilde{v}_{\tilde{z}} \\ 0 & \epsilon\tilde{w}_{\tilde{y}} & \tilde{w}_{\tilde{z}} \end{bmatrix}, \qquad (3.24) $$

from which all subsequent quantities (e.g. $\tilde{D}$, $\tilde{\sigma}$) can be computed. The dimensionless boundary conditions are

$$ \epsilon\tilde{v}\tilde{s}_{\tilde{y}} = \tilde{w}, \quad \tilde{z} = \tilde{s}, \qquad (3.25) $$

$$ \epsilon\tilde{\sigma}^{xy}\tilde{s}_{\tilde{y}} = \tilde{\sigma}^{xz}, \quad \tilde{z} = \tilde{s}, \qquad (3.26) $$

$$ \epsilon\tilde{\sigma}^{yy}\tilde{s}_{\tilde{y}} = \tilde{\sigma}^{yz}, \quad \tilde{z} = \tilde{s}, \qquad (3.27) $$

$$ \epsilon\tilde{\sigma}^{zy}\tilde{s}_{\tilde{y}} = \tilde{\sigma}^{zz}, \quad \tilde{z} = \tilde{s}, \qquad (3.28) $$

and the conservation of the total flux (3.9) becomes

$$ \int_{-\tilde{W}/2}^{\tilde{W}/2} \int_{\tilde{b}(\tilde{y})}^{\tilde{s}(\tilde{y})} \tilde{u}(\tilde{y}, \tilde{z}) \, d\tilde{z} \, d\tilde{y} = 1. \qquad (3.29) $$

Finally, the inertial number is

$$ \tilde{I} = \frac{\tilde{\gamma}}{\sqrt{\tilde{p}}}. \qquad (3.30) $$

### 4. Asymptotic solution for a transverse shallow flow

Here, the flow is assumed to be transverse shallow, i.e. to have typical depth $\mathcal{H}$ much smaller than its width $\mathcal{W}$, such that $\epsilon = \mathcal{H}/\mathcal{W} \ll 1$. First, under the shallow-flow assumption, the local flow basis can be taken as fixed: the dominant flow direction is $x$, the dominant gradient direction is $z$, and the dominant vorticity direction is then $y$. In this case, the scaled first normal stress difference (between flow and gradient direction) is

$$ N_1 = \frac{\sigma^{xx} - \sigma^{zz}}{p} = \frac{\tilde{\sigma}^{xx} - \tilde{\sigma}^{zz}}{\tilde{p}}, \qquad (4.1) $$

and the scaled second-normal stress difference (between gradient and vorticity direction) is

$$ N_2 = \frac{\sigma^{zz} - \sigma^{yy}}{p} = \frac{\tilde{\sigma}^{zz} - \tilde{\sigma}^{yy}}{\tilde{p}}. \qquad (4.2) $$

Note that here, the normal stress differences are scaled by the local pressure to be consistent with the previous literature (e.g. Nagy et al. 2020; Srivastava et al. 2021; Kim & Kamrin 2023).

Second, the transverse shallowness of the flow is exploited to perform an asymptotic expansion of the system of equations. The flow is then assumed to be described by a streamwise base flow plus a small correction as

$$ \tilde{u} = \tilde{u}^{(0)} e_x + \epsilon \begin{bmatrix} \tilde{u}^{(1)} \\ \tilde{v}^{(1)} \\ \tilde{w}^{(1)} \end{bmatrix} + \epsilon^2 \begin{bmatrix} \tilde{u}^{(2)} \\ \tilde{v}^{(2)} \\ \tilde{w}^{(2)} \end{bmatrix} + O(\epsilon^3), \qquad (4.3) $$

where the superscripts (i) denote quantities of order $i$, respectively. The leading order of the mass balance (3.20) implies directly that $\tilde{w}^{(1)} = 0$. Similarly, other quantities can be
expanded, including the rheological coefficients, as

$$ (\tilde{\sigma}, \tilde{p}, \tilde{I}) = (\tilde{\sigma}^{(0)}, \tilde{p}^{(0)}, \tilde{I}^{(0)}) + \varepsilon(\tilde{\sigma}^{(1)}, \tilde{p}^{(1)}, \tilde{I}^{(1)}) + O(\varepsilon^2), \quad (4.4) $$

$$ (\tilde{s}, \tilde{W}) = (\tilde{s}^{(0)}, \tilde{W}^{(0)}) + \varepsilon(\tilde{s}^{(1)}, \tilde{W}^{(1)}) + O(\varepsilon^2), \quad (4.5) $$

$$ (\mu_1, \mu_2, \mu_3) = (\mu_1^{(0)}, \mu_2^{(0)}, \mu_3^{(0)}) + \varepsilon(\mu_1^{(1)}, \mu_2^{(1)}, \mu_3^{(1)}) + O(\varepsilon^2), \quad (4.6) $$

$$ (N_1, N_2) = (N_1^{(0)}, N_2^{(0)}) + \varepsilon(N_1^{(1)}, N_2^{(1)}) + O(\varepsilon^2). \quad (4.7) $$

These are only expanded up to the linear order, as higher orders will not be involved in the expressions derived later.

### 4.1. Base state velocity field

Solutions for the base state are derived by keeping only the terms of order zero. The shear stress tensor reduces to

$$ \tilde{\sigma}^{(0)} = \begin{bmatrix} -A\tilde{p}^{(0)} & 0 & \mu_1^{(0)}\tilde{p}^{(0)} \\ 0 & -B\tilde{p}^{(0)} & 0 \\ \mu_1^{(0)}\tilde{p}^{(0)} & 0 & -C\tilde{p}^{(0)} \end{bmatrix}, \quad (4.8) $$

where

$$ A = 1 - \frac{1}{3}(6\mu_3^{(0)} - \mu_2^{(0)}), \quad (4.9) $$

$$ B = 1 - \frac{1}{3}(2\mu_2^{(0)}), \quad (4.10) $$

$$ C = 1 + \frac{1}{3}(6\mu_3^{(0)} + \mu_2^{(0)}). \quad (4.11) $$

At the leading order, the scaled first- and second-normal stress differences are then

$$ N_1^{(0)} = \frac{\tilde{\sigma}^{xx}(0) - \tilde{\sigma}^{zz}(0)}{\tilde{p}^{(0)}} = C - A = 4\mu_3^{(0)}, \quad (4.12) $$

$$ N_2^{(0)} = \frac{\tilde{\sigma}^{zz}(0) - \tilde{\sigma}^{yy}(0)}{\tilde{p}^{(0)}} = B - C = -(\mu_2^{(0)} + 2\mu_3^{(0)}). \quad (4.13) $$

The mass balance equation and the y momentum equation are trivially satisfied. The remaining x and z momentum equations then reduce to

$$ \mu_1^{(0)} \tilde{p}_z^{(0)} = -\tan \theta, \quad (4.14) $$

$$ -C\tilde{p}_z^{(0)} = 1. \quad (4.15) $$

As $\mu_1^{(0)}, \mu_2^{(0)}$ and $\mu_3^{(0)}$ are functions of $\tilde{I}^{(0)}$, taking the ratio of (4.14) and (4.15) leads to an implicit equation for $\tilde{I}^{(0)}$,

$$ \frac{\mu_1^{(0)}}{C} = \mu_1^{(0)} \left[ 1 + \frac{1}{3}(6\mu_3^{(0)} + \mu_2^{(0)}) \right]^{-1} = \tan \theta. \quad (4.16) $$

This shows that, at leading order, the inertial number is controlled by the slope $\theta$ and is constant across the whole flow, which is typical for dense granular flows on inclined planes (e.g. GDR MiDi 2004).

The shear rate, which simplifies here to $\tilde{u}_z^{(0)}$, is then selected from its relationship with the pressure (3.30) as

$$ \tilde{u}_z^{(0)} = \tilde{I}^{(0)} \sqrt{\tilde{p}^{(0)}}. \quad (4.17) $$
<figure><img src="image_9.png" /><figcaption>Figure 3. Base streamwise velocity $\tilde{u}^{(0)}$ computed from (4.19) in the case of a parabolic channel. (a) Velocity in colourscale in the plane $(\tilde{y}, \tilde{z})$. The red and brown lines represent the flow surface and channel base, respectively. (b,c) Horizontal and normal velocity profiles, as annotated by the solid and dashed lines in (a). In (b), the red curve shows the streamwise velocity along the surface, $\tilde{u}^{(0)}(\tilde{y}, \tilde{s}^{(0)})$. Here, $\tilde{I}^{(0)} = 1.01$, $\mu_2 = 0.2$, $\mu_3 = -0.02$ and $\tilde{b} = 4\tilde{y}^2$.</figcaption></figure>

The pressure can be obtained from (4.15),

$$ \tilde{p}^{(0)} = \frac{1}{C}(\tilde{s}^{(0)} - \tilde{z}), \quad (4.18) $$

taking into account the dynamic boundary condition in the z-direction (3.28), which reduces to a vanishing pressure at the surface, $\tilde{p}(\tilde{s}^{(0)}) = 0$. By integrating vertically (4.17) and taking into account the no-slip condition at the bottom (2.12), a Bagnold velocity profile is recovered,

$$ \tilde{u}^{(0)} = \frac{2}{3\sqrt{C}} \tilde{I}^{(0)}[(\tilde{s}^{(0)} - \tilde{b})^{\frac{3}{2}} - (\tilde{s}^{(0)} - \tilde{z})^{\frac{3}{2}}]. \quad (4.19) $$

At the leading order, the flow corresponds to a juxtaposition of non-interacting Bagnold velocity profiles, scaled by the cross-stream varying flow depth; that is, cross-slope stresses do not play a role (see figure 3). The normal stress differences contribute to decelerating the flow by modifying the vertical pressure distribution, since C is larger than unity according to the typical values of $\mu_2$ and $\mu_3$ (or, respectively, $N_1$ and $N_2$) found in DEM simulations (see Appendix B). In practice, based on these simulations, $1 < C < 1.04$ such that this effect is unlikely to be measurable.

To fully describe the base state, an equation describing the flow surface $s^{(0)}$ is required, which will come from looking at the equations at the next order.
## 4.2. Linear correction and secondary flows

As the system of equations is now balanced at the leading order, only terms of order $O(\epsilon)$ are kept. The relevant first-order corrections to the shear stress tensor are

$$ \tilde{\sigma}^{xz(1)} = \mu_1^{(1)} \tilde{p}^{(0)} + \mu_1^{(0)} \tilde{p}^{(1)}, \quad (4.20) $$

$$ \tilde{\sigma}^{yz(1)} = \frac{\tilde{p}^{(0)}}{\tilde{u}_{\tilde{z}}^{(0)}} (N_2^{(0)} \tilde{u}_{\tilde{y}}^{(0)} + \mu_1^{(0)} \tilde{v}_{\tilde{z}}^{(1)}), \quad (4.21) $$

$$ \tilde{\sigma}^{zz(1)} = (1-C)\tilde{p}^{(0)} - C\tilde{p}^{(1)}, \quad (4.22) $$

which lead to the momentum equations

$$ \tilde{\sigma}_{\tilde{y}}^{xy(0)} + \tilde{\sigma}_{\tilde{z}}^{xz(1)} = [\mu_1^{(1)} \tilde{p}^{(0)} + \mu_1^{(0)} \tilde{p}^{(1)}]_{\tilde{z}} = 0, \quad (4.23) $$

$$ \tilde{\sigma}_{\tilde{y}}^{yy(0)} + \tilde{\sigma}_{\tilde{z}}^{yz(1)} = -B\tilde{p}_{\tilde{y}}^{(0)} + \left[ \frac{\tilde{p}^{(0)}}{\tilde{u}_{\tilde{z}}^{(0)}} (\mu_1^{(0)} \tilde{v}_{\tilde{z}}^{(1)} + N_2^{(0)} \tilde{u}_{\tilde{y}}^{(0)}) \right]_{\tilde{z}} = 0, \quad (4.24) $$

$$ \tilde{\sigma}_{\tilde{y}}^{zy(0)} + \tilde{\sigma}_{\tilde{z}}^{zz(1)} = [C\tilde{p}^{(1)} + (1-C)\tilde{p}^{(0)}]_{\tilde{z}} = 0, \quad (4.25) $$

and the mass balance equation

$$ \tilde{v}_{\tilde{y}}^{(1)} + \tilde{w}_{\tilde{z}}^{(2)} = 0. \quad (4.26) $$

Noting that $\tilde{p}_{\tilde{y}}^{(0)} = \tilde{s}_{\tilde{y}}^{(0)}/C$ is independent of $\tilde{z}$ and that the pressure vanishes at the surface, the y-momentum (4.24) can be integrated vertically from the surface to obtain

$$ \tilde{v}_{\tilde{z}}^{(1)} = \frac{1}{\mu_1^{(0)}} \left[ \frac{B}{C} \tilde{s}_{\tilde{y}}^{(0)} (\tilde{z} - \tilde{s}^{(0)}) \frac{\tilde{u}_{\tilde{z}}^{(0)}}{\tilde{p}^{(0)}} - N_2^{(0)} \tilde{u}_{\tilde{y}}^{(0)} \right]. \quad (4.27) $$

This can be integrated vertically from the base, where the velocity vanishes, leading to

$$ \tilde{v}^{(1)} = \frac{\tilde{I}^{(0)}}{3\mu_1^{(0)}\sqrt{C}} \left( 3N_2^{(0)}\sqrt{\tilde{s}^{(0)} - \tilde{b}(\tilde{z} - \tilde{b})(\tilde{b}_{\tilde{y}} - \tilde{s}_{\tilde{y}}^{(0)})} - 2C\tilde{s}_{\tilde{y}}^{(0)}[(\tilde{s}^{(0)} - \tilde{b})^{\frac{3}{2}} - (\tilde{s}^{(0)} - \tilde{z})^{\frac{3}{2}}] \right). \quad (4.28) $$

Combining this expression for the cross-stream velocity and the mass balance (4.26) gives an expression for the normal velocity

$$ \tilde{w}^{(2)} = \frac{\tilde{I}^{(0)}}{\mu_1^{(0)} \sqrt{C(\tilde{s}^{(0)} - \tilde{b})}} \left( \frac{N_2^{(0)}}{4} A + \frac{C}{15} B \right), \quad (4.29) $$

where

$$ \begin{aligned} \mathcal{A} ={}& \tilde{b}[\tilde{b}(\tilde{b}_{\tilde{y}} - \tilde{s}_{\tilde{y}}^{(0)})^2 + 2(\tilde{b} - \tilde{s}^{(0)})(\tilde{b}[\tilde{b}_{\tilde{y}\tilde{y}} - \tilde{s}_{\tilde{y}\tilde{y}}^{(0)}] + 2\tilde{b}_{\tilde{y}}[\tilde{b}_{\tilde{y}} - \tilde{s}_{\tilde{y}}^{(0)}])] \\ & + \tilde{z}[(\tilde{z} - 2\tilde{b})(\tilde{b}_{\tilde{y}} - \tilde{s}_{\tilde{y}}^{(0)})^2 + 2(\tilde{b} - \tilde{s}^{(0)})(2\tilde{b}_{\tilde{y}}[\tilde{s}_{\tilde{y}}^{(0)} - \tilde{b}_{\tilde{y}}] + [\tilde{z} - 2\tilde{b}][\tilde{b}_{\tilde{y}\tilde{y}} - \tilde{s}_{\tilde{y}\tilde{y}}^{(0)}])], \end{aligned} \quad (4.30) $$

and

$$ \begin{aligned} \mathcal{B} ={}& \sqrt{\tilde{s}^{(0)} - \tilde{b}}(4\tilde{s}_{\tilde{y}\tilde{y}}^{(0)}[\tilde{s}^{(0)} - \tilde{z}]^{5/2} + 10(\tilde{s}_{\tilde{y}}^{(0)})^2(\tilde{s}^{(0)} - \tilde{z})^{3/2} \\ & - 5\tilde{z}\sqrt{\tilde{s}^{(0)} - \tilde{b}}[3\tilde{b}_{\tilde{y}}\tilde{s}_{\tilde{y}}^{(0)} + 2\tilde{s}_{\tilde{y}\tilde{y}}^{(0)}(\tilde{b} - \tilde{s}^{(0)}) - 3(\tilde{s}_{\tilde{y}}^{(0)})^2]) \\ & - (\tilde{s}^{(0)} - \tilde{b})(2[\tilde{s}^{(0)} - \tilde{b}][\tilde{s}_{\tilde{y}\tilde{y}}^{(0)}(2\tilde{s}^{(0)} + 3\tilde{b}) + 5(\tilde{s}_{\tilde{y}}^{(0)})^2] + 15\tilde{b}\tilde{s}_{\tilde{y}}^{(0)}[\tilde{s}_{\tilde{y}}^{(0)} - \tilde{b}_{\tilde{y}}]). \end{aligned} \quad (4.31) $$
The first non-zero corrections to the horizontal and normal velocities, $\tilde{v}^{(1)}$ and $\tilde{w}^{(2)}$, only depend on the quantities at the leading order. Computing first-order corrections to quantities that are already non-zero at the leading order, such as the pressure or the streamwise velocity, is then left for future work.

### 4.3. Flow surface deformation

At the linear order, the kinematic boundary condition (3.25) becomes

$$ \tilde{v}^{(1)} \tilde{s}_{\tilde{y}}^{(0)} = \tilde{w}^{(2)}, \quad \tilde{z} = \tilde{s}^{(0)}. \qquad (4.32) $$

Integrating the mass balance (4.26) and using Leibniz's integral rule leads to

$$ \frac{\partial}{\partial \tilde{y}} \left( \int_{\tilde{b}}^{\tilde{s}^{(0)}} \tilde{v}^{(1)} d\tilde{y} \right) - \left[ \tilde{v}^{(1)} \tilde{s}_{\tilde{y}}^{(0)} - \tilde{w}^{(2)} \right]_{\tilde{b}}^{\tilde{s}^{(0)}} = 0. \qquad (4.33) $$

Using the kinematic boundary condition (4.32) and no-slip boundary condition (2.12) simplifies the previous equation into

$$ \frac{\partial}{\partial \tilde{y}} \left( \int_{\tilde{b}}^{\tilde{s}^{(0)}} \tilde{v}^{(1)} d\tilde{y} \right) = 0, \qquad (4.34) $$

implying that the cross-stream flux is constant in the cross-stream direction. As there is no flux through the lateral edges, at $\tilde{y} = \pm \tilde{W}/2$, the cross-stream flux is zero everywhere,

$$ \int_{\tilde{b}}^{\tilde{s}^{(0)}} \tilde{v}^{(1)} d\tilde{y} = 0. \qquad (4.35) $$

Combining this with the expression of $\tilde{v}^{(1)}$ obtained previously, (4.28), leads to a linear relation between the surface and the bottom gradients, which can be integrated as

$$ \tilde{s}^{(0)} = \frac{5N_2^{(0)}}{5N_2^{(0)} + 4C} [\tilde{b} - \tilde{b}(0)] + \tilde{s}^{(0)}(0). \qquad (4.36) $$

Finally, the problem can be closed by determining $\tilde{s}^{(0)}(0)$ and $\tilde{W}^{(0)}$ from the mass flux conservation and the boundary condition at the lateral edges,

$$ \int_{-\frac{\tilde{W}^{(0)}}{2}}^{\frac{\tilde{W}^{(0)}}{2}} \int_{\tilde{b}}^{\tilde{s}^{(0)}} \tilde{u}^{(0)} d\tilde{z} d\tilde{y} = 1, \qquad (4.37) $$

$$ \tilde{s}^{(0)} = \tilde{b}, \quad \tilde{y} = \pm \frac{\tilde{W}}{2}. \qquad (4.38) $$

Hence, (4.28), (4.29) and (4.36) show that the free surface deformation and the secondary vortices occur in the system if both the base gradient $b_y$ and the second-normal stress difference $N_2$ are non-zero. Note that, unlike the secondary vortices, the free surface deformation is a leading-order effect in the flow aspect ratio $\epsilon$. The sign of $N_2$ then determines if the surface deforms upward or downward. For granular material, $N_2$ is negative, leading to an upward deformation of the surface (see figures 3 and 4) in agreement with the experimental observations (see figure 1 and § 6), numerical simulations (see § 5) and previous studies (McElwaine et al. 2012; Kim & Kamrin 2023). Importantly, this also determines the direction of rotation of the vortices. The latter can be investigated by looking at the surface velocity from (4.28). It shows that, in the shallow
<figure><img src="image_10.png" /><figcaption>Figure 4. Secondary flows induced by the second-normal stress difference in the case of a parabolic channel. (a) Cross-stream velocity $\tilde{v}^{(1)}$ using (4.28). (b) Normal velocity $\tilde{w}^{(2)}$ using (4.29). In both panels, the red and brown lines represent the flow surface and channel base, respectively. The black lines are streamlines, oriented anticlockwise (plain lines) and clockwise (dashed lines). Here, $\tilde{I}^{(0)} = 1.01$, $\mu_2 = 0.2$, $\mu_3 = -0.02$ and $\tilde{b} = 4\tilde{y}^2$.</figcaption></figure>

limit, the vortices are always downwelling in the centre and going up on the sides if the surface is curving upward (figure 4), and inversely if the surface is curving downward.

Finally, the absolute value of $N_2$ determines the magnitude of the surface deformation and the strength of the vortices. According to DEM simulations available in the literature (Srivastava et al. 2021; Kim & Kamrin 2023), these effects are therefore expected to become stronger as the inertial number of the flow increases, for example, by increasing the channel inclination $\theta$. This is observed in laboratory experiments, as shown and discussed in § 6.

The dimensionless model derived in § 4 therefore allows a prediction of the surface deformation and flow structure upon specifying the base shape $\tilde{b}(\tilde{y})$ and the functions $\mu_1(I), \mu_2(I)$ and $\mu_3(I)$, from which $\tilde{I}^{(0)}$ can be computed using (3.12) and (4.16). Note that, as the inertial number is constant across the whole flow, it is also possible to specify values for $\tilde{I}^{(0)}, \mu_1, \mu_2$ and $\mu_3$, providing those are consistent with the value of $\theta$ (see (4.16)).

## 5. Model validation against DEM simulations

Due to the opacity of granular materials, measuring internal velocity fields remains a challenge in laboratory experiments. To verify the continuum solution derived in § 4, grain-scale numerical simulations are carried out using a DEM from the open-source molecular dynamics code LAMMPS (Thompson et al. 2022).

### 5.1. *Numerical set-up*

The simulations consist of a three-dimensional dense flow of frictional spheres down an inclined rough parabolic channel. Details of the model, including contact forces and
chosen parameters, are given in Appendix A. To achieve a proper comparison with the
continuum modelling, the flow height in the centre is kept around 25d to minimise finite
size effects while keeping computation time reasonable. Additionally, the flow is kept
shallow by having $H/W \simeq 0.09$. A typical DEM simulation is shown in movie 2 of the
Supplementary material.

The DEM simulations provide instantaneous information about the grain velocities and
mechanical stresses, which are turned into averaged continuum fields following the method
described in Appendix A. In particular, this gives access to the velocity and pressure, but
also to the spatial distribution of the inertial number, rheological coefficients and normal
stress differences, allowing direct comparison with the predictions of the model derived in
§§ 2 and 4.

5.2. Results

The continuum model in §§ 2 and 4 is derived under the assumption of a constant volume
fraction across the whole domain. As shown by figure 5(a), this is verified in the DEM
simulations, where the volume fraction fluctuates in the domain around a constant value
$\Phi \in [0.53, 0.55]$ (depending on the slope angle), supporting the use of an incompressible
rheology for dense granular flows. These fluctuations mostly come from layering at the
subparticle scale of particles close to the base, an effect present despite the polydispersity
and the randomness of the roughness of the base, as observed in previous studies (Weinhart
et al. 2013). In addition, it also predicts the inertial number (and hence the friction
coefficients and normal stress differences) to be constant at leading order (see 4.16). As
shown by figure 5(b–e), this is also mainly the case in the DEM simulations, although
some discrepancies exist. The inertial number is slightly larger at the edges and close
to the surface, where the volume fraction is smaller and the shallow flow approximation
breaks down. The second-order rheological coefficients and the normal stress differences
fields exhibit a more significant spatial structure, with discrepancies between the top and
the base of the flow.

The volume fraction, inertial number and rheological coefficients can therefore be
averaged spatially. The resulting values are then used together with the values of *d*, $\rho_p$,
*g*, $\theta$, *Q* and the base shape *b*(*y*) (see Appendix A) as an input for the continuum model
to predict the surface height profile and the pressure and velocity fields. As shown by
figure 6, the model and the DEM data agree quantitatively, especially for the leading-
order quantities, the streamwise velocity and the pressure. Note that the cross-stream and
normal velocities predicted by the model are larger than the ones observed in the DEM
simulations (see discussion in § 7 for additional information). However, the model is still
able to reproduce accurately the spatial structures of the vortices. The predictions of the
leading-order quantities remain accurate when changing the slope angle from 23° to 29°.
However, the discrepancy between the predicted and observed velocities in the vortices
increases as the slope angle decreases. This is attributed to finite-size effects ($H/d$ not
large enough) and further discussed in § 7.

# 6. Laboratory experiments and inference of the second-normal stress difference

An interesting result of the solution derived in § 4 is the simple linear relationship between
the channel shape and the surface deformation (4.36). Based on the DEM simulations,
prediction of the surface shape using this relationship is especially accurate (see figure 7).
Hence, assuming *C* to be unity (see the end of § 4.1), the scaled second-normal stress
difference can be inferred from laboratory experiments by measuring the surface and basal
topographies.
<figure><img src="image_11.png" /><figcaption>Figure 5. Uniform fields from a DEM simulation with $\theta = 29^\circ$. (a) Solids volume fraction, $\Phi$. (b) Inertial number, $I$. (c) First-order friction coefficient, $\mu_1$. (d,e) Second-order rheological coefficients, $\mu_2$ and $\mu_3$. (f,g) Scaled normal stress differences, $N_1$ and $N_2$. In all panels, the red and brown lines represent the flow surface and channel base, respectively. Note that here, the aspect ratio has been increased for visualisation purposes, as the proper flow is much shallower ($\epsilon \approx 0.09$). A video of this DEM simulation is shown in movie 2 of the Supplementary material.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

Figure 6. Comparison between (a,c,e,g) dimensionless DEM results and (b,d,f,h) predictions from the model for $\theta = 29^\circ$ using as an input the values averaged from the uniform fields presented in figure 5. (a,b) Pressure, $\tilde{p}$, and (4.18). (c,d) Streamwise velocity $\tilde{u}$ and (4.19). (e,f) Cross-stream velocity $\tilde{v}$ and (4.28). (g,h) Normal velocity $\tilde{w}$ and (4.29). In all panels, the red and brown lines represent the flow surface and channel base, respectively. Note that here, the aspect ratio has been increased for visualisation purposes, as the proper flow is much shallower ($\varepsilon \approx 0.09$). A video of this DEM simulation is shown in movie 2 of the Supplementary material.
C. Gadal, C.G. Johnson and J.M.N.T. Gray

<figure><img src="image_24.png" /><figcaption>Figure 7. Comparison between dimensionless surface gradients $\tilde{s}_{\bar{y}}$ extracted from DEM simulations and predicted from the model using as an input the values averaged from the uniform fields presented in figure 5 together with the values of *d*, $\rho_p$, *g*, $\theta$, *Q* and the base shape *b*(*y*) used in the DEM (see Appendix A), for five different slopes. Each point is the value of the gradient at a single position along the *y*-axis. The black dashed line is the identity line.</figcaption></figure>

## 6.1. *Experimental set-up*

The experiments are carried out using a channel of curved parabolic cross-section, tilted at an angle $\theta \in [23^\circ, 30^\circ]$. Spherical glass beads, of diameter $d \in [125, 165]$ µm, are released from a funnel at the top of the chute behind a gate of controlled height. As they flow out of the funnel, they quickly pile up behind the door and reduce the flux coming out of the funnel, equilibrating with what comes out of the gate. This set-up ensures both a constant mass flux and a controlled height at the inlet, reducing the spatial relaxation length needed for the flow to reach the gravity–friction equilibrium. The mass flux is then self-selected by the slope and the gate height, which is varied from 0.75 to 1 cm.

Height measurements are performed using a laser scanner scanCONTROL 2950–100 from Micro-Epsilon at a frequency of 120 Hz, which has a horizontal spatial resolution of approximately 0.07 mm (depending on the experiments). The vertical resolution is of the order of 10 µm, but the precision of the measurements is dominated by temporal random fluctuations, of the order of the grain diameter. The measurements are averaged during the steady part (after the passage of the front and before the running out of the flow), which lasts from 60 to 100 s, corresponding to a total mass of approximately 5.7 kg of released particles. Similarly, the mass flux is recorded using a scale and averaged during the same temporal period.

## 6.2. *Surface height profiles*

As shown by figure 8(a), the measured surface height profiles all exhibit an upward curvature, as observed in the DEM simulations, and predicted by (4.36) in the case of negative $N_2$. However, profiles with different flow parameters are difficult to compare, particularly across different fluxes.

The profiles can be made dimensionless using the length scales derived in § 3.3. Their estimation requires a value for the first-order friction coefficient $\mu_1$, computed from the channel inclination following (4.16), and for the inertial number, which can be inferred from the mass flux measurements. The streamwise velocity is assumed to follow a Bagnold
<figure><img src="image_25.png" /><figcaption>Figure 8. Surface height profiles for various slopes and fluxes in dimensional coordinates. The shades of red represent the corresponding values of the inertial number. (a) Dimensional profiles. (b) Dimensionless profiles, with $H \in [2.1, 6.3]$ cm and $W \in [21, 34]$ cm. (c) Dimensionless profiles zoomed on the surface to highlight the evolution of the curvature with the inertial number. In all panels, the brown curves represent measurements of the empty channel section. Not all data are shown for the sake of clarity.</figcaption></figure>

profile scaled by the local flow depth as predicted by (4.19). Using the definition of the mass flux (3.9) and integrating, the inertial number can be estimated as

$$I = \frac{5}{2} d \frac{Q}{\Phi \rho_p} \frac{1}{\sqrt{\Phi g \cos(\theta) S}}, \quad (6.1)$$

where

$$S = \int_{-W/2}^{W/2} (s-b)^{5/2} dy, \quad (6.2)$$

is computed from the experimental data.

As shown by figure 8(b), all profiles collapse towards a single curve, indicating that most of the flow properties are encompassed in the gravity–friction balance used to derive the different scales. However, a careful examination of the surface shows a consistent increase of the surface curvature with the inertial number (see figure 8c).

### 6.3. Second-normal stress difference

The second-normal stress difference can then be inferred from the surface height profiles. For each experimental run, the flow surface $s^{(0)}$ and channel section $b$ can both be approximated by a parabola (figure 9a), and their gradients $s_y^{(0)}$ and $b_y$ are found to vary linearly with $y$ (figure 9b). Furthermore, their ratio is found to be constant, as predicted by (4.36) (figure 9c). Hence, both are fitted with a second-order polynomial,
<figure><img src="image_26.png" /><figcaption>Figure 9. Experimental measurements for $\theta = 23.4^\circ$ and $Q = 85 \text{ g s}^{-1}$. (a) Flow surface and channel section. The black dashed lines represent the fits of (6.3) and (6.4), leading to $B = 0.02 \pm 0.00005 \text{ mm}^{-1}$ and $S = -0.003 \pm 0.0003 \text{ mm}^{-1}$. (b) Cross-stream gradients of the flow surface and channel section. (c) Ratio $r = s_y/b_y$. The black dashed line represents the ratio between $S/B$ from the fits in (a).</figcaption></figure>

$$b(y) = \frac{B}{2}(y - y_0)^2 + b(0), \qquad (6.3)$$

$$s(y) = \frac{S}{2}(y - y_0)^2 + s(0), \qquad (6.4)$$

where $B$ and $S$ are the base and flow surface gradients, and $s(0)$, $b(0)$ and $y_0$ are constants accounting for spatial shifts. Assuming $C$ to be unity (see Appendix B), the scaled second-normal stress difference can then be computed from (4.36) as

$$N_2 = \frac{4r}{5(r-1)}, \qquad (6.5)$$

where $r = s_y/b_y = S/B$.

The first-order friction coefficient $\mu_1$ and the second-normal stress difference $N_2$ are shown as a function of the inertial number in figure 10. The results for $\mu_1$ are in agreement with friction laws obtained experimentally in the literature with inclined plane configurations (e.g. GDR MiDi 2004). Additionally, the measurements presented here reach higher $I$-values than those performed by GDR MiDi (2004). At such values, a linear trend is found rather than a saturation, in agreement with recent measurements performed on superstable granular heaps (Lloyd et al. 2025). This linear trend may result from the hypothesis of the local rheology under which all these measurements are inferred, as discussed by GDR MiDi (2004). To the authors' knowledge, no measurements of $N_2$ are available in the literature to compare with. Unfortunately, direct comparison with the results of DEM simulations is not relevant, as the results of the latter strongly depend on the particle-particle tangential friction coefficient $\mu_p$, which encompasses several physical quantities such as the particle surface roughness or geometry (Srivastava et al. 2021). However, $N_2$ is of the correct order of magnitude, and its variation with $I$ displays
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>Figure 10. (a,b) Measured first-order friction coefficient as a function of the inertial number. (c,d) Measured second-normal stress difference as a function of the inertial number. Note that the horizontal axes of (b) and (d) are in logarithmic scale. Errorbars represent the ±1σ uncertainty. The orange dashed line shows the fit of the power law (B2), leading to $N_2^* = 0 \pm 0.1$, $A_2 = -0.14 \pm 0.1$, $\alpha_2 = 0.1 \pm 0.1$.</figcaption></figure>

similar behaviour to the ones found in numerical simulations (Nagy *et al.* 2020; Srivastava *et al.* 2021; Kim & Kamrin 2023). Similar to numerical simulations (see Appendix B and Srivastava *et al.* (2021)), the shape of $N_2(I)$ can be represented with a power-law function (see figure 10c–d). However, additional data at smaller and larger inertial numbers are required to properly constrain the fit.

## 7. Conclusion and perspectives

In this work, the flow of a dense granular material down an inclined channel of arbitrary smooth cross-section is studied. Asymptotic solutions are derived from a continuum model with a second-order granular rheology accounting for the existence of normal stress differences. These solutions rely on the shallowness of the flow, which is used to expand the system of equations linearly. The resulting model can predict both the surface deformation and the existence of two counter-rotating vortices in the plane perpendicular to the flow direction, which are observed in numerical simulations and experiments in the existing literature. Importantly, the existence of these two features is linked to the presence of a cross-stream gradient of the streamwise velocity. Here, a non-flat topography is used to create this gradient, but similar results are expected in the case of a flat channel with frictional walls (Deboeuf *et al.* 2006; McElwaine *et al.* 2012; Kim & Kamrin 2023). The solution presented in this paper could be extended to this configuration, but this would require the computation of higher-order terms as the required velocity gradient would only come up at order $\epsilon$, pushing the secondary flows to higher orders.
C. Gadal, C.G. Johnson and J.M.N.T. Gray

The validity of these solutions is tested against DEM simulations. A quantitative match is found between the two, especially concerning the prediction of the surface deformation and the structure of the vortices. The velocities involved in the secondary vortices are, however, smaller in the DEM than predicted by the theory. This could be attributed to the aspect ratio of the flow, not being shallow enough, but also to finite-size effects. For a fixed aspect ratio $\varepsilon$, decreasing the size of the system (i.e. decreasing $H/d$) also reduces the strength of the vortices up to their disappearance. This may be induced by the breakdown of the continuum assumption and/or the increase (relative to the strength of the vortices) of the self-diffusion of particles in the plane perpendicular to the flow (e.g. Utter & Behringer 2004; Artoni et al. 2021), but this requires a dedicated study.

An interesting feature of the predicted and observed secondary flow is the direction of the vortices, always downwelling in the centre for a negative second-normal stress difference in the shallow limit. Kim & Kamrin (2023) have found vortices rotating in the opposite direction in a configuration deeper than it is wide. This suggests that the direction of the vortices could be driven not only by the sign of the normal stress difference but also by the flow aspect ratio (Lévay et al. 2025) or, more generally, by the shape of the container, as observed in the case of viscoelastic fluids (Yue, Dooley & Feng 2008). This interesting perspective is left for future work.

The continuum model also predicts a simple relationship between the cross-stream gradients of the surface and channel elevations. By measuring both in simple laboratory experiments, this relationship is used to infer measurements of the scaled second-normal stress difference in dense granular flows. This provides a simple method that could be used in the future to investigate the properties of the second-normal stress difference, such as its dependence on the particle properties. However, these experiments do not provide any constraint on the first normal stress difference. Hence, an alternative flow geometry providing complementary measurements of $N_1(I)$ would be required to provide a complete calibration of the second-order rheology in flow regimes where $N_1$ is not negligible. This rheology may then be used as an input in continuum simulations to predict the flow of granular materials in arbitrary configurations.

**Supplementary movies.** Supplementary movies are available at https://doi.org/10.1017/jfm.2025.11043.

**Acknowledgements.** The authors acknowledge the use of open-source numerical Python libraries, including Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020) and Uncertainties (https://pythonhosted.org/uncertainties/) as well as the open-source molecular dynamics code LAMMPS used to perform the DEM simulations (Thompson et al. 2022). They also acknowledge the contribution of Aarnav Panda to the preliminary experiments, and Dr J. Webb, Dr D. Schaeffer and the anonymous reviewers for useful suggestions. This work is licenced under CC BY 4.0. To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

**Funding.** C.G. acknowledges funding from the Royal Society through a Newton International Fellowship (NIF/R1/231983). 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). C.G.J. and J.M.N.T.G. acknowledge funding from NERC through grants NE/X013936/1 and NE/X00029X/1.

**Declaration of interests.** The authors report no conflict of interest.

**Data availability statement.** The experimental and numerical data used in this paper are available at https://doi.org/10.5281/zenodo.15639395.

## Appendix A. The DEM simulations

The DEM simulations are run using the molecular dynamics open-source code LAMMPS (29 August 2024 stable release with update 1) (Thompson et al. 2022). In the model, each
sphere i has a diameter $d_i$ sampled within a uniform distribution ranging from 0.8d to 1.2d to prevent crystallisation effects inside the flow, and a mass $m_i = \rho_s(4/3)\pi(d_i/2)^2$. Spheres interact only upon contact, which corresponds to overlapping in the model, i.e. when $\delta_{ij} < 0$, where $\delta_{ij}$ is the overlapping distance

$$ \delta_{ij} = r_{ij} - \frac{d_i + d_j}{2}, \qquad (A1) $$

with $r_{ij}$ the distance between the centres of the spheres. The contact forces are modelled using a combination of a viscoelastic normal force and Coulomb tangential friction, which has proven to be efficient in modelling various granular flows (e.g. Kim & Kamrin 2023).

The viscoelastic normal force corresponds to the classical spring-dashpot model (Cundall & Strack 1979)

$$ F_n^{i \to j} = \begin{cases} -(k_n \delta_{ij} + \gamma_n \delta_{ij}) & \text{for } \delta_{ij} < 0, \\ 0 & \text{else,} \end{cases} \qquad (A2) $$

where $k_n$ and $\gamma_n$ are the spring stiffness and damping coefficient, respectively. The damping coefficient is computed from the restitution coefficient $e$ as (Brilliantov et al. 1996)

$$ \gamma_n = \sqrt{\frac{4m_{eff}k_n}{1 + \frac{\pi^2}{\ln e}}} \qquad (A3) $$

where $m_{eff} = m_i m_j / (m_i + m_j)$ is the effective mass of the two colliding particles.

The tangential force corresponds to a regularised Coulomb friction, also taking into account the contact history (Luding 2008)

$$ F_t^{i \to j} = -\min(|\gamma_t u_t + k_t \xi|, |\mu_p| F_n^{i \to j}), \qquad (A4) $$

where $\gamma_t$ is the tangential damping coefficient, $u_t$ is the tangential component of the relative velocity of the two spheres, $\xi$ the tangential displacement accumulated during the entire duration of the contact and $\mu_p$ is the particle-particle friction coefficient. The tangential damping constant $k_t$ is fixed to (2/7)$k_n$ so that the frequencies of normal and tangential contact oscillations are similar, and the normal and tangential dissipation are comparable. The particles are also submitted to a constant gravitational acceleration of intensity $g$, which can be inclined to simulate different chute inclinations.

The channel is made of spheres of the same size distribution that are fixed on a regular array mapping the channel shape

$$ b(y) = ay^2, \qquad (A5) $$

where $a = 0.00213/d$. To ensure as much as possible a no-slip boundary condition at the base of the flow, the channel is made rough by initially separating the fixed spheres by a distance $0.7d$, and then perturbing this layer by randomly moving each sphere in every direction by a distance uniformly sampled in the range $[-0.4d, 0.4d]$. The details of this process do not matter providing the no-slip boundary condition at the base is verified.

For each particle, the particle motion is calculated by integrating Newton's second law using a velocity-Verlet scheme with a fixed time step $dt = 0.05t_R$, where $t_R$ is the Rayleigh time

$$ t_R = \pi \sqrt{\frac{m_{eff}}{k_n} \left[ 1 + \ln \left( \frac{e}{\pi} \right)^2 \right]}, \qquad (A6) $$
<table><caption>Table 1. Values of the parameters used for the DEM simulations.</caption><thead><tr><th>d [mm]</th><th>ρ<sub>p</sub> [kg m<sup>3</sup>]</th><th>k<sub>n</sub> [kg s<sup>-2</sup>]</th><th>e [-]</th><th>μ<sub>p</sub> [-]</th><th>g [m s<sup>-2</sup>]</th></tr></thead><tbody><tr><td>0.8</td><td>2650</td><td>263000</td><td>0.24</td><td>0.4</td><td>9.81</td></tr></tbody></table>

corresponding to half the period of oscillation of an undamped oscillator of the same properties. Table 1 contains the parameter values used for the simulations. The simulations are kept in the regime of rigid particles, where the rheology is unaffected by the particle stiffness, by ensuring that $k_n/(Pd) > 10^4$ (Favier de Coulomb et al. 2017).

To save computational time, the box is made periodic in the streamwise direction, with a length $\Delta_x = 20d$. The flow width and maximum depth are always approximately 200d and 30d, respectively, which corresponds to approximately $7 \times 10^4$ particles per simulation.

Simulations are run for a duration of $5 \times 10^7 dt$, corresponding to 25 s using the values of the parameters given in Table 1. Particle properties are saved every $2.4 \times 10^4 dt$, and are averaged using the methods described in Kim & Kamrin (2023) to obtain the continuum fields displayed in § 5. The particle mass flux can then be computed as

$$Q = \left\langle \frac{1}{\Delta_x} \sum_i m_i u_i \right\rangle_t, \qquad (A7)$$

where $m_i$ and $u_i$ are the mass and streamwise velocity of particle $i$.

## Appendix B. Normal stress differences in the DEM simulations

The rheological curves $N_1(I)$ and $N_2(I)$ extracted from the DEM simulations (see § 5) are shown in figure 11(a), together with the results from Kim & Kamrin (2023) obtained under linear shear with the same numerical parameters as used in this study. The results agree quantitatively despite the difference in geometrical configurations. Srivastava et al. (2021) have shown that $\mu_2(I)$ and $\mu_3(I)$ can be fitted using power-law functions (see figure 11a) as

$$\mu_i = \mu_i^* + A_i I^{2\alpha_i}, \qquad (B1)$$

where $\mu_i^*$, $A_i$ and $\alpha_i$ are fitting constants. Similarly, both $N_1(I)$ and $N_2(I)$ can be fitted using power-law functions (see figure 11b) as

$$N_i = N_i^* + B_i I^{2\beta_i}, \qquad (B2)$$

where $N_i^*$, $B_i$ and $\beta_i$ are fitting constants.

As shown by the blue points in figure 11(b), the constant C defined by (4.11) is slightly larger than unity (within 5%) in the numerical simulations. From (4.12) and (4.13), the constant C can also be expressed using the scaled normal stress differences as

$$C = 1 + \frac{1}{3}(N_1^{(0)} - N_2^{(0)}). \qquad (B3)$$

Then, to infer the value of $N_2$ from the experiments in § 6 using (4.36), one must either provide a value for $N_1$, or adopt an assumption removing the dependency in the first-normal stress difference. As $N_1 \ll N_2$ for $I < 10^{-1}$, one could simply set $N_1 = 0$ in (4.11). As shown by the orange points in figure 11(b), this provides a good approximation until $I \simeq 10^{-1}$. For larger values of the inertial number, $N_1$ is not negligible anymore but of the same order of magnitude as $N_2$, which leads to a large difference with the exact value
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

Figure 11. (a) Second-order rheological coefficients $\mu_2$ (blue) and $\mu_3$ (orange) obtained in DEM simulations as a function of the inertial number. The blue and orange dashed lines show the fits of the power law (B1), leading to $\mu_2^*=0.103 \pm 0.006$, $A_2=0.27 \pm 0.02$, $\alpha_2=0.40 \pm 0.04$ and $\mu_3^*=0 \pm 0.003$, $A_3=-0.05 \pm 0.005$, $\alpha_3=0.43 \pm 0.07$. (b) First (blue) and second (orange) normal stress differences obtained in DEM simulations as a function of the inertial number. The blue and orange dashed lines show the fits of the power law (B2), leading to $N_1^*=0.01 \pm 0.01$, $B_1=-0.26 \pm 0.03$, $\beta_1=0.46 \pm 0.08$ and $N_2^*=-0.1 \pm 0.01$, $B_2=-0.2 \pm 0.015$, $\beta_2=0.37 \pm 0.06$. (c) The C constant as a function of the inertial number computed using the exact expression (4.11) (blue) and by setting $N_1=0$ or equivalently $\mu_3=0$ (orange). In both panels, the filled circles are the simulations presented in this article, and the empty squares are data obtained under linear shear by Kim & Kamrin (2023) with the same numerical parameters as used in this study. The grey shaded area shows the range of inertial numbers corresponding to our experiments.

of C. Hence, a better approximation at larger inertial number is $(N_1 - N_2)/3 \ll 1$, hence $C \approx 1$, which keeps the error within 5%.

REFERENCES

ARTONI, R., LARCHER, M., JENKINS, J.T. & RICHARD, P. 2021 Self-diffusion scalings in dense granular flows. *Soft Matter* **17** (9), 2596–2602.

BARKER, T., RAUTER, M., MAGUIRE, E.S.F., JOHNSON, C.G. & GRAY, J.M.N.T. 2021 Coupling rheology and segregation in granular flows. *J. Fluid Mech.* **909**, A22.

BARKER, T., SCHAEFFER, D.G., SHEARER, M. & GRAY, J.M.N.T. 2017 Well-posed continuum equations for granular flow with compressibility and (I)-rheology. *Proc. R. Soc. A: Math. Phys. Engng Sci.* **473** (2201), 20160846.
BARKER, T., SCHAEFFER, D.G., BOHÓRQUEZ, P. & GRAY, J.M.N.T. 2015 Well-posed and ill-posed behaviour of the rheology for granular flow. *J. Fluid Mech.* **779**, 794–818.

BÖRZSÖNYI, T., ECKE, R.E. & MCELWAINE, J.N. 2009 Patterns in flowing sand: understanding the physics of granular flow. *Phys. Rev. Lett.* **103** (17), 178302.

BOUZID, M., IZZET, A., TRULSSON, M., CLÉMENT, E., CLAUDIN, P. & ANDREOTTI, B. 2015 Non-local rheology in dense granular flows: revisiting the concept of fluidity. *Eur. Phys. J. E* **38**, 1–15.

BOYER, F., POULIQUEN, O. & GUAZZELLI, É. 2011 Dense suspensions in rotating-rod flows: normal stresses and particle migration. *J. Fluid Mech.* **686**, 5–25.

BRILLIANTOV, N.V., SPAHN, F., HERTZSCH, J.-M. & PÖSCHEL, T. 1996 Model for collisions in granular gases. *Phys. Rev. E* **53** (5), 5382.

BRODU, N., DELANNAY, R., VALANCE, A. & RICHARD, P. 2015 New patterns in high-speed granular flows. *J. Fluid Mech.* **769**, 218–228.

CABRERA, M. & POLANÍA, O. 2020 Heaps of sand in flows within a split-bottom couette cell. *Phys. Rev. E* **102** (6), 062901.

CAMPBELL, C.S. 1989 The stress tensor for simple shear flows of a granular material. *J. Fluid Mech.* **203**, 449–473.

FAVIER DE COULOMB, A., BOUZID, M., CLAUDIN, P., CLÉMENT, E. & ANDREOTTI, B. 2017 Rheology of granular flows across the transition from soft to rigid particles. *Phys. Rev. Fluids* **2** (10), 102301.

COUTURIER, É., BOYER, F., POULIQUEN, O. & GUAZZELLI, É. 2011 Suspensions in a tilted trough: second normal stress difference. *J. Fluid Mech.* **686**, 26–39.

CUNDALL, P.A. & STRACK, O.D.L. 1979 A discrete numerical model for granular assemblies. *Geotechnique* **29** (1), 47–65.

DAI, S.-C., BERTEVAS, E., QI, F. & TANNER, R.I. 2013 Viscometric functions for noncolloidal sphere suspensions with Newtonian matrices. *J. Rheol.* **57** (2), 493–510.

DBOUK, T., LOBRY, L. & LEMAIRE, E. 2013 Normal stresses in concentrated non-brownian suspensions. *J. Fluid Mech.* **715**, 239–272.

DEBOEUF, S., LAJEUNESSE, E., DAUCHOT, O. & ANDREOTTI, B. 2006 Flow rule, self-channelization, and levees in unconfined granular flows. *Phys. Rev. Lett.* **97** (15), 158303.

DSOUZA, P.V. & NOTT, P.R. 2021 Dilatancy-driven secondary flows in dense granular materials. *J. Fluid Mech.* **914**, A36.

EDWARDS, A.N. & GRAY, J.M.N.T. 2015 Erosion-deposition waves in shallow granular free-surface flows. *J. Fluid Mech.* **762**, 35–67.

EINAV, I., ESCOBAR, A., BAKER, J., GUILLARD, F. & FAUG, T. 2024 *Experimental Confirmation of Secondary Flows in Granular Media*. Research Square.

FÉLIX, G. & THOMAS, N. 2004 Relation between dry granular flow regimes and morphology of deposits: formation of levées in pyroclastic deposits. *Earth Planet. Sc. Lett.* **221** (1–4), 197–213.

FISCHER, D., BÖRZSÖNYI, T., NASATO, D.S., PÖSCHEL, T. & STANNARIUS, R. 2016 Heaping and secondary flows in sheared granular materials. *New J. Phys.* **18** (11), 113006.

FORTERRE, Y. & POULIQUEN, O. 2001 Longitudinal vortices in granular flows. *Phys. Rev. Lett.* **86** (26), 5886.

FORTERRE, Y. & POULIQUEN, O. 2008 Flows of dense granular media. *Annu. Rev. Fluid Mech.* **40** (1), 1–24.

MIDI, G.D.R. 2004 On dense granular flows. *Eur. Phys. J. E* **14** (4), 341–365.

GUAZZELLI, É. & POULIQUEN, O. 2018 Rheology of dense granular suspensions. *J. Fluid Mech.* **852**, P11–P173.

HARRIS, C.R. et al. 2020 Array programming with NumPy. *Nature* **585** (7825), 357–362.

HEYMAN, J., DELANNAY, R., TABUTEAU, H. & VALANCE, A. 2017 Compressibility regularizes the $\mu(I)$-rheology for dense granular flows. *J. Fluid Mech.* **830**, 553–568.

HUNTER, J.D. 2007 Matplotlib: a 2D graphics environment. *Comput. Sci. Engng* **9** (03), 90–95.

JOP, P., FORTERRE, Y. & POULIQUEN, O. 2006 A constitutive law for dense granular flows. *Nature* **441** (7094), 727–730.

KAMRIN, K. 2019 Non-locality in granular flow: phenomenology and modeling approaches. *Front. Phys.* **7**, 116.

KIM, S. & KAMRIN, K. 2020 Power-law scaling in granular rheology across flow geometries. *Phys. Rev. Lett.* **125** (8), 088002.

KIM, S. & KAMRIN, K. 2023 A second-order non-local model for granular flows. *Front. Phys.* **11**, 1092233.

KOMATSU, T.S., INAGAKI, S., NAKAGAWA, N. & NASUNO, S. 2001 Creep motion in a granular pile exhibiting steady surface flow. *Phys. Rev. Lett.* **86** (9), 1757.

KRISHNARAJ, K.P. & NOTT, P.R. 2016 A dilation-driven vortex flow in sheared granular materials explains a rheometric anomaly. *Nat. Commun.* **7** (1), 10630.
LAGRÉE, P.-Y., STARON, L. & POPINET, S. 2011 The granular column collapse as a continuum: validity of a two-dimensional navier–stokes model with a $\mu$ (i)-rheology. *J. Fluid Mech.* **686**, 378–408.

LÉVAY, S., CLAUDIN, P., SOMFAI, E. & BÖRZSÖNYI, T. 2025 Transverse vortices induced by modulated granular shear flows of elongated particles. *Phys. Rev. E* **112** (1), 015402.

LLOYD, H.A., MAGUIRE, E.S.F., MISTRY, D., REYNOLDS, G.K., JOHNSON, C.G. & GRAY, J.M.N.T. 2025 On the formation of super-stable granular heaps. *J. Fluid Mech.* **1002**, A27.

LUDING, S. 2008 Cohesive, frictional powders: contact models for tension. *Granul. Matter* **10** (4), 235–246.

MAGUIRE, E.S.F., BARKER, T., RAUTER, M., JOHNSON, C.G. & GRAY, J.M.N.T. 2024 Particle-size segregation patterns in a partially filled triangular rotating drum. *J. Fluid Mech.* **979**, A40.

MAKLAD, O. & POOLE, R.J. 2021 A review of the second normal-stress difference; its importance in various flows, measurement techniques, results for various complex fluids and theoretical predictions. *J. Non-Newtonian Fluid* **292**, 104522.

MARTIN, N., IONESCU, I.R., MANGENEY, A., BOUCHUT, F. & FARIN, M. 2017 Continuum viscoplastic simulation of a granular column collapse on large slopes: $\mu(I)$ rheology and lateral wall effects. *Phys. Fluids* **29**, 013301.

MCELWAINE, J.N., TAKAGI, D. & HUPPERT, H.E. 2012 Surface curvature of steady granular flows. *Granul. Matter* **14** (2), 229–234.

MOHAMMADI, M., PUZYREV, D., TRITTEL, T. & STANNARIUS, R. 2022 Secondary flow in ensembles of nonconvex granular particles under shear. *Phys. Rev. E* **106** (5), L052901.

NAGY, D.B., CLAUDIN, P., BÖRZSÖNYI, T. & SOMFAI, E. 2020 Flow and rheology of frictional elongated grains. *New J. Phys.* **22** (7), 073008.

POULIQUEN, O., CASSAR, C., JOP, P., FORTERRE, Y. & NICOLAS, M. 2006 Flow of dense granular material: towards simple constitutive laws. *J. Statist. Mechan.: Theory Exp.* **2006** (07), P07020.

POULIQUEN, O. & FORTERRE, Y. 2009 A non-local rheology for dense granular flows. *Phil. Trans. R. Soc. A: Math. Phys. Engng Sci.* **367** (1909), 5091–5107.

RIVLIN, R.S. & ERICKSEN, J.L. 1997 Stress-deformation relations for isotropic materials. In *Collected Papers of RS Rivlin: Volume I and II*, pp. 911–1013.

ROCHA, F.M., JOHNSON, C.G. & GRAY, J.M.N.T. 2019 Self-channelisation and levee formation in monodisperse granular flows. *J. Fluid Mech.* **876**, 591–641.

SCHAEFFER, D.G., BARKER, T., TSUJI, D., GREMAUD, P., SHEARER, M. & GRAY, J.M.N.T. 2019 Constitutive relations for compressible granular flow in the inertial regime. *J. Fluid Mech.* **874**, 926–951.

SILBERT, L.E., ERTAŞ, D., GREST, G.S., HALSEY, T.C., LEVINE, D. & PLIMPTON, S.J. 2001 Granular flow down an inclined plane: bagnold scaling and rheology. *Phys. Rev. E* **64** (5), 051302.

SINGH, A. & NOTT, P.R. 2003 Experimental measurements of the normal stresses in sheared stokesian suspensions. *J. Fluid Mech.* **490**, 293–320.

SRIVASTAVA, I., SILBERT, L.E., GREST, G.S. & LECHMAN, J.B. 2021 Viscometric flow of dense granular materials under controlled pressure and shear stress. *J. Fluid Mech.* **907**, A18.

STARON, L., LAGRÉE, P.Y. & POPINET, S. 2014 Continuum simulation of the discharge of the granular silo: a validation test for the $\mu$ (i) visco-plastic flow law. *Eur. Phys. J. E* **37**, 1–12.

TAKAGI, D., MCELWAINE, J.N. & HUPPERT, H.E. 2011 Shallow granular flows. *Phys. Rev. E – Statist. Nonlinear Soft Matter Phys.* **83** (3), 031306.

THOMPSON, A.P. et al. 2022 LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comput. Phys. Commun.* **271**, 108171.

UTTER, B. & BEHRINGER, R.P. 2004 Self-diffusion in dense granular shear flows. *Phys. Rev. E* **69** (3), 031308.

VIRTANEN, P. et al. 2020 SciPy 1.0: fundamental algorithms for scientific computing in python. *Nat. Methods* **17**, 261–272.

WEINHART, T., HARTKAMP, R., THORNTON, A.R. & LUDING, S. 2013 Coarse-grained local and objective continuum description of three-dimensional granular flows down an inclined surface. *Phys. Fluids* **25** (7), 070605.

WORTEL, G., BÖRZSÖNYI, T., SOMFAI, E., WEGNER, S., SZABÓ, B., STANNARIUS, R. & VAN HECKE, M. 2015 Heaping, secondary flows and broken symmetry in flows of elongated granular particles. *Soft Matter* **11** (13), 2570–2576.

WU, W.-T., AUBRY, N., ANTAKI, J.F. & MASSOUDI, M. 2017 Normal stress effects in the gravity driven flow of granular materials. *Intl J. Nonlinear Mech.* **92**, 84–91.

YUE, P., DOOLEY, J. & FENG, J.J. 2008 A general criterion for viscoelastic secondary flow in pipes of noncircular cross section. *J. Rheol.* **52** (1), 315–332.