Rarefied particle motions on hillslopes – Part 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., 2021a), is based on a description of the kinetic energy balance of a cohort of particles treated as a rarefied granular gas, as well as a description of particle deposition that depends on the energy state 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 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.

1 Introduction

As described in our first companion paper (Furbish et al., 2021a), 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., 2011, 2013; DiBiase and Lamb, 2013; DiBiase et al., 2017; Doane et al., 2018, 2019), and the motions of rockfall 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 rockfall 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 (e.g., ravel involving many particles) but where during the collective motions of many particles the effects of particle–surface interactions dominate over effects of particle–particle interactions in determining the behavior of the particles – akin to granular shear flows at high Knudsen number (Kumaran, 2005, 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., 2021a) 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. 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 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 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., 2021a). In Sect. 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 basis of the theory leading to the generalized Pareto distribution of particle travel distances. In Sect. 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 new laboratory experiments designed to clarify how the size and shape of particles influence their motions and disentrainment based on high-speed imaging. In Sect. 4 we compare the formulation with the field-based measurements of travel distances reported by DiBiase et al. (2017) and Roth et al. (2020).

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., 2021b) we show that the generalized Pareto distribution in this problem is a maximum entropy distribution (Jaynes, 1957a, b) 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 arrangement. In the fourth companion paper (Furbish and Doane, 2021) 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

2.1 Probabilistic description of disentrainment

The problem of describing rarefied particle motions on hillslopes is motivated by the entrainment forms of the flux and the 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^{\prime }=x-r$, the flux q(x) may be written as

$$\begin{array}{}\text{(1)}& q\left(x\right)=\underset{-\mathrm{\infty }}{\overset{x}{\int }}{E}_{\mathrm{s}}\left(x^{\prime }\right)R_r(x-x^{\prime };x^{\prime })\mathrm{d}{x}^{\prime },\end{array}$$

where E_s(x) denotes the volumetric entrainment rate at position x. In turn, letting ζ(x,t) denote the local land-surface elevation, the entrainment form of the Exner equation is

$$\begin{array}{}\text{(2)}& c_{\mathrm{s}}{\displaystyle \frac{\partial \mathit{\zeta }(x,t)}{\partial t}}=-E_{\mathrm{s}}\left(x\right)+\underset{-\mathrm{\infty }}{\overset{x}{\int }}{E}_{\mathrm{s}}\left(x^{\prime }\right)f_r(x-x^{\prime };x^{\prime })\mathrm{d}{x}^{\prime },\end{array}$$

where $c_{\mathrm{s}}=\mathrm{1}-{\mathit{\varphi }}_{\mathrm{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

$$\begin{array}{}\text{(3)}& P_r(r;x)={\displaystyle \frac{f_r(r;x)}{R_r(r;x)}},\end{array}$$

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, Eqs. (1) and (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 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=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$ 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

$$\begin{array}{}\text{(4)}& f_k\left(k\right)=p_k\prod _{i=\mathrm{1}}^{k-\mathrm{1}}(\mathrm{1}-p_i).\end{array}$$

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 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. (2021a) 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 kinetic energy balance of rarefied particle motions and how this balance determines the deposition of particles in relation to their energy state.

Figure 1 Definition diagram of surface inclined at angle θ and control volume with edge length dx through which particles move. Figure reproduced from companion paper (Furbish et al., 2021a).

2.2 Energy and mass balances

Consider a rough, inclined surface with uniform slope angle θ (Fig. 1). At this juncture we simplify the notation and consider the motions of particles entrained at a single position x=0. Now the particle travel distance r→x and the probability density function $f_r(r;x)\to f_x\left(x\right)$. 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., 2021a, Appendix B). That is, let N(x) denote the number of particles that enter the control volume, and let N(x+dx) denote the number that leaves the volume. The number of particles deposited within the control volume is $\mathrm{d}N=N(x+\mathrm{d}x)-N\left(x\right)$. The objective is then to determine the rate of particle deposition, $\mathrm{d}N/\mathrm{d}x$, based on the energy state of the cohort of particles.

In particular, let $E_{\mathrm{p}}=(m/\mathrm{2})u^{\mathrm{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_{\mathrm{h}}=\mathrm{1}/\langle \mathrm{1}/E_{\mathrm{p}}\rangle$. The total energy E=NE_a. The formulation presented in the first companion paper (Furbish et al., 2021a) 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

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}E\left(x\right)}{\mathrm{d}x}}=Nmg\sin \mathit{\theta }-Nmg\mathit{\mu }\cos \mathit{\theta }-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}Nmg\mathit{\mu }\cos \mathit{\theta }.\end{array}$$

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

$$\begin{array}{}\text{(6)}& \mathit{\mu }={\displaystyle \frac{\langle {\mathit{\beta }}_x\rangle }{\mathrm{4}\tan \mathit{\varphi }}},\end{array}$$

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

$$\begin{array}{}\text{(7)}& L_{\mathrm{c}}={\displaystyle \frac{\mathit{\alpha }{E}_{\mathrm{h}}}{mg\mathit{\mu }\cos \mathit{\theta }}},\end{array}$$

and it is a function α=f(Ki) of the dimensionless Kirkby number Ki defined by

$$\begin{array}{}\text{(8)}& \mathit{Ki}={\displaystyle \frac{mg\sin \mathit{\theta }}{mg\mathit{\mu }\cos \mathit{\theta }}}={\displaystyle \frac{S}{\mathit{\mu }}},\end{array}$$

which represents the ratio of gravitational heating to frictional cooling.

Conservation of particle mass is given by

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}N\left(x\right)}{\mathrm{d}x}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }{E}_{\mathrm{h}}}}Nmg\mathit{\mu }\cos \mathit{\theta }=-{\displaystyle \frac{N}{L_{\mathrm{c}}}},\end{array}$$

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

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{E}_{\mathrm{a}}\left(x\right)}{\mathrm{d}x}}& =mg\sin \mathit{\theta }-mg\mathit{\mu }\cos \mathit{\theta }\\ & +{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}mg\mathit{\mu }\cos \mathit{\theta }\left({\displaystyle \frac{E_{\mathrm{a}}}{E_{\mathrm{h}}}}-\mathrm{1}\right),\end{array}\end{array}$$

where the arithmetic and harmonic averages are related as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{E_{\mathrm{a}}}{E_{\mathrm{h}}}}=\mathit{\gamma }\ge \mathrm{1}.\end{array}$$

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.

2.3 Generalized Pareto distribution of particle travel distances

Simultaneous solution of Eqs. (5), (9) and (10) using Eq. (11) leads to the disentrainment rate function

$$\begin{array}{}\text{(12)}& P_x\left(x\right)={\displaystyle \frac{\mathrm{1}}{Ax+B}}.\end{array}$$

In turn the probability density function f_x(x) of travel distances x for particles starting at x=0 is

$$\begin{array}{}\text{(13)}& f_x\left(x\right)={\displaystyle \frac{B^{\mathrm{1}/A}}{(Ax+B{)}^{\mathrm{1}+\mathrm{1}/A}}},\end{array}$$

and the exceedance probability function is

$$\begin{array}{}\text{(14)}& R_x\left(x\right)=\left\{\begin{array}{ll}\frac{B^{\mathrm{1}/A}}{(Ax+B{)}^{\mathrm{1}/A}}& A\ne \mathrm{0}\\ e^{-x/B}& A=\mathrm{0}.\end{array}\right.\end{array}$$

The shape and scale parameters A and B are

$$\begin{array}{}\text{(15)}& A={\displaystyle \frac{\mathit{\alpha }}{\mathit{\gamma }}}\left[{\displaystyle \frac{S}{\mathit{\mu }}}-\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}(\mathit{\gamma }-\mathrm{1})\right]\mathrm{and}\end{array}$$

$$\begin{array}{}\text{(16)}& B={\displaystyle \frac{\mathit{\alpha }}{\mathit{\gamma }}}{\displaystyle \frac{E_{\mathrm{a}\mathrm{0}}}{mg\mathit{\mu }\cos \mathit{\theta }}}.\end{array}$$

The mean of the distribution is

$$\begin{array}{}\text{(17)}& {\mathit{\mu }}_x={\displaystyle \frac{B}{\mathrm{1}-A}}A<\mathrm{1}.\end{array}$$

