## Abstract

We propose a mathematical model which suggests that the two main geological observations about shingle beaches, i.e. the emergence of predominant pebble size ratios and strong segregation by size, are interrelated. Our model is based on a system of ordinary differential equations (ODEs) called the *box equations* that describe the evolution of pebble ratios. We derive these ODEs as a heuristic approximation of Bloore's partial differential equation (PDE) describing collisional abrasion and verify them by simple experiments and by direct simulation of the PDE. Although representing a radical simplification of the latter, our system admits the inclusion of additional terms related to frictional abrasion. We show that non-trivial attractors (corresponding to predominant pebble size ratios) only exist in the presence of friction. By interpreting our equations as a Markov process, we illustrate by direct simulation that these attractors may only be stabilized by the ongoing segregation process.

## 1. Introduction

The shape of pebbles has been a matter of discussion since at least the time of Aristotle (Krynine 1960). In general, the central question is whether particular pebble shapes emerge from the abrasion and transport processes. Aristotle himself claimed that spherical shapes dominate. However, as Aristotle also observed, abrasion is a complex interaction between the abraded particle and for the abrading environment, represented by ‘other objects’ (i.e. other pebbles) where not only local properties (e.g. curvatures) but also semi-global effects caused by particle shape as well as global effects caused by particle transport play an important role. In this process pebbles mutually abrade each other, defining the time evolution both for the abraded pebble and for the abrading environment, represented by other particles subject to particle transport. In the simplest approach, one neglects the latter effects and regards the abrasion of a single pebble in a constant environment. Aristotle's model for individual abrasion may be translated as
1.1where *v* is the speed with which the pebble's surface moves along the inward normal, *R* is the radial distance from the centre of gravity of the pebble being abraded, *f* is the monotonically increasing function of *R* only, and, in particular, *f* is independent of time. Note that (1.1) is a partial integro-differential equation, since the location of the centre of gravity is determined by time-dependent integrals. The modern partial differential equation (PDE)-based theory of individual abrasion appears to have started with the pioneering work, both experimental and theoretical, of Lord Rayleigh^{1} (Rayleigh 1942, 1944*a*,*b*). For some earlier works see Palmer (1834) and Dobkins & Folk (1970). Rayleigh mainly considered axisymmetric pebbles and observed that the ultimate shape was not necessarily, indeed seldom, spherical. He found that some pebbles' shapes are far from ellipsoidal, being much more discoid in shape. He asserted that abrasion cannot be a simple function of the Gauss (or *specific*) curvature *K*.

On the mathematical side, Firey (1974) initiated a rigorous study by adopting a PDE model rejected by Rayleigh in which the shape evolved according to what Andrews (1999) calls the *flow by the Gauss curvature*; that is, he studied the PDE
1.2where *c* is a constant, and proved that all convex shapes ultimately converge to the sphere under the action of (1.2). Note that the word *flow* is being used in the sense in dynamical systems theory and should not be confused with physical fluid flow. Later, this proof was substantially amplified by Andrews (1999). Recently, Durian *et al.* (2006) investigated the statistical distribution of Gaussian curvature on pebble shapes.

The physical assumption underlying Firey's model is that the abraded particle (pebble) undergoes a series of small collisions with a very large, smooth abrader, and this might be the case when pebbles are carried by a fast river and collide repeatedly with the riverbed, a process called bedload transport. This concept of *collisional flows* was substantially generalized by Bloore, who studied abrasion by spherical abraders of radius *r* (Firey's model corresponds to ). With brilliant intuition, Bloore (1977) arrived at the PDE
1.3where *b* and *c* represent, as Bloore describes, the average or ‘effective’ values for *r* and *r*^{2}, respectively, and *H* is the mean curvature. Bloore also showed that a spherical shape of radius *R* abraded by an abrader of radius *r* is stable under the action of (1.3) if
1.4Apparently, (1.3) consists of three terms: the constant first term that corresponds to the so-called Eikonal equation, the mean curvature term 2*bH* and the Gaussian term *cK*. The latter was studied by Firey, while for the mean curvature flow Huisken (1984) showed that it also converges to the sphere. In the case of the Eikonal, we know that the *outward* flow (which can also be imagined as time-reversed abrasion) converges to the sphere, so, simply put, the question is whether the three terms can ‘balance’, i.e. homothetic solutions can exist.

