Research article 26 Mar 2021
Research article  26 Mar 2021
Particle size dynamics in abrading pebble populations
 ^{1}MTABME Morphodynamics Research Group, Budapest University of Technology and Economics, Műegyetem rakpart 1–3, Budapest, Hungary
 ^{2}Department of Mechanics, Materials and Structures, Budapest University of Technology and Economics, Műegyetem rakpart 1–3, Budapest, Hungary
 ^{3}Department of Theoretical Physics, Budapest University of Technology and Economics, Budafoki út 8, Budapest, Hungary
 ^{1}MTABME Morphodynamics Research Group, Budapest University of Technology and Economics, Műegyetem rakpart 1–3, Budapest, Hungary
 ^{2}Department of Mechanics, Materials and Structures, Budapest University of Technology and Economics, Műegyetem rakpart 1–3, Budapest, Hungary
 ^{3}Department of Theoretical Physics, Budapest University of Technology and Economics, Budafoki út 8, Budapest, Hungary
Correspondence: András A. Sipos (siposa@eik.bme.hu)
Hide author detailsCorrespondence: András A. Sipos (siposa@eik.bme.hu)
Abrasion of sedimentary particles in fluvial and eolian environments is widely associated with collisions encountered by the particle. Although the physics of abrasion is complex, purely geometric models recover the course of mass and shape evolution of individual particles in low and middleenergy environments (in the absence of fragmentation) remarkably well. In this paper, we introduce the first model for the collisiondriven collective mass evolution of sedimentary particles. The model utilizes results of the individual, geometric abrasion theory as a collision kernel; following techniques adopted in the statistical theory of coagulation and fragmentation, the corresponding Fokker–Planck equation is derived. Our model uncovers a startling fundamental feature of collective particle size dynamics: collisional abrasion may, depending on the energy level, either focus size distributions, thus enhancing the effects of sizeselective transport, or it may act in the opposite direction by dispersing the distribution.
1.1 Geological observations
Probably the most fundamental observation on pebbles is that they appear to be segregated both by size and shape, and it is broadly accepted that the dynamics are driven by two physical processes: transport and abrasion. Which of these processes dominates may depend on the geological location and also on timescales; however, geologists appear to agree that, in general, neither process should be ignored.
In coastal environments, one of the most remarkable accounts of pebble size and shape distribution is provided by Carr (1969) based on the measurement of approximately 100 000 pebbles on Chesil Beach, Dorset, England. In summarizing his results, Carr provides mean values and sample variations for maximal pebble size and pebble axis ratios along lines orthogonal to the beach. These plots reveal pronounced segregation by maximal size and shape; i.e., on shingle beaches pebbles of roughly similar maximal sizes and with roughly similar axis ratios appear to be spatially close to each other. Size and shape segregation has been broadly observed in various settings (Bird, 1996; Gleason et al., 1975; Hansom and Moore, 1981; Kuenen and Migliorini, 1950; Neate, 1967), and it was mostly attributed to the global transport of pebbles by waves (Lewis, 1931; Carr, 1969) but, in some settings, may also be related to abrasion. Indeed, a detailed account of the interaction of abrasion and transport is given by Landon (1930), who investigated the beaches on the west shore of Lake Michigan. He attributes size and shape variation to a mixture of abrasion and transport. Kuenen (1964) discusses Landon's observations but disagrees with the conclusions and attributes size and shape variation primarily to transport. Carr (1969) observes dominant sizes and shape ratios emerging as a result of abrasion and size grading, while Bluck (1967) describes beaches in South Wales where equilibrium distributions of size and shape are reached primarily by transport and abrasion plays a minor role. Which of the two processes (transport or abrasion) dominate may well depend on the timescales they operate on. While abrasion appears in some scenarios to act much more slowly than transport, a recent study (Bertoni et al., 2016) verified mass losses on the order of 50 % on a pebble beach over a 13month period, indicating that in some settings the two processes may indeed compete in determining size and shape distributions.
In fluvial environments, while downstream fining of sediment has been often attributed to transport (Paola et al., 1992; Ferguson et al., 1996; Fedele and Paola, 2007; Whittaker et al., 2011), other authors have pointed to the significance of attrition (Brewer and Lewin, 1993; Attal and Lavé, 2006; Dingle et al., 2017). In Miller et al. (2014) the authors, using field data, provide quantitative assessment of the significance of selective transport with respect to attrition in downstream fining. Beyond the evolution of smooth size and shape distributions, there is yet another common phenomenon in fluvial geomorphology where the interaction of transport and attrition could be far from trivial. The often observed presence of isolated large boulders in rivers (Huber et al., 2020) may be explained solely by transport, as these large pieces are often not carried by the river; rather, they move by a different process (e.g., landslide or debris flow). On the other hand, these large rocks could also be interpreted as outliers emerging spontaneously in a pebble size distribution on which collisional abrasion certainly has strong impact in upper reaches of rivers.
As we can see, both in coastal and fluvial environments it is a generally accepted fact that the two processes (transport and attrition) appear to compete in shaping the evolution of pebble shape and pebble mass distributions. How exactly this competition may play out and in what manner attrition may contribute to this process is the subject of our paper.
We also remark that while all available observations indicate that attrition could be a relevant factor in the evolution of shape and mass distributions, so far, in the absence of any predictive theory, no datasets have been collected which would allow verifying any theoretical predictions. We will point out potential strategies for verification in Sect. 4.
1.2 Existing theory
1.2.1 Individual abrasion
Individual abrasion is a theory describing the mass and shape evolution of one individual particle (abraded particle) under the impacts of many incoming particles (abraders) (see Fig. 1a). In the mean field theory for the geometry of individual abrasion only the mass and shape of the abraded particle is recorded; the effect of impacts is averaged and the evolution is determined by the size of the abraded particle compared to the average size of the abrading particles.
Since the seminal papers by Firey (1974) and Bloore (1977), the mean field geometric theory of individual abrasion (i.e., shape evolution) for sedimentary particles under collisions has been well understood and validated (Szabó et al., 2013, 2015; NovákSzabó et al., 2018). Still, despite the success of the Firey–Bloore geometric theory of shape evolution, it is clear (Domokos and Gibbons, 2018) that it is not suited to predict the evolution of size: in stark contrast with geological observations summarized in Sternberg's law (Sternberg, 1875), predicting exponential decay of particle mass and an infinite lifetime for all particles, geometric abrasion theory predicted a finite lifetime for all particles. On the other hand, Sternberg's broadly accepted theory of mass evolution (Sternberg, 1875) had nothing to offer regarding the evolution of shape. Recognizing this challenge, in Domokos and Gibbons (2018) a unified theory, called volumeweighted shape evolution, has been proposed which, on one hand, reproduces all the geometric features of the Firey–Bloore geometric theory and, on the other hand, also predicts mass evolution in accordance with Sternberg's law.
1.2.2 Binary abrasion
The first stepping stone between the theory of individual abrasion and collective abrasion is the model for mutual, binary abrasion, where two particles mutually abrade each other, and we track both evolutions (see Fig. 1b). In this case one can still write mean field equations by averaging over many collisions, and the mass and shape evolution of both particles are recorded. For any binary abrasion model of size evolution, we postulate the following requirements:

size evolution should follow Sternberg's law,

mass loss in a collision should be a monotonically increasing function of collision energy, and

