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

Abstract

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.

1 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 Iverson et al. (2010) showed that segregation effects significantly increase the run-out distance of granular flows over inclined planes (chute flows). On the laboratory scale, experiments also show that granular flows of mixtures behave differently than flows of uniform materials, exhibiting interesting fingering effects (Pouliquen & Vallance 1999; Vallance & Savage 2000; Baker, Johnson & Gray 2016).

Size segregation can have a very strong effect on a mixture’s composition, and often leads to a nearly complete separation of the different particle phases. For granular chute flows this tends to un-mix particles of different sizes into layers, creating an inversely graded structure; the larger particles tend to segregate to the upper free surface, while the smaller particles sink to the base (Middleton 1970). Although there is usually some degree of diffusive remixing (Gray & Chugunov 2006), in the bidisperse granular chute flows we consider here, layers of pure large and pure small particles develop, with only a thin layer of mixed material in between (Savage & Lun 1988).

There are many mechanisms of size segregation in granular flows (Dippel & Luding 1995; Luding et al. 1996; Jenkins 1998; Mullin 2000; Hong, Quinn & Luding 2001; Mobius et al. 2001; Jenkins & Yoon 2002; González et al. 2015; Windows-Yule et al. 2016), but kinetic sieving and squeeze expulsion are the mechanisms that are often used to describe segregation in dense granular chute flows (Middleton 1970; Savage & Lun 1988; Vallance & Savage 2000). In this description, void spaces open up between particles, allowing for percolation downwards. Small particles are more likely to fall into these gaps, which leads to a movement of small particles towards the base of the flow. At the same time, force imbalances squeeze all particles out of their layers towards the free surface, which leads to an upwards movement. The combination results in a net migration of small particles to the base and large particles towards the free surface. More recently, alternative mechanisms for size segregation, in dense chute flows, have been proposed, based on e.g. buoyancy effects (Khakhar, McCarthy & Ottino 1999), differences in fluctuating kinetic energy (Fan & Hill 2011; Staron & Phillips 2015b ) and kinetic theory (Larcher & Jenkins 2013).

In this study, the particles not only differ in their sizes, but also in their microscopic friction coefficients; larger particles are ‘rougher’ than smaller particles, which is expressed as a larger microscopic friction coefficient. This difference leads to a larger macroscopic friction coefficient for flows consisting of larger particles than for smaller particles. This is consistent with large, geophysical particulate flows, where the smaller particles usually have a higher flowability than larger particles (Iverson 2003). It also induces segregation based on microscopic friction, where smoother particles tend towards the bottom of the flow (Gillemot, Somfai & Börzsönyi 2017; van der Vaart et al. 2017). Thus the larger, rougher particles rise towards the free surface, while the smaller, smoother particles percolate towards the bottom of the flow.

The large rough particles in the higher layers of the flow are transported towards the front, due to the higher velocity near the free surface. This creates an additional downstream (‘horizontal’) segregation structure, with the front of the flow consisting of mostly large particles. Larger particles have a larger macroscopic friction coefficient than smaller particles for the same (dimensional) flow height (Weinhart et al. 2012; Thornton et al. 2012a ; Staron & Phillips 2015a ); this causes a reduced mobility of the flow in the large particle front compared to the bidisperse tail. In our case, this effect is even larger because, as stated before, we also use a higher microscopic friction for the large particles, compared to the small particles. This configuration, with larger frictional resistive grains at the front and finer more mobile material 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).

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.

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.

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; Bokhove & Thornton 2012). 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 advection–diffusion-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 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.

2 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.

2.1 Definitions