With reference to the presentation in the first companion paper (Furbish et al., 2021a) and for the purpose of presenting results below, let E_a0 denote the initial average particle energy at x=0 and let N₀ denote the initial number of particles at x=0. In turn we define a characteristic cooling distance $X=E_{\mathrm{a}\mathrm{0}}/mg\mathit{\mu }\cos \mathit{\theta }$. For plotting purposes, and to highlight the role of the Kirkby number, we now define the following dimensionless quantities denoted by circumflexes:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}x& =X\widehat{x},N=N_{\mathrm{0}}\widehat{N},E=N_{\mathrm{0}}{E}_{\mathrm{a}\mathrm{0}}\widehat{E},\\ & E_{\mathrm{a}}=E_{\mathrm{a}\mathrm{0}}{\widehat{E}}_{\mathrm{a}}\mathrm{and}{E}_{\mathrm{h}}=E_{\mathrm{a}\mathrm{0}}{\widehat{E}}_{\mathrm{h}}.\end{array}\end{array}$$

With these definitions in place, Eqs. (5), (9), (10) and (11) are rewritten as

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\mathrm{d}\widehat{E}\left(\widehat{x}\right)}{\mathrm{d}\widehat{x}}}=\left[\mathit{Ki}-\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\right)\right]\widehat{N},\text{(20)}& {\displaystyle \frac{\mathrm{d}\widehat{N}\left(\widehat{x}\right)}{\mathrm{d}\widehat{x}}}=-{\displaystyle \frac{\widehat{N}}{\mathit{\alpha }{\widehat{E}}_{\mathrm{h}}}},\text{(21)}& {\displaystyle \frac{\mathrm{d}{\widehat{E}}_{\mathrm{a}}\left(\widehat{x}\right)}{\mathrm{d}\widehat{x}}}=\mathit{Ki}-\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\left({\displaystyle \frac{{\widehat{E}}_{\mathrm{a}}}{{\widehat{E}}_{\mathrm{h}}}}-\mathrm{1}\right)\mathrm{and}\text{(22)}& {\displaystyle \frac{{\widehat{E}}_{\mathrm{a}}}{{\widehat{E}}_{\mathrm{h}}}}=\mathit{\gamma }\ge \mathrm{1}.\end{array}$$

The expressions involving energy, Eqs. (19) and (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,

$$\begin{array}{}\text{(23)}& {\mathit{Ki}}_{*}=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}(\mathit{\gamma }-\mathrm{1}).\end{array}$$

If $\mathit{Ki}<{\mathit{Ki}}_{*}$ then net cooling occurs, and if $\mathit{Ki}>{\mathit{Ki}}_{*}$ then net heating occurs. The condition $\mathit{Ki}={\mathit{Ki}}_{*}$ implies isothermal conditions such that $\mathrm{d}{\widehat{E}}_{\mathrm{a}}/\mathrm{d}\widehat{x}=\mathrm{0}$.

Figure 2 Plot of dimensionless probability density $f_{\widehat{x}}\left(\widehat{x}\right)$ versus dimensionless travel distance $\widehat{x}$ for scale parameter b=1 and different values of the shape parameter a for (a) a<0 and (b) a≥0 with associated exceedance probability plots (insets). Figure reproduced from companion paper (Furbish et al., 2021a). Compare with Fig. 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 Fig. 2.

¹ Truncation occurs at dimensionless distance $\widehat{x}=b/\left|a\right|$. ² Triangular with $a=-\mathrm{1}/\mathrm{2}$.

The dimensionless disentrainment rate is

$$\begin{array}{}\text{(24)}& P_{\widehat{x}}\left(\widehat{x}\right)={\displaystyle \frac{\mathrm{1}}{a\widehat{x}+b}}.\end{array}$$

The dimensionless probability density function of travel distances is

$$\begin{array}{}\text{(25)}& f_{\widehat{x}}\left(\widehat{x}\right)={\displaystyle \frac{b^{\mathrm{1}/a}}{(a\widehat{x}+b{)}^{\mathrm{1}+\mathrm{1}/a}}},\end{array}$$

and the exceedance probability function is

$$\begin{array}{}\text{(26)}& R_{\widehat{x}}\left(\widehat{x}\right)=\left\{\begin{array}{ll}\frac{b^{\mathrm{1}/a}}{(a\widehat{x}+b{)}^{\mathrm{1}/a}}& a\ne \mathrm{0}\\ e^{-\widehat{x}/b}& a=\mathrm{0}.\end{array}\right.\end{array}$$

The shape and scale parameters a and b are

$$\begin{array}{}\text{(27)}& a=A\mathrm{and}b={\displaystyle \frac{\mathit{\alpha }}{\mathit{\gamma }}}{\widehat{E}}_{\mathrm{a}\mathrm{0}}.\end{array}$$

The distribution $f_{\widehat{x}}\left(\widehat{x}\right)$ 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 Fig. 2, for a<0 this distribution decays more rapidly than an exponential distribution and is bounded at the position $\widehat{x}=b/\left|a\right|$. For a=0 it becomes an exponential distribution. For $\mathrm{0}<a<\mathrm{1}/\mathrm{2}$ the distribution $f_{\widehat{x}}\left(\widehat{x}\right)$ decays more slowly than an exponential distribution, but it possesses a finite mean and a finite variance. For $\mathrm{1}/\mathrm{2}\le a<\mathrm{1}$ the distribution possesses a finite mean but its variance is undefined. For a≥1 the mean and variance of $f_{\widehat{x}}\left(\widehat{x}\right)$ are both undefined, even though this distribution properly integrates to unity. For a>0 the tail of $f_{\widehat{x}}\left(\widehat{x}\right)$ decays as a power function, namely, $f_{\widehat{x}}\left(\widehat{x}\right)\sim {\widehat{x}}^{-(\mathrm{1}+\mathrm{1}/a)}$. The exceedance probability decays as $R_{\widehat{x}}\left(\widehat{x}\right)\sim {\widehat{x}}^{-\mathrm{1}/a}$. These results are summarized in Table 1. If the shape and scale parameters a and b are redefined as $a_{\mathrm{L}}=\mathrm{1}/a$ and $b_{\mathrm{L}}=b/a=b{a}_{\mathrm{L}}$, then Eq. (25) becomes

$$\begin{array}{}\text{(28)}& f_{\widehat{x}}\left(\widehat{x}\right)={\displaystyle \frac{a_{\mathrm{L}}{b}_{\mathrm{L}}^{a_{\mathrm{L}}}}{(\widehat{x}+b_{\mathrm{L}}{)}^{\mathrm{1}+a_{\mathrm{L}}}}}{a}_{\mathrm{L}},b_{\mathrm{L}}>\mathrm{0},\end{array}$$

which is a Lomax distribution with mean

$$\begin{array}{}\text{(29)}& {\mathit{\mu }}_{\widehat{x}}={\displaystyle \frac{b_{\mathrm{L}}}{a_{\mathrm{L}}-\mathrm{1}}}{a}_{\mathrm{L}}>\mathrm{1}.\end{array}$$

With a<0 the bounded form of $f_{\widehat{x}}\left(\widehat{x}\right)$ (Fig. 2) represents rapid thermal collapse with net frictional cooling. For a=0 the exponential form of this distribution represents isothermal conditions where frictional cooling is matched by gravitational heating and the apparent heating due to deposition. For a>0 the heavy-tailed form of $f_{\widehat{x}}\left(\widehat{x}\right)$ represents net gravitational heating.

For reference to data fitting presented below, a binomial expansion of Eq. (13) gives

$$\begin{array}{}\text{(30)}& f_x\left(x\right)={\displaystyle \frac{\mathrm{1}}{B}}[\mathrm{1}-{\displaystyle \frac{A}{B}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{A}}\right)x-\mathrm{\dots }].\end{array}$$

The expansion of an exponential distribution with mean μ_x is

