Time Scales for Rounding of Rocks through Stochastic Chipping

Various types of mechanisms act to weather or transform the shapes of rocks by eroding away their volume. Many of these processes, such as collision induced fragmentation, are inherently stochastic. Salient questions related to the emerging shapes as more and more volume is removed by random interactions with neighboring stones include the degree to which rock surfaces become smooth on a macroscopic scale, as well as the extent to which the accumulating random fractures yield spherical profiles, and the rate at which mass is removed through weathering. We investigate these concerns with large scale Monte Carlo simulations in which primordial convex polyhedral rocks (i.e. proto-clasts) with a small number of faces are subject to a sequence of planar fractures which successively carve away material from the parent stones.

With a combination of quantitative measures of the deviation from a perfect spherical shape, the overall oblateness and prolateness, and smoothness in terms of how evenly the total surface area is distributed over the individual facets, we characterize the evolution of the shapes of ensembles of rocks as stochastic fractures accumulate and volume is worn away. To statistical Monte Carlo statistical error, we follow the shape trajectories of at least 1000 stones, which accumulate as many as $2\times 10^{6}$ planar slices, yielding polyhedra with up to 3,600 facets on average.

Using a two-pronged approach, we find that stones subject to stochastically driven collisional weathering ultimately converge to spherical shapes. On the one hand, we consider a steady state scenario where observables (apart from the volume, which decreases monotonically in time) stabilize at constant values and thus in a sense equilibrate. By tuning a parameter $\gamma$ controlling the coarseness of the shapes at equilibrium (i.e. the mean relative area $\langle\Delta A\rangle/A_{\Sigma}$ of newly formed facets), we find that the equilibrium mean number $\langle n\rangle$ of planar faces scales as $\gamma$ while quantitative measures of the deviation from spherical shapes and/or a perfectly smooth surface tend systematically to zero with increasing $\gamma$. On the other hand, we also calculate the time dependence (i.e. with the number of sustained slices $N_{\mathrm{sust}}$ serving as a proxy for time) of observables in the scenario in which the latter eventually equilibrate. We use the latter to obtain corresponding quantities for a broad range of non-steady state scenarios, a technique we validate with direct Monte Carlo simulation in the context of the fixed velocity scheme, envisaged for rocks impelled by a current at a mean constant speed independent of their mass. In this way, we glean quantitative results for time scales for the attainment of milestones such as the removal of a specific fraction of the primordial volume, finding that such times vary depending on the specific erosion scenario under consideration but do not diverge, a point we underscore with approximate closed form expressions also yielding finite times while serving as upper bounds for simulation results.

Taking into consideration that in practice cohorts of primordial stones prior to weathering may be structurally diverse, we consider protoclasts in the form of irregularly shaped polyhedra to account for aggressive fragmentation events prior to the collision induced weathering. To mimic structures resulting from dramatic fragmentation processes, initial cohorts of structurally disordered shapes are generated with a sequence of random fracturing planes, invariably accepted independent of the orientation and depth of the prospective slice. With protoclast cohorts generated in this manner, we find relevant observables to vary smoothly with time. As companion calculations, we also examine structurally homogeneous primordial stones, including regular cubes, which ultimately are smoothed into spherical shapes just as occurs in the case of irregular protoclasts. However in contrast to structurally disordered stones, observables in the case of cubes exhibit two continuous phase transitions with all of the hallmarks of a second order phase transition, including power law scaling for singularities in relevant observables as we show explicitly in this work. The first transition is a structural transformation in which the primordial cube facets are sheared away, while in a subsequent transition the stones derived from regular solids revert to spherical shapes with concomitant singularities in quantities sensitive to deviations from a perfect spherical profile. A structural phase transition of the former type has been proposed in previous studies based on laboratory experiments and numerical studies Domokos2 ; Domokos3 . We apply the tools of finite size scaling to locate the transition and calculate associated critical exponents.

In addition to regular shapes, we consider a homogeneous ensemble of geometrically irregular solids finding asynchronous transitions for the erosion of primordial facets instead of a single structural transition for the solid as a whole, though again we find a single second order transition to a round spherical shape at a later time after all of the primordial facets have been removed. Thus, only in the case of highly symmetric shapes is there a simultaneous elimination of the facets inherited from the parent solid.