We consider a plane at an angle $\unicode[STIX]{x1D703}$ 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, $\boldsymbol{g}$ , has the direction $(\sin \unicode[STIX]{x1D703},0,-\!\cos \unicode[STIX]{x1D703})$ 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, $\bar{u}(x,t)$ , and the depth-averaged small particle concentration, $\bar{\unicode[STIX]{x1D719}}(x,t)$ , at all positions, $x$ , at any time, $t>0$ . Here, the small particle concentration, $\unicode[STIX]{x1D719}$ , is defined as the volume fraction of the solid phase that consists of small particles, so the large particle concentration equals $(1-\unicode[STIX]{x1D719})$ . The depth-averaged velocity and small particle concentration can be derived from their full two-dimensional equivalents $u(x,z,t)$ and $\unicode[STIX]{x1D719}(x,z,t)$ and the flow height, $h(x,t)$ , as

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}(x,t)=\frac{1}{h(x,t)}\int _{0}^{h(x,t)}u(x,z,t)\,\text{d}z, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D719}}(x,t)=\frac{1}{h(x,t)}\int _{0}^{h(x,t)}\unicode[STIX]{x1D719}(x,z,t)\,\text{d}z. & \displaystyle\end{eqnarray}$$

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 variables are non-dimensionalised by the large particle diameter, $d^{\ast L}$ , and gravitational constant, $g$ , as follows: the height and length of the flow are non-dimensionalised by $d^{\ast L}$ , time by $\sqrt{d^{\ast L}/g}$ and consequently the depth-averaged velocity by $\sqrt{d^{\ast L}~g}$ . We also non-dimensionalise the particle diameter for both types of particles with $d^{\ast L}$ , and thus $d^{S}=d^{\ast S}/d^{\ast 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.

Figure 2. Schematic drawing of a bidisperse chute flow. The granular material flows down a chute inclined at an angle $\unicode[STIX]{x1D703}$ 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, $\boldsymbol{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.

2.2 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

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D719}(x,z,t)=\left\{\begin{array}{@{}ll@{}}1\quad & \text{if }0<z<h(x,t)\bar{\unicode[STIX]{x1D719}}(x,t),\\ 0\quad & \text{if }h(x,t)\bar{\unicode[STIX]{x1D719}}(x,t)\leqslant z<h(x,t),\end{array}\right.\end{eqnarray}$$

where $h\bar{\unicode[STIX]{x1D719}}$ 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

(2.4) $$\begin{eqnarray}u(z)=\unicode[STIX]{x1D6FC}_{s}\bar{u}+2(1-\unicode[STIX]{x1D6FC}_{s})\bar{u}\frac{z}{h},\quad 0\leqslant \unicode[STIX]{x1D6FC}_{s}\leqslant 1.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FC}_{s}$ is the amount of basal slip; for $\unicode[STIX]{x1D6FC}_{s}=1$ , equation (2.4) describes a plug flow, while for $\unicode[STIX]{x1D6FC}_{s}=0$ it describes a linear velocity profile with no slip at the base. The choice of the constant value of $\unicode[STIX]{x1D6FC}_{s}$ is discussed in § 2.4. Note that $\unicode[STIX]{x1D6FC}_{s}$ is often denoted by $\unicode[STIX]{x1D6FC}$ , 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 $\unicode[STIX]{x1D6FC}$ , 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, $\unicode[STIX]{x1D6FC}$ , as

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\frac{\overline{u^{2}}}{\bar{u}^{2}}=\frac{1}{3}(1-\unicode[STIX]{x1D6FC}_{s})^{2}+1.\end{eqnarray}$$

This is where our model differs from that of Woodhouse et al. (2012), who set $\unicode[STIX]{x1D6FC}=1$ independent of $\unicode[STIX]{x1D6FC}_{s}$ . We choose to retain the dependency of $\unicode[STIX]{x1D6FC}$ on $\unicode[STIX]{x1D6FC}_{s}$ to keep the model consistent. Moreover, choosing $\unicode[STIX]{x1D6FC}\neq 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 $\unicode[STIX]{x1D6FC}=1$ also shows the bulbous-head profile. However, we prefer to keep the model internally consistent and hence choose $\unicode[STIX]{x1D6FC}$ 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

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{u})=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{u})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D6FC}h\bar{u}^{2}+\frac{1}{2}h^{2}\cos \unicode[STIX]{x1D703}\right)=hS(h,\bar{u},\bar{\unicode[STIX]{x1D719}}), & \displaystyle\end{eqnarray}$$

where the source term, $S$ , consists of the down-slope component of the gravitational acceleration and the basal friction respectively:

(2.8) $$\begin{eqnarray}S=\sin \unicode[STIX]{x1D703}-\unicode[STIX]{x1D707}(h,\bar{u},\bar{\unicode[STIX]{x1D719}})\frac{\bar{u}}{|\bar{u}|}\cos \unicode[STIX]{x1D703}.\end{eqnarray}$$

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, $\unicode[STIX]{x1D707}$ , which is the ratio of the shear and normal stress at the base. Traditionally, $\unicode[STIX]{x1D707}$ 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 $\unicode[STIX]{x1D707}$ 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 )

(2.9) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{\unicode[STIX]{x1D719}})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{\unicode[STIX]{x1D719}}\bar{u}-(1-\unicode[STIX]{x1D6FC}_{s})h\bar{\unicode[STIX]{x1D719}}\bar{u}(1-\bar{\unicode[STIX]{x1D719}}))=0,\end{eqnarray}$$

where it must be noted that $h\overline{\unicode[STIX]{x1D719}u}=(h\bar{\unicode[STIX]{x1D719}}\bar{u}-(1-\unicode[STIX]{x1D6FC}_{s})h\bar{\unicode[STIX]{x1D719}}\bar{u}(1-\bar{\unicode[STIX]{x1D719}}))$ when using the linear velocity profile with basal slip (2.4).

Note, that the macroscopic friction coefficient, $\unicode[STIX]{x1D707}$ , 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, $\bar{\unicode[STIX]{x1D719}}$ . To make the macroscopic friction coefficient, $\unicode[STIX]{x1D707}$ , dependent on the depth-averaged small particle concentration, $\bar{\unicode[STIX]{x1D719}}$ , we follow the approach of Pouliquen & Vallance (1999): the total basal friction is written as a concentration-weighted sum of the macroscopic friction coefficients $\unicode[STIX]{x1D707}^{S}$ and $\unicode[STIX]{x1D707}^{L}$ for monodispersed flows of small and large particles respectively:

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\bar{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D707}^{S}+(1-\bar{\unicode[STIX]{x1D719}})\unicode[STIX]{x1D707}^{L}.\end{eqnarray}$$

This provides a simple way of increasing the frictional resistance to motion as the proportion of larger particles increases at the flow front (with $\unicode[STIX]{x1D707}^{L}>\unicode[STIX]{x1D707}^{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

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D708}}=\tan \unicode[STIX]{x1D6FF}_{1}^{\unicode[STIX]{x1D708}}+\frac{\tan \unicode[STIX]{x1D6FF}_{2}^{\unicode[STIX]{x1D708}}-\tan \unicode[STIX]{x1D6FF}_{1}^{\unicode[STIX]{x1D708}}}{{\displaystyle \frac{\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D708}}h}{A^{\unicode[STIX]{x1D708}}d^{\unicode[STIX]{x1D708}}F}}+1}\quad \unicode[STIX]{x1D708}=S,L,\end{eqnarray}$$

where the Froude number $F=|\bar{u}|/\sqrt{h\cos \unicode[STIX]{x1D703}}$ can be expressed in terms of the non-dimensional velocity, $\bar{u}$ , and height, $h$ , of the flow, and the $z$ -component of gravity, $\cos \unicode[STIX]{x1D703}$ . 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 $\unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D6FF}_{1}^{\unicode[STIX]{x1D708}}$ is the minimum chute angle required for flow, $\unicode[STIX]{x1D6FF}_{2}^{\unicode[STIX]{x1D708}}$ is the maximum chute angle for which steady flows are possible, $A^{\unicode[STIX]{x1D708}}$ is a characteristic dimensionless scale and $\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D708}}$ 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 $\unicode[STIX]{x1D707}(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.

2.4 Closure relations

In this section, closure relations for the macroscopic friction coefficient, $\unicode[STIX]{x1D707}$ , and the shape factor, $\unicode[STIX]{x1D6FC}$ , 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, $\unicode[STIX]{x1D6FC}$ , 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, $\unicode[STIX]{x1D6FC}$ 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 $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}(h^{in},\unicode[STIX]{x1D703})$ constant throughout the domain, depending only on the inflow height and the chute angle. For example, for inflow height $h^{in}=15$ and angle $\unicode[STIX]{x1D703}=24^{\circ }$ , a constant value $\unicode[STIX]{x1D6FC}\approx 1.23$ is used throughout the whole domain. Combined with (2.5), this results in $\unicode[STIX]{x1D6FC}_{s}\approx 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.

