2.1. The granular Saint-Venant equations and the associated constitutive relations

In analogy with the shallow water approximation for normal fluids, and justified by the fact that the wavelength of the waves under consideration is orders of magnitude larger than the thickness of the sheet, which in turn is markedly larger than the wave amplitude (cf. figure 1), we describe the dynamics of the flowing granular sheet in terms of its height $h(x,t)$ and depth-averaged velocity $\bar {u}(x,t)$. The latter is defined as

(2.1)\begin{equation} \bar{u}(x,t) = \frac{1}{h(x,t)}\int_{0}^{h(x,t)} u(x,z,t) \,{\rm d}z, \end{equation}

where $u(x,z,t)$ is the detailed velocity profile, which also depends on the ordinate $z$. In the description, any velocity variations in the $z$-direction are averaged out, while flow variations in the cross-wise direction are ignored altogether, hence both $h(x,t)$ and $\bar {u}(x,t)$ depend only on $x$ and $t$, cf. figure 1. Furthermore, we treat the granular material as an incompressible fluid (Savage & Hutter Reference Savage and Hutter1989), assuming constant density throughout the flowing sheet.

The height $h(x,t)$ and the depth-averaged velocity $\bar {u}(x,t)$ are governed by a system of two coupled, nonlinear PDEs (Gray & Edwards Reference Gray and Edwards2014; Razis et al. Reference Razis, Edwards, Gray and van der Weele2014; Edwards & Gray Reference Edwards and Gray2015; Razis et al. Reference Razis, Kanellopoulos and van der Weele2018; Edwards et al. Reference Edwards, Russell, Johnson and Gray2019; Razis et al. Reference Razis, Kanellopoulos and van der Weele2019; Edwards et al. Reference Edwards, Viroulet, Johnson and Gray2021; Kanellopoulos Reference Kanellopoulos2021), namely, the mass conservation

(2.2)\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial}{\partial x} (h \bar{u}) = 0, \end{equation}

and the momentum balance

(2.3)\begin{equation} \frac{\partial}{\partial t} (h \bar{u}) + \frac{\partial}{\partial x} (h \bar{u}^{2}) = gh\sin\zeta - \frac{\partial}{\partial x}\left(\frac{1}{2}gh^{2}\cos\zeta\right) - \mu(h, \bar{u})gh \cos\zeta + \frac{\partial}{\partial x}\left(\nu h^{3/2}\frac{\partial \bar{u}}{\partial x}\right). \end{equation}

The left-hand side of (2.3) is the Stokesian derivative of the depth-averaged momentum per unit volume, divided by the constant density. The terms on the right-hand side of (2.3) represent the various depth-integrated forces (per unit volume and divided by the density) acting on the flowing material: (i) the gravity component along the $x$-direction (the direction of the flow); (ii) the force arising from variations in $h(x,t)$, i.e. the negative gradient of the depth-averaged pressure; (iii) the frictional force exerted by the chute bed, modelled as a (height- and velocity-dependent) friction coefficient $\mu (h,\bar {u})$ multiplied by the normal reaction force acting on the granular sheet; and (iv) the viscous-like diffusive term arising from depth averaging the in-plane stresses in the sheet.

The last two terms of (2.3) contain constitutive relations that are specific for granular matter, and it is here that the above equations deviate from their counterparts for Newtonian fluids such as water. For the friction coefficient $\mu (h,\bar {u})$ we adopt the modified expression introduced by Edwards et al. (Reference Edwards, Viroulet, Kokelaar and Gray2017) that is valid for granular matter in the fully dynamic regime, i.e. without stopping regions

(2.4)\begin{equation} \mu(h,\bar{u}) = \tan \zeta_1 + (\tan \zeta_2 - \tan \zeta_1)\left(1 + \frac{\beta h}{\mathcal{L} (Fr+\varGamma)}\right)^{{-}1}. \end{equation}

Here, $\zeta _1$, $\zeta _2$, $\mathcal {L}$, $\beta$ and $\varGamma$ are experimental parameters obtained from independent experiments (Edwards et al. Reference Edwards, Viroulet, Johnson and Gray2021). The angles $\zeta _{1}$ and $\zeta _{2}$ denote the two inclinations delimiting the interval in which uniform granular flows are possible: for $\zeta < \zeta _{1}$ the granular sheet remains at rest, whereas for $\zeta > \zeta _{2}$ it flows downwards in an accelerated fashion. The length scale $\mathcal {L}$ has its origins in the functional form of the experimental fit to the critical angle curves that determine the friction force, as described by Pouliquen & Forterre (Reference Pouliquen and Forterre2002). It represents a characteristic flow thickness (usually ranging from $1$ to $2$ particle diameters) over which the transition between $\zeta _{1}$ and $\zeta _{2}$ takes place in the above friction law, (2.4). Its exact value depends on the properties of the grains as well as on the bed roughness (Pouliquen & Forterre Reference Pouliquen and Forterre2002; Edwards et al. Reference Edwards, Viroulet, Kokelaar and Gray2017). The parameter $\beta$ and the Froude number offset $\varGamma$ are closely related to each other and are given by (Forterre & Pouliquen Reference Forterre and Pouliquen2003)

(2.5)\begin{equation} Fr = \beta \frac{h}{h_{{stop}}(\zeta)}-\varGamma, \end{equation}

where $h_{{stop}}(\zeta )$ denotes the thickness of the deposition layer that is left on the chute (inclined at an angle $\zeta$) after a uniform flow of height $h$ has passed over it. A first version of (2.5) (without the $\varGamma$ offset) was derived by Pouliquen (Reference Pouliquen1999) by experimentally measuring the Froude number as a function of the ratio $h/h_{{stop}}$ for four different systems of glass beads, revealing that this ratio is independent of the bead size and the roughness of the bed. The Froude offset $\varGamma$ was introduced by Forterre & Pouliquen (Reference Forterre and Pouliquen2003) in order to fit the experimental results concerning the mean velocity for different materials (sand and glass) into a single law, thus formulating (2.5). As stated in GDR-MiDi (2004), one must be cautious in using this relation for $h$-values close to $h_{stop}$, since the fit involved may lose its accuracy in this limit. In the present study this poses no problem, however, since for travelling waves without stopping regions one has to ensure anyhow that $h$ remains everywhere well above $h_{stop}$; so any possible ambiguities in the value of $\varGamma$ close to $h_{stop}$ can hardly affect our results. For an accurate comparison with the literature, we note that the factor $1/\sqrt {\cos \zeta }$ was absent from the definition of the Froude number $Fr$ in many of the pioneering studies.

In Edwards et al. (Reference Edwards, Viroulet, Kokelaar and Gray2017) a parameter $\beta _{*}$ was also introduced, indicating the value of the Froude number corresponding to a sheet thickness $h=h_{*}(\zeta )$ with $h_{{stop}}(\zeta )< h_{*}(\zeta )< h_{{start}}(\zeta )$, where $h_{{start}}(\zeta )$ denotes the height for which a static layer is mobilized when the inclination is increased to the angle $\zeta$. Physically, $h_{*}(\zeta )$ is the height of the thinnest (and therefore slowest) possible steady uniform flow which leaves a deposit of smaller thickness $h_{{stop}}(\zeta )$. The parameter $\beta _{*}$ is also affected by the Froude offset $\varGamma$, since it must satisfy $\beta _{*}>\beta -\varGamma$ (Edwards et al. Reference Edwards, Russell, Johnson and Gray2019). The granular flow is guaranteed to remain fully dynamic as long as $Fr(x,t) \geqslant \beta _{*}$ everywhere along the chute (Edwards et al. Reference Edwards, Viroulet, Kokelaar and Gray2017).

Another parameter that is characteristic for granular flow is the diffusive coefficient $\nu$ in the fourth term of the right-hand side of the momentum balance (2.3). We adopt the analytical expression derived by Gray & Edwards (Reference Gray and Edwards2014)

(2.6)\begin{equation} \nu = \nu(\zeta) = \frac{2 \mathcal{L} \sqrt{g} \sin\zeta }{9 \beta \sqrt{\cos\zeta}}\gamma(\zeta),\quad \textrm{where}\ \gamma(\zeta) = \frac{\tan \zeta_2 - \tan\zeta}{\tan\zeta - \tan \zeta_1}. \end{equation}

