The kinematics of bidisperse granular roll waves

Small perturbations to a steady uniform granular chute flow can grow as the material moves downslope and develop into a series of surface waves that travel faster than the bulk flow. This roll wave instability has important implications for the mitigation of hazards due to geophysical mass flows, such as snow avalanches, debris flows and landslides, because the resulting waves tend to merge and become much deeper and more destructive than the uniform flow from which they form. Natural flows are usually highly polydisperse and their dynamics is significantly complicated by the particle size segregation that occurs within them. This study investigates the kinematics of such flows theoretically and through small-scale experiments that use a mixture of large and small glass spheres. It is shown that large particles, which segregate to the surface of the flow, are always concentrated near the crests of roll waves. There are different mechanisms for this depending on the relative speed of the waves, compared to the speed of particles at the free surface, as well as on the particle concentration. If all particles at the surface travel more slowly than the waves, the large particles become concentrated as the shock-like wavefronts pass them. This is due to a concertina-like effect in the frame of the moving wave, in which large particles move slowly backwards through the crest, but travel quickly in the troughs between the crests. If, instead, some particles on the surface travel more quickly than the wave and some move slower, then, at low concentrations, large particles can move towards the wave crest from both the forward and rearward sides. This results in isolated regions of large particles that are trapped at the crest of each wave, separated by regions where the flow is thinner and free of large particles. There is also a third regime arising when all surface particles travel faster than the waves, which has large particles present everywhere but with a sharp increase in their concentration towards the wave fronts. In all cases, the significantly enhanced large particle concentration at wave crests means that such flows in nature can be especially destructive and thus particularly hazardous.

The kinematics of bidisperse granular roll waves

Introduction
Large-scale debris flows spontaneously develop wave-like disturbances that move downstream faster than the material flow (Li et al. 1983;McArdell et al. 2003;Zanuttigh & Lamberti 2007).These waves arise from flow instabilities that cause small perturbations to grow into a series of large-amplitude roll waves as they travel along a channel (Jeffreys 1925;Dressler 1949;Needham & Merkin 1984;Balmforth & Liu 2004).These roll waves can grow further through merging events until they reach a maximum size (Chang et al. 1996;Balmforth & Mandre 2004;Razis et al. 2014).As a result the destructive power of each wave can be significantly higher than a uniform flow of the same average mass flux (Razis et al. 2014;Köhler et al. 2016), and therefore understanding the flow dynamics is of crucial importance for debris-flow hazard mitigation (Hu, Cui & Zhang 2012;Jenkins et al. 2015).
Similar roll waves occur in experimental dry granular chute flows (Savage 1979;Forterre & Pouliquen 2003), where they likewise arise from instability of uniform flows.Steady uniform flows of monodisperse grains flowing on a rough bed are stable only if the Froude number Fr (the ratio of flow speed to surface wavespeed) is smaller than a threshold Fr ≈ 0.57 (Forterre & Pouliquen 2003).Above this threshold, the flows are unstable to low-frequency perturbations below a cutoff frequency, but remain stable to higher-frequency perturbations (Forterre & Pouliquen 2003).
A stability analysis of the depth-averaged avalanche equations (e.g.Savage & Hutter 1989) with the basal friction law of Pouliquen (1999a) predicts that this instability occurs above a critical Froude number Fr = 2/3, but this model does not predict the cutoff frequency, instead predicting that the growth rate tends to a positive constant for high-frequency perturbations (Forterre & Pouliquen 2003).Both the critical Froude number and the cutoff frequency are predicted by a stability analysis of the µ(I)rheology (Forterre 2006), which is a full constitutive law for dense granular flows in which the friction µ is dependent on the inertial number I (GDR MiDi 2004; Jop, Forterre & Pouliquen 2005, 2006).Depth integration of this µ(I)-rheology results in depth-averaged longitudinal viscous stresses, derived by Forterre (2006) and Gray & Edwards (2014) using two different approaches.When added into the depth-averaged avalanche equations these viscous stresses provide a prediction for the cutoff frequency for instability, although in the formulation of Forterre (2006) quantitative agreement is obtained only with an adjustable parameter.Gray & Edwards (2014) used their equations to construct the explicit thickness and velocity profiles for granular roll waves.These roll waves cannot be modelled by the full µ(I)-rheology due to its underlying ill-posedness at high and low inertial numbers (Barker et al. 2015), but are observed in numerical solutions of a partially regularised form of the µ(I)-rheology (Barker & Gray 2017).
The monodisperse granular roll waves studied previously are an idealisation of the roll waves that occur in highly polydisperse geophysical flows.Such polydisperse granular flows have a tendency to segregate according to size.Several micro-mechanical explanations have been offered as to the cause of this segregation (Middleton 1970;Savage & Lun 1988;Gray & Thornton 2005;Hill & Tan 2014;van der Vaart et al. 2015;Jing, Kwok & Leung 2017) and, whilst the exact mechanism remains unclear (Staron & Phillips 2015), the common phenomenological effect is that large particles migrate towards the surface of a granular avalanche.This process is well captured by advection-segregation-diffusion equations (e.g.Bridgwater, Foo & Stephens 1985;Dolgunin & Ukolov 1995;Gray & Chugunov 2006;Gray & Ancey 2011;Fan et al. 2014;Gray 2018) when the underlying bulk velocity field is known.

S. Viroulet and others
In shearing granular avalanches the segregation of large grains to the surface causes an effective segregation in the streamwise direction.Larger particles are initially segregated to the flow surface, where the velocity is highest, and are then transported rapidly downstream towards the flow front.When this process is combined with frictional differences between the large and small grains there is a very rich variety of behaviour, for example the spontaneous self-channelisation of the flow and the formation of coarse-grained lateral levees (Félix & Thomas 2004;Johnson et al. 2012;Kokelaar et al. 2014) and lobate finger-like channels (Pouliquen, Delour & Savage 1997;Pouliquen & Vallance 1999;Woodhouse et al. 2012;Baker, Johnson & Gray 2016b), which increase the distances that geophysical mass flows travel.In shallow flows, this increased streamwise transport rate of large particles is captured by a depth-integrated segregation model (Gray & Kokelaar 2010a,b), derived by integrating the segregation equation of Gray & Chugunov (2006) through the flow depth.When combined with a depth-averaged avalanche model for the bulk mass and momentum -and a concentration-weighted basal friction law producing greater friction in coarse-rich regions (Pouliquen & Vallance 1999) -this depth-integrated segregation model predicts the formation of fingers and levees (Woodhouse et al. 2012).However, the system of Woodhouse et al. (2012) is mathematically ill posed, leading to numerical solutions that do not converge under grid refinement (Woodhouse et al. 2012).The equations are regularised by including a two-dimensional extension of the viscous terms of Gray & Edwards (2014) in the momentum balance (Baker, Barker & Gray 2016a), leading to a well-posed predictive model for granular fingering (Baker et al. 2016b).
In this paper the flow kinematics of bidisperse roll waves are studied.In § 2 the results of small-scale laboratory experiments are shown, which demonstrate an increased concentration of large particles at the crest of roll waves.A depth-averaged model for these bidisperse flows is presented in § 3, which is similar to the fingering model of Baker et al. (2016b).In § 4 this model is used to construct travelling-wave solutions for the roll waves, as well as solutions for the segregation kinematics in § 5.In § 6 the kinematics in more complex aperiodic roll waves are discussed.

