Rarefied particle motions on hillslopes: 2. Analysis

We examine a theoretical formulation of the probabilistic physics of rarefied particle motions and deposition on rough hillslope surfaces using measurements of particle travel distances obtained from laboratory and field-based experiments, supplemented with high-speed imaging and audio recordings that highlight effects of particle-surface collisions. The formulation, presented in a companion paper (Furbish et al., 2020a), is based on a description of the kinetic 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 5 of the particles. Both laboratory and field-based measurements are consistent with a generalized Pareto distribution of travel distances and predicted variations in behavior associated with the balance between gravitational heating due to conversion of potential to kinetic energy and frictional cooling due to particle-surface collisions. For a given particle size and shape these behaviors vary from a bounded distribution representing rapid thermal collapse with small slopes or large surface roughness, to an exponential distribution representing approximately isothermal conditions, to a heavy-tailed distribution representing net 10 heating of particles with large slopes. The transition to a heavy-tailed distribution likely involves an increasing conversion of translational to rotational kinetic energy leading to larger travel distances with decreasing effectiveness of collisional friction. This energy conversion is strongly influenced by particle shape, although the analysis points to the need for further clarity concerning how particle size and shape in concert with surface roughness influence the extraction of particle energy and the likelihood of deposition. 15


Introduction
As described in our first companion paper (Furbish et al., 2020a), we are focused on rarefied motions of particles which, once entrained, travel downslope over the land surface. This notably includes the dry ravel of particles down rough hillslopes following disturbances (Roering and Gerber, 2005;Doane, 2018;Doane et al., 2019;Roth et al., 2020) or upon their release from obstacles (e.g., vegetation) following failure of the obstacles (Lamb et al., 2011DiBiase and Lamb, 2013;DiBiase et al., 20 2017; Doane et al., 2018Doane et al., , 2019, and the motions of rock fall material over the rough surfaces of talus and scree slopes (Gerber and Scheidegger, 1974;Kirkby and Statham, 1975;Statham, 1976;Tesson et al., 2020). By "rarefied motions" we are referring to the situation in which moving particles may frequently interact with the surface, but rarely interact with each other. Thus, rarefied particle motions are distinct from granular flows. Although this idea is most applicable to processes such as rock fall and the subsequent motions of the rock material over talus or scree slopes, our description of the motions of individual particles nonetheless may be entirely relevant to conditions that are not strictly rarefied, but where during the collective motions of many particles the effects of particle-surface interactions dominate over effects of particle-particle interactions in determining 5 the behavior of the particles -akin to granular shear flows at high Knudsen number (Kumaran, 2005(Kumaran, , 2006. We note that laboratory experiments (Kirkby and Statham, 1975;Gabet and Mendoza, 2012) and field-based experiments (DiBiase et al., 2017;Roth et al., 2020) designed to mimic particle motions and travel distances on hillslopes effectively focus on rarefied conditions. The formulation of rarefied particle motions presented in the first companion paper (Furbish et al., 2020a) is based on a 10 description of the kinetic 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 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 as this distribution varies downslope. The Kirkby number Ki -the ratio 15 of gravitational heating to frictional cooling -sets the basic deposition behavior and the form of the probability distribution f r (r) of particle travel distances r. For isothermal conditions where frictional cooling matches gravitational heating plus the apparent heating due to deposition, the distribution f r (r) is exponential. With non-isothermal conditions and small Ki this distribution is bounded and represents rapid thermal collapse. With increasing Ki the distribution f r (r) takes the form of a heavy-tailed Pareto distribution. It may possess a finite mean and finite variance with moderate Ki , or the mean and variance 20 may be undefined with large Ki .
The purpose of this second companion paper is to present an analysis of several data sets concerning particle motions on rough surfaces, as viewed through the lens of the theory presented in the first companion paper (Furbish et al., 2020a). In Section 2 we summarize the context for our work provided by recent probabilistic descriptions of the flux and the Exner equation (Furbish and Haff, 2010;Furbish and Roering, 2013), and then step through essential elements of the mechanical 25 basis of the theory leading to the generalized Pareto distribution of particle travel distances. In Section 3 we compare the formulation with the laboratory measurements of particle travel distances on rough surfaces reported by Gabet and Mendoza (2012) and Kirkby and Statham (1975). We also report laboratory experiments designed to clarify how the size and shape of particles influence their motions and disentrainment based on high-speed imaging. In Section 4 we compare the formulation with the field-based measurements of travel distances reported by DiBiase et al. (2017) and Roth et al. (2020). 30 Particle travel distances from both the laboratory and field-based experiments are consistent with the generalized Pareto distribution and provide compelling evidence for the full range of predicted behaviors, from rapid thermal collapse to approximately isothermal conditions to net heating of particles. Nonetheless, the analysis points to the need for further clarity concerning 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 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 5 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.
2 Key elements of theoretical formulation

Probabilistic description of disentrainment
The problem of describing rarefied particle motions on hillslopes is motivated by the entrainment forms of the flux and the 10 Exner equation. Namely, let f r (r; x) denote the probability density function of the travel distances r of particles whose motions start at x, and let R r (r; x) denote the associated exceedance probability function. Assuming motions are only in the positive x direction and noting that x = x − r, the flux q(x) may be written as where E s (x) denotes the volumetric entrainment rate at position x. In turn, letting ζ(x, t) denote the local land-surface eleva- 15 tion, the entrainment form of the Exner equation is where c b = 1 − φ s is the volumetric concentration of the surface with porosity φ s . The central elements of Eq. (1) and Eq. (2) are the probability density function f r (r; x) and the associated exceedance probability function R r (r; x). These are related to the disentrainment rate function defined as 20 P r (r; x) = f r (r; x) R r (r; x) , which, when multiplied by dr, 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, Eq, (3), connects the descriptions of the flux and its divergence, Eq. (1) and Eq.
(2), to the physics of particle motions and disentrainment.
For completeness we note that the formulation above involving continuous functions can be recast into a discrete form that is 25 useful for considering situations in which conditions influencing particle motions, for example the surface slope and roughness texture, change in the downslope direction. Let k = 1, 2, 3, ... denote a set of discrete intervals of length dr. Let p k denote the probability that a particle, having not been disentrained before the kth interval, then becomes disentrained within this interval.
The probability mass function of particle positions is then The probability p k , like its continuous counterpart P r (r; x), connects the descriptions of the flux and its divergence to the physics of particle motions and disentrainment. This discrete formulation opens the possibility of describing effects of varying disentrainment rates in response to changing downslope conditions in a manner intrinsic to particle-based treatments of 5 transport (Tucker and Bradley, 2010), but not readily incorporated in probabilistic descriptions. That is, if surface conditions change in the downslope direction, for example, giving net cooling followed by heating or vice versa, then particles whose travel distances are large enough "see" this change and their behavior concomitantly changes.
As summarized next, the analysis presented in Furbish et al. (2020a) describes the mechanical basis for the disentrainment rates P r (r; x) and p k , and the associated probability distributions f r (r; x) and f k (k). This involves a consideration of the 10 kinetic energy balance of rarefied particle motions and how this balance determines the deposition of particles in relation to their energy state.

Energy and mass balances
Consider a rough, inclined surface with uniform slope angle θ (Figure 1). At this juncture we simplify the notation and consider  (Furbish et al., 2020a).
the motions of particles entrained at a single position x = 0. Now the particle travel distance r → x and the probability density 15 function f r (r; x) → f x (x). 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 deposited within it. We now imagine collecting this great number of particles and treat them as a cohort, independent of time (Furbish et al., 2020a, Appendix B). 20 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. The number of particles deposited within the control volume is dN = N (x + dx) − N (x). The objective is then to determine the rate of particle deposition, dN/dx, based on the energy state of the cohort of particles.
In particular, let E p = (m/2)u 2 denote the translational kinetic energy of a particle with mass m and downslope velocity u. Letting angle brackets denote an ensemble average of a great number N of particles, then we denote the arithmetic average energy as E a = E p and the harmonic average energy as E h = 1/ 1/E p . The total energy E = N E a . The formulation presented in the first companion paper (Furbish et al., 2020a) then leads to four equations with four unknowns involving energy and mass.
Namely, conservation of the total energy of the particle cohort is given by The first term on the right side of Eq. (5) represents gravitational heating with the uniform conversion of potential to kinetic energy, the second term on the right side represents frictional cooling due to particle-surface collisions, and the last term on the right side represents a loss of energy associated with particle deposition. The friction factor µ is where β x is the expected proportion of translational energy extracted during a particle-surface collision, and φ is the expected reflection angle of particles following collision. The factor α modulates the particle deposition length scale L c , which represents an e-folding distance over which deposition occurs. This length scale is given by 15 and is a function α = f (Ki ) of the dimensionless Kirkby number Ki defined by which represents the ratio of gravitational heating to frictional cooling.
Conservation of particle mass is given by 20 which illustrates that deposition is proportional to frictional cooling depending on the particle energy state E h , modulated by the factor α.
The ensemble averaged energy satisfies the expression, where the arithmetic and harmonic averages are related as As with the total energy described by Eq. (5), the first term on the right side of Eq. (10) represents gravitational heating and the second term on the right side represents frictional cooling. Because the inequality in Eq. (11) must be satisfied, the last term in Eq. (10) represents an apparent heating due to particle deposition whose effect is to cull lower energy particles, thereby selecting higher energy particles for continued downslope motion. Brilliantov et al. (2018) describe an analogous behavior of granular gases due to particle aggregation. Simultaneous solution of Eq. (5), Eq. (9) and Eq. (10) using Eq. (11) leads to the disentrainment rate function, In turn the probability density function f x (x) of travel distances x for particles starting at x = 0 is and the exceedance probability function is The shape and scale parameters A and B are 15 The mean of the distribution is With reference to the presentation in the first companion paper (Furbish et al., 2020a) and for the purpose of presenting results below, 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 θ. For plotting purposes, and to highlight the role 20 of the Kirkby number, we now define the following dimensionless quantities denoted by circumflexes: With these definitions in place, Eq. (5), Eq. (9), Eq. (10) and Eq. (11) are rewritten as The expressions involving energy, Eq. (19) and Eq. (21), reveal that the Kirkby number Ki has a key role in the energy balance and therefore particle deposition. In addition we can define a transitional Kirkby number, If Ki < Ki * then net cooling occurs, and if Ki > Ki * then net heating occurs. The condition Ki = Ki * implies isothermal conditions such that dÊ a /dx = 0.
The dimensionless disentrainment rate is 15 The dimensionless probability density function of travel distances is and the exceedance probability function is The shape and scale parameters a and b are 20 The distribution fx(x) defined by Eq. (25) is a generalized Pareto distribution with position parameter equal to zero (Pickands, 1975;Hosking and Wallis, 1987). To summarize with reference to Figure 2, for a < 0 this distribution decays more rapidly than an exponential distribution and is bounded at the positionx = b/|a|. For a = 0 it becomes an exponential distribution. For 0 < a < 1/2 the distribution fx(x) 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 are summarized in Table 1. If the shape and scale parameters a and b are redefined as a L = 1/a and b L = b/a = ba L , then Eq.

5
(25) becomes   (Furbish et al., 2020a). Compare with Figure 1 in Hosking and Wallis (1987). Table 1. Behavior of the generalized Pareto distribution associated with the shape parameter a and Kirkby number Ki as illustrated in Figure   2.
Finite mean and variance Undefined mean and variance a ≥ 1 which is a Lomax distribution with mean With a < 0 the bounded form of fx(x) (Figure 2) represents rapid thermal collapse with net frictional cooling. For a = 0 10 the exponential form of this distribution represents isothermal conditions where frictional cooling is matched by gravitational For reference to data fitting presented below, a binomial expansion of Eq. (13) gives 15 The expansion of an exponential distribution with mean µ x is These two expansions indicate that unless the travel distance data span a significant proportion of the x domain, then at lowest order the fit of the generalized Pareto distribution looks like an exponential distribution with mean B. This result also is obtained from the disentrainment rate function, Eq. (12), in which for small x this function becomes P Moreover, if the travel distance data sample over a majority of the probability contained in the distribution, and if the tail of 5 the distribution is not "too" heavy, then B is an approximation of the mean (where A < 1 such that the mean exists).
Also note that in fitting the data to the generalized Pareto distribution, Eq. (13), we use the dimensional form of the exceedance probability, Eq. (14). Specifically, we estimate values of the exceedance probability as R where r x is the ascending rank of each datum. We then visually fit Eq. (14) to these estimated values to obtain values of the parameters A and B. This involves iteratively examining the data and theoretical lines in arithmetic, semi-log and log-log 10 plots, noting that semi-log plots generally highlight deviations in the tails whereas log-log plots highlight deviations near the origin. We mostly pay attention to the main body of data in the semi-log plots, avoiding over-fitting of the outer part of the tails given the inherent variability of estimates of the exceedance probability associated with the tails, notably with small sample sizes (Appendix A). We also examine the data using quartile-quartile (QQ) plots to ensure that these are consistent with the generalized Pareto distribution. Here we note that our objective is to demonstrate that the data are consistent with the several 15 forms of the generalized Pareto distribution, where in semi-log space the exceedance plots either: 1) have negative concavity (representing rapid thermal collapse); 2) are approximately straight (representing isothermal conditions); or 3) have positive concavity (representing net heating). We are aimed at reasonable estimates of the shape and scale parameters in order to achieve this objective, but we do not need refined estimates of these quantities. For reasons that are fully explained in Appendix A, we therefore use estimates of A and B obtained from visual fitting, avoiding known limitations of quantitative estimates (e.g., 20 maximum likelihood estimates) associated with small sample sizes, particularly in the presence of censorship. Assuming a value of γ (see description below), we then have two equations, Eq. (15) and Eq. (16), with two unknowns, µ and α. Thus, the fitting of A and B provides estimates of µ and α for subsequent consideration. In particular, we first estimate µ as and then obtain α as In turn the Kirkby number Ki is calculated from Eq. (8).