This expression is physically relevant in the interval $\zeta _1 < \zeta < \zeta _2$, decreasing from infinity at $\zeta = \zeta _1$ to zero at $\zeta = \zeta _2$. For inclination angles outside this interval the expression (2.6) becomes negative and cannot be used.

For completeness, we note that the advective term $\partial _x (h \bar {u}^{2})$ on the left-hand side of the momentum balance (2.3) is often multiplied by a shape factor $\alpha$ (of order $1$) to obtain an optimal correspondence with experiments (Saingier, Deboeuf & Lagrée Reference Saingier, Deboeuf and Lagrée2016; Lagrée et al. Reference Lagrée, Saingier, Deboeuf, Staron and Popinet2017). However, for shallow granular flows with relatively small Froude number, as in the current paper, the solutions are quite insensitive to the value of this shape factor (Saingier et al. Reference Saingier, Deboeuf and Lagrée2016; Viroulet et al. Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017) and therefore we choose to suppress it. Moreover, for shallow granular flows a certain amount of slipping motion at the bottom is almost unavoidable, even on a rough chute; this means that the difference in velocity between the lower and the upper layers is less pronounced, or in other words, that the velocity profile is not too far from plug flow. In this context, we note that in a study by Forterre (Reference Forterre2006) (see in particular figure 1 in that paper) it was found that a simple plug flow with shape factor $\alpha =1$ fitted the experimental data better than the Bagnold velocity profile, for which $\alpha$ would be $5/4$ (GDR-MiDi 2004; Börzsönyi, Hasley & Ecke Reference Börzsönyi, Hasley and Ecke2005; Gray & Edwards Reference Gray and Edwards2014).

2.2. Non-dimensional form of the equations

At this point, it is opportune to non-dimensionalize the above set of equations. We will take the height $h_0$ of the ‘uniform base flow’ as the characteristic length scale. This is the thickness of the sheet if the material released from the inlet is left to distribute itself evenly along the entire span of the chute. Practically, in the absence of obstacles, this height will be attained at some point not far from the inlet. The corresponding depth-averaged velocity $\bar {u}_0$ defines the characteristic time scale $h_0/\bar {u}_0$. So our non-dimensional variables are $\tilde {x}= x /h_{0}$ and $\tilde {t} = t / (h_{0}/\bar {u}_{0})$, and the height and depth-averaged velocity are rescaled as follows (we use tildes throughout to denote the non-dimensional quantities)

(2.7a,b)\begin{equation} \tilde{h}= \frac{h}{h_{0}},\quad \tilde{\bar{u}}=\frac{\bar{u}}{\bar{u}_{0}}. \end{equation}

The continuity equation (2.2) then becomes

(2.8)\begin{equation} \frac{\partial \tilde{h}}{\partial \tilde{t}} + \frac{\partial}{\partial \tilde{x}} (\tilde{h} \tilde{\bar{u}}) = 0 ,\end{equation}

and the momentum balance (2.3) takes the form

(2.9)\begin{equation} Fr_{0}^{2} \left(\frac{\partial}{\partial \tilde{t}} (\tilde{h} \tilde{\bar{u}}) + \frac{\partial}{\partial \tilde{x}} (\tilde{h} \tilde{\bar{u}}^{2})\right) = \tilde{h}\left(\tan\zeta-\mu(\tilde{h},\tilde{u})\right)- \tilde{h}\frac{\partial \tilde{h}}{\partial \tilde{x}} + \frac{Fr_{0}^{2}}{Re}\frac{\partial}{\partial \tilde{x}}\left(\tilde{h}^{3/2}\frac{\partial \tilde{\bar{u}}}{\partial \tilde{x}}\right), \end{equation}

where $Fr_0$ is the Froude number of the uniform base flow that is supposed to exist far upstream

(2.10)\begin{equation} Fr_{0}=\frac{\bar{u}_{0}}{\sqrt{h_{0}g\cos\zeta}}, \end{equation}

and $Re$ is given by the following dimensionless expression, being the granular analogue of the Reynolds number (introduced by Gray & Edwards (Reference Gray and Edwards2014) on the basis of the so-called $\mu (I)$-rheology), again defined in terms of the base flow conditions

(2.11)\begin{equation} Re=\frac{\bar{u}_{0}\sqrt{h_{0}}}{\nu(\zeta)}. \end{equation}

For consistency, also the friction coefficient $\mu$ is now expressed as a function of the new dimensionless parameters

(2.12)\begin{equation} \mu(\tilde{h},\tilde{\bar{u}}) = \tan \zeta_1 + (\tan \zeta_2 - \tan \zeta_1)\left(1 + \frac{\beta h_{0}\tilde{h}}{\mathcal{L} \left( \varGamma+Fr_{0} ( \tilde{\bar{u}}/\sqrt{\tilde{h}} ) \right)}\right)^{{-}1}.\end{equation}

Here, we note that the height $h_0$ and the velocity $\bar {u}_0$ appearing in the above expressions, are in fact related. For the uniform base flow all derivatives with respect to $x$ and $t$ vanish, and the momentum equation (2.3) reduces to a balance between gravity and friction, $\tan \zeta = \mu (h,\bar {u})$, or equivalently (in dimensionless parameters) $\tan \zeta = \mu (\tilde {h},\tilde {\bar {u}})$. Inserting $\tilde {h}=1$ and $\tilde {\bar {u}}=1$ in the latter equality, using the expression (2.12) for $\mu (\tilde {h},\tilde {\bar {u}})$ (and (2.10) for the Froude number $Fr_0$ appearing there), we get the following relation between $\bar {u}_0$ and $h_0$:

(2.13)\begin{equation} \bar{u}_0 = \frac{\beta \sqrt{g \cos \zeta}}{\mathcal{L} \gamma(\zeta)}\,h_0^{3/2} - \varGamma \sqrt{g \cos \zeta}\,h_0^{1/2}. \end{equation}

With this, the Froude number $Fr_0$ and the granular Reynolds number $Re$ can be expressed as functions of $h_0$ alone, and are in fact seen to be related to each other (Kanellopoulos Reference Kanellopoulos2021)

(2.14a,b)\begin{equation} Fr_0 = \frac{\beta h_0}{\mathcal{L} \gamma(\zeta)}-\varGamma,\quad Re =\frac{\mathcal{L} \gamma(\zeta) \sqrt{g \cos \zeta}}{\beta \nu(\zeta)}\,Fr_0(Fr_0+\varGamma). \end{equation}

The above expression for $Fr_0$, which is fully compatible with (2.5) for $h=h_0$ and $h_{stop}(\zeta ) = \mathcal {L} \gamma (\zeta )$, shows that the Froude number of the base flow is an increasing function of $h_0$. That is, for fixed $\zeta$, thicker base flows correspond to higher Froude numbers.

2.3. Travelling waves: reducing the problem to an ODE

The periodic train of granular roll waves is a travelling waveform, which propagates without change of shape and at a constant shock speed $\tilde {c}$, and therefore its height and depth-averaged velocity will have the general form $\tilde {h}(\tilde {x},\tilde {t}) = \tilde {h}(\tilde {x}-\tilde {c}\tilde {t}) = \tilde {h}(\tilde {\xi })$ and $\tilde {\bar {u}}(\tilde {x},\tilde {t}) = \tilde {\bar {u}}(\tilde {x}-\tilde {c}\tilde {t}) = \tilde {\bar {u}}(\tilde {\xi })$. Inserting this form in the governing equations (2.8)–(2.9), these become ODEs. Specifically, the dimensionless mass balance equation (2.8) takes the form

(2.15)\begin{equation} -\tilde{c}\tilde{h} + \tilde{h}\tilde{\bar{u}} ={-}K, \end{equation}

where the non-negative integration constant $K=\tilde {h}(\tilde {c}-\tilde {\bar {u}})$ corresponds to the dimensionless flux of material in the co-moving frame per unit width of the channel. Given the fact that for the base flow, which will be established at some distance from the inlet (depending among other things on the inclination angle $\zeta$), we have $\tilde {h}=1$ and $\tilde {\bar {u}}=1$, this constant is equal to $K=\tilde {c}-1$ along the whole chute.