Small-scale experiments
Experiments are carried out on a 3.3 m long chute, inclined at ζ = 29 • to the horizontal, with glass sidewalls 7.8 cm apart (figure 1).The base of the chute is roughened by attaching a single layer of glass beads, diameter 750-1000 µm, with double-sided tape.The flow consists of a bidisperse mixture of glass ballotini, consisting of small white (75-150 µm) and large green (200-250 µm) grains.These are both smaller than the beads on the bed, to produce no-slip conditions at the base, during flow, for both sets of particles (Pouliquen 1999a).An initially well-mixed sample is loaded into a hopper and released under two gates, one fixed at 3 mm above the bed and one movable to start and stop the flow.
As this latter gate is opened a green coarse-rich flow head rapidly forms as large particles are segregated to the surface and then preferentially transported to the front (Gray & Kokelaar 2010a,b).Behind the large rich front is an inversely graded region, with large particles concentrated at the top of the flow and small particles nearer the base.Near the gate this layer has uniform thickness, but approximately 1 m downstream (for the experiments shown here) it develops spatial and temporal instabilities.These small perturbations grow and merge into granular roll waves, which are free-surface waves that propagate faster than the bulk flow.Consequently, these waves catch up with the flow head, causing it to advance rapidly as each pulse reaches the front and more slowly between subsequent waves.Material is free to flow off the end of the chute, and, once the head has done this, the succeeding flow continues to develop roll waves (figure 1, supplementary movies 1 and 2 are available at https://doi.org/10.1017/jfm.2018.348)that are still merging and coarsening by the end of the channel (Razis et al. 2014).The slope angle and frictional properties of the grains are such that all material in the chute remains in motion.The resulting pulses are thus roll waves, as opposed to erosion-deposition waves (Edwards & Gray 2015), which are characterised by stationary regions between adjacent pulses, and occur at lower slope angles or with more frictional grains.
For an initial mixture consisting of 80 % small particles, 20 % large particles by volume, referred to as 80/20 from this point onwards (figure 1 and movie 1) the roll waves are easily observable by the end of the chute, with each wave crest appearing green due to higher concentrations of large particles near the free surface than in the corresponding white troughs.Roll waves still develop in identical experiments using 40 % fine and 60 % coarse material (40/60, see supplementary movie 2), but their amplitude is much smaller.Both the initial flow head and the waves propagate slightly slower for the mixture with a higher proportion of large particles.This is consistent with the µ(I)-rheology, which suggests that for two steady uniform flows of otherwise 840 S. Viroulet and others similar material properties, the smaller particles will move faster than the larger grains.This is because to leading order the downslope and normal momentum balances (GDR MiDi 2004;Gray & Edwards 2014) imply that µ(I) = tan ζ , where ζ is the slope angle, and hence that the inertial number I is equal to the same constant for both the large and small particles.The inertial number is defined as I = γ d/ √ p/ρ * , where γ is the shear rate, d is the particle diameter, p is the pressure and ρ * is the density of the grains.It follows, that if the shear rate γ is approximated by the ratio of the depth-averaged velocity ū to the flow depth h, and the pressure p is assumed to be the same for both flows, the small particle velocity ūsmall is related to that of the large particles ūlarge by ūsmall = (d large /d small )ū large .Since the large particles are bigger than the fines, d large > d small , it follows that ūsmall > ūlarge .
To investigate further, a high-speed colour camera (iPad Pro, Apple) is mounted perpendicular to the bed and used to capture images of the upper surface of the flow at 120 fps.A 1.05 m section of chute upstream from the outflow is recorded over a period of 25 s, and the central column of pixels from each image is extracted.These columns are then aligned next to each other to construct space-time plots of the flows (figure 2).The experiments are lit from directly above the chute to avoid introducing shadows, meaning the colours in figure 2 correspond directly to the colours of grains within a few grain diameters of the surface.The green lines in figure 2 represent the coarse-rich regions at the crest of waves, which are clearly visible in both of the mixtures shown.These are approximately straight lines, indicating that waves travel at a constant velocity u w .In both flows there is slight variation in the speed of different waves, and faster waves may catch up with slower ones.This leads to merging events, which are visualised as the intersection of two lines.Averaging the speed of 20 non-coarsening waves, the speed u w = 0.42 ± 0.035 ms −1 is obtained for the 80/20 mixture (figure 2a) and u w = 0.34 ± 0.02 ms −1 for the 40/60 mixture (figure 2b).Thus flows with higher proportions of large particles produce waves that, on average, propagate slower.This qualitative difference has also been observed for a wide range of compositions from 100 % small to 100 % large particles.In both of the flows shown here the average wavelength is ∼0.2 m, though this varies significantly in space and time.
The flow front, visible on the left of figure 2, travels at less than half the speed of the waves behind it (approximate front speeds are 0.18 ms −1 and 0.15 ms −1 for the 80/20 and 40/60 mixes, respectively), meaning the waves catch up with the front and cause it to advance in a series of pulses.The coarse-rich nature of the flow head is particularly apparent for the 40/60 mixture, with a thick green band travelling downslope, growing in size as large particles are continuously sheared to the front.
To understand the kinematics responsible for the formation of the coarse rich wave fronts present in both mixtures, another experiment is performed by releasing a small quantity of large red beads (300 µm) on top of existing waves that have developed in an otherwise identical experiment.These red grains are similar in diameter to the green ones, and so act as a tracer for large particles in the flow.These red tracers move more slowly than the roll waves for both 80/20 and 60/40 mixtures, as evidenced by space-time plots (figure 3).This is consistent with Dressler's description of roll waves as waves travelling faster than all fluid particles in the flow (Dressler 1949).In the experiments presented in this paper, the high concentration of large particles at the wave fronts is not caused by the individual particles accumulating there, because these particles move backwards relative to the waves (figure 3).Instead, large particles travel slowly backwards (relative to the wave) through the crest and are then rapidly transported down the leeward side of the wave.This results in a flux of large grains relative to the roll waves that is converging at a wave crest, increasing the concentration of large grains there, but diverging between waves, which creates a concertina effect in the concentration.The tracer grains dissipate as they move downslope due to diffusion and dispersion, particularly for the 80/20 mixture (figure 3a).The flow interior is examined by splitting the flow along its centreline with a 0.1 mm thick microscope cover slip at the end of the chute, allowing half of the flow to fall away and the other half to continue on an extended chute (see Barker & Gray 2017).Space-time plots filmed through this cover slip at 300 frames per second, using a Teledyne DALSA Genie HM1400 camera, allow the interior grain size distribution to be observed with minimal wall effects (figure 4).A coarse-rich front can clearly be observed for both the 80/20 and 40/60 mixtures, and is much larger for the flow with more large particles (figure 4b).Immediately behind the front is a region of mixed grains, a breaking size-segregation wave (Thornton & Gray 2008;Gray & Ancey 2009;Johnson et al. 2012;Gajjar et al. 2016).Roll wave instabilities follow behind the flow front, of amplitude ∼1 mm for the 80/20 mixture and ∼0.5 mm for the 40/60 mixture.For both mixtures large green particles are clearly concentrated at the flow surface (inverse grading).Although the higher concentration of large particles at the wave crests is robustly observed in the oblique views (figure 1) and aerial space-time plots (figures 2 and 3), it is less clear from the internal space-time plots.

A depth-averaged particle size-segregation and bulk flow model
The experimental flows of § 2 have a vertical length scale of ∼1 mm, much smaller than the downslope length scale of ∼1 m.This shallowness is now exploited by applying a depth-averaged modelling approach to the flows.In both cases a coarse-rich flow head is observed, followed by a mixing region and breaking size-segregation wave.Wave-like disturbances follow the flow front, with large green particles segregated to the surface, although these waves are much larger for higher concentrations of small particles (panel a).

3.1.
A depth-averaged model for particle size segregation A Cartesian coordinate system Oxz is defined at an angle ζ to the horizontal, with the x-axis pointing in the downslope direction and the z-axis aligned with the upward-facing normal (figure 5).A mass of granular material of thickness h(x, t) is assumed to lie between a flat, rigid base at z = 0 and free surface z = h(x, t).The concentration of small particles per unit granular volume is denoted φ, and the corresponding concentration of large particles is (1 − φ).The concentration φ is governed by an advection-segregation-diffusion equation (e.g.Bridgwater et al. 1985;Dolgunin & Ukolov 1995;Gray & Chugunov 2006;Gray 2018) where u = (u, w) are the bulk velocity components in the downslope and normal directions.The first term on the left-hand side is the time rate of change of the small particle concentration, the second and third terms describe advection by the bulk flow and the fourth term accounts for slope normal segregation.The quadratic flux in the segregation term ensures that the segregation stops when either the large or the small particles reach a pure phase and the factor q determines the strength of the segregation.More complicated flux functions are possible, for example, the asymmetry between large and small particle segregation velocities can be included by using a cubic flux (Gajjar & Gray 2014;van der Vaart et al. 2015) and the segregation rate q may depend explicitly on the grain size ratio and the shear rate (Dolgunin & Ukolov 1995;Marks, Rognon & Einav 2012;Schlick et al. 2015) or inertial number (Gray 2018).
FIGURE 5. Schematic diagram of the coordinate axes Oxz which are inclined at an angle ζ to the horizontal, so that the x-axis points downslope and the z-axis is aligned with the upward-facing normal.The granular material is assumed to have thickness h(x, t), downslope velocity u(x, z, t) and consists of large green particles lying above small white particles.The interface between the two regions is located at height η = h φ, where φ denotes the depth-averaged concentration of small particles.Roll wave disturbances move at a speed u w and are faster than the bulk flow.
The right-hand side of equation (3.1) describes the process of diffusive remixing of the grains through the depth of the flow.The diffusivity D may in general depend on the flow variables but is assumed constant here.
The segregation (3.1) can be integrated through the avalanche thickness using Leibniz's rule to interchange the order of integration and differentiation, and the condition that there is no flux of either large or small particles across the surface and base of the flow (see e.g.Gray & Kokelaar 2010a,b).This implies that where by definition the depth-averaged concentration and the depth-averaged flux of small particles are respectively.In the depth-integration process the segregation and diffusion terms disappear in (3.2), however, their physical effect is still incorporated via the integrals (3.3), because φ still evolves according to the full segregation equation (3.1).At this stage (3.2) is still exact and no approximations have been made.To turn the integro-differential equation (3.2) into a partial differential equation, The kinematics of bidisperse granular roll waves 845 Gray & Kokelaar (2010a,b) approximated the integrals (3.3) by assuming that (i) the segregation process dominated over the diffusion and (ii) that the segregation was sufficiently rapid that it could be considered to be instantaneous.As a result, the concentration φ could be approximated by a perfectly inversely graded profile, i.e. all the large particles over all the small grains.These assumptions are consistent with the internal space-time plots (figure 4), which show a relatively sharp interface between large and small particles, as well as sharp changes in the concentration as the roll waves pass by.Following, Gray & Kokelaar (2010a,b) the concentration profile is therefore assumed to be where η is the height of the sharp interface between the large and small particles as shown in figure 5.For simplicity Gray & Kokelaar (2010a,b) assumed that the velocity profile was linear with depth where ū is the depth-averaged velocity and the parameter A ∈ [0, 1] controls the degree of shear and basal slip.The advantage of assumptions (3.4) and (3.5) is that the integrals (3.3) are particularly simple and can be explicitly evaluated to give where the flux function is quadratic in the depth-averaged concentration.Assuming the thickness h is constant, the first term in (3.7) is the time rate of change of the depth-averaged concentration of small particles φ, the second term is the lateral transport of small particles due to the bulk flow and the final term is a reduction in the transport rate of fines due to velocity shear.Physically, equation (3.7) describes the process in which small particles are rapidly segregated to the bottom of the flow, where the velocity is lowest, and are therefore transported downslope slower than average.Since the depth-averaged concentration of large particles is 1 − φ, equation (3.7) may also be viewed as an equation that describes the preferential transport of large particles downslope.The physical mechanism is simple, large grains are rapidly segregated to the surface of the flow, where the velocity is greatest, and therefore move downslope faster than average.The final term in (3.7) has a quadratic flux φ(1 − φ) that is similar to the φ(1 − φ) structure in the full segregation (3.1).Here, however, the segregation is in the lateral x-direction, rather than through the depth of the flow.It is surprising that particle segregation through the depth of the flow, combined with velocity shear, effectively generates a secondary lateral segregation mechanism.This lateral segregation is, however, a very strong effect that gives rise to commonly observed features such as the formation of bouldery fronts in geophysical mass flows, segregation induced fingering instabilities and large particle rich levees (Pouliquen et al. 1997;Pouliquen & Vallance 1999;Félix & Thomas 2004;Johnson et al. 2012;Woodhouse et al. 2012;Kokelaar et al. 2014;Baker et al. 2016b).
Gray & Kokelaar's (2010a, b) derivation of the large particle transport equation (3.7) is very simple and captures the key physical effect.It is, however, possible to explore different approximations for the integrals.In their study of segregation-induced fingering, Baker et al. (2016b) derived the equivalent large particle transport equation assuming Bagnold flow (GDR MiDi 2004).Instead of the flux function G being symmetric about φ = 1/2 and having a maximum there, as in the quadratic case, for ), which is asymmetric with a maximum skewed towards lower concentrations of fines.This reflects the fact that the velocity gradient is stronger at the bottom of the avalanche where the small particles accumulate.In figure 7 of Baker et al. (2016b) a comparison is shown between the linear velocity (3.5) and the Bagnold profile.For A = 6/7 the linear velocity profile is close to that of Bagnold, and there is surprisingly little difference between the resulting flux curves.For simplicity, Baker et al. (2016b) therefore used the linear profile in their simulations, which made the code more stable at high concentrations, since it was no longer necessary to take the square root close to zero.The linear velocity profile (3.5) may therefore be thought of as an approximation to a more complex velocity profile with no slip at the base.Even more complex representations are possible.For instance, if one set of particles is considerably more frictional than the others, this may have an important feedback on the Bagnold-like velocity profile that develops in the two segregated media (Rognon et al. 2007).This can, in principle, also be built into the model, but it will become increasingly complex and harder to solve.In the interest of simplicity, in this paper the linear velocity profile is retained and A is treated as a control parameter with higher values corresponding to more highly sheared and segregated flows, with consequently greater preferential downslope transport of large particles.

A depth-averaged model for the bulk flow
To solve for the depth-averaged concentration φ, the large particle transport (3.7) is combined with a shallow-water model for the bulk thickness h and depth-averaged velocity ū.Following Gray & Edwards (2014), conservation of mass and momentum are given by the equations where g is the constant of gravitational acceleration, ν is a coefficient in the depthaveraged granular viscosity νh 1/2 /2, discussed later, and the shape factor χ = u 2 /ū 2 depends on the velocity shear profile.For simple shear (A = 1) this implies χ = 4/3, while for plug flow (A = 0) the shape factor χ = 1.For simplicity, χ is assumed to be equal to unity in this paper, as, in The kinematics of bidisperse granular roll waves 847 The source term where sgn is the sign function, represents the downslope component of gravity driving the flow, and the effective basal friction opposing the direction of motion.For a monodisperse granular material the basal friction coefficient µ b can be measured empirically, with the dynamic law of Pouliquen & Forterre (2002), where the Froude number is the ratio of the speed of the bulk flow to the speed of surface gravity waves.
The parameters µ 1 = tan ζ 1 and µ 2 = tan ζ 2 are constants, where angles ζ 1 and ζ 2 correspond to the minimum and maximum slope angles for which steady uniform flows are observed.The length scale L and dimensionless constant β depend on the granular material properties as well as the bed composition (Pouliquen 1999a;Goujon, Thomas & Dalloz-Dubrujeaud 2003).Interestingly, Gray & Edwards (2014) showed by depth averaging the µ(I)-rheology (GDR MiDi 2004; Jop et al. 2005Jop et al. , 2006)), subject to a no-slip condition at the base, the classical inviscid shallow-water-like avalanche equations (3.9)-(3.10)(with ν = 0) emerge at leading order with an effective basal friction given by (3.12).For rough beds, with no slip at the base, the friction law (3.12)may therefore be thought of as an integrated effect of the internal rheology, rather than a Coulombic basal sliding friction.The effective basal friction (3.12) is vital in determining the shape of the granular roll waves (Gray & Edwards 2014;Razis et al. 2014;Edwards & Gray 2015) as well as the critical Froude number Fr c = 2/3 for the instability (Forterre & Pouliquen 2003;Forterre 2006).Although the viscous terms in (3.9)-(3.10)are much smaller in magnitude, they are also needed in order to predict the correct cutoff frequency/wavenumber (Gray & Edwards 2014;Barker & Gray 2017) and obtain the right coarsening dynamics.
There is currently no widely accepted law for the effective basal friction of bidisperse flows on a rough bed.Clearly the frictional properties of the mixture should strongly depend on those of its individual constituents (large and small particles), as well as the relative amount of each of these constituents in the flow.The dependence on concentration is evident from the experimental results of § 2, where the bulk wave properties (amplitude and wavespeeds) depend on the composition of the mixture.In order to simplify the dynamics of the problem, this paper uses the monodisperse friction law of Pouliquen & Forterre (2002), given in (3.12).The changing composition in different experiments is reflected in the coefficients by using a weighted average of the large and small particle friction coefficients based on the initial concentration of the particles in the hopper φ0 , i.e.  1. Material parameters that will remain constant throughout this paper.
where the superscripts S, L denote the parameter values for pure small and large particles, respectively.These have been estimated for the laboratory set-up of figures 1-4, and are given in table 1, along with the other parameters that are kept constant in this paper.This basal friction law reduces to the monodisperse law when φ0 equals zero and unity.The chute and friction angles satisfy which allows for steady uniform flows of either species in a pure phase and also captures the higher effective friction of the larger grains.The final expression on the right-hand side of the momentum balance (3.10) is a viscous term, after Gray & Edwards (2014).The coefficient ν is given for a monodisperse flow by Gray & Edwards (2014) as where the constant With these definitions of the effective friction (3.12) and the coefficient in the viscosity (3.19), the mass and momentum balance equations (3.9) and (3.10) form a closed system for h and ū.This system is very similar to the depth-averaged µ(I)-rheology of Gray & Edwards (2014), but with a dependence on the constant φ0 .With h and ū determined, the large particle transport equation (3.7) can be used to solve for the evolution of the depth-averaged concentration of small particles φ.The model is uncoupled in the sense that there is no dependence of the bulk flow (3.9) and (3.10) on the local concentration φ.Motivated by the experimental observations of constant wavespeeds, steady travelling-wave solutions are now constructed to the bulk mass and momentum equations, (3.9) and (3.10).As noted previously, the transport equation (3.7) decouples and is solved a posteriori.In-plane deviatoric stresses, characterised by the viscous term in the momentum equation, are critical for predicting the correct growth rate and cutoff frequency of granular roll waves (Gray & Edwards 2014), but they are relatively small terms that do not change the essential shape of the fully developed waves.To gain greater insight the viscous terms are therefore neglected by setting ν = 0 in (3.10).This leads to a hyperbolic system of equations closely resembling the classical shallow-water theory investigated by Dressler (1949), who used a Chezy drag term and constructed periodic roll waves by piecing together discontinuous segments through suitable shock conditions.A similar approach is adopted here, but it is also shown how to construct the equivalent viscous solutions (ν > 0) in appendix A.

