Closure Relations for Shallow Granular Flows from Particle Simulations

The Discrete Particle Method (DPM) is used to model granular flows down an inclined chute. We observe three major regimes: static piles, steady uniform flows and accelerating flows. For flows over a smooth base, other (quasi-steady) regimes are observed where the flow is either highly energetic and strongly layered in depth for small inclinations, or non-uniform and oscillating for larger inclinations. For steady uniform flows, depth profiles of density, velocity and stress have been obtained using an improved coarse-graining method, which allows accurate statistics even at the base of the flow. A shallow-layer model for granular flows is completed with macro-scale closure relations obtained from micro-scale DPM simulations of steady flows. We thus obtain relations for the effective basal friction, shape factor, mean density, and the normal stress anisotropy as functions of layer thickness, flow velocity and basal roughness. For collisional flows, the functional dependencies are well determined and have been obtained.

of 1000 m 3 of material; to the flow of sinter, pellets and coke into a blast furnace for iron-ore melting; down to the flow of sand in an hour-glass. The dynamics of these flows are influenced by many factors such as: polydisperity; variations in density; non-uniform shape; complex basal topography; surface contact properties; coexistence of static, steady and accelerating material; and, flow obstacles and constrictions.
Discrete Particle Methods (DPMs) are an extremely powerful way to investigate the effects of these and other factors. With the rapid recent improvement in computational power the full simulation of the flow in a small hour glass of millions of particles is now feasible. However, complete DPM simulations of large-scale geophysical mass flow will, probably, never be possible.
One primary goal of the present research is to simulate large scale and complex industrial flows using granular shallow-layer equations. In this paper we will take the first step of using the DPM [CS79,SEG + 01,SGPL02,SLG03, Lud08] to simulate small granular flows of mono-dispersed spherical particles in steady flow situations. We will use a refined and novel analysis to investigate three particular aspects of shallow chute flows: i) how to obtain meaningful macro-scale fields from the DPM simulation, ii) how to asses the flow dependence on the basal roughness, and iii) how to validate the assumptions made in depth-averaged theory.
The DPM simulations presented here will enable the construction of the mapping between the micro-scale and macroscale variables and functions, thus enabling construction of a closed set of continuum equations. These mappings (closure relations) can then be used in continuum shallow-layer models and compared with full DPM simulations (DPMs). For certain situations, precomputed closures should work; but, in very complicated situations the pre-established relations may fail. Heterogeneous, multi-scale modelling (HMM) is then an alternative [WEL + 07], in which the local consititute relations are directly used in the continuum model. In HMM, continuum and micro-scale models are dynamically coupled with a two-way communication between the different models in selective regions in both space and time, thus reducing computational expense and allowing simulation of complex granular flows.
1.2 Shallow-layer models Shallow granular continuum models are often used to simulate geophysical mass flows, including snow avalanches [CGJ07], dense pyroclastic flows, debris flows [DI01], block and ash flows [DPP + 08] and lahars [WSS08]. Such shallowlayer models involve approximations reducing the properties of a huge number of individual particles to a handful of averaged quantities. Originally these models were derived from the general continuum incompressible mass and momentum equations, using the long-wave approximation [SH89, HSSN93,GWH99,WGH99,GTN03,BT12] for shallow variations in the flow height and basal topography. Despite the massive reduction in degrees of freedom made, shallow-layer models tend to be surprisingly accurate, and are thus an effective tool in modelling geophysical flows. Moreover, they are now used as a geological risk assessment and hazard planning tool [DPP + 08]. In addition to these geological applications, model applications involve smallscale laboratory chute flows containing obstacles [GTN03], wedges [HH05,GC07] and contractions [VATK + 07], showing good quantitative agreement between theory and experiment.
In fluid dynamics, the Navier-Stokes equations are established with full constitutive equations. Nonetheless, the shallow-layer equations or Saint-Venant equations are often used in large scale situations where it is unpractical to solve the full Navier-Stokes equations. Our present aim is to directly investigate the validity of the assumptions of granular shallow-layer models first from discrete particle simulations, before obtaining fully 3D 'kinetic theory'-style constitutive relations and simplifying these via the depth-integration process. A discussion of the full three-dimensional properties of our particle simulations will be undertaken later.
Here we restrict attention to the closures required for twodimensional shallow granular flow equations.
A key difference between shallow-layer fluid models and granular ones is the appearance of a basal friction coefficient, µ, being the ratio of the shear to normal traction at the base. In early models, a dry Coulomb-like friction law was used [SH89]. It implies µ to be constant, given by the tangent of the friction angle between the material and the base, δ , i.e., µ = tan δ . As a consequence constant uniform flow is only possible in such a model at that angle δ , independent of height. There is a considerable amount of experimental evidence, e.g., [DD99,GDR04], that suggests that such a simple Coulomb law does not hold on rough beds or for moderate inclination angles. Furthermore, detailed experimental investigations using glass beads [Pou99] lead to an improved empirical 'Pouliquen' friction law characterised by two angles: the angle at which the material comes to rest, δ 1 , where friction dominates over gravity and the angle, δ 2 , above which gravity dominates over friction and the material accelerates. Between these two angles steady flow is possible, and in the limit δ 1 → δ 2 = δ the original Coulomb style model is recovered.
Since its formulation a lot of work has been performed on extending and understanding this Pouliquen law. The original law was obtained by retarding flowing material and measuring the angle at which the material stopped as a function of height h stop (θ ), or equivalently, by inverting this relation, θ stop (h). For most materials, granular included, a greater angle is required to initiate stationary than to retard flowing material. Pouliquen and Forterre [PF02], by measuring the angle required to start motion, measured θ start (h), i.e., the friction law for initially stationary material. As expected θ start was greater than θ stop and this information was used to extend the friction law to all values of the height and velocity within the steady regime. Borzsonyi & Ecke performed a series of additional experiments: firstly, in [BE06] they looked at higher angles were the mean flow rates are close to the terminal velocity of a single particle, and found regions were the Pouliquen law is not valid and the Froude number becomes inversely proportional to the height, as opposed to the linear relationship observed for steady flow. Borzsonyi and Ecke, and Pouliquen and Forterre [BE07,FP03] have all worked on extending the original law to be valid for more complicated non-spherical materials like sand and metallic materials. Also, the effect of basal surface roughness has been systemically studied in [GTDD03] by varying the size of both the free flow and fixed basal particles. For convenience, we define λ to be the size ratio of the fixed and the free particles. They observed a peak in roughness at a certain diameter ratio, λ c , which depends on the compactness of the basal layer. Measured values of λ c in [GTDD03] ranged between 1 and 3 for a monolayer of fixed particles. For fixed particles with smaller size and λ < λ c , the range of angles where steady flow was observed decreased, and eventually the steady flow regime completely vanished, i.e., δ 2 − δ 1 → 0 as λ → 0 (yielding Coulomb type behaviour). For smaller flow particle diameters, i.e., with λ > λ c , there was also a reduction in friction, but weaker than in the small λ case. For much larger λ , the friction saturated to a constant value, which they contributed to free particles that filled the holes in the basal surface and effectively created a stable basal surface of free particles. In a later publication [GDDT07], they extended this investigation to flows containing two particle sizes.
Louge and Keast [LK01] modified the kinetic theory presented in [Jen93] by modelling enduring contacts via a fric-tional rate-independent stress component in order to obtain steady flow on flat frictional inclines. This work was later extended to bumpy inclines [Lou03]. Jenkins [Jen06] took a different approach and theoretically formulated a phenomenological modification of granular kinetic theory to account for enduring particle contacts. His idea is that enduring contacts between grains, forced by the shearing, reduce the collision rate of dissipation. Therefore a modification to the dissipation is introduced, which does not affect the stress. It leads to a law very similar to the one experimentally obtained by Pouliquen. He extended the theory in [Jen07] to very dissipative frictional particles, with a coefficient of restitution less than 0.7. Later, a detailed comparison with new experiments was performed, showing agreement with flows of low inclination [JB10].
Silbert et al.
[SEG + 01] used DPMs to simulated chute flow of cohesionless particles. They found that a steady-state flow regime exists over a wide range of inclination angles, heights and interaction parameters, in confirmation of the experiments of Pouliquen [Pou99]. They found for steadystate flows that the volume fraction is constant throughout the flow, in agreement with the assumptions of shallow-layer theory [SH89]. They also observed that the shear stress is proportional to the square of the shear and the flow velocity scales with the height to the power 3/2. This result coincides with Bagnold's analysis of dilute binary collisions flows [Bag54]. They also observed small systematic deviations from isotropic stress, which shows a deviation from fluid-like behaviour. However, normal stresses do not approach a Coulomb-yield criterion structure at the angle of repose except near the surface, hinting that the failure of flow starts near the surface. In [SGPL02], they looked at the effect of different basal types and found that for an ordered chute base the steady state regime splits into three distinct flow regimes: at smaller angles, the flowing system self-organizes into a state of low-dissipation flow consisting of in-plane ordering in the bulk; at higher angles, a high-dissipation regime similar to that for a rough base but with considerable slip at the bottom is observed; and, between these two subregions they observe a transitional flow regime characterized by large oscillations in the bulk averaged kinetic energy due to the spontaneous ordering and disordering of the system as a function of time. Finally, [SLG03] investigated the initiation and cessation of granular chute flow more carefully and computed both θ stop and θ start . For inclinations θ ≫ θ stop they observed a Bagnold rheology, for θ > ∼ θ stop a linear profile, and for θ ≈ θ stop avalanching behaviour.