Now, by solving (2.15) with respect to $\tilde {\bar {u}}$, we find the depth-averaged velocity $\tilde {\bar {u}} = \tilde {c}-K\tilde {h}^{-1} = \tilde {c}-(\tilde {c}-1)\tilde {h}^{-1}$ and its derivatives $\tilde {\bar {u}}_{\tilde {\xi }},\tilde {\bar {u}}_{\tilde {\xi }\tilde {\xi }}$ in terms of $\tilde {h}$. Substituting these expressions into the non-dimensional momentum balance (2.9), the velocity is eliminated and we obtain the following second-order ODE for $\tilde {h}(\tilde {\xi })$:

(2.16)\begin{equation} \frac{{\rm d}^{2} \tilde{h}}{{\rm d} \tilde{\xi}^{2}} - \frac{1}{2 \tilde{h}} \left(\frac{{\rm d} \tilde{h}}{{\rm d} \tilde{\xi}} \right)^{2} + \frac{ Re \tilde{h}^{3/2} }{ Fr_0^{2} (\tilde{c}-1) } \left[ \left( \frac{ Fr_0^{2} (\tilde{c}-1)^{2} }{\tilde{h}^{3}} -1 \right) \frac{{\rm d} \tilde{h}}{{\rm d} \tilde{\xi}} + \tan \zeta - \mu(\tilde{h}) \right]=0, \end{equation}

where also the friction coefficient $\mu$ is now expressed exclusively in terms of $\tilde {h}$

(2.17)\begin{equation} \mu(\tilde{h}) = \tan \zeta_1 + \frac{\tan \zeta_2 - \tan \zeta_1}{1 + \displaystyle\frac{\gamma(\zeta)\tilde{h}(Fr_0+\varGamma)}{\varGamma+Fr_0(\tilde{c} \tilde{h}-\tilde{c}+1) \tilde{h}^{{-}3/2}}}. \end{equation}

Equation (2.16) first appeared in a study by Gray & Edwards (Reference Gray and Edwards2014) on granular roll waves and also featured in Razis et al. (Reference Razis, Kanellopoulos and van der Weele2019), where the transition from granular monoclinal waves to roll waves was discussed. The distinguishing difference here is that the friction coefficient (2.17) is enriched by the inclusion of the Froude offset parameter $\varGamma$ (Edwards et al. Reference Edwards, Viroulet, Kokelaar and Gray2017; Kanellopoulos Reference Kanellopoulos2021).

Roll waves are nonlinear, finite-amplitude travelling waveforms that emerge as an instability to the uniform flow. Close to their birth, when their amplitude is still small, these waves are well described by a linear approximation of the governing equations (2.8)–(2.9) around the uniform solution $(\tilde {h}_0,\tilde {\bar {u}}_0) = (1,\tilde {\bar {u}}_0)$. This linear analysis yields the wave speed $\tilde {c} = 1 \pm 1/Fr_0 + \textrm {h.o.t.}$, (where h.o.t. stands for higher order terms), or in dimensional form $c = u_0 \pm \sqrt {h_0 g \cos \zeta } + \textrm {h.o.t.}$, in which one may recognize the familiar expression for the velocity (in the laboratory frame) of ‘gravity waves’ propagating on the surface of a fluid moving at speed $u_0$. As the roll waves gain height, as a result of the aforementioned instability, the nonlinear terms become increasingly important, in such a way as to mitigate their growth and eventually arrest it altogether, producing a fully ripened roll wave pattern. This sequence of events is in fact at the core of the present study, and will be treated in detail in the next sections.

2.4. The dynamical system formulation

Equation (2.16) can equivalently be written as two coupled first-order ODEs, forming a dynamical system

(2.18a)\begin{gather} \frac{{\rm d}\tilde{h}}{{\rm d} \tilde{\xi}} = \tilde{s} = f(\tilde{h},\tilde{s}), \end{gather}

(2.18b)\begin{gather}\frac{{\rm d} \tilde{s}}{{\rm d} \tilde{\xi}} = {\frac{\tilde{s}^{2}}{2\tilde{h}}} - \frac{Re \tilde{h}^{3/2}}{Fr_0^{2}(\tilde{c}-1) } \left[ \left( \frac{Fr_0^{2} (\tilde{c}-1)^{2}}{\tilde{h}^{3}} -1 \right) \tilde{s} + \tan \zeta - \mu(\tilde{h}) \right] = g(\tilde{h},\tilde{s}). \end{gather}

The fixed points of the above dynamical system lie at the intersections of the nullclines in the $(\tilde {h},\tilde {s})$ phase space. The nullclines are defined as the curves where the direction field ${\rm d}\tilde {s}/{\rm d}\tilde {h} = g(\tilde {h},\tilde {s})/f(\tilde {h},\tilde {s})$ attains infinite and zero slope, respectively,

(2.19a)\begin{gather} f(\tilde{h},\tilde{s}) = \tilde{s} = 0, \end{gather}

(2.19b)\begin{gather}g(\tilde{h},\tilde{s}) = {\frac{\tilde{s}^{2}}{2\tilde{h}}} - \frac{Re \tilde{h}^{3/2}}{Fr_0^{2}(\tilde{c}-1) } \left[ \left( \frac{Fr_0^{2} (\tilde{c}-1)^{2}}{\tilde{h}^{3}} -1 \right) \tilde{s} + \tan \zeta - \mu(\tilde{h}) \right] = 0. \end{gather}

In order to clearly distinguish between the two, we will call $f(\tilde {h},\tilde {s})=0$ the ‘vertical-motion nullcline’ and $g(\tilde {h},\tilde {s})=0$ the ‘horizontal-motion nullcline’.

Wherever the nullclines intersect, we have a fixed point. From (2.19a) we immediately see that any fixed point necessarily corresponds to a flat profile of the flowing sheet, with $\tilde {s}=0$ (zero slope). Substituting this (2.19a) into (2.19b), we recover the condition for uniform flow $\tan \zeta = \mu (\tilde {h})$, which we already discussed in connection with (2.13). This condition yields two different values of $\tilde {h}$, namely $\tilde {h}=1$ and $\tilde {h} =\tilde {h}_{+} < 1$; the latter is not a constant, but rather a function of the parameters $\tilde {c}$, $Fr_0$ and $\varGamma$, as can be seen from the following expression that relates these quantities:

(2.20)\begin{equation} \tilde{c} = \frac{1}{1-\tilde{h}_{+}}\,\frac{1}{Fr_0} \left((Fr_0+\varGamma)(1-\tilde{h}_{+}^{5/2})-\varGamma(1-\tilde{h}_{+}^{3/2}) \right). \end{equation}

The associated fixed points $(\tilde {h},\tilde {s})=(1,0)$ and $(\tilde {h}_{+},0)$ (see figure 2) correspond to a uniformly flowing sheet with constant thickness and velocity. The fixed point $(1,0)$ represents the familiar base flow near the inlet, whereas the second fixed point $(\tilde {h}_{+},0)$ represents a thinner uniform flow. Both these uniform flows can occur simultaneously on the chute, linked by a travelling shock front, forming a ‘monoclinal’ flood wave or an undular variation of this waveform (Razis et al. Reference Razis, Kanellopoulos and van der Weele2018, Reference Razis, Kanellopoulos and van der Weele2019; Kanellopoulos Reference Kanellopoulos2021).

Figure 2. The fixed points $(\tilde {h}_{+},0)$ and $(1,0)$ (indicated by open dots) in the $(\tilde {h},\tilde {s})$ phase plane for $Fr_0=0.7$ and $\tilde {c}=2.42457 (=\tilde {c}_{H}-0.004)$, together with the vertical-motion nullcline $\tilde {s}=0$ (blue) and the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (magenta, given by (2.26)). The black dotted line depicts the almost vertical contour of zero divergence of the vector field, which is discussed in § 3.3. The grey background arrows denote the direction field of the dynamical system and the large black arrows signify this field evaluated on the nullclines. The system parameters are as in table 1 and the inclination angle is chosen to be $\zeta =36^{\circ }$.

Table 1. Values of the system parameters as measured by Edwards et al. (Reference Edwards, Viroulet, Kokelaar and Gray2017) for a flow of carborundum particles over a bottom layer of glass beads glued on a chute. The inclination angle of the chute $\zeta$ must be chosen between $\zeta _1$ and $\zeta _2$ to get a steady, non-accelerated flow; in the present paper we take $\zeta = 36^{\circ }$ as a typical case, but for completeness we will also investigate other values on the interval ($\zeta _1,\zeta _2$), see § 3.5. For the given parameter values, the critical Froude number where the uniform flow gives way to roll waves is $Fr_{cr} = (2/3) (1- \varGamma ) = 0.40$, cf. (1.2).