In spite of the stochastic nature of the erosion scheme, we find in the regime that fragments cleaved away through random fracture events are small in volume relative to that of the stone as a whole that shapes evolve deterministically, with reproducible changes on scales small in comparison with the overall size of the stone but large relative to freshly exposed faces. As in a recent study, we find a universal dependence of observables on the remaining volume fraction Domokos1 . We take advantage of this characteristic to calculate time dependent observables for a variety of erosion schemes using results from the scenario of fixed relative area for typical fresh facets created by fracture events.

Aside from measures of deviation from a perfectly spherical shape, we also examine more specific macroscopic structural characteristics, including the degree to which stones are oblate or prolate as volume is chipped away. While both the former and the latter tend to zero as rocks are rounded, the oblateness measure is non-monotonic in time, initially rising and then peaking after approximately 50% of the original volume has been removed, and declining and tending to zero thereafter.

As a novel feature of the calculations in this work, no restrictions are imposed as to how many vertices and facets are truncated or removed by a sustained slice, and our calculation is compatible with the elimination of an arbitrary number of vertices per fracture. For a plausible physical basis for fracture events, measures such as the area of a newly exposed face are considered to determine if a candidate fracture plane is accepted. To set the scale of new facet areas, we implement an energy based criterion for accepting prospective slices, taking the kinetic energy input needed to cleave away a slice to be proportional to the area of the new face. In this vein, in the spirit of statistical mechanical treatments, we assume the fracture probability to be exponentially suppressed in the exposed area $\Delta A$, and given by $P(\Delta A)=e^{-\gamma\Delta A/\theta}$ where $\gamma$ is a constant related to the toughness of the material comprising the rocks and $\theta(V)$ is a volume dependent measure of the mean kinetic energy or “temperature”, with volume serving as a proxy for the mass in the case of stones of uniform composition as assumed in this manuscript. In addition, since we consider initial cohorts of stones to be mono-dispersed with respect to mass (though structurally diverse in the case of irregular protoclasts), we operate in terms of the relative volume $\tilde{v}\equiv V/V_{0}$ with $V_{0}$ the volume prior to the erosion process. With analytical arguments and direct simulation results, we show that characteristic times scale quadratically in the toughness parameter $\gamma$.

In this work to be definite we assume a power law dependence $\theta=\theta_{0}\tilde{v}^{\alpha}$ (i.e. with $\alpha>0$), reflecting kinetic energy dependence on volume as a rock is borne along (e.g. by water currents in a river or stream as collisions chip away material). Defining $\gamma^{{}^{\prime}}\equiv\gamma/\theta_{0}$ for the sake of brevity, the acceptance probability for a prospective fracture plane is $P(\Delta A)=e^{-\gamma^{{}^{\prime}}\Delta A\tilde{v}^{-\alpha}}$. Assuming $\gamma^{{}^{\prime}}\langle\Delta A\rangle\tilde{v}^{-\alpha}\sim 1$ for the mean area $\langle\Delta A\rangle$ of a new facet, one would argue on dimensional grounds that the mean number of faces is $\langle n\rangle\sim A_{\Sigma}/\langle\Delta A\rangle=\gamma^{{}^{\prime}}\tilde{v}^{-\alpha}A_{\Sigma}$, where $A_{\Sigma}$ is the total polyhedron surface area. For the 3D case, taking the total area to scale as $\tilde{v}^{2/3}$, $\alpha_{c}$ marks a boundary between stones which in principle could become smoother over time ($\alpha>\alpha_{c}$) and shapes which become coarser ($\alpha<\alpha_{c}$) as fractures accumulate; whereas relative areas $\langle\Delta A\rangle/A_{\Sigma}$ of new facets decrease for the former, newly exposed facets encompass an increasingly large share of the surface in the case of the latter. The boundary case $\alpha=\alpha_{c}$ thus offers the possibility for the attainment of a steady state for the mean relative area of new faces, as well as other variables of interest.

Another merit of obtaining steady state or relative area results is the fact that the universal dependence of observables on volume fraction permits one to map remaining stone volumes onto time for a scheme distinct from the relative area case, an technique bolstered with results from large-scale Monte Carlo simulations. We show in this manuscript that both the former and latter are identical up to random Monte Carlo statistical error. Time scales obtained in this manner are then used to validate a closed form expression serving as an approximation and strict upper bound for times associated with the chipping away of a specific volume fraction or the structural phase transition for cohorts of regular protoclasts in which primordial facets disappear.