Elements of collisional friction
The quantities µ and α defined in relation to Eq. (6) and Eq. (7) merit further description. We start by noting that the formulation summarized in the preceding section is based on the assumption that a change in translational kinetic energy ∆E p associated 10 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 of the first companion paper (Furbish et al., 2020a), 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 15 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. (34) represents losses associated with particle and surface deformation as well as work performed against friction during collision impulses (thence converted to heat and sound). 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 in Furbish et al., 2020a). However, the complexity of particle-surface collisions 25 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.
This hypothesis formally enters the formulation via Eq. (6). Namely, from this relation we may write µ ∼ β x , highlighting that µ is associated with the cooling rate. In turn, particle collision mechanics (Appendix E in Furbish et al., 2020a) suggest, 30 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. These points are examined below in relation to experimental data.

5
In turn, Furbish et al. (2020a) suggest that the quantity α is related to the Kirkby number Ki as where α 0 denotes the value of α associated with a flat surface (Ki = 0) and µ 1 is a factor of order unity. Substitution into Eq.
(7) gives 10 Viewed in this manner, α represents a direct effect of heating described by mgµ 1 sin θ, 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 relative to the cooling length scale l c = E h /mgµ cos θ. In this view, µ goes with the cooling rate (not the deposition rate). But we also may write Eq. (36) Viewed in this manner, we may define an apparent friction factor as µ 0 = µ(1 − µ 1 Ki ) associated with deposition. Here again, µ is associated with the cooling rate but is then modulated by heating. We suggest below 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 -both decreasing the likelihood of stopping and not represented in just the factor µ. We also 20 suggest that for the same slope and surface roughness, α increases with increasing particle size, decreasing the likelihood of frictional loss with increasing rotational energy.  The experiments reported by Gabet and Mendoza (2012) involved launching spherical, sub-angular one-cm particles with initial velocity u 0 = 0.7 m s −1 onto an inclined surface with fixed roughness elements embedded within concrete (see Figure 2 therein). The experiments stepped through slope angles of θ = 0 degrees to θ = 30 degrees in increments of three degrees. The travel distances of 100 particles were measured for each slope angle. All 100 particles stopped within the 3 m flume for angles 0 to 15 degrees. For angles 18, 21, 24 and 27 degrees, 92, 58, 26 and 4 particles stopped, respectively. No particles stopped on the 30 degree slope.
Because the same particle is launched each time with the same initial velocity u 0 , the initial arithmetic and harmonic averages of the particle energy are the same, that is, E a0 /E h0 = γ 0 ≈ 1. Over some unknown distance the particle motions experience randomization via collisions such that E a /E h = γ > 1. With reference to Eq. (9) where E h = E a /γ, this means that because 5 γ is initially one and then increases, the expected disentrainment rate likely increases initially over a short distance. Indeed, in preliminary plots (not shown) of the exceedance probability R x (x) versus x, an inflection occurs in some of the data close to x = 0, which we suspect represents a delay in the onset of deposition of the lowest energy particles. In fitting the data to the exceedance probability R x (x) we assume that E a /E h = γ > 1 (and fixed) over the entirety of the travel distances. For this reason we truncate distances shorter than the inflection position then recalculate the travel distances and the exceedance 10 probabilities, while avoiding large truncation given the limited data size. We cannot know for certain the appropriate truncation position, so this point should be kept in mind. Note, however, that the formulation does not care where motions start, so this truncation is just a resetting of the staring position (x = 0) with fewer data, assuming γ remains approximately fixed beyond the adjusted starting position. These points are examine further in Appendix B with reference to experiments described in Section 3.3. 15 We choose γ = 1.5 in this and subsequent analyses of travel distance data. Note first that this choice has no influence on the estimates of the parameters A and B. However, it does influence the estimated values of µ and α via Eq. (32) and Eq. (33).
The assumption that γ is fixed may be incorrect. However, to rigorously constrain γ would require solving the Fokker-Planck equation describing the evolving number density n Ep (E p , x) of particle energy states E p as described in Section 3.3.2 of the first companion paper (Furbish et al., 2020a); this is beyond the scope of our objective of demonstrating the basic behavior 20 of the particle energy balance. The effect of fixed γ = 1.5 on estimates of µ and α is systematic throughout the analyses, but whether this choice is correct remains an open question.