To determine the local stability properties of either fixed point we consider the Jacobian matrix

(2.21)\begin{equation} J(\tilde{h},\tilde{s}) =\left.\begin{pmatrix} \dfrac{\partial f}{\partial \tilde{h}} & \dfrac{\partial f}{\partial \tilde{s}} \\ \dfrac{\partial g}{\partial \tilde{h}} & \dfrac{\partial g}{\partial \tilde{s}} \end{pmatrix}\right\vert_{(\tilde{h},\tilde{s})} = \left.\begin{pmatrix} 0 & 1 \\ \dfrac{\partial g}{\partial \tilde{h}} & \dfrac{\partial g}{\partial \tilde{s}} \end{pmatrix}\right\vert_{(\tilde{h},\tilde{s})}, \end{equation}

where the elements of the matrix must be evaluated at the fixed point in question. The corresponding eigenvalues are

(2.22)\begin{equation} \lambda_{a,b}(\tilde{h},\tilde{s}) = \frac{1}{2} \left( \textrm{Tr}\, J \pm \sqrt{(\textrm{Tr}\,J)^{2}-4\,\textrm{det} \,J} \right) =\frac{1}{2} \left( \frac{\partial g}{\partial \tilde{s}} \pm \sqrt{ \left( \frac{\partial g}{\partial \tilde{s}} \right)^{2} + 4 \frac{\partial g}{\partial \tilde{h}} }\,\right). \end{equation}

The above eigenvalues can be evaluated analytically, yielding four expressions (two for each fixed point) depending on $\tilde {c}$, $Fr_0$, $\varGamma$ and the angles $\zeta _1$, $\zeta$ and $\zeta _2$. The corresponding formulas are