Solution procedure
Introducing a travelling coordinate ξ = x − u w t moving at constant speed u w and seeking steady travelling-wave solutions, the bulk governing equations reduce to where the depth-averaged velocity ū is assumed to be positive.The mass balance (3.9) integrates directly to give where the constant Q 1 , corresponding to the upstream flux of particles in the frame moving with the wave, is assumed positive to ensure that waves travel faster than the bulk flow.Substituting (4.3) into the momentum balance equation (4.2) and rearranging gives the first-order ordinary differential equation (ODE) where the Froude number dependence in the coefficient µ b (3.12) can be written in terms of h, u w and Q 1 using (3.13) and (4.3).Equation (4.4) therefore defines an autonomous ODE for the flow thickness h, although both constants u w and Q 1 are unknown at this stage.A single roll wave profile can be constructed by integrating (4.4) from the wave trough to its peak, and periodic wavetrains are then formed by applying jump conditions (see e.g.Chadwick 1976) for h and ū to join the peak of one wave to the trough of the next.For a shock moving at velocity u w , the depth-averaged mass and momentum jump conditions imply where '+' denotes quantities at the forward side of the shock and '−' at the backward side, which are assumed to represent the wave trough and peak respectively.For the travelling waves considered here, (4.5) implies that the constant Q 1 in (4.3) is the same on either side of the shock.All the velocities can therefore be eliminated in the momentum jump condition (4.6) by using (4.3) and (4.5).Neglecting the trivial root h + = h − , it follows that the thicknesses satisfy the quadratic equation Taking the positive roots to ensure that thicknesses remain positive everywhere, the wave heights on either side of the shock are related by The smooth part of the solution to (4.4) for the roll wave profile is the one in which the forward thickness is less than the backward thickness, i.e. h + < h − .Using this inequality in (4.8) and (4.9) it follows that the thickness must pass through the critical point where the denominator of the ODE (4.4) is zero, since In order to obtain smooth solutions in the neighbourhood of h * , the numerator must also be zero at this critical point.This implies a balance between the downslope component of gravity and basal friction at the critical point, where Fr * = ū * / √ gh * cos ζ is the corresponding Froude number and ū * the velocity at this critical point.From the friction law (3.12), it follows that where the constant γ is defined in (3.20), and hence that Evaluating (4.3) at the critical point and using the definition (4.10) allows u w and Q 1 to be written in terms of h * as The kinematics of bidisperse granular roll waves (4.15) Equation ( 4.3) may then be rearranged to write the depth-averaged velocity as Substituting the friction law (3.12)into the ODE (4.4) and using the definition of the Froude number and equations (4.13) and (4.16) to eliminate ū * and ū implies that the ODE becomes where Since both the numerator and denominator of (4.17) are zero at the critical point h = h * , the gradient is evaluated here using L'Hôpital's rule, which is positive when Fr * > 2/3, as required for roll waves.This is the same condition for instability of a steady uniform flow of thickness h * (Forterre & Pouliquen 2003; Gray & Edwards 2014) and, using (4.12), defines a minimum critical thickness for roll wave solutions to be possible, In fact, the gradient dh/dξ must remain positive for all values of h + < h < h − , not just at the critical point, in order to construct appropriate solutions.Examining the functional form of (4.17), it can be seen that this depends on the two expressions (4.18) and (4.19) making up parts of the numerator and denominator, respectively.
Clearly f 1 (h * ) = 0, and straightforward algebra shows that there is also second root of f 1 that is strictly less than h * when Fr * > 2/3.Let this root be denoted h min , which can be reliably found with standard root-finding algorithms.Again, simple algebra reveals that f 2 (h) > 0 for h h min , and combining these results implies that dh/dξ > 0 for all h > h min .Rearranging equation (4.16) shows that the depth-averaged velocity is only positive for Since f 1 (h zero ) > 0 it follows that h zero < h min , meaning that both ū and dh/dξ are positive for h > h min .The value h min therefore provides a lower thickness bound for roll Crosses represent the critical point on both graphs.For clarity only one period is shown for each wave, but periodic trains are formed by connecting many such profiles.
wave profiles with a given h * , and the construction procedure can be summarised as follows: For a given critical thickness h * satisfying (4.21), pick a minimum (trough) height h + in the range h min < h + < h * and use the jump condition (4.9) with (4.15) to calculate the corresponding maximum (peak) height h − .Then integrate (4.17) in ξ , starting from initial condition h(0) = h + and stopping when h = h − , using (4.20) to integrate through the critical point.Finally, the depth-averaged velocity ū(ξ ) is obtained by using (4.16).
Figure 6 shows an example family of roll wave solutions, all computed using a critical thickness h * = 0.002 m and mean concentration φ0 = 0.5.All of these waves travel at the same velocity u w , determined by (4.14), but have differing amplitudes and wavelengths depending on the forward flow thickness h + chosen.For h + close to h * the waves have very small amplitudes and short wavelengths, but these both increase as h + decreases.As the trough reaches its minimum value h + → h min , the amplitude tends to the constant h max − h min , where h max is the thickness calculated by substituting h min into the jump condition (4.7).The kinematics of bidisperse granular roll waves 853 4.2.Relation to inflow conditions The method described above provides a systematic way to construct families of roll wave solutions that all travel at the same speed.This is not, however, what is observed experimentally, where constant inflow conditions lead to different velocity waves that consequently merge and coarsen as they move downstream.To directly relate the travelling-wave solutions of § 4.1 to the disturbances that form in chute-flow experiments solutions are now considered that have the same flux in a stationary (laboratory) frame.Assuming that there is a steady uniform upstream flow of thickness h 0 and velocity ū0 , which are related by replacing ū * and h * by ū0 and h 0 in (4.13), the upstream flux is constant in space and time and given by This flux depends on the mean concentration φ0 through the coefficients β, L and γ .
A periodic train of roll waves forming downstream of this flow must have the same average flux, which for the travelling waves is given by the average flux across one wavelength Λ, Considering this average flux as a function of the critical and minimum wave thicknesses, the problem reduces to finding the pairs h * and h + such that q(h * , h + ) = q 0 .For a fixed h * , q is a monotonically increasing function of the trough thickness h + , so smaller-amplitude waves, despite having lower peak fluxes than larger waves, actually have greater average fluxes.The wave amplitude decreases to zero as h + → h * , and the mean flux tends to its maximum value for a given h * , say q max (h * ), which corresponds to a steady uniform flow of thickness h * , i.e. (4.25) In the opposite limit h + → h min the mean flux approaches a finite lower limit Since q is also found to increase with increasing critical thickness (for all h + ), equations (4.25) and (4.26) imply that the critical thickness must lie in the range where h * max satisfies q min (h * max ) = q 0 and the lower bound is required for the initial steady uniform flow to become unstable.For each h * in the range (4.27) the corresponding forward shock thickness h + can then be found iteratively by ensuring that the mean flux q of the resulting wave is equal to q 0 .This defines a new family of roll waves, each having different amplitudes but all with the same flux as a given steady uniform inflow.Figure 7 shows one such family of waves, calculated using h 0 = 2 mm and φ0 = 0.5.The individual wavespeeds are determined by h * through (4.14), meaning that each wave now travels at a different speed and larger-amplitude, longer wavelength waves travel faster.Figure 8 . Family of roll wave solutions resulting from a steady uniform inflow of thickness h 0 = 2 mm for a mean concentration φ0 = 0.5.Crosses denote the critical thickness h * , and dashed lines the theoretical absolute minimum and maximum thicknesses, h min and h max respectively, for each h * .
wavespeed, wavelength and amplitude in more detail by considering different inflow concentrations φ0 for the same h 0 = 2 mm.The maximum possible wavelength and amplitude of waves with more small particles is larger, and for a given amplitude these small-rich waves have shorter wavelengths.In general, more small particles results in faster moving waves (for a given amplitude), which is consistent with the experimental observations.Note that the pure large particles case φ0 = 0 is not shown on figure 8 because h 0 < h * min in this more frictional regime, and so the steady uniform flow cannot develop roll waves.