Results
The data are reasonably well fit by the theoretical curves, where we plot the data twice ( Figure 3 and Figure 4) in order to highlight several points. Estimated parametric values are provided in Table 2, where estimates of the variability in µ, α, Ki and 25 Ki * are based on a Monte Carlo analysis as described in Appendix C.  Table 2. Fitted and estimated values of the parameters for the data shown in Figure 3 and Figure 4 with coefficients of variation in parentheses. b Estimated from data, recognizing that these do not incorporate distances of particles that moved beyond measurement distance of 3 m.

Slope
c Notation (−) means undefined or coefficient of variation is less than 0.01.
For data involving an estimated shape parameter A < 0 ( Figure 3a, Figure 4a), the relatively rapid decrease in the exceedance probability R x (x) with little indication of an asymptotic approach to R x (x) = 0 strongly suggests that the data represent bounded forms of the distribution f x (x) (Figure 3a). Nonetheless, one might on empirical grounds reasonably fit the data to, say, an exponential distribution (although quartile-quartile plots, not shown, would advise against this). However, the negative concavity of the fits are unambiguous in Figure 4a, strongly reinforcing the point that the data represent bounded forms of the distribution. This result is consistent with the idea of rapid thermal collapse for data involving θ ≤ 9 degrees and Ki ≤ 0.73.
Note that Ki < Ki * in all cases. For unknown reasons (and not attributable to truncation) the average particle motion on the flat surface is larger than the next three slope angles (3, 6 and 9 degrees). This may simply reflect the stochasticity of the disentrainment process at these small angles. Also note that several data sets share "kinks" in their estimated exceedance probabilities at similar distances, for example, around 0.8 m and 1.3 m. This could be due to chance, or it may reflect a persistent effect of the structure of the surface roughness, specifically the occurrence of relatively large roughness elements. Roth et al.
(2020) note this behavior in their field-based measurements of particle travel distances on vegetated hillslopes (Section 4.2), 5 where vegetation acts as roughness elements.
For data involving an estimated shape factor A ≥ 0 ( Figure 3b, Figure 4b), the first two sets (12 and 15 degrees) are approximately exponential with Ki ≈ Ki * . This is reflected in the close correspondence between the values of the scale parameter B and the estimated average travel distances µ x ( Table 2). The data show a clear asymptotic appearance of the exceedance probability R x (x) (Figure 3b) with essentially straight line fits in Figure 4b. This result is consistent with the idea that these two 10 experiments involved approximately isothermal conditions. For larger shape factor A, the fitted lines decrease in an exponentiallike manner in Figure 3b and they appear as essentially straight lines over the domain of measured travel distances in Figure   4b. The fits therefore cannot be used to distinguish between an exponential distribution and a generalized Pareto distribution.
Note, however, that regardless of the distribution, whereas the data for 18 degrees span about 90% of the probability of the distribution, the data for 21 degrees sample only about 50% of the probability contained in the distribution, and the data for 15 24 degrees sample only about 15% of the probability. This directly points to the difficulty of working with tail-censored data, especially if the underlying distribution is heavy-tailed (Appendix A). That is, at large Ki the experimental flume is sampling only a fraction of the distribution, representing just the shortest travel distances associated with the specific surface-roughness conditions. Nonetheless, the mechanical basis of the distribution combined with its consistency with data for A < 0 reinforces the merit of the hypothesis of a heavy-tailed behavior for A > 0. 20 Further note that for all cases less than 24 degrees the estimated values of A and B suggest that the distributions have finite moments. These moments are undefined (A > 1) for the case of 24 degrees. Also recall that only four particles stopped on the flume at a slope angle of 27 degrees, and no particles stopped at an angle of 30 degrees. These points are consistent with the idea that gravitational heating is systematically increasing relative to frictional cooling with increasing Kirkby number. Using the largest estimated value of µ (Table 2), the values of the Kirkby number would be Ki = 1.2 and Ki = 1.3 for the slope 5 angles of 27 and 30 degrees, respectively.
Here is an interesting sidebar. Following Samson et al. (1999) we plot the reciprocal, 1/µ x , of the average travel distance µ x versus slope angle θ ( Figure 5). For spheres rolling bumpety-bump down an inclined surface roughened with a quasi-random monolayer of small particles, Samson et al. (1999) show a linear decline in 1/µ x with increasing slope angle (see their Figure   3) associated with trapping (deposition) related to collisional friction. This reciprocal then smoothly transitions to values close 10 to zero with further increase in slope angle representing continuing motions without significant trapping. The plot in Figure   5 also reinforces the idea that if one assumes an exponential distribution to calculate average travel distances, the values will be similar to those associated with the generalized Pareto distribution for small to moderate Kirkby numbers, but then deviate significantly with increasing heaviness of the distribution tail. (1975) 15

Experiments
The experiments reported by Kirkby and Statham (1975) involved dropping particles onto an inclined surface with fixed roughness composed of particles of different sizes, thus giving different ratios of particle to roughness size. For different slopes, the particles were dropped from different initial heights h at the crest of the surface, and travel distances were measured. Here we focus on average travel distances reported in their Figure 1 involving a slope angle of 35 degrees and three dropped particle sizes with minor-to-intermediate axis ratios of 0.40, 0.53 and 0.77. Kirkby and Statham (1975) report that the distributions of travel distances are exponential, but do not present the travel distance data.
For A and B defined by Eq. (15) and Eq. (16) we may write the average travel distance as 5 where B * = (α/γ) sin 2 θ/µ cos θ or B = B * 2 h. The form of Eq. (38) requires a straight line fit between the average distance µ x and the drop height h, but it does not provide a unique fit. We must ensure that the estimated shape parameter A < 1 to yield a finite average; otherwise the comparison is not meaningful. Based on the results described in the previous section we assume that the Kirkby number is close to unity, and for the given slope angle we choose a friction factor of µ = 0.7. Whereas the coefficient of restitution may be relatively large for natural rock material, this coefficient typically decreases significantly on average for irregular particles due to the high probability that collisions are not collinear with particle centers of mass (see next section). For illustration we choose = 0.5. As before we fix γ and then vary α to estimate A and B * .

Results
The data are well fit by Eq. (38) ( Figure 6). Estimated parametric values are provided in Table 3.  Figure 1 of Kirkby and Statham (1975) for three different particle sizes. Note that we are assuming the largest size is 21.5 mm rather than the reported value of 0.215 mm. Table 3. Fitted values of the parameters for the data shown in Figure 6. We emphasize that other choices of the quantities , µ and α would provide equally good fits, given that Eq. (38) does not 15 provide a unique solution for A and B * . With this in mind, the estimated values of A suggest that the data are consistent with a heavy-tailed form of the generalized Pareto distribution with finite mean and finite variance. The dropped particles experience different rates of frictional cooling, manifest in the increasing value of α with increasing particle size (Table 3); and the data are consistent with the idea that the average travel distance is directly proportional to the initial energy state determined by the drop height h. 20 16 https://doi.org/10.5194/esurf-2020-99 Preprint. Discussion started: 8 December 2020 c Author(s) 2020. CC BY 4.0 License.

Experiments
We conducted two sets of experiments. The first set was aimed at demonstrating the basis for treating the proportion of energy extraction, β x , as a random variable. To do this we focused on the analogous quantity β z , which is the proportion of energy extraction associated with particle collision on a horizontal surface following vertical free fall. This allowed us to show, and calculate, the partitioning of translational energy into deformational friction and rotational energy during collisions. We used particles of varying size and angularity. The experiments involved dropping the particles onto a smooth rigid surface of slate, and onto a rough surface. The rough surface, made of concrete, had a granular roughness smaller than the particles (Figure 7).
We recorded particle motions using a Lightning RDT monochrome camera (DRS Technologies) operating at 800 frames per ligible rotational motion and with a prescribed surface-parallel velocity. We used two sizes (D ≈ 1 cm and D ≈ 0.5 cm) of irregular, rounded particles and one size (D ≈ 1 cm) of angular particles. We recorded motions of particles launched from the catapult onto the rough surface with high-speed imaging at a resolution of 640 × 640 pixels. In addition we made audio recordings of particle-surface interactions during their downslope motions. Audio recordings were made and processed using the GarageBand application on a 6th generation Apple iPad. Figure 8. Image of pendulum catapult. A particle is placed on the low friction (glossy cardboard) cradle at the base of the pendulum arms (∼ 20 cm); using a wand the arms and particle are pushed back to a preset position as one would a toddler on a playground swing (albeit not with a wand), then released; the arms are arrested by the front bar, whence the momentum of the particle launches it onto the surface. The cradle is about 2 mm above the rough surface. A stone is placed in the base of the catapult for stability.