(2.23)\begin{align} \lambda_{a,b}(\tilde{h},\tilde{s}) &= \frac{1}{2Fr_0^{2}(\tilde{c}-1)\tilde{h}^{3/2}} \left[ Re \tilde{h}^{3} + Fr_0^{2} (\tilde{c}-1) \left( \tilde{h}^{1/2}\tilde{s}-Re(\tilde{c}-1) \right)\right.\nonumber\\ &\quad \pm \left\{Re^{2} \tilde{h}^{6} + Fr_0^{4} (\tilde{c}-1)^{2} \left( Re^{2} (\tilde{c}-1)^{2}+4Re(\tilde{c}-1) \tilde{h}^{1/2}\tilde{s} -\tilde{h}\tilde{s}^{2}\right)\right.\nonumber\\ &\quad +2Fr_0^{2} Re (\tilde{c}-1) \tilde{h}^{3} \left(\left(4\tilde{h}^{1/2}\tilde{s}-Re(\tilde{c}-1)\right)\right.\nonumber\\ &\quad \left.\left.\left.+\tilde{h}^{1/2}\left(3\mu(\tilde{h})-3\tan\zeta+2\tilde{h}\mu'(\tilde{h})\right) \right) \right\}^{1/2} \right], \end{align}

where $\mu (\tilde {h})$ is given by (2.17). Since these expressions are rather long, we will simply evaluate their numerical values whenever we address the issue of the stability, with specific system parameters ($\varGamma, \zeta _1, \zeta _2$). For definiteness, in the rest of this study we adopt the parameter values measured experimentally by Edwards et al. (Reference Edwards, Viroulet, Kokelaar and Gray2017), concerning the flow of carborundum particles on a bed of glass beads glued to the chute. The values are given in table 1. This constitutes a realistic setting to demonstrate our results. Other sets of parameters are also possible (see e.g. table for the flowing sand studied by Rocha et al. Reference Rocha, Johnson and Gray2019), yet the changes this produces are only quantitative, while the qualitative trend is preserved.

Before we proceed, let us provide the explicit expressions of the trace and the determinant of the Jacobian matrix at the fixed point $(\tilde {h},\tilde {s})=(1,0)$, which we will use in the analysis to follow

(2.24)\begin{gather} \textrm{Tr} \,J(1,0) \equiv \left.\frac{\partial g}{\partial \tilde{s}}\right\vert_{(1,0)} = \frac{Re \left( 1 - Fr_0^{2} (\tilde{c}-1)^{2} \right)}{Fr_0^{2} (\tilde{c}-1)} , \end{gather}

(2.25)\begin{gather}\displaystyle \textrm{det} \,J(1,0) \equiv\left. -\frac{\partial g}{\partial \tilde{h}}\right\vert_{(1,0)} ={-}\frac{(\tan\zeta -\tan\zeta_1)(\tan\zeta_2-\tan\zeta) Re \left( \left( \tilde{c}-\frac{5}{2} \right) Fr_0-\varGamma \right)}{(\tan\zeta_2-\tan\zeta_1) (Fr_0+\varGamma) Fr_0^{2} (\tilde{c}-1) } . \end{gather}

For the same reason, we give the analytical expression for the horizontal-motion nullcline associated with (2.19b). This is a quadratic equation for $\tilde {s}$, with the following two roots:

(2.26)\begin{equation} \tilde{s}_{N} = \frac{Re\tilde{h}^{5/2}}{Fr_0^{2} (\tilde{c}-1)} \left\{\frac{Fr_0^{2}(\tilde{c}-1)^{2}}{\tilde{h}^{3}}-1 \pm \sqrt{ \left( \frac{Fr_0^{2}(\tilde{c}-1)^{2}}{\tilde{h}^{3}} -1 \right)^{2} + \frac{2Fr_0^{2} (\tilde{c}-1)}{Re \tilde{h}^{5/2}} \left( \tan \zeta - \mu(\tilde{h}) \right) } \right\}. \end{equation}

The nullcline in question consists of two branches, as can be seen in figure 2, where $\tilde {s}_{N}$ as given by (2.26) is represented by the magenta curves. These branches will play a central role in our analysis.

In closing, we note that the condition (2.19b) for the horizontal-motion nullcline corresponds to (2.16) without the second-order derivative ${\rm d}^{2}\tilde {h}/{\rm d}\tilde {\xi }^{2}$. This means that on the nullcline $\tilde {s}_{N}$ the second-order ODE governing the travelling waveforms reduces locally to a first-order ODE. Or equivalently, since the term ${\rm d}^{2}\tilde {h}/{\rm d}\tilde {\xi }^{2}$ stems from the diffusive term in the momentum balance (2.9), on the horizontal-motion nullcline the full problem locally reduces to its inviscid counterpart.

3.1. On the parameter range where roll waves may be expected

The local stability of both fixed points with $\varGamma =0$ has been addressed in Razis et al. (Reference Razis, Kanellopoulos and van der Weele2019). As we will see, the inclusion of the $\varGamma$ offset makes a marked quantitative difference yet it does not give rise to any qualitative changes. We will show this in detail for $\varGamma = 0.4$ (the value of table 1). In § 3.5 we will also discuss other values of $\varGamma$, corresponding to different experimental set-ups, which lead to the same conclusion.

Regarding the fixed point $(\tilde {h},\tilde {s})=(\tilde {h}_{+},0)$ given by (2.22) we find that the eigenvalues $\lambda _{a,b}(\tilde {h}_{+},0)$ are real and of opposite sign for all relevant parameter values, showing that $(\tilde {h}_{+},0)$ is always a saddle point.

The stability properties of the fixed point $(\tilde {h},\tilde {s})=(1,0)$ are more intricate. In the interval $\beta _{*}\leqslant Fr_0 < Fr_{cr}$ (where the uniform flow is stable, and where also the monoclinal flood wave is found) it is an unstable node or, in the immediate vicinity of $Fr_{cr}$, an unstable spiral (Razis et al. Reference Razis, Kanellopoulos and van der Weele2019; Kanellopoulos Reference Kanellopoulos2021). The monoclinal flood wave manifests itself as a heteroclinic connection in phase space between the two fixed points $(1,0)$ and $(\tilde {h}_{+},0)$. For $Fr_0 \geqslant Fr_{cr}$ (i.e. the regime where roll waves exist) the fixed point $(1,0)$ can be either an unstable spiral, a centre, a stable spiral or a stable node, depending on the parameter values. Figure 3, for the typical value $Fr_0=0.7$, shows the changing characteristics of the fixed point $(1,0)$ for varying values of the non-dimensional wave speed $\tilde {c}$.

Figure 3. The trajectories in the complex plane, for decreasing $\tilde {c}$ and at a fixed Froude number $Fr_0=0.7$, of the two eigenvalues $\lambda _{a,b}(1,0)$ of the fixed point $(\tilde {h},\tilde {s})=(1,0)$. In the insets above the main plot we sketch how the local stability characteristics of the fixed point $(1,0)$ change as the eigenvalues vary. The five schematic pictures below the main plot depict how, simultaneously, the orbits surrounding the fixed point $(1,0)$ interact with those around the neighbouring saddle point $(\tilde {h}_{+},0)$ denoted by ‘S’. The stable and unstable manifolds of the latter are indicated in red and blue, respectively. The letter ‘H’ stands for Hopf bifurcation (where $(1,0)$ is momentarily a centre point and gives birth to a limit cycle) and ‘hom’ denotes the critical event where the stable and unstable manifolds of the saddle connect to form a closed homoclinic loop. In the region between ‘H’ and ‘hom’ (dashed) the fixed point $(1,0)$ is surrounded by a stable limit cycle, see the middle picture below the main plot; this limit cycle corresponds to a periodic train of roll waves.

From the discussion below (2.15) we know that $\tilde {c}=1+K>1$. Physically, this means that the travelling waveforms of the granular sheet are always propagating faster than the carrying base flow. For the specific case of roll waves, we can also give an upper bound for $\tilde {c}$. Since roll waves manifest themselves as limit cycles in phase space, and since a limit cycle cannot be formed around a saddle point, it follows that the determinant of the Jacobian matrix $J(1,0)$, given by (2.25), must be positive in the roll wave regime. This means that the factor $(\tilde {c}-5/2)Fr_0-\varGamma$ in the numerator of this determinant must be negative, since all other factors in the expression are positive. This gives us the upper bound $\tilde {c}<5/2+\varGamma /Fr_0$. So roll waves must be sought for $\tilde {c}$ values within a sub-interval of

(3.1)\begin{equation} 1<\tilde{c}<\frac{5}{2}+\frac{\varGamma}{Fr_0}. \end{equation}

Now, the positive determinant alone does not ensure the existence of a stable limit cycle. Also needed is the presence of a Hopf bifurcation, where the fixed point $(1,0)$ turns from a stable spiral into an unstable one (momentarily and locally taking the form of a centre point). The Hopf bifurcation can be found by setting the trace of $J(1,0)$ equal to zero, see (2.24): $1-Fr_0^{2} (\tilde {c}-1)^{2} =0$. Thus one finds that the Hopf bifurcation takes place at the following value of $\tilde {c}$, cf. the analyses by Gray & Edwards (Reference Gray and Edwards2014) and Razis et al. (Reference Razis, Kanellopoulos and van der Weele2019):

(3.2)\begin{equation} \tilde{c}=\tilde{c}_{H}=1+\frac{1}{Fr_0}. \end{equation}

It may be noted that the equation $1-Fr_0^{2} (\tilde {c}-1)^{2} =0$ also gives a second solution $\tilde {c}=1-1/Fr_0$, yet this one is discarded since it lies outside the interval (3.1) of admissible $\tilde {c}$ values. The local stability properties of the fixed point $(1,0)$ in the neighbourhood of the Hopf bifurcation at $\tilde {c}_{H}$, while keeping $Fr_0$ fixed, show that $(1,0)$ is an unstable spiral for $\tilde {c}<\tilde {c}_{H}$ and a stable one for $\tilde {c}>\tilde {c}_{H}$. This means that a stable limit cycle (i.e. a roll wave solution) exists only for $\tilde {c}<\tilde {c}_{H}$. In fact, as we will see, they appear only in a very narrow interval of $\tilde {c}$-values directly below $\tilde {c}_H=1+1/Fr_0$ (cf. figure 3).

Before we proceed, let us comment upon a subtle point concerning the dimensionless wave speed $\tilde {c}_{H}$ in (3.2). This relation may give the erroneous impression that roll waves originating from a base flow with large $Fr_0$ have small wave speeds (and vice versa) since $\tilde {c}_{H}$ is a decreasing function of $Fr_0$. In reality, however, the opposite is true, as was shown by Razis et al. (Reference Razis, Edwards, Gray and van der Weele2014).

To get the correct picture, one has to consider the dimensional wave speed $c_{H}$, which reads as follows (with $\bar {u}_0$ given by (2.13) and using the expression $h_0 = (\mathcal {L} \gamma (\zeta )/\beta ) (Fr_0+\varGamma )$, cf. (2.14a,b)):

(3.3)\begin{equation} c_H=\tilde{c}_{H}\bar{u}_0=(1+Fr_0)\sqrt{(Fr_0+\varGamma) \frac{\mathcal{L} \gamma(\zeta) }{\beta}g \cos\zeta}, \end{equation}

which is an increasing function of $Fr_0$. This tells us that roll waves stemming from a base flow with higher $Fr_0$ (and hence higher $h_0$, cf. (2.14a,b)) have a larger amplitude, and therefore a larger shock speed. This is also observed in practice.

In the non-dimensional context, however, if we fix the value of $Fr_0$ (and hence $h_0$) one has to follow a path through $(Fr_0, \tilde {c})$ parameter space with decreasing $\tilde {c}$ to witness the birth of roll waves, and to see their amplitude gradually grow. This will be the type of path that we follow in our analysis.

3.2. The birth and evolution of the limit cycle

Figure 3, based on the analytical expressions of (2.22), depicts the trajectories of the eigenvalues – also known as the root locus – of the fixed point $(1,0)$ for decreasing values of $\tilde {c}$, while the Froude number $Fr_0$ is kept constant.

The sketches above the central plot of figure 3 show the changing character of the fixed point in question. We also depict, in the five figures below the main plot, the simultaneously changing interplay between $(1,0)$ and the saddle $(\tilde {h}_{+},0)$, here denoted by the letter ‘S’. The roll waves occur in the narrow dashed interval to the right of the imaginary axis, where $(1,0)$ has just become an unstable spiral. In this interval, which is only a small part of the section where $(1,0)$ is an unstable spiral, the fixed point $(1,0)$ happens to be surrounded by a stable limit cycle, corresponding to a periodic train of roll waves. Figure 4 shows that, in the $(Fr_0,\tilde {c})$ parameter space, the region of stable limit cycles (grey zone) is indeed remarkably narrow: the curves along which the Hopf bifurcation and the homoclinic event take place are seen to lie very close to each other.

Figure 4. The $(Fr_0,\tilde {c})$ parameter plane, with the region in which roll waves occur shaded in grey. The upper bound of this region (denoted by ‘H’) is given by (3.2) where the Hopf bifurcation takes place, cf. figure 3, and the lower bound corresponds to the homoclinic event (‘hom’). For Froude numbers below $Fr_{cr} = (2/3) (1-\varGamma ) = 0.4$ (hatched region) no roll waves exist; the corresponding critical value of $\tilde {c}$ follows from (3.1) and equals $3.5$, so the curves ‘H’ and ‘hom’ both emanate from the critical point $(Fr_0,\tilde {c})=(Fr_{cr}, 5/2 + \varGamma /Fr_{cr})=(0.4, 3.5)$. On the scale of the main figure the two curves ‘H’ and ‘hom’ seemingly coincide, therefore we also show an inset in which we magnify the section around $Fr_0=0.7$ (which is one of the typical values used in the text). The system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$.

In figure 5 we present the limit cycle for three different values of $\tilde {c}$ and for fixed Froude number $Fr_0=0.7$. In figure 5(a,b), with $\tilde {c}=\tilde {c}_{H}-10^{-4}$, the limit cycle – the blue closed curve surrounding $(1,0)$ – has just been born and hence the corresponding roll wave train (plot to the right) exhibits an almost sinusoidal behaviour with relatively small amplitude. Somewhat further away from the Hopf bifurcation, at $\tilde {c}=\tilde {c}_{H}-10^{-3}$, figure 5(c,d), the limit cycle is seen to have increased in size. The amplitude and the wavelength of the corresponding roll waves have increased accordingly, and the wave pattern starts to visibly deviate from a sinusoidal one. In figure 5(e,f), at $\tilde {c}=\tilde {c}_{H}-4\times 10^{-3}$ (again further away the Hopf bifurcation), the amplitude and wavelength have grown further and the asymmetry of the limit cycle has become very pronounced, corresponding to the characteristic roll wave pattern observed in practice.

Figure 5. Phase space portraits and corresponding wave profiles for three successive values of $\tilde {c}$ at fixed $Fr_0=0.7$, as shown in the accompanying section of the $(Fr_0, \tilde {c})$ parameter plane: (a,b,g(i)) at $\tilde {c}=\tilde {c}_{H}-10^{-4}=2.42847$, the limit cycle has just been born from $(1,0)$ and is still quite symmetric. The waves are nearly sinusoidal and of small amplitude. (c,d,g(ii)) At $\tilde {c}=\tilde {c}_{H}-10^{-3}=2.42757$, the limit cycle has become bigger and starts to develop an asymmetry. The roll waves are therefore of larger amplitude and wavelength, and their shape deviates from the harmonic. (e,f,g(iii)) At $\tilde {c}=\tilde {c}_{H}-4\times 10^{-3}=2.42457$, the limit cycle has grown further and now exhibits a marked asymmetry. A practically straight segment appears on its topmost part, which is responsible for the nearly exponential shape of the rising flanks of the roll waves. The latter are now seen to have an even larger amplitude and wavelength. The background grey arrows denote the direction field of the dynamical system. The system parameters are as in table 1 and the inclination angle is $\zeta =36^{\circ }$.

A striking feature of the limit cycle is the almost straight segment on its topmost part, which (as we shall see in the next subsection, cf. figure 6) can be attributed to the horizontal-motion nullcline $\tilde {s}_N$. This almost straight segment of the limit cycle is responsible for the nearly exponentially rising flank of the roll waves. This can be seen by writing down the mathematical expression for the observed straight segment above the fixed point $(1,0)$

(3.4)\begin{equation} \tilde{s} = \frac{{\rm d}\tilde{h}}{{\rm d}\tilde{\xi}} = \tilde{s}_1 + a (\tilde{h}-1), \end{equation}

where $a$ is the small positive slope of the segment, and $\tilde {s}_1$ denotes the $\tilde {s}$-value of the segment at $\tilde {h}=1$. The above equation comprises a first-order linear ODE for $\tilde {h}(\tilde {\xi })$, which is readily solved to give

(3.5)\begin{equation} \tilde{h}(\tilde{\xi}) = 1 - \frac{\tilde{s}_1}{a} + C {\rm e}^{a \tilde{\xi}}, \end{equation}

where $C$ is the constant of integration. This expression (3.5) is exponentially rising with $\tilde {\xi }$, demonstrating that the straight segment in the $(\tilde {h},\tilde {s})$ plane translates to an exponential portion of the wave profile $\tilde {h}(\tilde {\xi })$.

Figure 6. The path traced out by the limit cycle (solid blue curve) in the upper half-plane is seen to closely follow a path in the immediate vicinity of the two branches of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_N$ (magenta dashed curves, given by (2.26)). In fact, the vector field forces all trajectories in this region through a narrow channel around the left branch until they cross the contour of zero divergence (dotted), after which they start to drift apart. Subsequently, the limit cycle crosses the right branch of the horizontal-motion nullcline and curves downwards. The background grey arrows denote the direction field of the dynamical system; the parameters are as in table 1 and the inclination angle is $\zeta =36^{\circ }$.

One of the advantages of the dynamical systems approach presented here, is that it allows for the analytical determination of the wavelength of roll waves at the moment of their birth. Close to the Hopf bifurcation (i.e. when $\tilde {c} \simeq \tilde {c}_{H}$) the periodic wave train is almost perfectly sinusoidal, as seen in figure 5(b). Its wavelength is

(3.6)\begin{equation} \varLambda \simeq \varLambda_H = \frac{2{\rm \pi}}{|\lambda_{H}(1,0)|}, \end{equation}

where $|\lambda _{H}(1,0)|$ stands for the magnitude of the eigenvalues $\lambda _{a,b}(1,0)$ at the Hopf bifurcation. These eigenvalues are purely imaginary. For the parameters used here we have $\lambda _H(1,0) = \pm 0.506 i$. At the Hopf bifurcation itself, i.e. for $\tilde {c} = \tilde {c}_{H}$, the selected wavelength is precisely $\varLambda _H$, but at that moment the size of the limit cycle – and hence the amplitude of the roll waves – is still nil. In the case of figure 5(a,b), slightly beyond the Hopf bifurcation, we find numerically that the wavelength is $\varLambda = 12.45$; this is indeed close to the value at the exact Hopf bifurcation, $\varLambda _H = 12.40$, given by (3.6).

3.3. The influence of the nullclines on the shape of the limit cycle

Having examined the birth and evolution of the roll waves, we now focus on their characteristic shape, showing that this is dictated to a large extent by the form of the nullclines. Figure 5(e,f) shows that the roll waves are comprised of an elongated rising flank followed by a relatively abrupt shock front. The former corresponds to the upper part of the limit cycle, lying in the positive $\tilde {s}$ half-plane, while the latter is associated with the parabolic-like section in the negative half-plane. Along the upper part (the straight segment) the dynamics is slow, whereas the lower part is traversed very quickly.

Here, we will show that the overall shape of the limit cycle is heavily influenced by the nullclines of the dynamical system. The nullclines are the curves in phase space where the slope ${\rm d}\tilde {s}/{\rm d}\tilde {h}$ is either zero or infinite. As we have shown in § 2.4, the vertical-motion nullcline (where ${\rm d}\tilde {s}/{\rm d}\tilde {h}$ diverges) coincides with the horizontal axis $\tilde {s} =0$, see (2.19a). The horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (where ${\rm d}\tilde {s}/{\rm d}\tilde {h}=0$), which is given by (2.26), has a rather more intricate form and consists of two branches. In figure 6 the nullclines are shown for the specific parameter values $Fr_0=0.7$ and $\tilde {c}=\tilde {c}_{H}-0.004406$. The dotted line in this figure denotes the contour of zero divergence, i.e. the curve along which the expression

(3.7)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \begin{pmatrix} f(\tilde{h},\tilde{s}) \\ g(\tilde{h},\tilde{s}) \end{pmatrix} = \frac{\partial f(\tilde{h},\tilde{s})}{\partial \tilde{h}}+\frac{\partial g(\tilde{h},\tilde{s})}{\partial \tilde{s}} \end{equation}

becomes zero. To the left of this line the divergence (3.7) is negative, whereas to its right it is positive.

Let us first concentrate on the aforementioned nearly straight segment in the upper part of the limit cycle. This is seen to trace out a path adjacent to the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$. The limit cycle enters the upper half of the phase plane vertically, when it crosses the $\tilde {h}$-axis, but soon it has to re-adjust its orientation – always following the arrows of the vector field – when it approaches the left branch of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$.

Along this left branch, the divergence of the vector field is negative, so all nearby trajectories are squeezed together. As we see in figure 6, they converge into a narrow channel in the proximity of the horizontal-motion nullcline $\tilde {s}_{N}$, following its course until they are catapulted into the gap between the two branches and pass through the contour of zero divergence (the dotted line in figure 6). After this point, they start to drift away from each other, with some trajectories heading upwards and others going downwards. The limit cycle belongs to the latter category. The reason for this is that, for the current parameter values, the unstable manifold of the saddle in $(\tilde {h}_{+},0)$ (which embraces the limit cycle, cf. the central sketch in the lower part of figure 3) curls downwards. This will be discussed further in the next subsection.

When it starts to curve downwards, the trajectory of the limit cycle quickly intersects the right branch of the horizontal-motion nullcline. Here, the second derivative ${\rm d}^{2}\tilde {h}/{\rm d}\xi ^{2}$ momentarily vanishes (which corresponds to an inflection point of the associated wave profile) and immediately afterwards, the trajectory bends towards the South and makes its way to the negative half-plane. It crosses the $\tilde {h}$-axis vertically, so at this point $\tilde {h}(\tilde {\xi })$ attains its maximum value, corresponding to the crest of the roll wave.

The above description, in which the limit cycle has a pronounced straight segment, holds when the limit cycle approaches its maximal size and is restrained in its growth by the horizontal-motion nullcline (and by the unstable manifold emanating from the saddle point). We illustrate this in figure 7 for three successive values of the Froude number: $Fr_0=0.7$, $Fr_0=0.9$ and $Fr_0=1.1$, and corresponding values of $\tilde {c}$ close to $\tilde {c}_{hom}(Fr_0)$, where the subscript ‘hom’ refers to the homoclinic event that takes place at this value, as we will see shortly. The three limit cycles in the left column are very close to their maximal size and the role of the nullclines is seen to be profound.

Figure 7. (a,d,g) The maximal size limit cycles in phase space (blue solid curves) along with the two branches of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (magenta curves) for Froude numbers $Fr_0=0.7$, $Fr_0=0.9$ and $Fr_0=1.1$, respectively. The corresponding wave speeds are $\tilde {c}_{hom}(0.7)=\tilde {c}_{H}(0.7)-0.004406 \simeq 2.424165$, $\tilde {c}_{hom}(0.9)=\tilde {c}_{H}(0.9)-0.002679 \simeq 2.108432$ and $\tilde {c}_{hom} (1.1)=\tilde {c}_{H}(1.1)-0.00179804 \simeq 1.907293$. (b,e,h) The associated fully ripened wave profiles. (c,f,i) ‘Escape’ of the limit cycles (through the gap between the two branches of the horizontal-motion nullcline) for $\tilde {c}$-values that are just below the aforementioned critical ones, namely: $\tilde {c}= 2.424161$ for $Fr_0=0.7$, $\tilde {c}= 2.108431$ for $Fr_0=0.9$ and $\tilde {c}= 1.907292$ for $Fr_0=1.1$. The above choices of $Fr_0$ and $\tilde {c}$ are indicated in the $(Fr_0,\tilde {c})$ diagram in (j–m); the regions marked by the arrows in the first plot are given in more detail in the three magnifications to the right. The other parameters are invariably taken as in table 1 and the inclination angle is $\zeta =36^{\circ }$.

In the positive $\tilde {s}$ half-plane, the two branches of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (magenta dashed curves) together constitute the nearly straight upper border of the limit cycle. The orbit is forced to follow the left branch, where the divergence of the vector field is negative, and starts to drift away on the right hand side, where the divergence of the vector field is positive.

In the negative $\tilde {s}$ half-plane, where $\tilde {s}=\tilde {s}_{N}$ runs almost vertically, this nullcline is responsible for redirecting the trajectory of the limit cycle: at its lowest point it turns from South-West to North-West. This specific point corresponds to the spot of steepest descent in the shock front of the roll wave. (In fact, it is one of the two inflection points of the wave profile, where ${\rm d}^{2}\tilde {h}/{\rm d}\xi ^{2}=0$. These inflection points are located at the intersections of the right branch of the horizontal-motion nullcline with the limit cycle. One of these is the present point of steepest descent, while the other one is the aforementioned point close to the wave crest.)

On the other hand, the intersections of the limit cycle with the vertical-motion nullcline $\tilde {s}=0$ (where the slope of the direction field diverges) correspond to the highest and lowest point of the roll wave pattern. The nullcline in question is merely the horizontal axis. That is to say, the limit cycle (as well as any other trajectory) always crosses the $\tilde {h}$-axis vertically, towards the South for $\tilde {h}>1$ and towards the North for $\tilde {h}_{+}<\tilde {h}<1$, cf. figure 6.

Among the three cases shown in figure 7, the roll wave depicted in the top row (where $Fr_0=0.7$) has the smallest amplitude and the largest wavelength. When the Froude number increases, as in the second and the third row, where $Fr_0=0.9$ and $Fr_0=1.1$, respectively, the two fixed points of the dynamical system are seen to move further away from each other. As a result, the limit cycle follows the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ to larger positive values of $\tilde {s}$, i.e. it reaches larger slopes, which means that the rising flank of the roll wave has more space to develop itself. This is the reason for the continuing increase in amplitude observed in the second and third rows.

The depicted limit cycles are of ‘maximal size’ for the given Froude numbers: for values of $\tilde {c}$ below the critical ones of figure 7, and with the same numerical precision, the limit cycles break down. In the rightmost column of figure 7, the phase space trajectories are now seen to escape to infinity in the North-East direction, so the dynamical system ceases to give physically meaningful solutions for the Froude numbers in question. Physically this means that the profiles shown in figure 7 are the largest possible roll waves for the given $Fr_0$ values. This will be discussed in more detail in § 3.4.

As a side note, we observe (from the values of $\tilde {c}$ given in the caption of figure 7) that the difference between $\tilde {c}_{H}(Fr_0)$ and $\tilde {c}_{hom}(Fr_0)$ gradually diminishes for growing $Fr_0$. In the limit of large Froude numbers $Fr_0$, the interval of $\tilde {c}$-values for which roll waves exist becomes truly narrow, in agreement with the findings of Gray & Edwards (Reference Gray and Edwards2014).

3.4. On the maximal size of the roll waves

In this subsection we explain, from the dynamical systems perspective, why roll waves cannot become arbitrarily large for fixed parameter values and a given base flow. More specifically, we will address the question: Why does the limit cycle for a given Froude number $Fr_0$ exhibit a maximal size? To answer this, we will study what happens when the parameter $\tilde {c}$ is pushed below its threshold value. In practice, the system will never spontaneously push itself beyond the threshold, and the present discussion will illustrate why this is so: the system would then enter a regime where no bounded travelling wave solutions exist.

In the rightmost column of figure 7 we present the breakdown of the limit cycle when the value of $\tilde {c}$ is decreased below the critical value $\tilde {c}_{hom}(Fr_0)$, for each of the three cases $Fr_0=0.7$, $0.9$ and $1.1$. Responding to the change in $\tilde {c}$, the limit cycle tends to expand, but in the process it runs into the ‘escape channel’ in the gap region between the two branches of the horizontal-motion nullcline $\tilde {s}_N$. In these plots, we also see that the branches of $\tilde {s}_{N}$ reach higher positive values of $\tilde {s}$ for growing Froude number, as mentioned earlier, allowing the spiralling trajectories to expand further in that direction.

The reason behind the limit cycle's escape lies in the interplay between the stable and unstable manifolds of the saddle point $(\tilde {h}_{+},0)$. In figure 8 we illustrate this for the case of $Fr_0=0.7$, where we follow a vertical path through a small part of the $(Fr_0,\tilde {c})$ parameter space.

Figure 8. The changing interplay between the stable (red) and unstable (blue) manifolds of the saddle point $(\tilde {h}_{+},0)$ for $Fr_0=0.7$, at three values of $\tilde {c}$ as shown in the depicted section of the $(Fr_0, \tilde {c})$ parameter plane. In (b,d,f), a close-up of the area of interest is given, showing how the interplay determines the fate of the limit cycle. (a,b,g(i)) For $\tilde {c}_{hom}(0.7) < \tilde {c} < \tilde {c}_{H}(0.7)$, the unstable manifold (blue curve) is seen to dive below the stable manifold (red curve) and ends on the stable limit cycle surrounding $(1,0)$. (c,d,g(ii)) For $\tilde {c} = \tilde {c}_{hom}(0.7)$, the manifolds join precisely together and form a homoclinic orbit. The maximal size limit cycle has been indicated as a dotted curve. (e,f,g(iii)) For $\tilde {c} < \tilde {c}_{hom}(0.7)$, the unstable manifold sweeps over the stable one and flies off, whereas the stable manifold now forms a heteroclinic connection between $(1,0)$ and the saddle $(\tilde {h}_{+},0)$. The background grey arrows denote the direction field in phase space; the system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$.

In figure 8(a,b) (top row), when the value of the wave speed is $\tilde {c}_{hom}(0.7) < \tilde {c} < \tilde {c}_{H}(0.7)$, the unstable manifold of the saddle (blue line) curls itself below the saddle's stable manifold (red line) and closes in upon the stable limit cycle around $(1,0)$. Meanwhile, the stable manifold comes down from the positive $\tilde {s}$ half-plane; it describes a wide curve around the region containing the limit cycle and finally arrives in the saddle $(\tilde {h}_{+},0)$.

By contrast, when $\tilde {c} < \tilde {c}_{hom}(0.7)$, as depicted in figure 8(e,f) (third row), the stable manifold (red) begins at the unstable spiral point $(1,0)$. It spirals outwards and ends, like before, in the saddle $(\tilde {h}_{+},0)$ (establishing a heteroclinic connection). The unstable manifold (blue) now sweeps over its stable counterpart, flying off in the North-East direction, and – importantly – opening the aforementioned ‘escape channel’. To see this more clearly, in the right column of figure 8 we have expanded the relevant strip of the phase plane around $\tilde {s}=0$. Thus the limit cycle ceases to exist. All the orbits in the region of the former limit cycle (with the exception of the saddle's stable manifold itself) eventually escape to the North-East.

In the intermediate threshold case, $\tilde {c} = \tilde {c}_{hom}(0.7)$, the unstable and stable manifolds join precisely onto each other, thus forming a homoclinic loop, see figure 8(c,d) (second row). This homoclinic loop corresponds to a granular solitary wave with an unusually high amplitude and infinite wavelength (Razis et al. Reference Razis, Kanellopoulos and van der Weele2019). One cannot expect to witness this solitary wave in practice, however, since it exists only at this exact value $\tilde {c}_{hom}(Fr_0)$, which requires an unrealistic amount of fine tuning. Moreover, the homoclinic orbit is unstable to any small perturbations.

At this critical state, the limit cycle has attained its maximal size, and its upper segment coincides with the corresponding part of the homoclinic loop. In figure 8(c,d) the limit cycle is represented by the dotted blue line, and in figure 9(a) we depict it as a solid blue line. Evidently, even at this maximum size, the limit cycle is still distinctly smaller than the homoclinic loop. This can also be seen from the associated wave profiles (the roll wave train and the solitary wave, respectively) depicted in figure 9(b). Interestingly, despite their difference in amplitude, both of these waveforms have the same wave speed $\tilde {c}_{hom}(Fr_0)$. In fact, the rising flanks of the roll wave pattern (solid blue) have exactly the same form as the corresponding part of the large solitary wave (black dashed curve). The instability of the homoclinic loop implies that also the roll wave pattern becomes unstable at this critical event.

Figure 9. (a) The maximal size limit cycle (solid blue line) and the unstable homoclinic loop (black dashed line) in phase space, cf. figure 8(c,d). (b) The corresponding roll wave train of maximal amplitude (solid blue curve) and the unstable solitary waveform (black dashed curve). The rising flanks of the roll wave train have precisely the same shape as the associated part of the solitary wave. Both these waveforms travel at dimensionless wave speed $\tilde {c}_{hom}$. The solitary waveform only exists at this very specific value of $\tilde {c}_{hom}$. For the roll wave train, this same value of $\tilde {c}_{hom}$ is associated with the maximum amplitude the roll waves can attain for the given Froude number $Fr_0$. The system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$; in view of the extreme fine tuning of $\tilde {c}$ needed to create the homoclinic loop, and since the associated wave profiles are very unstable, the waves presented here have been generated using a value of $\tilde {c}$ marginally larger than $\tilde {c}_{hom}(Fr_0)$.

From the physical point of view, the breakdown of the limit cycle implies that for a specific value of $Fr_0$, the amplitude of the granular roll wave cannot grow beyond a certain size, with a well-defined wave speed value $\tilde {c}_{hom}(Fr_0)$. This is in agreement with the findings of Razis et al. (Reference Razis, Edwards, Gray and van der Weele2014), who determined the maximum amplitude of roll waves in the absence of viscous-like effects. To witness larger roll waves, one should open the inlet at the top of the chute further or increase the inclination angle $\zeta$, thereby raising the value of $Fr_0$. The point $(1,0)$ then corresponds to a higher uniform flow and, with the new set of parameter values, a higher roll wave will be formed. This is described below in § 4.

3.5. The influence of varying the system parameters $\zeta$ and $\varGamma$

Up to now, we have varied the parameters $\tilde {c}$ and $Fr_0$. It is natural to ask to what extent the above scenario (concerning the birth of the roll waves and the subsequent ripening to their maximal size) is affected by a change in the other system parameters. So here we discuss the inclination angle $\zeta$ and the Froude offset $\varGamma$. It will turn out that changes in these two parameters cause only quantitative changes; the scenario as a whole remains qualitatively the same.

As far as $\zeta$ is concerned, figure 10 shows how the thickness of the region between ‘H’ (the birth of the roll waves) and ‘hom’ (where they attain their maximal size) depends on the inclination angle. At the boundaries of the inclination interval, i.e. for $\zeta =\zeta _1$ and $\zeta =\zeta _2$, the thickness becomes zero. This stands to reason, because at $\zeta =\zeta _1$ the sheet is still stagnant and hence the waves have no opportunity to establish themselves. For $\zeta =\zeta _2$, on the other hand, the uniform base flow breaks down (since the flow is on the verge of turning into an uncontrollable avalanche), so here the waves cannot develop because they lose the substrate on which to do so. In the intermediate interval of $\zeta$ values, however, it is seen that the roll waves always follow the same sequence from birth (zero amplitude) to maximal size. The three detailed plots in figure 10 show that the thickness of their region of existence is largest for $\zeta \approx 42^{\circ }$. Physically, this means that (for an experimental set-up with the parameters of table 1) roll waves with the largest possible amplitude may be anticipated to be found around this inclination angle.

Figure 10. The three-dimensional $(Fr_0,\tilde {c},\zeta )$ parameter space, being an extension of the $(Fr_0,\tilde {c})$ plane shown in figure 4, with the region in which roll waves occur shaded in grey. The upper bound of this region (the green surface denoted by ‘H’, where the roll waves are born via a Hopf bifurcation) is given by (3.2); the lower bound ‘hom’ is where the roll waves attain their maximal size, which is determined by the homoclinic event as shown in figure 7. The three plots in the lower row depict the region in question – with the numerically calculated values of ‘hom’ – for three cuts through the three-dimensional plot at $Fr_{0}=0.7$, $0.9$, and $1.1$; note the different scaling along the $\tilde {c}$-axis in these plots. The system parameters are as in table 1.

To illustrate the influence of the parameter $\varGamma$, in figure 11 we show three stages of the sequence between ‘H’ and ‘hom’ for an entirely different set of system parameters (see table 2) corresponding to the experimental set-up used by Rocha et al. (Reference Rocha, Johnson and Gray2019), who studied the flow of sand on a chute with ballotini beads glued to the bottom. The inclination angle of the chute, which must be chosen between $\zeta _1$ and $\zeta _2$, is taken to be $\zeta = 34^{\circ }$, just as in the actual experiment. For the given parameter values, the critical Froude number that must be exceeded to get roll waves is $Fr_{cr} = (2/3) (1- \varGamma ) = 0.107$, cf. (1.2). The wave profiles in figure 11 have been generated for a base flow with $Fr_0=0.7$, which amply fulfils the condition $Fr_0>Fr_{cr}$. The values of the dimensionless wave speed in the three plots are (a) $\tilde {c} = \tilde {c}_H - 0.5 \times 10^{-4}$, (b) $\tilde {c} = \tilde {c}_H -7.0 \times 10^{-4}$ and (c) $\tilde {c} = \tilde {c}_H -24.28 \times 10^{-4}$, where $\tilde {c}_H=2.42857$ marks the birth of the roll waves. It is seen that the ripening of the roll waves from birth to maximal size follows the same trend as described in the previous subsection (see especially figure 5) for the system associated with the parameter values of table 1.

Figure 11. Roll wave profiles for the experimental set-up employed by Rocha et al. (Reference Rocha, Johnson and Gray2019), concerning the flow of sand on a chute. The system parameters of this set-up are given in table 2. The inclination angle of the chute is $\zeta = 34^{\circ }$ and the Froude number of the incoming flow is $Fr_0=0.7$, well above the critical value $Fr_{cr} = 0.107$ for this set-up. The depicted wave profiles correspond to three successive values of the wave speed: (a) $\tilde {c} = \tilde {c}_H - 0.5 \times 10^{-4}$, (b) $\tilde {c} = \tilde {c}_H -7.0 \times 10^{-4}$ and (c) $\tilde {c} = \tilde {c}_H -24.28 \times 10^{-4}$, where $\tilde {c}_H=2.42857$ marks the birth of the roll waves.

Table 2. The system parameters of the experimental set-up used by Rocha et al. (Reference Rocha, Johnson and Gray2019) for the flow of red sand over a bottom layer of ballottini beads glued on a chute.

For completeness, we have also checked other sets of system parameters, such as those corresponding to the flow of glass beads from the aforementioned study by Rocha et al. (Reference Rocha, Johnson and Gray2019) (for which $\varGamma =0$). In all cases studied, we found that the roll waves undergo the same qualitative sequence of events. So all the evidence suggests that the scenario is generic and quite insensitive to the precise values of $\zeta$ or $\varGamma$.