the model should be fully compatible with the geometric evolution model.
The unified theory in Domokos and Gibbons (2018) offers a model satisfying all three requirements: by extending the Firey–Bloore equations and Sternberg's theory and using the kinetic energy of collision, models for binary shape evolution and for binary mass evolution of two mutually abrading particles were put forward. These two models have been merged in Domokos and Gibbons (2018) into a unified volumeweighted theory of binary abrasion, compatible both with the Firey–Bloore and with the Sternberg theory. The volumeweighted model for binary mass evolution, describing the time evolutions for the masses X(t), Y(t) of two particles with respective material properties m_{1}, m_{2} can be written as
where the subscript t refers to differentiation with respect to time and the constant prefactors c_{12} and c_{21}, which we call the binary abrasion parameters, depend simultaneously on the materials m_{1} and m_{2} of the X and Y particles, respectively.
We also note that in the case of two identical particles (e.g., two particles with identical masses X=Y and identical material properties c_{XY}=c_{YX}) the system (Eqs. 1 and 2) predicts mass evolution according to Sternberg's law. In the case of different masses or properties we still have an infinite lifetime, with one of the particles approaching zero mass asymptotically as time goes to infinity and the other particle approaching a finite mass.
1.2.3 Collective size dynamics
Independently of individual (and binary) abrasion theory there exists broad interest in collective shape and size evolution models tracking mutually colliding populations of N particles (see Fig. 1c). Similar problems arise in particular in the context of coagulation (da Costa, 2015) and dynamic fragmentation processes (Cheng and Redner, 1988). In such collective evolution models the main question is how the size distribution of particles, starting from an initial distribution, evolves in time due to the mutual collisions. These models use a standard framework relying on a socalled collision kernel. In a more general setting, the collision kernel is referred to as the interaction kernel. Our choice of terminology is motivated by the fact that in our case the only interactions are collisions.
The collision kernel can be derived from the binary equations (the physical model of the N=2 case) by incorporating statistical effects, i.e., that collision probability may depend on particle speed or mass. In Domokos and Gibbons (2018) the binary model (Eqs. 1 and 2) was extended to a kernel by introducing an additional scalar parameter r (to which we will also refer as the environmental parameter of the evolution), representing the assumption that on average, only the collision probability depends on particle size and the collision speed is independent of mass:
Note that these equations are identical to the formulas (118) and (119) in Domokos and Gibbons (2018) with α=0 in their notation and taking r=ν, X=V_{X}, Y=V_{Y}, ${c}_{\mathrm{12}}={c}_{\mathrm{21}}=c$. Alternative interpretations of r are also possible; we will discuss the role of the environmental parameter r in detail in Sect. 2.4. Henceforth, in the main body of this paper (apart from Appendix A) we assume that the pebble population is homogeneous, i.e., that the material for all pebbles is identical so we have $c={c}_{\mathrm{12}}={c}_{\mathrm{21}}$ and the sole role of the constant c is to set the timescales. We will incorporate this into the time variable t, and, henceforth, for homogeneous pebble populations, we set c≡1. We will discuss the role and identification of material constants in heterogeneous pebble populations in Appendix A.
Once the kernel has been established, we make the assumption that for large N the collective size evolution is a stochastic process driven by many binary events among the particles, implying that the core of the collective process is still the abovementioned collision kernel. This allows for the construction of the master equation, also known as the Fokker–Planck equation, which describes the time evolution of the particle size distribution. Although the collective abrasion is a stochastic process, in the N→∞ limit the collision kernel will uniquely determine the global evolution of the continuous size distribution. The master equation (or Fokker–Planck equation) expresses this evolution. Determining the master equation based on the collision kernel is the second step in the statistical model.
1.3 Our model
1.3.1 Relationship to earlier models
The aboveoutlined structure is characteristic of coagulation fragmentation models (da Costa, 2015), in particular for nonlinear fragmentation, which describe fragmentation processes triggered by binary collisions of particles. Our model may be regarded as a special case of the nonlinear fragmentation models (Cheng and Redner, 1988) since, in addition to the standard framework adopted in these models, we also make two simplifying assumptions:

we only consider collisions where the relative mass loss is small (i.e., the particles lose only fragments with small relative mass), and

the small fragments generated in the collisions are not considered further in the evolution.
By implementing these two assumptions into the statistical model based on the collision kernel (Eqs. 3 and 4), we take the first step towards establishing the statistical theory of collective size and shape evolution of sedimentary particles. This approach offers multiple methodological advantages. On one hand, by using Eq. (4) as the collision kernel, our statistical model will be compatible with Sternberg's law, so we can expect the collective evolution also to observe this theory, albeit in a statistical sense. On the other hand, we can also expect all our results to be compatible with an extended (future) theory which also describes collective shape evolution based on the unified, volumeweighted geometric theory in Domokos and Gibbons (2018).
1.3.2 Basic notations
To describe our construction we will need to address both the size evolution of individual particles (under the collision kernel) as well as the evolution of size distributions. While particle size appears in both settings, we need to distinguish carefully: in individual and binary models particle size evolves in time; in collective models size distribution evolves in time. As a consequence, in the individual setting the variable denoting size may be differentiated with respect to time; in the collective setting this is not the case. We will use X, Y to denote individual particle sizes (either volume or mass) and we will use x, y to denote the independent variables of size distributions. The time evolution of individual particle size will be denoted by X(t), Y(t) with time derivatives X_{t}(t), Y_{t}(t) (arguments of a function written in subscript will refer to differentiation throughout the paper). The time evolution of size densities will be denoted by f(x,t), f(y,t) with time derivatives f_{t}(x,t), f_{t}(y,t) and size derivatives f_{x}(x,t), ${f}_{y}(y,t).$ We denote the expected value and variance of these size distributions, respectively, by E(t) and W(t) and we will primarily use the relative variance $R\left(t\right)=W\left(t\right)/E(t{)}^{\mathrm{2}}$ to characterize the evolution of the distributions.
1.3.3 Main results
The collision kernel (Eqs. 3 and 4) for mass evolution in Domokos and Gibbons (2018) has one single environmental parameter r, which is inherited by the corresponding Fokker–Planck equation (shown in Sect. 2.2). As we will describe in Sect. 2, the environmental parameter r may, depending on interpretation, represent either the size dependence of the number of collisions or, alternatively, the size dependence of collision energy. Regardless of the interpretation, in Sect. 2.3 we find that the value r=0.5 is critical as it separates two regimes of collective abrasion with a qualitatively different evolution R(t) of the relative variance:

For r>0.5 we find focusing processes with decreasing R(t), approaching R(t)=0 in the limit as time approaches infinity. Here the size distribution converges to a Diracdelta function. This parameter range corresponds to lower energy levels. Natural abrasion processes belonging to this regime will thus amplify the segregating effects of sizeselective transport.

For r<0.5 we find dispersing processes with increasing R(t), thus counteracting sizeselective transport processes. This corresponds to collisional abrasion at higher energy levels.
As collisional abrasion may occur within a broad range of energies, these two basic scenarios of the model (illustrated in Fig. 2) offer an explanation for the broad range of geological observations (Bluck, 1967; Landon, 1930; Carr, 1969; Kuenen, 1964) in relating the relative significance of transport and abrasion in various scenarios. Our model also reflects the universality of Sternberg's law by predicting, regardless of the environmental parameter r, exponential decay as the universal evolution E(t) of the expected value.
In general, the evolution equations generated by Eqs. (3) and (4) for the mean E(t) and the variance W(t) are integrodifferential equations which are hard to solve analytically. To support our claims, we will use three types of approximations:
 a.
We approximate the kernel (Eqs. 3–4) by its truncated Taylor series expansion and investigate the evolution of general initial density functions. This is found in Appendix C1.
 b.