In Section II, we discuss principal quantitative measures of deviation from sphericity as well as key methodological elements. Section III contains results for rocks in the case of the steady state scenario, while Section IV broadens the discussion to non-steady state situations. In Section V, we calculate time scales for a range of erosion schemes, with a closed form expression obtained as an upper bound for time scales for salient structural milestones. Conclusions are found in Section VI.

Two salient quantities which permit one to address in a direct manner the degree to which a shape departs from a perfect spherical profile include the configuration averaged sphericity $\phi_{3}$ and the ratio $r^{\mathrm{min}}_{\mathrm{max}}$ of the minimum and maximum distance of surface points for the center of mass. Both measures distinguish among a spherical shape and a stone which is not perfectly round. The sphericity benefits from the fact that surface area to volume ratio of a solid attains a global minimum for spheres, and is given by $\phi_{3}\equiv[6\pi^{\frac{1}{2}}V]^{\frac{2}{3}}/A_{\Sigma}$ Wadell , where $A_{\Sigma}$ is the surface area and $V$ the volume. Since $\phi_{3}<1$ (except in the case of a perfect sphere where $\phi_{3}=1$), the complement $1-\phi_{3}$ serves as a measure of the departure from a perfectly spherical shape. An alternative measure, sensitive to local features, is $r^{\mathrm{min}}_{\mathrm{max}}\equiv d_{\mathrm{min}}/d_{\mathrm{max}}$, the ratio of the minimum and maximum distances, respectively, of points on the rock surface from its center of mass. In the context of our calculations involving faceted objects, $d_{\mathrm{min}}$ is the distance to the closest planar facet and $d_{\mathrm{max}}$ is the distance to the farthest vertex. Again, $r^{\mathrm{min}}_{\mathrm{max}}$ peaks at unity only in the case of a perfect sphere, with $1-r^{\mathrm{min}}_{\mathrm{max}}$ serving as a measure of the deviation from a perfect spherical shape.

Related to whether there is a systematic convergence to a spherical shape is if erosion through the stochastic chipping mechanism yields a stone with a smooth surface. Although sphericity implies smoothness, the converse is not true. In fact, in this work we find an example of a perfectly smooth non-spherical shape at the phase transition for structurally homogeneous cohorts where primordial facets disappear. Hence, as a measure of smoothness not anchored to a specific geometry, we use the Inverse Participation Ratio (IPR), where $\textrm{IPR}=\langle A_{\Sigma}^{-2}\sum_{i=1}^{N}A_{i}^{2}\rangle$ with $A_{\Sigma}$ being the total surface area, $A_{i}$ the area of the $i$th of $N$ facets, and angular brackets indicated a configurational average. The IPR tends to zero with increasing $N$ for an even distribution of the area of the polyhedron faces (i.e. for smooth shapes), but converges to a finite value if a small number of facets contain a macroscopic fraction of $A_{\Sigma}$ (i.e for rough stones), and thus serves as a quantitative measure of smoothness. Appealing to the IPR in this manner is analogous to its usage in the context of charge transport discussions to distinguish among extended and localized carrier states Wegner .

In this work, we simulate collisional weathering with large scale Monte Carlo simulations involving sequences of randomly oriented planar fractures, where rocks begin as regular or irregular protoclasts with a small number of faces [i.e. 6 and a mean of 9.03(1) sides for cubic and irregular polyhedra respectively], ultimately yielding polyhedra with as many as 3,600 facets. Although there are qualitative differences for the time dependence of observables in the cases of cube shaped and irregular protoclasts (e.g. the structural phase transition for cube shaped protoclasts not seen in ensembles of irregular shapes), both cohorts ultimately converge to spherical shapes in the long time regime. Moreover, we obtain an approximate analytical expression for time scales for salient stages in the attainment of a perfectly spherical shape. The fact that these times are finite, bounding corresponding numerical results from above is compatible with the eventual conversion of rocks into spherical shapes for a wide range of erosion scenarios.