Overview of this study
Our present research is novel on the following three counts: Firstly, we compute more meaningful macro-scale fields from the DPM simulations than before by carefully choos-ing the coarse graining function. In order to homogenize the DPM data, the micro-scale fields need to be coarse-grained to obtain macroscopic fields. Coarse-grained micro-scale fields of density, momentum and stress were derived directly from the mass and momentum balance equations by Goldhirsch [Gol10]. The quality of the statistics involved depends on the coarse graining width w, which defines the amount of spatial smoothing. For small coarse-graining width w, micro-scale variations remain visible, while for large w these smooth out in the macro-scale gradients. Since one of the objectives is to obtain the value of µ at the base, we use a novel adaptation of Goldhirsch' statistics [Gol10] near boundaries. This is a non-trivial issue neglected in the literature; often continuum fields are simply not computed within a distance w of the boundary. Our new approach [WTLB11] is consistent with the continuum equations everywhere, enabling the construction of continuum fields within one course-graining width of the boundary.
Secondly, we follow the approach of [GTDD03] and vary the basal particle diameter to achieve different basal conditions. For particles with smaller basal than flowing diameter, λ < 1, the flow becomes more energetic and oscillatory behaviour is observed. This phenomena has previously been reported in [SGPL02], but was achieved by changing the basal particles to a more regular, grid-like configuration. By investigating flow over fixed particles of different size than the free ones, we are able to quantify the roughness and numerically investigate the transition from rough to smooth surfaces. For smoother surfaces, we show that the parameter space can be split into to two types of steady flow, and we obtain a general friction law.
Finally, we test the assumptions made in depth-averaged theory and determine the required closure laws. For shallow granular flows, the flow can be described by depth-averaged mass and momentum-balance equations [GTN03]. Solving the depth-averaged equations requires a constitutive relation for the basal friction, a way to account for mean density variations, the shape of the velocity profile and the pressure anisotropy. We extract such data from DPMs obtained for steady uniform flows, and establish a novel, extended set of closure equations. Also, the depth-averaged equations are obtained under the assumptions that a) the density is constant in space and time and does not vary through the flow; b) the ratio between mean squared velocity and the squared mean velocity is known; c) the downward normal stress is lithostatic, i.e., balances the gravitational forces acting on the flow; and, d) the ratio between the normal stresses is known. Gray et al. [GTN03] assumed the latter ratio to be one. The depth profiles of these quantities are discussed by Silbert et al. in [SEG + 01,SGPL02,SLG03] for steady flow. We extend these measurements to validate our DPM simulations, using our superior statistical procedure. Hence, we establish the range in which the shallow-layer model is valid, and the closure relations required for shallow-layer continuum simulations, for a much wider range of flow regimes than had been considered before.

Outline
We introduce the force model used in the DPM in Sect. 2, and the statistical method used to obtain macroscopic density, velocity and stress profiles in Sect. 3. In Sect. 4, we discuss the continuum shallow-layer equations for modelling granular flow including some macro-scale closures. The set up of the simulations is discussed in Sect. 5, and the steadystate regime is mapped for flows over a rough basal surface in Sect. 6. Depth profiles of the flow are introduced in Sect. 7, which are then used to characterize the steady flow over smoother surfaces in Sect. 8. Finally, the closure relations for the shallow-layer model are established in Sect. 9, before we conclude in Sect. 10.

Contact law description
A Discrete Particle Method (DPM) is used to perform the simulation of a collection of identical granular particles. Boundaries are created by special fixed particles, which generally will have different properties than the flow particles. Particles interact by the standard spring-dashpot interaction model [CS79,WB86,Lud08], in which it is assumed that particles are spherical and soft, and that pairs have a single contact point.
Each particle i has a diameter d i , density ρ i , position r i , velocity v i and angular velocity ω i . For pairs of two particles {i, j}, we define the relative distance vector r i j = r i − r j , their separation r i j = |r i j |, the unit normaln i j = r i j /r i j and the relative velocity v i j = v i − v j . Two particles are in contact if their overlap, is positive. A single contact point c at the centre of the overlap is assumed, which is a valid assumption as long as the overlap is small. For our simulations the overlap between two particles is always below 1% of the particle radius; the simple contact model we employ thus captures the key process as there are no multiple contact points. The force acting on particle i is a combination of the pairwise interaction of two particles. The force f i j represents the force on particle i from the interaction with particle j and can be decomposed into a normal and a tangential component, We assume that the particles are viscoelastic; therefore, they experience elastic as well as dissipative forces in both normal and tangential directions. The normal force is modelled as a spring-dashpot with a linear elastic and a linear dissipative contribution, with spring constant k n , damping coefficient γ n and the normal relative velocity component, For a central collision, no tangential forces are present, and the collision time t c between two particles can be calculated as with the reduced mass m i j = m i m j /(m i + m j ). The normal restitution coefficient r c (ratio of relative normal speed after and before collision) is calculated as We also assume a linear elastic and a linear dissipative force in the tangential direction, with spring constant k t , damping coefficient γ t , elastic tangential displacement δ δ δ t i j (which is explained later), and the relative tangential velocity at the contact point, with b i j = − (d i − δ n i j )/2 n i j the branch vector from point i to the contact point; for equal size particles b i j = −r i j /2.
The elastic tangential force is modelled and derives from the roughness of the particle surface. Near the contact point, small bumps on a real particle would stick to each other, due to the normal force pressing them together, and elongate in the tangential direction resulting in an elastic force proportional to the elastic tangential displacement. It is defined to be zero at the initial time of contact, and its rate of change is given by where the first term is the relative tangential velocity at the contact point, and the second term ensures that δ δ δ t i j remains normal ton i j . The second term is always orthogonal to the spring direction and, hence, does not affect the rate of change of the spring length: it simply rotates it, thus keeping it tangential.
When the tangential to normal force ratio becomes larger than the particle contact friction coefficient, µ c , for a real particle the bumps would slip against each other and their elongation is shortened until the bumps can stick to each other again. This is modelled by a static yield criterion, truncating the magnitude of δ δ δ t i j as necessary to satisfy |f t i j | ≤ µ c |f ,n i j |. Thus, the contact surfaces are treated as stuck while |f t i j | < µ c |f n i j | and as slipping otherwise, when the yield criterion is satisfied 1 . The total force on particle i is a combination of contact forces f i j with other particles and external forces such as gravity g. The resulting force f i and torque q i acting on particle i are (10) Finally, using these expressions we arrive at Newton's equations of motion for the translational and rotational degrees of freedom, with m i the mass and I i the inertia of particle i. We integrate (11) forward using Velocity-Verlet [AT93], formally second order in time, with an adequate time step of ∆t = t c /50. The collision time t c is given by (5), while (9) is integrated using first-order forward Euler. Hereafter, we distinguish between identical free flowing and identical fixed basal particles. Base particles are modelled as having an infinite mass and are unaffected by body forces: they do not move. This leaves two distinct types of collision: flow-flow, and flow-base. Model parameters for each of these collision types are changed independently.