Table 1. Summary of parameters for the continuum model used in this paper. The friction parameters, $\unicode[STIX]{x1D6FF}_{1},\unicode[STIX]{x1D6FF}_{2},A$ and $\unicode[STIX]{x1D6FD}$ , depend on the granular materials used, while $\unicode[STIX]{x1D703}$ , $h^{in}$ and $\bar{\unicode[STIX]{x1D719}}^{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, $\unicode[STIX]{x1D6FC}_{s}$ , depends on both the material and the geometry.

3 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 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\in [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\bar{u}(x,0)=h\bar{\unicode[STIX]{x1D719}}(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\bar{u}(0,t)=h^{in}\bar{u}^{in}$ and $h\bar{\unicode[STIX]{x1D719}}(0,t)=h^{in}\bar{\unicode[STIX]{x1D719}}^{in}$ . These inflow values are chosen such that the flow at $x=0$ is non-accelerating, i.e. $S(h^{in},\bar{u}^{in},\bar{\unicode[STIX]{x1D719}}^{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 $\unicode[STIX]{x1D703}=24^{\circ }$ 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$ , $\bar{\unicode[STIX]{x1D719}}^{in}=0.5$ . For a flow with uniform inflow height, it then follows that $\bar{u}^{in}\approx 3.4$ , see § B.1. For DGFEM simulations, we have chosen a time step of $\unicode[STIX]{x0394}t=10^{-4}$ and $N=250$ elements. All parameter values used in the numerical solutions are summarised in table 1.

Figure 3. 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\bar{\unicode[STIX]{x1D719}}$ . Both $x$ and $z$ are scaled by the large particle diameter, $d^{\ast 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.

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 noteworthy that as a result of choosing $\unicode[STIX]{x1D6FC}\neq 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\approx 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\approx 1000$ (figure 3 c), the bulbous-head height stays constant. Near the inflow, the flow is uniform, with the prescribed height, velocity and small 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 3 c–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 $\unicode[STIX]{x1D6FC}\neq 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.

4 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 $\bar{u}$ are continuous, but there is a steady moving shock in $\bar{\unicode[STIX]{x1D719}}$ . 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

(4.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=t,\quad \unicode[STIX]{x1D709}=x-u_{f}t,\end{eqnarray}$$

such that the frame moves with speed $u_{f}$ in the positive direction with $\unicode[STIX]{x1D709}=0$ being the fixed position of the shock. Applying this coordinate transformation leads to the system

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}(h(\bar{u}-u_{f}))=0, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}(h\bar{u})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\left(h\bar{u}(\unicode[STIX]{x1D6FC}\bar{u}-u_{f})+\frac{1}{2}h^{2}\cos \unicode[STIX]{x1D703}\right)=hS(h,\bar{u},\bar{\unicode[STIX]{x1D719}}), & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}(h\bar{\unicode[STIX]{x1D719}})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}(h\bar{\unicode[STIX]{x1D719}}(\bar{u}-u_{f})-h\bar{\unicode[STIX]{x1D719}}\bar{u}(1-\unicode[STIX]{x1D6FC}_{s})(1-\bar{\unicode[STIX]{x1D719}}))=0. & \displaystyle\end{eqnarray}$$

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 $\unicode[STIX]{x1D70F}$ -derivatives to zero, so the system becomes a set of coupled ordinary differential equations:

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}\unicode[STIX]{x1D709}}(h(\bar{u}-u_{f}))=0, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}\unicode[STIX]{x1D709}}\left(h\bar{u}(\unicode[STIX]{x1D6FC}\bar{u}-u_{f})+\frac{1}{2}h^{2}\cos \unicode[STIX]{x1D703}\right)=hS(h,\bar{u},\bar{\unicode[STIX]{x1D719}}), & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}\unicode[STIX]{x1D709}}(h\bar{\unicode[STIX]{x1D719}}(\bar{u}-u_{f})-h\bar{\unicode[STIX]{x1D719}}\bar{u}(1-\unicode[STIX]{x1D6FC}_{s})(1-\bar{\unicode[STIX]{x1D719}}))=0. & \displaystyle\end{eqnarray}$$