The rigorous foundation of the collisional model requires a probabilistic approach, as was already conjectured by Hilbert & Cohn-Vossen (1932). Based on the results from stochastic and integral geometry (Schneider & Weil 2008), one can not only justify Bloore's equations but also identify the physical meaning of the coefficients *b* and *c* in (1.3) (Várkonyi & Domokos 2011)
1.5where *M* is the integrated mean curvature and *A* is the area of the *abrading* particles, assumed convex and identical. In the case of spherical abraders, this is identical with Bloore's interpretation of the coefficients. By *Minkowski's inequality* for convex bodies (Minkowski 1897), the quantities *b* and *c* may not be given arbitrarily; rather, based on (1.5), they must satisfy
1.6The main question now is whether collisional flows defined in (1.3) and (1.5) can explain the observation of spherical and non-spherical pebble shapes, i.e. whether (1.3) and (1.5) admit homothetic solutions and what they look like.

Independently of the mathematical work that we have briefly reviewed, there has been considerable geological interest in pebble shape, which has partly focused on individual pebble shapes and partly on global understanding of pebble transport and spatial distribution on large shingle beaches. An extensive list of references about individual pebble shapes up to 1970 is contained in Dobkins & Folk (1970).

One of the most remarkable accounts of pebble shape and distribution is provided by Carr (1969) based on the measurement of approximately 100 000 pebbles on Chesil Beach, Dorset, UK. In summarizing his results, Carr first notes that *a pebble appears to retain its proportions irrespective of size and angularity*. Apparently, the same phenomenon was described by Landon (1930), who claimed that coastal pebbles evolve towards what he calls *equilibrium shapes*. The subject was taken up in the 1990s in *Nature* (Ashcroft 1990; Lorang & Komar 1990; Yazawa 1990) following a note by Wald (1990) claiming that pebbles with three planes of symmetry typically had the form of flattened ovoids in the proportion 7 : 6 : 3. In the absence of pebble transport, these observations would imply the existence of stable (attracting) homothetic (self-similar) solutions of the governing equations. As we will point out, in the absence of transport such solutions do not exist.

However, Carr also provides mean values and sample variations for maximal pebble size along lines orthogonal to the beach. These plots reveal pronounced segregation by maximal size, i.e. on shingle beaches pebbles of roughly similar maximal sizes appear to be spatially close to each other. This phenomenon, also confirmed by Bird (1996), Gleason *et al.* (1975), Hansom & Moore (1981), Kuenen & Migliorini (1950), Landon (1930) and Neate (1967), is closely related to the global transport of pebbles by waves (Lewis 1931; Carr 1972). Grading not only by size but also by shape has also been reported; the most notable description is due to Bluck (1967), reporting on what he calls *equilibrium states* of pebble collections corresponding to stationary shape distributions (see also Landon 1930; Orford 1975; Williams & Caldwell 1988).

Apparently, the distributions of pebble shapes and size are controlled by two fundamental geological processes: particle abrasion and particle transport. Which of these processes dominates may depend on the geological location and also on time; however, geologists appear to agree that the role of transport may not be neglected. Another general observation is that this process approaches equilibria in the form of stationary pebble-shape distributions with dominant peaks. A detailed account of the interaction of abrasion and transport is given by Landon (1930) investigating the beaches on the west shore of Lake Michigan. He attributes shape variation to a mixture of abrasion and transport. Kuenen (1964) discusses Landon's observations, but disagrees with the conclusions and attributes the shape variation primarily to transport. Carr (1969) observes dominant shape ratios emerging as a result of abrasion and size grading, whereas Bluck (1967) describes beaches in South Wales where equilibrium distributions of shape are reached primarily by transport and abrasion plays a minor role.

These observations call for a mathematical model in which both abrasion and transport are included and one may study their complex interaction in the entire range from pure abrasion to pure segregation. Our goal in this study is to make one step towards such a general model which is based on the classical theory of particle abrasion but also admits the inclusion of global transport. We derive a system of representative ordinary differential equations (ODEs) called the *box equations* that provide a huge simplification of Bloore's PDEs but that, nevertheless, allow us to study the spatio-temporal process. In this framework we can show that in accord with the observations of Landon, Carr and Bluck, in the absence of global transport and segregation, stationary distributions with sharp peaks (centred around homothetic solutions of the deterministic equations) *may not emerge* in a stable manner. After verifying our analytical results both by direct simulation of Bloore's PDE and by simple experiments we also show that abrasion and size grading alone (most clearly stated by Carr) might be sufficient for producing persistent dominant basic proportions. Strong shape grading (as described by Bluck) enhances this process even further. Our model creates a framework in which both extreme cases (pure abrasion in the absence of transport and pure transport in the absence of abrasion), as well as their combination, may be studied.