In stochastic chipping sequences, each sustained slice exposes a new face while truncating one or more vertices. The likelihood of accepting a prospective slice hinges solely on the area $\Delta A$ and volume of the parent solid with no restriction on the number of truncated vertices, raising the possibility of the elimination of an entire face or faces if each of the constituent vertices are cleaved away. For the sake of an efficient implementation, updates of lists of vertices and planes making up the polyhedron include the addition of new features as well as the deletion of those eliminated by the most recent fracture event. Indeed, in the steady state scenario, as many vertices and facets are removed on average as accumulate with each sustained slice. Previous studies in two and three dimensional geometries have involved the truncation of a finite number of vertices or facets per fracture event Durian2 ; Krapivsky ; Domokos2 ; Domokos1 . A calculation involving two dimensional shapes allowed for the removal of an arbitrary number of vertices Durian1 To our knowledge, our calculations are the first to examine erosion phenomena due to a sequence of planar fractures allowing for the elimination of an arbitrary number of vertices and/or facets for three dimensional geometries.

To facilitate the examination of polyhedra with a large number of faces, a variety of measures are employed to optimize the efficiency of locating vertices of prospective facets. While one could in principle examine all possible intersections of a fracture plane with the faces comprising the stone, in general only a small subset of vertices identified in this manner populate the new face, comprising a share of the total which diminishes as the number of planar faces increases. We avoid considering spurious vertices by only examining faces which contain vertices both above and below the prospective fracture plane, as facets not meeting this condition are either sheared away altogether or are not truncated by slice at all and hence do not yield any new constituent vertices. As a further constraint, we only seek intersections of the slicing plane with pairs of faces which have an edge in common. Finally, we avoid computational costs associated with slices unlikely to be accepted by considering a sphere inscribed in the polyhedron and concentric with its center of mass, for which circular regions exposed by a fracture plane bound the area of a prospective face from below, with the corresponding acceptance probability overestimating that of the prospective new facet. In this manner, we avoid consideration of fracture planes for which the acceptance probability is no greater than $10^{-8}$, events with negligible incidence over the course of an erosion sequence. Operating in this manner, the computational burden of a sustained slice scales no more rapidly than the number of facets of the solids we consider.

The areas of individual facets and the total volume are key quantities in finding the probability of a prospective fracture event as well as many of the observables we calculate. For individual faces, the arithmetic mean of each of the member vertex coordinates provides a convenient interior point for subdividing the polygonal facet into triangles, whose areas are calculated and summed to yield the combined face area. On the other hand, to find the total volume of the stone one operates in an analogous fashion, partitioning the polyhedron into constituent tetrahedra whose combined volume is the shape total volume. To this end, we again subdivide faces into triangles; the three vertices of the latter along with an interior point (given by the mean of all of the polyhedron member vertices) define the component tetrahedron.

Complementary to the sphericity measures and the IPR are quantities which bear more specifically on the shape of the stone, namely parameters which characterize the degree to which a stone is oblate or prolate. Measures of the latter are gleaned from the principle moments of inertia of the equivalent ellipsoid, obtained by diagonalizing the moment of inertia tensor of the stone. As in the calculation of total volume, the polyhedron is partitioned into constituent tetrahedra with elements $I_{ij}$ of the moment of inertia tensor given by Tonon

where $V_{\mathrm{tet}}$ is the volume of the constituent tetrahedron, while $s_{i}\equiv\sum_{l=1}^{4}x_{li}$ and $t_{ij}\equiv\sum_{l=1}^{4}x_{li}x_{lj}$ with, e.g., $x_{li}$ being the $i$th component of the $l$th vertex of the tetrahedron. Combining moments for each component tetrahedron yields the moment of inertia tensor for the shape as a whole. To find the moment of inertia relative to the center of mass, as appropriate for a bona fide ellipsoid, we appeal to the Steiner parallel axis theorem Marion while taking advantage of the fact that the center of mass for a tetrahedron of uniform composition is the arithmetic mean of its four vertices to find the center of mass of the entire shape. In terms of semi-major axes $\{a,b,c\}$ of the equivalent ellipsoid (i.e. with $a>b>c$), the oblateness and prolateness measures $\psi_{\mathrm{O}}$ and $\psi_{\mathrm{P}}$ are taken to be $\psi_{\mathrm{O}}=(b-c)/a$ and $\psi_{\mathrm{P}}=(a-b)/a$, respectively.