Travelling-wave solutions for the concentration profile
Having found the family of wave thickness and bulk velocity profiles that can form from a given steady uniform inflow, the large particle transport equation (3.7) can now be solved to find the resulting distribution of large and small particles within the wave.Switching to the moving-frame coordinate ξ and seeking steady travelling-wave solutions, (3.7) reduces to which can be integrated directly using (4.3) to give where Q 2 is the upstream flux of small particles in the frame moving with the wave.This is a non-negative constant, equalling zero only when φ = 0. Equation (5.2) is quadratic in φ (from the definition of G( φ) in (3.8)), and has two roots (5.4) For a given bulk wave, ū is a monotonically increasing function of h, determined using (4.16), and the constant Q 1 is also known explicitly via (4.15).Consequently, (5.3) and (5.4) are functions of h, with Q 2 acting as a control parameter (figure 9).To understand the different regions it is useful to consider the discriminant (5.5) For Q 2 < Q 1 , the upstream flux of large particles relative to the wave Q 1 − Q 2 is positive, and it follows that (5.6) and therefore the two roots (5.3), (5.4) are real and distinct for all values of h (figure 9, solid lines).It is also straightforward to show that, when Conversely, when Q 2 > Q 1 there is a net downstream flux of large particles relative to the wave and the discriminant (5.5) is no longer positive for all values of h, becoming zero when (5.9) The discriminant is therefore positive for Ahū > D + and Ahū < D − , meaning the two concentration roots (5.3) and (5.4) are real and distinct in these cases, for (5.11) these roots are physically admissible only for Ahū > D + .Finally, when Q 1 = Q 2 and there is no net flux of large particles relative to the wave, equations (5.3) and (5.4) reduce to (5.12) (5.13) shown by bold solid lines in figure 9.The kinematics of bidisperse granular roll waves 857 There is a qualitative change in behaviour when Ahū = Q 1 , and, substituting for (4.15) and (4.16), this bifurcation point occurs at thickness (5.14) The physical significance of the bifurcation point can be seen by considering the freesurface velocity, u s ≡ u(z = h).From the definition of the shear profile (3.5) and the fact that (5.15) and therefore surface particles at h = h b travel at precisely the velocity of the waves.
It follows that u s < u w for h < h b and u s > u w for h > h b , meaning the bifurcation point represents the divide between surface particles travelling slower or faster than the waves themselves.The relative position of h b compared to the wave thickness range h + < h < h − plays an important role in determining admissible concentration profiles φ ∈ [0, 1].There are three different cases to consider as shown in figure 10.
5.1.Case 1: h < h b everywhere If h < h b at all points in the wave, meaning all particles travel more slowly than the wave, then the only permissible root is φ = φ(1) , assuming that Q 2 < Q 1 as shown in figure 9.In this case there is a family of concentration profiles that are determined by the specific choice of Q 2 ∈ [0, Q 1 ], which can in turn be determined by evaluating (5.2) at a given concentration φ = φp ∈ [0, 1] and bulk thickness h = h p .A selection of these profiles are shown in figure 10(a) and corresponding interfaces η in figure 10(b).Each profile φ(ξ ) is continuous along the wave, but must necessarily experience a jump at the peak where the thickness h and velocity ū are discontinuous.The jump condition relating the forward ( φ+ ) and backward ( φ− ) concentrations here is (5.16) which is equivalent to saying that the value of Q 2 must be the same on either side of the discontinuity.Consequently, it is not possible to jump between different permissible solutions (corresponding to different values of Q 2 ) at any point along the wave.This class of solution is therefore referred to as 'continuous', even though the profiles are still discontinuous at the end points.Figure 10(a) shows that φ decreases as h increases, implying higher concentrations of large particles at wave crests, as observed experimentally.