We regard the full kernel; however, we only investigate density functions obtained as a small perturbation of the Dirac delta (i.e., populations of almost identical particles). This is done in Appendices C2 and C3.
 c.
We numerically compute both the discrete and the continuum models. For details see Sect. 3.
We will briefly refer to the first two approximations as the continuum model. In the case of the third approximation we do direct, discrete simulations of finite particle populations; we use the full kernel and we call this the discrete model. One startling feature of the latter (as compared with the former) is the appearance of outliers, i.e., particles substantially larger than the vast majority (illustrated in Fig. 2). As we can observe, the bulk of the density function closely mimics the evolution in the continuum model. The quantitative analogy in the evolution of the relative variation R can also be recovered if we consider a reduced density function ${f}^{*}({t}_{\mathrm{1}},x)$ by omitting the outliers, i.e., by applying an upper cutoff in size, omitting bins containing only one particle. The reduced density function ${f}^{*}({t}_{\mathrm{1}},x)$ is characterized by the reduced relative variation R^{*}, which will decrease in a focusing process; however, in contrast to the continuum model, it will not approach zero but a positive constant.
1.4 Testing of the model for homogeneous pebble populations
As outlined above, our model is defined on two levels: the collision kernel (Eqs. 3–4) we will briefly refer to as the input level as it defines the basic physics of the underlying collisions. The Fokker–Planck equation we will briefly refer to as the output level as it defines the evolution of the mass density function based on the collision kernel. One may test the model at both levels. Below we discuss the case of homogeneous pebble populations where the evolution of the mass distribution is controlled by the single material parameter c and the single environmental parameter r:
 a.
One may test the model at the input level, by fitting the kernel (Eqs. 3–4) to laboratory tests where the abrasion rate is plotted as a function of particle size. Such an experiment could be used to determine the material parameter c for a given homogeneous population. Also, if the laboratory test imitates the environment of the natural process, the environmental parameter r may also be obtained in this manner. We also note that the functional relationship between particle size and abrasion rate will not only depend on the parameters but also on particle size. For details, see Appendix A.
 b.