Coarse-graining
The main aims of this paper are to use discrete particle simulations to both confirm the assumptions of and provide closure rules required for the depth-averaged shallow-water equations. Hence, continuum fields have to be extracted from the discrete particle data. There are many papers in the literature on how to go from the discrete to the continuum, binning micro-scale fields into small volumes [IK50,SH82, Lud04,Lud09,LLV + 01], averaging along planes [TED95], or coarse graining spatially and temporally [Bab97,SA04, Gol10]. Here, this will be achieved using a coarse-graining approach first described in [Bab97] based on a later derivation of Goldhirsch [Gol10], extended further by us [WTLB11] to account for boundary forces due to the presence of the base.
The method has several advantages over other methods because: (i) the fields produced automatically satisfy the equations of continuum mechanics, even near the flow base; (ii) it is neither assumed that the particles are rigid nor spherical; and, (iii) the results are even valid for single particles as no averaging over groups of particles is required. The only assumptions are that each particle pair has a single point of contact (i.e., the particle shapes are convex), the contact area can be replaced by a contact point (i.e., the particles are not too soft), and that collisions are not instantaneous.

Notation and basic ideas
Vectorial and tensorial components are denoted by Greek letters in order to distinguish them from the Latin particle indices i, j. Bold vector notation will be used when convenient.
Assume a system given by N p flowing particles and N b fixed particles. Since we are interested in the flow, we will calculate macroscopic fields pertaining to the flowing particles only. From statistical mechanics, the microscopic mass density of the flow, ρ mic , at a point r at time t is defined by where δ (r) is the Dirac delta function and m i is the mass of particle i. The following definition of the macroscopic density of the flow is used thus replacing the Dirac delta function in (12) by an integrable 'coarse-graining' function W whose integral over space is unity. We will take the coarse-graining function to be a Gaussian with width or variance w. Other choices of the coarse-graining function are possible, but the Gaussian has the advantage that it produces smooth fields and the required integrals can be analyzed exactly. According to [Gol10], the answer depends only weakly on the choice of function, and the width w is the key parameter. It is clear that as w → 0 the macroscopic density defined in (14) reduces to the one in (13). The coarse-graining function can also be seen as a convolution integral between the micro and macro definitions,

Mass balance
Next we will consider how to obtain the other fields of interest: the momentum vector field and the stress tensor. As stated in Sect. 3.1 the macroscopic variables will be defined in a way compatible with the continuum conservation laws. The coarse grained momentum density p(r,t) is defined by where the v iα 's are the velocity components of particle i. The macroscopic velocity field V (r,t) is then defined as the ratio of momentum and density fields, It is straightforward to confirm that equations (13) and (16) lead to the continuity equation with the Einstein summation convention for Greek letters.

Momentum balance
Finally, we will consider the momentum conservation equation with the aim of establishing the macroscopic stress field. In general, the desired momentum balance equations are written as, where σ αβ is the stress tensor, and g α a component of the acceleration vector of gravity.
Expressions (16) and (17) for the momentum p and the velocity V have already been defined. The next step is to compute their temporal and spatial derivatives, respectively, and reach closure. Taking the time derivative of (16) gives Using (11), the first term in (20) can be expressed as In the simulations presented later the force on each particle contains three contributions: particle-particle interactions, particle-base interactions, and the gravitational body force. Hence, where f i j is the interaction force between particle i and j, and f b ik the interaction between particle i and base particle k, or base wall if the base is flat. Therefore, we rework (21) as (23) can be simplified to ρg α by using (13). From Newton's third law, the contact forces are equal and opposite, such that where in the first step we interchanged the dummy summation indices. It follows from (24) that (23) can be written as Next, we will write A α as the divergence of a tensor in order to obtain a formula for the stress tensor. The following identity holds for any smooth function W where r i j = r i − r j ; we used the chain rule and differentiation to the full argument of W (·) per component. The next step extends Goldhirsch' analysis near a boundary. To obtain a similar expression for the interaction with base particles, we write which holds because W i decays towards infinity. Substituting identities (26), (27) and (13)

into (25) leads to
From [Gol10], it follows that the second term in (20) can be expressed as follows where v ′ i is the fluctuating velocity of particle i, with components given by Substituting (28) and (29) into momentum balance (19) yields Therefore the stress is given by The kinetic component of the stress tensor as well as the contact stress coming from normal forces are symmetric; only the contact stress from tangential forces can contribute to the antisymmetric part of the stress tensor. In our simulations the tangential forces contribute less than 5% to the total stress in the system, such that the stress is almost symmetric. Equation (32) differs from the results of [Gol10] by an additional term that accounts for the stress created by the presence of the base. The contribution to the stress from the interaction of two flow particles i, j is spatially distributed along the contact line from r i to r j , while the contribution from the interaction of particles i with a fixed particle k is distributed along the line from r i to r k , extending further beyond r k . We explain the situation as follows, see Fig. 1. Stress and density profiles are calculated using (15) and (32) Fig. 1: Stress and density profiles are shown for two onedimensional two-particle systems, each with two particles of unit mass at positions x = ±1, and repelling each other (so with d > 2 for our granular case). In the top figure, both particles are flowing, while in the bottom figure the left particle is fixed and the right one flowing.
for two one-dimensional two-particle systems, each with two particles of unit mass at positions x = ±1, repelling each other with a force | f | = 1 and with w = 0.2. In the top figure, both particles belong to the flowing species, so the density is distributed around the particles' center of mass and the stress along the contact line. In the bottom figure, the left particle is a fixed base particle and the right particle is a free flowing one, so density is distributed around the flowing particle's center of mass and the stress along the line extending from the flowing particle to negative infinity.
The strength of this statistical method is that the spatial coarse graining satisfies the mass and momentum balance equations exactly at any given time, irrespective of the choice of the coarse graining function. Further details about the accuracy of the stress definition (32) are discussed in [WTLB11]. The expression for the energy is also not treated in this publication, see [Bab97].

Mathematical background
In this section, we briefly outline the existing knowledge to develop a continuum shallow-layer theory for granular flow.