Case 2: h
= h b at an internal point Next, if h = h b at an internal point in the wave, say ξ = ξ b , then surface particles move slower than the waves for ξ < ξ b and faster for ξ > ξ b .There are two different classes of possible solutions for the concentration profile in this regime.Similarly to case 1, continuous profiles can be constructed by specifying the concentration φ = φp ∈ [0, 1] at thickness h = h p h b and selecting the first root (5.3).This again corresponds to Q 2 < Q 1 , and example profiles are shown in figure 10(c,d).For h > h b , there is a region where both φ(1) and φ( 2 In panels (c) and (d) the bifurcation point h = h b is present at ξ = ξ b and shown by a black marker, and possible shocks between φ− s = 1 and φ+ s = Q 1 /(Ahū) for ξ > ξ b are marked with dash-dotted black lines.In panels (e) and ( f ) the dashed lines show the roots φ(1) (blue) and φ(2) (red) for Q 2 > Q 1 , and dash-dotted lines represent possible shock solutions between φ− s = φ(2) and φ+ s = φ(1) in the special case (5.21) when The shaded green regions show specific solutions for the region occupied by large grains when the mean flux of small particles across one wave q S is equal to that at the inflow q S 0 .All cases use the same bulk wave, calculated using h * = 2.05 mm, h + = 1.38 mm and φ0 = 0.8, but different shear parameters A = 0.2 for (a,b), A = 0.6 for (c,d) and A = 0.9 for (e,f ).
(figure 9).However, for a particular choice of Q 2 in this region neither concentration profile remains real across the full thickness range, and hence these solutions are not permissible.The boundary case Q 2 = Q 1 is important because the two roots coalesce at h = h b and interchange through φ = 1.For h > h b , (5.12) and (5.13) imply that φ(1) < 1 and φ(2) = 1, meaning that both roots are permissible for the same choice of Q 2 .This raises the possibility of shock solutions transitioning between φ(1) and φ(2) at an internal point in the wave, where the thickness and velocity remain continuous.Such a solution would automatically satisfy the shock condition (5.16) but must also The kinematics of bidisperse granular roll waves 859 satisfy the causality condition to ensure that sufficient information is propagated into the shock (see e.g.Ockendon et al. 2004).This is equivalent to the Lax entropy condition (Lax 1957), and can be formulated in terms of the characteristic lines of the transport equation (3.7) (5.17) The causality condition requires that the characteristics on either side of an internal concentration shock must travel into it, i.e. dx/dt > u w for ξ < ξ s and dx/dt < u w for ξ > ξ s , where ξ s ξ b is the assumed internal shock position.Substituting in for the two roots (5.3) and (5.4), and using (4.3), implies and it necessarily follows that φ− s = φ(2) and φ+ s = φ(1) , where φ± s are the values of φ on either side of the internal shock.Furthermore, the causality condition should also be applied at the wave end points, where there are also shocks in h and ū.In this case one of the concentration characteristics (5.17) must travel into the shock and another travel out, because there are already three of the bulk characteristics travelling in and one travelling out (see e.g.Viroulet et al. 2017).Choosing φ− = φ(1) and φ+ = φ(2) at the end points satisfies this requirement, and, due to the two roots swapping over at the bifurcation point, is also consistent with the internal shock values.'Full internal shock' solutions can therefore be constructed as where 'full' refers to the fact that the backward side of the shock fully consists of small particles.Figure 10(c,d) shows some example solutions of this type, with different possible shock positions ξ s indicated with dash-dotted lines.Note that there is again a region with a higher concentration of large particles towards the wavefront, as seen in the 'continuous' concentration profiles and experimental results.However, the pure small particle region φ = 1 behind the wavefront is qualitatively different.This will prevent transport of large particles backwards through the wave, which was determined to be the dominant mechanism from the tracer particle experiments.This shock regime represents large particles travelling downslope with the wave itself and is indicative of a breaking size-segregation wave in a non-depth-averaged framework (Thornton & Gray 2008;Gray & Ancey 2009;Johnson et al. 2012;Gajjar et al. 2016), where large particles are continuously segregated, sheared and recirculated inside the crest of the wave.
5.3.Case 3: h > h b everywhere Finally, consider the case where h > h b at all points along the wave profile, meaning all surface particles travel faster than the waves.Continuous solutions can again be constructed by picking a thickness h p and concentration φp , using (5.2) to calculate Q 2 , and then substituting into the root (5.3).Note that Q 2 is only less than Q 1 if φp < Q 1 /(Ah p ū(h p )) at this point.Specifying φp Q 1 /(Ah p ū(h p )) places the solution on a branch where Q 2 Q 1 , and figure 9 shows that both φ(1) and φ(2) are valid at the specified location.However, they may become complex if the discriminant (5.5) is zero, which occurs when Ahū = D + .To avoid such a point lying on the wave profile one can specify the concentration at the minimum wave thickness h + .In this case Q 2 Q 1 can be determined by picking the forward concentration φ+ in the range Q 1 /(Ah + ū(h + )) φ+ 1, and both roots (5.3) and (5.4) are then real, valid solutions throughout the entire wave.These are shown as dashed lines in figure 10(e,f ).
Since there are two permissible concentration profiles for the same value of Q 2 , it should also be investigated whether internal shock solutions connecting the two, similar to those described in § 5.2, are possible in this case.The causality condition again implies that an internal shock should have φ− s = φ(2) and φ+ s = φ(1) , and the end-point shock have φ− = φ(1) and φ+ = φ(2) .It is only possible to satisfy both criteria if the concentration can switch between the two roots without becoming discontinuous.Previously, the bifurcation point h = h b provided the means to achieve this, but such a point is not present in these wave profiles.However, there is precisely one pair of solutions that do coalesce at the point h = h + where the discriminant D(h + ) = 0, and this allows for the desired interchange.From (5.3) and (5.4) the concentration at this point is φ+ = (1 + Q 1 /(Ah + ū+ ))/2 and, substituting into (5.2), the corresponding value of Q 2 is (5.21) One can therefore construct solutions of the form (5.22) which have an internal shock at position ξ s ∈ [0, Λ], where the profiles φ(1) and φ(2) are calculated by substituting (5.21) into (5.3) and (5.4).These are referred to as 'partial internal shock' solutions because the backward side of the wave is only partially saturated with small particles ( φ− s < 1), in contrast to § 5.2 where φ− s = 1. Figure 10(e,f ) shows some example solutions of this type, with different possible shock positions ξ s being indicated with dash-dotted lines.As before, the concentration of large particles is significantly higher towards the wavefront, but some large grains are present at all points in the wave.

Relation to inflow conditions
For a given bulk wave profile, figure 10 illustrates the families of possible travellingwave solutions for the concentration φ.To understand which of these profiles are likely to form downstream in chute-flow experiments, the flux of small particles can be considered in an analogous way to using the total flux for the bulk in § 4.2.Assuming the waves form from a steady uniform flow of thickness h 0 , velocity ū0 and mean concentration φ0 , the inflow has constant small particle flux q S 0 = h 0 ū0 φ0 − h 0 ū0 G( φ0 ).The average flux across one travelling wave is The kinematics of bidisperse granular roll waves 861 and the waves that will be realised are those for which q S = q S 0 .For the continuous concentration profiles, such as those shown on figures 10(a) and 10(b), the profile with the correct mean small particle flux can be found by iteratively choosing Q 2 .This profile is indicated in figures 10(a) and 10(b) by the interface between the large green region and the underlying white small particle region.The same approach can also be applied when internal shock solutions are possible (either full or partial, figure 10c-f ), providing the desired small particle flux is sufficiently low.For higher values of q S 0 , however, concentration profiles require internal shocks in order to incorporate enough small particles.In this case the correct profile is selected by iteratively choosing the shock position ξ s so that q S = q S 0 (indicated by the interface between the large particle green region and the underlying white small particle region in figure 10c-f ).In the third regime (figure 10e,f ), even higher small particle fluxes may require choosing the second concentration root (5.4) across the whole wave, with appropriate iterative choice of Q 2 .This would result in an alternative type of continuous solution with higher concentrations of small particles at the wavefront, in disagreement with the experimental results and other theoretical solutions.
The three different cases in figure 10 are all computed using the same bulk profile and varying the effective shear parameter A, which controls the relative position of the bifurcation point through (5.14).The continuous solutions in case 1 correspond to low shear regimes and, for higher shear, full internal shock solutions become possible (case 2).If the degree of shear becomes even higher then these full internal shock solutions disappear but alternative partial internal shock solutions may occur (case 3).Now, figure 7 shows that each inflow condition (h 0 , φ0 ) gives rise to a one-parameter family of bulk waves with the same mean flux q 0 as the inflow, and figure 11 applies the concentration-flux matching approach described above to each wave in this family.Here, the bulk profiles are parameterised by their frequency f , which can be related to the wavespeeds and wavelengths of figure 8 using f = u w /Λ, with smaller-amplitude, slower, shorter waves corresponding to higher frequencies.The phase diagram figure 11 shows the different classes of solution that would form from a given steady uniform inflow, depending on f and the shear parameter A. Below a minimum frequency no travelling-wave solutions are possible for the bulk, but above this frequency the bulk waves may be augmented with any of the continuous concentration profiles, full internal shock concentration profiles or partial internal shock solutions for the concentration.Note that travelling-wave solutions for φ are unique for a given bulk profile.

Time-dependent numerical simulations
The full system of partial differential equations (PDEs) (3.7), (3.9) and (3.10) is now solved numerically using the shock-capturing central scheme of Kurganov & Tadmor (2000), whose semi-discrete formulation is combined with a Runge-Kutta-Chebyshev adaptive time stepper (Medovikov 1998) and weighted essentially non-oscillatory (WENO) flux limiter (detailed in Noelle 2000).Similar numerical methods have been employed for related systems of conservation laws governing granular flows, for example monodisperse roll waves (Razis et al. 2014), erosion-deposition waves (Edwards & Gray 2015) and flow over variable basal topography (Viroulet et al. 2017), as well as segregation-induced finger formation in bidisperse flows (Baker et al. 2016b).
To simulate bidisperse roll waves in an inclined chute, a numerical domain 0 x 5 m is discretised over 50 000 grid points, and initial conditions solutions that theoretically form from a steady uniform inflow of thickness h 0 = 2 mm and mean concentration φ0 = 0.8, depending on the perturbation frequency f and effective shear parameter A. Boundaries are calculated by matching both the mean total flux q and small particle flux q S to the inflow values, q 0 and q S 0 respectively.
representing steady uniform flow are enforced along its length.Note that these conditions do not capture the initial front propagation down the empty chute, but have the advantage of avoiding the degeneracy of the system as h → 0. The variables at the inflow of the chute are set to be h(0, t) = h 0 + δh 1 (t), ū(0, t) = ū0 , φ(0, t) = φ0 , (6.2a−c) which constitute the same steady uniform flow as the initial conditions but with a time-dependent perturbation to the flow thickness, designed to trigger the roll wave instability.All simulations are carried out with the same perturbation magnitude δ = h 0 /100, but different forms of the function h 1 (t) are subsequently employed.The numerics are computed using the viscous form of the momentum equation (3.10) with ν = 0.001 m 3/2 s −1 .

Periodic inflow perturbation
The inflow perturbation is initially taken to be a sinusoidal function of the form h 1 (t) = sin(2πft), where f = 2 Hz is the perturbation frequency.This periodicity is motivated by the desire to produce well-defined regular waves that can be directly related to the ODE travelling-wave solutions derived in § § 4 and 5.
Figure 12 shows the results of two numerical simulations at time t = 20 s, each computed with steady uniform thickness h 0 = 2 mm and concentration φ0 = 0.8.The resulting thickness profiles are thus identical, and the bulk waves grow as they move downslope before their amplitudes eventually saturate (as shown in movies 3 and 4), forming a periodic train of steady travelling waves.Differences in the interface profiles η(x, t) arise due to different shear parameters A being used in the two cases, with figure 12(a) showing a low shear value (A = 0.1) where the interface largely follows the wave.This corresponds to the 'continuous' concentration profiles described in § 5. Figure 12(b), on the other hand, shows a higher shear value (A = 0.5) which leads to the formation of the second class of solutions with full internal concentration shocks.Numerical simulations have also been carried out in the third regime where h > h b everywhere, leading to concentration profiles with partial internal shocks, but these are omitted here for brevity.

S. Viroulet and others
These simulations suggest that all classes of travelling-wave solutions predicted by the inviscid ODE are stable and manifest themselves in the full system of viscous PDEs. Figure 12(c,d) extends this idea further by calculating the ODE solutions that have mean total flux q = q 0 and small particle flux q S = q S 0 , and overlaying the profiles on the final downstream wave extracted from the PDE solutions.There is excellent agreement for both the thickness and concentration profiles in both regimes, indicating that solutions of the viscous PDEs are well approximated by inviscid travelling waves.Further comparisons between the inviscid and viscous solutions are described in the Appendix.
To further investigate the kinematics of these two classes of solution, a tracer region of large red particles is now introduced into the simulations, which is designed to mimic the experimental approach of § 2. To achieve this, the mass and momentum balance and transport equations, (3.9), (3.10) and (3.7) are solved as above until t = 20 s, allowing a well-developed periodic wavetrain to form.The final thickness, velocity and concentration profiles are then taken as initial conditions for a second set of simulations, where the same equations are solved for h, ū and φ (with identical inflow perturbation) but an additional transport equation is solved for a new variable, say φR .This represents the relative amount of large red tracer particles, and is governed by the same transport equation (3.7) but with initial and inflow conditions where x 0 and x 1 are the up-and downslope limits of the tracer region.For a sufficiently narrow range (x 0 , x 1 ), solving for φR thus represents the evolving concentration of a (large particle rich) red tracer region, whose boundary is determined by the interface η R = h φR .Note that this approach works because the transport equation is decoupled from the bulk, meaning that it can be solved without altering the overall wave properties.
It is also insightful to study the paths of individual tracer particles on the surface of the flow, whose downslope position is determined by dx dt = u s , (6.5) where the surface velocity u s (x, t) is calculated by substituting z = h into the shear profile (3.5). Figure 13 and supplementary movies 5 and 6 show the results of these tracer region/particle simulations for the two different shear regimes.In the low shear case (A = 0.1, figure 13a,c,e,g and movie 5) the waves always travel faster than the tracer region, consequently catching up and then passing through the red shaded section.A compression/dilation concertina effect is also apparent, with the tracer region quickly being compressed as the wavefront initially passes through and then slowly dilating as it retreats relative to the leeward side of the wave.This is confirmed by the trajectories of surface particles starting at the lateral boundaries of the tracer region x 0 and x 1 .These stay a finite distance apart from each other, but the distance decreases/increases during the wave cycle.The kinematics of the high shear regime (A = 0.5, figure 13b,d,f,h and movie 6) are completely different.In this case a red tracer region starting just behind the wavefront initially travels faster than the waves and is sheared to the front.It then proceeds to

S. Viroulet and others
The breaking wave ensures that large particles that are over-run by the advancing front are re-segregated upwards into the faster moving layers allowing the large-rich frontal region to grow in size.Small perturbations to the inversely graded upstream flow grow as they travel downslope, eventually developing into fully formed roll waves that travel faster than the bulk flow.These waves have approximately constant velocity (figure 2), although waves exist with a range of wavespeeds, leading to merging events as faster waves engulf slower moving ones.The inversely graded flow, large particle flow head, internal breaking size-segregation wave, and roll waves are directly observed though internal visualisation of the experimental flows (figure 4).
This basic roll wave instability mechanism and subsequent coarsening dynamics are understood for monodisperse flows (Forterre & Pouliquen 2003;Gray & Edwards 2014;Razis et al. 2014;Edwards & Gray 2015).In bidisperse flows, there are higher proportions of large particles at the wave fronts compared to the rest of the flow, which consequently form dark green bands in figure 2. Large red tracer particles seeded onto the flow surface travel more slowly than the waves and pass through the wave crests (figure 3), indicating that the increased large particle concentration at the wave fronts is not caused by large particles becoming trapped and accumulating in the wave crests.
Instead, it is caused by spatial variations in the large particle flux in a frame moving with the waves.This results in the layer of large particles on the surface of the flow being compressed in the streamwise direction near the wave crest, and consequently thickening this layer.This compression and thickening is less pronounced in the small particle layer at the base of the flow, meaning that the depth integrated concentration of large particles is increased in the wave crests.This mechanism is observed experimentally; as a wavefront catches up with a region of red tracer particles, these tracers become compressed in the streamwise direction, before dilating after the wave has passed (figure 3).
The formation of bidisperse roll waves and increased concentration of large particles at roll wave crests are predicted by a depth-averaged model for the flow that combines particle size segregation (Gray & Kokelaar 2010a,b) with bulk mass and momentum balance equations (Gray & Edwards 2014).Motivated by the experimental observations of constant wavespeeds, inviscid travelling-wave solutions are constructed for the bulk thickness and velocity profiles by switching to a moving reference frame.Equating the mean flux of material across one wave with that at the inflow, the model predicts that a given steady uniform flow can develop into a family of different steadily travelling waves, each with different wavespeeds, amplitudes and wavelengths (figure 7).Fixing any one of these properties determines a unique bulk wave profile.Larger-amplitude waves are predicted to travel faster and have longer wavelength (figure 8), which is in agreement with experimental observations of larger-amplitude waves travelling faster than (and consequently merging with) smaller waves.For a given thickness of inflow, mixtures with higher proportions of small particles are predicted to produce waves that travel faster and reach higher amplitudes (figure 8), again consistent with experimental measurements (figure 2).
This agreement occurs despite the fact that frictional differences between large and small particles are accounted for in the model using a basal friction law that depends only on the mean composition of the mixture.The agreement suggests that any increase in the local basal friction caused by increased large particle concentration at wave crests is not central to the formation of roll waves.This is despite such segregation-mobility feedback playing an essential role in other bidisperse flows such as segregation-induced finger formation (Pouliquen et al. 1997 The kinematics of bidisperse granular roll waves 869 2012; Baker et al. 2016b).A natural extension to our modelling would be to account for this feedback, by using the local depth-averaged concentration φ(x, t) in our friction law (3.13)-(3.16).This coupling adds significant mathematical complexity to the model (since the flow dynamics now depend on the predicted concentration profiles throughout the wave), but is simple to implement in numerical codes, and is likely necessary for quantitative prediction of the dynamics, as well as the kinematics, of bidisperse roll waves.
Having obtained these bulk travelling-wave profiles, travelling-wave solutions can then be constructed for the relative concentration of particles throughout a wave.Three distinct classes of solution are found which all have higher concentrations of large particles towards the front but have key qualitative differences (figure 10).The first 'continuous' class of solution has an effective interface between layers of large and small particles that mostly follows the thickness profile, whereas the second 'full internal shock' solutions have a region of pure small particles at the rear of the wave transitioning rapidly to a coarser-rich zone at the front across a shock.The third 'partial internal shock' solutions are similar to case 2, except that some large particles are present at all points along the wave.By matching the average flux of small particles across a wave to the inflow conditions, one can predict the types of solution that may develop in chute-flow experiments (figure 11).It is found that the continuous concentration profiles typically correspond to low amounts of shear in the flow, the full internal shock solutions occur at higher degrees of shear and the partial internal shock solutions appear for the highest shear values.Fully time-dependent numerical simulations confirm the existence and stability of the different types of concentration profile solutions (figure 12 and movies 3 and 4), which spontaneously form by perturbing a steady uniform inflow by a known frequency and allowing these perturbations to grow into a well-defined wavetrain.The resulting waves agree extremely well with the theoretical travelling waves, both in terms of the bulk and concentration profiles.The kinematics are also investigated numerically by simulating the evolution of a region of tracer particles, as well as individual grains (figures 13 and 14 and movies 5-8).In the low shear continuous regime, it can be seen that the particles move backwards relative to oncoming waves, and the tracer region dilates and compresses in the same concertina effect observed in the experiments.The kinematics are very different in the higher shear partial and full internal shock regimes, where surface tracer particles either move faster than the waves or move with the velocity of the wavefronts.At the moment there is no experimental evidence of this type of behaviour, but it may well occur in highly sheared geophysical flows.
Irrespective of which type of kinematics dominates, hazardous geophysical mass flows spontaneously develop waves that move significantly faster than depth-averaged flow, have wave peak heights that are significantly higher than an equivalent steady uniform flow and have high concentrations of large particles at their wave crests.All of these factors strongly enhance the impact pressures that such flows can exert on structures in their path, making these flows significantly more destructive than current design criteria may allow for.It is therefore anticipated that modelling roll waves and surges in hazard assessments will become much more important in future.smooths the shape of the wave and leads to a less sharp shock, but the differences are otherwise minimal.

FIGURE 1 .
FIGURE 1. Oblique view showing a bidisperse roll wave experiment in a chute inclined at ζ = 29 • .An initially homogeneous mixture consisting of 80 % small glass ballotini (white, 75-150 µm diameter), 20 % large ballotini (green, 200-250 µm) flows steadily through a hopper gate raised 3 mm from the bed.The steady flow near the inflow develops wave-like pulses further downstream, with large green particles rising to the surface of the flow and accumulating in higher concentrations at the fronts of the waves and at lower concentrations in the troughs.

FIGURE 2 .
FIGURE 2. Aerial space-time plots for bidisperse mixtures consisting of (a) 80 % small white ballotini (75-150 µm) 20 % large green ballotini (200-250 µm) and (b) 40 % small 60 % large particles.The plots are obtained by using a high-speed camera to capture a section of the chute from the outflow to approximately 1.05 m upstream.The central columns of pixels from each image are then aligned next to each other to give the resulting figures.Straight lines represent the green coarse-rich crests of roll waves travelling at roughly constant velocities, faster than the initial flow head.Merging events correspond to locations where two lines intersect.Distances are given relative to the outflow.

FIGURE 3 .
FIGURE 3. Aerial space-time plots for mixtures consisting of (a) 80 % small white ballotini (75-150 µm), 20 % large green ballotini (200-250 µm) and (b) 40 % small, 60 % large, particles with red large particles (300 µm) released on the surface of the waves.The red particles travel backwards relative to the waves, and are compressed together when passing through the wave crests before being stretched out at the back of each wave. https://doi.org/10.1017/jfm.2018.

FIGURE 4 .
FIGURE 4. Internal space-time plots for bidisperse mixtures consisting of (a) 80 % small white ballotini (75-150 µm), 20 % large green ballotini (200-250 µm) and (b) 40 % small, 60 % large particles.The plots are obtained by using a splitter plate and high-speed camera at the outflow to capture a vertical slice of the flow interior at each moment in time.Images are then aligned next to each other to give the resulting figures.In both cases a coarse-rich flow head is observed, followed by a mixing region and breaking size-segregation wave.Wave-like disturbances follow the flow front, with large green particles segregated to the surface, although these waves are much larger for higher concentrations of small particles (panel a). https://doi.org/10.1017/jfm.2018.
.20) is positive for slope angles ζ 1 < ζ < ζ 2 where monodisperse steady uniform flows are possible.As in the basal friction coefficient, the definition (3.19) is extended to depend on the mean concentration by substituting (3.14)-(3.17)into (3.19) and (3.20), which ensures ν > 0 and the equations are well posed for all φ0 and angles in the range (3.18).

