Rarefied particle motions on hillslopes: 1. Theory

We describe the probabilistic physics of rarefied particle motions and deposition on rough hillslope surfaces. The particle energy balance involves gravitational heating with conversion of potential to kinetic energy, frictional cooling associated with particle-surface collisions, and an apparent heating associated with preferential deposition of low energy particles. Deposition probabilistically occurs with frictional cooling in relation to the distribution of particle energy states whose spatial evolution is described by a Fokker-Planck equation. The Kirkby number Ki — defined as the ratio of gravitational heating 5 to frictional cooling — sets the basic deposition behavior and the form of the probability distribution fr(r) of particle travel distances r, a generalized Pareto distribution. The shape and scale parameters of the distribution are well-defined mechanically. For isothermal conditions where frictional cooling matches gravitational heating plus the apparent heating due to deposition, the distribution fr(r) is exponential. With non-isothermal conditions and small Ki this distribution is bounded and represents rapid thermal collapse. With increasing Ki the distribution fr(r) becomes heavy-tailed and represents net particle heating. 10 It may possess a finite mean and finite variance, or the mean and variance may be undefined with sufficiently large Ki . The formulation provides key elements of the entrainment forms of the particle flux and the Exner equation, and it clarifies the mechanisms of particle-size sorting on large talus and scree slopes. Namely, with conversion of translational to rotational kinetic energy, large spinning particles are less likely to be stopped by collisional friction than are small or angular particles for the same surface roughness. 15