$$\begin{array}{}\text{(31)}& f_x\left(x\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_x}}\left(\mathrm{1}-{\displaystyle \frac{x}{{\mathit{\mu }}_x}}+\mathrm{\dots }\right).\end{array}$$

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_x\approx \mathrm{1}/B\approx \mathrm{1}/{\mathit{\mu }}_x$. Moreover, if the travel distance data sample over a majority of the probability contained in the distribution, and if the tail of 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_x\left(x\right)=\mathrm{1}-r_x/(N+\mathrm{1})$, 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 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 quantile–quantile (Q–Q) 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 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., 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, Eqs. (15) and (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

$$\begin{array}{}\text{(32)}& \mathit{\mu }=S-{\displaystyle \frac{E_{\mathrm{a}\mathrm{0}}(A-\mathrm{1}+\mathrm{1}/\mathit{\gamma })}{Bmg\cos \mathit{\theta }}},\end{array}$$

and then we obtain α as

$$\begin{array}{}\text{(33)}& \mathit{\alpha }={\displaystyle \frac{B\mathit{\gamma }mg\mathit{\mu }\cos \mathit{\theta }}{E_{\mathrm{a}\mathrm{0}}}}.\end{array}$$

In turn the Kirkby number Ki is calculated from Eq. (8).

2.4 Elements of collisional friction

The quantities μ and α defined in relation to Eqs. (6) and (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 with a particle–surface collision can be expressed as $\mathrm{\Delta }{E}_{\mathrm{p}}=-{\mathit{\beta }}_x{E}_{\mathrm{p}}$ so that ${\mathit{\beta }}_x=-\mathrm{\Delta }{E}_{\mathrm{p}}/E_{\mathrm{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., 2021a), in general we may write the energy balance of a particle as

$$\begin{array}{}\text{(34)}& \mathrm{\Delta }{E}_{\mathrm{p}}=-\mathrm{\Delta }{E}_{\mathrm{r}}-f_{\mathrm{c}}-f_y.\end{array}$$

Here, a positive change in rotational energy ΔE_r is seen as an extraction of translational energy. This loss of translational energy with the onset of rotation may be relatively large if a collision involves stick following initial sliding due to a large normal impulse, and such a loss also may occur due to the imposed torque of friction during a collision that does not necessarily involve stick. The term f_c in Eq. (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, due to momentum exchanges associated with the sputtering of loose surface particles during collision. (The videos published as Supplement to DiBiase et al., 2017, nicely illustrate this sputtering as well as the onset of rotational motion.) The term f_y in Eq. (34) represents the energy loss associated with glancing collisions that produce transverse translational motions 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., 2021a). However, the complexity of particle–surface collisions on natural hillslopes precludes explicitly demonstrating such a relation for all possible scenarios. Nonetheless, it is entirely defensible to assume that energy losses can be related to the energy state E_p if the elements involved are formally viewed as random variables. Then, the simple relation $\mathrm{\Delta }{E}_{\mathrm{p}}=-{\mathit{\beta }}_x{E}_{\mathrm{p}}$ is to be viewed as a 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., 2021a) suggest, for example, that $\mathit{\mu }\sim \langle {\mathit{\beta }}_x\rangle \sim M\left(\mathit{\theta }\right)$, 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; Brach and Dunn, 1992, 1995). Moreover, M(θ) is independent of particle size. These points are examined below in relation to experimental data.

In turn, Furbish et al. (2021a) suggest that the quantity α is related to the Kirkby number Ki as

$$\begin{array}{}\text{(35)}& \mathit{\alpha }={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}}{\mathrm{1}-{\mathit{\mu }}_{\mathrm{1}}\mathit{Ki}}},\end{array}$$

where α₀ denotes the value of α associated with a flat surface (Ki=0) and μ₁ is a factor of order unity. Substitution into Eq. (7) gives

$$\begin{array}{}\text{(36)}& L_c={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}{E}_{\mathrm{a}}}{mg\mathit{\mu }\cos \mathit{\theta }-mg{\mathit{\mu }}_{\mathrm{1}}\sin \mathit{\theta }}}.\end{array}$$

Viewed in this manner, α represents a direct effect of heating described by mgμ₁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_{\mathrm{c}}=E_{\mathrm{h}}/mg\mathit{\mu }\cos \mathit{\theta }$. In this view, μ goes with the cooling rate (not the deposition rate). But we also may write Eq. (36) as

$$\begin{array}{}\text{(37)}& L_{\mathrm{c}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}{E}_{\mathrm{h}}}{mg\cos \mathit{\theta }\mathit{\mu }(\mathrm{1}-{\mathit{\mu }}_{\mathrm{1}}\mathit{Ki})}}.\end{array}$$