Equations (4.5) and (4.7) can be integrated directly:

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle h(\bar{u}-u_{f})=C_{1}, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle h\bar{\unicode[STIX]{x1D719}}(\bar{u}-u_{f})-h\bar{\unicode[STIX]{x1D719}}\bar{u}(1-\unicode[STIX]{x1D6FC}_{s})(1-\bar{\unicode[STIX]{x1D719}})=C_{2}, & \displaystyle\end{eqnarray}$$

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\bar{u}$ , and small particle height, $h\bar{\unicode[STIX]{x1D719}}$ , go asymptotically to a fixed value for $\unicode[STIX]{x1D709}\rightarrow -\infty$ , so $h\rightarrow h^{in}$ , $\bar{u}\rightarrow \bar{u}^{in}$ and $\bar{\unicode[STIX]{x1D719}}\rightarrow \bar{\unicode[STIX]{x1D719}}^{in}$ as $\unicode[STIX]{x1D709}\rightarrow -\infty$ for some constants $h^{in}>0$ , $\bar{u}^{in}>0$ and $\bar{\unicode[STIX]{x1D719}}^{in}\in [0,1]$ . Using (4.8) and (4.9), the integration constants, $C_{1}$ and $C_{2}$ , are then given by

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle C_{1}=h^{in}(\bar{u}^{in}-u_{f}), & \displaystyle\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle C_{2}=h^{in}\bar{\unicode[STIX]{x1D719}}^{in}(\bar{u}^{in}-u_{f})-h^{in}\bar{\unicode[STIX]{x1D719}}^{in}\bar{u}^{in}(1-\unicode[STIX]{x1D6FC}_{s})(1-\bar{\unicode[STIX]{x1D719}}^{in}). & \displaystyle\end{eqnarray}$$

Substituting (4.10) into (4.8) gives the relation for the height depending on the depth-averaged velocity

(4.12) $$\begin{eqnarray}h(\bar{u})=\frac{h^{in}(\bar{u}^{in}-u_{f})}{\bar{u}-u_{f}}.\end{eqnarray}$$

Similarly, substitution of (4.11) into (4.9) and rearrangement of the terms shows that $\bar{\unicode[STIX]{x1D719}}$ is a solution to the quadratic equation

(4.13) $$\begin{eqnarray}(1-\unicode[STIX]{x1D6FC}_{s})\bar{u}\bar{\unicode[STIX]{x1D719}}^{2}+(\unicode[STIX]{x1D6FC}_{s}\bar{u}-u_{f})\bar{\unicode[STIX]{x1D719}}-\frac{C_{2}}{h}=0.\end{eqnarray}$$

This equation can be solved with the quadratic formula, so if $\bar{u}$ is known, then $h(\bar{u})$ follows from (4.12), and $\bar{\unicode[STIX]{x1D719}}(\bar{u})$ can be found from

(4.14) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D719}}(\bar{u})=\frac{-(\unicode[STIX]{x1D6FC}_{s}\bar{u}-u_{f})+\sqrt{(\unicode[STIX]{x1D6FC}_{s}\bar{u}-u_{f})^{2}+4(1-\unicode[STIX]{x1D6FC}_{s})\bar{u}C_{2}/h(\bar{u})}}{2(1-\unicode[STIX]{x1D6FC}_{s})\bar{u}}.\end{eqnarray}$$

Note that the other solution to the quadratic equation (4.13) gives a negative value for $\bar{\unicode[STIX]{x1D719}}$ , which is unphysical. Using these results to eliminate $h$ and $\bar{\unicode[STIX]{x1D719}}$ from (4.6) leads to a first-order nonlinear ordinary differential equation for $\bar{u}$ :

(4.15) $$\begin{eqnarray}\left(2\unicode[STIX]{x1D6FC}\bar{u}-u_{f}-\frac{\bar{u}(\unicode[STIX]{x1D6FC}\bar{u}-u_{f})}{\bar{u}-u_{f}}-\frac{h^{in}(\bar{u}^{in}-u_{f})\cos \unicode[STIX]{x1D703}}{(\bar{u}-u_{f})^{2}}\right)\frac{\text{d}\bar{u}}{\text{d}\unicode[STIX]{x1D709}}=S(h(\bar{u}),\bar{u},\bar{\unicode[STIX]{x1D719}}(\bar{u})).\end{eqnarray}$$