FIGURE 6 .
FIGURE 6.An example family of roll waves all having the same critical thickness h * = 0.002 m and mean concentration φ0 = 0.5 but varying minimum thicknesses h + .All of the waves shown travel at the same velocity u w .(a) Shows the thickness profile, with dotted lines denoting the absolute minimum and maximum thicknesses, h min and h max respectively, for this value of h * .(b) Shows the corresponding depth-averaged velocity ū.Crosses represent the critical point on both graphs.For clarity only one period is shown for each wave, but periodic trains are formed by connecting many such profiles. https://doi.org/10.1017/jfm.2018.

FIGURE 8 .
FIGURE 8.The relationship between the amplitude and (a) wavespeed u w and (b) wavelength Λ for families of waves resulting from different steady uniform inflows, all with h 0 = 2 mm but different mean concentrations φ0 .

FIGURE 9 .
FIGURE 9. Plots of the depth-averaged concentration profiles (5.3) (blue) and (5.4) (red) as functions of flow thickness h for different values of the constant Q 2 .Solid lines denote cases where Q 2 < Q 1 and dashed lines where Q 2 > Q 1 .Bold solid lines represent the boundary Q 2 = Q 1 , and the black marker shows the bifurcation point h = h b , determined by (5.14).All profiles are calculated using h * = 2 mm, φ0 = 0.5 and A = 0.25.Note permissible concentrations lie in the range φ ∈ [0, 1].