## 2. The collisional box equations

Motivated partly by the above-described geological observations, instead of the numerical study of (1.3) we propose representing shapes by orthogonal bounding boxes (defined by six orthogonal bounding planes) with sizes 2*u*_{1}≤2*u*_{2}≤2*u*_{3}, so that the attrition speed *v* is replaced by ; the local quantities *K* and *H* we obtain from the ellipsoid that fits inside the box
2.1Since we are interested in the existence of self-similar solutions, we introduce the box ratios *y*_{1}=*u*_{1}/*u*_{3} and *y*_{2}=*u*_{2}/*u*_{3}, we denote the origin (*y*_{1},*y*_{2})=(0,0) by *O* and the sphere at (*y*_{1},*y*_{2})=(1,1) by *S* and to simplify the notation we write *y*_{3}=*u*_{3}. This approach has the advantage that most of the available geological data are expressed in terms of these variables; in fact, the representation in the (*y*_{1},*y*_{2}) plane is often referred to as the *Zingg triangle* (Zingg 1935). Now the PDE (1.3) is reduced to the system of ODEs in the (*y*_{1},*y*_{2},*y*_{3}) system
2.2and
2.3where
2.4As can be seen from (2.2) and (2.3), the full flow is three-dimensional; however, the component flows (Eikonal, mean curvature and Gaussian) have autonomous planar components (illustrated in figure 1):
2.5where ()′=d/d*y*_{1}. By introducing the vector notation **y**=[*y*_{1},*y*_{2},*y*_{3}]^{T}, ** F**=[

*F*

_{1},

*F*

_{2},

*F*

_{3}]

^{T}, (2.2) and (2.3) can be rewritten as 2.6Our main concern is the detection of fixed points and self-similar (homothetic) solutions. According to (2.3), , so (2.2) and (2.3) do not have fixed points. Homothetic solutions satisfying 2.7may exist; in fact, spheres at

*S*are a homothetic solution of (2.2) and (2.3). Apart from spheres, (2.2) and (2.3) will only admit homothetic solutions if (2.2) is autonomous and has fixed points; we explore this case in §2

*b*. If (2.2) is not autonomous, then homothetic solutions do not exist; nevertheless, one may still look at solutions of 2.8as ‘temporary’ homothetic solutions; we explore this in §2

*a*. Whether permanent or temporary, the existence of homothetic solutions requires that the (

*y*

_{1},

*y*

_{2}) projection of the flow should be balanced and this observation admits the formulation of a

*necessary*condition. By considering that the flow is a linear combination of three components (2.5), two of which converge to the sphere, we can observe that if either (2.7) or (2.8) is to be satisfied, then the slope of the Eikonal should be ‘bracketed’ by two other slopes 2.9Based on (2.5), we can readily verify that in the interior

*y*

_{1}<

*y*

_{2}<1 of the Zingg triangle this is indeed the case since 2.10i.e. the necessary condition (2.9) is always fulfilled (see figure 4

*a*).

Box flows represent only a heuristic approximation of the PDE and we may use several alternative definitions of the bounding box when measuring the box dimensions. A selection of special approximations based on extremal properties has been studied and a famous example is the unique John ellipsoid with maximal volume (John 1948). Box equations are the basic tool enabling us to follow Aristotle's ideas and to include semi-global and global effects in our model. In contrast, however, box equations are Aristotelian also in the local sense, since in the box approximations of the Bloore flows (1.3) attrition speeds are a monotonically growing function of the distance *u*_{i}. When deriving the box ODEs from the PDE (1.3), the abrading environment (represented by the coefficients *b* and *c*) can also be approximated by its box dimensions *v*_{1}≤*v*_{2}≤*v*_{3} and box ratios *z*_{1}=*v*_{1}/*v*_{3} and *z*_{2}=*v*_{2}/*v*_{3} and maximal box size *z*_{3}=*v*_{3}. According to (1.5), the coefficients *b* and *c* depend on the surface integrals *M* and *A*, and the latter we derive directly from the abrading orthogonal box, with ratios *z*_{i},
2.11The same quantities can be expressed for the unit cube as *M*_{1}=6*π*, *A*_{1}=24, so we have
2.12We remark that (1.6) remains valid also in the box approximations. Whereas the PDE model (1.3) is describing the collisional abrasion of an *individual* pebble in the abrading environment represented by the coefficients *b* and *c*, the box equations suggest a natural generalization which tracks the *mutual* abrasion of two pebbles **y** and ** z**, each of which represents the abrading environment for the other. (Alternatively, one may think of two pebble

*types*abrading each other. We will return to the proper statistical interpretation of these equations in §4.) Based on (2.6) and (2.12), we obtain the simultaneous ODEs 2.13and 2.14describing the (approximate) evolution of box ratios under a purely collisional abrasion process. At first sight, the mutual abrasion of two pebbles may not seem to be a realistic model for the shape evolution of vast pebble collections on shingle beaches; nevertheless, we will show that (2.13) and (2.14) open up the possibility of describing the global interaction between the abraded pebble and the abrading environment, including segregation. Another important feature of the box equations is that they can be easily extended to model abrasion beyond collisional processes. First we look at two special collisional processes.

### (a) The stationary case

First we study the case where pebble ** z** is infinitely harder than pebble

**y**, i.e. we have

**(**

*z**t*)≡

**(0), so (2.13) and (2.14) are reduced to the three-dimensional system 2.15which we call the**

*z**steady-state flow*. Since here we have constant coefficients

*b*and

*c*, this case is directly comparable to the PDE (1.3). The infinitely hard pebble

**could represent either the average value of a very hard pebble population in which the relatively soft pebble**

*z***y**is being abraded or, alternatively, the average value of a pebble population which remains invariant not because every pebble remains invariant but because it has a constant source and a constant sink. In the latter interpretation

**y**is also a member of this population. We are interested in the geometric properties of these flows. Although the explicit formulae in (2.13) and (2.14) admit a rather straightforward, rigorous mathematical study, here we merely state the main geometric facts as observations.

As stated in §2, the flow (2.15) cannot have any genuine fixed points, and the only homothetic solution is the sphere *S*. Nevertheless, to explore the geometry of the flow we can still ask whether ‘temporary’ homothetic solutions satisfying (2.8) may exist; thus we seek the simultaneous solution of the quadratic equations
2.16Two of the four roots of (2.16) are strictly positive and two are strictly negative. The former define two surfaces *λ*_{i}(*y*_{1},*y*_{2},*b*,*c*) (*i*=1,2), where *b* and *c* are substituted from (2.12). Above these surfaces (*y*_{3}>*λ*_{i}) the flow moves away from *S*, towards *O*; *below* the surfaces (*y*_{3}<*λ*_{i}) the opposite happens; and *between* the surfaces the derivatives and have opposite signs. Points along the intersection lines defined by *λ*_{1}=*λ*_{2} satisfy (2.8); trajectories hitting these lines slow down and turn back at *cusps*. Since the (*y*_{1},*y*_{2}) projection of such trajectories *loiters* near the cusp, we call these points *loitering points*, and the surfaces *λ*_{1} and *λ*_{2} we call loitering surfaces. The distance between the two loitering surfaces is rather small, so most trajectories take a sharp U-turn. Both loitering surfaces escape to infinity at *O* and both have coinciding, global minima at *S*: . Based on this and by using (1.6) and (2.12) we obtain a global, sufficient condition for the attractivity of the sphere
2.17In the case of spherical objects with radius *y*_{3}=*R* being abraded by spherical abraders with radius *r* (*z*_{1}=1,*z*_{2}=1,*z*_{3}=*r*), we arrive at the stability result of Bloore (1.4). In general, (2.17) is a global extension of this result (in the frame of the box equations); this condition does not depend on the box ratios *y*_{1} and *y*_{2} of the abraded object. In addition to the previous observations, one can also show that trajectories of (2.15) transversally intersect each loitering surface *λ*_{i} at most once. Combining these facts, one can globally describe the steady-state flow (2.15). The (*y*_{1},*y*_{2}) projection of all trajectories of (2.15) are either of ‘type I’: *O*→*S*, or of ‘type II’: *S*→*S*, and they have at most one extremum in either variable, i.e. type II trajectories are ‘simple loops’ originating and ending at the sphere (see figure 2*a*). Although the PDE (1.3) may have a more complex structure, simple loops appear to be present. Figure 2*b* illustrates two trajectories originating at an ellipsoid with axis ratios *y*_{1}=0.83 and *y*_{2}=0.91 evolved under (1.3) and (2.15), respectively. Whereas steady-state flows correspond to collisional pebble abrasion among well-mixed pebbles and this is an uncommon geological situation, there exist observations (Nielsen 1983) that strongly suggest that non-monotonic evolution of box ratios may be present in some geological settings. A detailed description of this phenomenon is provided by Kuenen (1964).

We return to the experimental verification of loops in §3. We also remark that type I trajectories may correspond to the evolution of pebbles carried by rivers where the abrading environment is the riverbed, represented by a very hard and very large abrader.

### (b) The self-dual case

We now study the case when two identical pebbles **y** and ** z** mutually abrade each other, i.e. we have

**y**(

*t*)≡

**(**

*z**t*), so (2.13) and (2.14) are again reduced to a three-dimensional system 2.18which we call the

*self-dual flow*with no free parameters and an autonomous direction field in the (

*y*

_{1},

*y*

_{2}) plane 2.19where and are defined in (2.12) by a

**=**

*z***y**substitution and the superscript

*C*stands for ‘collisional’. Since (2.18) has time-dependent coefficients, it is difficult to compare it directly with (1.3). Nevertheless, it represents the abrasion of pebbles in an environment similar to itself, i.e. it is the simplest model of the abrasion of segregated pebble populations. The geometry of self-dual flows is simple (visually rather similar to the Gaussian flows illustrated in figure 1); in fact

*all*trajectories are of type I and converge globally to the sphere and in (2.18) we have 2.20on all trajectories, at all times. This is not very surprising, not only because it follows from (2.19) but also because it intuitively matches (2.17). The latter condition was derived for constant

**however, we can see that it is fulfilled in the self-dual flow on every trajectory, at any time. Since self-dual flows are the simplest model including segregation, they are closer to the physical process of pebble abrasion on shingle beaches than stationary flows. Whereas self-dual flows have an autonomous two-dimensional component as in (2.10), we showed that such a system might produce homothetic solutions; nevertheless, in the collisional box equations there is no such indication.**

*z*;## 3. Friction

### (a) Basic assumptions and a semi-local PDE model

From §2, we learnt that *collisional abrasion*, as modelled by (1.3), is dominated by Gaussian and mean curvature flows (2.4); thus it converges ultimately to the sphere. It is not entirely surprising that collisional abrasion is determined by local quantities since in this model we assume that abrading particles arrive from uniformly distributed directions. This model is based on the assumption that in the abrasion process the energy level is high enough to support the free flight of mutually abrading particles.

Frictional abrasion is rather different: here the kinetic energy can be rather low and the orientation (and thus the global shape) of the particle play a key role, so we look for a *semi-local* PDE model depending on both local and global quantities. Adopting Aristotle's approach, we regard the radial distance *R* (measured from the centroid) as our primary variable but we also assume frictional abrasion to depend on the global minimum and maximum of *R* described by the variables , , respectively
3.1

Thus our model may be regarded as a generalization of that of Aristotle (1.1). Friction influences abrasion either by sliding or by rolling, and we proceed by making specific assumptions about these components, starting with sliding. In shear band experiments with granular materials, it has been observed (Pena *et al.* 2008) that the long axis of elongated particles gets aligned with the direction of shear. This result indicates that, in the case of *sliding*, abrasion will be concentrated on the flat sides of the abraded pebble. This implies that for *sliding friction* we can assume that
3.2In the case of *rolling abrasion*, the situation is radically different. Rolling is essentially a two-dimensional motion, the plane of which is determined by the kinetic energy. If this energy is so high that the pebble can even roll over the point furthest from the centroid, then we are at the transition to collisional abrasion. There is also a lower energy threshold below which rolling is not possible: this corresponds to a trajectory with minimal maximum distance. Rolling occurring between these two critical energy levels results in a non-monotonic abrasion law, i.e. we assume that abrasion speed is small for radii to the minimum, zero at the maximum and larger in between, i.e. we assume the existence of a radius *R*_{☆}, , for which we have
3.3The non-monotonicity of the abrasion law could be explained by observing that surface points corresponding to very large radii are hardly affected, because they are not part of the typical trajectory. Points corresponding to very small radii are part of typical trajectories; however, they will abrade less than points with larger radii along the same trajectories. In other words: rolling is essentially a planar problem; the plane of rolling prefers smaller radii because of energy minimization; however, from among radii in the plane of rolling, larger radii will be abraded faster (in accordance with Firey's (1974) argument).

The above ideas can be translated into what we shall call a semi-local PDE model, i.e. we propose a definite form for *f* in (3.1). Here we introduce separate terms for sliding and rolling with independent coefficients *ν*_{s} and *ν*_{r}, respectively, and use the dimensionless ratios and
3.4For sufficiently high values of *n*, this model appears to capture the most essential physical features of the processes we are aiming to describe. The PDE model is clearly not unique but provides a simple basis for a qualitative analysis.

### (b) Unified model for collisional and frictional box flows

In the box approximations, assumptions (3.2) and (3.3) translate into
3.5respectively. In the limit, the PDE defined by (3.1) and (3.4) yields for the box flows
3.6which is equivalent to
3.7where
3.8These equations satisfy the conditions (3.5). Similarly to the self-dual flows in the collisional model, we can observe that, while the full flow is three-dimesnional, both component flows (sliding and rolling) can be reduced to simple planar flows
3.9and in (3.7) we have
3.10We can now simply add collisional and frictional flows (2.13)–(2.14) and (3.7) to obtain the unified equations (denoted by superscript *u*) for the two-body problem
3.11and
3.12

### (c) The steady-state flow in the presence of friction: experimental evidence for simple loops

We study the steady-state flows in the presence of friction, defined by the three-dimensional system
3.13The basic geometry is similar to the collisional case described in §2*a*: we can only observe type I (*O*→*S*) and type II (*S*→*S*) trajectories. To verify the latter, we conducted a simple table-top experiment in the spirit of Kuenen (1956). In each experiment, we rolled five rounded chalk pieces of approximately 9.8 mm diameter surrounded by approximately 400 hard, plastic balls of diameter 1.8 mm in a horizontal glass cylinder at 5 r.p.m. for approximately 24 h. The box dimensions *u*_{1},*u*_{2} and *u*_{3} were measured every hour. (*u*_{3} was the largest diameter, *u*_{2} the largest diameter orthogonal to *u*_{3}, and *u*_{1} was the diameter orthogonal to both *u*_{3} and *u*_{2}.) The data for the five chalk pieces were averaged. For the computations, we used the unified frictional equations (3.13) and the coefficients *b* and *c* were computed from (2.12). In the case of abrading plastic balls with radius *r*=0.9, we had *b*=*r*=0.9 and *c*=*r*^{2}=0.81, respectively. Initial conditions were set identical to the experiments and the friction coefficients *ν*_{S} and *ν*_{R} were selected to achieve a good match.

As one can observe in figure 3, both experiments and numerical simulations show fair quantitative agreement, with both displaying the simple loop predicted in §2*a*.

### (d) The self-dual flow in the presence of friction and the existence of stable, non-trivial attractors

Analogously to the self-dual collisional flows (2.18), we define the self-dual frictional flows as the mutual abrasion of two identical pebbles in the presence of friction
3.14which, similarly to (2.19), also have an autonomous two-dimensional direction field, but with two free parameters *ν*_{s} and *ν*_{r}. We look for homothetic solutions (2.7): as fixed points of (3.14), i.e. we solve , which is a 2×2 linear system for *ν*_{s} and *ν*_{r}, yielding
3.15where *F*_{i} are from (2.2) and are from (3.8). By substituting (*y*_{1}=*y*_{10},*y*_{2}=*y*_{20}) into (3.15), we obtain the friction coefficients in the presence of which (3.14) will have a homothetic solution at (*y*_{1}=*y*_{10},*y*_{2}=*y*_{20}).

Friction will produce physically observable, stable attractors if both critical parameters are positive. The condition for this is analogous to the necessary condition (2.9) defined for purely collisional flows. Here we regard the vector sum of the self-dual collisional flow with the frictional flows. To obtain equilibrium, the former must be ‘bracketed’ by the latter; this can be expressed based on (2.19) and (3.9) as (see figure 4*b*)
3.16Condition (3.16) is satisfied by our equations (3.7) and (3.8) for friction; in fact, the stricter condition
3.17is also met and therefore the proposed PDE (3.4) is capable of predicting the existence of non-trivial attractors in the (*y*_{1},*y*_{2}) plane. Since both *F*_{1} and *F*_{2} appear to be monotonically decreasing in *y*_{1} and *y*_{2}, the fixed points defined by (3.15) are stable attractors, illustrated in figure 5*a*.

## 4. Stochastic process: multi-body simulations

The unified collisional–frictional equations (3.11) and (3.12) describe the mutual abrasion of two individual pebbles; however, they also define a *Markov process* where **y** and ** z** are random vectors with

*identical*distributions since they represent two random samples of the same pebble population. The evolution of this Markov process (and thus the time evolution of the pebble size and ratio distributions) is of prime interest since it determines the physical relevance of the stable attractors identified in §3

*d*. Whereas the analytical investigation of the Markov process is beyond the scope of this study, direct simulations are relatively straightforward. We consider

*N*pebbles out of which we randomly draw two with coordinates

**y**

^{0}and

**z**

^{0}and run equations (3.11) and (3.12) for a short time period Δ

*t*on these initial conditions to obtain the updated vectors

**y**

^{1}and

**z**

^{1}. In the simplest linear approximation we have the recursive formula 4.1and 4.2Such an iterative step can be regarded as the cumulative, averaged effect of several collisions between the two selected pebbles. Apparently, the

*N*=2,Δ

*t*→0 case is identical to (3.11) and (3.12). Previously, we investigated the behaviour of the deterministic flows in the special cases of steady-state and self-dual flows. Multi-body simulations allow the numerical study of the statistical stability of the flows, i.e. one can assess the stability of the aforementioned special cases.

### (a) Stochastic, multi-body simulation of collisional flows

The characteristic U-turn identified in the continuous flow for type II trajectories appears to be very robust. In multi-body simulations one can observe the fuzzy zigzag geometry of the trajectory near the turning point, indicating that a larger amount of time is spent near the loitering points. This behaviour is illustrated in figure 6*a* for *N*=4. The self-dual collisional flow (2.18) proves to be remarkably stable and in multi-body simulations the trajectories of the continuous flow are often followed rather closely. However, there are exceptions. Figure 6*b* illustrates an *N*=9 body simulation with random initial conditions, where eight trajectories converge to the sphere (similar to the continuous flow); however, one trajectory turns back. This deviation from the continuous case can be expected in the vicinity of the sphere, since here the derivatives and converge to zero and so small perturbations can change their signs.

### (b) Stochastic, multi-body simulations of the unified collisional and frictional model

In §3*d*, we showed that, in the presence of friction, stable, non-trivial attractors appear in the self-dual flows (3.14). Whereas these fixed points collectively attract the abrading objects, the latter mutually repel each other, as we illustrate in figure 7, where two almost identical pebbles initially follow the same trajectory; however, at some critical time they split. Here we used (4.1) and (4.2) with *N*=2 and *almost* identical initial conditions to simulate (3.14).

This shows that self-dual flows (describing the mutual abrasion of a homogeneous pebble population consisting of identical pebbles) are structurally unstable in the sense that the homogeneity of initial data decays (pebble ratios may mutually repel each other in the abrasion process). Based on this, we expect that in multi-body random simulations according to (4.1) and (4.2) the attractors do not appear. The geological motivation for the study of self-dual flows was the observation that pebbles are often segregated by size by the wave transport. If this process is happening on the same (or faster) time scale as abrasion, then the self-dual flows and their attractors may be stabilized by sustaining the homogeneity of the initial data, i.e. by selecting only nearly similar pebbles for the collision process. This could be implemented by introducing a correlation *r* between the random vectors **y** and ** z**; our observation of the instability corresponds to the uncorrelated

*r*=0 case. The fully correlated

*r*=1 case is essentially identical to the deterministic evolution of (3.14). Another option to implement segregation (possibly closer to the geological process) is to omit pebbles from the simulation the maximal size

*y*

_{3,i}of which differs from

*all other*pebble sizes

*y*

_{3,j}by a prescribed factor. For each pebble, the coefficient 4.3was computed and pebbles with

*σ*

_{i}>

*σ*were omitted. Increasing segregation effect corresponds to decreasing values of

*σ*; this is illustrated in figure 8. Observe that, for weak segregation

*σ*=2, the stable fixed point of the deterministic flow vanishes, while for stronger segregation it persists.

## 5. Conclusions

Our goal was to create a framework in which the interaction of mutual particle abrasion with global transport and segregation can be studied.

### (a) Mutual abrasion versus individual abrasion

As a first step, we radically simplified Bloore's classical PDE, describing the general collisional abrasion process of an individual pebble. Our system of ODEs, called the box equations, attempts to track the evolution of maximal pebble size and the main pebble proportions. Whereas the box equations are merely a heuristic simplification of the PDE describing individual abrasion, they admit the global, analytic study of *mutual abrasion*. To achieve this, we incorporated the results of Schneider & Weil (2008) and Várkonyi & Domokos (2011), providing the general interpretation of the coefficients in Bloore's PDE. We arrived at the system of six-coupled ODEs, defining the mutual abrasion of two pebbles. We investigated two special cases of this rich process when the six-dimensional system reduces to a three-dimensional system: the steady-state flows where one pebble is regarded as constant and the self-dual flows where both pebbles and abraders are treated identically.

### (b) Special cases: steady-state and self-dual flows

Steady-state flows may serve as a skeleton model of abrasion of *well-mixed* pebble populations where the abrading environment has stationary size distributions. In the box equations, we extended Bloore's local stability result (1.4) for nearly spherical objects abraded by spherical objects to obtain the sphere's global attractivity criterion. Equation (2.17) provides a sufficient, global criterion for the convergence to the sphere, depending only on the size ratio of the abraded and the abrading object and the box proportions of the latter. We found that the steady-state box flows admit only two types of trajectories, one of which appears as *simple loops* in the plane of the two box proportions, originating and ending at the sphere. We confirmed the existence of this ‘boomerang effect’ by direct simulation of the PDE, in simple table-top experiments and in stochastic simulations. Whereas loops may cause some statistical accumulation of pebble ratios in the vicinity of the turning point (which we called ‘loitering points’), they certainly do not account for the very clear and striking appearance of dominant pebble ratios.

We found that the self-dual flows, modelling the abrasion of pebbles in an environment identical to the abrading pebble, uniformly and globally converge to the sphere. These flows represent the box equation model for pebble abrasion in segregated environments.

### (c) Friction

Since none of the collisional equations accounted for the existence of homothetic solutions, we added frictional terms, corresponding to sliding and rolling, with coefficients *ν*_{s} and *ν*_{r}, respectively. We postulated a semi-local PDE (3.4) based on simple physical observations and derived box equations describing the unified, self-dual collisional and frictional process. We found that even an arbitrarily small amount of friction can stabilize non-trivial shapes as homothetic solutions in the vicinity of the sphere. We also observed that sliding friction stabilizes discoid shapes, whereas rolling friction stabilizes elongated rod-like shapes. At larger values of the coefficients *ν*_{s} and *ν*_{r}, shapes further away from the sphere can be stabilized. Physically, the increase of the coefficients might be due either to the increased importance of frictional efficiency or, more importantly, to the decrease of total energy of the abrasion process. In the box equations, we determined the coefficients *ν*_{s} and *ν*_{r} corresponding to fixed points with given locations. The latter data could be extracted from statistical measurements on shingle beaches and our equations could be used to identify the associated coefficients. Our proposed semi-local PDE model is certainly not unique; alternative PDEs, also meeting the same criteria, can predict multiple attractors for the same value of the coefficients; therefore, several dominant box ratios could co-exist in a pebble environment governed by mutual abrasion by collisions and friction. As an example we mention the following model as an alternative to (3.4)
5.1which, similarly to (3.4), meets the criteria (3.2) and (3.3); the corresponding box flow (defined in the limit) is illustrated in figure 5*b* with a saddle-type critical point.

### (d) Statistical approach: segregation by transport

Whereas the existence of attractors in the presence of friction is promising, the model calls for statistical verification. Our box equations in the self-dual case describe the mutual abrasion of two identical pebbles, i.e. the abrasion process in which the abraded pebble and abraders are treated completely symmetrically. To analyse the statistical stability of the results, we described the Markov process based on the box equations. This process defines the evolution of pebble size and ratio distributions under the combined action of collisional and frictional abrasion and segregation. We ran direct simulations starting with *almost identical* initial conditions. In this process, the colliding pair of pebbles is drawn randomly from the same distribution. First we modelled abrasion without transport and segregation. To this end, the draw was uncorrelated and we found that the attractors disappeared: this matched the geological observations of Carr (1969), Landon (1930) and Bluck (1967) claiming that global transport played an important role in the emergence of equilibrium shapes. Next, we introduced transport and segregation by size, so the maximal size of the colliding pairs was highly correlated and we found that stable attractors were sustained in the random process. Our model thus predicts that dominant pebble ratios as stationary equilibrium shapes may robustly appear in abrasion processes where, in addition to collisions, friction also plays an important role, and an ongoing global transport process segregates pebbles based on size. Segregation by shape (also included in our model) would enhance the stability of these dominant ratios even further. These predictions match well geological observations; also, they confirm Aristotle's intuition that the basic mechanisms of pebble abrasion can only be understood by including semi-global and global effects.

## Acknowledgements

The authors thank András A. Sipos for his valuable help with the numerical simulation of the PDE (figure 2) and Ottó Sebestyén for builiding the equipment and conducting the chalk experiments. G.D. was a Visiting Fellow Commoner at Trinity College, Cambridge, during part of this collaboration. This research was supported by the Hungarian National Foundation (OTKA) grant K104601.

## Footnotes

↵1 Son of the Nobel Prize winner Lord Rayleigh.

- Received September 14, 2011.
- Accepted April 12, 2012.

- This journal is © 2012 The Royal Society