4.1 Boundary conditions

As in § 3, we prescribe the inflow height, $h^{in}$ , and the inflow concentration of small particles, $\bar{\unicode[STIX]{x1D719}}^{in}$ . We also want the height to go to the inflow height asymptotically for $\unicode[STIX]{x1D709}\rightarrow -\infty$ , so $(\text{d}h/\text{d}\unicode[STIX]{x1D709})^{in}=0$ is prescribed. Using (4.5) and the chain rule, it follows that

(4.16) $$\begin{eqnarray}\left(\frac{\text{d}\bar{u}}{\text{d}\unicode[STIX]{x1D709}}\right)^{in}=\frac{-C_{1}}{(h^{in})^{2}}\left(\frac{\text{d}h}{\text{d}\unicode[STIX]{x1D709}}\right)^{in}=0.\end{eqnarray}$$

Equations (4.15) and (4.16) now imply that $\bar{u}^{in}$ must satisfy $S(h^{in},\bar{u}^{in},\bar{\unicode[STIX]{x1D719}}^{in})=0$ . For the current friction law (2.10)–(2.11), this is given by:

(4.17) $$\begin{eqnarray}\bar{u}^{in}=f(\bar{\unicode[STIX]{x1D719}}^{in})(h^{in})^{3/2},\end{eqnarray}$$

where $f(\bar{\unicode[STIX]{x1D719}}^{in})$ is a function of $\bar{\unicode[STIX]{x1D719}}^{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 $\bar{u}^{out}$ at the outflow position of the domain, $\unicode[STIX]{x1D709}\rightarrow \infty$ . We look for a flow that is steady for $\unicode[STIX]{x1D709}>0$ , with no small particles in this region, so $\bar{\unicode[STIX]{x1D719}}=0$ if $\unicode[STIX]{x1D709}>0$ . To determine the values of $h^{out}$ and $\bar{u}^{out}$ , we use the global mass balance,

(4.18) $$\begin{eqnarray}h^{in}(\bar{u}^{in}-u_{f})=h^{out}(\bar{u}^{out}-u_{f}).\end{eqnarray}$$

Furthermore, we prescribe that the flow is non-accelerating and uniform with only large particles in the front:

(4.19) $$\begin{eqnarray}S(h^{out},\bar{u}^{out},0)=0.\end{eqnarray}$$

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 $\bar{u}=\bar{u}^{out}$ for $\unicode[STIX]{x1D709}>0$ . They can be computed with Matlab’s solve routine from

(4.20) $$\begin{eqnarray}(\bar{u}^{out})^{5/3}-u_{f}(\bar{u}^{out})^{2/3}-\left(\frac{\tan \unicode[STIX]{x1D703}-\tan \unicode[STIX]{x1D6FF}_{1}^{L}}{\tan \unicode[STIX]{x1D6FF}_{2}^{L}-\tan \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x1D6FD}^{L}\sqrt{\cos \unicode[STIX]{x1D703}}}{A^{L}}\right)^{2/3}h^{in}(\bar{u}^{in}-u_{f})=0.\end{eqnarray}$$

The value of $h^{out}$ then follows directly from (4.18).

4.2 Shock properties

To compute the frame speed, $u_{f}$ , consider a plane perpendicular to the chute at a certain $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{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) $$\begin{eqnarray}\int _{0}^{\bar{\unicode[STIX]{x1D719}}h}(u(z)-u_{f})\,\text{d}z=0.\end{eqnarray}$$

Evaluating this integral at $\unicode[STIX]{x1D709}\rightarrow -\infty$ with the linear velocity profile (2.4), it can be shown that the shock speed is given by