One may test the model at the output level by measuring the time evolution of full mass distributions and fitting the respective material and environmental parameters c and r to this dataset. While we are not aware of any such public dataset, this could be performed in a laboratory either in a flume or in a drum experiment. In the field the optimal solution appear to be radiotagged pebbles (Bertoni et al., 2016).
The above simple procedures apply only for homogeneous populations. We lay out the procedures for the testing of the model for heterogeneous populations in Appendix A, where we also perform partial testing for the laboratory data obtained by Attal and Lavé (2009).
2.1 General form of the collision kernel
The first simplification described in Sect. 1.3.1 implies that the limit where relative fragment mass approaches zero offers a good approximation; thus it permits a collision kernel of the type used in Ernst and Pagonabarraga (2007), describing continuous mass evolution via coupled ordinary differential equations for the evolution of particles with masses X(t) and Y(t):
where ψ^{1}(X,Y) and ψ^{2}(X,Y) are differentiable (C^{1}) functions, with positive values (i.e., ${\mathbb{R}}^{+}\times {\mathbb{R}}^{+}\to {\mathbb{R}}^{+}$). Symmetry of the binary process implies ${\mathit{\psi}}^{\mathrm{1}}(X,Y)={\mathit{\psi}}^{\mathrm{2}}(Y,X)$, so often superscripts are suppressed and the kernel is simply referred to as $\mathit{\psi}(.,.)$. Selection of the kernel encapsulates not only the physics of binary collisions, it also may include the massdependent probability of collision between two particles. We will discuss the identification of physically sound kernels in Sect. 2.3.
2.2 General form of the master equation
The second simplification in Sect. 1.3.1 permits the construction of the master equation solely based on the collision kernel (by omitting additional terms for the remainder of the fragmented material). These simplifying assumptions also set our model apart from general fragmentation models in another respect: in the latter, constant mass is prescribed as a global time invariant while the (integer) number of particles changes, whereas in our model total mass decreases while the number of particles remains constant and serves as a global invariant.
Using these considerations, for our problem the master equation is found to be
where subscripts stand for partial derivatives. Without loss of generality, the evolution starts at t=0 and we consider the initial distribution of the volume $f(\mathrm{0},x)\equiv {f}_{\mathrm{0}}\left(x\right)$ to be known a priori. Note that contrary to the majority of Fokker–Planck models, our model contains solely the advection term, which readily follows from the deterministic nature of the kernel. Here we aim to determine the collective behavior implied by Eq. (7). Nonetheless, a stochastic kernel would produce diffusion in the master equation. Such a generalization would inevitably reduce the analytic transparency and thus the qualitative predictive capability of the model. Whether or not it is justified from the quantitative point of view can be decided based on extensive testing campaigns.
We aim to understand some scenarios characteristic of pebble populations by investigating the Cauchytype initial value problem associated with Eq. (7), starting at the distribution f_{0} with mean value E_{0}, variance W_{0}, and relative variance ${R}_{\mathrm{0}}:={W}_{\mathrm{0}}/{E}_{\mathrm{0}}^{\mathrm{2}}$.
2.3 Collision kernels
Detailed physical modeling of the collisional event can make the interaction kernel highly complex; for a recent review on kernels see Meyer and Deglon (2011). On the other hand, mathematical studies tend to prefer simple expressions for ψ(X,Y), permitting rigorous, analytical conclusions. Our goal is to find a kernel which has a strong physical basis yet permits an analytical approach; thus it offers a tradeoff between physical and mathematical preferences.
We first consider two simple kernels which satisfy the mathematical requirement of leading to analytically soluble Fokker–Planck equations. However, as we will show, these very analytical results highlight that these kernels are physically not admissible. Next, we investigate the parameterdependent compound kernel suggested in Domokos and Gibbons (2018), which grabs the essential physics of the investigated process, yet the corresponding Fokker–Planck equation still permits analytical conclusions.
First, we consider the summation kernel (denoted by $\mathit{\{}.{\mathit{\}}}^{+}$), where the mass loss rate is proportional to the sum of the masses of the colliding particles:
stating that the rate of mass loss in binary collisions is proportional to the total mass of the two colliding particles. Appendix B1 demonstrates that the relative variance of the mass in the case of the summation kernel follows ${R}^{+}\left(t\right)={R}_{\mathrm{0}}{e}^{\mathrm{2}t}$; hence it is a dispersive process regardless of the initial distribution f_{0}.
In the very same manner let us investigate the product kernel distinguished by the sign $\mathit{\{}.{\mathit{\}}}^{*}$. The product kernel is defined via
According to Appendix B2, the relative variance in this case is constant as ${R}^{*}\left(t\right)={R}_{\mathrm{0}}$ for all t≥0, which means that the model is neither focusing nor dispersing. Note that the time invariance of R^{*}(t) under the product kernel does not imply the invariance of the probability density function (PDF) f(t,x) per se. In addition, we see a polynomial decay in the mass as ${E}^{*}\left(t\right)=(t+{E}_{\mathrm{0}}^{\mathrm{1}}{)}^{\mathrm{1}}$, which contradicts Sternberg's law (Sternberg, 1875), which postulates an exponential decay.
In order to be in accordance with Sternberg's law and to have a control on the evolution of the relative variance, following the lead of Domokos and Gibbons (2018) we investigate the interaction law (Eqs. 3 and 4), which we call a compound kernel, and using the introduced general notation for kernels, we distinguish it with the {.}^{c} sign:
where $\mathrm{0}\le r\le \mathrm{1}$ is a fixed parameter. Henceforth we investigate the evolution of mass density functions under the Fokker–Planck equation derived from Eq. (10). In Appendix C1 we show analytical results for the evolution if the kernel (Eq. 10) is replaced by its truncated Taylor expansion. In Appendices C2 and C3 we show analytical results for the evolution under Eq. (10) using a Dirac delta as initial distribution. The evolution under Eq. (10) with no restrictions for the initial condition is studied numerically. The essential properties of the three investigated kernels are summarized in Fig. 3.
2.4 Interpretation of the parameter r
In natural events, both velocity and collision probability (cross section) may depend on particle size: in laminar flows relative velocity and collision probability is proportional to linear size, while in a turbulent flows velocity could be inversely proportional to linear size and collision probability could be proportional to projected area. In the collision kernel (Eq. 10) both effects (dependence of velocity and dependence of collision probability on speed) are represented by the single scalar parameter r, so one may freely assign various interpretations to this parameter. In Domokos and Gibbons (2018) one particular interpretation was used: the compound kernel was derived using the assumption that particle velocity is independent of the size (e.g., instead determined by the surrounding fluid), but the collision probability works as a power law with particle size, i.e., X^{r}. The effective mass combined with the collision probability gives the kernel in Eq. (10). However, alternative interpretations are possible; the only essential underlying assumption is that we regard a oneparameter family of scenarios. In this family, if velocity is proportional to X^{a} and collision probability is proportional to X^{b}, then we have $r\simeq a+b$.
To have a global view, it may be of interest to estimate the parameter r in two extreme (limiting) scenarios. Laminar flows are characterized by a linear velocity profile. The particles hit each other if their trajectories intersect. The integration of the linear velocity profile combined with a spherical particle shape yields a collision probability proportional to $\sim {X}^{\mathrm{2}/\mathrm{3}}$ or, alternatively, $r\simeq \mathrm{2}/\mathrm{3}$. The other extreme case corresponds to turbulent flows, where we have equipartition; i.e., the kinetic energy of the particles is independent of their size (see, e.g., Uberoi, 1957), implying that particle velocity is proportional to ${X}^{\mathrm{1}/\mathrm{2}}$. Since the area of the cross section is proportional to ${X}^{\mathrm{2}/\mathrm{3}}$ we arrive at a collision probability ${X}^{\mathrm{1}/\mathrm{6}}$ or, alternatively, $r\simeq \mathrm{1}/\mathrm{6}$. As we can see, both extreme scenarios yield r values far away to either side of the critical value ${r}_{\mathrm{crit}}=\mathrm{1}/\mathrm{2}$, so these estimates suggest that smooth steady conditions should result in a focusing and turbulent gaslike behavior in a dispersing process. For a detailed derivation see Appendix D.
In order to examine the validity of these assumptions we made discrete element simulations using the eventdriven method (Lubachevsky, 1991). In eventdriven dynamics, collisions are considered instantaneous and resolved accordingly, which is best suited to obtaining proper collision statistics. We emulated the abovementioned processes by choosing an artificial mass for the particles and simulating a chaotic system. The artificial mass was used to obtain different volumevelocity relations in different scenarios. We found that in chaotic or turbulent systems relative velocities were proportional to $v\sim {X}^{\mathrm{1}/\mathrm{2}}{X}^{\mathrm{1}/\mathrm{3}}$ and the system behaved as the continuum model with $r\sim \mathrm{1}/\mathrm{6}\mathrm{1}/\mathrm{3}$. On the other hand, if velocities were proportional to $v\sim {X}^{\mathrm{1}/\mathrm{6}}{X}^{\mathrm{0}}$, then the system was similar to a continuum model with $r=\mathrm{1}/\mathrm{2}\mathrm{2}/\mathrm{3}$. Thus the discrete element simulations fully support the results of the compound kernel.
2.5 Fluvial abrasion
Here we interpret the intuitive picture of fluvial abrasion in the context of our statistical model. In our model a fluvial environment may be represented by a fluvial population, consisting of N+1 particles: a very large number (N) of small particles X^{i} (i=1, 2, … N) representing the pebbles carried by the river and one very large particle Y representing the riverbed. Such a scenario cannot be explored directly in the context of our continuum model; however, as we will discuss in detail in Sect. 3.3, the discrete model can capture this situation even in the limit of N→∞.
To make a meaningful characterization of geologically relevant scenarios, we will regard two extreme cases which represent brackets on geological processes. In both cases we assume that the mass evolution is driven by binary collisions and we regard the limit as N, Y→∞ (while the masses X^{i} of the small particles remain finite). Since we are interested in the mass evolution of pebbles (and in the current paper we are not interested in the mass evolution of the riverbed), we will denote the relative variance of the pebble population (i.e., all X^{i} particles, the riverbed Y not included) by R(t). Our aim is to establish the sign of R_{t}(t) as the main qualitative feature of collective dynamics, as R_{t}(t)<0 and R_{t}(t)>0 imply focusing and dispersing processes, respectively.
In the first extreme scenario we assume that particles are chosen uniformly from the full fluvial population: i.e., the riverbed has no special role. In this case almost all collisions will happen among a pair of small particles (X^{i}, X^{j}); thus the presence of the riverbed has no impact on the evolution of R(t). For this extreme case all predictions of our continuum model remain valid: r=0.5 will be a critical parameter value above which we see focusing (R_{t}<0) and below which we see dispersing (R_{t}>0) behavior. At the critical value r=0.5 our model predicts neutral behavior with R_{t}=0.
In the second extreme scenario we assume that the small particles exclusively collide with the riverbed (large particle); i.e., we only have (X^{i}, Y)type collisions. This means that the evolution for each of the small particles is an identical, independent twoparticle process governed by the model (Eqs. 1–2) for binary collisional mass evolution. In this process, in the Y→∞ limit each individual small particle X^{i} will thus evolve as
and thus follow Sternberg's law. It is easy to show that for any initial distribution for the masses X^{i}(0), in this process we have R_{t}=0. The large Y particle (riverbed) will lose some mass as well, but in this publication we are not interested in that part of the process.
Intuitively it is clear that any geologically relevant process is in between the above two extreme cases, and, although we do not deliver a rigorous proof, it appears plausible that in a geologically relevant setting R_{t} will also be bounded by the two evolutions predicted for the two extreme scenarios. As for the second extreme scenario we have R_{t}=0; we expect that for any intermediate scenario the sign of R_{t} will agree with the sign of R_{t} based on the first extreme scenario. Our results show that the focusing behavior of the particle size distribution, or lack thereof, depends on interparticle interactions and not on the collisions between the particles and the riverbed. This would imply that all our qualitative predictions remain valid in fluvial environments.
Here we perform computations to illustrate the main results presented in Sect. 1.3.3 by discretizing time with a fixed time step Δt. The discrete model has been simulated with custommade codes in Matlab and Python performing ${M}^{*}[N/\mathrm{2}]$ collisions between pairs during one time step Δt, where M is fixed model parameter and N is the size of the population. The simulation starts with the creation of N particles whose volumes are randomly sampled from the initial distribution f_{0}. Binary collisions are performed on uniformly selected pairs; i.e., all particles have an equal chance of being selected irrespective of their volume. Once a pair is selected, the collision kernel ψ^{c} is applied and the volume decrement is computed with time step $\mathrm{\Delta}t/M$. After the binary collision event both particles with a reduced volume are replaced into the sample. In the presented simulations we set the population size to be N=5000, the time step Δt=0.01, and M=10. In the continuum setting f(t,x) evolves under Eq. (7) with some initial value f_{0}(x). This code uses the operator exponential syntax of a the Chebfun toolbox (Driscoll et al., 2014) in Matlab.
3.1 Focusing and dispersing regimes
The evolution of a pebble population under the compound kernel was simulated both in the frame of discrete and the continuum model, i.e., by direct eventbased simulation and by discretizing the partial differential equation. (see Sect. 1.3.3). The results show excellent agreement with our analytical predictions: $r=\mathrm{1}/\mathrm{2}$ does indeed appear to be a critical parameter in the model. This is illustrated in Fig. 4, where a lognormal distribution is used as an initial value for the evolution.
3.2 Fitted lognormal distribution
Although the lognormal distribution is certainly not invariant under the compound kernel (i.e., an initially lognormal density function does not remain lognormal in the evolution), mass distributions in later time steps highly resemble lognormal distributions. To test this visual observation we fitted lognormal distributions to the computed mass distributions in the discrete simulations. The evolution of the two parameters (respectively, denoted μ and σ) of the lognormal distribution is given in Fig. 5 at values of the parameter r. The criticality of r=0.5 is obvious in this setting, too: while the initially lognormal distribution is almost invariant under the evolution at $r=\mathrm{1}/\mathrm{2}$, the evolution of the parameters μ and σ takes an opposite direction in the parameter space for r=0.0 and r=1.0, respectively. The 95 % confidence levels of the fit confirm the visual intuition: the evolved distributions are close to lognormal: in practical applications an approximation with a lognormal distribution produces an acceptable error.
3.3 Outliers: anomalies in smaller samples
The continuum model describes the N→∞ limit of the system. In the computations shown in Sect. 3.1 and 3.2 we either showed results based on the continuum model or in the direct, discrete simulations we treated large (N=5000) populations. However, if we look at the discrete simulations on smaller samples we may observe unexpected phenomena not recorded in the previous computations. In Fig. 6 we show the mass distribution of a system at r=0.6 with N=2000 particles. The bulk of the histograms can be approximated well with a lognormal distribution. However, there are 12 particles with somewhat larger volume than predicted by the lognormal distribution and one approximately 150 times the median volume (5.3 times the radius). Thus inside the focusing regime we may observe a situation where we have a welldefined narrow distribution which describes the bulk of the particles, but a few might escape from this process and may be left behind, at larger mass. This effect is persistent and it was observed also for the parameter value of r≃0.7.
In order to estimate the robustness of this scenario we use a simple approximation by assuming that all but one particle have volume X and one single, exceptional particle, called the “outlier”, has a volume aX with a≫1. As is demonstrated in Appendix E, the outlier can coexist with the population of the small particles. In the N→∞ limit the condition of such a coexistence reads
The numerical solution of Eq. (12) for equality yields the critical curve a_{crit}(r) on the [r, a] parameter plane, separating systems where outliers may coexist with the population from systems where they may not. While we computed the a_{crit}(r) critical curve for the case of infinitely large populations (the N→∞ limit), we stress the fact that the illustrated phenomenon is inherently discrete and does not arise in the continuum model. We may explain this curious phenomenon in the following manner. Assume that we start from a narrow distribution. Then random fluctuations in the discrete system may create particles with large relative mass (i.e., a large parameter value a). If these fluctuations are sufficiently large to create particles above the critical curve a_{crit}(r), then these outliers will be sustained; otherwise their mass will again approach the average mass of the majority. The critical curve in Fig. 7 shows that in the vicinity of the critical value r=0.5 almost any such fluctuation will be sustained and outliers are likely to survive. However, as the parameter r increases, it becomes increasingly less likely to see sustained outliers. Another observation is that as the likelihood for the existence of outliers decreases, their expected relative size increases, which matches the commonsense observation that the larger the outlier, the less frequently it may be observed. We also note that the relationship between the collection of small particles and the large particle is essentially asymmetrical. While the evolution of the latter is strongly influenced by both the factor a and the control parameter r, the evolution of the density function for the small particles is solely controlled by the latter. In other words, adding one (or a few) very large particles to a collection of many small particles will not alter the fate of the latter, as long as the collisions between a pair of particles are based on a uniform choice.
In this paper we presented the first statistical model for the collective mass evolution of pebble populations under collisional abrasion. While our model is certainly not unique, it is compatible with
 a.