Shallow-layer model
Shallow-layer models have been shown to be an effective tool in modelling many geophysical mass flows. Early avalanche models were formulated by adding gravitational acceleration and Coulomb basal friction to shallow-layer models [GEY67,KE73]. Similar dry granular models have been derived using the long-wave approximation [SH89, HSSN93,Ive97,GWH99,WGH99] for shallow variations in the flow height and slope topography and included a Mohr-Coulomb rheology by the use of an earth pressure coefficient. The key to these theories is to depth-integrate general three-dimensional equations in the shallow direction, resulting in a system of two-dimensional equations which still retains some information about variations in thickness.
Let Oxyz be a coordinate system with the x-axis downslope and the z-axis normal to a channel with mean slope θ . For simplicity, we further consider boundaries, flows, and external forcing to be (statistically) uniform in y. The continuum macro-scale fields are thus independent of y, while the DPM simulations remain three-dimensional and will be periodic in y. The free-surface and base location are z = s(x,t) and z = b(x), respectively. The thickness of the flow is thus , and the bulk density and velocity components are ρ and u = (u, v = 0, w) t , respectively, as functions of x, y, z and t.
The three-dimensional flow viewed as continuum is described by the mass and momentum balance equations (18) and (19). At the top and bottom surface, kinetic boundary conditions are satisfied: D(z − s)/Dt = 0 and D(z − b)/Dt = 0 at their respective surfaces, and with material time derivative . Furthermore, the top surface is traction-free, while the traction at the basal surface is essentially Coulomb-like. We decompose the traction t = t t + t nn in tangential and normal components, with normal traction t n = −n · (σn), wheren is the outward normal at the fixed base. The Coulomb Ansatz implies that t t = −µ|t n |u/|u| with friction factor µ > 0. Note that µ generally can be a function of the local depth and the flow velocity. Its determination is essential to find a closed system of shallow-layer equations.
We consider flows that are shallow, such that a typical aspect ratio ε between flow height and length, normal and alongslope velocity, or normal and alongslope variations in basal topography, is small, of order O(ε). Furthermore, the typical friction factor µ is small enough to satisfy µ = O(ε γ ) with γ ∈ (0, 1). We follow the derivation of the depth-averaged swallow layer equations for granular flow by [GTN03] without assuming that the flow is incompressible. Instead we start the asymptotic analysis from the dimensionless form of the mass and momentum conservation equations (18) and (19), assuming only that the density is independent of depth at leading order. Density, velocity, and stress are depth averaged as follows In the end, one retains the normal stress ratio K =σ xx /σ zz , the shape factor α = u 2 /ū 2 , and the friction µ as unknowns.
The goal is to investigate whether these unknowns can be expressed as either constants or functions of the remaining shallow flow variables, to leading order in O(ε). The latter variables are the flow thickness h = h(x,t) and the depthaveraged velocityū =ū(x,t). At leading order, the momentum equation normal to the base yields that the downward normal stress is lithostatic, Depth-averaging the remaining equations, while retaining only terms of order O(ε µ), yields the dimensional depthaveraged shallow-layer equations, cf., [VATK + 07,BT12], To demarcate the dimensional time and spatial scales, we have used starred coordinates. These scales differ from the ones used before in the particle dynamics and the dimensionless ones used later in the DPM simulations. The shallowlayer equations (34) consist of the continuity equation (34a) and the downslope momentum equation (34b). The system arises also via a straightforward control volume analysis of a column of granular material viewed as continuum from base to the free surface, using Reynold-stress averaging and a leading order closure with depth averages. While the mean densityρ can be modelled as a system variable by considering the energy balance equation, we will assume that it can be expressed as a function of height and velocityρ(h,ū). Thus, the closure to equations (34) is determined when we can find the functionsρ(h,ū), K(h,ū), α(h,ū), and µ(h,ū). In Section 9.2, we will analyze if and when DPM simulations can determine these functions.

Granular friction laws for a rough basal surface
The friction coefficient, µ, was originally [HM84] taken to be a simple Coulomb type µ = tan δ , where δ is a fixed friction angle. Note that in steady state for a flat base with b = 0, the shallow-layer momentum equation (34b) then yields µ = tan θ . Pure Coulomb friction implies that there is only one inclination, θ = δ , at which steady flow of constant height and flow velocity exists. That turns out to be unrealistic. Three parameterizations for µ have been proposed in the literature.
Firstly, Forterre and Pouliquen [FP03] found steady flow in laboratory investigations for a range of inclinations concering flow over rough basal surfaces. They measured the thickness h stop of stationary material, left behind when a flowing layer was brought to rest, with the following fit where δ 1 is the minimum angle required for flow, δ 2 the maximum angle at which steady uniform flow is possible, d the particle diameter, and A a characteristic dimensionless length scale over which the friction varies. Note that h stop diverges for θ = δ 1 and is zero for θ = δ 2 . For h > h stop , steady flow exists in which the Froude number, the aspect ratio between flow speed and surface gravity-wave speed , is a linear function of the height, where β and γ are constants independent of chute inclination and particle size. Provided one uses the steady state assumption µ = tan θ to hold (approximately) in the dynamic case as well, one can combine it with (35) and (36) to find an improved empirical friction law It is a closure for µ in terms of the flow variables, which has been shown its practical value. Note that in the limit δ 1 → δ 2 = δ the Coulomb model is recovered.
Secondly, in an earlier version [Pou99], another, exponential fitting was proposed for h stop , as follows with the same limiting behaviour, and primes used to denote the difference in the fit. It yields the friction factor Equation (35) did, however, prove to be a better fit to experiments and is computationally easier to evaluate. Finally, an additional correction to (36) was made in [Jen06] to include dependence on the inclination accounting for enduring particle contacts in very dissipative frictional flows, for which we can use any appropriate fit for h stop . It leads subsequently to a more complicated evaluation of the friction law for µ. We omit further details. We will compare our DPM simulations against these rules and fits for the rough basal surface, and extend it to smoother surfaces in the upcoming sections.

Simulation description
In this section, DPM simulations are used to simulate monodispersed granular flows. Parameters have been nondimensionalised such that the flow particle diameter d = 1, mass m = 1 and the magnitude of gravity g = 1. The normal spring and damping constants are k n = 2 · 10 5 and γ n = 50; thus the contact duration is t c = 0.005 and the coefficient of restitution is ε = 0.88. The tangential spring and damping constants are k t = (2/7)k n and γ t = γ n , such that the frequency of normal and tangential contact oscillation and the normal and tangential dissipation are equal. The microscopic friction coefficient was taken to be µ c = 1/2. The interaction parameters are chosen as in [SEG + 01] to simulate glass particles of 0.1 mm size; this corresponds to a dimensional time scale of d/g = 3.1 ms and dimensional velocity scale √ dg = 0.031 ms −1 . The above parameters are identical to the simulations of Silbert et al., except that dissipation in tangential direction, γ t , was added to damp rotational degrees of freedom in arresting flow. Adding of such tangential damping removes all vibrational energy for flows otherwise arrested. The differences in the simulations with γ t = γ n reported here are minor relative to the case with γ t = 0. Silbert et al. also investigated the sensitivity of the results to the particle interaction parameters t c , ε, the ra- tio k n /k t , and µ c ; they found that while the density of the bulk material is not sensitive to these interaction parameters, the flow velocity increased with decreasing friction µ c . Nonetheless, the qualitative behaviour of the velocity profiles did not change.
The chute is periodic and of size 20 × 10 in the x-and y-directions and has a layer of fixed particles as a base. The bottom particles are monodispersed with (nondimensional) diameter λ . Various basal roughnesses are investigated by taking λ = 0, 1/2 to 2 in turn, with λ = 0 as flat base. This bottom particle layer is obtained by performing a simulation on a horizontal, smooth-bottom chute. It is filled with a randomly distributed set of particles and we simulate until a static layer about 12 particles thick is produced. Then the layers of particles at height z ∈ [9.3, 11]λ are selected, remaining particles deleted, and the selected ones moved downwards by 11λ . The layer is thick enough to ensure that no flowing particles can fall through the base. Their positions are fixed.
Defining a 'filling height' H, an integer amount of about 200 H particles are inserted into the chute. To insert a particle, a random location (x, y, z) ∈ [0, 20] × [0, 10] × [0, I] is chosen, where I = H. If the particle at this position overlaps other particles, the insertion is rejected, and the insertion domain is enlarged by increasing I to I + 0.01 to ensure that there is enough space for all particles. The initial packing fraction is about ρ/ρ p = 0.3. Thus the particles initially compact to an approximated height H, giving the chute enough kinetic energy to initialize flow. Dimensionless time is integrated between t ∈ [0, 2000] to allow the system to reach a steady state. A screen shot of the system in steady state is given in Fig. 2.
To ensure that the size of the periodic box does not influence the result, we compared density and velocity profiles of the flow at an angle θ = 24 • and height H = 17.5 for domain sizes L x = 10, 20, 40, L y = 10 and L x = 10, 20, 40, L y = L x /2, and found no significant changes.