For the $\alpha=\alpha_{c}$ case where the acceptance probability for a prospective slice is $e^{-\gamma\Delta A/A_{s}}$ ($A_{s}=[6\pi^{\frac{1}{2}}V]^{\frac{2}{3}}$ being the area of the volume equivalent sphere), a steady state may in principle develop with the mean number of faces $\langle n\rangle$ converging to a fixed value since $\langle\Delta A\rangle\sim A_{\Sigma}/\gamma$ per the acceptance probability relation. Images of stones sampled from the steady state ensemble appear in Fig. 3, showing the emergence of smooth spherical shapes for toughness parameters ranging from $\gamma=10$ to $\gamma=2000$.

A salient question related to the attainment of equilibrium is how long is needed (in terms of the number of sustained slices $N_{\mathrm{sust}}$) for an ensemble of stones to reach steady state. In the context of numerical simulations, we insist that all observables of interest remain constant with respect to doubling of $N_{\mathrm{sust}}$. For $\gamma=100$, 3,000 sustained slices satisfies this criterion. As we now argue from geometric considerations, for $\gamma\gg 1$ (i.e. certainly for $\gamma\geq 100$) time scales for the removal of equivalent fractions of the original volume of an ensemble of rocks (and hence salient stages such as the achievement of steady state) scale quadratically with the toughness parameter $\gamma$.

Although the acceptance probability sets the scale $\langle\Delta A\rangle$ for the typical area of a new facet, fractures leaving behind faces with the same area need not cleave away the same volume. One anticipates more volume to be removed from sharply peaked features than from regions with less curvature by the same number of sustained slices Jerolmack1 . Thus the edges and vertices of the primordial polyhedron, are quickly worn down in the initial stages of the erosion process.

In the $\gamma\gg 1$ regime, fracture planes are constrained to be shallow, unable to penetrate deeply due to the exponential penalty on $\Delta A$ in the acceptance probability. We introduce the cleaving plane, just below the rock surface, as well as the parallel plane of tangency with a single point of contact with the surface, while superimposing Cartesian axes with the $xy$ plane coinciding with the tangent plane and the $z$ axis directed toward the interior of the stone. Hence, for the rock surface one would generically have have $f(x,y)\approx Ax^{2}+2Bxy+Cy^{2}$ near its minimum at the point of tangency, where $A$, $B$, and $C$ are locally determined constants. Seeking to eliminate the diagonal term, we have in matrix form

where $\vec{v}_{\mathrm{plane}}$ is the component in the tangent plane. The eigenvalues of the symmetric matrix are: $\lambda_{1,2}=\frac{1}{2}[(A+C)\pm\sqrt{(A-C)^{2}+4B^{2}}]$ with $\lambda_{1,2}>0$ if $AC>B^{2}$, as must be the case for the convex objects we consider. In terms of the orthonormal eigenvectors $\hat{u}_{1}$ and $\hat{u}_{2}$, one may write $\vec{v}_{\mathrm{plane}}=\alpha_{1}\hat{u}_{1}+\alpha_{2}\hat{u}_{2}$ and thus

The cleaving plane hence exposes an elliptical region for which $z=\lambda_{1}\alpha_{1}^{2}+\lambda_{2}\alpha_{2}^{2}$ for a fixed depth $z$ beneath the point of tangency, where the area is $A(z)=\pi ab=\pi z/l$ with $l\equiv(\lambda_{1}\lambda_{2})^{-1/2}$ being a locally defined length scale. Thus, the depth of a typical fracture plane (i.e. with an area $\langle A\rangle$) is $D=\langle\Delta A\rangle(\pi l)^{-1}$, and the typical volume $\Delta V$ removed is $\Delta V=\int_{0}^{D}A(z)dz=\langle\Delta A\rangle^{2}/(2\pi l)\sim\gamma^{-2}A_{\mathrm{\Sigma}}/(2\pi l)$ since $\langle\Delta A\rangle\sim A_{\Sigma}\gamma^{-1}$ for the relative area scenario. Thus, time scales for the removal of volume scale as $\gamma^{2}$ in the steady state scenario and more broadly as the square of the toughness parameter in non-steady state schemes as well.