existing geological observations,
 b.
existing geometrical theory of individual and binary abrasion of pebbles,
 c.
existing theory for individual mass evolution of pebbles (Sternberg's law), and
 d.
exiting statistical theory of coagulation and fragmentation.
In the spirit of standard statistical theory for collective evolution, our model is based on two components: (i) the binary collision kernel and, based on that, (ii) the governing equation for the evolution of probability density functions for mass distribution. Regarding (i) we used the model from Domokos and Gibbons (2018), which incorporates the existing theory for individual and binary abrasion; regarding (ii) we used the Fokker–Planck equation, which is broadly used in the theory of coagulation in fragmentation.
Our collision kernel includes the single scalar parameter r which can be associated with the energy level of the collective collisional evolution process. We found that r=0.5 is critical, separating two regimes with fundamentally different behavior: for r>0.5 (lowenergy regime) we found focusing behavior with decreasing relative variance R(t), and for r<0.5 (highenergy regime) we found dispersing behavior with increasing relative variance R(t). In geological terms, this result suggests that in lowenergy environments collisional abrasion acts on mass distributions in unison with sizeselective transport, while in highenergy environments the opposite happens and the two processes counteract each other. In accordance with prevailing geological observations and Sternberg's law, our models predicts exponential decay of particle mass in both energy regimes.
We investigated our model on two levels: (i) as a continuum model by regarding the evolution of the Fokker–Planck equation and (ii) as a discrete model by running discrete eventbased simulations. In the case of the continuum model we derived our results analytically and also from numerical simulation of the Fokker–Planck equation, while in the discrete model we relied on numerical computations. With regard to the existence of the critical parameter r=0.5 and the existence of the focusing and dispersing regimes, the two approaches yielded quantitatively matching results.
Among many small pebbles, large boulders are often visible in mountain ranges or rivers. While this phenomenon is commonly attributed to transport, our model suggests that under some conditions, here again transport and abrasion may act in unison: we identified a curios phenomenon not present in the continuum model but present in the discrete model (even in the N→∞ limit). If the parameter r was in the focusing r>0.5 range but not very far from the critical value r=0.5, the bulk of the distribution narrowed (in accordance with our analytical predictions); however, we could also observe a few particles with substantially larger mass (outliers), escaping the bulk of the distribution. We characterized the mass ratio of outliers versus the mean of the bulk distribution by the parameter a, and we derived a critical curve a_{crit}(r) separating systems where outliers may be observed from those where this may not happen. Our result predicts that larger outliers are less likely to be observable.
While our paper only dealt with size distributions, there exist also related observations on shape: sharp peaks in distributions of axis ratios (also referred to as equilibrium shapes) are mentioned in Bluck (1967), Dobkins and Folk (1970), Landon (1930), Orford (1975), Williams and Caldwell (1988), Ashcroft (1990), Lorang and Komar (1990), Yazawa (1990), and Wald (1990). In Domokos and Gibbons (2012) a plausible argument was presented that equilibrium shapes may emerge on shingle beaches as the result of the interaction of abrasion and transport. We hope that the extension of the statistical theory presented in this paper may be capable of verifying these observations.
In the context of the binary evolution model (Eqs. 1–2) we introduced the binary abrasion parameters c_{12} and c_{21} and for simplicity (since we only aimed to treat homogeneous populations) we used the same notation in the collision kernel (Eqs. 3–4). Here we refine this concept in the statistical setting for heterogeneous populations where we regard the collective evolution of N particles with M≤N different materials m_{i}, (i=1, 2, … M). (The binary case corresponds to N=2; if the two pebbles are made from different material, then we have M=2, and for pebbles with identical materials we have M=1. In the latter case in (Eqs. 1–2) we have ${c}_{\mathrm{12}}={c}_{\mathrm{21}}=c$.)
In the statistical setting the binary abrasion parameters can be organized into an M×M matrix M with entries c_{ij}, (i, j=1, 2, … M). The binary parameter c_{ij} is defined as the constant coefficient in the collision kernel (Eqs. 3–4) associated with the abrasion rate of particles with material m_{i}, bombarded by particles with material m_{j}. Needless to say, the matrix M is not symmetrical; in general we have c_{ij}≠c_{ji}. In particular, if material m_{i} is much harder than material m_{j}, then we expect c_{ij}≪c_{ji}.
Based on the above considerations, the statistical model is controlled by the $M\times M={M}^{\mathrm{2}}$ binary abrasion parameters and the single environmental parameter r. Testing this model can be done along the strategies outlined in Sect. 1.4 for homogeneous populations; however, more detail has to be observed.
 a.
One may test the model at the input level, by fitting the kernel (Eqs. 3–4) to laboratory tests for pairwise selected materials m_{i}, m_{j}. In such a test the abrasion rate of particles of material m_{i} under abrasion from particles of material m_{j} is plotted as a function of particle size of the abraded particle (with material m_{i}). Such experiments can be used to determine the binary abrasion parameters c_{ij} for a given heterogeneous population. If the laboratory test imitates the environment of the natural process, the environmental parameter r may also be obtained in this manner. We will show such an example below.
 b.
One may test the model at the output level by measuring the time evolution of full mass distributions and fitting the M^{2} material parameters c_{ij} and the environmental parameter r to this dataset.
Next we show an example for testing the model at the input level by using the data obtained in Attal and Lavé (2009). Here the authors report on flume experiments where they measured the abrasion rate E_{d} of individual limestone gravels with a diameter between D=9 and D=39 mm mixed in approximately 400 g of 10–18 and 18–28 mm granitic gravel. In our terminology, we have M=2 (two materials) and we will use m_{1} for the limestone and m_{2} for the granite. The joint evolution of such a heterogeneous population is described by M^{2}=4 binary material constants: c_{11}, c_{12}, c_{21}, and c_{22}. Attal and Lavé (2009) were primarily interested in the abrasion rates for limestone and they produced the E_{d}(D) plots for these particles. In this experiment we may assume that the abrasion of the limestone pebbles was exclusively due to collisions with the granitic gravel (i.e., we disregard limestone–limestone collisions). Thus the only relevant collisions are between limestone and granite, and for the mass loss X(t) of the limestone we will use Eq. (3) with X and Y denoting the volumes of the colliding limestone and granitic particles, respectively, and c_{12} denoting the binary abrasion parameter associated with limestone being abraded by granite (not reported, but we may assume c_{21}≪c_{12}). If we replace the volume of the granitic particles by their average, the abrasion rate as function of its diameter can be calculated numerically. Note that the abrasion rate E_{d} in our notation reads
We fitted Eq. (3) to the dataset provided in Attal and Lavé (2009). We minimized the mean square error (with respect to the results in Attal and Lavé, 2009) for the parameters r and c_{12} and obtained r=0.19 and c_{12}=0.28. Our fitted curves are illustrated in Fig. A1 showing fair agreement between the data and the fitted model. The value of the environmental parameter is in the range where we expect dispersing behavior, as we discussed in Appendix D, which is in accordance with the target of the original experiment which simulated abrasion in fluvial environments. We note that the same parameter pair r=0.19 and c_{12}=0.28 is valid for both limestone experiments (i.e., these parameters do not depend on the size of the particle). Our fit appears to be consistent in this respect.
B1 Summation kernel
Differential equations governing the time evolution of the first and second moments can be readily obtained; hence the mean E^{+}(t) and variance W^{+}(t) follow the following initial value problems (IVPs):
It follows that both the expectation and the variance exhibit exponential decay, namely ${E}^{+}\left(t\right)={E}_{\mathrm{0}}{e}^{\mathrm{2}t}$ and ${W}^{+}\left(t\right)={W}_{\mathrm{0}}{e}^{\mathrm{2}t}$. It is straightforward to show that the relative variance R^{+}(t) increases exponentially:
B2 Product kernel
In the case of the product kernel the IVPs describing the evolution of the mean E^{*}(t) and variance W^{*}(t), respectively, read
Here the decay of the mean and the variance are polynomial as we find
which result in a steady relative variance R^{*}(t), determined by the initial distribution f_{0}. Specifically
C1 Truncated compound kernel
The truncated compound kernel is obtained from the compound kernel as the truncated Taylor polynomial computed at y=x with an (O(y−x)^{2}:
Using the master equation, the following Cauchy problems are found that define the evolution of the mean and the variance:
Solution of these ordinary differential equations (ODEs) yields the evolution of the relative variance as
C2 A population of identical particles preserved
Here we show that a population of identical particles, characterized by a Diracdelta function as input PDF is preserved in the model with the compound kernel regardless of the value of parameter r. Without loss of generality, we investigate the evolution from the f_{0}=δ(1) initial condition, where δ(x) denotes the Diracdelta function at x. Obviously, E_{0}=1 and W_{0}=0. We show that now $f(t,x)=\mathit{\delta}\left(c\right(t\left)\right)$ holds for any t>0. Let us assume that at some ${t}^{*}\ge \mathrm{0}$, the distribution is $f({t}^{*}*)=\mathrm{0}$. Observe that
The time derivative of the mean can be computed via
where we used Eq. (7), applied integration by parts, and employed Eq. (C5). Similarly, the evolution of the variance is found to follow
This shows that the variance of the distribution is constant, and it started at W_{0}=0; i.e., it vanishes entirely in its the evolution. In other words, we have a Diracdelta (degenerate) distribution at any t≥0. Employing Eq. (C6) we find that the location c(t) follows the initial value problem ${c}_{t}\left(t\right)=\frac{\mathrm{1}}{\mathrm{2}}c\left(t\right)$ with c(0)=1; hence $c\left(t\right)=\mathrm{exp}(\frac{t}{\mathrm{2}})$.
C3 Dispersing and focusing behavior identified in the population of almost identical particles
As the model lacks diffusion, the behavior of a degenerate distribution with all the mass concentrated at a single value is worth studying because longterm existence of a set consisting of identical particles can take place in the model. In Appendix C2 we show that a population of identical particles remains identical in our model. In other words, the time invariance of the Diracdelta distribution holds in our model, regardless of the value of the parameter r. Nevertheless, the value of r affects the stability of that Dirac delta: next we show that the evolution for a population of almost identical particles (i.e., a perturbed version of the Diracdelta distribution) is either focusing or dispersing, depending on the value of r. To see this, we define a perturbed distribution. Let ε>0 be a fixed parameter and define
It is straightforward to show that $\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\widehat{f}}_{\mathrm{0}}\left(y\right){\mathit{\psi}}^{\mathrm{c}}(x,y)\mathrm{d}y={\mathit{\psi}}^{\mathrm{c}}(x,\mathrm{1})$, ${E}_{\mathrm{0}}^{\mathrm{c}}=\mathrm{1}$, and ${M}_{\mathrm{2}}^{\mathrm{c}}\left(t\right):=\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}f(t,x){x}^{\mathrm{2}}\mathrm{d}x$ with ${M}_{\mathrm{2}}^{\mathrm{c}}\left(\mathrm{0}\right)=\mathrm{1}+\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\epsilon}}^{\mathrm{3}}$. We aim to investigate the sign of ${R}_{t}^{\mathrm{c}}$ at t=0. Since ${R}^{\mathrm{c}}\left(t\right)={M}_{\mathrm{2}}^{\mathrm{c}}\left(t\right){E}^{\mathrm{c}}(t{)}^{\mathrm{2}}\mathrm{1}$, we need to study the sign of ${M}_{\mathrm{2},t}^{\mathrm{c}}\left(\mathrm{0}\right){E}^{\mathrm{c}}\left(\mathrm{0}\right)\mathrm{2}{E}_{t}^{c}\left(\mathrm{0}\right){M}_{\mathrm{2}}^{\mathrm{c}}\left(\mathrm{0}\right)$. Integration by parts yields
where algebraic manipulations leads the last equality. In accordance with the results on the truncated model, we found that $r=\mathrm{1}/\mathrm{2}$ is critical. At $r<\mathrm{1}/\mathrm{2}$ the relative variance ${R}_{t}^{\mathrm{c}}$ is positive, it increases for any ε>0; i.e., the population of identical particles is unstable and small perturbations disperse the mass distribution. At $r>\mathrm{1}/\mathrm{2}$ the relative variance ${R}_{t}^{\mathrm{c}}<\mathrm{0}$, which shows that the population of identical particles is stable; the model is focusing.
In the paper we assumed that the particle collision probability depends on the volume of the particles as
Here we investigate two extreme scenarios, associated with the collision probabilities P_{smooth}(X) and P_{turbulent}(X) where we expect r to assume its extreme values.
The first is the smooth gradient flow. In such a case the driving fluid has a strong, but on a particle size scale constant, velocity gradient in one of the spatial directions. Such situations may arise, e.g., in shallow water layers. Here the relative velocity of the particles grows with the distance. If we are at distance u from the center of the particle in the direction of the flow velocity gradient, the collision probability P_{smooth}(X) can be estimated by the product of the velocity difference and the linear cross section of the particles (note that $R\equiv {X}^{\mathrm{1}/\mathrm{3}}$ is the linear size of the particle):
Based on Eq. (D1), this gives us an estimate for high $r=\mathrm{2}/\mathrm{3}$.
The other extreme case is a fully chaotic motion where equipartition takes place (Uberoi, 1957). Thus the kinetic energy of the particles ($\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho}X{v}^{\mathrm{2}}$) is independent of their volume. Thus the speed of the particles must be proportional to ${X}^{\mathrm{1}/\mathrm{2}}$. If we disregard correlations, the particles have a cross section proportional to their projected area which is proportional to ${X}^{\mathrm{2}/\mathrm{3}}$. Combining the two gives us
and based on Eq. (D1) we obtain $r=\mathrm{1}/\mathrm{6}$. Thus it is possible to have physical scenarios apparent in nature where the value of r falls to either side of the critical value of r_{c}=0.5 with a large enough margin.
Let us have a sample with N particles with (N−1) having the identical volume X. The last particle is an outlier with volume aX, where a≪1. In a single binary collision, a hit between particles with volume X is called an Atype event, while a collision with the outlier being involved is a Btype event. Based on discrete probabilistic considerations, the probability of an Atype event equals $(N\mathrm{2})/N$ and a Btype event is $\mathrm{2}/N$. In the Atype event the average size $\stackrel{\mathrm{\u203e}}{X}$ of the particles with volume X after the collision that lasts for Δt reads
Computing $aX/\stackrel{\mathrm{\u203e}}{X}$ and truncating the Taylor series expansion in Δt after linear terms around the value Δt=0 yields the time derivative of the parameter a associated with an Atype event:
In the case of the Btype event both the outlier and one of the small particles follow the compound kernel via
The second equation is employed to compute the average volume of the small particles (i.e., $\stackrel{\mathrm{\u203e}}{X}$ associated with this event). Now we need to truncate the Taylor series of $\frac{aX(aX{)}_{t}\mathrm{\Delta}t}{X\stackrel{\mathrm{\u203e}}{X}}$ at Δt=0. After algebraic manipulations we find
Considering the probabilities of events A and B we arrive at
Note that an increase in the value of a, i.e., a_{t}>0 implies that the outlier is getting further from the population. In the case of the N→∞ limit we find
Here the sign of the expression in the brackets determines the sign of a_{t}, which coincides with Eq. (12) in the text. One can also show that if there exist a_{crit}>1 such that a_{t}=0 at a_{crit}, then a_{t}>0 for any a>a_{crit}. Hence, we need the a_{crit}>1 that makes the expression in brackets vanish. Existence of such a critical value can be shown for the case with finitely many particles, too. As we can see, sufficiently large outliers may coexist with the population in the long run. The control parameter r determines how large an outlier needs to be for sustained coexistence.
The archived code, referred to in this paper, is available at the following public repository: https://doi.org/10.5281/zenodo.4634368 (Sipos, 2021).
GD proposed the problem and supervised the research. AAS carried out the analytical and numerical study of the continuous model. TJ developed the discrete numerical model. GD, AAS, and JT wrote the paper.
The authors declare that they have no conflict of interest.
The authors are indebted to David J. Furbish, Sebastien Carretier, and Duccio Bertoni for reviews and feedback. Responding to their comments helped to improve the paper substantially. The authors sincerely thank Jérôme Lavé and Mikaël Attal for sharing the original dataset of their abrasion experiments. The research reported in this paper was supported by the BME Water Sciences & Disaster Prevention TKP2020 Institution Excellence Subprogram, grant no. TKP2020 BMEIKAVIZ and the NKFIH grant K134199.
This research has been supported by the Nemzeti Kutatási Fejlesztési és Innovációs Hivatal (grant no. K134199) and the Institution Excellence Subprogram (grant no. TKP2020 BMEIKAVIZ).
This paper was edited by Simon Mudd and reviewed by Duccio Bertoni and Sebastien Carretier.
Ashcroft, W.: Beach pebbles explained, Nature, 346, 227, 1990. a
Attal, M. and Lavé, J.: Changes of bedload characteristics along the Marsyandi River (central Nepal): Implications for understanding hillslope sediment supply, sediment load evolution along fluvial networks, and denudation in active orogenic belts, in: Tectonics, Climate, and Landscape Evolution, Geological Society of America, Boulder, CO, 2006. a
Attal, M. and Lavé, J.: Pebble abrasion during fluvial transport: Experimental results and implications for the evolution of the sediment load along rivers, J. Geophys. Res.Earth, 114, F04023, https://doi.org/10.1029/2009JF001328, 2009. a, b, c, d, e, f, g
Bertoni, D., Sarti, G., Grottoli, E., Ciavola, P., Pozzebon, A., Domokos, G., and NovákSzabó, T.: Impressive abrasion rates of marked pebbles on a coarseclastic beach within a 13month timespan, Mar. Geol., 381, 175–180, 2016. a, b
Bird, E.: Lateral Grading of Beach Sediments: A Commentary, J. Coast. Res., 12, 774–785, 1996. a
Bloore, F. J.: The Shape of Pebbles, Math. Geol., 9, 113–122, 1977. a, b
Bluck, B. J.: Sedimentation of beach gravels; examples from South Wales, J. Sediment. Res., 37, 128–156, 1967. a, b, c
Brewer, P. A. and Lewin, J.: InTransport Modification of Alluvial Sediment: Field Evidence and Laboratory Experiments, John Wiley & Sons, Ltd, Hoboken, NJ, 23–35, 1993. a
Carr, A. P.: Size grading along a pebble beach; chesil beach, England, J. Sediment. Res., 39, 297–311, 1969. a, b, c, d, e
Cheng, Z. and Redner, S.: Scaling Theory of Fragmentation, Phys. Rev. Lett., 60, 2450–2453, 1988. a, b
da Costa, F. P.: Mathematical Aspects of CoagulationFragmentation Equations, in: Mathematics of Energy and Climate Change, edited by: Bourguignon, J.P., Jeltsch, R., Pinto, A. A., and Viana, M., Springer International Publishing, Cham, 83–162, 2015. a, b
Dingle, E. H., Attal, M., and Sinclair, H. D.: Abrasionset limits on Himalayan gravel flux, Nature, 544, 471–474, https://doi.org/10.1038/nature22039, 2017. a
Dobkins, J. E. and Folk, R. L.: Shape development on TahitiNui, J. Sediment. Res., 40, 1167–1203, 1970. a
Domokos, G. and Gibbons, G. W.: The evolution of pebble size and shape in space and time, P. Roy. Soc. A, 468, 3059–3079, 2012. a
Domokos, G. and Gibbons, G. W.: The Geometry of Abrasion, in: New Trends in Intuitive Geometry, edited by: Ambrus, G., Bárány, I., Böröczky, K. J., Fejes Tóth, G., and Pach, J., Springer, Berlin, Heidelberg, 125–153, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m
Driscoll, T. A., Hale, N., and Trefethen, L. N.: Chebfun Guide, Pafnuty Publications, Oxford, 2014. a
Ernst, M. H. and Pagonabarraga, I.: The Nonlinear Fragmentation Equation, J. Phys. A, 40, F331–F337, 2007. a
Fedele, J. J. and Paola, C.: Similarity solutions for fluvial sediment fining by selective deposition, J. Geophys. Res.Earth, 112, F02038, https://doi.org/10.1029/2005JF000409, 2007. a
Ferguson, R., Hoey, T., Wathen, S., and Werritty, A.: Field evidence for rapid downstream fining of river gravels through selective transport, Geology, 24, 179–182, 1996. a
Firey, W. J.: Shapes of worn stones, Mathematika, 21, 1–11, 1974. a, b
Gleason, R., Blackley, M. W. L., and Carr, A. P.: Beach stability and particle size distribution, Start Bay, J. Geol. Soc., 131, 83–101, 1975. a
Hansom, J. D. and Moore, M. P.: Size Grading along a Shingle Beach in Wicklow, Ireland, J. Earth Sci., 4, 7–15, 1981. a
Huber, M. L., Lupker, M., Gallen, S. F., Christl, M., and Gajurel, A. P.: Timing of exotic, fartraveled boulder emplacement and paleooutburst flooding in the central Himalayas, Earth Surf. Dynam., 8, 769–787, https://doi.org/10.5194/esurf87692020, 2020. a
Kuenen, P. H.: Experimental Abraison: 6. Surf Action, Sedimentology, 3, 29–43, 1964. a, b
Kuenen, P. H. and Migliorini, C. I.: Turbidity Currents as a Cause of Graded Bedding, J. Geol., 58, 91–127, 1950. a
Landon, R. E.: An Analysis of Beach Pebble Abrasion and Transportation, J. Geol., 38, 437–446, 1930. a, b, c
Lewis, W. V.: The Effect of Wave Incidence on the Configuration of a Shingle Beach, Geogr. J., 78, 129–143, 1931. a
Lorang, M. and Komar, P.: Pebble shape, Nature, 347, 433–434, 1990. a
Lubachevsky, B. D.: How to simulate billiards and similar systems, J. Comput. Phys., 94, 255–283, 1991. a
Meyer, C. J. and Deglon, D. A.: Particle collision modeling – A review, Miner. Eng., 24, 719–730, 2011. a
Miller, K. L., Szabó, T., Jerolmack, D. J., and Domokos, G.: Quantifying the significance of abrasion and selective transport for downstream fluvial grain size evolution, J. Geophys. Res.Earth, 119, 2412–2429, 2014. a
Neate, D. J. M.: Underwater pebble grading of Chesil Bank, Proc. Geol. Assoc., 78, 419–426, 1967. a
NovákSzabó, T., Sipos, A. Á., Shaw, S., Bertoni, D., Pozzebon, A., Grottoli, E., Sarti, G., Ciavola, P., Domokos, G., and Jerolmack, D. J.: Universal characteristics of particle shape evolution by bedload chipping, Sci. Adv., 4, aao4946, https://doi.org/10.1126/sciadv.aao4946, 2018. a
Orford, J. D.: Discrimination of particle zonation on a pebble beach, Sedimentology, 22, 441–463, 1975. a
Paola, C., Parker, G., Seal, R., Sinha, S. K., Southard, J. B., and Wilcock, P. R.: Downstream fining by selective deposition in a laboratory flume, Science, 258, 1757–1760, 1992. a
Sipos, A. A.: aasipos/ESD_collective_paper_code: First release (Version v.1.0.2), Zenodo, https://doi.org/10.5281/zenodo.4634368, 2021. a
Sternberg, H.: Untersuchungen uber Langen und Querprofilgeschiebefuhrender Flusse, Z. Bauwes., 25, 486–506, 1875. a, b, c
Szabó, T., Fityus, S., and Domokos, G.: Abrasion model of downstream changes in grain shape and size along the Williams River, Australia, J. Geophys. Res.Earth, 118, 2059–2071, https://doi.org/10.1002/jgrf.20142, 2013. a
Szabó, T., Domokos, G., Grotzinger, J. P., and Jerolmack, D. J.: Reconstructing the transport history of pebbles on Mars, Nat. Commun., 6, 8366, https://doi.org/10.1038/ncomms9366, 2015. a
Uberoi, M. S.: Equipartition of energy and local isotropy in turbulent flows, J. Appl. Phys., 28, 1165–1170, 1957. a, b
Wald, Q. R.: The form of pebbles, Nature, 345, 211, 1990. a
Whittaker, A. C., Duller, R. A., Springett, J., Smithells, R. A., Whitchurch, A. L., and Allen, P. A.: Decoding downstream trends in stratigraphic grain size as a function of tectonic subsidence and sediment supply, GSA Bull., 123, 1363–1382, 2011. a
Williams, A. T. and Caldwell, N. E.: Particle size and shape in pebblebeach sedimentation, Mar. Geol., 82, 199–215, 1988. a
Yazawa, T.: More pebbles, Nature, 348, 398, 1990. a
 Abstract
 Introduction
 Modeling collective size dynamics
 Numerical results
 Conclusions
 Appendix A: Testing the model for heterogeneous pebble populations
 Appendix B: Some properties of the kernels in Sect. 2.3
 Appendix C: Approximate investigation of the compound kernel
 Appendix D: Estimating physically possible values of r
 Appendix E: Investigation of outliers in finite samples
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Modeling collective size dynamics
 Numerical results
 Conclusions
 Appendix A: Testing the model for heterogeneous pebble populations
 Appendix B: Some properties of the kernels in Sect. 2.3
 Appendix C: Approximate investigation of the compound kernel
 Appendix D: Estimating physically possible values of r
 Appendix E: Investigation of outliers in finite samples
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References