Arrested, steady, and accelerating flow
From the experiments of Pouliquen [Pou99], granular flow over a rough base is known to exist for a range of heights and inclinations. DPMs by [SEG + 01] also showed that steady flows arose for a range of flow heights and (depth-averaged) velocities or Froude numbers. Their simulations did, however, provide relatively few data points near the boundary of arrested and steady flow to allow a more adequate fit of the stopping height. In this section, we therefore perform numerous DPMs at heights and angles near the separatrix between the steady flow regime and the regime with static piles. To study the full range of steady flow regimes, simulations were performed for inclinations θ varying between 20 • and 60 • and approximated or 'filling heights' H = 10, 20, 30 and 40. In Section 8, we will repeat (some of) these simulations for varying base roughness.
Whether a steady flow regime has been reached, is quantified here by assessing the time whereafter the ratio of ki-netic energy normalized by the mean elastic potential energy becomes time independent. This is shown in Fig. 3, where we plot such an energy ratio for a rough base, constant height, and varying chute angle. The elastic potential energy is averaged over t ∈ [1000, 2000] to minimize fluctuations after start-up, but any interval larger than 100 appears sufficient. For chute angles at most 20.5 • the kinetic energy vanishes after a short time, thus the flow arrests; for chute angles between 21 • − 29 • , a constant value is reached, indicating steady flow; and, for inclinations above 29 • the energy keeps increasing: thus flow steadily accelerates. If the energy ratio remained constant within t ∈ [1800, 2000], the flow was deemed steady, otherwise the flow was deemed to be either accelerating or stopping.
Unlike fluids, the free surface of granular flows, and thus the flow height, are not well defined. In [SEG + 01], the height of the flow was estimated by H, which is equivalent to assuming a constant packing fraction of ρ/ρ p = π/6. However, the exact height h = s − b of the flow varies from the approximated height H due to compaction of the flow and H is typically an overestimate. In [VATK + 07], the surface of the flow was defined by the time-average of the maximum vertical position of all flow particles. One could also define the free surface of the flow as the height where the density vanishes. The latter two methods, however, have the disadvantage that saltating particles can lead to slightly overestimated flow heights.
Instead, we will define the height via the downward normal stress. For steady uniform flows the downward normal stress is lithostatic, i.e., balances the gravitational weight, such that This is a direct consequence of the momentum balance equations. Thus, σ zz (z) has to decrease monotonically; the base and free surface are the heights at which σ zz (z) reaches its maximum and minimum value, respectively. However, in order to avoid effects of coarse graining or single particles near the boundary, we cut off the stress σ zz (z) on either boundary by defining threshold heights with κ = 1%. We subsequently linearly extrapolate the stress profile in the interval (z 1 , z 2 ) to define the base b and surface height s as the height at which the linear extrapoation reaches the maximum and minimum values of σ zz , respectively, The variable most sensitive to these height choices isρ. It shows well-defined functional behaviour for our definition of height, shown later. This is not the case if we define height by the density or the method in [VATK + 07]. The threshold κ = 1% was chosen because the results in Fig. 12 were relatively insensitive to the choice of κ at or above 1%.
To determine the demarcation line h s (θ ; λ ) between arrested and steady flow with good accuracy, we performed a set of simulations with initial conditions determined by the following algorithm. Starting from an initial 'filling height' H = 4 and inclination θ = 21 • , the angle was increased in steps of 1 • until eventually a flowing state was reached. From this initial flowing state, the height was increased by 2 particle heights, if the flow arrested, or else the angle decreased by 1/2 • , assuming the curve is monotonically decreasing. Flow was defined to be arrested when the energy ratio E kin /E ela fell below 10 −5 within 500 time units, otherwise the flow was classified as flowing. In simulations in which such arrested flows were continued after t = 500, a further decrease of kinetic energy was observed, thus validating the approach. This procedure yields intervals of the inclination angle for each height and, vice versa, height intervals for each angle, between which the demarcation line lies. The values presented in [SEG + 01] deviate at most 0.5 • from our observations, perhaps due to the preparation of the chute bottom, or the slightly different dissipation used. A demarcating curve between steady and arrested flow was fitted to equations (35) and (38) by minimizing the horizontal, respectively vertical, distance of the fit to these intervals, see Fig. 4. Fitting h stop (θ ) yields better results than h ′ stop (θ ) for all roughnesses and only the fit (35) will be used hereafter. Similar fits will be made in Section 8 for varying basal roughness. It leads us to a study of the depth profiles for steady state flow in the following section.