Although we discuss a variety of non-steady state schemes, choose the constant velocity scenario (e.g. for stones borne along at constant speed in a stream or river) for direct Monte Carlo simulations. With kinetic energy scaling as $\frac{1}{2}mv^{2}$, and $m=\rho V$ in the case of stones of uniform composition, typical energy would be proportional to the volume. That the $\alpha=1$ exponent for the fixed velocity scenario exceeds the critical $\alpha_{c}=2/3$ for the relative area scheme implies an ever diminishing relative area with time for typical sustained slices and a concomitant rise in the mean number of facets. Nevertheless, although the relative area and fixed velocity scenarios are distinct, the universal dependence of stone profiles on fractional remaining volume Domokos1 dictates that key structural milestones occur at finite values of the reduced time $\tilde{\tau}$, borne out in the case of the fixed velocity scheme.

Fixed velocity erosion trajectories are calculated for irregular protoclasts with 2000 distinct realizations of disorder. A range of toughness parameter values from $\gamma^{{}^{\prime}}=20$ to $\gamma^{{}^{\prime}}=160$ are considered with the mean number of facets $\langle n\rangle$ exceeding 4000 for the latter. Representative examples of relevant variables are shown in Fig. 22 with the sphericity complements $1-\phi_{3}$ and $1-r_{\mathrm{max}}^{\mathrm{min}}$ and the IPR shown in panel (a) and (b) while oblateness and prolateness measures are shown in panel (c) with the primordial facet survival probability $f_{\mathrm{sur}}$ plotted in panel (d). As in the case of the time dependent relative area results, $f_{\mathrm{sur}}$ and ellipticity measures are converged with respect to the $\gamma^{{}^{\prime}}$ parameter with minor discrepancies for the oblateness measure near its maximum likely due to Monte Carlo statistical error. Quantities such as measures of deviation from a spherical shape and the IPR do diverge on the semilogarithmic scale, but at later times with increasing $\gamma$ suggesting convergence with respect to the latter. In a qualitative sense, in comparison to results in the case of the relative area scenario, fixed velocity variables evolve more rapidly initially with structural change slowing as sustained slices accumulate. This trend is compatible with the typical area of freshly exposed facets decreasing relative to the total area with decreasing remaining volume fraction $\tilde{v}$.

An additional merit of discussing the steady state scenario in the equilibrated regime is its relevance to non-steady state scenarios after a sufficient number of fractures have accumulated. We define an instantaneous or “effective” $\gamma$ with $\gamma_{\mathrm{eff}}=A_{\Sigma}/A_{\textrm{char}}$ where $A_{\textrm{char}}$ is the characteristic area exposed by a fracture event obtained by setting the exponent in the prospective slice acceptance probability to 1 and solving for the area. In the case of the fixed velocity erosion scheme, one has $\gamma_{\mathrm{eff}}=\tilde{v}^{-\frac{1}{3}}(6\pi^{\frac{1}{2}})^{-1}\gamma^{{}^{\prime}}$, which increases with decreasing remaining volume fraction. For representative variables of interest in Fig. 23, we show on the same plots equilibrated observables from steady state calculations and fixed velocity simulation results (solid or broken curves) plotted with respect to $\gamma_{\mathrm{eff}}$. The mean number of sites normalized with respect to $\gamma_{\mathrm{eff}}$, shown in panel (a), ultimately saturates at the large $\gamma$ equilibrium value of $1.82$. The IPR and complements $1-\phi_{3}$ and $1-r_{\mathrm{max}}^{\mathrm{min}}$ obtained from fixed velocity scheme simulations decrease monotonically with time in panels (b), panel (c), and panel (d) of Fig. 23, converging with steady state equilibrium results (open or filled circles). The merging of the fixed velocity simulation and equilibrium results suggests that even in non-steady state situations, stones reach a condition of quasi-equilibrium, in the sense that variables correspond closely with the steady-state counterparts calculated for $\gamma_{\mathrm{eff}}$ after enough time has elapsed.