FIGURE 10 .
FIGURE 10.Permissible travelling-wave solutions for the concentration profile φ(ξ ) (a,c,e) and corresponding interface height η(ξ ) (b,d,f ) in the three cases described in § § 5.1, 5.2 and 5.3.Solid black lines represent φ = 1 on the left-hand plots and the flow thickness z = h on the right-hand plots, and solid blue lines represent the root φ(1) forQ 2 < Q 1 .Bold blue lines in (c-f ) show φ(1) when Q 2 = Q 1 .In panels (c) and (d) the bifurcation point h = h b is present at ξ = ξ b and shown by a black marker, and possible shocks between φ− s = 1 and φ+ s = Q 1 /(Ahū) for ξ > ξ b are marked with dash-dotted black lines.In panels (e) and ( f ) the dashed lines show the roots φ(1) (blue) and φ(2) (red) for Q 2 > Q 1 , and dash-dotted lines represent possible shock solutions between φ− s = φ(2) and φ+ s = φ(1) in the special case (5.21) whenQ 2 = Q * 2 .The shaded green regions show specific solutions for the region occupied by large grains when the mean flux of small particles across one wave q S is equal to that at the inflow q S 0 .All cases use the same bulk wave, calculated using h * = 2.05 mm, h + = 1.38 mm and φ0 = 0.8, but different shear parameters A = 0.2 for (a,b), A = 0.6 for (c,d) and A = 0.9 for (e,f ). https://doi.org/10.1017/jfm.2018. https://doi.org/10.1017/jfm.2018.

FIGURE 11 .
FIGURE 11.Phase diagram showing the different classes of inviscid travelling-wave (TW)solutions that theoretically form from a steady uniform inflow of thickness h 0 = 2 mm and mean concentration φ0 = 0.8, depending on the perturbation frequency f and effective shear parameter A. Boundaries are calculated by matching both the mean total flux q and small particle flux q S to the inflow values, q 0 and q S 0 respectively.

FIGURE 12 .
FIGURE 12. Numerical simulations showing the thickness h (thick solid lines) and interface η (thin solid lines) profiles at time t = 20 s for flows beginning from h 0 = 2 mm and φ0 = 0.8, subject to a periodic inflow perturbation of frequency f = 2 Hz.The green shaded region therefore corresponds to the region with large grains and the white region underneath it to small grains.The shear parameter in (a) is A = 0.1 and in (b) is A = 0.5.Panels (c) and (d) show close ups of the final wave extracted from the PDE solutions (black lines), with the predictions from the inviscid ODE solutions superimposed on top (red lines).Supplementary movies 3 and 4 are available online showing the time evolution of both states.