Bulbous head formation in bidisperse shallow granular flow over an inclined plane

Rapid shallow granular flows over inclined planes are often seen in nature in the form of avalanches, landslides and pyroclastic flows. In these situations the flow develops an inversely graded (large at the top) particle-size distribution perpendicular to the plane. As the surface velocity of such flows is larger than the mean velocity, the larger material is transported to the flow front. This causes size segregation in the downstream direction, resulting in a flow front composed of large particles. Since the large particles are often more frictional than the small, the mobility of the flow front is reduced, resulting in a so-called bulbous head. This study focuses on the formation and evolution of this bulbous head, which we show to emerge in both a depth-averaged continuum framework and discrete particle simulations. Furthermore, our numerical solutions of the continuum model converge to a travelling wave solution, which allows for a very efficient computation of the long-time behaviour of the flow. We use small-scale periodic discrete particle simulations to calibrate (close) our continuum framework, and validate the simple one-dimensional (1-D) model with full-scale 3-D discrete particle simulations. The comparison shows that there are conditions under which the model works surprisingly well given the strong approximations made; for example, instantaneous vertical segregation.


Introduction
For accurate zonation and risk assessment, it is important to accurately predict the distance to which hazardous particulate natural flows such as debris flows, pyroclastic flows or snow slab avalanches might travel (Dalbey et al. 2008). These flows are heavily influenced not only by the basal topography, but also by the mixture properties, such as water saturation and the size distribution of the particulate components (e.g. snow, pumice, rock, sand). Large-scale experiments performed by FIGURE 1. Experimental debris-flow deposit at the United States Geological Survey (USGS) flow flume, Oregon, USA, August 2008 (Logan & Iverson 2013). The back of the flow has a constant height, while the front shows evidence of a bulbous head; the flow is higher near the front than at the back of the flow. Picture courtesy of USGS.
behind, can be unstable and leads to a lateral flow instability, commonly called the fingering instability (Pouliquen & Vallance 1999). This instability leads to the formation of levees, that channelise the flow, and therefore can have a significant influence on the run-out distance (Iverson 1997). When the slope is steep enough that the flow of large particles does not arrest, the downstream segregation structure leads to an accumulation of material in the front, leading to the so-called bulbous head, see figure 1. This bulbous head is seen both in geological deposits (Kokelaar et al. 2014;Takahashi 2014) and large-scale experiments (Iverson et al. 2010;Johnson et al. 2012;Haas et al. 2015). The difference in friction in the front of these flows is not only due to segregation effects, but also due to differences in water saturation: all experimental flows mentioned above are partially water saturated, for which it is known that the front contains less water than the tail (Takahashi 2014).
The characterisation of the behaviour of water-saturated granular flows is a topic of ongoing study. Some early experimental studies were done by Davies (1990), who focused on the characterisation of granular 'waves' on moving belts, and Iverson & LaHusen (1993), who studied the main origin of the frictional behaviour of water-saturated granular flows on large-scale inclined planes. They were followed by many publications, see e.g. Haas et al. (2015) and references therein. At the same time, many models for (partially) water-saturated monodisperse granular flows were developed. Examples include the early model of Iverson (1997), who developed a depth-averaged model based on mixture theory, and the model of Berzi & Jenkins (2009), which uses dense kinetic theory for the granular phase. For a recent overview, we refer the interested reader to Delannay et al. (2017). To fully understand large particulate geophysical flows, it is important to understand both the influence of segregation and the influence of water saturation in these flows; however, here we are interested in a leading-order model that does not explicitly take into account the presence of water. In § 3 we show that even with this strong simplification, we are able to capture much of the dynamics of these flows.
Generally, the debris-flow models discussed above do not take into account segregation effects inside granular chute flows; however, particles in such flows typically have a size distribution that spans many orders of magnitude and will therefore segregate into layers of different sizes (Dunning & Armitage 2011). Here, we will investigate the influence of particle segregation on dry granular flows. We focus on bidisperse granular materials, i.e. materials consisting of just two constituents of different sizes. This greatly simplifies the problem, while maintaining the leading-order segregation behaviour (Pouliquen & Vallance 1999;Goujon, Dalloz-Dubrujeaud & Thomas 2007;Gray & Ancey 2009). This is a good model system to investigate the fundamental dynamics in the more complex wet geophysical problem. 266 I. F. C. Denissen and others Discrete particle simulations, which model the motion of each grain, are a logical choice for studying and understanding particulate flows (Campbell, Cleary & Hopkins 1995;Luding 2008). However, large geophysical flows often consist of several cubic kilometres of material; combined with the average grain size of less than a cubic centimetre (Dunning & Armitage 2011), this means that there can be of the order O(10 15 ) particles in such a flow. Previously, discrete particle simulations have been shown to be very accurate for simulating segregating laboratory-scale chute flows (Thornton et al. 2012b); but, are currently limited to O(10 8 ) particles, which is seven orders of magnitude less than needed for the simulation of geophysical flows.
Alternatively, to overcome the limits of computation time, one can model the flow as a continuum, expressing it in terms of the height, velocity and small/large particle concentration of the flow at any given position. Geophysical flows are often modelled using single-phase shallow layer models, that ignore effects due to segregation (Savage & Hutter 1989;Gray, Tai & Noelle 2003;. These models start with the general incompressible mass and momentum balances, which are then depth averaged to reduce the complexity of the model. A source term is introduced to prescribe the difference between the driving force (gravitational acceleration along the slope) and the resistive force (such as basal friction).
To account for effects of segregation in granular chute flows, an advectiondiffusion-like model was developed by Bridgwater, Foo & Stephens (1985) and later Savage & Lun (1988) derived an equivalent model using statistical arguments. Subsequently, models based on mixture theory have been developed, which led to the same advection-diffusion structure, e.g. Gray & Thornton (2005), Thornton, Gray & Hogg (2006), Gray & Chugunov (2006), Fan & Hill (2011), Marks, Rognon & Einav (2012), Tunuguntla, Bokhove & Thornton (2014). For the interested reader, Tunuguntla, Weinhart & Thornton (2016b) gives an overview of segregation models, and uses a unified notation to compare and contrast them. Moreover, an extension to include segregation based on gradients in granular temperature has been given by Hill & Fan (2016). Note that this model reduces to the model of Gray & Chugunov (2006) when kinetic stress gradients are negligible. Alternatively, one can use the more complicated models based on kinetic theory, e.g. Larcher & Jenkins (2013) and Larcher & Jenkins (2015), which use different assumptions, for example that collisions between particles are instantaneous, binary and uncorrelated.
While the aforementioned segregation models are useful for looking at evolving segregation profiles, they can be further simplified by looking at either only the mean diameter at a certain position (x, y) (Takahashi et al. 1992), or only some depth-averaged quantities (Gray & Kokelaar 2010a). Here, we concentrate on the depth-averaged concentration of each constituent, which requires additional assumptions on the segregation behaviour. Gray & Kokelaar (2010a) assumed that the flow rapidly segregates into fully separated inversely graded layers; whereas, Edwards & Vriend (2016) used a more complicated three layer approach. It is noteworthy that, under Gray & Kokelaar (2010a) assumption of instantaneous complete vertical segregation, all models compared by Tunuguntla et al. (2016b) reduce to the same depth-averaged segregation model that is used in this paper. Moreover, this depth-averaged model can be coupled to single-phase shallow layer models to investigate the mobility feedback that segregation provides (Woodhouse et al. 2012). This simple model represents the leading-order behaviour of the flow, as demonstrated by the fact that Woodhouse et al. (2012) successfully showed it can reproduce the fingering instability of bidisperse granular flow fronts.
In this study, we adapt the model of Woodhouse et al. (2012) by including a non-unity shape factor for the velocity profile in the momentum balance, and using a Bulbous head formation in bidisperse shallow granular flow 267 different choice for the computation of the macroscopic friction coefficient. We then show that the bulbous-head profile is also a solution of this kind of depth-averaged continuum framework when using simple initial and boundary conditions. Furthermore, we demonstrate that the bulbous head also develops in three-dimensional discrete particle simulations, with reasonable agreement between the simple one-dimensional continuum model and the fully three-dimensional discrete particle simulations. This paper is structured as follows: our model adapted from Woodhouse et al. (2012) and our choices for the constitutive relations are discussed in § 2. The results of the numerical solutions of this model are described in § 3, where it is shown that the bulbous head emerges. Furthermore, it is shown in § 4 that our numerical solutions to the continuum model converge to a simple travelling wave solution, which makes computation of the long-time flow properties very efficient. Finally, in § 5 the results from the one-dimensional continuum model are compared to fully three-dimensional discrete particle simulations, with reasonable agreement between both approaches, despite the simplifications in the continuum model.

Continuum model
In this section, we describe the size bidisperse shallow granular model that is used in both the time-dependent numerical solutions of § 3 and the travelling wave solution of § 4. The model is based on Woodhouse et al. (2012), with some generalisations that are discussed later in this section.

Definitions
We consider a plane at an angle θ to the horizontal with the x-axis pointing downslope, the y-axis pointing across the slope and the z-axis being aligned with the upward pointing normal. In this coordinate system, the gravity vector, g, has the direction (sin θ , 0, −cos θ ) and magnitude g. Particles of two different sizes are released onto the inclined plane, namely small particles and large particles. These have the same density but different radii. Where appropriate, the different sizes are denoted by superscripts L and S for large and small particles, respectively. A schematic overview of this geometry is given in figure 2.
It is assumed that gradients and velocities in y-direction are neglected. We are interested in the height, h(x, t), the depth-averaged downslope velocity,ū(x, t), and the depth-averaged small particle concentration,φ(x, t), at all positions, x, at any time, t > 0. Here, the small particle concentration, φ, is defined as the volume fraction of the solid phase that consists of small particles, so the large particle concentration equals (1 − φ). The depth-averaged velocity and small particle concentration can be derived from their full two-dimensional equivalents u(x, z, t) and φ(x, z, t) and the flow height, h(x, t), asū ( 2.2) The length scale for the non-dimensionalisation is chosen to be the large particle diameter, instead of a typical flow height that is usually used with this kind of model. This leads to an easier comparison with discrete particle simulations later on. All 268 I. F. C. Denissen and others x g oe y z u FIGURE 2. Schematic drawing of a bidisperse chute flow. The granular material flows down a chute inclined at an angle θ to the horizontal. A coordinate system is taken with the x-axis aligned with the downslope direction and the z-axis normal to the chute's surface. The flow is assumed to be uniform in the cross-slope y-direction. Due to the gravity, g, pointing down, the avalanching material flows in the positive x-direction with a downslope velocity, u(x, z, t). We assume particle-size segregation completely separates the two particle sizes, where large particles are on top of small particles.
variables are non-dimensionalised by the large particle diameter, d * L , and gravitational constant, g, as follows: the height and length of the flow are non-dimensionalised by d * L , time by d * L /g and consequently the depth-averaged velocity by d * L g. We also non-dimensionalise the particle diameter for both types of particles with d * L , and thus d S = d * S /d * L and d L = 1, respectively. The bulk density of the flow is assumed constant, and hence it does not appear in the non-dimensionalised system of equations in § 2.3.

Assumptions
The full two-dimensional fields for the small particle concentration and velocity of the flow are needed for the derivation of the one-dimensional depth-averaged system of equations in § 2. 3 (Woodhouse et al. 2012). This reconstruction is not unique, and needs assumptions on both the segregation behaviour and the velocity profile.
Regarding the segregation behaviour, it is assumed that the flow is already segregated at the inflow boundary of the domain, where the large particles flow above the small particles. Furthermore, they are assumed to stay segregated, so there is no diffusive remixing (Gray & Chugunov 2006). While this is a rather strong assumption, it is fairly accurate and validated by experiments and discrete particle simulations of bidisperse chute flows, where it is shown that the segregation rate is 10-30 times higher than the diffusive-remixing rate (Wiederseiner et al. 2011;Thornton et al. 2012b). The combination of instantaneous segregation and the absence of diffusive remixing is equivalent to the assumption that the flow is fully segregated at all times and positions. This can be summarised as where hφ can also be referred to as the height of the small particle layer. Note, that the assumption of a fully segregated flow leaves the depth-averaged segregation equation heavily dependent on the velocity profile, i.e. on how much faster the large particles in the top layers are moving compared to the small particles in the bottom layers.
The velocity profile for bidisperse flows over a rough base is complicated and cannot be described by a simple linear or Bagnold velocity profile (Rognon et al. 2008;Tripathi & Khakhar 2011). However, for very shallow monodisperse flows over sufficiently rough bases, a linear profile is a good approximation (Silbert, Landry & Grest 2003;GDR-MiDi 2004;Weinhart et al. 2012), while a Bagnold velocity profile is more appropriate for thicker monodisperse flows (>20 particle diameters) (Silbert et al. 2001;Tunuguntla et al. 2014) and monodisperse flows over smooth bases (Brodu, Richard & Delannay 2013;Faug et al. 2015). For the shallow bidisperse flows studied in this paper, we use a simplified linear velocity profile with basal slip, of the form Here, α s is the amount of basal slip; for α s = 1, equation (2.4) describes a plug flow, while for α s = 0 it describes a linear velocity profile with no slip at the base. The choice of the constant value of α s is discussed in § 2.4. Note that α s is often denoted by α, see e.g. Gray & Thornton (2005), Woodhouse et al. (2012) and Baker et al. (2016). However, this causes a conflict in nomenclature with the shape factor, which we here denote by α, see (2.5), as is used by many authors e.g. Forterre & Pouliquen (2003), Weinhart et al. (2012) and Saingier, Deboeuf & Lagrée (2016). It must be noted that the linear velocity profile (2.4) is deeply embedded in the depth-averaged segregation equation (2.9), and that this equation will change drastically for different velocity profiles. For a detailed derivation of the depth-averaged segregation equation with other velocity profiles, we refer the interested reader to Baker et al. (2016).
Using the velocity profile (2.4), we can compute the shape factor, α, as This is where our model differs from that of Woodhouse et al. (2012), who set α = 1 independent of α s . We choose to retain the dependency of α on α s to keep the model consistent. Moreover, choosing α = 1 is important to preserve the influence of inertia terms, especially for rapid flows (Saingier et al. 2016). The downside of using a non-unity shape factor is that the slope of the height at the front goes to zero, which causes numerical difficulties. Note, that the assumption α = 1 also shows the bulbous-head profile. However, we prefer to keep the model internally consistent and hence choose α as given by (2.5).
2.3. System of equations With the above assumptions on the segregation behaviour and the velocity profile, we can now formulate our depth-averaged, non-dimensionalised governing equations, based on the model of Woodhouse et al. (2012). The mass and momentum balances are respectively given by where the source term, S, consists of the down-slope component of the gravitational acceleration and the basal friction respectively: The gravitational acceleration drives the flow in the positive x-direction, while the basal friction opposes the avalanche motion. It can be described with a basal friction coefficient, µ, which is the ratio of the shear and normal stress at the base. Traditionally, µ is taken either as a constant, e.g. Savage & Hutter (1989), or as a function of the height and velocity of the flow, e.g. Pouliquen (1999b). Here, we generalise the computation of µ to also include a dependence on the depth-averaged concentration of small particles, as discussed later in this section. Furthermore, the depth-averaged conservation of small particle mass can be expressed as (Gray & Kokelaar 2010a) where it must be noted that hφu = (hφū − (1 − α s )hφū(1 −φ)) when using the linear velocity profile with basal slip (2.4). Note, that the macroscopic friction coefficient, µ, is the only feedback mechanism affecting the segregation behaviour of the bulk; the depth-averaged mass and momentum balances (2.6)-(2.7) do not involve gradients in the concentration of small particles,φ. To make the macroscopic friction coefficient, µ, dependent on the depth-averaged small particle concentration,φ, we follow the approach of Pouliquen & Vallance (1999): the total basal friction is written as a concentration-weighted sum of the macroscopic friction coefficients µ S and µ L for monodispersed flows of small and large particles respectively: This provides a simple way of increasing the frictional resistance to motion as the proportion of larger particles increases at the flow front (with µ L > µ S ). Other friction laws are possible, for example the basal friction could be dependent on the local concentration of large and small particles at the base of the flow, or there may be a complex nonlinear dependence, as shown in the discrete particle simulations of Rognon et al. (2008). For the basal friction coefficients of each pure individual constituent, we use the non-dimensionalised version of the friction law for rough beds of Weinhart et al. (2012), which is closely related to Forterre & Pouliquen (2003), and has the form where the Froude number F = |ū|/ √ h cos θ can be expressed in terms of the non-dimensional velocity,ū, and height, h, of the flow, and the z-component of gravity, cos θ . The definition of the Froude number is different in Weinhart et al. (2012) compared to Forterre & Pouliquen (2003), which results in a slightly different form of the friction law. Here we use the definition of Weinhart et al. (2012).
Concerning the non-dimensional friction parameters for each constituent ν, δ ν 1 is the minimum chute angle required for flow, δ ν 2 is the maximum chute angle for which steady flows are possible, A ν is a characteristic dimensionless scale and β ν is an empirical constant. The friction law (2.11) is the second point where our model differs from Woodhouse et al. (2012), who used an exponential version of this friction law (Pouliquen 1999b). The advantages of the reciprocal friction law (2.11) are that it both gives a better fit, and is computationally faster compared to the exponential friction law (Weinhart et al. 2012). Moreover, it has been shown by Jop, Forterre & Pouliquen (2005) that the friction law (2.11) is closely related to the popular µ(I)-rheology (GDR-MiDi 2004) evaluated at the base.
One of the major advantages of the model described in this section is that it consists entirely of non-dimensional quantities; one can use this model across many different scales, from the laboratory scale to full particulate geophysical flows. Recently, Kokelaar et al. (2018) presented field evidence that these models can describe large-scale features on the moon, in addition to the laboratory-scale phenomena they have already been shown to capture.

Closure relations
In this section, closure relations for the macroscopic friction coefficient, µ, and the shape factor, α, are presented for incorporating the material-dependent properties into the model (2.6)-(2.11). To determine the friction parameters in (2.11), one can use small-scale periodic discrete particle simulations, as shown by Thornton et al. (2012a) and Weinhart et al. (2012).
For glass beads with a volume ratio V L /V S = 2 on a rough bottom of small beads, the friction parameters have been determined by Voortwis (2013), using the open-source software package MercuryDPM (Thornton et al. 2013a,b;Weinhart et al. 2017).
Similarly, the shape factor, α, can be measured from small-scale discrete particle simulations (Weinhart et al. 2012). For the shallow bidisperse flows in this study, we observed a height dependence for this shape factor, similar to the relation for monodisperse flows seen by Weinhart et al. (2012): with increasing height, α decreases asymptotically to a fixed value, in our case approximately 1.2. This is lower than the value of 1.25, observed by Weinhart et al. (2012), for monodisperse flows with a Bagnold velocity profile. For simplicity, we take α = α(h in , θ ) constant throughout the domain, depending only on the inflow height and the chute angle. For example, for inflow height h in = 15 and angle θ = 24 • , a constant value α ≈ 1.23 is used throughout the whole domain. Combined with (2.5), this results in α s ≈ 0.17 as the constant for the basal slip in the velocity profile (2.4). All relevant parameters for the continuum model used to obtain the results presented in this paper are summarised in table 1.

Time-dependent numerical solution
In this section, we solve the model outlined in § 2 numerically in order to show that it leads naturally to the bulbous head, where the flow is thicker downstream near the front than further upstream. Since the system is hyperbolic, it can be simulated with a discontinuous Galerkin finite element method (DGFEM) (Reed & Hill 1973). This type of method combines the numerical fluxes of finite volume methods with the high-order polynomial approximations through basis functions of finite element methods. DGFEMs are widely used for computational fluid dynamics, see for example 272 I. F. C. Denissen and others

Friction parameters
Other parameters Small particles Large particles . Summary of parameters for the continuum model used in this paper. The friction parameters, δ 1 , δ 2 , A and β, depend on the granular materials used, while θ, h in andφ in depend on the geometry. The procedure we used for determining the friction parameters, for a given granular material, is based on the experiments of Pouliquen (1999b) and is fully detailed in Thornton et al. (2012a). The basal slip, α s , depends on both the material and the geometry. Shu (2013) and the references therein. One of the advantages of DGFEM over conforming FEM is that it can deal with discontinuities in the solution, without producing spurious oscillations. We follow the general DGFEM framework for a second-order scheme with the Harten, Lax and van Leer (HLL) flux (Harten, Lax & van Leer 1983), the total variation bounded (TVB) limiter of Cockburn & Shu (1989) and the wetting/drying treatment of Bunya et al. (2009). The details of the numerical method used can be found in appendix A. The method is implemented in the hpGEM software package (Pesch et al. 2007;Nurijanyan 2013), which provides a framework for DGFEM simulations. Due to the wetting and drying treatment and the small but finite gradients of the height in the front, the accuracy of the method is reduced to first order.
A domain x ∈ [0, x end ] is constructed with an appropriate choice of x end , such that this end boundary does not influence the flow profile. The initial conditions are chosen to represent an empty chute, so h(x, 0) = hū(x, 0) = hφ(x, 0) = 0. At the inflow boundary x = 0, Dirichlet boundary conditions are prescribed to fix the inflow values, h(0, t) = h in , hū(0, t) = h inūin and hφ(0, t) = h inφin . These inflow values are chosen such that the flow at x = 0 is non-accelerating, i.e. S(h in ,ū in ,φ in ) = 0. At the outflow boundary x = x end , outflow boundary conditions are prescribed (Bristeau & Coussin 2001). Since the bulk flow never reaches this boundary, the exact boundary conditions do not matter, as long as they keep the height, momentum and small particle height non-negative and nearly zero.
For the results presented in this paper, a granular flow over a chute at angle θ = 24 • has been simulated until t end = 2500. For this end time, an appropriate choice for the length of the domain is x end = 10 000. In order to be able to compare the DGFEM solutions with the discrete particle simulations of Voortwis (2013), the inflow height and small particle concentration are prescribed as h in = 15,φ in = 0.5. For a flow with uniform inflow height, it then follows thatū in ≈ 3.4, see § B.1. For DGFEM simulations, we have chosen a time step of t = 10 −4 and N = 250 elements. All parameter values used in the numerical solutions are summarised in table 1.
The height and concentration profiles that emerge from these DGFEM solutions can be found in figure 3. Since the initial and boundary conditions are exactly the same as in a dam-break structure, the profiles are initially very similar to this structure, see figure 3(a) (Hogg & Pritchard 2004). We see that the bulbous head has not formed yet, and the profile of the small particle concentration is smooth. Moreover, it is (a) . Height profiles at various times generated by DGFEM solutions. The thick black curve denotes the height, h, of the flow, the grey area is bounded by the height of the small particles, hφ. Both x and z are scaled by the large particle diameter, d * L . Initially, the bulbous head shape develops (a,b), until the head reaches its maximum height (c). It then advects downwards, with the faster large particle front growing longer at a constant rate. The dotted red and dashed blue lines illustrate the different speeds of the shock position and large particle front, respectively. Note, the axis has been significantly compressed in the x-direction, in order to fit the page. An animation of this DGFEM solution is included in the online supplementary material available at https://doi.org/10.1017/jfm.2019.63.
noteworthy that as a result of choosing α = 1, the height tends asymptotically to zero and a thin precursor layer arises. The smooth profile with gradient continuously decaying to zero is consistent with monodisperse dam breaks and chute flows, for which analytical solutions have been derived by Hogg & Pritchard (2004) and Saingier et al. (2016). Around t ≈ 250 the flow starts to develop a clear front of only large particles, separated by a shock in the small particle concentration, see figure 3(b). Furthermore, the friction causes the large particle front to have a lower speed than the inflowing material. This makes the height deviate from the dam-break structure, towards a structure in which the height at the front is higher than at the inflow. Once the head is fully developed, around t ≈ 1000 (figure 3c), the bulbous-head height stays constant. Near the inflow, the flow is uniform, with the prescribed height, velocity and small 274 I. F. C. Denissen and others particle concentration. The bulbous head at the front is now clearly higher than the inflow height, with a very pronounced large particle front.
When comparing the profiles at t = 1000, t = 1500, t = 2000 and t = 2500 (figure 3c-f ), we see that the shock in small particle concentration moves with a constant rate, as illustrated by the dotted red line. In the region just behind this shock, the height profile and small particle concentration profile stay constant with respect to the shock. Near the inflow, the height and small particle concentration stay constant at their prescribed values. So up to the shock position, the profiles of the height and small particle concentration essentially stay the same, they are simply advected in the x-direction and have a growing 'tail'. This indicates that there may be a travelling wave solution of the model we presented in § 2. In § 4 we derive such a travelling wave solution.
The large particle front grows longer with time, and moves with a constant speed. This speed is illustrated by the dashed blue line in figure 3(c-f ). The speed of the large particle front is higher than the shock speed, which is represented by the different slopes of the dotted red and dashed blue lines in figure 3(c-f ). The shape of the front stays the same, suggesting that this can be described by another travelling wave solution of this model. In § 4, it is shown that the large particle front indeed fits a travelling wave solution, namely the monodisperse front solution derived by Saingier et al. (2016), which is a generalisation for α = 1 of the solutions of Pouliquen (1999a) and Gray & Ancey (2009).
In the next section, the travelling wave solution for the region behind the shock will be derived. Furthermore, we show that the numerical solutions of this section converge to a combination of two travelling wave solutions with different speeds, one with a shock in the small particle concentration and one for the monodisperse front.

Travelling wave solution
In this section, a travelling wave solution is derived based on the observations from the numerical solution of § 3. That is, we seek a steady solution in a reference frame moving at constant speed. So if we can capture this wave in a reference frame, there exists a frame speed, u f , for which the travelling wave solution is a steady-state solution in this moving frame.
Based on our observations in § 3, we postulate a travelling wave solution to system (2.6)-(2.9) in which the profiles of h andū are continuous, but there is a steady moving shock inφ. This shock speed is the same as the frame speed, u f , as we want the shock to stand still in the moving frame. We can introduce a new coordinate system such that the frame moves with speed u f in the positive direction with ξ = 0 being the fixed position of the shock. Applying this coordinate transformation leads to the system Note, that the source term, S, given by (2.8) remains unchanged by this transformation since it does not directly depend on x or t.
As we are interested in a travelling wave solution, we assume that the solution of the system is in steady state with respect to the new coordinate system. Therefore, we set the τ -derivatives to zero, so the system becomes a set of coupled ordinary differential equations: Equations (4.5) and (4.7) can be integrated directly: where C 1 and C 2 are arbitrary integration constants. Based on the numerical results in § 3, we assume that the height, h, momentum, hū, and small particle height, hφ, go asymptotically to a fixed value for ξ → −∞, so h → h in ,ū →ū in andφ →φ in as ξ → −∞ for some constants h in > 0,ū in > 0 andφ in ∈ [0, 1]. Using (4.8) and (4.9), the integration constants, C 1 and C 2 , are then given by (4.11) Substituting (4.10) into (4.8) gives the relation for the height depending on the depthaveraged velocity Similarly, substitution of (4.11) into (4.9) and rearrangement of the terms shows that φ is a solution to the quadratic equation This equation can be solved with the quadratic formula, so ifū is known, then h(ū) follows from (4.12), andφ(ū) can be found from (4.14) Note that the other solution to the quadratic equation (4.13) gives a negative value for φ, which is unphysical. Using these results to eliminate h andφ from (4.6) leads to a first-order nonlinear ordinary differential equation forū:

Boundary conditions
As in § 3, we prescribe the inflow height, h in , and the inflow concentration of small particles,φ in . We also want the height to go to the inflow height asymptotically for ξ → −∞, so (dh/dξ ) in = 0 is prescribed. Using (4.5) and the chain rule, it follows that dū dξ Equations (4.15) and (4.16) now imply thatū in must satisfy S(h in ,ū in ,φ in ) = 0. For the current friction law (2.10)-(2.11), this is given by: where f (φ in ) is a function ofφ in , the friction parameters and the geometry. For the exact form, see § B.1. Furthermore, it is beneficial to know the values h out andū out at the outflow position of the domain, ξ → ∞. We look for a flow that is steady for ξ > 0, with no small particles in this region, soφ = 0 if ξ > 0. To determine the values of h out andū out , we use the global mass balance, Furthermore, we prescribe that the flow is non-accelerating and uniform with only large particles in the front: Note, that conditions (4.17) and (4.19) only hold for steady flows; in particular, in the initial phase of the flows described in § 3, these conditions do not hold. Conditions (4.18) and (4.19) imply that the height and depth-averaged velocity have constant values h = h out andū =ū out for ξ > 0. They can be computed with Matlab's solve routine from (4.20) The value of h out then follows directly from (4.18).

Shock properties
To compute the frame speed, u f , consider a plane perpendicular to the chute at a certain ξ = ξ 0 < 0. Since we are looking at a steady-state solution in the moving frame and there are no small particles at the front, we know that the net volume flux of small particles through this plane is zero. Because we assume that the flow is fully segregated, this can be expressed as (4.21) Evaluating this integral at ξ → −∞ with the linear velocity profile (2.4), it can be shown that the shock speed is given by Bulbous head formation in bidisperse shallow granular flow

277
To compute the small particle concentration at the left side of the shock, we realise that for a general hyperbolic system of partial differential equations of the form ∂u/∂t + ∂f (u)/∂x = 0, the Rankine-Hugoniot jump condition states that the shock speed s satisfies s(u 2002). Applying this to the relation for the small particle concentration (2.9) and shock speed equal to the frame speed, u f , we find where (φ) − = lim ξ ↑0φ and (φ) + = lim ξ ↓0φ . Since h andū are continuous at the shock position, it holds that (h) − = (h) + = h out and (ū) − = (ū) + =ū out . Working out the parentheses and substituting these relations for the height-and depth-averaged velocity in (4.23) gives (4.24) Assuming the small particle concentration at the downstream position of the shock disappears, (φ) + = 0, it follows from (4.24) that is the value of the small particle concentration at the left side of the shock.

ODE solution
The last missing step from the travelling wave solution is the exact shape of the height, velocity and concentration profiles. In order to compute these profiles, equation (4.15) is solved with MATLAB's ode45 routine using the parameters of table 1. The boundary conditions at the shock position, ξ = 0, are h = h out ,ū =ū out and (φ) − =φ out and (φ) + = 0. Equation (4.15) is integrated in both positive and negative ξ -direction starting from ξ = 0 with the different initial conditions forφ.
The resulting profiles are displayed in figure 4. For ξ > 0, the flow is nonaccelerating, since the height, velocity and small particle concentration stay constant, in accordance with constraint (4.19). There are no small particles present, sinceφ = 0 for all ξ > 0. The flow for ξ < 0 is more interesting, where it is important to note that h → h in ,ū →ū in andφ →φ in asymptotically as ξ → −∞. All curves are smooth, except for the shock inφ at ξ = 0.
Comparing the inflow and outflow heights and velocities, we see that the flow in the front is deeper and slower than the flow in the back. This is caused by the higher friction in the front, since there are only large, more frictional, particles here. The higher friction coefficient causes more momentum to be dissipated, but the mass must be conserved; therefore, the flow of pure large particles must be deeper to conserve mass. Of course, the magnitude of these differences depend on the difference in friction of both constituents. The larger the difference in friction between the bidisperse and large particle phases, the deeper and slower the flow in the front. This influence was verified by changing various combinations of friction parameters in (4.15), e.g. by making A, δ L 1 and δ L 2 larger or smaller. The exact relation between the macroscopic friction coefficients and flow height for a bidisperse granular flow within the bulbous head needs further investigation. However, a detailed investigation of this The shock speed of u f ≈ 1.9 is significantly lower than the outflow speed, u out ≈ 3.1, which might look surprising at first. However, we assume that the velocity of the flow at the surface is higher than at the base, varying linearly as described by (2.4). The small particles are at the base, and therefore moving much slower than the large particles at the free surface. Thus the shock in the small particle concentration must be located at the base of the flow, since that is where all small particles reside. It follows that the shock speed must be the same as the average velocity of the flow from the base to the small particle height, which can be expressed as This is exactly the shock speed that we see, which is significantly lower than the outflow speed. This difference in u f andū out leads to a net transport of large particles to the front, causing the bulbous head to grow in volume. Note, that the net transport of large particles towards the front also directly follows from the travelling wave solution, see § B.2 for the exact expression. Next, this travelling wave solution is compared with the time-dependent DGFEM solution of § 3.

Comparison with time-dependent solution
To verify that this travelling wave solution describes the time-dependent solution of § 3, both solutions are compared by shifting the DGFEM solution such that the shock in small particle concentration is located at ξ = 0 for all times, see figure 5. Moreover, the monodisperse front is compared to the solution of Saingier et al. (2016), which satisfies the ordinary differential equation wherex is the distance from the point where the height vanishes. This is a generalisation for α = 1 of the monodisperse front solutions derived by Pouliquen (1999a) and Gray & Ancey (2009). By shiftingx, we can fix the mass such that the front solution agrees with the DGFEM solutions. Figure 5 shows the comparison between the DGFEM solution at t = 2500 and both ODE solutions for the height and height of small particles. We see similar agreement for the depth-averaged velocity of the flow. For ξ < 0, the time-dependent solution matches the travelling wave solution of this section. For ξ > 0, the time-dependent solution matches the solution of ODE (4.28). It should be noted that both ODE solutions travel with different velocities: the large particle front travels with a speed u out , which is larger than the shock speed, u f , leading to an increasing volume in the monodisperse head. Furthermore, the volume flow rate of large particles is identical for the travelling wave solution and the DGFEM solution, see § B.2. We therefore conclude that the DGFEM solution can be seen as a combination of two travelling wave ODE solutions, travelling at different velocities. In the next section we validate the model by comparing the DGFEM solution to discrete particle simulations.

Comparison with discrete particle simulations
To analyse the quality of the model described in § 2 for this kind of flow, the DGFEM solutions are compared with two preliminary sets of discrete particle simulations of a chute flow. The discrete particle simulations are performed using MercuryDPM (Thornton et al. 2013a,b;Weinhart et al. 2017), where the simulation details can be found in appendix C. The chute has the same inclination as in the DGFEM solutions and the particles are the same as used for the calibration of the friction parameters. To model the inflow boundary conditions, a so-called maser boundary has been developed: for x ∈ [−20, 0], a small periodic chute is simulated.

I. F. C. Denissen and others
Once the flow is fully developed, the periodic boundary is removed and all particles flowing through the boundary at x = 0 are both entering the non-periodic part of the chute (x > 0) and duplicated at the beginning of the domain, x ≈ −20. This way, the flow has a fully developed, steady velocity and segregation profile at the inflow. See § C.5 for full details of the maser inflow boundary. Figure 6 shows a time-evolving height profile from our initial discrete particle simulations. Initially, at t = 0, we only have the maser particles between x = −20 and x = 0, otherwise the chute is empty. By t = 100, the flow has already developed a front of purely large particles, that is higher than inflow height; a bulbous head is forming. As time progresses (t = 200 until t = 500), the bulbous head gets longer and higher, eventually saturating in height but always growing in length.

Height profiles
In figure 6, the discrete particle simulations are overlaid with the DGFEM solution at the same time. Note, that there is no fitting involved, since the friction parameters and shape factor were calibrated, i.e. measured with small-scale periodic discrete particle simulations using the same particle properties. The procedure we used for calibrating the friction parameters, for a given granular material, is based on the experiments of Pouliquen (1999b) and is fully detailed in Thornton et al. (2012a). Since the discrete particle simulations and continuum model of § 2 both use the large particle diameter and gravity for scaling, we can directly compare the results of both methods. Note that the x-axes start at x = −20 to show the maser particles.
The similarity between the predicted shape of the flow front and the discrete particle data is remarkably good, especially for larger t. Initially there are some discrepancies: the head is already forming at t = 100 in the discrete particle simulations and not in the continuum model. We hypothesise that this is because the DFGEM prescribed inflow conditions do not exactly match the outflow conditions of the maser, in particular, the effective friction between the flow and surface. The travelling distance of the flow is very well predicted for all t, even though the only condition that is prescribed in both simulation methods is that the flow is steady and non-accelerating at the inflow.
Still, the depth-averaged continuum theory does not fully capture all of the details of the height profile. There are two main differences between the shape of the head in the continuum model and discrete particle simulations. The first is at the very front of the flow, which in the discrete particle simulations shows a gaseous behaviour with a few saltating particles being ejected from the flow. This gaseous behaviour cannot be predicted by the continuum model, since a constant bulk density is assumed. However, the main discrepancy for the height profile between the continuum model and the discrete particle simulations is at the back of the bulbous head; e.g. around x = 1000 and t = 500. The discrete particle simulations display smoother behaviour than the continuum model, where the mass of the particles is also slightly more in the back compared to the numerical solution of the continuum model. We suspect the disagreement in the continuum and discrete particle simulations stems from the fact that the segregation structure in the front is approximated by a shock rather than the full structure of the breaking size-segregation wave, that is expected at the back of the bulbous head (Thornton & Gray 2008;Gajjar et al. 2016;van der Vaart et al. 2018). This difference modified the effective basal friction; and, hence, the bulk dynamics and shape of the front. To remove the discontinuity in the first derivative of the height profile, at the top of the bulbous head, one would need a more complex model that FIGURE 6. Snapshots of discrete particle simulation of a bidisperse chute flow over a rough bottom at various times. The inflow height is 15 particle diameters, resulting in supercritical inflow and the mixture is 50/50 large and small particles by volume. Red are the large particles, green are the small particles. The black lines denote the DGFEM solution for the height. Both x and z are scaled by the large particle diameter d * L , in both the discrete particle simulations and DGFEM solutions. Note, that the x-axis is squeezed by a factor 100 compared to the z-axis, so the individual particles appear as very thin vertical stripes in this plot.
retains more details of the vertical segregation structure. A more detailed discussion of the segregating structure can be found in § 5.2. Figure 7 shows a second comparison with particle simulations for the case h in = 7.4, α s = 0.08. These parameters result in subcritical inflow conditions; however the structure is still very similar to a bulbous head, forming and growing longer over time. In this case, the height is slightly under-predicted, in contrast to the super-critical FIGURE 7. Snapshots of discrete particle simulation of a bidisperse chute flow over a rough bottom at various times. The inflow height is 7.4 particle diameters, resulting in subcritical inflow and the mixture is 50/50 large and small particles by volume. Red are the large particles, green are the small particles. The black lines denote the DGFEM solution for the height. Both x and z are scaled by the large particle diameter d * L , in both the discrete particle simulations and DGFEM solutions. Note, that the x-axis is squeezed by a factor 100 compared to the z-axis, so the individual particles appear to be very thin in this plot. . Note the breaking size-segregation wave in the discrete particle simulations between x ≈ 500 and x ≈ 1500.
case, see figure 6. The height profile seems to be less steep in the continuum model compared to the discrete particle simulations. Furthermore, the continuum model overpredicts the flow velocity, which can partially be explained by the fact that the flow is subcritical. In the particle simulations the maser boundary is constructed such that it can interact with the main flow particles. The flow front is more resistive, due to an evolving segregation profile, which feeds back on the maser boundary lowering the inflow velocity produced by the maser. Therefore, we see that the inflow velocity is decreasing with time, in the discrete particle simulations, and lower than the steadyflow relation, equation (4.17), used by the continuum model. This is similar to the complex feedback between segregation profile, friction and hence velocity, seen in the breaking size-segregation structures studied by van der Vaart et al. (2018) and is not captured by the current continuum models. Further study is needed to investigate the limits of the current model with respect to the inflow height, Froude number and chute angle. Moreover, the inflow conditions for subcritical flow configurations should be examined.

Segregation profiles
Looking at the segregation profiles from the discrete particle simulations in more detail, we see that the shock in the small particle concentration is unphysical and that the region with only large particles is much smaller than predicted by the continuum model. In general, the segregation profiles in the discrete particle simulations are not as sharp due to diffusive remixing, and possibly some deposition of large particles at the front with re-circulation. In the front of the flow, we see a thin band consisting of only large particles. Behind that, there is an expanding structure around the shock position, see figure 8. The structure observed here in the discrete particle simulations 284 I. F. C. Denissen and others has been predicted from a non-depth-averaged continuum model (Thornton & Gray 2008;Gray & Ancey 2009;Gajjar et al. 2016) and in both experiments and discrete particle simulations of a moving belt (van der Vaart et al. 2018). This structure was named the breaking size-segregation wave and its study is ongoing. Nevertheless, these results indicate that the assumption of full segregation needs to be relaxed in order to obtain a better continuum model for the segregation behaviour in this region.

Conclusion and discussion
Granular flows containing a particle-size distribution, in both laboratory experiments and the natural environment, often show a bulbous head structure at their flow front. Using three-dimensional discrete particle simulations, and a simple one-dimensional depth-averaged continuum model, we have shown that this bulbous head structure is predicted at both the discrete and the continuum level. Furthermore, our long-time numerical solutions of the continuum model converge to a novel combination of two travelling wave solutions. This allows for an efficient computation of key features of the flow, such as the maximum flow depth and the mass flux in steady state.
The simple one-dimensional depth-averaged model, that we present for this complex problem, is calibrated using only small-scale, steady-state, discrete particle simulations, i.e. there are no fitting parameters. We performed a preliminary comparison (one super-and one subcritical) between the computationally cheap depth-averaged model (requiring minutes of computation on a single core) with two large expensive (requiring several months of computation on a single core) fully three-dimensional particle simulations. For the super-critical flow case the depth-averaged model was able to capture key flow properties e.g. maximum flow depth, and propagation speed; however, it was not able to accurately capture the details of the segregating profile. For the subcritical case discrepancies arose and increased with time; the problem, lying with the inflow boundary condition for the depth-averaged continuum model. A more detailed comparison between the continuum model and the full-scale particle simulations should be undertaken and the full range of validity of the simple continuum model determined. However, given the efficiency of the model it could be an elegant tool for investigating the leading-order behaviour of complex segregating geophysical flows in the future.
As stated above, despite capturing various bulk flow properties the current continuum model does not capture all the details of the segregation behaviour, since it assumes full segregation at every position. There are many ways the model could be improved with respect to the segregation profile. Firstly, one could insert a diffusion layer between the pure large and pure small phases as done by Edwards & Vriend (2016). Alternatively, one could fit the 'correct' two-dimensional breaking size-segregation wave profile (Thornton & Gray 2008;Gray & Ancey 2009;Gajjar et al. 2016;van der Vaart et al. 2018) at the one-dimensional shock position. This structure is known to appear in similar systems, both in discrete particle simulations and solutions of the full two-dimensional segregation equation. Furthermore, the breaking size-segregation wave becomes a simple shock when the segregation equation is depth averaged (Gray & Kokelaar 2010a,b). Finally, we could even couple the two-dimensional segregation equation to the one-dimensional depth-averaged bulk flow model. For example, one could construct a two-dimensional domain which is bounded by the flow base, inflow boundary and the one-dimensional height profile. On this domain the two-dimensional segregation model can then be solved numerically, using a velocity profile reconstructed from the computed depth-averaged velocity via (2.4). This formulation would allow the diffusion terms to be retained within the segregation model and hence diffusive-remixing effects could also be captured.
From a mathematical point of view, it would be interesting to check under which conditions our depth-averaged model is well posed. The closely related model of Woodhouse et al. (2012) is ill posed when the characteristic from the segregation equation coincides with one of the characteristics of the bulk flow model, i.e. when the model fails to be strictly hyperbolic. However, it becomes unconditionally well posed if a viscosity term is added to the momentum balance (Baker et al. 2016). Thus, another possible step is to see what happens to the travelling wave solution if such a viscous term is added to our model. The addition of this term has recently been investigated for monodisperse flows by Gray & Edwards (2014) and by Baker et al. (2016) for the model of Woodhouse et al. (2012).
Special attention should also be paid to the use of the shape factor, α = 1, as this makes the complexity of the computational algorithm for the continuum method a lot higher; requiring the use of a more sophisticated wetting and drying treatment (see § A.4). The inclusion of the non-unity shape factor creates a thin precursor layer which causes a strong restriction on the time step of the numerical DGFEM solutions. The numerical solution of the continuum model is closer to the particle simulations with the inclusion of the 'correct' shape factor; but, the general structure of the height and velocity of the bulbous head are all captured by the approximate case α = 1. A more detailed comparison between the DGFEM and particle simulations forms the central theme of future work where the cost versus error of including a non-unity shape factor will be further investigated. For obtaining accurate closure relations and performing a detailed comparison between discrete and continuum models, it is helpful to coarse grain the discrete particle simulations. Coarse graining allows us to extract three-dimensional continuum fields from the discrete particle data, for e.g. the momentum and stresses (Babic 1997;Goldhirsch 2010;Weinhart et al. 2013), both for the flow as a whole and for each constituent separately (Tunuguntla, Thornton & Weinhart 2016a). This, in turn, can show us which assumptions in the continuum model are most restrictive and help us to improve full three-dimensional as well as shallow granular segregation models. Currently, coarse-grained fields have not been constructed for the system studied here.
To make the depth-averaged shallow granular flow framework more general, it can be extended from bidisperse to tri-and polydisperse mixtures, for which segregation frameworks have been developed by  and Marks et al. (2012), respectively. Both models have the same structure as the current model for bidisperse flows, and it would therefore be relatively straightforward to perform the same analysis as is done in this paper. It has also been shown by Schlick et al. (2016) that such models can predict segregating patterns in the formation of heaps, when first calibrated with discrete particle simulations. Due to of the weak coupling between the bulk flow model and the segregation model, one could even think of coupling a depth-averaged segregation model, such as the one of Gray & Kokelaar (2010a), to a model for water-saturated granular chute flows, for example the model described in Kowalski & McElwaine (2013).
Currently, the calibration of the constitutive relations for our chute flow model has to be repeated for each different combination of materials and particle sizes. To make the continuum model more applicable, the friction laws for bidisperse flows should be validated and extended to varying size ratios and mixture compositions. This can be done by determining a more general friction law, or by designing efficient discrete particle simulations that can determine the friction coefficient on the fly. The 286 I. F. C. Denissen and others same holds for the a priori unknown velocity profile, which generally depends on many parameters, such as base geometry and height. Baker et al. (2016) have already developed a framework for the depth-averaged segregation equation with a general velocity profile, and extending this work is an interesting avenue of future research.
Bulbous head formation in bidisperse shallow granular flow 287 A.2. Discretisation We now put the system into a discrete weak form. Multiply each equation in (A 1) with an arbitrary test function v ∈ V h and integrate over each element, which results in Integrating the second term by parts yields: There are discontinuities in U h allowed on the element faces ∂K, therefore the flux F(U h ) is replaced by a numerical fluxF(U L h , U R h ), where U L h and U R h are the values at ∂K of the left and right elements respectively, giving the discretisation Here we use the numerical flux of Harten et al. (1983), which is positivity preserving and works well for the similarly structured one-dimensional shallow water equations (Kesserwani et al. 2008;Toro 2013). For completeness, this is repeated here. Define a L = αū − √ (α 2 − α)ū 2 + h cos θ = λ 1 and a R = αū + √ (α 2 − α)ū 2 + h cos θ = λ 2 as the characteristic wave speeds of the bulk flow and define f L = F(U L h ) and f R = f (U R h ), respectively, the HLL flux is then given bŷ A.3. Slope limiting Since the system is hyperbolic, shocks in the solution can form over time. This can lead to severe oscillations near the shock, which can be prevented by using a slope limiter. Here, the TVB limiter of Cockburn & Shu (1989) is used.
A.4. Wetting and drying treatment Physically, the height of the particle flow can never become negative, so this should also be enforced in the numerical solutions of the continuum model. In order to do so, we use the wetting and drying treatment of Bunya et al. (2009), which conserves mass and, if the average height in an element is large enough, momentum. It is a shallow layer approach, in the sense that if the height at a node is below a certain threshold, the slope of the height and momentum are changed such that the height is above that threshold everywhere. If the average height in an element is below that threshold, the height in that element is set to its average and the momentum is set to 0. Note that we expect a very thin precursor layer in front of the flow, since α = 1 (Saingier et al. 2016). This layer can be much smaller than the numerical precision, so it is important to correct the height every time step. As the gradient of the height also goes to zero, we must not be too aggressive with our limiting, as this causes numerical artefacts where the gradient of the height gets the wrong sign for some elements. Here, the threshold is chosen to be 10 −10 , which is well above numerical precision, but low enough that numerical artefacts are minimal. 288 I. F. C. Denissen and others A.5. Time integration With some basic manipulation, the system of (A 6) can be written as where U are the coefficients and L is some nonlinear operator. We use the forward Euler method (Press 2007) to compute all coefficients at the next time step: where Υ is the operator that denotes the combination of slope and non-negativity limiting and t is the time step. From experience we know that the most severe restriction on the time step comes from the wetting and drying treatment, so a higherorder time integration method is not needed for a better stability and accuracy of the method.
Appendix B. Details of the travelling wave solution

B.1. Derivation ofū in
In this appendix, the exact formula forū in such that S(h in ,ū in ,φ in ) = 0, i.e. µ(h in ,ū in ,φ in ) = tan θ, will be derived. It is assumed that h in ,φ in and all friction and system parameters are known beforehand, so thatū in is the only unknown. The friction law is given by Introducing the following notation: one can rewrite (B 2) as Rearrangement of the terms then shows that F can be expressed as B.2. Conservation of large particle volume It must be noted that in the travelling wave solution of § 4, the values of h andū for ξ > 0 are constant, so the volume of the large particles in the front is infinite. However, there is a net flow of large particles towards the front, as can be seen in the time-dependent solutions of § 3. This transport rate for large particles Q L can be computed by taking the difference of the volumetric flow rate at the inflow ξ = −∞ and outflow ξ = ∞: using the assumptions of a fully segregated flow and applying velocity profile (2.4). For the parameters in table 1, this gives the value Q L = 22.0. Since Q L > 0 in this situation, we expect the volume of the head to grow over time based on this analysis.
The front of the time-dependent DGFEM solutions are integrated numerically and evaluated at t = 2000, 3000 and 4000, see table 2. We see that the volume of the large particle front increases at a rate of Q L · t, which shows that the global mass balance and hence the travelling wave solution accurately predicts the growth of the head.

Appendix C. Details of discrete particle simulations
In this appendix, details are given regarding the discrete particle simulations used to determine the friction parameters in table 1 and to produce figure 6. All simulations are performed with MercuryDPM (Thornton et al. 2013a,b;Weinhart et al. 2017), an open-source package designed for performing discrete particle simulations. It has been used for bidisperse chute-flows before, see e.g. (Thornton et al. 2012b). 290 I. F. C. Denissen and others

C.1. Definitions
The system consists of two types of soft, spherical particles that have the same density, but different diameters and therefore different contact properties. On a microscopic scale, both types of particles have the same restitution coefficient and dimensionless contact time, but the large particles have a larger microscopic friction than the small particles, see § C.2. Each particle i has diameter d i and density ρ p . The system variables are, for each particle i, its position r i , velocity v i and angular velocity ω i .
Each contact between two particles is treated as acting at a single point. The distance between two particles i and j is given by r ij = |r i − r j |, and the overlap between these particles is δ n ij = max(0, (d i + d j )/2 − r ij ). Note, that it is assumed that this overlap is small, otherwise contacts cannot be treated as being point-like. If two particles are in contact, one can define the unit normal vectorn ij = (r i − r j )/r ij and the vector from the centre of the particle to the contact point, the branch vector, b ij = −(d i − δ n ij )n ij /2. The relative velocity v ij = v i − v j can be decomposed in a normal and tangential component as follows: The tangential displacement δ t ij can be evaluated as (Weinhart et al. 2012 (C 4) which is integrated with the first-order forward Euler method (Press 2007), starting from the initial time of contact t 0 .
C.2. Contact law There are many different contact laws to evaluate the forces between colliding particles. Here, the normal and tangential forces are modeled with linear elastic and linear dissipative contributions (Cundall & Strack 1979;Luding 2008). The contact law does not strongly affect the macroscopic friction coefficient (Thornton et al. 2012a), which justifies choosing the simple linear contact model. The total contact force on particle i due to particle j is then given by f ij = f n ij + f t ij , where the normal and tangential forces are given by f n ij = k n δ n ijn ij − γ n v n ij , (C 5) with normal and tangential spring constants k n and k t , and damping constants γ n and γ t respectively. The contact time t c for a central collision can be related to the contact properties as where the reduced mass is given by m ij = m i m j /(m i + m j ). The restitution coefficient is given by = exp(−t c γ n /2m ij ). (C 8) It should be noted, that our definition of γ n differs with a factor two from the definition of Silbert et al. (2001), but is in line with the definitions of Luding (2008) and Weinhart et al. (2012). The contact-stiffness and damping are chosen such that the restitution coefficient, which is a material property, is the same for all types of collisions, LL = SS = LS , (C 9) where the superscripts LL, SS and LS denote large-large, small-small and large-small contacts respectively. Furthermore, the ratio between collision time and gravitational time is kept constant, which means that . (C 10) For each type of contact, the tangential stiffness is related to the normal stiffness as k t = 2/7 k n and the the tangential damping coefficient equals the normal damping coefficient.
Furthermore, Coulomb friction on the contact level has been assumed, which means that the magnitude of δ t ij is capped if necessary such that | f t ij |/| f n ij | µ p . For this study, the value of µ p depends on the size of the particles, where it is larger for large particles than for small particles. The relevant parameter values can be found in table 3. For a more detailed description of this contact law we refer to (Weinhart et al. 2012).
The total force F i and torque q i acting on particle i are now given by (C 12) where g is the gravitational vector, see figure 2.
C.3. Dynamics The dynamics of a particle i with mass m i and moment of inertia I i can now be described by These equations of motion are integrated with a second-order velocity Verlet (leapfrog) scheme (Allen & Tildesley 1989) to compute the position, velocity and angular velocity of each particle at each time step. All necessary parameters for the discrete particle simulations are given in table 3, non-dimensionalised such that for a large particle d L = 1, m L = 1 and g = 1. From these expressions, all properties for collisions between small particles and mixed collisions can be derived, see the thesis of Voortwis (2013) for details.

Parameter
Value Description g 1 Gravitational constant d L 1 Large particle diameter d S 3 √ 1/2 ≈ 0.79 Small particle diameter ρ p 6/π Particle density (m L = 1, m S = 0.5) k n,L 2 × 10 5 Normal spring constant large particles k t,L 2 7 k n,L Tangential spring constant large particles γ n,L 25 Normal damping coefficient large particles γ t,L γ n,L Tangential damping coefficient large particles µ p,LL 0.5 Microscopic particle friction large-large contacts µ p,SS 0.25 Microscopic particle friction small-small contacts µ p,LS 0.375 Microscopic particle friction large-small contacts t 10 −4 (≈ t LL c /50) Time step for time integration schemes TABLE 3. Nondimensionalised parameter values used for the particle simulations.
C.4. Chute geometry For all simulations, the base of the flow consists of a few layers of fixed large particles, to imitate a rough surface. To determine the friction parameters for the continuum model, a short periodic chute of size (x, y) ∈ [0, 20] × [0, 10] is constructed such that it tilts in x-direction. It is filled with 1000 large particles and 2000 small particles, so that the volume fraction of each is the same. For a description of the algorithm to get the closure relations from a short periodic chute flow, we refer the interested reader to Weinhart et al. (2012). For the long chute simulations, the base of the short periodic chute is copied in x-direction to a total length of 2000, while keeping the width in y-direction constant. The inflow is regulated through a maser boundary, see § C.5 for details.

C.5. Maser inflow boundary condition
For the long chute simulations, the inflow must be non-accelerating and segregated in order to compare them to the continuum model. One of the options is to do this with a hopper with a large reservoir of particles. However, this means that at all times, the positions and velocities of all the particles that will ever appear in the simulation have to be computed. Furthermore, it can neither be assumed that a flow from a hopper is already segregated at the inflow of the chute, nor that the mass fraction is constant over time, which makes comparison with continuum models more difficult.
As an alternative, we designed a maser inflow boundary, which generates a fully developed inflow for a chute. Initially, it only consists of a periodic box for (x, y) ∈ [−20, 0] × [0, 10]. After the periodic flow reaches a steady state, the maser boundary starts generating particles: for each particle that passes the line x = 0, the particle is both moved back to x = −20 as in a normal periodic boundary, and copied as inflow particle to its original position at the right, becoming an inflowing particle for the chute, see figure 9. The particles around x = −20 do not exert any forces on the particles around x = 0; otherwise, the particles around x = 0 would experience similar forces twice. However, the particles at x = 0 do exert forces at the particles around x = −20 by the means of ghost particles. This way, a fully developed inflow with nearly constant height is realised. It should be noted, that the velocity inside the maser boundary is influenced by the velocity in the outflow domain, since those particles do