(4.22) $$\begin{eqnarray}u_{f}=\bar{u}^{in}(\unicode[STIX]{x1D6FC}_{s}+(1-\unicode[STIX]{x1D6FC}_{s})\bar{\unicode[STIX]{x1D719}}^{in}).\end{eqnarray}$$

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 $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}t+\unicode[STIX]{x2202}f(u)/\unicode[STIX]{x2202}x=0$ , the Rankine–Hugoniot jump condition states that the shock speed $s$ satisfies $s(u^{+}-u^{-})=f(u^{+})-f(u^{-})$ (LeVeque 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

(4.23) $$\begin{eqnarray}\displaystyle u_{f}((h\bar{\unicode[STIX]{x1D719}})^{+}-(h\bar{\unicode[STIX]{x1D719}})^{-}) & = & \displaystyle ((h\bar{\unicode[STIX]{x1D719}})^{+}(\bar{u})^{+}-(h\bar{\unicode[STIX]{x1D719}})^{+}(\bar{u})^{+}(1-\unicode[STIX]{x1D6FC}_{s})(1-(\bar{\unicode[STIX]{x1D719}})^{+}))\nonumber\\ \displaystyle & & \displaystyle -\,((h\bar{\unicode[STIX]{x1D719}})^{-}(\bar{u})^{-}-(h\bar{\unicode[STIX]{x1D719}})^{-}(\bar{u})^{-}(1-\unicode[STIX]{x1D6FC}_{s})(1-(\bar{\unicode[STIX]{x1D719}})^{-})),\end{eqnarray}$$

where $(\bar{\unicode[STIX]{x1D719}})^{-}=\lim _{\unicode[STIX]{x1D709}\uparrow 0}\bar{\unicode[STIX]{x1D719}}$ and $(\bar{\unicode[STIX]{x1D719}})^{+}=\lim _{\unicode[STIX]{x1D709}\downarrow 0}\bar{\unicode[STIX]{x1D719}}$ . Since $h$ and $\bar{u}$ are continuous at the shock position, it holds that $(h)^{-}=(h)^{+}=h^{out}$ and $(\bar{u})^{-}=(\bar{u})^{+}=\bar{u}^{out}$ . Working out the parentheses and substituting these relations for the height- and depth-averaged velocity in (4.23) gives

(4.24) $$\begin{eqnarray}u_{f}=\unicode[STIX]{x1D6FC}_{s}\bar{u}+\bar{u}(1-\unicode[STIX]{x1D6FC}_{s})((\bar{\unicode[STIX]{x1D719}})^{+}+(\bar{\unicode[STIX]{x1D719}})^{-}).\end{eqnarray}$$

Assuming the small particle concentration at the downstream position of the shock disappears, $(\bar{\unicode[STIX]{x1D719}})^{+}=0$ , it follows from (4.24) that

(4.25) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D719}}^{out}=(\bar{\unicode[STIX]{x1D719}})^{-}=\frac{u_{f}-\unicode[STIX]{x1D6FC}_{s}\bar{u}^{out}}{(1-\unicode[STIX]{x1D6FC}_{s})\bar{u}^{out}},\end{eqnarray}$$

is the value of the small particle concentration at the left side of the shock.

4.3 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, $\unicode[STIX]{x1D709}=0$ , are $h=h^{out}$ , $\bar{u}=\bar{u}^{out}$ and $(\bar{\unicode[STIX]{x1D719}})^{-}=\bar{\unicode[STIX]{x1D719}}^{out}$ and $(\bar{\unicode[STIX]{x1D719}})^{+}=0$ . Equation (4.15) is integrated in both positive and negative $\unicode[STIX]{x1D709}$ -direction starting from $\unicode[STIX]{x1D709}=0$ with the different initial conditions for $\bar{\unicode[STIX]{x1D719}}$ .

Figure 4. Profiles for the travelling wave solution for the height $h$ , height of small particles $h\bar{\unicode[STIX]{x1D719}}$ and depth-averaged velocity $\bar{u}$ of the flow. Note, that this solution is valid for $\unicode[STIX]{x1D709}\in (-\infty ,\infty )$ , with $h\rightarrow h^{in}$ , $\bar{\unicode[STIX]{x1D719}}\rightarrow \bar{\unicode[STIX]{x1D719}}^{in}$ and $\bar{u}\rightarrow \bar{u}^{in}$ as $\unicode[STIX]{x1D709}\rightarrow -\infty$ and $h=h^{out}$ , $\bar{\unicode[STIX]{x1D719}}=0$ and $\bar{u}=\bar{u}^{out}$ for $\unicode[STIX]{x1D709}>0$ . Therefore the total travelling wave solution has infinite mass.