Viewed in this manner, we may define an apparent friction factor as ${\mathit{\mu }}_{\mathrm{0}}=\mathit{\mu }(\mathrm{1}-{\mathit{\mu }}_{\mathrm{1}}\mathit{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 suggest that for the same slope and surface roughness, α increases with increasing particle size, decreasing the likelihood of frictional loss with increasing rotational energy.

3 Laboratory measurements

3.1 Gabet and Mendoza (2012)

3.1.1 Experiments

The experiments reported by Gabet and Mendoza (2012) involved launching spherical, sub-angular 1 cm particles with initial velocity u₀=0.7 m s⁻¹ onto an inclined surface with fixed roughness elements embedded within concrete (see Fig. 2 therein). The experiments stepped through slope angles of θ=0 to θ=30^∘ in increments of 3^∘. 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^∘. For angles 18, 21, 24 and 27^∘, 92, 58, 26 and 4 particles stopped, respectively. No particles stopped on the 30^∘ slope.

Because the same particle is launched each time with the same initial velocity u₀, the initial arithmetic and harmonic averages of the particle energy are the same, that is, $E_{\mathrm{a}\mathrm{0}}/E_{\mathrm{h}\mathrm{0}}={\mathit{\gamma }}_{\mathrm{0}}\approx \mathrm{1}$. Over some unknown distance the particle motions experience randomization via collisions such that $E_{\mathrm{a}}/E_{\mathrm{h}}=\mathit{\gamma }>\mathrm{1}$. With reference to Eq. (9) where $E_{\mathrm{h}}=E_{\mathrm{a}}/\mathit{\gamma }$, this means that because γ 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_{\mathrm{a}}/E_{\mathrm{h}}=\mathit{\gamma }>\mathrm{1}$ (and fixed) over the entirety of the travel distances. For this reason we truncate distances shorter than the inflection position and then recalculate the travel distances and the exceedance 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 examined further in Appendix B with reference to experiments described in Sect. 3.3.

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 Eqs. (32) and (33). The assumption that γ is fixed may be incorrect. However, rigorously constraining γ would require solving the Fokker–Planck equation describing the evolving number density $n_{E_{\mathrm{p}}}(E_{\mathrm{p}},x)$ of particle energy states E_p as described in Sect. 3.3.2 of the first companion paper (Furbish et al., 2021a); this is beyond the scope of our objective of demonstrating the basic behavior 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.

Figure 3 Plot of exceedance probability R_x(x) versus travel distance x for experiments described by Gabet and Mendoza (2012) where (a) A<0 and (b) A≥1. In (a), slope angles are 0^∘ (top open), 3^∘ (light gray), 6^∘ (bottom open) and 9^∘ (black); in (b), slope angles are 12^∘ (bottom black), 15^∘ (bottom light gray), 18^∘ (open), 21^∘ (top black) and 24^∘ (top light gray).

Figure 4 Plot of logarithm of exceedance probability R_x(x) versus travel distance x for experiments described by Gabet and Mendoza (2012) where (a) A<0 and (b) A≥0, together with fitted distributions (lines). In (a), slope angles are 0^∘ (top open), 3^∘ (light gray), 6^∘ (bottom open) and 9^∘ (black); in (b), slope angles are 12^∘ (bottom black), 15^∘ (bottom light gray), 18^∘ (open), 21^∘ (top black) and 24^∘ (top light gray).

3.1.2 Results

The data are reasonably well fit by the theoretical curves, where we plot the data twice (Figs. 3 and 4) in order to highlight several points. Estimated parametric values are provided in Table 2, where estimates of the variability in μ, α, Ki and 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 Figs. 3 and 4 with coefficients of variation in parentheses.

^a Estimated from parameters A and B. ^b Estimated from data, recognizing that these do not incorporate distances of particles that moved beyond measurement distance of 3 m. ^c Notation (–) means undefined or coefficient of variation is less than 0.01.

For data involving an estimated shape parameter A<0 (Figs. 3a, 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) (Fig. 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 is unambiguous in Fig. 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^∘ and Ki≤0.73. Note that $\mathit{Ki}<{\mathit{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^∘). 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 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 (Sect. 4.2), where vegetation acts as roughness elements.

For data involving an estimated shape factor A≥0 (Figs. 3b, 4b), the first two sets (12 and 15^∘) are approximately exponential with $\mathit{Ki}\approx {\mathit{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) (Fig. 3b) with essentially straight line fits in Fig. 4b. This result is consistent with the idea that these two experiments involved approximately isothermal conditions. For larger shape factor A, the fitted lines decrease in an exponential-like manner in Fig. 3b and they appear as essentially straight lines over the domain of measured travel distances in Fig. 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^∘ span about 90 % of the probability of the distribution, the data for 21^∘ sample only about 50 % of the probability contained in the distribution, and the data for 24^∘ 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.

Further note that for all cases less than 24^∘ 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^∘. Also recall that only four particles stopped on the flume at a slope angle of 27^∘, and no particles stopped at an angle of 30^∘. 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 angles of 27 and 30^∘, respectively.

Figure 5 Plot of reciprocal of average travel distance $\mathrm{1}/{\mathit{\mu }}_x$ versus slope angle θ calculated from Eq. (17) using estimated values of A and B (open circles) and calculated directly from data (black circles).

Here is an interesting sidebar. Following Samson et al. (1999) we plot the reciprocal, $\mathrm{1}/{\mathit{\mu }}_x$, of the average travel distance μ_x versus slope angle θ (Fig. 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 $\mathrm{1}/{\mathit{\mu }}_x$ with increasing slope angle (see their Fig. 3) associated with trapping (deposition) related to collisional friction. This reciprocal then smoothly transitions to values close to zero with further increase in slope angle representing continuing motions without significant trapping. The plot in Fig. 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.

3.2 Kirkby and Statham (1975)

3.2.1 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 Fig. 1 involving a slope angle of 35^∘ 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. Here we briefly examine the effect of the initial energy state determined by the drop height h.

For A and B defined by Eqs. (15) and (16) we may write the average travel distance as

$$\begin{array}{}\text{(38)}& {\mathit{\mu }}_x={\displaystyle \frac{B^{*}}{\mathrm{1}-A}}{\mathit{ϵ}}^{\mathrm{2}}h,\end{array}$$

where $B^{*}=(\mathit{\alpha }/\mathit{\gamma }){\sin }^{\mathrm{2}}\mathit{\theta }/\mathit{\mu }\cos \mathit{\theta }$ or $B=B^{*}{\mathit{ϵ}}^{\mathrm{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 yields 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^*.

3.2.2 Results

The data are well fit by Eq. (38) (Fig. 6). Estimated parametric values are provided in Table 3.

Figure 6 Plot of average travel distance μ_x versus drop height h based on data in Fig. 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 Fig. 6.

We emphasize that other choices of the quantities ϵ, μ and α would provide equally good fits, given that Eq. (38) does not 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); 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.

3.3 Vanderbilt data

3.3.1 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 (Fig. 7). We recorded particle motions using a Lightning RDT monochrome camera (DRS Technologies) operating at 800 frames per second. The camera was mounted on a tripod and oriented parallel to the horizontal surface. The image resolution was 1280 × 640 pixels.

Figure 7 Image of rounded (left), small (center) and angular (right) test particles on the concrete surface with granular roughness texture.

In the second set of experiments we launched particles of varying size and angularity onto the rough surface, then measured their travel distances for several slope angles. The slopes were S= 0, 0.09, 0.15, 0.18, 0.25 and 0.28. The launching device consisted of a pendulum catapult (Fig. 8) configured so that particles were delivered to the rough surface with negligible 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 sixth-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.

3.3.2 Results

As a point of reference, the vertical rebounds of ordinary spherical glass marbles following their impacts on the hard slate reveal no surprises. These collinear collisions give a normal coefficient of restitution of $\mathit{ϵ}=\mathrm{0.81}\pm \mathrm{0.017}$ yielding ${\mathit{\beta }}_z=(\mathrm{1}-{\mathit{ϵ}}^{\mathrm{2}})=\mathrm{0.34}\pm \mathrm{0.028}$ (m=0.014 kg, N=5), and $\mathit{ϵ}=\mathrm{0.81}\pm \mathrm{0.018}$ yielding ${\mathit{\beta }}_z=\mathrm{0.35}\pm \mathrm{0.030}$ (m=0.0033 kg, N=5). The variation in ϵ is likely attributable to small differences in the marble-surface deformation mechanics during collision. Rebounds from the rigid granular surface give $\mathit{ϵ}=\mathrm{0.26}\pm \mathrm{0.59}$ yielding ${\mathit{\beta }}_z=\mathrm{0.93}\pm \mathrm{0.031}$ (m=0.014 kg, N=10), and $\mathit{ϵ}=\mathrm{0.41}\pm \mathrm{0.038}$ yielding ${\mathit{\beta }}_z=\mathrm{0.83}\pm \mathrm{0.032}$ (m=0.030 kg, N=10). Although the collisions are collinear, the granular texture of the surface leads to some variation in the reflection angles. The smaller values of ϵ relative to the slate surface indicate that, although the concrete surface is rigid, its granular texture gives more dissipative collisions, likely involving deformation of micro-asperities. This effect evidently is more pronounced with the larger marbles.

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 ${\mathit{\beta }}_z=\mathrm{1}-{\mathit{ϵ}}_z^{\mathrm{2}}$ is the proportion of the normal (vertical) component of kinetic energy extracted, analogous to β_x. That is, ${\mathit{\beta }}_z=-\mathrm{\Delta }{E}_{\mathrm{p}}/E_{\mathrm{p}}$. This proportion involves a mechanical loss due to particle and surface deformation during the collision, with conversion of energy to heat and sound. But a significant part can go into transverse components of translational energy and rotational energy. Thus, this is a simple demonstration of the idea that β_z (and β_x) must be treated as a random variable rather than a fixed quantity as in the case of the normal coefficient of restitution ϵ associated with spherical particles in a granular gas, although recent efforts have treated this quantity as a random variable (Gunkelmann et al., 2014; Serero et al., 2015).

To illustrate this we write a normalized form of β_z as ${\mathit{\beta }}_z^{*}=({\mathit{\beta }}_z-{\mathit{\beta }}_z^{\min })/(\mathrm{1}-{\mathit{\beta }}_z^{\min })$, where ${\mathit{\beta }}_z^{\min }=\mathrm{1}-{\mathit{ϵ}}^{\mathrm{2}}$ is the minimum value associated with a collinear collision. Although we cannot offer a theoretical basis for the distribution of ${\mathit{\beta }}_z^{*}$, we nonetheless on empirical grounds fit ${\mathit{\beta }}_z^{*}$ to the beta distribution because of its versatile bounded form over [0,1], then transform the fitted distribution back to the original values of β_z. For comparison we fit a Gaussian distribution to the values β_z measured for the spheres.

Figure 9 Plot of (a) cumulative distributions of energy extraction quantity β_z for glass spheres fit to a Gaussian distribution and rounded and angular particles fit to a beta distribution, and (b) associated probability density functions of fitted distributions. Collisions are on hard slate.

Figure 10 Plot of (a) cumulative distributions of the energy extraction quantity β_z for glass spheres fit to a Gaussian distribution and rounded and angular particles fit to a beta distribution, and (b) associated probability density functions of fitted distributions. Collisions are on rough surface.

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 (Fig. 9). For the rough surface, spheres show a larger variance, and there is a stronger redistribution of probability toward β_z=1 with increasing angularity (Fig. 10). We cannot directly map this result 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 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.

Table 4 Average energy partitioning as a proportion of initial energy E_p0= mgh for estimated coefficient of restitution ϵ.

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 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 energy is recovered in vertical motion. In contrast, on average much of the initial energy is dissipated by friction on the rough surface, and only about 10 % or so goes to transverse motion and rotation. About 5 % of the initial energy is recovered in vertical motion. These results are consistent with those reported by Williams and Furbish (2021) involving a larger data set.

Table 5 Slope-parallel launch velocities u₀ measured from high-speed imaging.

High-speed imaging of particle motions launched by the catapult onto the rough surface provides estimates of initial surface-parallel particle velocities u₀ (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.

Figure 11 Plots of exceedance probability versus travel distance for the Vanderbilt experiments over six values of slope S showing angular (open circles), rounded (black circles) and small (gray circles) particles together with fitted distributions (lines).

For all surface slopes, the large rounded particles systematically travel farther than the large angular particles, and the small particles typically travel distances similar to those of the large angular particles (Fig. 11). The data are reasonably fit by the theoretical curves, notably at small and large Ki. At intermediate Ki, several cases involve a systematic deviation from the curves. Estimated parametric values are provided in Table 6, where estimates of the variability in μ, α, Ki and Ki_* are based on a Monte Carlo analysis as described in Appendix  C. As with the data of Gabet and Mendoza (2012), we displace the starting position (x=0) to the first inflections in the raw exceedance probability plots and then recalculate distances and exceedance probabilities. A specific example is provided in Appendix B.

Table 6 Fitted and estimated values of the parameters for the data shown in Fig. 11 with coefficients of variation in parentheses.

^a Estimated from parameters A and B. ^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.

Note that, for reference below, two values of μ, α, Ki and Ki_* are provided. The first value is based on the measured launch velocity u₀ (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 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 α.

For the lower slopes (S=0 to S=0.15), all particles experience rapid thermal collapse (A<0). However, between S=0.15 and S=0.18, the rounded particles appear to transition to a heavy-tailed behavior. By S=0.25, no rounded particles stop on the surface; gravitational heating far exceeds frictional cooling. At this slope the angular and small particles exhibit heavy-tailed behavior with net heating. At S=0.28, only angular particles stopped on the surface with net heating, and these are strongly censored (126 of 210 particles stopped).

For a given slope, values of the friction factor μ for the angular and small particles are systematically larger than values for the rounded particles. This result is consistent with the expectation that angular particles on average experience a greater energy loss during collisions than rounded particles and that larger rounded particles are less likely than are small particles (both rounded and angular) to be influenced by the surface roughness texture during collisions. Whereas values of the factor α generally increase with slope S (and Kirkby number Ki), no pronounced differences between particle size or shape appear.

Video and audio recordings (Supplement, Vanderbilt University Institutional Repository, https://ir.vanderbilt.edu/handle/1803/9742, last access: 9 June 2021) provide clear evidence in support of the results above. In particular the files “Rounded_colinear.avi” and “Angular_colinear.avi” show examples of collinear collisions on the hard slate, with negligible rotation following collision and maximum vertical rebound. 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 Mark Schmeeckle) in that the point of impact involves a particle corner that becomes aligned directly beneath the center of mass at the instant of impact.

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 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 “Angular_rotational_to_vertical.avi” shows conversion of rotational to vertical motion with the second collision. The file “Semiangular_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 at relatively small incident angle. Nonetheless, the strong conversion of translational to rotational motion associated with the former is analogous to the behavior of a particle that experiences stick during a small incident angle collision with conversion of translational to rotational energy (Appendix E in Furbish et al., 2021a).

The files “Angular_18%slope.avi” and “Angular_28%slope.avi” show examples of angular particles launched from the catapult onto the rough surface. Although the surfaces in these videos appear flat because of camera alignment, the slopes are S=0.18 (10.2^∘) and S=0.28 (15.6^∘), so gravitational heating starts immediately. The particles are launched with 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 it continues with heating. For contrast the file “Rounded_0slope.avi” shows an example of a rounded particle that rapidly decelerates and then “dies” when launched onto the flat rough surface (S=0). The increase in free-flight distances (before initial collisions) with increasing slope is 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 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.

4 Field measurements

4.1 DiBiase et al. (2017)

4.1.1 Experiments

The field-based experiments reported by DiBiase et al. (2017) involved launching three different sizes of particles down a natural hillslope surface. The particle size classes included diameters D= 2–3, D= 4–6 and D= 9–12 cm, involving 58, 93 and 43 particles, respectively. Of these, 53, 61 and 14 particles stopped within a 14 m measurement distance with approximately uniform steepness of 38^∘. The distributions of travel distance systematically varied with the different particle sizes. Further details of the experiments, including measurements of surface roughness, are provided by DiBiase et al. (2017).

As with the data of Gabet and Mendoza (2012), we again fit the parameters A and B and then calculate values of μ and α 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, and then we recalculate distances and exceedance probabilities.

Figure 12 Plot of exceedance probability versus travel distance for experiments described by DiBiase et al. (2017) showing small (black circles), medium (gray circles) and large (open circles) particles together with fitted distributions (lines).

4.1.2 Results

The data are reasonably well fit by the theoretical curves (Fig. 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 Fig. 12 with coefficients of variation in parentheses.

^a Estimated from parameters A and B. ^b Estimated from data, recognizing that these do not incorporate distances of particles that moved beyond measurement distance of 14 m. ^c Notation (–) means coefficient of variation is less than 0.01.

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 ($\mathrm{1}/\mathrm{2}\le A<\mathrm{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 (Fig. 4b), the concavity in the semi-log exceedance probability plot (Fig. 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. 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 $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$, $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{8.5}\times {\mathrm{10}}^{-\mathrm{4}}$ J. Estimated values of the average kinetic energies E_a measured over the total travel distances are $\mathrm{8.0}\times {\mathrm{10}}^{-\mathrm{6}}$, $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{3}}$ J. The changes in average energies ΔE_a are $-\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{6}}$, 0.0×10⁰ and $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{3}}$ J. These changes are qualitatively consistent with net cooling, isothermal conditions and net heating, although in all three cases the estimated parametric values suggest net heating ($\mathit{Ki}\ge {\mathit{Ki}}_{*}$). Using the largest value of the friction factor μ (Table 7), the Kirkby number of the 40^∘ surface immediately downslope from the measurement site is Ki≈1.1. The transition Kirkby number ${\mathit{Ki}}_{*}=\mathrm{1}$. If particles on average experienced net heating on the upper 38^∘ surface, then this result is consistent with the reported observation that particles reaching the steeper slope continued to the base of the hillslope without stopping.

4.2 Roth et al. (2020)

4.2.1 Experiments

The field-based experiments reported by Roth et al. (2020) involved dropping three different sizes of particles on eight natural hillslope surfaces in the Oregon Coast Range. Five of the hillslopes were covered with natural vegetation (designated by V) and included slope angles of 0, 14, 20, 24 and 39^∘. Three of the hillslopes had recently been burned (designated by B) and included slope angles of 17, 20 and 28^∘. Particle size classes involved average diameters of 1.7, 4.5 and 7.3 cm. These were dropped from a height of h≈0.2 m onto each hillslope surface. The distributions of travel distances systematically varied with slope angle, particle size and surface roughness conditions.

The surfaces of the vegetated hillslopes had a layer of duff, woody debris and banana slugs beneath small plants (e.g., ferns) and trees. The surfaces of the burned hillslopes had little vegetal cover and were markedly smoother than the vegetated sites. Further details of the experiments, including measurements of surface roughness, are provided by Roth et al. (2020). Banana slugs, whose locomotive energetic costs are constrained by the shear-thinning rheology of their mucus (Lauga and 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 and 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, and then we 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 meters 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 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, and then use the average downslope reflection 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.

Figure 13 Plot of logarithm of exceedance probability R_x(x) versus travel distance x for experiments described by Roth et al. (2020). These examples are for sites V14 (a) and V24 (b) showing data for small (black circles), medium (gray circles) and large (open circles) particle sizes, together with fitted distributions (lines).

Table 8 Fitted and estimated values of the parameters for the data reported by Roth et al. (2020) as shown in Figs. 13 and 14 with coefficients of variation in parentheses.

^a Estimated from parameters A and B. ^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.

Figure 14 Plot of logarithm of exceedance probability R_x(x) versus travel distance x for experiments described by Roth et al. (2020). These examples are for sites B17 (a) and B20 (b) showing data for small (black circles), medium (gray circles) and large (open circles) particle sizes, together with fitted distributions (lines).

4.2.2 Results

The data for the V sites are reasonably well fit by the theoretical curves (Fig. 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 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 (Fig. 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), 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.

Figure 15 Plot of exceedance probability versus travel distance for experiment B28M described by Roth et al. (2020) showing data (circles) fit to mixed distribution (black line) composed of sum of two distributions (gray lines) weighted by proportions p₁ and $p_{\mathrm{2}}=\mathrm{1}-p_{\mathrm{1}}$. Note the effect of increased frictional cooling after slope inflection at ∼33 m; data at x=50 m are censored but included for reference. Plots for B28S and B28L (Table 9) are similar in appearance.

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 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 (Fig. 15) – in effect giving a mixed distribution. Namely, let x₁ and x₂ denote the travel distances of the two groups, and let p₁ denote the proportion represented by first group such that $p_{\mathrm{2}}=\mathrm{1}-p_{\mathrm{1}}$ is the proportion of the second group. The simplest form of a mixed distribution is

$$\begin{array}{}\text{(39)}& f_x\left(x\right)=p_{\mathrm{1}}{f}_{x_{\mathrm{1}}}\left(x_{\mathrm{1}}\right)+(\mathrm{1}-p_{\mathrm{1}})f_{x_{\mathrm{2}}}\left(x_{\mathrm{2}}\right),\end{array}$$

and the cumulative distribution is

$$\begin{array}{}\text{(40)}& F_x\left(x\right)=p_{\mathrm{1}}\underset{\mathrm{0}}{\overset{x}{\int }}{f}_{x_{\mathrm{1}}}\left(x_{\mathrm{1}}^{\prime }\right)\mathrm{d}{x}_{\mathrm{1}}^{\prime }+(\mathrm{1}-p_{\mathrm{1}})\underset{\mathrm{0}}{\overset{x}{\int }}{f}_{x_{\mathrm{2}}}\left(x_{\mathrm{2}}^{\prime }\right)\mathrm{d}{x}_{\mathrm{2}}^{\prime },\end{array}$$

where primes denote variables of integration. As above, the exceedance probability is $R_x\left(x\right)=\mathrm{1}-F_x\left(x\right)$.

Table 9 Fitted and estimated values of the parameters for the data shown in Fig. 15.

For the three particle sizes, exceedance probabilities are well matched by the weighted sum of a nearly exponential distribution (A₁≈0) reflecting isothermal conditions and a heavy-tailed distribution (A₂>0) reflecting net heating (Table 9). Note, however, that the parametric values A₁, B₁, A₂ and B₂ cannot be combined to estimate associated factors such as μ and α. 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 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 (Fig. 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 (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.

5 Analysis

We emphasize at the outset a key point in comparing travel distances measured in experiments (laboratory or field) with theoretical distributions. By definition a sample of measured values drawn from a distribution possesses a finite sample mean and variance, regardless of the form of the underlying distribution. If the underlying distribution possesses a finite mean and variance (e.g., an exponential distribution), then the calculated sample average and variance are unbiased estimates of the underlying parametric values. If the mean or variance of the underlying distribution is undefined (e.g., the generalized Pareto distribution for $A\ge \mathrm{1}/\mathrm{2}$), then the calculated sample average and variance have no meaningful relation to the underlying (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, quantile–quantile 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 that the calculated values do not meaningfully characterize a property (e.g., absence of central tendencies) of the underlying distribution.

In contrast, estimates of the parameters A and B are less sensitive to this uncertainty when these values are used to calculate moments (if they exist) – but only if the selected form of the distribution is the correct choice. The mechanical basis of the generalized Pareto distribution lends confidence, but does not guarantee, that it is the correct choice. Moreover, uncertainty in the estimated values of A and B increases when a decreasing proportion of the tail of the distribution is sampled. Indeed, one can never know the form of the censored tail (Ballio et al., 2019).

Figure 16 Plot of modified exceedance probability R_* versus dimensionless travel distance x_* and line with log–log slope of −1 for (a) experiments described by Gabet and Mendoza (2012) (gray circles) and experiments described by DiBiase et al. (2017) (open circles), (b) Vanderbilt experiments, and (c) experiments described by Roth et al. (2020) for V sites (gray circles) and B sites (open circles). Data for A<0 fall to the left of $x_{*}=\mathrm{1}$ with values in the tails represented by smaller values of x_*. Data for A>0 fall to the right of $x_{*}=\mathrm{1}$ with values in the tails represented by larger values of x_*. Total data numbers are (a) N=813, (b) N=2980 and (c) N=1878.

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:

$$\begin{array}{}\text{(41)}& R_{*}=R_x^A\mathrm{and}{x}_{*}={\displaystyle \frac{A}{B}}x+\mathrm{1}.\end{array}$$

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 (Fig. 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 3 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 variety of slope and surface roughness conditions.

The bounded form of the generalized Pareto distribution (A<0) must not be viewed as involving a “hard” boundary. Because of the stochasticity of motions associated with varying sizes and shapes, some particles by chance “leak” beyond the position $x=B/\left|A\right|$. This aspect of the formulation is necessarily simplified. What is clear, nonetheless, is the rapid thermal collapse reflected by the (approximately) bounded form of the distribution in the laboratory measurements of Gabet and Mendoza (2012) (Figs. 3, 4) and the measurements at Vanderbilt (Fig. 11), as well as the field-based measurements of Roth et al. (2020) (Figs. 13, 14).

From an empirical point of view the data are consistent with the generalized Pareto distribution and reflect the predicted variation in behavior from rapid thermal collapse to approximately isothermal conditions to net heating of particles. Nonetheless we proceed by asking whether the estimated values of the quantities μ and α make physical sense, while recognizing that these quantities do not readily map to established formulations of friction, and represent a complexity that cannot be encapsulated in idealized collision mechanics (Appendix E in Furbish et al., 2021a).

The laboratory measurements with zero slope merit particular attention. The Kirkby number Ki is known and zero. The effect of heating, and thus the influence of heating on α, is removed. Focusing on the length scale L_c given by Eq. (37), motions are mass independent and the initial velocity is approximately fixed. Assuming fixed γ and comparing the angular and rounded particles in the Vanderbilt data (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 (Figs. 9, 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 particles; 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 remain 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 (Fig. 17). Note that the Vanderbilt data in this figure are based on the reduced velocity calculations (Table 6). For completeness this figure is reproduced in Appendix E using the initial launch velocities u₀ (Table 5).

Figure 17 Plot of friction factor μ versus slope S for laboratory experiments described by Gabet and Mendoza (2012) (black circles) and Vanderbilt data (open squares), and field-based experiments of DiBiase et al. (2017) (gray circles) and Roth et al. (2020) for V sites (open circles) and B sites (open triangles) together with 1:1 line.

Values of μ for Ki<1 systematically fall above the 1:1 line (Fig. 17) and 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). 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 Sect. 2.4, scaling suggests that the factor μ∼M(θ) is independent of particle size (Furbish et al., 2021a). 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 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., 2021a). However, this probably is insufficient to explain the relation in Fig. 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 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. (33) versus the Kirkby number Ki (Fig. 18) together with the function in Eq. (35). Note that the Vanderbilt data in this figure are based on the reduced velocity calculations (Table 6). For completeness this figure is reproduced in Appendix E using the initial launch velocities u₀ (Table 5).

Figure 18 Plot of factor α versus Kirkby number Ki for experiments described by Gabet and Mendoza (2012) (black circles), Vanderbilt data (open squares), DiBiase et al. (2017) (gray circles) and Roth et al. (2020) for V sites (open circles) and B sites (open triangles) together with function $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{0}}/(\mathrm{1}-{\mathit{\mu }}_{\mathrm{1}}\mathit{Ki})$ with α₀=1 and μ₁=0.98 (line). Only data for A<1 are included.

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) 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, Fig. 6) rather than being estimated from A and B^*, we do not plot these values in Fig. 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

$$\begin{array}{}\text{(42)}& L_{\mathrm{c}}={\displaystyle \frac{\mathit{\alpha }{E}_{\mathrm{a}}}{\mathit{\gamma }mg\mathit{\mu }\cos \mathit{\theta }}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}{E}_{\mathrm{a}}}{\mathit{\gamma }mg\mathit{\mu }\cos \mathit{\theta }(\mathrm{1}-{\mathit{\mu }}_{\mathrm{1}}\mathit{Ki})}}.\end{array}$$

With $E_{\mathrm{a}}=(m/\mathrm{2})\langle u^{\mathrm{2}}\rangle$, the effect of particle mass m does not explicitly appear. This means that any effect of particle size must 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 (Fig. 18). Although not explicit in the formulation, we suspect that the effect 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 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 (Fig. 11a). When pooled these data can be approximately fitted (not shown) to a generalized Pareto distribution. However, the data are well fit using the mixed distribution defined by Eq. (39) (Fig. 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.

Figure 19 Plot of exceedance probability versus travel distance for Vanderbilt University experiment with S=0 showing data (circles) fit to mixed distribution (M) composed of sum of distributions of angular particles (A) and rounded particles (R) depicted in Fig. 11.

6 Discussion and conclusions

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 point made in the first companion paper (Furbish et al., 2021a). 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; it 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.

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; 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 surface-parallel motions, a tumbling angular particle is more likely than is a rounded tumbling particle to experience a noncollinear 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 (Fig. 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 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., 1998, 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 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 terms of particle properties and surface-roughness conditions, as well as 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 parameter 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 Sect. 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 process, specifically the likelihood of converting translational to rotational energy and the decreasing extraction of energy by collisional friction (Williams and Furbish, 2021).

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 soon after launch from the catapult as particle motions become randomized with the onset of particle–surface collisions, so there is uncertainty in choosing the truncation distance and the associated particle energy state. Similarly, little is known about the distribution $f_{E_{\mathrm{p}}}\left(E_{\mathrm{p}}\right)$ 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 (Fig. 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 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 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 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 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 and possibly for fresh ideas concerning particle motions more generally.

Appendix A: Parameter estimation

Here we demonstrate the basis for using visual fits of the exceedance probability plots to illustrate the behavior of the generalized Pareto distribution, and we provide context for interpreting these fits. We work with the dimensionless form of the distribution for comparison with Fig. 2 and pursue a straightforward Monte Carlo analysis.

First, let ${\widehat{x}}_{\mathrm{u}}$ denote a random number drawn from a uniform distribution with support [0,1] and cumulative distribution function $F_{{\widehat{x}}_{\mathrm{u}}}\left({\widehat{x}}_{\mathrm{u}}\right)={\widehat{x}}_{\mathrm{u}}$. In turn, the cumulative distribution function of the generalized Pareto distribution is

$$\begin{array}{}\text{(A1)}& F_{\widehat{x}}\left(\widehat{x}\right)=\mathrm{1}-{\displaystyle \frac{b^{\mathrm{1}/a}}{(a\widehat{x}+b{)}^{\mathrm{1}/a}}}.\end{array}$$

Equating $F_{{\widehat{x}}_{\mathrm{u}}}\left({\widehat{x}}_{\mathrm{u}}\right)$ and Eq. (A1) leads to

$$\begin{array}{}\text{(A2)}& \widehat{x}={\displaystyle \frac{b}{a}}\left[{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-{\widehat{x}}_{\mathrm{u}}{)}^a}}-\mathrm{1}\right],\end{array}$$

which provides an algorithm for generating values of $\widehat{x}$ drawn from the generalized Pareto distribution with shape and scale parameters a and b, starting with values of ${\widehat{x}}_{\mathrm{u}}$ selected by a uniform random-number generator.

Second, among the methods for estimating the values of a and b are the method of moments and the maximum likelihood estimation (MLE) method. Both are unsuitable for censored data, and the method of moments is unsuitable when either the first or second moment is undefined. Nonetheless, for the purpose of this appendix we focus on the MLE method, noting that an MLE is the same as a Bayesian estimate assuming a uniform (maximum entropy) prior distribution of the parametric values. The MLE method is a popular, standard choice for estimating parametric values of distributions, notably heavy-tailed distributions, because of its asymptotic properties of consistency and efficiency. (However, see the delightful review by Cam (1990) regarding maximum likelihood estimates, in particular his nine principles on p. 165.) The shape and scale parameters a and b are not orthogonal. The MLE of a is $\stackrel{\mathrm{̃}}{a}=\mathrm{1}/{\stackrel{\mathrm{̃}}{a}}_{\mathrm{L}}$ where

$$\begin{array}{}\text{(A3)}& {\stackrel{\mathrm{̃}}{a}}_{\mathrm{L}}={\displaystyle \frac{N}{{\sum }_{i=\mathrm{1}}^N\ln (\mathrm{1}+{\widehat{x}}_i/{\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}})}},\end{array}$$

and the MLE estimate of b is $\stackrel{\mathrm{̃}}{b}={\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}}/{\stackrel{\mathrm{̃}}{a}}_{\mathrm{L}}$ where ${\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}}$ is obtained from an iterative solution of

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\displaystyle \frac{N}{{\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}}{\sum }_{i=\mathrm{1}}^N{\widehat{x}}_i/({\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}}^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}}{\widehat{x}}_i)}}\\ & -{\displaystyle \frac{N}{{\sum }_{i=\mathrm{1}}^N\ln (\mathrm{1}+{\widehat{x}}_i/{\stackrel{\mathrm{̃}}{b}}_{\mathrm{L}})}}-\mathrm{1}=\mathrm{0}.\end{array}\end{array}$$

These are biased estimates (Giles et al., 2013a), 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 ${\mathit{\mu }}_{\widehat{x}}=b$. The MLE of the mean ${\mathit{\mu }}_{\widehat{x}}$ of an exponential distribution is just the sample average,

$$\begin{array}{}\text{(A5)}& \stackrel{\mathrm{̃}}{b}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\widehat{x}}_i,\end{array}$$

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 and then calculate and plot exceedance probabilities for varying values of the shape parameter a, holding the scale parameter b fixed for convenience. We also calculate the MLE of a using Eq. (A3) for each sample. The MATLAB/GNU Octave code for doing this is available in the Supplement (Vanderbilt University Institutional Repository, https://ir.vanderbilt.edu/handle/1803/9742, last access: 9 June 2021) and includes an animation of the results.

Plots of estimated exceedance probabilities $R_{\widehat{x}}\left(\widehat{x}\right)$ for all samples provide a visual sense of the variability in these probabilities associated with the inherent randomness in drawing samples of $\widehat{x}$ from the known distribution (Fig. A1). The animation mentioned above shows the outcome of successive draws and nicely illustrates that many draws, by chance, bear little resemblance 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 variability in the MLE of a increases with increasing a (Fig. 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=-\mathrm{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; 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 – that is, these values converge to the true values – only when N approaches 10 000 or more (Fig. A3). Fourth, with increasing censorship of the data, the MLEs based on the uncensored values become strongly biased (Fig. 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 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 (Figs. A1, A2), we calculate sample exceedance probabilities and $\stackrel{\mathrm{̃}}{b}$ for the exponential distribution. According to the central limit theorem, values of $\stackrel{\mathrm{̃}}{b}$ are approximately normally distributed with variance $\sim {\mathit{\sigma }}_{\widehat{x}}^{\mathrm{2}}/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 so-called confidence estimates associated with these known values of a and b. These specifically give information regarding the 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 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 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).

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 Eqs. (A3) and (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) estimates for $A\mathit{\lesssim }-\mathrm{1}/\mathrm{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 (Fig. 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, 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., 2013a, b; Pak 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 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 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 scaling relations may be particularly sensitive to the presence or absence of extreme values in the data sets used to estimate the exponents. To quote the sixth of nine delightful principles offered by 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.

Figure A1 Plots of exceedance probability $R_{\widehat{x}}\left(\widehat{x}\right)$ versus dimensionless travel distance $\widehat{x}$ for different values of the shape parameter a assuming scale parameter b=1. Each plot shows 1000 samples, each of size n=100, together with theoretical exceedance probability function (line).

Figure A2 Histograms of maximum likelihood estimates of shape parameter a assuming scale parameter b=1, and maximum likelihood estimate of scale parameter b for a=0. Each histogram is based on 10 000 samples.

Figure A3 Plots of exceedance probability $R_{\widehat{x}}\left(\widehat{x}\right)$ versus dimensionless travel distance $\widehat{x}$ for shape parameter a=0.5 assuming scale parameter b=1, showing convergence to theoretical exceedance probability function (lines) with increasing sample size N. Examples involve (a) 100 samples each of size N=1000 and (b) 20 samples each of size N=10 000.

Figure A4 Histograms of maximum likelihood estimates of shape parameter a=1 assuming scale parameter b=1 with censorship at $\widehat{x}=\mathrm{50}$ based on 10 000 samples. The bias increases as the censorship distance decreases. Compare with Fig. A2.

Table A1 Fitted and estimated values of the parameters for the data reported by Gabet and Mendoza (2012), the Vanderbilt experiments, DiBiase et al. (2017), and Roth et al. (2020).

¹ Estimated visually; values reproduced from Tables 2, 6, 7 and 8. ² Estimated from MLE algorithm; asterisk denotes problematic estimate due to A≲−0.5, A→0 or censored data. ³ A – large angular, R – large rounded, S – small; V – vegetated, B – burned.

Figure A5 Example of fit (gray line) based on MLE values of A and B versus visual fit (black line). This example coincides with the smallest particle size reported by DiBiase et al. (2017) (Table 7).

Appendix B: Particle launching conditions

Consider as an example the initial exceedance probability plot for the angular particles on a flat surface (Fig. B1a), which shows a clear inflection at about 5 cm. In this example, high-speed imaging of particles launched from the catapult reveals 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 an initial disentrainment rate that is smaller than after randomization. However, we note that the inflection also may involve other effects.

Figure B1 Plot of exceedance probability R_x(x) versus travel distance x for the example of angular particles on a flat surface showing (a) initial data set and (b) truncated data set with fitted distribution (line).

For these reasons we truncate the plots at the inflection position and then recalculate exceedance probabilities with reduced N (Fig. 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 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/\sqrt{N}$ and $B/\sqrt{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 $\sim \mathrm{1}/\sqrt{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 (Sect. 3.3) we also incorporate the uncertainty in the launch velocities u₀ provided in Table 5.

We perform a straightforward Mont Carlo analysis. Assuming values of A, B and u₀ are approximately Gaussian, we successively solve Eqs. (32), (33), (8) and (23) 10 000 times and then calculate the associated coefficients of variation of μ, α, Ki and Ki_*. Because A and B are obtained by visual fitting (as opposed to being based on MLE values) where the number of censored data are included in calculations of exceedance probabilities, we include censored data in setting N.

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 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 order of 10 % or more for smaller slopes in the experiments reported by Gabet and Mendoza (2012) and in the Vanderbilt experiments. Coefficients of variation generally are on the order of 1 % or smaller in the field-based experiments of DeBiase et al. (2017) and Roth et al. (2020) involving steep slopes.

Appendix D: Particle energy balance

Let $E_{\mathrm{p}\mathrm{0}}=mgh=(\mathrm{1}/\mathrm{2})m{w}_{\mathrm{0}}^{\mathrm{2}}$ denote the initial impact energy of a particle with mass m falling from height h onto a horizontal surface, where w₀ is the vertical impact velocity. Then let w₁ and u₁ denote the vertical and horizontal rebound velocity components. Assuming negligible rotational energy during the initial free fall, the energy balance may be written as

$$\begin{array}{}\text{(D1)}& E_{\mathrm{p}\mathrm{0}}=f_{\mathrm{c}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m{u}_{\mathrm{1}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m{w}_{\mathrm{1}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}I{\mathit{\omega }}^{\mathrm{2}},\end{array}$$

where f_c is the frictional loss due to particle–surface deformation, I is the moment of inertia and ω is the angular velocity. If we set the initial and final vertical positions of the rebounding motion at $z\left(\mathrm{0}\right)=z\left(T\right)=\mathrm{0}$, then from Newton's second law,

$$\begin{array}{}\text{(D2)}& w_{\mathrm{1}}={\displaystyle \frac{g}{\mathrm{2}}}T,\end{array}$$

where T is the travel time to the second collision. We assume as an approximation that $f_c=(\mathrm{1}-{\mathit{ϵ}}^{\mathrm{2}})E_{\mathrm{p}\mathrm{0}}$, 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

$$\begin{array}{}\text{(D3)}& {\mathit{ϵ}}^{\mathrm{2}}{E}_{\mathrm{p}\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}m{g}^{\mathrm{2}}{T}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m{u}_{\mathrm{1}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}I{\mathit{\omega }}^{\mathrm{2}}.\end{array}$$

We can experimentally determine ϵ, E_p0 and T for individual particles. However, we generally cannot determine u₁ from side 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,

$$\begin{array}{}\text{(D4)}& E_{\mathrm{c}}={\mathit{ϵ}}^{\mathrm{2}}{E}_{\mathrm{p}\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}m{g}^{\mathrm{2}}{T}^{\mathrm{2}}.\end{array}$$

Appendix E: Effect of launch velocity

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₀ rather than a reduced velocity in the calculations as presented in Figs. 17 and 18.

Values of μ, notable at smaller slopes S (Fig. E1), are noticeably larger than those calculated in Fig. 17. The values appear to converge to the 1:1 line with increasing slope as Ki→1, as in Fig 17.

Figure E1 Plot of friction factor μ versus slope S for data shown in Fig. 17 using the initial launch velocity u₀ in the calculations.

Values of α and associated Kirkby numbers tend to be smaller (Fig. E2) than those calculated in Fig. 18. However, the overall variation between α and Ki is similar.

Figure E2 Plot of factor α versus Kirkby number Ki for data shown in Fig. 18 using the initial launch velocity u₀ in the calculations.

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=\mathrm{1}/{\mathit{\mu }}_x$. We then write the conditional distribution as

$$\begin{array}{}\text{(F1)}& f_{x|P_x}\left(x\right|P_x)=P_x{e}^{-P_{x}x}.\end{array}$$

We may now treat the rate P_x as a random variable that is distributed as a gamma distribution, namely,

$$\begin{array}{}\text{(F2)}& f_{P_x}(P_x;A_{\mathrm{g}},B_{\mathrm{g}})={\displaystyle \frac{B_{\mathrm{g}}^{A_{\mathrm{g}}}}{\mathrm{\Gamma }\left(A_{\mathrm{g}}\right)}}{P}_x^{A_{\mathrm{g}}-\mathrm{1}}{e}^{-B_{\mathrm{g}}{P}_x},\end{array}$$

with shape parameter A_g and scale parameter B_g. The unconditional distribution of travel distances x is then obtained as a gamma weighting of the conditional exponential distribution, namely,

$$\begin{array}{}\text{(F3)}& f_x\left(x\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_{x|P_x}\left(x\right|P_x\left)f_{P_x}\right(P_x;A_{\mathrm{g}},B_{\mathrm{g}})\mathrm{d}{P}_x.\end{array}$$

Substituting Eqs. (F1) and (F2) into Eq. (F3) and evaluating the integral then leads to

$$\begin{array}{}\text{(F4)}& f_x\left(x\right)={\displaystyle \frac{A_{\mathrm{g}}{B}_{\mathrm{g}}^{A_{\mathrm{g}}}}{(x+B_{\mathrm{g}}{)}^{\mathrm{1}+A_{\mathrm{g}}}}}.\end{array}$$

This is a Lomax distribution (compare with Eq., 28) with shape parameter $A_{\mathrm{g}}=\mathrm{1}/A$ and scale parameter $B_{\mathrm{g}}=B/A=B{A}_{\mathrm{g}}$.

The expected value ${\mathit{\mu }}_{P_x}=A_{\mathrm{g}}/B_{\mathrm{g}}=\mathrm{1}/B$ and the variance is ${\mathit{\sigma }}_{P_x}^{\mathrm{2}}=A_{\mathrm{g}}/B_{\mathrm{g}}^{\mathrm{2}}=A/B^{\mathrm{2}}$. This immediately implies that an experimental estimate of B provides an estimate of the expected disentrainment rate $\mathrm{E}\left(P_x\right)={\mathit{\mu }}_{P_x}$, 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

$$\begin{array}{}\text{(F5)}& P_x={\displaystyle \frac{\mathit{\gamma }mg\mathit{\mu }\cos \mathit{\theta }}{\mathit{\alpha }{E}_{\mathrm{a}\mathrm{0}}}},\end{array}$$

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_{{\mathrm{a}}_{\mathrm{0}}}$, 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_{P_x}(P_x;\mathit{\alpha },\mathit{\beta })\mathrm{d}{P}_x$ is the probability that particles possess the value P_x. These particles are deposited exponentially with mean ${\mathit{\mu }}_x=\mathrm{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.

Figure F1 Plot of (a) probability density $f_{P_x}(P_x;A_{\mathrm{g}},B_{\mathrm{g}})$ of disentrainment rate P_x and (b) probability density of mean travel distance ${\mathit{\mu }}_x=\mathrm{1}/P_x$ for $A=\mathrm{0.01},\mathrm{0.5},\mathrm{0.95}$ $(A_{\mathrm{g}}=\mathrm{100},\mathrm{2},\mathrm{1.05})$ with B=1.

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 $\mathrm{1}/B$ (Fig. E1a) and approaches a Dirac function in the limit of A→0 as the variance ${\mathit{\sigma }}_{P_x}^{\mathrm{2}}=A/B^{\mathrm{2}}\to \mathrm{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 distribution of the reciprocal ${\mathit{\mu }}_x=\mathrm{1}/P_x$ is given by the inverse gamma distribution (Fig. 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 surface roughness – then Eq. (F5) suggests the possibility of formulating a multiplicative relation between these properties and the shape parameter $B=\mathrm{1}/\mathrm{E}\left(P_x\right)$. An initial effort to this effect is reported by Roth et al. (2020).