Introduction
Sediment transport on steepland hillslopes involves a great range of scales of particle motions. These vary from relatively small motions that collectively produce the slow en masse motion of disturbance driven creep (Culling, 1963;Roering et al., 1999Roering et al., , 2002Gabet, 2000;Anderson, 2002;Gabet et al., 2003;Furbish, 2003;Roering, 2004;Furbish et al., 2009bFurbish et al., , 2018a in concert with athermal granular creep (Houssais and Jerolmack, 2017;BenDror and Goren, 2018;Ferdowsi et al., 2018;Deshpande 20 et al., 2020) to the long-distance and relatively fast en masse motions of landsliding and the rarefied motions associated with rockfall and ravel (Kirkby and Statham, 1975;Statham, 1976;Dorren, 2003;Gabet, 2003;Roering and Gerber, 2005;Luckman, 2013;Tesson et al., 2020). Particularly in relation to long-distance motions, there is a growing interest in non-continuum Figure 1. Image of talus slope at the base of cliffs of the Bandelier Tuff showing downslope sorting of particle sizes, with the largest particles preferentially accumulating near the base of the slope. The largest boulders in the foreground are about 1 m in diameter. As described in the text, we suspect that with conversion of translational to rotational kinetic energy, large spinning particles are less likely to be stopped by collisional friction than are small or angular particles for the same surface roughness, thus contributing to the sorting in this image. Image location is at the confluence of the Rito de los Frijoles river canyon with the Rio Grande River canyon on the eastern boundary of the Bandelier National Monument, New Mexico, USA.
to a heavy-tailed form representing net particle heating. In Section 5 we compare the formulation with previous mechanical descriptions of disentrainment, showing both similarities and dissimilarities with these descriptions.
We emphasize that this initial phase of our work on rarefied particle motions is aimed at clarifying how particle disentrainment works. With this in place we will be positioned to consider effects of rarefied transport over time scales spanning many transport events, including ensemble-averaged particle fluxes and changes in land-surface elevation as described by 5 formulations of nonlocal transport. As a step in this effort we show in the second companion paper (Furbish et al., 2020a) that the theory in this first paper is entirely consistent with data from laboratory and field-based experiments involving measurements of particle travel distances on rough surfaces. These include data reported by Kirkby and Statham (1975), Gabet and Mendoza (2012), DiBiase et al. (2017) and Roth et al. (2020), and new travel distance data from laboratory experiments supplemented with high-speed imaging and audio recordings that highlight effects of particle-surface collisions. Outstanding questions concern how particle size and shape in concert with surface roughness influence the extraction of particle energy and the likelihood of deposition. In the third companion paper (Furbish et al., 2020b) we show that the generalized Pareto distribution in this problem is a maximum entropy distribution (Jaynes, 1957a(Jaynes, , 1957b constrained by a fixed energetic "cost" -the total cumulative energy extracted by collisional friction per unit kinetic energy available during particle motions. That 5 is, among all possible accessible microstates -the many different ways to arrange a great number of particles into distance states where each arrangement satisfies the same fixed total energetic cost -the generalized Pareto distribution represents the most probable arrangement. In the fourth companion paper (Furbish et al., 2020c) we step back and examine the philosophical underpinning of the statistical mechanics framework for describing sediment particle motions and transport.

Continuous form
Following the presentations of Furbish and Haff (2010) and Furbish and Roering (2013), let f r (r; x) denote the probability density function of particle travel distances r whose motions begin at position x. By definition the cumulative distribution function is 15 where the prime denotes a variable of integration. In turn, the exceedance probability, also referred to as the survival function, is R r (r; x) = 1 − F r (r; x) = ∞ r f r (r ; x) dr .
With these definitions in place we now define the spatial disentrainment rate as 20 which is a conditional probability per unit distance. Namely, upon multiplying both sides of Eq. (3) by dr, then P r (r; x)dr = f r (r; x)dr/R r (r; x) is interpreted as the probability that a particle will become disentrained within the small interval r to r + dr, given that it "survived" travel to the distance r. The disentrainment rate P r (r; x) also may be interpreted as an inhomogeneous Poisson rate (Feller, 1949). Now, using the fact that f r (r; x) = −dR r (r; x)/dr, one may deduce from Eq. (3) that the probability density f r (r; x) is given by 25 f r (r; x) = P r (r; x)e − r 0 Pr(r ;x) dr .
Thus, according to Eq. (4), the disentrainment rate P r (r; x) completely determines the probability density f r (r; x) of travel distances r.
Assuming particle motions occur only in the positive x direction, the entrainment form of the volumetric particle flux is where E s (x) denotes the volumetric entrainment rate at position x. In turn, letting ζ(x, t) denote the local land-surface elevation, the entrainment form of the Exner equation is (Tsujimoto, 1978;Nakagawa and Tsujiomoto, 1980) 5 where c b = 1 − φ s is the particle volumetric concentration of the surface with porosity φ s . These probabilistic formulations of the flux and the Exner equation have three lovely properties. They are mass conserving, they are nonlocal in form, and they are scale independent. They illustrate that the probability density f r (r; x) of particle travel distances r and its related survival function R r (r; x) form the centerpiece of describing mass conservation and the particle flux. In turn, the significance of the disentrainment rate P r (r; x) becomes clear. This rate connects Eq. (5) and Eq. (6) to the physics of particle motions on a hillslope. That is, this rate, together with the entrainment rate E s (x), represent the elements in the formulation that can be elucidated by physics.
To date, previous formulations of the disentrainment rate P r (r; x) have envisioned a friction dominated behavior in which the land-surface slope S(x) = |∂ζ(x)/∂x| has a primary role (Furbish and Haff, 2010;Furbish and Roering, 2013;Doane, 2018;Doane et al., 2018a;Section 5.3). The disentrainment rate is specified as a function of the land-surface slope at the position 15 of entrainment, with the idea that the slope changes over a distance much larger than the average particle travel distance.
That is, P r (r; x) is assumed to be a determined by the slope S(x) at position x such that the distribution of travel distances of particles entrained at x is exponential with mean µ r [S(x)]. As the land-surface slope S varies with increasing downslope distance x, the mean µ r [S(x)] changes. The disentrainment rate is qualitatively consistent with limiting cases, namely, it yields a fixed small average travel distance at zero slope, and it approaches zero in the limit of a steep critical slope beyond which 20 disentrainment does not occur. However, the mechanical elements of the disentrainment rate P r (r; ) are otherwise not explicitly specified. We also note that Kirkby and Statham (1975) first pointed out the relation between the distribution of travel distances and the disentrainment rate function. These authors defined a posteriori the disentrainment rate from an assumed exponential distribution of travel distances whose mean value is expressed in terms of a Coulomb-like description of particle friction (Section 5.1).

Discrete form
It is valuable to recast the ideas of disentrainment above in discrete form. The motivation is this. Instead of trying to formulate a continuous disentrainment rate function that is generally applicable to the entirety of a hillslope, we instead break it into discrete spatial intervals, where certain physics may be more or less important in some intervals than in others. This gets us closer to the physical ingredients of disentrainment that are occurring at different locations on a hillslope, where the mechanical 30 behavior at a location transitions to another behavior in the downslope direction. We may then combine the intervals together as a whole.
Let k = 1, 2, 3, ... denote a set of discrete intervals of length dr. Let p denote the probability that a particle is disentrained within the first interval (k = 1). If N denotes a great number of particles, then the number of particles n(1) disentrained within the first interval is n(1) = N p. Because q = 1 − p is the probability that a particle is not disentrained within the first interval, 5 then the number of particles moving beyond the first interval is N q = N (1 − p). That is, this is the number of particles that "survived" without being disentrained within the first interval. In turn, of the number of particles that survived, the number that is disentrained within the second interval is n(2) = N (1 − p)p. More generally, n(k) = N (1 − p) k−1 p. Dividing this by N then gives the probability mass function 10 which defines the well-known geometric distribution with mean µ k = 1/p. Note that the probability p is taken here as being fixed. That is, in this formulation, the probability that a particle survives the kth interval is (1 − p) k−1 , so the disentrainment probability is constant, namely, P k (k) = p.
The geometric distribution, Eq. (7), is the discrete counterpart of the exponential distribution. Here we relate the two. The cumulative distribution function of Eq. (7) is F k (k) = 1 − (1 − p) k . We may thus write F r (r = kdr) = F k (k) = 1 − q k . The 15 quantity q k is a memoryless geometric series, and because q ≤ 1 we may write q = e −dr/µr , where µ r is a characteristic distance. In turn, then, F r (r) = 1 − (e −dr/µr ) k = 1 − e −r/µr . Finally, f r (r) = dF r (r)/dr = (1/µ r )e −r/µr , where it becomes clear that µ r is the mean of the exponential distribution, analogous to µ k for the discrete counterpart. Also note that the disentrainment rate P r (r) = 1/µ r is fixed. Below we show that the exponential and geometrical distributions represent isothermal conditions, where gravitational heating of particles is balanced by frictional cooling. 20 In contrast, suppose that the probability of disentrainment p varies from one interval k to another. Here we generalize the ideas above. Let p k denote the probability that a particle, having not been disentrained before the kth interval, then becomes disentrained within this interval. Similar to the formulation above, the number of particles n(1) within the first interval is n(1) = N p 1 and the number moving beyond the first interval is N (1 − p 1 ). In turn the number of particles disentrained within the second interval is n(2) = N (1 − p 1 )p 2 , the number disentrained within the third interval is N (1 − p 1 )(1 − p 2 )p 3 , and so 25 on. In general, n(k) = N (1 − p 1 )(1 − p 2 )(1 − p 3 )...(1 − p k−1 )p k . Dividing this expression by N then gives Note that if p k = p is fixed, then Eq. (8) reduces to Eq. (7).
This generalization has a lovely property. Namely, by definition it conserves probability, and it therefore is mass conserving.
That is, the sum of f k (k) over all possible k is equal to unity, regardless of how p k might vary with k. As alluded to above, 30 the physics of each p k may be treated differently if desired. Moreover, like its continuous counterpart presented above, this discrete formulation of mass conservation is nonlocal and scale independent.
We now set these results aside. Our next objective is to illustrate the mechanical elements of disentrainment, which we then use to elaborate the continuous and discrete cases described above.
3 Mechanical interpretation of disentrainment

Conservation of mass
Consider a rough, inclined surface with uniform slope angle θ (Figure 2). At this juncture we simplify the notation and consider Figure 2. Definition diagram of surface inclined at angle θ and control volume with edge length dx through which particles move.

5
the motions of particles entrained at a single position x = 0. Now the particle travel distance r → x and the probability density . Consider a control volume with edge length dx parallel to the mean particle motion. Over a period of time a great number of particles enters the left face of the control volume. Some of these particles move entirely through the volume, exiting its right face, and some come to rest within the control volume. Many, but not necessarily all, of the particles interact with the surface one or more times in moving through the volume or in being disentrained within it. 10 We now imagine collecting this great number of particles and treat them as a cohort, independent of time (Appendix B).
That is, let N (x) denote the number of particles that enter the control volume, and let N (x+dx) denote the number that leaves the volume. We may imagine for the purpose of visualizing the problem that the N (x) particles enter the control volume at the same time, but this actually is not essential. Similarly, we may imagine that the N (x + dx) particles exit the control volume at approximately the same time, but again, this reference to time only is a means to envision particle motions (Appendix B). The 15 number of particles disentrained within the control volume then is If N (0) denotes the great number of particles whose motions started at position x = 0, then the exceedance probability Our objective is to determine the derivative dN (x)/dx in relation to particle energy, as this derivative represents disentrainment. Here we summarize the essence of this problem before turning to a description of conservation of energy. Let E p = (m/2)u 2 denote the translational kinetic energy of a particle with mass m and downslope velocity u. Here we are assuming that the total translational energy is dominated by downslope motion. Let f Ep (E p , x) denote the probability density function of particle energies E p as these vary with position x. For a great number N of particles the number density denote the probability that a particle at energy state E p will become disentrained within the small interval x to x + dx. Because N f Ep (E p , x)dE p is the number of particles within the small interval E p to dE p is the number of particles in this energy interval that becomes disentrained. The total number of particles that becomes disentrained within the interval x to x + dx is then Letting angle brackets denote an ensemble average, then according to the Law of the Unconscious Statistician, Eq. (10) is simply dN (x) = −N (x) p(E p , x) . Below we introduce the expected number of particle-surface collisions per unit dis-10 tance n x = 1/λ, where λ is the expected travel distance between successive collisions. We then show that dN (x)/dx = −N (x)n x p(E p , x) . Thus, the essence of the problem is to determine the averaged probability p(E p , x) as this depends on particle energy E p . This in turn requires specifying the particle energy as this varies with position x.

Particle energy
We start our formulation with a general statement concerning conservation of the kinetic energy of a system of particles. 15 Because of its familiarity in relation to studies of granular gas systems, we initially consider changes with respect to time, then return to changes with respect to space as in the preceding section. Namely, let E p denote the kinetic energy of a particle, and let E p denote the expected energy state, where angle brackets represent an ensemble average over a great number N of moving particles. The total energy of the system is E = N E p . Neglecting transport of energy, the rate of change in the total energy of the system with respect to time is then The first term on the right side of Eq. (11) represents the rate of change in the average energy state of N moving particles, and thus describes either a net heating (d E p /dt > 0) or cooling (d E p /dt < 0) of the system, depending on the relative contribution of the sources of each. The second term on the right side represents the rate of change in the number of moving particles with average energy state E p , and thus describes the rate of change in the total energy due to either the addition or 25 loss of moving particles. For a closed system, this represents either a net sublimation (dN/dt > 0) or net deposition (dN/dt < 0) of particles, depending on the relative contribution of each.
The first term on the right side of Eq. (11) has been studied extensively for granular gas systems, specifically in relation to the "homogeneous cooling state" of a closed system as described by Haff's cooling law (Haff, 1983;Brilliantov and Pöschel, 2004;Dominguez and Zenit, 2007;Volfson et al., 2007;Brilliantov et al., 2018;Yu et al., 2020). In what follows, we start with 30 similar concepts of particle energy; but the formulation is designed to be independent of time and focused on changes in energy and particle disentrainment over space.
Reconsider a control volume with edge length dx parallel to the mean motion of particles over a rough, inclined surface ( Figure 2). Analogous to Eq. (11) we write where now the angle brackets formally denote a Gibbs ensemble average over a cohort of particles (Appendix B). As described below, the first term on the right side of Eq. (12) represents the spatial rate of change in energy due to the sum of gravitational 5 heating and frictional cooling. The second term on the right side represents the rate of change in energy due to deposition, that is, disentrainment. In this problem, we assume that sublimation (entrainment) does not occur over x > 0. Eq. (12) provides a basic starting point. However, it is not particularly useful in this form. If in fact the probability of deposition varies with energy state, then in general the derivative dN/dx contributes to the derivative d E p /dx, as removal of energy by deposition affects the average energy of the remaining particles. We note that Brilliantov et al. (2018) demonstrate an analogous effect, 10 as described below, associated with aggregation of particles in a granular gas. We therefore must be careful in formulating a statement of conservation of particle energy, as deposition preferentially involves particles at low energy states.
Here is a sidebar concerning our focus on conservation of particle energy versus momentum. Particle motions down a rough hillslope surface involve numerous details that control momentum exchanges during particle-surface interactions. As a scalar quantity, energy forces us to blur our eyes appropriately, focusing on the essence of these complex interactions rather 15 than attempting to describe details of momentum exchanges that ultimately cannot be constrained given the stochastic nature of the phenomenon. As an example, below we introduce the random variable β x to represent the proportion of downslope kinetic energy extracted during a particle-surface collision. This quantity blurs over many details (e.g., differences between collisions during rolling, tumbling and bouncing motions, rotational versus translational motion, and the roles of normal and tangential coefficients of restitution), yet β x is entirely meaningful when treated as a random variable. (We provide a description 20 (Appendix E) of how the energy-centric quantity β x is related to momentum exchanges during collisions, and in the companion paper we illustrate the elements of β x using high-speed imaging.) In contrast, when describing the collisional behavior of an ideal granular gas, one can at lowest order appeal to a single coefficient of restitution because of the relative simplicity of the particles and their collisions (e.g., Haff, 1983;Jenkins and Savage, 1983). This simplicity is not possible here. The focus on energy thus offers tractable and defensible simplicity amidst the messiness of natural hillslopes. Focusing just on slope parallel motions, let E p = (m/2)u 2 denote the translational kinetic energy of a particle with mass m and downslope velocity u. Then let f Ep (E p , x) denote the probability density function of particle energies E p as these vary with downslope position x (Appendix B). For a great number N of particles the number density is n Ep (E p , x) = N f Ep (E p , x). 30 9 https://doi.org/10.5194/esurf-2020-98 Preprint. Discussion started: 8 December 2020 c Author(s) 2020. CC BY 4.0 License.
The average particle energy is The total energy E(x) = N E p , so We now take the derivative of Eq. (14) with respect to x using Leibniz's rule to give The derivative within the integral of Eq. (15) satisfies a Fokker-Planck equation (see next section and Appendix C), the solution of which represents the evolution of the distribution n Ep (E p , x) of particle energy states E p with distance x. In particular this 10 derivative has three parts. The first part, denoted below by K h (E p , x), is associated with a change in the density n Ep (E p , x) due to gravitational heating. The second part, K c (E p , x), is associated with a change in this density due to frictional cooling.
The third part, K d (E p , x), is associated with a loss of energy due to deposition (which does not involve the analogue of release of latent heat; but see below). We thus write 15 and then rewrite Eq. (15) as The next task consists of showing the correspondence of K h (E p , x), K c (E p , x) and K d (E p , x) to terms in the Fokker-Planck equation, then describing the physical elements of these terms. This is followed by evaluating each of the integral quantities in Eq. (17). There are a lot of moving parts in this formulation, so bear with us.

Fokker-Planck-like equation
The density n Ep (E p , x) within Eq. (15) and Eq. (16) satisfies a Fokker-Planck equation (Appendix C), which describes the evolution of this density with increasing distance x. Namely,

10
The first term on the right side of Eq. (18) represents advective gravitational heating, where k 1h (E p , x) is a drift speed, the average spatial rate of change in particle energy over the energy domain due to heating. The second term on the right side represents advective frictional cooling, where k 1c (E p , x) is a drift speed, the average spatial rate of change in particle energy due to cooling. The third term represents diffusive frictional cooling, where k 2c (E p , x) is a diffusion coefficient. The last term represents a loss of energy due to deposition, where for now we have retained the notation from above. Explicitly, for 15 K h (E p , x) and K c (E p , x) we now have and In the next section we step through gravitational heating, frictional cooling and deposition, in each case unfolding the mechan- Here is a didactic sidebar if the formulation above seems counterintuitive. Notice that Eq. (18) effectively represents an advection-diffusion equation with two advective terms, a diffusive term and a sink term. Normally we think of an advection-25 diffusion equation as involving space and time, that is, where the rate at which a quantity changes with respect to time at a given position is equal to the sum of an advective term and a diffusive term involving derivatives of the quantity with respect to space. Indeed, imagine replacing E p with x, and x with t, in Eq. (18). The result looks like a familiar advection-diffusion equation with a sink term (albeit involving two advective terms rather than one). The basic idea of Eq. (18) is the same. It just describes the rate of change in n Ep with respect to position x (rather than time t) in relation to advection and diffusion of n Ep occurring over the energy coordinate E p (rather than x). A consideration of the rate of change with respect to position x as in Eq. (18) is perhaps unusual, but the idea of advection and diffusion of a quantity occurring over a domain other than a spatial coordinate (e.g., a velocity coordinate) is common in statistical physics, of which examples pertaining to sediment motions include those presented in Furbish et al. (2012Furbish et al. ( , 2018aFurbish et al. ( , 2018b.

Gravitational heating
We start by noting that the rate at which the potential energy of a particle is converted to kinetic energy per unit distance x is mg sin θ. To be clear, between collisions a particle that is not in contact with the inclined surface beneath it accelerates vertically at a rate of −g, independently of the orientation of the surface. The factor sin θ therefore is a geometrical constraint on the magnitude of the potential energy that is accessible for net heating when viewed with respect to x. This means that so that Eq. (19) becomes We now write the first integral in Eq. (17) as Because ∂(E p n Ep )/∂E p = E p ∂n Ep /∂E p + n Ep , Eq. (23) may be written as 20 Assuming n Ep (∞, x) → 0, the second integral in Eq. (24) vanishes and the first integral in Eq. (17) becomes Note that the form of the density n Ep (E p , x) is immaterial in this formulation.
If for illustration we assume that no cooling or deposition occurs, then dE(x)/dx = N mg sin θ. The solution of this is E(x) = E(0) + N mg sin θx, where E(0) denotes the starting energy at x = 0. That is, the total energy E(x) increases 25 linearly with downslope distance x. Moreover, for reference below, no particle can be heated to an energy greater than E p (0) + mg sin θx, representing a complete conversion of gravitational to kinetic energy without any loss due to particlesurface collisions. This ensures that the density n Ep (E p , x) is bounded with finite mean and variance, a point that becomes useful below.

5
We start by assuming that a change in the downslope energy of a particle associated with a collision is ∆E p = −β x E p , so that β x = −∆E p /E p is the proportion of energy extracted by the collision (Appendix E). By definition β x is a random variable.
(Note that the negative sign above is by convention. As a random variable we are assuming that 0 ≤ β x ≤ 1. The sign associated with β x will be clear from the context in the developments below.) The change ∆E p includes frictional loss, any conversion of translational to rotational energy, and any apparent change when downslope incident motion is reflected to transverse motion 10 during a glancing particle-surface collision. Note that ∆E p generally is a negative quantity. But strictly speaking it could be positive, albeit with low probability, if transverse incident motion is reflected to downslope motion during a collision.
Because E p and β x are random variables, ∆E p is a random variable. As a point of reference, in granular gas theory where the total translational energy is considered rather than just the energy associated with one coordinate direction, the proportion β x = 1 − 2 where is the normal coefficient of restitution (Haff, 1983). Moreover, is treated as a fixed deterministic quantity 15 rather than a random variable. Here, in contrast, collision mechanics theory suggests that the constitution of β x is far more complicated in relation to normal, tangential and rotational impulses during particle-surface collisions (Appendix E).
Let q = E p (x + dx) − E p (x) denote a change in the energy of a particle over the small distance dx. Then as described in Appendix C, the drift speed k 1c (E p , x) = dq/dx ≈ n x β x E p and the diffusion coefficient k 2c = dq 2 /dx ≈ n x β 2 x E 2 p , where the overline denotes an average over particles at the energy state E p (rather than an ensemble average), and n x = 1/λ denotes 20 the expected number of particle-surface collisions per unit distance where λ is the expected travel distance between collisions.
Scaling (Appendix D) shows that where φ is the expected reflection angle of a particle with energy E p following a surface collision. We now assume that 25 and that Now Eq. (20) becomes We now use these results to write the second integral in Eq. (17) as Upon applying the product rule to the derivative ∂(E p β x n Ep )/∂E p , the first integral in Eq. (30) may be written as

10
Assuming that n Ep (∞, x) → 0, the second integral in Eq. (31) vanishes and the first integral becomes N β x , where the angle brackets now represent an ensemble average.
In turn, upon applying the product rule to the derivative ∂[E p ∂(β 2 x E p n Ep )/∂E p ]/∂E p , the second integral in Eq. (30) may be written as Assuming that n Ep (∞, x) → 0 and ∂n Ep /∂E p | Ep→∞ → 0, the integrals in Eq. (32) reduce to (mg cos θ/8 tan φ)β 2 x E p n Ep (0, x) with β 2 x E p = 0 when evaluated at E p = 0. Thus, whereas the diffusive term in Eq. (18) redistributes energy by modifying the density n Ep (E p , x) (see below), it does not contribute to the total energy balance. The second integral in Eq. (17) is thus We return to these results below.

Energy loss with deposition
For illustration, suppose initially (unrealistically) that deposition is independent of the particle energy state E p . This means that the number of particles disentrained within any small energy interval E p to E p + dE p is a fixed proportion k 3 of the particles within this interval. Thus, K d (E p , x) = −k d n Ep (E p , x) and the third integral in Eq. (17) becomes 10 paring this result with Eq. (12), the situation in which deposition is independent of the particle energy state is consistent with isothermal conditions wherein the average energy state is unchanging, that is, d E p /dx = 0.
More generally, deposition is unlikely to be independent of the particle energy state, as particles with small energy are on average more likely to become disentrained than are particles with large energy. Thus, K d (E p , x) likely possesses a more 15 complicated form than in the example above. Whereas early work on granular gases focused on their behavior in the absence of deposition, the phenomenon of thermal collapse, condensation and freezing in a gravitational field now is receiving significant attention (Volfson et al., 2006;Kachuck and Voth, 2013). We can lean on insight from this work, but because energy dissipation in a granular gas is dominated by particle-particle collisions rather than particle-boundary collisions, the rarefied problem considered here is quite different. As with approaches used in the study of condensation and freezing of granular gases, our 20 analysis at this stage is aimed at lowest order behavior.
For any position x, we do not know the ensemble distribution f Ep (E p , x) of particle energy states E p with expected value Ep . Because no particle can be heated to an energy greater than E p (0) + mg sin θx (representing a complete conversion of gravitational to kinetic energy without any loss due to particle-surface collisions), we know only that 0 ≤ E p ≤ E p (0) + mg sin θx. Most energies likely are significantly smaller than the upper limit due to collisions.

25
Collecting results from above, the density n Ep (E p , x) satisfies a Fokker-Planck-like equation, namely, where we are assuming for simplicity that β x and β 2 x are fixed. As a reminder, the first term on right side of Eq. (36) represents 5 gravitational heating, and the second and third terms on the right side represent frictional cooling. The term −K d (E p , x), which describes the loss of energy associated with deposition, is defined below.
Let n Ep0 and E p0 denote suitable characteristic values of the density n Ep and the energy E p , and let X denote a characteristic length scale. We now define the following dimensionless quantities denoted by circumflexes: 10 Upon substituting these quantities into Eq. (36), we may identify three characteristic length scales, namely, The first of these, X h , represents the distance required to heat a particle to the energy state E p0 in the absence of frictional cooling. The second, X cA , represents the distance over which thermal collapse by advective cooling occurs. The third, X cD , represents a distance over which diffusive cooling occurs.
We now define two dimensionless numbers, the Kirkby number 1 , and a cooling Péclet-like number, The Kirkby number Ki is the ratio of gravitational heating to advective cooling. The Péclet-like number Pe c is the ratio of advective cooling to diffusive cooling. Choosing X cA as the characteristic length scale and neglecting the deposition term in Eq. (36), we now rewrite it as Note that with β x 1, then Pe c 1 according to Eq. (42), such that the diffusive term in Eq. (43) becomes insignificant relative to the advective cooling term.
With reference to Figure 3, imagine a great number of particles whose initial energy states at x = 0 are described by the and (c) net cooling. The triangular region represents an idealized "window" of increasing likelihood of deposition with decreasing particle energy Ep. Note that an effect of deposition is to increase the average energy Ep by culling lower energy particles, thereby selecting higher energy particles for continued travel with increasing distance. density n Ep (E p , 0). With just gravitational heating, this distribution is advected to higher energy values at a fixed rate mg sin θ.

10
With just frictional cooling, but in the absence of diffusion, the distribution is advected to lower energy values at a fixed rate mgβ x cos θ/4 tan φ. If gravitational heating is balanced by advective cooling (Ki = 1), the form of the distribution remains fixed with increasing distance x. With diffusive cooling, advective cooling of the density n Ep (E p , x) to lower energy values involves smoothing of this density. When these effects are combined, whether heating is greater than advective cooling (Ki > 1) or vice versa (Ki < 1), no value of E p is larger than E p (0)+mg sin θx, and most values are significantly less than this maximum due to the increasing likelihood of particle-surface interactions (cooling) within increasing x. When the magnitude of the term in Eq. (43) involving Ki is greater than the sum of the magnitudes of the two cooling terms, then net heating occurs. When the magnitudes of the cooling terms are larger than the heating term, then net cooling occurs. For particles reaching relatively small energy states, there is an increasing likelihood of deposition (see below). As a reminder, this description does not pertain 5 to the energy states of a great number of particles during an interval of time. Rather, this description pertains to an ensemble of particles reaching any position x over a long period of time when treated as a cohort. That is, n Ep (E p , x) is the density of particle energies at any x representing the great number of particles that occupied this position while in motion at many previous instants in time.
We now offer a simple hypothesis describing the loss of energy associated with deposition. Recall that X cA is a measure of 10 the distance over which particles with energy E p0 thermally collapse by frictional cooling. We may imagine, for example, a sudden removal of the source of heating such that X cA is a measure of the distance of relaxation to a total loss of energy. For particles with energy E p , this length scale can be expressed more generally as which becomes unbounded only in the limit of θ → π/2. Because thermal collapse involves deposition, we then assume at 15 lowest order that where the subscript d denotes that the derivative refers to a change in the density n Ep (E p , x) just associated with deposition.
We emphasize that Eq. (45) pertains to the imagined situation in which gravitational heating is not involved. This is the same as assuming a spatial Poisson process of deposition, that is, a fixed disentrainment rate keyed to the specific energy state E p .

20
In the presence of heating, however, the length scale of deposition increases relative to l c . That is, heating suppresses the disentrainment rate. The factor α thus modulates the length scale l c so the product αl c is a net e-folding length in the presence of heating. As described below, the factor α is assumed to be a function of the Kirkby number.
Substituting Eq. (45) into Eq. (17) and evaluating the integral then yields where we now redefine the Kirkby number as assuming that β x is independent of E p . Comparing this result with Eq. (33), the energy loss rate due to deposition is the same as the advective cooling rate, but modulated by the factor α.

Conservation of mass revisited 5
The preceding provides the basis for separately calculating the disentrainment rate, consistent with the deposition rate. Because n Ep (E p , x)dE p represents the number of particles within the small energy interval E p to E p + dE p , using Eq. (44) and Eq.
(45) the total disentrainment rate is therefore Thus, the deposition rate is proportional to the cooling rate, as it should be. Here it is important to note that the expected value 1/E p = 1/ E p . In fact, 1/E p is the reciprocal of the harmonic mean (Appendix F). This means that E p 1/E p ≥ 1. 15 Only in the limit where n Ep (E p , x) has zero variance does 1/E p → 1/ E p . To simplify the notation, hereafter we denote the arithmetic mean as E p = E a and the harmonic mean as As a point of reference we may now define an ensemble averaged deposition length as The factor α has a key role in the formulation. As described above, this factor modulates the length scale L c in the presence of gravitational heating. Note that Eq. (48) is equivalent to For a given value of α the length scale L c is set by the cooling rate, and this length scale increases with increasing slope angle θ. But gravitational heating also increases with θ, the effect of which is to suppress the rate of deposition and increase L c .
That is, the deposition length scale is not the same as the cooling length scale. As described below, whereas l c is a measure of the rate of extraction of translational energy, this includes its conversion to rotational energy whose effect is to decrease the likelihood of stopping. On dimensional grounds an inspired guess suggests that this effect is a function f (Ki ) of the Kirkby number Ki . For example, suppose that where µ 1 is a coefficient of order unity. This leads to where α → α 0 as µ 1 Ki → 0. Now, This example suggests that L c → ∞ as µ 1 Ki → 1. That is, µ 1 Ki → 1 sets an upper limit above which deposition is insignifi-10 cant. More generally we may write to indicate the possibility of other dependencies of α on Ki . Note that we provide evidence for this behavior in the companion paper, including the form of Eq. (52) based on experiments. For notational simplicity in subsequent sections, we use α with the understanding that this implies α = α 0 f (Ki ).

15
Here is a key sidebar for reference in our descriptions below of related formulations. We emphasize that according to Eq.
(45) and Eq. (48) the deposition rate is proportional to the advective cooling rate rather than the net cooling rate (the difference between the rates of heating and cooling), where the rate of heating then modulates the deposition rate, therein increasing the deposition length scale. Moreover, the deposition rate explicitly depends on the energy state of the particles. Consider a thought experiment. Let us imagine a system consisting of a box containing a finite number of particles. Suppose that we mechanically 20 add energy to the system such that some proportion of the particles becomes a rarefied granular gas, and suppose that the gas achieves a non-equilibrium steady state with a specific average energy state (Appendix G). This means that the rate of (mechanical) heating is equal to the rate of cooling due to dissipative particle-box collisions, and sublimation (entrainment) matches deposition (disentrainment). That is, depending on the energy state of the particles, deposition occurs even though the difference between the rate of heating and cooling is zero. Now imagine that when a particle is deposited, it cannot become re-25 entrained. The rate of heating and cooling of the remaining gas particles is still the same, yet the deposition process continues for those particles which, by chance, cool to sufficiently low energies for deposition to occur -just as deposition of these particles would have occurred before re-entrainment was "turned off." Furthermore, the average energy of the gas (active) particles can remain fixed while their total energy decreases due to irreversible deposition. Thus, we are assuming that the deposition rate is proportional to the cooling rate rather than the net cooling rate, depending on the energy state of the particles. 30 The effect of heating is to decrease the likelihood of deposition by decreasing the proportion of particles that cool to sufficiently low energies for deposition to occur -which translates to suppressing the disentrainment rate and increasing the length scale of deposition.

Energy and mass balances
With µ = β x /4 tan φ we now collect results from above. The total energy balance is given by To summarize, the first term on the right side of Eq. (55) is due to gravitational heating, the second term is due to frictional cooling, and the last term represents a loss of energy due to deposition. Note that none of these terms explicitly involves the energy E(x). In turn, conservation of mass is given by This is coupled with Eq. (55) via the relation between the total energy E(x), the average energy E a and the harmonic average energy E h (see below), and the explicit appearance of N in both of these equations.
At this point we emphasize that the quantity µ = β x /4 tan φ is not to be interpreted as Coulomb-like dynamic friction coefficient. Indeed, the product mgµ cos θ in Eq. (55) and Eq. (56) looks like an ordinary Coulomb friction force (e.g., Kirkby 15 and Statham, 1975;Gabet and Mendoza, 2012). Recall, however, that cos θ enters from the geometry of particle motions, and does not represent the angle needed to specify the normal component of the weight mg. Similarly, tan φ is an expected reflection angle, not a friction angle. We elaborate these points below.
To close the circle in reference to our stating point, Eq. (12), we now combine Eq. (12), Eq. (55) and Eq. (56) to give This balance involving the average energy E a rather than the total energy E reveals an important behavior associated with deposition, centered on the parenthetical part of the last term. Namely, it is straightforward to show (Appendix F) that E a /E h − 1 ≥ 0. The last term in Eq. (57) therefore represents an apparent heating associated with deposition. With reference to Figure   25 3, a net advective cooling uniformly lowers all particle energy states, thus lowering the average energy E a as well as the total energy E. As this cooling lowers all energy states, some particles enter the range where deposition occurs, and the deposition rate therefore is proportional to the net advective cooling rate. In the absence of a net advective cooling, particles with small energy nonetheless are preferentially disentrained, so the average energy state increases. When cooling and deposition are combined, the average energy decreases more slowly than it otherwise would in the absence of deposition. This effect increases with increasing variance in the distribution of energies (Appendix F), and it vanishes as the variance goes to zero. The balance described by Eq. (57) thus provides a formal description of what we intuitively know: deposition culls lower energy particles, thereby selecting higher energy particles for continued travel with increasing distance. We note that Brilliantov et al. (2018) 5 demonstrate an analogous unexpected behavior of granular gases, namely, the heating of a granular gas associated with particle aggregation with continued loss of total energy. This occurs when the rate of loss of particles by aggregation exceeds the rate of loss of total energy, such that by definition the average particle energy increases.
The balance described by Eq. (57) also reveals an important constraint on particle energies. Namely, if we imagine the special situation of isothermal conditions (dE a /dx = 0), then frictional cooling given by the second term on the right side 10 of Eq. (57) must balance two sources of heating, namely, the first and third terms on the right side. This requires that either Because this latter condition is highly unlikely, an isothermal condition generally requires that Ki < 1. Conversely, net heating must occur with Ki > 1.
According to Eq. (55) or Eq. (57), for a given slope angle θ the spatial rate of net cooling (or net heating) of the ensemble is 15 a fixed quantity in which this slope angle has a dual role. Namely, an increasing slope decreases the rate of frictional cooling by decreasing the expected occurrence of particle-surface collisions, and it simultaneously increases the rate of gravitational heating. With θ = 0, heating vanishes and frictional cooling occurs at a maximum rate of µmg. In turn, as θ → π/2, which represents a vertical cliff, frictional cooling vanishes and heating matches that of free-fall motion. This transition from small to large slopes nicely illustrates what virtually every undergraduate student learns intuitively from the sport of boulder rolling 20 (or "trundling" (Forrester, 1931)), and why this sport is so spectacular and satisfying in steep terrain. Moreover, recall that the Kirkby number Ki = S/µ is the ratio of gravitational heating to advective cooling. If these are balanced, Ki = 1 and Qualitatively, this is the slope at which an undergraduate student may expect that boulder rolling starts to become particularly interesting.

25
The formulation also nicely illustrates that if the heating and cooling rates are matched, this does not imply an absence of deposition, as the last terms in Eq. (55) and Eq. (57) may be finite with Ki = 1. Moreover, because this is a probabilistic phenomenon, some particles are likely to become disentrained even on steep, rough slopes where heating on average exceeds cooling. Experienced undergraduates indeed inform us that some boulders just do not make it all the way to the bottom of the hillslope despite their best efforts to select conditions satisfying Ki > 1.
4 General behavior

Effects of energy balance
There is value in restating Eq. (55), Eq. (56) and Eq. (57) in dimensionless form. Let E a0 denote the initial average particle energy at x = 0 and let N 0 denote the initial number of particles at x = 0. In turn we define a characteristic cooling distance X = E a0 /mgµ cos θ so that E a0 = mgµ cos θX. We now define the following dimensionless quantities denoted by circum-5 flexes: With these definitions we write Eq. (55), Eq. (56) and Eq. (57) as 15 Because the dimensionless disentrainment ratePx(x) = −(1/N )dN/dx, notice that Eq. (61) provides the basis for determining the distribution of travel distances using Eq. (4). This requires specifying howÊ h varies withx for given values of α and Ki . At this point, however, we must confront the fact that we have four unknowns,N ,Ê,Ê a andÊ h , and three equations, one of which is nonlinear in the ratioÊ a /Ê h . Here we add a fourth equation by assuming that this ratio remains fixed, namely, We do not know the distribution fÊ p (Ê p ) required to determine γ (Appendix F). Nonetheless, Eq. (63) essentially assumes that the form of fÊ p (Ê p ) remains similar with distancex. This allows us to illustrate key elements of the formulation. With the assumption of Eq. (63) we note that Eq. (61) becomes and Eq. (62) becomes Focusing initially on Eq. (65), isothermal conditions exist if dÊ a /dx = 0. We then rearrange Eq. (62) or Eq. (65) to define a transition value of the Kirkby number, namely, in this case cooling occurs with Ki < 1 and heating occurs with Ki > 1. The ratioÊ a /Ê h = γ generally increases with the variance of E p , thus decreasing Ki * . That is, as this variance increases, the transition between cooling and heating occurs at a smaller value of the Kirkby number. This represents a stronger culling (deposition) of lower energy particles. The largest 10 possible transition value is Ki * = 1.
We now start with an idealized example that illustrates key elements of the formulation, including the coupling between Eq.
(60), Eq. (61) and Eq. (62). Assume that the Kirkby number Ki is fixed, and assume isothermal conditions. Thus dÊ a /dx = 0 withÊ a =Ê a0 so that Eq. (62) leads tô Eq. (61) and Eq. (67), In turn, using Eq. (4) this yields an exponential distribution of travel distances with mean 20 so that Note that an increasing value of γ in Eq. (69) represents an increasing proportion of lower energy particles available for deposition relative to this availability with γ → 1, the effect of which is to decrease the mean travel distance.
The total energyÊ also declines exponentially withx. Namely, substituting Eq. (70) into Eq. (60) leads to This example of isothermal conditions illustrates that withÊ a =Ê a0 , then according to Eq. (69) the average travel distance µx is directly proportional to the initial average energy. However, isothermal conditions are unlikely, because according to Eq. (66), such a condition requires a specific value of Ki for the ratioÊ a /Ê h = γ. We now consider the more general case 5 involving either net cooling or net heating.
As above we assume that the ratioÊ a /Ê h = γ is fixed, although the averagesÊ a andÊ h otherwise are unconstrained. Net cooling or net heating is not prescribed; either condition is allowed. Using Eq. (60) and Eq. (61) the disentrainment rate is Note that as a → 0 the disentrainment rate goes to a fixed value equal to 1/b = γ/αÊ a0 , and the distribution fx(x) of travel 15 distancesx goes to an exponential distribution with mean µx More generally, the distribution of travel distances is a generalized Pareto distribution with position parameter equal to zero (Appendix H), namely, 20 where now a ∈ is interpreted as a shape parameter and b > 0 is a scale parameter (Pickands, 1975;Hosking and Wallis, 1987). The cumulative distribution iŝ and the exceedance probability iŝ For a < 1 the mean is which is independent of the ratio γ =Ê a /Ê h . This is the same as Eq. (69) for isothermal conditions, although the denominator in Eq. (77) generally is not equal to γ. In turn, Eq. (77) requires that 5 which provides the upper limit of Ki for which the mean µx is defined. Because α > 0, this limit may be greater than one. For a < 1/2 the variance is 10 Unlike the mean, the variance depends on γ =Ê a /Ê h . In turn, for a ≥ 1 such that the mean of fx(x) is undefined. Moreover, for a ≥ 1/2 the variance is undefined. These results reflect the heavy-tailed behavior of the generalized Pareto distribution.
As a point of reference in the second companion paper (Furbish et al., 2020a), the generalized Pareto distribution defined 15 by Eq. (74) also may be considered a generalized Lomax distribution. This distribution can be rewritten as an ordinary Lomax distribution (Appendix H). Namely, if we define the shape parameter a L = 1/a and the scale parameter b L = b/a, then Eq.
(74) becomes which is a Lomax distribution with mean For a L > 0 (a > 0) the forms and behaviors of Eq. (74) and Eq. (81) are identical. Notice, however, that if a < 0 then a L = 1/a < 0 and b L = b/a < 0 for positive b. This means that we cannot use the form of the Lomax distribution given by Eq. (81) to examine conditions involving a < 0. Yet these conditions are mechanically meaningful, so we proceed using the generalized Pareto distribution given by Eq. (74). To be clear, the ordinary Pareto distribution that is normally referred to in the literature is 25 a special case of the generalized Pareto distribution. In turn the Lomax distribution is a special case of the Pareto distribution (and therefore of the generalized Pareto) with position parameter equal to zero. This represents a condition of rapid thermal collapse. Specifically, when a < −1 this distribution monotonically increases and becomes asymptotically unbounded atx = b/|a|. In the limit of a → −1 it becomes a uniform distribution. When a = −1/2 this distribution is triangular. For −1/2 < a < 0 this distribution decays more rapidly than an exponential distribution and is bounded at the positionx = b/|a|. For a = 0, fx(x) becomes an exponential distribution, representing an isothermal condition 5 as described above. For a > 0 the distribution fx(x) is heavy-tailed. This represents a condition of net heating. Specifically, for 0 < a < 1/2 this distribution decays more slowly than an exponential distribution, but it possesses a finite mean and a finite variance. For 1/2 ≤ a < 1 the distribution possesses a finite mean but its variance is undefined. For a ≥ 1 the mean and variance of fx(x) are both undefined, even though this distribution properly integrates to unity. For a > 0 the tail of fx(x) decays as a power function, namely, fx(x) ∼x −(1+1/a) . The exceedance probability decays as Rx(x) ∼x −1/a . These results 10 are summarized in Table 1. We provide evidence of all three behaviors -rapid thermal collapse, isothermal conditions, and net heating -in our second companion paper (Furbish et al., 2020a).
The formulation above assumes uniform surface conditions, specifically, uniform slope angle and roughness texture. We show below (Section 6) how it may be adapted to varying downslope conditions. We also note that the distribution fx(x) given by Eq. (74) can be incorporated into a mixed distribution. Indeed, a mixed distribution is the natural choice for describing 15 the travel distances of a mixture of particle sizes, each involving a different frictional cooling behavior for a given surface roughness (Roth et al., 2020). Table 1. Behavior of the generalized Pareto distribution associated with the shape parameter a and Kirkby number Ki as illustrated in Figure   4.

Elements of the average travel distance
The average travel distance given by Eq. (77) for Ki < 1 + 1/α contains all of the elements that influence particle motions except the quantity γ. Thus, whereas the average by itself does not reveal the source of variations in the form of distribution of travel distances, Eq. (74), the average nonetheless provides a focal point. Here we rewrite this average in its dimensional form, then step through the significance of its elements. Namely, with E a0 = (1/2)m u 2 0 and Ki = S/µ = (1/µ)mg sin θ/mg cos θ, For an ensemble of particles whose motions start at x = 0, the average travel distance µ x increases directly with the average starting energy E a0 ∝ u 2 0 . This is entirely akin to the formulation by Kirkby and Statham (1975) (see below), and it highlights 10 the significance of the initial particle energy conditions at x = 0 in setting their travel distances. The archetypal example involves rock fall from cliffs followed by their motions over talus and scree slope surfaces (Figure 1), where fall heights and initial rebounds set the initial average downslope energy. This also is a key element in experiments where initial energies are set by the choice of drop height (Kirkby ana Statham, 1975) or launch speed (Gabet and Mendoza, 2012;DiBiase et al., 2017).
This aspect of the formulation also points to the significance of energetics associated with the entrainment rate E s (x) in Eq.

15
(5) and Eq. (6) at hillslope positions that are not necessarily as well-defined as, say, the base of a cliff (see Section 6).
The average travel distance µ x is inversely proportional to the rate of frictional cooling represented by µgµ cos θ. Here we reemphasize that despite its form, this expression does not represent a Coulomb-like friction. Rather, this expression enters the formulation via the characteristic length λ in setting the expected number of collisions per unit distance, n x . As described below, the surface-normal component of the particle weight does not set collisional friction; this is set by dynamic forces during collisions. Moreover, the appearance of the Kirkby number Ki in the denominator of Eq. (83) indicates that as Ki increases, the denominator becomes smaller (subject to the conditions that Ki < 1 + 1/α), so the average travel distance increases. We also note that, except for purely bouncing motions, it is incorrect to interpret the length λ strictly as a saltation-like distance.
This is a scaling approximation (Appendix D) to show that n x must involve the average energy (∝ u 2 ) and the geometrical factor cos θ at lowest order.

5
Notice that Eq. (83) indicates that with E a0 = (1/2)m u 2 0 the average travel distance µ x is independent of the particle mass m. Viewed in isolation, this suggests that large particles should on average travel no farther than small particles. However, this is inconsistent with what is observed in laboratory and field-based experiments (Kirkby and Statham, 1975;DiBiase et al., 2017;Roth et al., 2020) and with downslope size sorting on natural talus and scree slopes (Statham, 1976). We examine this topic in the second companion paper (Furbish et al., 2020a); here we offer a synopsis, which centers on the interpretation and 10 significance of the quantities µ and α.
Recall that the formulation is based on the assumption that a change in translational kinetic energy ∆E p associated with a particle-surface collision can be expressed as ∆E p = −β x E p so that β x = −∆E p /E p is the proportion of the energy extracted during the collision. Both ∆E p and β x are random variables. As described in Appendix E, in general we may write the energy balance of a particle as Here, a positive change in rotational energy ∆E r is seen as an extraction of translational energy. This loss of translational energy with the onset of rotation may be relatively large if a collision involves stick following initial sliding due to a large normal impulse, and such a loss also may occur due to the imposed torque of friction during a collision that does not necessarily involve stick. The term f c in Eq. (84) represents losses associated with particle and surface deformation as well as work 20 performed against friction during collision impulses (thence converted to heat, sound, etc.). But this term also includes losses associated with deformation of the surface at a scale larger than that of an idealized particle-surface impulse contact, namely, and rotation oriented differently than any incident rotation. In some cases the change in energy ∆E p can be expressed directly in terms of the energy state E p (Appendix E). However, the complexity of particle-surface collisions on natural hillslopes precludes explicitly demonstrating such a relation for all possible scenarios. Nonetheless, it is entirely defensible to assume that energy losses can be related to the energy state E p if the elements involved are formally viewed as random variables. Then, the simple relation ∆E p = −β x E p is to be viewed as an hypothesis to be tested against data. 30 This hypothesis formally enters the formulation via the right side of Eq. (58). Namely, from this relation we may write µ ∼ β x , highlighting that µ is associated with the cooling rate. In turn, particle collision mechanics (Appendix E) suggest, for example, that µ ∼ β x ∼ M (θ), where M (θ) involves the coefficients associated with tangential and normal impulses contributing to energy losses during collisions, and depends on the slope angle θ in that the expected surface normal impact velocity varies with this angle. (In an idealized particle-surface collision these coefficients include the normal coefficient of restitution and a coefficient describing the ratio of tangential to normal impulses during the collision (e.g., Brach, 1991;Dunn, 1992, 1995)). Moreover, M (θ) is independent of particle size.
In turn, focusing on the definition of the deposition length scale L c , Eq. (51), α may be viewed as representing a direct effect of heating described by mg sin θ, namely, to decrease the likelihood of deposition by decreasing the proportion of 5 particles that cool to sufficiently low energies for deposition to occur -which translates to suppressing the disentrainment rate and increasing the length scale of deposition L c relative to the cooling length scale l c = E h /mgµ cos θ. Specifically, heating decreases the spatial rate of the Poisson deposition process below that which would occur in the absence of heating. In this view, µ goes with the cooling rate (not the deposition rate). But we also may write Eq. (51) as in Eq. (53). Viewed in this manner, we may define an apparent friction factor as µ 0 = µ(1 − µ 1 Ki ) associated with deposition. Here again, µ is associated 10 with the cooling rate but is then modulated by heating. We suggest in the second companion paper (Furbish et al. 2020a) that for the same particle size, α increases with increasing Ki , very likely due to a combination of increased heating and increased partitioning of translational energy to rotational motion (Dorren, 2003;Luckman, 2013) -both decreasing the likelihood of stopping and not represented in just the factor µ. We also suggest that for the same slope and surface roughness, α increases with increasing particle size, decreasing the likelihood of frictional loss with increasing rotation. 15 Turning to the factor γ (which does not appear in Eq. (83) large. It is possible to qualitatively explore possible values of γ based on assuming different forms of the distribution of energies (and we have done this), but in the absence of knowing the specific form of the distribution, this exercise is not particularly meaningful. In the companion paper we show that fits of experimental travel distances to the theoretical distribution f x (x) are relatively insensitive to the specific value of γ selected.

Related formulations 25
Here we briefly examine three related formulations of particle disentrainment, focusing on the mechanical basis of this work for comparison with the formulation above. (We examine associated experiments in the companion paper.) We start with the formulations of Kirkby and Statham (1975) and Gabet and Mendoza (2012). These begin with descriptions of particle motions over time rather than space, centered on consideration of a combination of momentum and energy. We then consider elements of the probabilistic formulation presented by Furbish and Haff (2010), Furbish and Roering (2013) and Doane et al. (2018).

Kirkby-Statham formulation
In their study of particle motions on scree surfaces, Kirkby and Statham (1975) start with a statement of conservation of energy for a particle falling from height h at x = 0 onto a rough surface inclined at an angle θ. Namely, if w 0 = √ 2gh is the vertical impact velocity, then it is assumed that the initial downslope velocity on average is u 0 = w 0 sin θ. (This actually should be u 0 = w 0 sin θ with the normal coefficient of restitution .) The initial downslope particle energy therefore is (1/2)mu 2 0 = 5 2 mgh sin 2 θ = E p0 . In turn, because work is W = F x l, where F x is the downslope force and l is a displacement, then for a fixed force F x the displacement is l = W/F x . Assuming that W must be equal to the initial kinetic energy E p0 (that is, this initial energy is dissipated over the distance l), then l = E p0 /F x . Assuming a Coulomb-like friction behavior, F x = mg sin θ − µ d mg cos θ, where µ d is a dynamic friction coefficient. Upon asserting that the length l represents the expected travel distance where it is assumed that |µ d cos θ| ≥ sin θ. As described below, this is equivalent to assuming that particle energy decreases linearly with distance.
This formulation correctly describes the motion of an individual particle that experiences a fixed Coulomb-like friction force, but it cannot represent the rarefied behavior of an ensemble of particles. Nonetheless, it shares elements of the preceding 15 formulation. Namely, in comparing Eq. (85) with Eq. (83), let us momentarily set aside the fact that E p0 cannot represent E a0 except in the limit of zero variance of initial energy states (γ = 1), and that the friction factor µ in Eq. (83) and the dynamic friction coefficient µ d in Eq. (85) have different interpretations. These two descriptions of the average travel distance µ x then converge in the limit of α → ∞. Inasmuch as Eq. (52) correctly describes the behavior of α, this limit coincides with Ki → 1.
More generally, Eq. (85) implies that the deposition rate is independent of the extant energy state of particles. If Eq. (85) 20 denotes the average of a distribution of travel distances with fixed disentrainment rate, then this fixed rate P x = 1/µ x . In dimensionless form this is That is, Eq.(86) cannot allow for the possibility of variations in the cooling rate or heating rate that give spatial variations in the disentrainment rate, as in Eq. (72). The resulting distribution f x (x) of travel distances therefore is exponential for all Ki < 1.

25
Interestingly, the formulation of Kirkby and Statham (1975) is equivalent to (Appendix I) This is like Eq. (57), but lacks the heating effect of deposition. Like Eq. (85), Eq. (87) implies that deposition is independent of the extant energy state; and when Ki = 1 the energy E p remains fixed as x → ∞ (again noting that E p0 cannot represent E a0 except in the limit of zero variance of initial energy states). Dorren (2003) provides a review of efforts to elaborate the Kirky-Statham description of particle motions in relation to hazard assessment. These mostly appeal to a Coulomb-like frictional behavior and are focused on predicting rockfall runout distances.

Gabet-Mendoza formulation
In support of their experimental work involving particle motions on a rough, inclined surface, Gabet and Mendoza (2012) 5 appeal to ideas from Quartier et al. (2000) and Samson et al. (1998) and assume that particle acceleration is described as a linear combination of the gravitational force, a Coulomb-like friction and collisional friction, namely, As written, the dimensions of the coefficient κ depend on the value of the exponent ψ, which is thought to vary between one and two based on experiments. The principal significance of this formulation is that it points to the idea of collisional friction, 10 thus representing an important step beyond the Coulomb-like model of Kirkby and Statham (1975). However, because there is confusion in the literature regarding the form and interpretation of Eq. (88), we summarize the basis of this formulation in Appendix J. The essence is this: The Coulomb term and the collisional term as written are not additive for an individual particle. The collisional term is a stochastic quantity and applies to an averaged behavior, not to the instantaneous behavior of an individual particle. If this term is involved, the velocity u must be considered a time-averaged or ensemble-averaged 15 velocity, or Eq. (88) must be recast as a Langevin-like equation. Parts of this formulation are appropriate for describing the behavior of particles that roll bumpety-bump over a surface roughened with a monolayer of particles, but this formulation is problematic in its description of the mechanics involved in rarefied motions over the roughness of natural hillslopes.
In both formulations above the idea of a Coulomb-like friction with a dynamic friction coefficient is mechanically unsound (Appendix J). For particles that tumble, bounce and skitter down a rough surface, the static normal weight of a particle, 20 mg cos θ, does not set the particle-surface friction. Rather, dynamic forces during collision impulses matter (Brach, 1991;Stronge, 2000). This includes the dynamic Coulomb-like friction force associated with conversion of translational to rotational kinetic energy during collisions (Appendix E). Formulating a dynamic friction coefficient would require ensemble averaging of the ratio of tangential to normal momentum exchanges, both of which are random variables. A Coulomb-like friction is appropriate for solid body and dense granular motions, but not for the rarefied conditions described here.

Furbish-Haff-Roering-Doane formulation
The probabilistic formulation presented by Furbish and Haff (2010), Furbish and Roering (2013) and Doane et al. (2018) assumes that travel distances are described by an exponential distribution whose mean µ x is a function of the local slope S.
Namely, the mean increases with S and becomes unbounded as S approaches a critical value S c . This formulation is equivalent to setting the mean µ x ∼ L c . Here we consider the behavior of L c over small S then as S → S c .
Starting with Eq. (49) we write If α is described by Eq. (52), and neglecting the factor µ 1 for simplicity, then this is A binomial expansion of Eq. (89) gives and Eq. (90) gives 10 From Furbish and Haff (2010), If we interpret the length scale λ 0 ∼ αE h /mgµ, then for small to modest slopes S, Eq. (91) and Eq. (93)  in Section 2. Recall that in this frame of reference the particle travel distance is denoted by r and the starting position may involve any position x. Then, with reference to Eq. (4), Eq. (5) and Eq. (6), the disentrainment rate is P r (r; x), the distribution of travel distances r is f r (r; x) and the exceedance probability (survival function) is R r (r; x) = 1 − F r (r; x). In turn, for particles starting at x, the mean travel distance is µ r (x).
To use the results of Section 2.1 in specifying the exceedance probability R r (r; x) and the probability density f r (r; x) in the entrainment forms of the flux and the Exner equation, Eq. (5) and Eq. (6), requires a key assumption. Namely, one must assume 5 that the factors controlling the disentrainment rate on a hillslope change sufficiently slowly over x such that these factors defined at any position x correctly determine the conditions for the downslope motions of particles starting at x (Furbish and Roering, 2013;Doane et al., 2018). This is equivalent to assuming that during its downslope motion a particle "sees" conditions similar to those at its starting position. However, in actuality particles may see new conditions during their motions that change their behavior relative to what was "predicted" by the conditions at their starting positions. Let λ S denote a characteristic distance 10 over which conditions persist. For example, focusing on the Kirky number Ki , Thus, a rapid change in Ki over position x implies that λ S is small, and if Ki changes slowly then λ S is large. Uniform conditions imply that λ S → ∞. We may then assume that if µ x λ S , conditions change sufficiently slowly that use of the continuous forms of R r (r; x) and f r (r; x) with Eq. (5) and Eq. (6) provides a reasonable approximation of collective particle 15 behavior.
This strategy might be acceptable for an exponential-like distribution with finite moments, but it is problematic if particle travel distances r involve a heavy-tailed distribution or if conditions transition along x between net cooling and net heating, or vice versa. Herein resides the merit of the discrete form of the disentrainment rate and the distribution of travel distribution as summarized in Section 2.2. Recall that this formulation is aimed at describing the ingredients of disentrainment that are 20 occurring at different locations on a hillslope, where the mechanical behavior at a location transitions to another behavior in the downslope direction. In this formulation we let p k denote the probability that a particle, having not been disentrained before the kth interval, then becomes disentrained within this interval.
Let dr denote a finite (rather than infinitesimal) interval. Then the kth interval begins at r and ends at r + dr. Letting N k = N (r) denote the number of particles reaching the kth interval, then based on Eq. (56) the probability that a particle will 25 be disentrained within this interval is This will be recognized as the setup for a simple finite-difference scheme, to be coupled with a similar finite-difference expres- where both the Kirky number Ki and the elements of the transition value of the Kirkby number Ki * = 1 − γ/α + 1/α may vary from one interval to the next as conditions vary in the downslope direction. The proportion ofN (0;x) particles starting from positionx is then recovered from We note that, although different in form and implementation, this description is similar to the particle-based scheme of Tucker and Bradley (2010) in which particle behavior adjusts to newly encountered conditions during downslope motion.

10
Consider for illustration a situation in which the Kirkby numbers Ki and Ki * systematically vary with positionx, relative to uniform conditions ( Figure 5). This may be due, for example, to variations in steepness or in the friction µ with increasing Figure 5. Cartoon of hillslope surfaces with downslope variations in steepness leading to concomitant variations in heating, cooling and deposition; this is in contrast to a planar slope with uniform µ that produces either net heating or net cooling, or isothermal conditions. travel distance. Also recall that Ki < Ki * implies cooling whereas Ki > Ki * implies heating.
In these examples we let α vary with the Kirkby number Ki according to Eq. (52) in anticipation of results presented in the companion paper. A decreasing rate of cooling associated with, for example, steepening in the downslope direction generally 15 increases the heaviness of the tail of the distribution relative to the tail associated with a fixed rate of cooling ( Figure 6). We present evidence in the companion paper that this occurs in the field-based experiments reported by DiBiase et al. (2017).
Specifically, particles were launched down a rough hillslope surface, and then their travel distances were measured over a 14 m interval. Particles reaching the steeper slope below the measurement interval continued to the base of the hillslope without stopping. In turn, an increasing rate of cooling (e.g., with decreasing slope in the downslope direction) generally lightens the tail, and may lead to truncation of the distribution if the rate increases rapidly enough. Moreover, a condition involving initial heating followed by cooling (e.g., with a concave hillslope surface) can lead to a distribution with a finite mode (Figure 7).
These examples represent situations where particle travel distances cannot necessarily be approximated by a distribution whose parametric values are set by the hillslope conditions at the position where particle motions start. We defer further examination of this behavior, including use of the convolutions in Eq. (5) and Eq. (6), for a later time.

Discussion and conclusions
Our formulation of rarefied particle motions is based on a description of the energy balance of a cohort of particles treated as a rarefied granular gas, and a description of particle deposition that depends on the energy state of the particles. The net heating lead to a heavy-tailed distribution of travel distances. We provide compelling evidence of all three behaviors in our companion paper (Furbish et al., 2020a). Here we emphasize that we do not choose the generalized Pareto distribution in the empirical manner of selecting a distribution based on goodness-of-fit criteria applied to data sets. Rather, this distribution is dictated by the physics of the problem, just as, for example, the Boltzmann distribution (an exponential distribution) emerges in classical statistical mechanics from consideration of the accessible energy microstates of a gas system. We elaborate this 15 point in the third companion paper (Furbish et al., 2020b).
Two of the most important elements of the formulation are the deposition length scales l c (E p ) and L c (E h ), the former being keyed to the specific particle energy state E p and the latter being keyed to the harmonic average energy E h of the particle cohort. Indeed, these lengths provide the essential connection between particle deposition and the energy balance of the particle cohort. We assume that l c is set by the advective cooling length scale in the Fokker-Planck equation, that is, Eq. (44). This is 20 a natural choice in that deposition must go with cooling. The energy specific deposition rate in the absence of heating is then specified as if deposition proceeds as a spatial Poisson process. We emphasize that this represents a maximum (information) entropy choice in the sense that it is faithful to what we think we know, namely, the connection between deposition and cooling, as well as to what we do not know (Jaynes, 1957a(Jaynes, , 1957b, namely, any detailed physics that would produce a different rate (for example, involving a nonlinear dependence on energy state) but which cannot be specified or constrained with available 25 information. This description then leads to the interesting result, Eq. (46), that the loss of total energy due to deposition appears to be independent of the energy state. In particular, the loss of large energy states occurs at a relatively slow rate whereas the loss of small energy states occurs at a relatively fast rate. In effect the rate of loss of energy per energy interval is fixed across energy states. The result is that the energy E p cancels with substitution of Eq. (45) into the integral in Eq. (46) such that the total loss becomes independent of the energy state. That is, the loss of total energy goes simply with the loss of particles (and 30 the energy they possess).
In turn, the total deposition rate is energy dependent. This rate, defined by the length scale L c , is obtained by integrating the number density of particles over all possible energy states as in Eq. (48). Because l c is keyed to the energy state E p , but the integral in Eq. (48) does not involve this energy in the numerator, the result involves the reciprocal of the harmonic average energy E h . In general, the harmonic average diverges from the arithmetic average E a with increasing variance of the distribution of energy states. With E a /E a = γ, the resulting ratio γ/E a in Eq. (56) (with dimensionless form given by Eq. (61) or Eq. (64)) reflects an increasing proportion of lower energy particles available for deposition, relative to this availability with γ → 1. This effect is directly apparent in the expression of the mean travel distance, Eq. (69), associated with isothermal 5 conditions.
Note that the formulation does not involve specifying a threshold energy for deposition. Such an idea is mechanically irrelevant. Whereas low energy particles are on average more likely to become disentrained than are high energy particles, a set of particles with precisely the same low energy will for probabilistic reasons not be disentrained simultaneously. Each particle experiences a unique set of conditions that disentrain it; and because of this uniqueness of conditions a particle with energy 10 below an arbitrarily assigned threshold can with finite probability be gravitationally reheated to a higher energy state. For given particle and surface roughness conditions, the formulation treats this aspect of disentrainment as a probabilistic process. In effect, this aspect is incorporated into the deposition lengths l c and L c as these are related to the distribution of particle energy states and the probabilistically expected extraction of energy during collisions.
Frictional cooling is formulated in terms of extraction of translational kinetic energy associated with particle-surface col- 15 lisions. This involves the random variable β x = −∆E p /E p whose energy specific average β x is the expected proportion of energy extracted from particles with energy E p . In detail the change in energy ∆E p may be partitioned between a frictional loss, any conversion of translational to rotational energy, and any apparent loss associated with downslope incident motion reflected to transverse motion during a glancing collision. Our treatment of β x as a random variable does not distinguish the details involved in collisions. Yet these details may be important in terms of effects of different particle sizes and shapes, 20 specifically the likelihood that the partitioning of energy losses differs between sizes or shapes. Herein the quantity α has a dualistic role. As incorporated in Eq. (51), this quantity represents the effect of heating, namely, to decrease the likelihood of deposition by decreasing the proportion of particles that cool to sufficiently low energies for deposition to occur -which translates to suppressing the disentrainment rate and increasing the length scale of deposition L c . As incorporated in Eq. (53), this quantity modulates the frictional cooling described by µ ∝ β x to give an apparent decrease in friction associated with 25 deposition.
Whereas particles that are small relative to the surface roughness texture are on average more likely to experience near collinear collisions with surface bumps and be "captured" within divots and pockets, particles that are large relative to the roughness texture are less likely to experience direct collisions with, or strong deflections by, smaller surface bumps. In addition, large particles are more likely to experience conversion of their translational energy into rotational energy with less 30 loss during collisions. In particular, large spherical particles are more likely to roll or spin with increased heating, and large spinning particles are less likely than are smaller particles to be frictionally cooled. These points are reflected in the laboratory experiments of Samson et al. (1998Samson et al. ( , 1999) (Appendix J), the laboratory experiments of Kirkby and Statham (1975) and the field experiments of DiBiase et al. (2017) and Roth et al. (2020) (see the second companion paper (Furbish et al., 2020a)). This also implies that for a given slope angle and surface roughness, some particle sizes may experience net cooling while some sizes experience net heating (Roth et al., 2020), likely contributing to the size sorting observed on many talus and scree slopes (Kirkby and Statham, 1975;Statham, 1976;Luckman, 2013). We suspect the noticeable sorting in Figure 1 is due to these effects.
The formulation readily accommodates the idea of a mixed distribution composed of different distributions associated with different particle sizes or mechanical behaviors. This amounts to forming a sum of distributions, each weighted in proportion to 5 the size classes involved in transport. As with individual sizes, the formulation assumes rarefied conditions -that particles of different sizes do not interact during their downslope motions, or that such interactions negligibly influence the particle energy balance relative to particle-surface interactions. We provide an example in the second companion paper (Furbish et al., 2020a).
With rockfall and subsequent particle motions over talus and scree surfaces, the initial energy state E a0 can be approximated in terms of the fall height (Kirkby and Statham, 1975). But this is a special situation in which the initial energy can be reasonably 10 constrained. More generally, and with reference to the entrainment forms of the flux and the Exner equation, Eq. (5) and Eq.
(6), we are concerned with entrainment of particles from many if not all positions on a hillslope in relation to disturbances. This points to the idea that entrainment, if followed by long distance motions, requires sufficient initial heating to keep particles moving downslope. This in turn echoes the conclusion of Doane et al. (2018a), that correctly specifying the entrainment rate is a key part of implementing formulations of nonlocal transport and mass conservation. Because of the significance of sediment 15 capacitors (e.g., vegetation) in trapping and storing sediment on hillslopes Doane, 2018a), there is merit in clarifying the initial energetics of particles upon their release (i.e., entrainment) from storage. There also is a need to examine re-entrainment and transport associated with particle collisions, analogous to work on particle splash during aeolian transport and the energetics of collective entrainment (Ancey et al., 2008) by collisions during bed load transport (Lee and Jerolmack, 2018). 20 That the energy and mass balances are expressed in the form of coupled differential equations opens the possibility of describing effects of varying disentrainment rates in response to changing downslope conditions in a manner intrinsic to particlebased treatments of transport (Tucker and Bradley, 2010), but not readily incorporated in previous probabilistic descriptions.
Namely, if surface conditions change in the downslope direction, for example, giving net cooling followed by heating or vice versa ( Figure 5), then particles whose travel distances are large enough "see" this change and their behavior concomitantly 25 changes. In this case the coupled equations of energy and mass in principle can be solved to accommodate these changing conditions. Interestingly, as differential (or finite difference) equations these have a local form, yet they intrinsically represent nonlocal behavior in that information concerning the energy state E a and the mass N is cumulatively handed from one position to the next downslope. In turn, the forms of R r (r; x ) and f r (r; x ) associated with any position x in the expressions of the flux and its divergence, Eq. (5) and Eq. (6), must be based on information downslope from this position. 30 In this regard, here we offer further perspective on what is meant by local versus nonlocal transport on hillslopes. A transition of travel distances involving a distribution with a light tail to one with a heavy tail, as embodied in the generalized Pareto distribution, does not distinguish local from nonlocal transport. As fully explained in Furbish and Roering (2013) and in Furbish et al. (2016), the convolutions in Eq. (5) and Eq. (6) represent nonlocal transport regardless of the form of the probability density function f r (r; x) and its associated exceedance probability function R r (r; x). These scale independent expressions are just specialized forms of the Master equation used in probabilistic descriptions of particle motions over a remarkable range of scales (Einstein, 1905;von Smoluchowski, 1906;Chandrasekhar, 1943Risken, 1984. Nonlocal transport is a physical thing, and refers to the idea that attributes of particle motions used in defining the rheology, the flux or its divergence at a position x depend on conditions "far" from this position (e.g., Bocquet et al., 2009;Brantov and Bychenkov, 2013;Henann and Kamrin, 2013). In contrast, local transport is a mathematical thing, not a physical thing, and refers to the idea that under 5 certain circumstances the convolution form of the Master equation can be approximated such that the flux or its divergence has the form of a local mathematical expression -for example, a Fokker-Planck equation -whose terms involve conditions associated with the local position x. As alluded to above, a local expression can be formulated when the distribution f r (r; x) has finite moments and is peaked near the origin (r = 0). A heavy-tailed behavior means that this is not justified. Rather, the full convolution or a fractional derivative approximation of it must be used (Schumer et al., 2009). Because of the generality 10 and scale independence embodied in the Master equation and the convolutions in Eq. (5) and Eq. (6), the use of "nonlocal" as a qualifier of "transport" in reference to hillslopes actually is redundant (Doane, 2018). Its use is merely a reminder that the flux or its divergence at a position x depends on things happening upslope.
The entrainment rate E s (x), the exceedance probability function R(r; x) and the distribution of travel distances f r (r; x) within the integrals in Eq. (5) and Eq. (6) (2013) and Furbish et al. (2016Furbish et al. ( , 2017Furbish et al. ( , 2018b. The formulation may have interesting implications for examining Martian landforms. For example, the appearance of the acceleration g in Eq. (51), Eq. (56) and Eq. (83) immediately suggests the possibility that particle travel distances are on average significantly longer on Mars than on Earth for otherwise similar particle sizes and surface-roughness conditions; and 30 we are confident in suggesting that future Martians likely will have far more fun than Earthlings in the sport of boulder rolling, notably on the crater rim of Olympus Mons. Nonetheless, we leave it to folks more familiar with Mars than we are to examine this. A key element of doing this is to either assume that the friction factor µ is similar to what occurs on Earth (which may be entirely reasonable) or further unfold the elements of this factor. We comment on this idea again in the second companion paper (Furbish et al., 2020a). We meanwhile note that a similar question arises in relation to the role of g in setting the friction of granular slopes on Mars. Atwood-Stone and McEwen (2013) address this question by examining dune slip-face angles on Mars, and suggest that the similarity of these angles with those observed on Earth weakens any argument for different granular behavior associated with g -consistent with independent assessments (Moore et al., 1987;Tesson et al., 2020) and the idea that this angle is set by the static granular force-chain network (Cates et al., 1998;Furbish et al., 2008).
Appendix A: Choice of terminology 5 The study of granular materials is concerned with the behavior of the phases of these materials and associated phase transitions (Jaeger et al., 1996;Baldassarri et al., 2005;Daniels and Behringer, 2006;Forterre and Pouliquen, 2008;Jerolmack and Daniels, 2019). These phases and transitions share attributes with ordinary materials -solids, liquids and gases -although granular materials often exhibit behavior that is much different than ordinary materials. Nonetheless, it has become customary in the study of granular materials to adopt terminology similar to that used to describe ordinary materials.

10
The ideas of heating and cooling of a granular material are straightforward, to mean a change in the granular temperature of the material, specifically the average translational kinetic energy of the particles (but see van Zon and MacKintosh (2004) and Baldassarri et al. (2005)). However, granular materials do not possess an internal energy in the sense that we attribute to the particles of an ordinary liquid or gas. This means that heating of a granular material requires a mechanical input of energy, whereas cooling is associated with dissipative (non-conservative) collisions of particles with each other and with boundaries. 15 In the problem at hand, gravitational heating occurs as particles move downslope, and their gravitational potential energy is converted to kinetic energy. Frictional cooling is associated with dissipative particle-surface interactions (e.g., collisions).
The ideas of melting and freezing of a granular material (Daniels and Behringer, 2006) pertain to the transition between a solid-like phase and a hydrodynamic (fluid-like) phase. However, in the problem at hand, we are concerned with rarefied particle conditions in which disentrainment from the rarefied state to the solid-like state or vice versa does not involve an 20 intermediate hydrodynamic phase (e.g., Haff, 1983;Jenkins and Savage, 1983;Jaeger et al., 1996). Entrainment is akin to sublimation, and disentrainment is akin to deposition (or desublimation). Phase transitions involving an intermediate hydro- perspectives on this emerging topic, notably in relation to transport by shear flows.

Appendix B: Particle cohort
In order to clarify the idea of a cohort of particles associated with a control volume with edge length dx (Figure 2), here we offer a straightforward thought experiment. As a point of reference, the study of granular gases typically involves consideration of the behavior of an individual system composed of many particles that are mechanically heated, where energy dissipation 30 is associated with particle-particle collisions. In contrast, our problem involves an unusual situation in that we must start by considering a system composed of one particle, where energy dissipation occurs with particle-surface collisions, and then in turn consider the behavior of an ensemble of such systems.
Imagine a box containing one particle. We mechanically shake the box and the particle is heated. At any instant the particle has kinetic energy E p . Each time the particle collides with the floor of the box it is re-heated, and each time it collides with a wall of the box energy is extracted. Eventually the particle by chance has sufficiently low energy that when it next encounters 5 a wall it becomes irreversibly deposited (disentrained) onto the wall. Then the box has no moving particle.
Like Gibbs (1902), we now imagine a great number N of nominally identical but independent single-particle systems, where each particle in each system (box) behaves according to the same laws of physics, each undergoing heating and collisional cooling, and occasionally being deposited ( Figure B1). We now choose one instant in time and examine the state of each Figure B1. Schematic diagram of surface inclined at angle θ and control volume with edge length dx through which particles move, with Gibbs-like ensemble of single-particle systems leading to definition of the cohort of N (x) particles starting at the left face of the control volume.
particle. Some particles previously have been deposited, so at this instant N refers to those systems whose particles are in 10 motion. At this instant each particle has kinetic energy E p , and we may define the ensemble probability density f Ep (E p ) of energy states E p . As a consequence we also may at this instant define the ensemble averaged kinetic energy E p and the total energy E = N E p . (Alternatively, we could imagine all N particles in a single box at one instant, but with the caveat that we must imagine them as not interacting with each other, only with the floor and walls of the box.) We now choose a successive instant in time, namely, t + dt. During dt the number N has decreased with deposition of some particles, the distribution 15 f Ep (E p ) has changed, and the average energy E p and the total energy E have changed.
More generally we can choose N different instants in time t, one instant for each box, and examine the state of each particle. Then, upon collecting the particles as a cohort independently of the selected times, like above we observe an ensemble distribution of energy states with specific average energy E p and total energy E. At this point we relax the idea of a box, and simply view particle-wall collisions more generally as particle-surface collisions during motions parallel to x; and instead of heating the particles via particle-floor collisions we imagine this to occur continuously by gravitational heating. We then 5 let the N selected instants in time coincide with those instants that each of the particles is located at a specified position x.
That is, these are the N (x) particles located at the left face of the interval x to x + dx ( Figure B1). We may then examine how the number N (x) and the ensemble distribution f Ep (E p , x) and its moments change over the interval x to x + dx as this particle cohort moves downslope. Note that each member of the cohort not deposited within this interval may arrive at position x + dx at a different instant in time. This is unimportant, however, as we are interested only in how the energy states of the 10 particles vary with position x. Similarly, upon choosing any subsequent downslope position x, we must recognize that the N (x) particles reaching this position do so at entirely different instants in time. Here is a final note: In this problem a particle ensemble average is identical to a Gibbs ensemble average.

Appendix C: The Fokker-Planck-like equation
Let q = E p (x + dx) − E p (x) denote a change in the energy of a particle over the small distance dx, and let f q (q; E p , x) denote 15 the probability density function of the changes q associated with the energy state E p and position x. If n Ep (E p , x) denotes the number density of particle energies E p , then according to the Master equation, 20 Assuming the density f q (q; E p , x) is peaked near q = 0 with finite first and second moments, we may expand the integrand in Eq. (C1) as a Taylor series to second order, subtract n Ep (E p , x) from both sides, then divide by dx and take the limit as dx → 0 to obtain a Fokker-Planck-like equation, namely, Here, k 1 (E p , x) is a drift speed and k 2c (E p , x) is a diffusion coefficient defined by and The drift speed k 1 (E p , x) has two parts, one associated with gravitational heating and one associated with frictional cooling.
Starting with gravitational heating, let h(x) denote the height of a particle within the gravitational field at position x. If E p (x) denotes the particle kinetic energy equal in magnitude to the potential energy mgh(x) at height h(x), then E p (x + dx) = 5 mgh(x + dx) is the magnitude of the particle kinetic energy at the height h(x + dx), assuming a complete conversion of gravitational to kinetic energy without loss. Thus, This indicates that q in Eq. (C3) is independent of the energy state E p and therefore may be removed from the integral. We thus write Eq. (C3) as This is the steady rate of gravitational heating.
The part of k 1 (E p , x) associated with frictional cooling is obtained as follows. With particle-surface collisions we may 15 assume that q is proportional to the expected value of ∆E p . In turn we let n x = 1/λ denote the expected number of collisions per unit distance, where λ is the expected travel distance between collisions. This leads to where the overline denotes an average over particles at the energy state E p (rather than a global average).
Because gravitational heating is a fixed quantity according to Eq. (C6), heating does not involve diffusion. In turn, the 20 diffusion coefficient k 2c (E p , x) associated with frictional cooling is given by Note that whereas k 1h is a fixed quantity, k 1c and k 2c must be viewed as statistically expected quantities.

Appendix D: Expected travel distance between collisions
Momentarily let v = u 2 1/2 , and then let v 0 denote the surface-parallel velocity of a particle rebounding with reflection angle 25 φ measured from the surface ( Figure D1). We then know that where U 0 denotes the horizontal velocity and W 0 denotes the vertical velocity. For a vertical change in elevation Z over a horizontal distance X associated with the surface-parallel distance λ, we know that Z = −SX = −SU 0 t 0 , where t 0 is the travel time. For a rebounding particle starting at position z 0 = 0 we may deduce from Newton's second law that which gives With λ = X/ cos θ = U 0 t 0 / cos θ, we then combine Eq. (D1), Eq. (D2) and Eq. (D3) to obtain λ = 2v 2 0 sin(φ − θ) cos(φ − θ) g cos 2 φ cos θ + 2Sv 2 0 cos 2 (φ − θ) g cos 2 φ cos θ .
With reference to Figure E1, consider an idealized collision of a spherical particle with a rigid planar surface with slope Figure E1. Definition diagram for idealized collision of a spherical particle with a rigid planar surface.
angle θ. Let u, w and ω respectively denote the surface parallel velocity, the surface normal velocity and the angular velocity of the particle with mass m and radius r = D/2, and let the subscripts 1 and 2 denote incident and reflection values. With appropriate modification of the coordinate and sign convention used by Brach (1991), the momentum components associated with impulses can be expressed as (Brach, 1991;Brach and Dunn, 1995;Brach, 1998) 15 u 2 = u 1 + µ c (1 + )w 1 + gτ (sin θ − µ c cos θ) and (E2) where w 1 < 0, is the normal coefficient of restitution attributed to Newton, µ c is the ratio of tangential to normal impulses 20 during the collision, and τ is the impulse duration. Note that µ c generally is not a coefficient of friction, although it may be equal to a coefficient of friction in special cases, for example, with sliding throughout the entire duration of the collision (Brach, 1991;Brach and Dunn, 1992). Also note that µ c = 0 if u 1 = 0.
The second term on the right side of Eq. (E2) represents the effect of tangential friction on the velocity u, increasing with the magnitude of the normal impulse associated with the velocity w 1 . This term may be considered the dynamic contribution to friction during τ . The term gτ sin θ represents the downslope contribution to the impulse associated with the weight of the particle, and the term gτ µ c cos θ represents an enhancement of friction associated with this weight. The impulse duration τ may be on the order of milliseconds for a hard particle impacting a hard surface. It may be longer for a hard particle impacting 5 a relatively soft surface (Brach, 1991). If the magnitude of w 1 is sufficiently large and τ is sufficiently short, the gravitational terms in Eq. (E2) and Eq. (E3) may be neglected. The second term on the right side of Eq. (E3) represents the effect of tangential friction in contributing to rotational motion, that is, the conversion of translational energy to rotational energy.
Collisions involving small incident angles begin with sliding during the impulse duration τ . If with a sufficient normal dynamic force this initial sliding gives way to stick prior to separation, then for a sphere with moment of inertia I = (2/5)mr 2 , 10 the velocity u 2 = rω 2 at separation. This leads to which represents the outcome of a conversion of translational to rotational motion with stick. Whereas the resultant velocity u 2 can be determined in this situation, the effect of sliding on u 2 cannot be analytically constrained. Nonetheless, Eq. (E3) indicates that collisions induce a conversion of translational to rotational motion in that tangential friction during an impulse 15 exerts a torque on the particle, thereby extracting translational kinetic energy that is in addition to work performed during particle deformation and by friction. We also note that low-angle collisions likely dominate in the problem at hand.
In order to recast the problem in terms of kinetic energy, we start by squaring Eq. (E1), Eq. (E2) and Eq. (E3) to give We may immediately neglect terms involving τ 2 , and for sufficiently large w 1 and small τ we may neglect terms involving τ .
The next task involves scaling the normal velocity w 1 in terms of the tangential velocity u 1 in relation to particle motions down an inclined surface. Hereafter we focus on lowest order effects. With reference to the analysis presented in Appendix 5 D, let W 0 denote the vertical reflection velocity of a particle following a collision. Assuming downslope motion, then for any finite horizontal reflection velocity U 0 and reflection angle φ, the magnitude of the vertical velocity at the next collision is given by where Z ≤ 0 is the vertical distance between the collisions. That is, From Appendix D, W 0 = v 0 sin(φ−θ)/ cos φ, Z = − 1 2 gt 2 0 +W 0 t 0 and t 0 = 2W 0 /g +2SW 0 /g, where v 0 is the surface parallel velocity associated with W 0 , t 0 is the travel time and S = tan θ. Using these relations with Eq. (E10) we obtain Expanding the trigonometric functions in Eq. (E11) as Taylor series and retaining the lowest order term in φ we obtain In effect the magnitude of w 1 is set by the gain in the magnitude of the vertical velocity associated with conversion of gravitational potential energy to translational energy with finite slope. This strengthens the normal impulse of the particle, but only 20 up to a slope (nominally 45 degrees) beyond which the surface normal component of the vertical velocity begins to decrease.

25
Omitting subscripts, this is Comparing this result with the assumption ∆E p = −β x E p , we may conclude that Note that Eq. (E15) and Eq. (E16) pertain to a highly idealized collision. In fact, the quantities µ c , and φ are each random variables. Moreover, on an irregular hillslope surface the angle θ also is a random variable when viewed at the particle-surface collision scale. Nonetheless, for the purpose of scaling we may view this angle as a locally averaged value, and we now take 5 the ensemble average of Eq. (E16) to give In turn, with µ = β x /4 tan φ we may write with 10 M (θ) = 2µ c (1 + ) 4 sin φ cos θ sin θ .
At lowest order, cos θ sin θ ∼ θ. We therefore may expect µ to systematically vary with the slope angle θ. Also note that µ is independent of particle size. We examine both of these points in the second companion paper (Furbish et al., 2020a).
Subtracting u 2 1 from both sides of Eq. (8), multiplying by m/2 and retaining the lowest order term, 1 2 m(u 2 2 − u 2 1 ) ≈ − 1 − 25 49 This is This result indicates that the onset of rotation with stick produces a large change in the slope-parallel kinetic energy. In this case, β x ≈ 0.5. Again notice that this result is independent of particle size. Nonetheless, the numerical factors in Eq. (4), Eq. (8) and Eq. (E20) are set by the moment of inertia of the particle, which means that these factors vary with irregular particles. Also 20 note that Eq. (E21) does not imply that half of the translational energy E p is converted entirely to rotational energy. Rather, half is converted to rotational energy and lost to work performed by friction prior to stick and by particle/surface deformation, thence dissipated as heat, vibrations and sound.
Subtracting ω 2 1 from both sides of Eq. (E22) and multiplying by I/2 = (1/5)mr 2 , 1 2 I(ω 2 2 − ω 2 1 ) = −mµ c r(1 + ) cos θ sin θ cos φ ω 1 u 1 which is This result suggests that in the absence of initial rotation (E r = 0), a change in rotational energy is directly related to the translational energy E p , where the proportion β x now represents the leading factors in the first term on the right side of Eq. 10 (E24). With extant rotational motion, a weaker conversion of translational to rotational energy occurs according to the second term on the right side. Both cases are slope dependent due to the connection between w 1 and u 1 ∼ v 0 implied by Eq. (E12).
Focusing on downslope motions, in general we may write the energy balance of a particle as Here, a positive change in rotational energy ∆E r is seen as an extraction of translational energy. Then, for example, this loss is 15 given explicitly by Eq. (E20) in the specific case of stick with the onset of rotation. An approximation of this loss is given by Eq.
(24) for a frictional collision that does not necessarily involve stick. The term f c in Eq. (E25) represents losses associated with particle and surface deformation as well as work performed against friction during collision impulses (converted to heat, sound, etc.). This is represented, for example, by Eq. (15). But this term also includes losses associated with deformation of the surface at a scale larger than that of an idealized particle-surface impulse contact, namely, due to momentum exchanges associated with glancing collisions that produce transverse translational motions and rotational motions oriented differently than that considered above ( Figure E1). In some cases, as described above, the change in energy ∆E p can be expressed directly in 25 terms of the energy state E p . However, the complexity of particle-surface collisions on natural hillslopes precludes explicitly demonstrating such a relation for all possible scenarios. Nonetheless, the examples above suggest that it is entirely defensible to assume that energy losses can be related to the energy state E p if the elements involved are formally viewed as random variables. Specifically, with the effect of slope angle θ on the impact velocity w 1 and its relation to u 1 ∼ v 0 via Eq. (E12), we can be confident that the loss ∆E p is functionally related to the energy state E p . The simple relation ∆E p = −β x E p thus is to 30 be viewed as an hypothesis to be tested against data, as elaborated in the second companion paper (Furbish et al., 2020a).
writing this inequality as This means that Whereas the left side of Eq. (F2) is the arithmetic average E a , the right side is the harmonic average E h . Thus, Because the arithmetic average of a set of positive numbers is always greater than or equal to the harmonic average of this set, this inequality is indeed satisfied. These averages are equal only if all values of the set are equal, that is, the variance of the set 10 is zero.
We do not know the form of the underlying distribution f Ep (E p , x). For physical reasons, however, it cannot be a distribution that supports E p → 0, as this coincides with particles at rest. For example, f Ep (E p , x) cannot be an exponential or Weibull distribution with support E p ∈ [0, ∞). In contrast, the lognormal and gamma distributions with support E p ∈ (0, ∞) are admissible, and the Pareto distribution with support E p ∈ [E pm , ∞) is admissible. 15 As a point of reference, for a density f Ep (E p ) with finite expected value E p , the density f y (y) of the reciprocal y = 1/E p may not have defined moments. This occurs, for example, if f Ep (E p ) is exponential with support E p ∈ [0, ∞). Interestingly, if with x = ln(E p ) the density f x (x) is lognormal with mean µ, then with y = 1/x the density f y (y) also is lognormal with mean −µ.
For a density f y (y) of y = 1/E p with undefined mean, the average y calculated from a sample nonetheless is finite, as the 20 probability of sampling precisely a value E p = 0 is identically zero. Moreover, as the variance of E p becomes small for finite mean E p , the product E p y = E p 1/E p → 1, as in the discrete case above.

Appendix G: Deposition rate
Our description of the deposition rate for a granular gas in a box has both similarities and dissimilarities with the processes of deposition (de-sublimation) and condensation. Here we briefly outline key points.

25
In a closed system involving two phases (solid/gas or solid/liquid) at thermodynamic equilibrium, the rates of deposition and sublimation (or condensation and evaporation) are equal. That is, the rate at which molecules move from the solid phase to the gas phase (or from the liquid phase to the gas phase) is balanced by the rate at which molecules move from the gas phase to the solid phase (or from the gas phase to the liquid phase). These rates, in each direction, depend only on the thermal state of the system. Because the system has specified internal energy involving conservative particle-particle collisions, we do not need to appeal to the idea of heating and cooling (although this could be occurring). For a granular gas involving dissipative collisions, however, a non-equilibrium steady state is achieved only if it is continuously mechanically heated, and the rate of heating is matched by the rate of cooling due to the collisions. (Note that we refer to a non-equilibrium steady state condition rather than thermal equilibrium, as unlike an ordinary gas, a granular gas can exhibit strong spatial correlations in the particle 5 number density (see Pöschel (2004, 2005) and Brilliantov et al. (2018) and references therein; and van Zon and MacKintosh, (2004)). However, this distinction is unimportant in relation to the behavior of particle motions on a hillslope envisioned as a rarefied granular gas.) Like an ordinary solid-gas system, the rate of sublimation (entrainment) is matched by the rate of deposition (disentrainment), and the total particle energy and the average particle energy are fixed. Moreover, like an ordinary solid-gas system, the deposition rate depends on the physics of disentrainment in relation to its thermal state, not 10 on the difference between the heating and cooling rates (which is zero at steady state). Heating modulates the deposition rate as described in the text.

Appendix H: Generalized Pareto distribution
Solving Eq. (65) giveŝ 15 Using Eq. (64) the disentrainment rate is then which we write aŝ Making use of Eq. (4) we then obtain the distribution of travel distances, namely, This is a generalized Pareto distribution with location parameter equal to zero.
To show how the generalized Pareto distribution is related to the ordinary Lomax distribution, we start by rewriting Eq. (H4) as

25
This is We now define the shape parameter a L = 1/a and the scale parameter b L = b/a. This gives a Lomax ditribution, namely, Thus, for a > 0 the behavior of the generalized distribution, Eq. (H4), is the same as that of a Lomax distribution. The mean is We work with the generalized Pareto distribution in the form of Eq. (H4) because of the clear connection between its parameters 5 and the disentrainment rate function, Eq. (H3), and because the condition a < 0 is physically meaningful.

Appendix I: Kirkby-Statham formulation
The formulation of Kirkby and Statham (1975) assumes that initial particle kinetic energy is dissipated in work performed by a fixed Coulomb-like friction to give an average travel distance. This idea can be formulated in terms of momentum and energy, then recast in terms of the rate of change in energy with respect to position x for comparison with the formulation presented in 10 the main text.
In appealing to a Coulomb-like friction behavior, Kirkby and Statham (1975) start with F x = mg sin θ − µ d mg cos θ. With particle velocity u we write this as du(t) dt = g sin θ − µ d g cos θ .
Note that u(t) must be envisioned as representing an idealized "average" velocity of a group of particles viewed over time. 15 This gives u(t) = (g sin θ − µ d g cos θ)t + u 0 .
For a total travel time T , u p (T ) = 0 = (g sin θ − µ d g cos θ)T + u 0 , so that 20 T = − u 0 g sin θ − µ d g cos θ . (I4) In turn we rewrite Eq. (I2) as dx(t) dt = (g sin θ − µ d g cos θ)t + u 0 , so that x(t) = 1 2 (g sin θ − µ d g cos θ)t 2 + u 0 t . The total travel distance X is thus Using the initial squared velocity u 2 0 = 2 gh sin 2 θ, This is the result that Kirkby and Statham (1975) offer as representing the average travel distance.

5
We now turn to kinetic energy. Let A = g sin θ − µ d g cos θ. Multiplying Eq. (I1) by mu then leads to d dt With u = At + u 0 from Eq. (I3), This leads to 10 E p (t) = 1 2 mA 2 t 2 + mAu 0 t + E p0 .
Substituting this into Eq. (I11) and doing algebra then yields E p (x) = mAx + E 0 . The derivative of this result with respect to x is 15 dE p (x) dx = mg sin θ − µ d mg cos θ .
This result is like Eq. (57), but absent the effect of deposition and the associated apparent heating, as it strictly applies to the motion of an individual particle or a group of particles acting like a rigid body. It does not describe an ensemble averaged motion.
Appendix J: Gabet-Mendoza formulation 20 Gabet and Mendoza (2012) appeal to ideas from Samson et al. (1998) and Quartier et al. (2000) and suggest that the motion of an individual particle can be described as du(t) dt = g sin θ − µ d g cos θ − κu ψ .
However, whereas the derivative term on the left side of Eq. (J1) and the first two terms on the right side pertain to the instantaneous motion of an individual sliding particle or group of particles acting like a rigid body, the third term on the right side, representing collisional friction, actually is relevant to time-averaged or ensemble-averaged behavior rather than the instantaneous behavior of an individual particle (Riguidel et al., 1994a(Riguidel et al., , 1994bSamson et al., 1998Samson et al., , 1999. These terms are not additive as written. The gravity and Coulomb friction terms are like those in the formulation of Kirkby and Statham (1975).
Because there is confusion in the literature regarding the collisional friction term, here we elaborate its form.
Let n t denote the expected number of particle-surface interactions (collisions) per unit time as a particle moves downslope, 5 and let β x denote the proportion of momentum parallel to x that is extracted during an individual collision involving the particle velocity u. Recognizing that both β x and u must be treated as random variables, and letting angle brackets denote an ensemble average, we may now assume that The first term on the right side of Eq. (J2) represents the uniform gravitational force, and the second term on the right side 10 represents a frictional force due to particle-surface collisions (compare with Eq. (2) in Riguidel et al., 1994). As a reminder, this term is entirely analogous to the dissipation term that Haff (1983) introduced (formulated in terms of energy rather than momentum), leading to Haff's cooling law (Brilliantov and Pöschel, 2004;Yu et al., 2020). The proportion of momentum extracted, β x , involves an appropriate coefficient of restitution depending on the geometrical details of the collision. We may now assume that n t ∼ u /l, where l denotes a characteristic length scale representing the expected distance between collisions. 15 This leads to which is close to the form of Eq. (J1) with ψ = 2 (neglecting the Coulomb friction term), but not quite.
We now focus on uniform, steady conditions such that u is unchanging with position or time, consistent with various experiments (Riguidel et al., 1994;Samson et al., 1998). This leads to 20 β x u u ≈ lg sin θ .
We now write β x = β x + β x and u = u + u , where primes denote deviations about the expected values. Substituting these expressions into Eq. (J4) and taking expected values then leads to β x u 2 + β x u u ≈ lg sin θ .
The product β x u 2 has the appearance of a nominal, nonlinear viscous term. Samson et al. (1998) suggest that this represents 25 a Bagnold-like friction based on analogy with the scaling provided by Bagnold (1954), preceding the critical assessment of Bagnold's experimental work presented by Hunt et al. (2002). The term β x u u , neglected at the outset by Riguidel et al. (1994), looks like a linear viscous term, where the "viscosity" is given by the covariance β x u .
Noting that Eq. (J5) is quadratic, we can solve for the velocity u and determine that at lowest order u ≈ lg sin θ β x so long as ( β x u / β x ) 2 < 4lg sin θ/ β x . If this inequality is satisfied, then Eq. (J3) becomes giving ψ = 2. Note that the squared average velocity in Eq. (J7) does not imply that collisional friction is scaled with kinetic energy rather than momentum. This result occurs because n t is initially scaled with u and l. Quartier et al. (2000) present an analogous formulation; see their Eq. (4) and explanation of the squared velocity term. Dippel et al. (1997) also discuss this 5 point.
In relation to their experiments involving particles of radius R moving down an inclined surface roughened with a quasirandom monolayer of particles with radius r m , Riguidel et al. (1994) and Samson et al. (1998) propose the hypothesis that u ∼ sin θ. This derives from a scaling analysis in which the magnitude of the collisional momentum extraction (i.e., β x ) is written as a function of the relative smoothness R/r m by introducing an unconstrained velocity quantity. These authors plot 10 measured values of u versus sin θ and suggest that the linear fit confirms a viscous-like behavior. Note, however, that because of the rather limited experimental range of sin θ (Figure 2 in Riguide et al., 1994; Figure 2 in Samson et al., 1998; Figure 4 in Samson et al., 1999), the data are equally well fit by a straight line in a plot involving √ sin θ ( Figure J1), consistent with the Figure J1. Plot of ensemble averaged particle velocity u versus √ sin θ involving a steel sphere (R = 2.5 mm) moving over glass beads (rm = 0.53 mm) giving R/rm = 4.7; data from Samson et al. (1998). collisional-based formulation, Eq. (J7). (Sampson et al. (1999) acknowledge this limitation of the range of sin θ.) In addition, we can scale the length l as l ∼ r m /c A , where c A = 0.67 is the areal concentration of the surface-roughness particles. For a 15 fixed velocity u , Eq. (J7) gives l ≈ β x u 2 /g sin θ. With R/r m = 4.7 and a coefficient of restitution of ≈ 0.8, we can estimate β x ≈ 0.05. This gives l ≈ 0.5 − 0.6 mm over the range of measured velocities in Figure J1, which is close to the experimental value of l = r m /c A = 0.8 mm, thereby reinforcing the collisional basis of Eq. (J7). Thus, a spherical particle that macroscopically rolls over a monolayer roughness is actually going bumpety-bump, colliding with monolayer particles during its motion. 20 https://doi.org/10.5194/esurf-2020-98 Preprint. Discussion started: 8 December 2020 c Author(s) 2020. CC BY 4.0 License.
Because the relative smoothness R/r m is not entirely adequate in scaling the collisional friction as c A varies (Samson et al., 1998), it remains unclear whether these experimental conditions involve an apparent viscous-like behavior where the effective viscosity depends only on roughness geometry (Dippel et al., 1997;Samson et al., 1999) or a squared-velocity behavior as in Eq. (J7). Dippel et al. (1997) note that there is an apparent cross-over in behavior for very large and heavy spheres. Effects of the covariance of β x and u in relation to roughness geometry and the details of motions, including transverse motions, likely are 5 important. Nonetheless, we emphasize that in the formulations of Riguidel et al. (1994) and Samson et al. (1998), a Coulomb friction behavior is not involved.
Returning to Eq. (J1), similarly there is no clear reason to include a Coulomb-like friction term, as natural irregular particles mostly do not slide down natural rough surfaces. In addition, if the starting point involves the derivative term on the left side and the gravitational term on the right side as written, then the collisional term on the right side should be a random quantity, 10 thus leading to a stochastic differential equation -that is, a Langevin-like equation (Riguidel et al., 1994) -not an ordinary differential equation. Moreover, the idea of a dynamic friction coefficient is misapplied in the situation where rarefied particles tumble, roll and skitter over the surface. A Coulomb model is appropriate for sustained contact, and even then a dynamic friction involves collisional friction at the surface asperity scale. Particle-surface contacts on natural granular surfaces are not smooth at a scale commensurate with a sliding Coulomb model. A rolling coefficient of friction works for spheres moving over 15 a relatively smooth surface, not for irregular tumbling particles involving non-collinear impacts. Moreover, the static normal weight of a particle, mg cos θ, does not set the particle-surface friction. Rather, dynamic forces during collision impulses matter (Brach, 1991;Stronge, 2000). This includes the dynamic Coulomb friction force associated with conversion of translational to rotational kinetic energy during collisions (Appendix E). Any resulting dynamic friction coefficient represents an ensemble averaged ratio of tangential to normal momentum exchanges, both of which are random variables. (This point currently is being 20 examined in studies of bed load and aeolian transport; see for example Pähtz and Duran (2018).) Finally, the experiments of Quartier et al. (2000) involved rolling a cylinder over an inclined row of cylinders in an experiment designed to remove the transverse degree of freedom of motion. The Coulomb-like term in their formulation (see their Eq. (5)) arises from trapping of the rolling cylinder between bumps, and is unrelated to sliding as in a conventional Coulomb model.
Author contributions. All authors contributed to the conceptualization of the problem and its technical elements. DJF wrote much of the paper with contributions by the other authors. Because AMA is deceased, we offer the following clarification. AMA made key observations regarding particle motions in initial laboratory experiments during her MS work, and worked closely with DJF in conceptualizing a description of these motions. Whereas we did not succeed in formalizing AMA's ideas in mathematical terms while she was pursuing her work, these ideas were essential to a proper description of particle motions and are now formalized in the paper. We also note that AMA's co-authorship carries full approval of her family.