Statistics for uniform steady flow
To obtain detailed information about steady flows, we use the statistics defined in Sect. 3. Since the flows of interest are steady and uniform in x and y, density, velocity and stress will be averaged over x, y and t. The resulting depth profiles will depend strongly on the coarse-graining width w, which needs to be carefully selected. Representative depth profiles for particular heights, inclinations and basal roughnesses will also be analyzed.
Depth profiles for steady uniform flow are averaged using a coarse graining width w over x ∈ (0, 20], y ∈ (0, 10] and t ∈  with χ w in turn the macroscopic field(s) of density, momentum and stress, as defined in (13), (16) and (32). We average in time with time snapshots taken every t c /2 units.
To determine an appropriate time averaging interval T , we calculate the rate of change in momentum from the density, velocity and stress fields by For steady flow, the temporal variations in mass and momentum should approach zero when averaged over a long enough time interval T . This is shown in Fig. 5, where we plot the depth-averaged norm of the momentum rate of change for varying time averaging interval. For T ≥ 100, the temporal fluctuations decrease to less than 2% of the largest term, ρg, in the momentum equation. In the remainder, we choose T = 100 as the averaging interval. The effect of varying coarse-graining width w is shown in Fig. 6, which shows the z-profile of particle volume fraction ρ/ρ p , where ρ p is the particle density. For small w we observe strong oscillations of about 0.9 particle diameters width, particularly at the base. The microscopic oscillations are increasingly smoothed out and finally vanish as we approach w = 0.5. For larger w, such as w ≥ 1, the macroscopic gradients at the base and surface are smoothed out, an unwanted effect of the coarse-graining. The same behaviour is observed in the stress and velocity fields. Smoothing over the microscopic effects makes it impossible to observe microscopic layering in the density, which we still wish to identify in our averaged fields. Hence, we choose w = 0.25 as the coarse-graining width, such that layering effects remain visible along with macroscopic gradients.
The microscopic oscillations at the base indicate a strong layering effect of particles near the boundary, despite the rough bottom surface. The macroscopic density throughout the flow is almost constant in the bulk and decreases slightly at the base. An approximately constant density profile is a feature of all steady flows and is an important assumption of depth-averaging. Non-zero stress components are plotted in Fig. 7. We observe that the stress components are nearly symmetric. Shear stresses σ yx and σ yz are zero since the flow velocity in y-direction vanishes. For steady flow, the downward normal stress σ zz (z) is lithostatic and satisfies equation (41) with a maximum error of 0.5%. Since the density is nearly constant, we obtain a linear stress profile, another assumption of depth-averaged theory. Applying the momentum balance (19) to steady uniform flow further yields that the shear @ @ @ @ @ I λ Fig. 8: Flow velocity profile of thick flow for H = 30, θ = 24 • , λ = 0, 2/3, 5/6, 1, 2 and for H = 30, θ = 22 • , 26 • , λ = 1/2. For a rough base with λ ≥ 5/6, we see a Bagnold velocity profile (dashed line), except near the surface. For a smooth bases with λ ≤ 2/3, the profile becomes more convex. For λ = 1/2, θ < 24 • , the flow velocity shows layering while still observing the Bagnold profile. For λ = 0, a considerable slip velocity is observed. For λ = 2, the basal shear is small due to flow particles trapped between basal particles so that the definition of the base b(x) is rather fuzzy. stress satisfies σ xz = ∞ z ρ(z ′ )g sin θ dz ′ . Thus, the macroscale friction µ satisfies µ = σ xz /σ zz = −g x /g z = tan θ . This relation is locally satisfied for all steady flow cases to an accuracy of |θ − tan −1 (µ)| < 0.4 • . The remaining normal stress components, σ xx and σ yy , are not constrained by this mass balance. We thus see in Fig. 7 significant anisotropy in the amplitude of the normal stresses, in particular in σ yy , showing that the confining stress is largest in the flow direction, except for very small inclinations. It is always weakest in the lateral or y-direction with fluctuations at the base that are in phase with the fluctuations of the density. Generally, the anisotropy increases with higher inclinations and smoother bases, as will be analyzed in future work.

Transition from rough to smooth base
Next, we study the effect of smoother bases on the range of steady flows by decreasing the diameter λ of the base particles, with the limiting case of a flat bottom wall for λ = 0. Such an extensive numerical study of the effects of changing bottom roughness appears to be novel. To that effect, the DPM simulations from Section 6 were extended such that results for basal roughnesses λ = 0, 1/2, 2/3, 5/6, 1, and 2 can be compared. Before we show the h stop -curves for these simulations, we investigate the extent to which changes in basal roughness lead to more complex density and velocity profiles.
We summarize the density profiles seen without explicitly showing the results. For decreasing basal roughness λ , we observe that the microscopic oscillations and the dip in density at the base increase, while the bulk density remains constant. For λ = 1/2 and λ = 2/3 and small inclinations, we see steady flow that is strongly layered throughout the flow. In contrast, for λ = 2 there is a low flow density in the basal region, since some of the free particles are small enough to sink a little into the base, forming a mixed layer of fixed and free particles.
Velocity profiles for H = 30 and θ = 24 • are shown for varying basal roughness in Fig. 8. For λ = 1, we observe the Bagnold profile [Bag54] for thick collisional flows, differing only at the surface. For very thin flows at H = 10 or inclinations near the arresting flow regime, the profile differs strongly from the Bagnold profile and becomes linear. For smoother bases, the flow velocity increases, and the profile becomes more concave. Weak to stronger slip velocities are observed for λ < 2/3. For λ = 0, thicker flows have constant velocity throughout the depth, almost corresponding to plug flow.
A family of demarcation curves h stop (θ ; λ ) between steady and arrested flow is shown in Fig. 9. The curve fits are based on in which the dependencies on λ are explicitly denoted. The fitting parameters δ 1,λ , δ 2,λ , A λ appearing in (45) are found in Table 1. Again, a fit based on the original equation (35) (or (45)) rather than Pouliquen's early fit (38) yields the best results. For a flat bottom, such that λ = 0, steady flow initiates and resides at or very tightly around an inclination θ = 12.5 • for all heights, see Fig. 9. It is in agreement with the angle found in the laboratory experiments of [GTDD03]. Hence, for a smooth base the flow is steady only at a single inclination, arrests for lower inclinations and accelerates for larger inclinations. Such behaviour is in line with laboratory observations and DPM simulations, e.g. [VATK + 07]. For 1/2 < λ ≤ 2, we observe Pouliquen-style behaviour in Fig. 9. The angle δ 1,λ of flow initiation is nearly constant with respect to λ . In contrast, the range of angles at which both steady and arrested flow is possible, δ 2,λ − δ 1,λ , is maximal for λ = 1 and decreases for smoother chutes with λ < 1, as follows from Table 1. This has been reported in [GTDD03] for laboratory experiments, who also observed a slight decrease of the interval δ 2,λ − δ 1,λ for λ > λ c ≈ 2. However, their λ c was measured for basal particles fixed at the same height and depended on the compactness of the base. We observe a slight decrease of δ 2,λ for λ = 2; however, the fitting curves in Fig. 9 do mildly overlap for λ ≥ 1.
We recall that δ 1,λ and δ 2,λ are fitting parameters for the h stop -curve (45) which does not necessarily imply, though it is expected, that the flow accelerates for angles greater than δ 2,λ . Surprisingly, while steady flow is observed exclusively for θ ∈ (δ 1,1 , δ 2,1 ) when λ = 1, the range of angles associated with steady flow for smoother chutes (i.e., when λ < 1) extends to greater inclinations with θ > δ 2,λ . For these latter cases, δ acc,λ > δ 2,λ is defined as the smallest angle at which accelerating flow is observed; the DPM simulations show that  Note that for the λ = 0 case, between angles 12.5 • and 25 • , the flow is steady and layered, because the friction factor is nonzero.
The above is illustrated in Fig. 10 for λ = 1/2, where one observes the two different steady state regimes. At higher angles, δ 3,1/2 < θ < δ acc,1/2 , a disordered regime similar to that for a rough base is observed. At smaller angles, δ 1,1/2 < θ < δ 3,1/2 , the flowing system self-organizes into a state of layered flow consisting of ordering in the x-y-plane for the bulk (bottom left panel of Fig. 10), except for a small intermediate region, θ ≈ δ 3,1/2 , where a transitional flow regime can be found. It is characterized by large oscillations in the ratio of bulk averaged kinetic to elastic energy due to a spontaneous ordering and disordering, or stop-and-go flow, of the system as a function of time (lower right panel). The same flow regimes have been observed in [SGPL02], where the smoother bottoms were achieved by arranging the base particles in a grid-like fashion. In contrast, we always use of a fully disorded base and vary the roughness by changing the basal particle size.

Closure relations for the depth-averaged model
The goal of this section is to close the shallow-layer equations (34) by a determination of the basal friction µ, the mean densityρ, the stress ratio K, and the velocity profile α, using our DPMs. A demarcation will be made of the flow regimes in which such a determination is possible.