Due to the universal dependence of relevant variables on volume fraction $\tilde{v}$, to find time scales for a salient event such as a structural phase transition, one need only determine the volume fraction (e.g. in the context of the relative area scenario) for the case of interest and calculate the elapsed time in terms of sustained slices using Equation 8. Since $f(\tilde{v})$ decreases monotonically, ultimately converting to the asymptotic value $f_{0}\equiv 0.216$, one has $f(\tilde{v})\leq f_{0}$ and thus time scales gleaned from Equation 8 are approximately given and bounded from above by $f_{0}\int_{\tilde{v}}^{1}A_{\mathrm{char}}^{2}\tilde{v}^{\frac{1}{3}}d\tilde{v}$, often amenable to exact calculation. For an power law scaling of the kinetic energy where $A_{\mathrm{char}}=(\kappa/\gamma)\tilde{v}^{\alpha}$, we obtain

Fig. 27 shows salient reduced time scales $\tilde{\tau}$ for a range of values of the exponent $\alpha$; solid curves are obtained from Monte Carlo simulation results in the case of the relative area scenario, while broken lines are an upper bound for time scales from Equation 9. In panel (a) of Fig. 27, pertaining to cube shaped protoclasts, the two time scales shown include the structural phase transition involving the loss of primordial facets as well as the attainment of spherical shapes (corresponding to equilibration in the relative area scenario) for volume fractions $\tilde{v}=0.57$ and $\tilde{v}=0.37$, respectively. On the other hand, in panel (b) of Fig. 27, for cohorts of irregular protoclasts, the three time scales shown correspond to $\tilde{v}=0.50$, $\tilde{v}=0.25$, and $\tilde{v}=0.10$. In spite of the monotonic rise in $\tilde{\tau}$ for the attainment of various volume fraction milestones, panel (a) and panel (b) of Fig. 24 do not diverge for a finite value of $\alpha$. Moreover, the upper bound in Eq. (i.e. the same upper bound for the cases of cubic and irregular protoclasts), albeit increasing with increasing $\alpha$, nevertheless remains finite for finite $\alpha$.

In conclusion, with large scale Monte Carlo simulations, we have examined the erosion of rocks through stochastic chipping, considering polyhedral stones with as many as 3,600 facets and averaging over at least 1000 realizations of disorder. Using an energy based criterion for accepting a randomly chosen slicing plane, our calculation is unique in placing no restrictions on the number of vertices and faces sheared away with each fracture event. We have argued on theoretical grounds and verified by direct simulation that time scales for the removal of a specific amount of volume or the attainment of a particular structural milestone scale quadratically in the toughness parameter $\gamma$ which specifies the amount of energy per area associated with a fracture; the largest time scales examined in our calculations exceed a million sustained slices per erosion sequence.

Cohorts of stochastically chipped protoclasts in the form of regular polyhedra undergo a structural phase transition in which all primordial facets are abruptly lost with concomitant singularities in other relevant observables, marking a genuine second order phase transition; in a subsequent structural transformation evident in measures of sphericity, the stones revert to spherical shapes. More broadly, however, ensemble averaged observables in the case of initially identical irregular shapes are punctuated by an elimination of facets in multiple stages instead of a single event. On the other hand, cohorts of irregular protoclasts exhibit no structural phase transitions as an aggregate, with individual eliminations of primordial surfaces being blurred by structural disorder.

We find direct measures of deviation from a spherical profile (i.e. the sphericity $\phi_{3}$ and ratio of minimum and maximum center of mass distances $r^{\mathrm{min}}_{\mathrm{max}}$) decrease monotonically with accumulating sustained slices. However, the oblateness measure actually rises at the initial stages of the erosion trajectory, reaching a peak when half of the original volume has been chipped away and then decreasing. This non-monotonic behavior of the oblateness is attributed to initially blade-like protoclasts which lose material more rapidly from the ends than from the edges, briefly becoming more oblate in shape before eventually being rounded into a spherical profile.

That data from distinct erosion scenarios (i.e. relative area and fixed velocity schemes) collapse onto common curves when plotted with respect to the remaining volume fraction confirms the universal dependence of observables on the latter Domokos1 . This phenomenon provides an avenue for using results in the relative area case to calculate time dependent observables for arbitrary erosion schemes, a technique we validate by reproducing results gleaned in the context of the fixed velocity scenario by direct Monte Carlo calculation. In addition, one also is afforded an efficient approach for calculating characteristic times (e.g. for specific remaining volume fractions) for a range of erosion scenarios. We have obtained closed form approximate expressions which bound these time scales from above, and which are in good quantitative agreement with the latter.