Results
As a point of reference, the vertical rebounds of ordinary spherical glass marbles following their impacts on the hard slate reveal  In contrast, the rebounds of natural particles from the hard slate reveal how noncollinear collisions strongly influence the rebound angle and normal rebound height. If z denotes the effective normal coefficient of restitution, then β z = 1 − 2 z is the proportion of the normal (vertical) component of kinetic energy extracted, analogous to β x . That is, β z = −∆E p /E p . This proportion involves a mechanical loss due to particle and surface deformation during the collision, with conversion of energy to is a simple demonstration of the idea that β z (and β x ) must be treated as a random variable rather than a fixed (deterministic) quantity as in the case of the normal coefficient of restitution associated with spherical particles in a granular gas.
To illustrate this we write a normalized form of β z as β value associated with a collinear collision. Although we cannot offer a theoretical basis for the distribution of β * z , we nonetheless on empirical grounds fit β * z to the beta distribution because of its versatile bounded form over For the hard slate surface, spheres show a small variance about the mean value. Probability is then redistributed toward β z = 1 with increasing particle angularity ( Figure 9). For the rough surface, spheres show a larger variance, and there is a stronger redistribution of probability toward β z = 1 with increasing angularity ( Figure 10). We cannot directly map this result 20 to an interpretation of β x because of differences in the geometrical conditions of collisions. Nonetheless, as described below, the effect of angularity appears in measurements of travel distances as an increasing likelihood of disentrainment with increasing angularity.
For the rough experimental surface, β z is on average larger than for the hard slate surface. The particle-surface impact is unlike that of an idealized rigid surface, and more like that of a quasi-rigid (deformable) granular material. Nonetheless, despite 25 the small effective coefficient of restitution, noncollinear collisions yield significant conversion of energy to transverse motions and rotation. The particles do not simply "die" on impact.
For each particle-surface combination, the largest rebound heights provide an estimate of the (ordinary) normal coefficient of restitution . These heights are associated with approximately collinear collisions as confirmed by the high-speed imaging. We can therefore estimate the partitioning of energy between the frictional loss f c (deformation, heat and sound) and reflectional 30 transverse motion and rotation E c (Appendix D). Results indicate that on average less than half of the initial energy is dissipated by friction on the hard surface, and slightly less goes to transverse motion and rotation (Table 4). About 20-30% of the initial  High-speed imaging of particle motions launched by the catapult onto the rough surface provides estimates of initial surfaceparallel particle velocities u 0 ( Table 5). The imaging reveals that the free-flight distances before first collisions increase with increasing surface slope. The particles then experience complex collisions with the surface that randomize their motions, including the onset of rotation.  Note that, for reference below, two value of µ, α, Ki and Ki * are provided. The first value is based on the measured launch velocity u 0 ( Table 5). The second value is based on a reduced velocity. Namely, we do not know the average energy state of the particles at the truncation position used in the fitting of A and B, although it most likely is smaller than that associated with 5 the launch velocity. To provide a sense of the uncertainty in the calculated values we thus assume that the particle energy is reduced by half of its initial launch energy, although this is less likely with increasing surface slope and gravitational heating.
At small slopes S this adjustment has a larger effect on µ than on α. At large slopes this adjustment has a larger effect on α. Table 5. Slope-parallel launch velocities u0 measured from high-speed imaging.
Particles Slope u0 ( These examples are used in estimating the coefficient of restitution for the particle-slate collisions. The angular particle collision is a low probability event (dubbed the "pogo stick" by M. Schmeeckle) in that the point of 5 impact involves a particle corner that becomes aligned directly beneath the center of mass at the instant of impact. Table 6. Fitted and estimated values of the parameters for the data shown in Figure 11 with coefficients of variation in parentheses. b Estimated from data, recognizing that these do not incorporate distances of particles that moved beyond measurement distance of 1.9 m.
c Notation (−) means undefined or coefficient of variation is less than 0.01.
One of the more compelling results appearing in several of the videos is when the translational kinetic energy of a particle at first impact is converted to translational kinetic energy involving transverse motion and rotational kinetic energy, then during a second or third collision, rotational energy is converted back to vertical motion thence to gravitational potential energy. The likelihood of this occurring increases with particle angularity, where noncollinear collisions are the rule rather than 10 exception, and pointy particle corners lead to unusual collision configurations. The file "Angular_all_rotational.avi" shows a particularly strong conversion of translational to rotational motion with the initial collision on hard slate. The file "Angu-lar_rotational_to_vertical.avi" shows conversion of rotational to vertical motion with the second collision. The file "Semian-gular_rotational_die.avi" shows rapid cessation of motion following conversion of rotational to vertical motion.
The geometry of a noncollinear collision following the vertical drop of an angular particle is different from that of a particle negligible initial rotation and the motions start to become randomized, including the onset of rotation, following free flight and initial surface collisions. Rather than decelerating, gravitational heating maintains velocities similar to the launch velocities. Indeed, the particle on S = 0.18 seems likely to stop, but then continues with heating. For contrast the file "Rounded_0slope.avi" shows an example of a rounded particle that rapidly decelerates then "dies" when launched onto the flat rough surface (S = 0).

25
The increase in free-flight distances (before initial collisions) with increasing slope are apparent in the three videos.
Audio recordings of particle-surface interactions during their downslope motions reveal the distinctive clickety-click sounds of collisions ("Bouncing.m4a"), which are markedly different from the sounds emitted by particles that are either gently or forcefully made to slide on the rough surface ("Sliding.m4a"). These clickety-click sounds occur with high frequency, particularly when particles are in a tumbling (nominally "rolling") mode, giving way to decreasing frequency when particles 30 undergo runaway bouncing motions. In contrast, sliding motions emit continuous scraping sounds. The key result of these recordings is to audibly reinforce the idea that motions involve collisional friction rather than a sliding Coulomb-like behavior, except in association with the brief collision impulses as described in collision mechanics theory (Brach, 1991;Stronge, 2000).
Further analyses of these detailed particle motions in relation to downslope and cross-slope motions are to be reported elsewhere. As with the data of Gabet and Mendoza (2012), we again fit the parameters A and B, then calculate values of µ and α 10 assuming a value of γ. We also displace the starting position (x = 0) to the first inflections in the raw exceedance probability plots for the smaller two particle sizes, then recalculate distances and exceedance probabilities.

Results
The data are reasonably well fit by the theoretical curves ( Figure 12). Estimated parametric values are provided in Table 7, where estimates of the variability in µ, α, Ki and Ki * are based on a Monte Carlo analysis as described in Appendix C. Table 7. Fitted and estimated values of the parameters for the data shown in Figure 12 with coefficients of variation in parentheses. b Estimated from data, recognizing that these do not incorporate distances of particles that moved beyond measurement distance of 14 m.

Particle size (m)
c Notation (−) means coefficient of variation is less than 0.01.

15
For all particle diameters the Kirkby number Ki ≈ 1, and the fits are insensitive to the value of γ. Moreover, estimated values of the friction factor µ are similar; these do not show a systematic change with particle size. The estimated values of A suggest that the smallest particles represent a distribution with finite mean but undefined variance (1/2 ≤ A < 1). The larger two particle sizes represent conditions with an undefined mean and undefined variance (A ≥ 1).
In contrast to the ambiguity of an exponential versus a Pareto fit to the data of Gabet and Mendoza (2012) for A ≥ 0 20 (Figure 4b), the concavity in the semi-log exceedance probability plot (Figure 12) for the smaller two particle sizes is readily apparent and consistent with a Pareto distribution. Certainly one could fit a straight line to these data, but the fit would degrade (as revealed, although not shown, by quartile-quartile plots). Nonetheless, we must be cautious. Inasmuch as the theoretical distribution is the correct choice, then the data represent only a fraction of the total probability. For the smallest particles about 15% of the tail is not sampled. For the intermediate particles about 35% is not sampled; and for the largest particles about 70% of the tail is not sampled. This reinforces the difficulty of working with tail-censored data with relatively small sample sizes. 5 Note also that for the smallest particle size the average travel distance µ x calculated from the data is similar to the estimated parameter B but is significantly smaller than the average estimated from the values of A and B (Table 7). This occurs because the data "see" probability distributed in a manner that is not dissimilar from that expected for an exponential distribution, whereas the values of A and B incorporate information contained in the tail of the Pareto distribution.
Based on reported particle velocities, values of the initial average kinetic energies E a0 of the three particle sizes are 1.0 × Banana slugs, whose locomotive energetic costs are constrained by the shear-thinning rheology of their mucus (Lauga and 30 Hosoi, 2006), appear as slow moving Dirac functions in the power spectra of surface elevation. None were injured during the experiments.
As above, we fit the parameters A and B, then calculate values of µ and α assuming a value of γ. We also displace the starting position (x = 0) to the first inflection in the raw exceedance probability plots, then recalculate distances and exceedance probabilities. This displacement is applied only for cases involving relatively small travel distances (typically the smallest particles), and is only about one or two particle diameters. For travel distances involving tens of m, we focus the fitting on the central part of the data, deemphasizing the fit for small values and for the extreme tails. In addition, changes in surface slope 5 occur on all sites, and we restrict the data fitting to positions upslope of these changes. These changes occur at: 2.7 m (V14), 3.5 m (V20), 11 m (V24), 16.6 m (V39), 34 m (B17), 31 m (B20) and 33 m (B28). For site V0, 20 travel distances were measured for each particle size class. Initial examination indicated that the distributions of travel distances were similar, so we pooled these data. By dropping (rather than launching) the particles, initial energies are less certain than those in the experiments reported above. We calculate the impact velocities, assume a coefficient of restitution, then use the average downslope reflection 10 velocities. We note that uncertainty in the initial energies affects the estimates of the quantities µ and α but does not influence the values of the parameters A and B obtained from the data fitting.

Results
The data for the V sites are reasonably well fit by the theoretical curves ( Figure 13). In these two examples as well as in cases not plotted, the estimated parameter A systematically increases with particle size (Table 8), and reflects a transition from Table 8. Fitted and estimated values of the parameters for the data reported by Roth et al. (2020) as shown in Figure 13 and Figure 14 with coefficients of variation in parentheses.

Slope Particle
Site b Estimated from data, recognizing that these do not incorporate distances of particles that moved beyond positions of noted slope changes.
c Notation (−) means undefined or coefficient of variation is less than 0.01.

15
rapid thermal collapse (A < 0) to approximately isothermal conditions (A ≈ 0) or net heating (A > 0). Interestingly, travel distances on the steep V24 slope are systematically larger than on V14. Yet evidently the surface roughness on V24 leads to approximately isothermal conditions for the larger particles, whereas the roughness on V14 leads to net particle heating. The data for the B sites ( Figure 14) similarly show that the estimated parameter A systematically increases with particle size (Table   8), transitioning from rapid thermal collapse to net heating. Note that the fitted shape and scale parametric, A and B (Table 8), 20 are consistent in sign and approximate magnitude with those presented by Roth et al. (2020) for the same data set. This paper also presents exceedance probability plots for all 21 experiments (excluding the zero slope case). Also note that estimates of the variability in µ, α, Ki and Ki * are based on a Monte Carlo analysis as described in Appendix C.
The steepest burned site, B28, offers a further interesting perspective on particle motions. On this steep and relatively smooth hillslope, exceedance probabilities associated with all three particle sizes cannot be fit with reasonable fidelity by individual 25 curves. Rather, the data suggest a mingling of two particle behaviors -rapid cooling for many particles, and runaway heating for a second group leading to a pronounced heavy tail ( Figure 15) -in effect giving a mixed distribution. Namely, let x 1 and x 2 denote the travel distances of the two groups, and let p 1 denote the proportion represented by first group such that p 2 = 1 − p 1 is the proportion of the second group. The simplest form of a mixed distribution is 30 and the cumulative distribution is where primes denotes variables of integration. As above, the exceedance probability is R x (x) = 1 − F x (x).
For the three particle sizes, exceedance probabilities are well matched by the weighted sum of a nearly exponential distribution (A 1 ≈ 0) reflecting isothermal conditions and a heavy-tailed distribution (A 2 > 0) reflecting net heating (Table 9). Note, Table 9. Fitted and estimated values of the parameters for the data shown in Figure 15. Although in these three experiments (B28) the particles in each size group are nominally similar, we nonetheless suspect that the steep slope and relatively smooth surface give conditions that "filter" the particles into two subsets. One subset consists of particles whose motions are strongly randomized and become disentrained over short distances. The other subset consists 5 of particles whose motions by chance quickly transition to rotation and travel much longer distances over the smooth surface.
This filtering likely includes effects of the dropping of non-spherical particles. Namely, particles that are by chance initially dropped onto their relatively flat faces are less likely to transition to rotation and thus are more likely to travel short distances.
We also suspect the existence of a similarly nonuniform behavior in the Vanderbilt data, manifest as systematic variations in several of the exceedance probability plots (Figure 11).
Using the same data set, Roth et al. (2020) directly calculate the disentrainment rate function P x (x) using finite-differencing of the empirical cumulative distribution and exceedance probability functions. Although noisy, the data clearly illustrate the forms of P x (x) representing rapid thermal collapse (A < 0), approximately isothermal conditions (A ≈ 0) and net heating 5 (A > 0). Of particular note is the result that the roughness of natural vegetation exerts a strong cooling effect, and that the spatial structure of roughness elements together with local changes in surface slope can contribute to noticeable variations in travel distances about those expected for nominally uniform conditions.

Analysis
We emphasize at the outset a key point in comparing travel distances measured in experiments (laboratory or field) with (undefined) moments. We can never know this, although it might be suggested, for example, by the absence of convergence of estimated moments with increasing sample size or from multiple samples. The best we can do is to infer the veracity of the form of the distribution from descriptive statistics (e.g., exceedance probability plots, quartile-quartile plots), but this generally requires large data sets to support the tails of heavy-tailed distributions. In some of the comparisons above, there is the real possibility that calculated averages are just numbers associated with a distribution whose mean is undefined, such With these points in mind, we suggest that the fits presented above are consistent with the idea that each of the data sets represents a specific case of the generalized Pareto distribution. To further illustrate this idea we calculate the following quantities: 5 Based on Eq. (14), values of the modified exceedance probability R * and the dimensionless travel distance x * should collapse to a single straight line in a log-log plot with slope of −1 (Figure 16). Note that these plots magnify the deviations in the tails of the distribution. Also note that these fits suggest that all data, if plotted together, would collapse to the same line spanning more than three orders of magnitude of the dimensionless travel distance x * . This does not prove, but nonetheless supports, the idea that the generalized Pareto distribution correctly describes the energetics of the behavior of rarefied particle motions for a  (Table 6), the effect of particle angularity evidently appears as a difference in µ, consistent with µ ∼ β x and the measurements indicating that angular particles on average extract more translational energy during collisions than rounded particles (Figure 9, Figure 10). This also is consistent with the measurements of Gabet and Mendoza (2012) for a flat surface in that all values of µ in the Vanderbilt data are larger than the value of µ for a spherical particle in the Gabet and Mendoza data. We suspect that the small particles "feel" the roughness texture more than the larger rounded and angular 30 particles; and because the small particles are a mixture of rounded and angular shapes, the value of µ is similar to the angular particles. These differences in the values of µ persist with larger slopes, where the values of µ for rounded particles remains less than the other values (although we must be cautious to not overinterpret these differences given the uncertainty of the calculations). At zero slope the values of α are similar across particle shape and size. No similar systematic variation in α with particle size and shape is apparent, although rounded particles appear to be more readily heated with increasing slope.
The field-based measurements of DiBiase et al. (2017, Table 7) and Roth et al. (2020, Table 8) suggest that the friction factor µ is insensitive to particle size for the same slope and roughness conditions, and these data together with the laboratory measurements of Gabet and Mendoza (2012, Table 8) and the Vanderbilt data suggest that µ systematically varies with surface slope S (Figure 17). Note that the Vanderbilt data in this figure are based on the reduced velocity calculations (Table 6). For 5 completeness this figure is reproduced in Appendix E using the initial launch velocities u 0 ( Table 5).
Values of µ for Ki < 1 systematically fall above the 1:1 line (Figure 17), then converge to this line as Ki → 1. Using Eq.
(32) to estimate µ, evidently near S = 0 the second term on the right side of this equation dominates and gives positive µ with negative A for the smallest four slope angles in the data of Gabet and Mendoza (2012, Table 2) and the Vanderbilt data ( Table   6). The magnitude of this term then decreases (for A > 0) with increasing S such that µ ∼ S as Ki → 1 according to Eq. (8).
10 Note that Eq. (32) does not provide a physical explanation of the factor µ; it is just an estimate of µ based on the parameters A and B.
As summarized in Section 2.4, scaling suggests that the factor µ ∼ M (θ) is independent of particle size (Furbish et al., 2020a). This consistency with the experimental data reinforces the idea that the elements of µ ∼ M (θ), despite the complexity of the collisions involved, are akin to results from idealized collision mechanics, namely, that these elements are determined 15 by the coefficients associated with tangential and normal impulses, where particle size is not involved for a given slope and surface roughness. Recall that the expected dependency µ ∼ M (θ) on the slope angle θ arises because the expected surface normal impact velocity varies with this angle (Appendix E in Furbish et al., 2020a). However, this probably is insufficient to explain the relation in Figure 17. Unfortunately, we cannot further unfold the physical basis of µ in relation to particle-surface interactions beyond observing that the values of A and B return estimates of µ that are consistent with the expectation of its 20 slope dependency, independent of particle size. On empirical grounds, as Ki → 1 the factor µ → S reflects an approximate balance between heating and cooling with respect to translational energy. That is, the rate of extraction of translational energy (partitioned to all other forms) increases to match the rate of heating. Presumably this balance is exceeded (Ki 1) with slopes that are so steep that cooling is insufficient for deposition to occur -as in several of the experiments described above.
Turning to the quantity α, we plot the estimated values of this factor based on Eq.  (Table 6). For completeness this figure is reproduced in Appendix E using the initial launch velocities u 0 ( Table   5).
The data from Gabet and Mendoza (2012) support the idea that α systematically increases with Ki and becomes unbounded near Ki ∼ 1. The Vanderbilt data similarly support this idea. The data for all three particle sizes from DiBiase et al (2017) 30 involve Ki ≈ 1 and thus only reinforce the unbounded behavior of α. Similarly, the data from Roth et al. (2020) support this behavior. Note that values with large Ki and A ≥ 1 are not meaningful, as the underlying deposition length scales L c are undefined. Also note that because the Kirkby number Ki and the factor µ are specified in the fits of the data from Kirkby and Statham (1975, Figure 6) rather than being estimated from A and B * , we do not plot these values in Figure 18.
Recall that α reflects a direct 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. This translates to suppressing the disentrainment rate and increasing the deposition length scale L c , rewritten here as With E a = (m/2) u 2 p , the effect of particle mass m does not explicitly appear. This means that any effect of particle size must 5 appear in α or µ, or both. Similarly, any effect of particle angularity must appear in one or both of these quantities. The results above, with particular reference to the flat rough surface in the Vanderbilt experiments, suggest that effects of angularity appear in µ, whereas the data sets together suggest that effects of particle size primarily appear in α. We emphasize that the functional relation of α to Ki given by Eq. (35) is not definitive. Other functional forms are possible, although the basic form of Eq. (35) seems to be reasonably consistent with the data (Figure 18). Although not explicit in the formulation, we suspect that the effect 10 of heating includes an increasing partitioning of energy into rotational motion that is amplified for larger particles for a given slope and roughness, giving a decreasing likelihood of stopping as reflected in increasing α. Further disentangling the effects of α and µ must a await clearer mechanical basis for these quantities.
To reinforce the idea of a mixed distribution, consider the example of angular and rounded particles from the Vanderbilt experiments for S = 0 (Figure 11a). When pooled these data can be approximately fitted (not shown) to a generalized Pareto 15 distribution. However, the data are well fit using the mixed distribution defined by Eq. (39) (Figure 19). Note that this mixture of generalized Pareto distributions is not a generalized Pareto distribution. The distribution of travel distances of a mixture of particle sizes and shapes therefore must be described empirically or formed as a weighted mixture of distributions characterizing the behavior of the individual particle size and shape groups involved.
6 Discussion and conclusions 20 The laboratory and field-based measurements of particle travel distances presented above provide clear evidence that these distances are well described by a generalized Pareto distribution, where the form of the distribution reflects variations in particle behavior associated with the balance between gravitational heating and frictional cooling by particle-surface collisions. These behaviors vary from a bounded distribution associated with rapid thermal collapse to an exponential distribution representing approximately isothermal conditions to a heavy-tailed distribution associated with net heating of particles. Here we reiterate a 25 point made in the first companion paper (Furbish et al., 2020a). Namely, 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 probabilistic physics of the problem, and is based on a description of the kinetic 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.

30
The experiments involving high-speed imaging of particle motions reinforce what we intuitively already understand. Relative to a spherical particle, a rounded non-spherical particle is more likely to experience a noncollinear collision that converts the translational energy of free fall into transverse motion and rotational energy; and an angular particle is more likely than is a rounded particle to experience such conversions. The effect of this behavior is a systematic increase in the proportion β z with increasing angularity. Moreover, following the first free fall collision, an angular particle is more likely than is a rounded particle to experience a noncollinear collision that extracts either rotational or translational energy, or both. Translating this to surfaceparallel motions, a tumbling angular particle is more likely than is a rounded tumbling particle to experience a noncollinear 5 collision that extracts either translational energy or rotational energy, or both. Although we did not directly measure changes in surface-parallel energy associated with collisions, we can infer that the proportion β x likely systematically increases with increasing particle angularity as reflected in systematically shorter travel distances of angular particles relative to those of rounded particles ( Figure 11) on a surface with only a granular roughness texture. These experiments also illustrate the value of treating β x as a random variable. Although this quantity is related to the normal coefficient of restitution as used in granular gas theory, the complexity and richness of collisions and associated conversions of energy among modes necessitates a probabilistic description in this problem.
The essence of rapid thermal collapse (A < 0) involves the situation in which gravitational heating is absent or is insufficient to replace frictional cooling, particularly with angular particles and small particles. That is, a small tumbling particle is more likely than is a large tumbling particle to "see" the bumps and divots of the roughness texture at its scale, and to experience 15 collisions that arrest its motion. Indeed, this is the basic lesson of experiments involving spheres rolling bumpety-bump over monolayer roughness elements (Dippel et al., 1997;Samson et al., 1998Samson et al., , 1999, the experiments of Kirkby and Statham (1975) involving particles moving down surfaces with different granular roughness, and the experiments of Roth et al. (2020) involving the different roughnesses of vegetated and burned hillslopes. With increasing gravitational heating the transition to a heavy-tailed distribution of travel distances likely involves an increasing conversion of translational to rotational kinetic energy 20 leading to larger travel distances with decreasing effectiveness of collisional friction. In this regard the analysis points to the need for further clarity concerning how particle size and shape in concert with surface roughness influence the extraction of particle energy and the likelihood of deposition.
Although not essential to the fitting of particle travel distances to the generalized Pareto distribution, it nonetheless is desirable to have a clearer mechanical interpretation of the quantities µ and α and their relation to the Kirkby number Ki in 25 terms of particle properties and surface-roughness conditions, and the modes of particle motions. Here we note that the Pareto distribution with positive shape parameter A can be obtained as a mixture of exponential distributions whose rate parameters are distributed as a gamma distribution (Appendix F). This result suggests an interesting physical interpretation of the Pareto distribution of particle travel distances, and it also may indicate a strategy for clarifying how particle size and shape in concert with surface roughness influence the extraction of particle energy and the likelihood of deposition, inasmuch as the scale pa-30 rameter B is equivalent to the reciprocal of the expected disentrainment rate, E(P x ). For example, Roth et al. (2020) show that B systematically varies with particle size and surface slope based on the data described in Section 4.2.
We suggest that, in designing and conducting particle launching experiments, we have a propensity to select pretty particles, and rounded (if not spherical) particles are pretty. This is not a bad thing. But it skews our view of particle motions toward the behavior of rounded particles. The experiments clearly demonstrate that particle angularity matters in the disentrainment 35 process, specifically the likelihood of converting translational to rotational energy and the decreasing extraction of energy by collisional friction (Williams and Furbish, 2020).
Here in essence are the shortcomings of the formulation and its application to the experimental data sets of particle travel distances. We do not understand the transient (probabilistic) physics of near-launch conditions as particle motions become randomized, so there is uncertainty in choosing the truncation distance and the associated particle energy state. Similarly, little 5 is known about the distribution f Ep (E p ) of particle energy states E p and how this distribution might vary in the downslope direction. Thus, the assumption that the ratio γ of the arithmetic and harmonic means of the particle energies remains fixed may be incorrect. This is a parsimonious choice to close the formulation analytically. The friction factor µ is tentative. Namely, its essential element, the expected proportion of energy extraction β x , is consistent with the experimental results as reflected in the behavior of rounded versus angular particles, but the mechanical reasons for its asymptotic behavior, µ → S as Ki → 1 (Figure 17), remain unclear. Similarly the factor α is tentative. We need a clearer understanding of the elements that lead to increasing α and the associated lengthening of the deposition length scale L c , notably in relation to particle size. Many of the individual fits between the data and the generalized Pareto distribution in the exceedance probability plots are reasonably close; but some exhibit systematic deviations about the fitted distribution. As is usual in this situation, it is difficult to fully assess whether such misfits are related to stochasticity associated with small sample sizes (Appendix A) or to inadequacy of 15 the experimental design or to underlying flaws in the formulation leading to the generalized Pareto distribution. Likely all of these are involved. Despite the conceptual simplicity of particles moving bumpety-bump down a rough inclined surface, this is a hard problem.
We reemphasize that the work reported here is aimed at a probabilistic description of expected particle travel distances. This is a part of a larger effort to understand and inform the essence of the ingredients of nonlocal formulations of transport. We 20 are not suggesting that the results presented here can be immediately cast as a nonlocal formulation of transport. But in order to progress beyond current formulations, the probabilistic physics of particle motions merits closer examination. For example, this level of understanding provides the basis for justifying a Taylor expansion of the convolution (Furbish and Haff, 2010) to form a local Fokker-Planck-like description of transport assuming an exponential-like distribution of travel distances -with clarity regarding the limitations of this description. Furthermore, we have focused on the energetics of particles in motion. But 25 this is one of two ingredients of nonlocal formulations. The other involves the probabilistic physics and energetics of particle entrainment -a particularly difficult ingredient to constrain because of the difficulty of observing the entrainment process and because we do not yet know how to properly simulate this process. For this we must rely on theory and measurements of tracer particles in ways that have yet to be designed.
We end with a philosophical point. We enjoy eating our favorite tortilla chips, and mostly we enjoy them with a well prepared 30 dip, for example, spicy guacamole. But let us be honest. The experience then is no longer about the chips, it's about the dip.
The chips are just the guacamole delivery system. (Yumm.) Similarly, these companion papers nominally concern particle motions on inclined rough surfaces. But these particles are just the delivery system. The dip consists of the coherent statistical mechanics framework for describing the particle motions, and a demonstration that such a framework, albeit with rough edges, is possible. This represents a solid basis for subsequent efforts aimed at replication, falsification and refinement or replacement, 35 and possibly for fresh ideas concerning particle motions more generally.
Code and data availability. Emmanuel Gabet provided the data described in Section 3.1. The data described in Section 4.1 and Section 4.2 are available from Dibiase et al. (2017) and Roth et al. (2020). The data described in Section 3.3, including video and audio files, and function Fx u (x u ) =x u . In turn, the cumulative distribution function of the generalized Pareto distribution is Equating Fx u (x u ) and Eq. (A1) leads tô which provides an algorithm for generating values ofx drawn from the generalized Pareto distribution with shape and scale The MLE of a isã = 1/ã L wherẽ and the MLE estimate of b isb =b L /ã L whereb L is obtained from an iterative solution of These are biased estimates , but they nonetheless provide useful information concerning parameter estimation. This bias increases with decreasing sample size and with increasing censorship of the distribution tail. Moreover, the MLE may not converge near a = 0 nor if the sample has a coefficient of variation less than one. For reference below, whereas the Lomax distribution requires that a L > 0, the MLE given by Eq. (A3) may be negative, and therefore provides an estimate of a < 0 for the generalized Pareto distribution if b is known. Note also that in the limit of a → 0 the Pareto distribution is replaced with the exponential distribution with mean µx = b. The MLE of the mean µx of an exponential distribution is just the sample average, which is an unbiased estimate.
Now consider a sample size of N = 100, consistent with the data sets of Gabet and Mendoza (2012) and Roth et al. (2020).
We draw 10,000 samples then calculate and plot exceedance probabilities for varying values of the shape parameter a, holding Plots of estimated exceedance probabilities Rx(x) for all samples provide a visual sense of the variability in these probabilities associated with the inherent randomness in drawing samples ofx from the known distribution ( Figure A1). The animation mentioned above shows the outcome of successive draws, and nicely illustrates that many draws, by chance, bear little re- 15 semblance to the theoretically expected exceedance probability function as well as, in particular, the squirrelly behavior of values in the tails. (We hope that the animation drives home the point to avoid over-fitting and over-interpreting data in the tail of a heavy-tailed distribution with small sample size N .) Here are key items to consider. First, the variability in calculated exceedance probabilities increases with a, that is, with increasing heaviness of the distribution tail. This is not surprising, as a finite sample size represents a decreasing proportion of the total probability in the distribution as a increases. Second, the 20 variability in the MLE of a increases with increasing a (Figure A2), reflecting the first point above. Aside from the bias of these estimates, the difference between estimates of a and the true value can be large, although proportional differences are similar across values of a. For example, with a = −0.5, about 32% of the estimated values exceed a difference of ±10% from the true value; with a = 0, about 32% exceed a difference of ±10% about the true value of b = 1; and with a = 1, about 31% of the estimated values exceed a ±10% difference. Third, the variability in the exceedance probability plots and the MLEs decreases 25 -that is, these values converge to the true values -only when N approaches 10,000 or more ( Figure A3). Fourth, with increasing censorship of the data, the MLEs based on the uncensored values become strongly biased ( Figure A4). Moreover, this simple demonstration of the inherent variability in estimates of a does not involve the collinear effects and added variability associated with simultaneously estimating b. Fifth, despite the variability in the exceedance probability plots, the sign of the concavity of the plots for large positive or negative a is clear. However, near isothermal conditions (a = 0), individual samples 30 could appear to represent net particle heating when in actuality conditions of net cooling exist, and vice versa. Note that in the example of a → 0 ( Figure A1, Figure A2), we calculate sample exceedance probabilities andb for the exponential distribution.
According to the central limit theorem, values ofb are approximately normally distributed with variance ∼ σ 2 x /N . Here is an important sidebar. In the presence of an exact theory that predicts the values of a and b, one can appeal to, for example, the central limit theorem or a Monte Carlo analysis (as above) or MLE methods or bootstrapping to assign socalled confidence estimates associated with these known values of a and b. These specifically give information regarding the 5 likelihood that a sample of size N will yield values of a and b that differ from the true values. In contrast, in the absence of an exact theory and a priori knowledge of the true values of a and b, one can construct similar confidence estimates based on values of a and b estimated from a single sample. These specifically give information regarding the likelihood that a second sample of size N will yield values of a and b that differ from those estimated from the first sample. But no method -besides making N → ∞ (Cam, 1990) -can provide information regarding how close the first set of estimated values (or the second 10 set) is to the true unknown values. Alas, the literature is awash with confidence estimates, based on a single sample, incorrectly interpreted as measures of likelihood of containing the unknown values (Amrhein et al., 2019).
We therefore reemphasize our objective. At this stage of our work we are aimed at reasonable estimates of the shape and scale parameters in order to demonstrate the existence of the behaviors -rapid thermal collapse, isothermal conditions, and net heating of particles -represented by the generalized Pareto distribution. Refined values of these parameters are not needed 15 until we possess a clearer understanding of the mechanics of deposition. Semi-log plots highlight deviations in the tails, and provide a clear sense of the concavity that discriminates between cooling and heating. Log-log plots highlight deviations near the origin, and provide a sense of the log-linear decay of the tails for heavy-tailed conditions. The variability in the tails of the distribution as outlined above emphasizes the importance of avoiding over-fitting of the tails in visual fitting (or in any other method of fitting). 20 For comparison with our fits, we return to dimensional quantities and compute the MLE values of A and B (Table A1).
The MLE is implemented in the "flomax" algorithm in the Renewal Method for Extreme Values Extrapolation library of the R Project for Statistical Computing; it is implemented in the "gpfit" algorithm of the MATLAB programming language; or it can be coded directly from Eq. (A3) and Eq. (A4). With reference to Table A1, the MLE algorithm converges in all cases using the MATLAB "gpfit" algorithm (but not the R "flomax" algorithm). However, it returns poor (sometimes nonsensical) 25 estimates for A −1/2 or near A ≈ 0. The MLE degrades with increasing A and increasing censorship. Also note that the MLE estimates do not necessarily improve the fits ( Figure A5), likely due to the relatively small sample sizes and the likelihood that the data represent samples that are strongly censored, that is, where a significant proportion of the distribution is contained in that part of the tail that is not sampled.
One alternatively can choose, say, a nonlinear least-squares fitting algorithm that weights various parts of data differently, 30 emphasizing or deemphasizing values near the origin or in the tails. We suggest, however, that this is just a rule-based version of visual fitting. We also note that visual fitting is not directly influenced by censorship, although the form of the censored tail can never be known (Ballio et al., 2019). Bringing more sophisticated techniques to bear (e.g., Hosking and Wallis, 1987;Castillo and Hadi, 1997;Cramer and Schmiedt, 2011;Giles et al., , 2016Pak and Mahmoudi, 2018) to refine estimates of A and B is premature. There is a need to collect larger data sets, avoiding censorship if possible, and only then aim at refined 35 estimates of the parametric values as their theoretical basis is improved. Moreover, any real data set is not immune from the possibility, by chance, of representing a misfit from the underlying distribution and yielding parametric estimates that markedly differ from this underlying distribution -just as with the numerical examples above. But no formal quantitative analysis can reveal or fix this misfit.
We end with a cautionary note: Parametric values of a heavy-tailed distribution estimated from a data set with N < 1000 (and 5 possible with N < 10, 000), if presented as being precise, merit a healthy skepticism, particularly if the tail of the distribution is censored. Korup et al. (2012) address a related point in demonstrating how the exponents in power functions involved in Table A1. Fitted and estimated values of the parameters for the data reported by Gabet and Mendoza (2012) Cam (1990), "Never trust an estimate which is thrown out of whack if you suppress a single observation." Stumpf and Porter (2012) suggest that more generally statistically fitted power laws have little more than anecdotal value in the absence of a theoretical basis. Consider as an example the initial exceedance probability plot for the angular particles on a flat surface ( Figure B1a), which shows a clear inflection at about 5 cm. In this example, high-speed imaging of particles launched from the catapult reveal that the particles consistently travel ∼ 2 cm horizontally before their first collisions with the surface, as expected from calculations using Newton's second law for measured initial velocities. (Free-flight distances increase with increasing surface slope.) These initial flights involve negligible rotational motion. The particles then experience widely varying changes in their motions over the next 2-3 cm following the first collisions, often with the onset of rotational motion. In the text we suggest that the inflection in the exceedance probability plot reflects the uniformity of the launch velocities followed by a finite distance over which randomization of the motions occurs. This is consistent with the idea that the factor γ ∼ 1 before randomization occurs, giving 5 an initial disentrainment rate that is smaller than after randomization. However, we note that the inflection also may involve other effects.
For these reasons we truncate the plots at the inflection position, then recalculate exceedance probabilities with reduced N ( Figure B1b). In this example the truncation distance is a significant proportion of the total travel distances. However, this truncation distance is less important when effects of initial (near-launch) conditions occur over a distance that is small relative 10 to total travel distances. Unfolding the details of the physics of particle motions over short distances is for a later time.

Appendix C: Uncertainty in calculated quantities
Of interest is how uncertainty in the estimates of the shape and scale parameters, A and B, propagates to uncertainty in the calculated values of µ, α, Ki and Ki * . Because A and B are obtained by visual fitting, for illustration we conservatively assume that the standard deviations of these values associated with a great number of samples of similar size N vary as A/ √ N and 15 B/ √ N based on Monte Carlo simulations as described in Appendix A. This effectively assumes the coefficient of variation formed by the sampling standard deviation is ∼ 1/ √ N . This provides the basis for gaining a sense of the relative magnitude of the variability in the calculations of µ, α, Ki and Ki * . For the Vanderbilt data (Section 3.3) we also incorporate the uncertainty in the launch velocities u 0 provided in Table 5.
We We emphasize that these calculations provide a sense of the variability only associated with that of A and B as this cascades through the successive calculations of µ, α, Ki and Ki * . For example, based on Eq. (32) the variability in µ reflects that in both 25 A and B. Based on Eq. (33), the variability in α reflects that in B and the variability in µ previously calculated. Also note that, starting with Eq. (32), as the slope S increases the relative contribution of this fixed term to calculated values of µ increases.
These calculations do not represent the natural variability in µ, α, Ki and Ki * if these quantities somehow could be measured directly, independently of A and B.
In general, calculated coefficients of variation decrease with increasing surface slope. Coefficients of variation are on the 30 order of 10% or more for smaller slopes in the experiments reported by Gabet and Mendoza (2012)  where T is the travel time to the second collision. We assume as an approximation that f c = (1 − 2 )E p0 , where is the normal coefficient of restitution. This effectively assumes that the energy loss due to particle-surface deformation is the same as that of a collinear collision. Now Eq. (D1) becomes We can experimentally determine , E p0 and T for individual particles. However, we generally cannot determine u 1 from side 15 imaging (except by chance with motion transverse to the camera). We also cannot readily determine I for irregular particles, nor ω from side imaging. Solving Eq. (D3) for the last two terms then represents the conversion E c of translational kinetic energy just prior to impact into translational energy associated with surface parallel motion and rotational energy. That is, Appendix E: Effect of launch velocity 20 For completeness, here we show plots of the friction factor µ versus the slope S and the factor α versus the Kirkby number Ki using the initial launch velocity u 0 rather than a reduced velocity in the calculations as presented in Figure 17 and Figure 18.
Values of µ, notable at smaller slopes S ( Figure E1), are noticeably larger than those calculated in Figure 17. The values appear to converge to the 1:1 line with increasing slope as Ki → 1, as in Figure 17.
Values of α and associated Kirkby numbers tend to be smaller ( Figure E2) than those calculated in Figure 18. However, the 25 overall variation between α and Ki is similar.
Appendix F: The Pareto distribution as a mixture of exponential distributions It is well known that a Pareto distribution with positive shape parameter can be obtained as a mixture of exponential distributions whose rate parameters are distributed as a gamma distribution. This result suggests an interesting physical interpretation of the Pareto distribution of particle travel distances, and it also may indicate a strategy for clarifying how particle size and shape in concert with surface roughness influence the extraction of particle energy and the likelihood of deposition. For completeness we therefore offer the following.
Recall that for an exponential distribution of travel distances x the fixed disentrainment rate is P x = 1/µ x . We then write 5 the conditional distribution as We may now treat the rate P x as a random variable that is distributed as a gamma distribution, namely, with shape parameter A g and scale parameter B g . The unconditional distribution of travel distances x is then obtained as a 10 gamma weighting of the conditional exponential distribution, namely, Substituting Eq. (F1) and Eq. (F2) into Eq. (F3) and evaluating the integral then leads to This is a Lomax distribution (compare with Eq. (28)) with shape parameter A g = 1/A and scale parameter B g = B/A = BA g .

15
The expected value µ Px = A g /B g = 1/B and the variance is σ 2 Px = A g /B 2 g = A/B 2 . This immediately implies that an experimental estimate of B provides an estimate of the expected disentrainment rate E(P x ) = µ Px ; and an estimate of A together with B provides an estimate of the variance of P x .
Because P x is a random variable, and because the exponential distribution, Eq. (F1), implies isothermal conditions, we now use Eq. (16) to write 20 P x = γmgµ cos θ αE a0 which is the reciprocal of the deposition length L c with specified average energy E a0 . For a given particle mass m, slope angle θ and energy E a0 , the quantities γ, µ and α are random variables. That is, we may envision an ensemble of combinations of these quantities, each of which yields isothermal conditions. In turn, envision a great number (cohort) of particles. Then f Px (P x ; α, β)dP x is the probability that particles possess the value P x . These particles are deposited exponentially with mean 25 µ x = 1/P x . When combined with all other exponential distributions with varying means (i.e., different combinations of γ, µ and α), the collective effect is a Pareto distribution. Each subset of particles with rate P x behaves isothermally, but collectively the downslope energy variation of the entire cohort involves net heating.
As an example, for a value of A = 0.01 representing nearly isothermal conditions, the gamma distribution of P x is centered about the value of 1/B ( Figure E1a) and approaches a Dirac function in the limit of A → 0 as the variance σ 2 Px = A/B 2 → 0.
This represents the exponential limit of a Pareto (or Lomax) distribution. As A increases, the distribution of the disentrainment rate P x becomes increasingly skewed toward P x = 0, which collectively gives a heavy-tailed Pareto distribution. In turn, the 5 distribution of the reciprocal µ x = 1/P x is given by the inverse gamma distribution ( Figure E1b). Again, in the limit of A → 0 the mean travel distances µ x of the mixture of exponential distributions converge to the mean value of the Pareto distribution, namely, B. With increasing A the mixture of mean values is distributed about the mean, notably incorporating increasingly larger values of µ x .
Inasmuch as the quantities γ, µ and α can be related to measurable quantities -for example, particle size, particle shape and 10 surface roughness -then Eq. (F5) suggests the possibility of formulating a multiplicative relation between these properties and the shape parameter B = 1/E(P x ). An initial effort to this effect is reported by Roth et al. (2020).
Author contributions. All authors contributed to the conceptualization of the problem and its technical elements. SGW led the VU experiments. DJF wrote much of the paper with contributions by the other authors.
Competing interests. We have no competing interests.