Friction µ of shallow-layer model
For the rough base several friction laws have been proposed, as detailed in Section 4.2. In the following, we will compare these friction laws for the rough base, λ = 1, as well as for varying basal ratios λ .
To obtain a function for the basal friction µ, we used the approach of Pouliquen, who found that for a rough base the Froude number is a linear function of h/h stop (θ ). A first approach was to fit the Froude number to a linear function of h/h stop (θ ; λ ) across the range of non-accelerating DPMs. While this does work for λ ≥ 5/6, a (linear or other) fit does not work for well for λ ≤ 2/3 because for the smoother bases layered and oscllating flows occur for δ 1,λ < θ stop (h, λ ) < θ < δ 3,λ . This is illustrated for λ = 1/2 in Fig. 10. Instead, the Froude number is fitted with h stop (θ ; 1) such that where δ 3,λ is the largest angle at which oscillating flow is observed. In Fig. 10, these angles δ acc,λ and δ 3,λ are shown for the case λ = 1/2. Overall, the simulations reveal that .5 • ± 0.5 • if λ = 1/2 and H = 10, 24.5 • ± 0.5 • if λ = 1/2 and H > 10, 24.5 • ± 0.5 • if λ = 2/3 and H = 10, θ stop (h; λ ) otherwise.
The results of such fits to the Pouliquen law for δ 3,λ < θ < δ acc,λ are shown in Fig. 11, with corresponding fitting parameters provided in Table 1. Shown is the Froude number F =ū/ √ g cos θ h against the ratio of flow and stopping heights h/h stop (θ ; 1). For the disordered steady flow regime, concerning angles δ 3,λ < θ < δ acc,λ , the data are seen to fit better with the stopping angle h stop (θ ; λ = 1), the one for basal surface λ = 1, rather than with the actual stopping height h stop (θ ; λ ). This is a key observation. It shows that the Froude number F increases as the roughness λ decreases, due to the lower dissipation at the base. The weaker Froude number dependence for λ = 2 seen in the right panel of Fig. 11 is in line with the zero shear observed at the base in Fig. 8. The full set of fitting parameters and the standard error for the fit to (47) are found in Table 1 with a standard error defined by We remark that a fit to equation (36) is marginally better than Jenkins' adaption (40), but the differences are too small to discriminate accurately. The situation for layered and oscillating flows is more complicated. We illustrate that for the case λ = 1/2. Two fits are shown in Fig. 11, one for the layered case (dotted line concerning the crosses), and one for the steady case (solid line concerning the circles). The oscillating flows seem to defy a sensible fit because the flow swings irregularly between the layered and disordered states. That oscillating behavior was also shown in Fig. 10 (bottom right panel).
9.2 Functionsρ, K, α of shallow-layer model DPM simulations of steady uniform flows are considered for disordered steady flow with δ λ 3 < θ ≤ δ λ acc , to determine closures forρ, K and α as functions of continuum fields u and h. The layered and oscillating flow regimes are thus excluded momentarily.
All steady disordered flows show a constant density profile in the bulk of the flow, cf. Fig. 6, while the density decreases near the base and the surface. The lower density region at the base spans about two particle diameters for λ > 0, and less than 9d for λ = 0, while the surface region spans always less than 4d. Thus, a mean bulk density can roughly be defined as In Fig. 12, the bulk volume fraction and the mean volume fraction are shown for roughness λ = 1 and varying height and inclination. The bulk volume fraction decreases with inclination θ , but is independent of flow height and roughness, whereas the mean volume fraction depends also on flow height and roughness. We fit the mean bulk density of all steady disordered flows with λ > 0 to an arbitrary function Secondly, the normal stress ratio K =σ xx /σ zz is determined. It describes the anisotropy of the stress tensor and is expected to be unity under isotropic conditions. The range of K for steady disordered flow is generally small, ranging from 0.98 to 1.07, except for λ = 0, where it can be as low as 0.68. The stress anisotropy generally increases with inclination. For λ > 0, K fits to a function linear in θ , with d 0 = 132 • and d 1 = 21.50 • . The model results give a small standard error of err(K − K f it ) = 0.013. Given that the dependence on inclination is small, we may as well take K ≈ 1. Finally, we develop a fit for the shape factor α(λ ) = u 2 /u 2 . The fit is based on a phenomenological model of the observed velocity profiles, as shown in Fig. 8. For rough bases λ ≥ 5/6, a Bagnold velocity profile, is observed in the bulk of the flow; a linear profile in the surface layer, which is about 5d thick; and, a convex profile with no slip in the base layer, whose thickness b λ increases as the height approaches the stopping height. No kinks occur at the intersection of the layers. Thus, we model the velocity  (54). Right: Shape factor α for λ = 1 and varying height h and inclination θ . Markers denote the simulation data, while dotted lines denote fits using (54) with corresponding coefficients from Table 2. Fitted values and simulation data are connected by a solid line.  Fig. 8. For rough bases, the strain is modelled by a Bagnold profile, except near the base and surface. For smoother bases, λ ≤ 2/3, the profile becomes linear. For λ = 0, a large slip velocity is observed, and the strain becomes inversely proportional to the depth.
A fit to the strain ∂ z u in Fig. 14 shows it is well approximated by this model. The parameter b λ decreases with increasing distance from the stopping height, and a simple fit reads where h stop (θ , 1) was chosen since h stop (θ , λ ) does not provide values for all inclinations for which steady flow is observed. Subsequently, the shape factor α(λ ) =ū 2 /u 2 can be computed numerically and compared to the measured values in Fig. 13. The coefficients b λ are given in Table 2. For λ ≤ 2/3, the dependence of the shape factor on height and inclination diminishes and can be approximated with a constant value α(λ ). The Bagnold profile disappears and the flow becomes more convex and plug-like, as shown in Fig. 8. Each velocity profile will be analyzed in turn next.
For λ = 0, the slip is so large that we can assume plug flow to hold. There is almost no slip for λ = 2/3 and a slip of approximately u(0)/ max z u(z) = 1/6 for λ = 1/2. We neglect the variations at the surface and the bulk and model. Thus, we model the velocity profiles as The corresponding coefficients α(λ ) are found in Table 2 and provide a good fit to the data. In summary, the functionsρ, α and K depend on the inclination θ and the height h. The inclination θ in turn can be written as a function of the friction coefficient µ such that θ = tan −1 (µ(h,ū)). This allows us to describe the parameters of the shallow-layer model in terms of the height h, roughness λ , and friction µ(h,ū) and thus provides a closure proper for the system. The different behavior for the varying λ 's remains an open issue, since we only provided emperical fits above.  Table 2: Fitting for the shape factor α = α(λ ) for λ ≤ 2/3 and α =ū 2 /u 2 , u = u(z; b λ ) for λ ≥ 5/6, and the standard error. Closure relations are fitted to all data sets of steady unordered flow, δ λ 3 < θ ≤ δ λ acc .