The resulting profiles are displayed in figure 4. For $\unicode[STIX]{x1D709}>0$ , the flow is non-accelerating, since the height, velocity and small particle concentration stay constant, in accordance with constraint (4.19). There are no small particles present, since $\bar{\unicode[STIX]{x1D719}}=0$ for all $\unicode[STIX]{x1D709}>0$ . The flow for $\unicode[STIX]{x1D709}<0$ is more interesting, where it is important to note that $h\rightarrow h^{in}$ , $\bar{u}\rightarrow \bar{u}^{in}$ and $\bar{\unicode[STIX]{x1D719}}\rightarrow \bar{\unicode[STIX]{x1D719}}^{in}$ asymptotically as $\unicode[STIX]{x1D709}\rightarrow -\infty$ . All curves are smooth, except for the shock in $\bar{\unicode[STIX]{x1D719}}$ at $\unicode[STIX]{x1D709}=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$ , $\unicode[STIX]{x1D6FF}_{1}^{L}$ and $\unicode[STIX]{x1D6FF}_{2}^{L}$ 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 relationship, for the similar constant height breaking wave structure can be found in van der Vaart et al. (2018).

The shock speed of $u_{f}\approx 1.9$ is significantly lower than the outflow speed, $\bar{u}^{out}\approx 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

(4.26) $$\begin{eqnarray}\displaystyle u_{f} & = & \displaystyle \unicode[STIX]{x1D6FC}_{s}\bar{u}^{out}+2(1-\unicode[STIX]{x1D6FC}_{s})\frac{\bar{u}^{out}}{h^{out}}\left(\frac{h^{out}\bar{\unicode[STIX]{x1D719}}^{out}}{2}\right),\end{eqnarray}$$

(4.27) $$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D6FC}_{s}\bar{u}^{out}+(1-\unicode[STIX]{x1D6FC}_{s})\bar{u}^{out}\bar{\unicode[STIX]{x1D719}}^{out}.\end{eqnarray}$$

This is exactly the shock speed that we see, which is significantly lower than the outflow speed. This difference in $u_{f}$ and $\bar{u}^{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.

4.4 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 $\unicode[STIX]{x1D709}=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

(4.28) $$\begin{eqnarray}\left((\unicode[STIX]{x1D6FC}-1)\frac{(\bar{u}^{out})^{2}}{h\cos \unicode[STIX]{x1D703}}+1\right)\frac{\text{d}h}{\text{d}\tilde{x}}=\tan \unicode[STIX]{x1D703}-\unicode[STIX]{x1D707}^{L}(h,\bar{u}^{out}),\end{eqnarray}$$

where $\tilde{x}$ is the distance from the point where the height vanishes. This is a generalisation for $\unicode[STIX]{x1D6FC}\neq 1$ of the monodisperse front solutions derived by Pouliquen (1999a ) and Gray & Ancey (2009). By shifting $\tilde{x}$ , 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 $\unicode[STIX]{x1D709}<0$ , the time-dependent solution matches the travelling wave solution of this section. For $\unicode[STIX]{x1D709}>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 $\bar{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.

Figure 5. Numerical height and concentration profiles at $t=2500$ obtained from DGFEM solutions (black lines) compared with the travelling wave solution of figure 4 (red circles) and monodisperse front solution of (4.28) (blue diamonds). The solid lines and objects denote the height, $h$ , the dotted line and open circles denote the small particle height, $h\bar{\unicode[STIX]{x1D719}}$ . The profiles obtained from the DGFEM solutions are shifted in $x$ -direction such that the shock is located at $\unicode[STIX]{x1D709}=0$ .

5 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\in [-20,0]$ , a small periodic chute is simulated. 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\approx -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.

5.1 Height profiles

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.

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.

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^{\ast 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.

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^{\ast 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.

Figure 8. Enlarged image of the $t=500$ snapshot from figure 6. Red are the large particles, green are the small particles. The black lines denote the DGFEM solution for the height (solid) and height of small particles (dashed). Note the breaking size-segregation wave in the discrete particle simulations between $x\approx 500$ and $x\approx 1500$ .

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 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$ , $\unicode[STIX]{x1D6FC}_{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 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 over-predicts 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 steady-flow 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.

5.2 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 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.

6 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, $\unicode[STIX]{x1D6FC}\neq 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 $\unicode[STIX]{x1D6FC}=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 Gray & Ancey (2011) 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 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.