Summary
In this article, an extensive series of DPM simulations was used to determine closure relations for the popular shallowlayer model of granular flows on inclined chutes. The latter model is a depth-averaged continuum model with a macroscale variable thickness h = h(x,t) and with a mean velocitȳ u =ū(x,t) as variables. For simplicity, we assumed uniformity in the lateral y-direction. The flow consisted of monodispersed particles of diameter d and the base of monodispersed particles of diameter λ d. The bottom roughness, or diameter ratio, λ was systematically varied. Simulations revealed the existence of a range of chute inclinations θ for which horizontal and temporal variations are small enough to produce an approximately steady and uniform flow. Particle flows with variations in height h and inclination θ were numerically investigated for varying basal roughness λ . We observed the following phenomenology: at small inclinations, the flow quickly forms a static pile, while at large inclinations, the flow continues to accelerate. Between these two regimes there was a range of inclinations in which steady flows occurred (cf., Fig. 4). The curve h stop (θ ; λ ), a function of height versus inclination, forms the separatrix between arrested and steady flow, as a function of basal roughness (Figs. 4 and 9). For smaller basal roughness, steady states arise at smaller inclinations and heights, and the range of angles shrinks for which steady flow is possible. Other types of steady flow were observed at small inclinations for small base particles, showing a strong layering in depth as well as disordered and oscillatory flows (cf., Fig. 10).
Depth profiles for density, velocity and stress were constructed using coarse-grained macroscopic fields. The coarsegraining width was carefully chosen to preserve some microscopic structure as well as macroscopic gradients (cf., Fig. 6). The assumptions of depth-averaged theory were confirmed in the simulations for a certain range of steady, uniform flows: the density was constant at depth, and the downward normal stress as well as shear stress were lithostatic. A often-used key assumption, often used, is that statistically steady DPM simulations, or laboratory experiments, are relevant to find the closure relations even for time-dependent continuum shallow-layer models.
Consequently, four closure relations could in principle be determined: for basal friction µ, stress ratio K, mean densityρ, and for the shape of the velocity profile α. Firstly, basal friction µ = tan θ was shown to be a function of height and flow velocity (Table 1). Pouliquen's approach was found to be valid, with the Froude number as a linear function of h/h stop (θ ; λ = 1). This fitting approach was extended to smoother cases with λ < 1, where the Froude number was fitted to h/h stop (θ ; λ = 1) instead of the h/h stop (θ ; λ ) for the actual λ or basal roughness. The stopping curve associated with the diameter λ = 1 of the flowing particles is more relevant than the stopping with the actual λ . One possible explanation is that there is a boundary layer of intermittently slow flow particles that originated in the bulk, and that shields the smoother base from the bulk flow. Closure relations for the mean densityρ, stress anisotropy K and shape factor α were also established as follows. For rough bases with λ ≥ 5/6, the determined closures were valid for θ stop (h; λ ) = δ 3,λ < θ < δ acc,λ , with θ stop (h; λ ) the inverse of the h stop (θ ; λ )-curve between arrested and dynamic flow. For smaller roughnesses with λ ≤ 2/3 and θ stop (h; λ ) < θ < δ 3,λ , layered and oscillating flows arose for which we are (as yet) unable to capture closures. For these smoother bases with λ ≤ 5/6, closures were obtained for the range δ 3,λ < δ < δ acc,λ .

Open questions
What does the granular shallow-layer model enable us to do, and what can we not do with it? In the range of steady flows, this continuum model can be used to predict steady and time-dependent flows. Strictly speaking, this is only allowed for steady flow in the established inclination range δ 3,λ < θ < δ acc,λ , but it can be expected to remain valid for the slowly-varying dynamic cases as well. It is often the case, however, that even rapidly-changing flows can be captured by models that should only be valid for the slowlyvarying cases. Consequently, a systematic study of the validity of the shallow-layer model is required. By respectively extending the "hydraulic" analysis for fluidized granular matter and water in Vreman et al. [VATK + 07] and Akers and Bokhove [AB08], granular flows within constrictions become a nice, analytically-treatable targets. Such flows in constrictions reach a steady state and appear (partially) accessible by direct DPM simulations.
What do these results enable us to do? Whether the steady DPM-based closures are valid across granular "hydraulic" jumps in such steady and constricting flows is of interest. Whether the steady DPM-based closures hold for (slow) transient routes towards such steady states is of interest, too.
What closures should be used outside the formal range of applicability for the smoother bases, so for the layered and oscillating flows for θ stop (h; λ ) < θ < δ 3,λ and the accelerating flows for θ > δ acc,λ , appears a tantalising, and as yet, open question.
What are we not able do? Although, we did observe layered and oscillating flows in our DPM simulations, it is doubtful as to whether the homogenization assumption that led to the shallow-layer model is sufficient. Nonetheless, the lithostatic balance relation is shown to hold for the DPM simulations, as expected from standard asymptotic analysis using the aspect ratio of normal to planar velocity and length scales.

Outlook
Alternatively, a multi-scale modelling approach might be adopted such as the heterogeneuous, multiscale methodology [WEL + 07], among others, in which closure relations for discretisations (e.g., [?]) depth-averaged shallow-layer models are coupled to DPM simulations in selected regions in space and time. Thus computational costs would be diminished while accurate closure relations are gathered intermittently in time and space.
For future work, we advocate the extension of our DPM simulations with investigation of the three-dimensional closure relations. We surmise that reduced lithostatic models for shallow granular flows could be more consistently derived from three-dimensional continuum models with stress closure determined from DPM simulations in combination with laboratory measurements. These new models would be reduced and therefore computationally still manageable for large-scale debris flows; for example, when the degrees of freedom in the vertical remain limited, but are extended beyond only one degree of freedom. Such reduced modelling is akin to hydrostatic modelling in water-wave and coastal hydrodynamics.

A System of Hamiltonian equations in the dissipation-free limit
The purpose of this appendix is to show that the particle system in the dissipation-and yield-free limit, i.e., γ n = γ t = 0 and µ c → ∞, is a Hamiltonian system. It subsequently facilitates the derivation of so-called conservative or symplectic discretization schemes, in time. These have been shown and are believed to provide better long-term statistics than classical time discretization schemes. Furthermore, analysis of the Hamiltonian limit allows one to clearly demarcate the transfer of energy between its kinetic, elastic, and internal components.
We show that the normal and tangential forces are elastic, that is the system does not dissipate energy; instead kinetic energy is converted into potential energy in the springs and vice versa. If the tangential spring is not fully unloaded when two particles loose contact, the potential energy stored in the tangential spring is converted into internal potential energy in each particle (vibrations).
We use the notation given in §2. In the dissipation-and yield-free limit, the contact force between particles i, j is given by where δ n i j and δ δ δ t i j are given by equations (1) and (9). The equations of motion for translational and angular momentum of particle i are given by, in three dimensiones, where r i is the position, α i is the angle, p i the momentum and φ i the angular momentum of particle i. To define the Hamiltonian system, we pair these generalized position and momentum vectors as follows Then the kinetic energy can be calculated using only the generalized momenta P as follows The potential energy is a combination of the potential of gravity, the potential of the normal and tangential springs and the internal potential energy in the particles, created from the remaining potential of the tangential spring at the time t e i j that a particle pair {i, j} looses contact, {t e i j } = {t : with r i j = r j (t) − r i (t). The potential can be expressed in terms of the position and the tangential springs at all times, which itself is a function of the previous positions of the particle pair, where the gravitational, elastic and internally stored potential energy is defined by Substituting (69) into (68) where we used the identity Thus, using (68) and that b i j and n i j are parallel, we obtain ∂ H Subsequently, we derive (62b) since and Finally, we show that the total energy is conserved. Since mass m i , radius r i and spring constants k n , k t are constant, H has no direct dependence on t and thus Using and (62) and (75) yields Fig. 15 shows the energy balance of two particles colliding noncollinear in the dissipation-and yield-free case. One can see the jump in energy at the end of contact, where potential tangential spring energy is converted into internal energy.  Fig. 15: Energy balance of two particles colliding noncollinear in the dissipation-and yield-free case. The kinetic and elastic potential energy are in balance, until the particles loose contact and potential spring energy is converted into internal energy.

B Orthogonality of the tangential spring
To show that the tangential spring is orthogonal to r i j , note that the tangential spring is initially of zero length and therefore orthogonal to r i j ; further Thus, we can integrate equation (77) to obtain a continuously orthogonal tangential spring with δ t i j · r i j = 0.

C Algorithms for time integration and the calculation of the tangential force
The algorithm for the time integration and the calculation of the tangential force is shown in Algorithms 1 and 2.