## Abstract

We have measured grain size distributions of the results of laboratory decompression explosions of volcanic rock. The resulting distributions can be approximately represented by gamma distributions of weight per cent as a function of *d* is the grain size in millimetres measured by sieving, with a superimposed long tail associated with the production of fines. We provide a description of the observations based on sequential fragmentation theory, which we develop for the particular case of ‘self-similar’ fragmentation kernels, and we show that the corresponding evolution equation for the distribution can be explicitly solved, yielding the long-time lognormal distribution associated with Kolmogorov's fragmentation theory. Particular features of the experimental data, notably time evolution, advection, truncation and fines production, are described and predicted within the constraints of a generalized, ‘reductive’ fragmentation model, and it is shown that the gamma distribution of coarse particles is a natural consequence of an assumed uniform fragmentation kernel. We further show that an explicit model for fines production during fracturing can lead to a second gamma distribution, and that the sum of the two provides a good fit to the observed data.

## 1. Introduction

Fragmentation is a ubiquitous phenomenon in many natural and engineering systems. It is the process by which an initially competent medium, whether solid or liquid, is broken up into a population of constituents. Examples occur in collisions and impacts of asteroids/meteorites, explosion-driven fragmentation of munitions on a battlefield, as well as of magma in a volcanic conduit causing explosive volcanic eruptions, rock explosions in mining, break-up of liquid drops, and many industrial processes, such as atomization and spraying [1–4]. Besides the mechanism of fragmentation the resulting frequency–size distribution of the generated constituents is of central interest. Initially, their distributions were fitted empirically using lognormal, Rosin–Rammler and Weibull distributions ([5] and references therein). The sequential fragmentation theory [2,5–7] and the application of fractal theory to fragmentation products [1,8–10] attempt to circumvent this empiricism by providing a more physical basis for the applied distribution. Both rely on an at least partially scale-invariant and thus self-similar random fragmentation process.

Our present interest stems from the mechanism by which explosive volcanic eruptions occur, and in particular from experiments in which vesicular volcanic rock samples explode naturally when suddenly depressurized from initial pressures of tens of MPa to ambient conditions [11,12]. The physics governing this fragmentation process has been successfully modelled and the observed fragmentation pattern can be numerically reproduced [13]. In these experiments, the fragmentation of natural rocks leads to grain size distributions (GSDs) which vary depending on the experimental starting conditions; typical examples of these are shown in figures 1 and 2. While our earlier theory provided an explanation for the fragmentation patterns observed, it did not consider the GSD.

Our purpose in this paper is to provide a stochastic model for the evolution of the grain size distribution during the explosion process, and thus to provide a theoretical description of the different distributions which are shown in figures 1 and 2. Our model combines a sequential model of the type outlined by Turcotte [1], but generalized to cater for the explosive process appropriate here, in particular by including in the description of the fracturing events in which the rock fragments a recipe for the production of fines, as observed in the experiments.

## 2. Experimental products and theory of explosive magma fragmentation

Constraining dynamics and products of explosive fragmentation of vesicular magma and rocks is key to better understanding volcanic eruptions and their associated products. Laboratory experiments have been conducted for almost 20 years with increasing vigour, investigating various aspects of fragmentation behaviour of magma analogues (e.g. [14–17]) as well as natural rocks and magmas (e.g. [10–12,18–26]). Experiments on rocks and magmas are conducted in a shock-tube apparatus where samples are exposed to rapid decompression from various starting conditions in the form of applied pressures ranging from 2 to 40 MPa and temperatures between 20°*C* and 900°*C*. This study uses experiments at ambient temperature only, however, very similar GSDs are also found for high-temperature experiments [19,20] owing to the brittle nature of these very fast processes.

Starting conditions are reached by slow gradual build-up of temperature and gas pressure using either argon or nitrogen gas; a dwell period prior to rapid decompression assures equilibrium conditions and negligible path effects. Rapid decompression is triggered by the bursting of a series of diaphragms, separating the high-pressure section with the sample from a large low pressure section, into which the generated clasts will be ejected. The low-pressure section is tailored to an efficient collection of the clasts, reducing the loss of very fine clasts to a minimum.

In this study, we analyse the products obtained from rapid decompression experiments of a BKB, a volcanic rock from Unzen Volcano, Japan. Cylindrical samples (diameter: 24 mm, length: 50–60 mm) were drilled out of a larger block. The connected porosity of each sample was determined from its geometrical volume and the matrix volume obtained by helium pycnometry. All samples exhibit very similar values of connected porosity between 45% and 49%. Isolated porosity was always estimated to be in the order of a few per cent only; thus, a significant role of isolated porosity on the evolution of GSD can be ruled out. The experiments were conducted at ambient temperature in a transparent autoclave and monitored with a high-speed video camera. Applied pressure was varied between 2.5 and 10 MPa, higher pressures were inhibited by the autoclave's strength. Figure 1 shows the individual grain size distributions of the clasts generated at increasing applied pressures. Size distributions were obtained by manual dry sieving at half-*ϕ* steps ([27], pp. 476–478 and references therein), where
*d* is the particle diameter and *d*_{1}=1 mm; thus, the equation above is often written as *ϕ*=−5 (equivalent to 32 mm) is assumed to be the coarsest fraction. The finest grain fraction used in this study corresponds to *ϕ*=4 (equivalent to 63 μm). Figure 2 comprises all four datasets as differential weight fraction of particles (top) or as cumulative weight fraction of particles (bottom), both plotted against *ϕ*. In both, the change of the distribution in shape as well as in position with increasing applied pressure can be seen. Further GSDs derived from rapid decompression experiments on volcanic rocks can be found for instance in [19] and [20]. These distributions are very similar to the BKB ones in shape, as well as in their responses to changes in the applied pressure. Such differences as do occur are mainly attributed to different samples being used, i.e. different connected porosity and permeability. Thus, we consider the GSDs presented in figures 1 and 2 to be representative for rapid decompression-induced fragmentation of vesicular rocks at increasing applied pressure. The advantage of the dataset presented here lies both in its cleanness, and in the additional monitoring of the fragmentation process which was obtained using a high-speed camera. This allowed further analysis of the individual fragmentation events which occurred in each sample.

The experimentally generated GSDs do not reflect any transport-related sorting, and in that they differ from the GSD of natural volcanic deposits. Their closest natural counterparts are the total grain size distributions (TGSD), which are harder to obtain, but more relevant to eruption dynamics [28,29]. We neglect transport-related grain size changes due to abrasion, collisions (particle–particle and particle–tank) and comminution, as these secondary effects scale with the fragmentation depth [30], i.e. the length of the experimental conduit, which is rather short in our laboratory setting (approx. 25 cm). However, these secondary grain size reductions also depend on (i) initial size distribution and (ii) fragmentation dynamics, as the ejected clasts and surrounding gas do not form a continuous granular flow with homogeneous density distribution, but a heterogeneous granular flow with variable density. The shape of the grain size distribution depends on the fragmentation process; for the dataset analysed here it depends in particular on the applied pressure: the higher the pressure, the more fragmentation events will occur before the pressure is reduced. In case of an applied pressure just slightly above the pressure needed to fragment the sample, only a few fragmentation events will occur, and a longer time lies between these events.

## 3. Sequential fragmentation theory

The basic mechanism of the theory of fragmentation was described by Kolmogorov [31]. He explained how a lognormal distribution of grain sizes could arise through the successive fragmentation of a rock sample. If the splitting position occurs at random, then the successive segments have diameters *d*_{i} given by *d*_{i+1}=*f*_{i}*d*_{i}, where *f*_{i} is drawn from a random distribution on (0,1). Thus,
*n*, *d*_{n} is lognormally distributed.

Brown [2] and Wohletz *et al*. [6] introduce a theory called sequential fragmentation theory, which they use to explain observed distributions in various rock fragment populations. Their theory is discussed further below. Its context is that of a stochastic process, in which the continuous-time evolution of the grain size distribution *n*(*x*,*t*), where *n* is number density and *x* is fragment volume (specifically, *n*(*x*,*t*) *dx* is the number of fragments having volume in the range (*x*,*x*+*dx*)), can be written in terms of a master equation which takes the form [32,33]
*p*(*x*) is the probability (per unit time) of fragmenting rocks of size *x* at time *t*, and *K*(*x*,*y*) is the fragmentation distribution of particles of size *x* from fragmentation of a rock of size *y*. Note that the time scale of the problem is just *p*^{−1}. If *x*_{0} is the volume of the initial block, then
*x*<*y*, so that

Much of the literature on fragmentation revolves around measurements of distributions having a power-law distribution in part of the range, and this is commonly thought to imply a self-similar process [34]. Two examples which we will focus on are Turcotte's own, where a cubic block is fractured into eight sub-blocks of equal size [1], and Brown & Wohletz's [5] power-law kernel, designed to produce a commonly assumed Weibull distribution, and popular among volcanologists (e.g. [7]).

Turcotte's model defines a kernel *K* in terms of a delta function (essentially a point source),
*λ*<1 is the volume fraction of the blocks which result from the fracturing event; in Turcotte's model *K*=*Cx*^{ν−1} (in [5], *ν*=1+*γ*∈(0,1)). Brown and Wohletz take *C* to be constant, but this is inconsistent with the requirement in (3.5), and the correct form for *K* is

### (a) Self-similar fragmentation

The purpose of this section is to show that for *self-similar fragmentation*, the integral in (3.2) takes a particularly simple form, of convolution type. We define a general self-similar kernel to be one of the form

We can also define a generalized Turcotte model (GT) which splits a block up into a number of different-sized fragments, of fractions *λ*_{1},*λ*_{2},… of the original block size. The resulting GT kernel is

It is convenient to relate particle size *x* to the sieve size variable *ϕ* introduced in (2.1); the two are related by
*d*_{1}=1 mm as before. The value of *A* is
*ϕ*. Additionally, we relate the number density *n* to the volume or weight fraction density *W* measured in terms of *ϕ*, as shown in figure 2. Details of the calculations are given in appendix A. There we show that *W* satisfies the simplified equation
*p*=1 (note that this defines the time scale of the process). The kernel *G* is determined by *g*, and *G*(*s*)=0 for *s*<0, and additionally
*W* is
*ϕ*. Since *W*=0 for *ϕ*<*ϕ*_{0}, the upper limit of the integral in (3.13) can equally be taken as *ϕ*−*ϕ*_{0} (and the lower limit taken as 0, since *G*=0 when *ψ*<0).

As we mentioned earlier, Brown & Wohletz [5] use their choice of kernel to motivate the derivation of the Weibull distribution. We also mentioned that their choice of kernel is inaccurate. Further, they derive the Weibull distribution by assuming the solution of (3.13) is a steady-state one. This makes no sense. And in any case, no equilibrium solution of (3.13) exists unless *G*(*ϕ*)=*δ*(*ϕ*), in which case any distribution *W* is a solution.

The solution of (3.13) is easily found using the Fourier transform
*w* is

Some useful insight can be obtained by solving the basic Turcotte model, for which we may take
*a* equal sub-blocks; Turcotte took *a*=8, but any integer will suffice, and in fact mathematically any positive number greater than one will do. Equation (3.13) then reduces to
*ϕ*_{0}+*α*, *ϕ*_{0}+2*α*, etc., and is determined by writing
*λ*_{1}=*λ*, *λ*_{2}=*λ*^{2}, etc. For example, this choice truncated at *λ*_{2} yields a second order such equation, providing *λ*=0.618 is the golden ratio. These problems are easily solved using generating functions, but useful closed-form solutions equivalent to (3.26) are not easily come by.

The general solution *w* whose transform is given by (3.20) is
*k*. Noting that *ϕ*) at speed *G*_{1}. This is the manifestation of the central limit theorem, and since *ϕ* is the logarithm of particle size, it gives a lognormal distribution in particle size.

Figure 3 shows the evolution of *w* with time when *G* given by (A 3) corresponds to a choice of the fragmentation kernel *g* as a uniform function, *g*=2, equivalent to the choice *ν*=1 in the Brown kernel (3.9). It can be seen that the solution approaches the Gaussian given by (3.31). A feature of note is that the value of *w* at *ϕ*=0 decreases exponentially with time, rather than descending instantly to zero. As we discuss in the following section, this is somewhat unrealistic in our experimental context. A further feature of note is that at large time the distribution decays faster than exponentially, unlike the fractal distribution (see appendix A) given by (A 7) (but note that for fixed time and large *ϕ*, the decay is weaker).

## 4. Application to rock fragmentation experiments

In considering the application of fragmentation theory to the experimental data, there are four issues to deal with, and we consider each of these in turn. These four issues are those of time evolution of the grain size distribution, the apparent self-similarity of the coarse fragment distribution, the apparent advection and truncation at coarse wavelengths of the distributions in figure 2, and the problem of the production of fines. Further, in order to account for some of these issues, we introduce what we call a reductive fragmentation theory.

### (a) Time evolution

One of the most striking features of the four graphs in figure 2*a* is the sense that they can be interpreted as a time evolution of a distribution, much as figure 3 shows a single evolving distribution; the distributions in figure 2 decay and move to the right as the initial confining pressure increases.

In order to understand why this should be so, we need to recall the mechanism of fragmentation, as discussed by Fowler *et al.* [13]. When a block fractures and separates from the column, its total stress *σ*_{tot}=(1−*ϕ*)*σ*−*ϕp* rapidly relaxes to equal the external pressure in the chamber, *p*_{c}, which itself is rapidly decreasing with time. Here, *ϕ* is the sample porosity, *σ* is the stress in the solid and *p* is the pore pressure in the gas. It follows that the effective stress in the block is
*σ*+*p*>*σ*_{Y}, this criterion in the block is thus
*t*_{c}≈1 ms is much smaller than the diffusive time scale *t*_{0}≈22 ms for the gas pore pressure to relax due to exhalation from the block, and therefore, at least for large blocks, the time to the next fracture is essentially controlled by the time from the preceding fracture for the chamber pressure to decrease by the quantity (1−*ϕ*)*σ*_{Y}. Therefore, the higher the initial confining pressure, the longer does the fragmentation process continue, and thus the final GSD represents the evolution of the distribution to a longer time. It is because of this that experiments with different initial confining pressures yield results which can be interpreted as representing distribution evolution over different times. In the simplest interpretation, the initial confining overpressure would simply be proportional to the number of fracture events, but in reality the correspondence will not be linear, because the diffusive gas relaxation time decreases with the square of the block size, so that at some small scale, block fragmentation will cease. In summary, we interpret the GSDs from experiments at different initial confining pressures as representing proxy time sequences of a single evolving GSD, much in the same way as proxy time sequences of, for example, evolving plant root systems are studied from a population whose members are sacrificed at different times.

### (b) Self-similarity

The self-similar nature of the long-term solution (3.31) suggests that we might seek a similar collapse of the experimental data on to a single curve, and figure 4 shows that such a collapse is possible, although the curves start to diverge at larger *Φ*; the reasons for this are discussed in §4e. Equally, it is clear that the chosen fitting gamma distribution decays too rapidly at larger *Φ*, but while one could choose other curves to provide a better fit to the heavy tail, the later discussion indicates that there is no purpose to this.

### (c) Advection and truncation

Given that our fragmentation model has an explicit solution in (3.20) or (3.27), and that we can interpret the results of different initial confining pressures as representing different times in the evolution of the distribution, it is tempting to seek to fit the data in figure 2 by the choice of a suitable kernel *G*. Of course, the problem with this is that we have no real idea what to choose for *G*. We might try the various generalized Turcotte kernels, or the Brown kernels, but initial efforts in this direction all yield figures similar to figure 3. A more enterprising possibility is to use the data to determine an approximate formula for *w*, and then use the inverse of (3.20), thus
*G*. Efforts to do this fail, apparently because the inverse problem is ill-posed: (3.20) only has a solution for *G* if *L*^{1} function *G*, as the Riemann–Lebesgue lemma implies that *f*, the inverse transform of *f*(*ϕ*)=*e*^{−αϕ}/*ϕ*, *ϕ*>0, *f*(*ϕ*)=0, *ϕ*<0, which is evidently nonsense.

The problem is that while there are some similarities in figure 3 to figures 1 and 2, there are some essential differences. Two which are critical are *truncation* and *advection*. Since the initial value of *ϕ*_{0}≈−5, we see that in all samples, the value of *W* at *ϕ*=*ϕ*_{0} goes to zero immediately; we call this truncation. Further, the first value of *ϕ* where *W*=0 advances with time; we call this advection. Neither of these properties is satisfied by the solution of the master equation (3.13).

### (d) Reductive fragmentation theory

The physical reason for this is in fact obvious. After the initial fracture of the sample, it is split into two pieces, each of them smaller than the original. So long as the blocks continue to fracture, both truncation and advection will occur. In order to describe this, we revisit the derivation of the master equation. The equation (3.2) is derived from an incremental difference equation
*p*(*x*)Δ*t* is the probability that a block of size *x* will fracture in a time Δ*t*. In our situation, if Δ*t* is the time interval between fractures, then we need to take *pΔt*=1, and this leads to the difference equation
*reductive* fragmentation theory, as the block sizes reduce at each step.

We proceed as in §3a to define self-similar fragmentation by means of the self-similar kernel *K*(*x*,*y*)=(1/*y*)*g*(*x*/*y*), and this leads to the reductive equivalent of (3.13),
*G* is as defined earlier.

Truncation and advection are determined by choosing a kernel *g*(*η*) such that *g*=0 for *η*>1−Δ*η*, where we imagine Δ*η* is small (though this is not necessary). This essentially says that any block fracture yields a following block of maximum relative size 1−Δ*η*, with the details of the grain sizes produced dependent on the kernel *g*.

We show in appendix A that the solution of (4.9) can then be written as
*χ*∝Δ*η*, approximately, and *U*_{n}(*ξ*)≡*U*(*ξ*,*nΔt*) satisfies the difference equation
*G**(*ψ*) is a phase shifted version of *G*, but satisfies the same integral constraint
*ψ*<0. Equation (4.10) demonstrates the fact that the solution for *W* is advected downstream in *ϕ*, at speed Δ*χ*/Δ*t*.

The solution for *U* can be written in terms of the Fourier transform. Explicitly,
*g* is constant, so that *G* and therefore *G**∝*e*^{−2Aχ}. Then *g*∝*η*^{λ}, *λ*>−1, indicating a weighting towards production of fines (*λ*<0) or coarse particles (*λ*>0), and this also leads to a gamma function (4.14), but with 2*A* replaced by *A*(2+*λ*). Note that the exponent *λ* is unrelated to its previous use in the Turcotte kernels.

### (e) The production of fines

It is clear in figures 1 and 2 that the distributions are heavy-tailed at small particle size (large *ϕ*). Both figures show that fine particles with *ϕ* is fairly mild. We associate this with the production of fine particles during the fracturing process. Figure 5 shows typical images of the vesicularity of the rock samples. We imagine that as a fracture passes across the sample, a range of particles of sizes

In order to describe this, we return to the reductive equation (4.8), but now we consider the particle size distribution to consist of a coarse population *n*_{c} and a fine population *n*_{f}, with the dividing size being *x*_{f}. Our assumption is that the fracture of coarse particles produces both coarse and fine particles, but only the coarse particles fracture, essentially because the fine particles are at the scale of the vesicularity, where the fracturing mechanism (exhalation of pore gas) does not apply. We thus take *K*=0 both for *x*>*y* and also for *y*<*x*_{f}. In consequence we have
*x*_{f} is sufficiently small that it does not affect the coarse fraction. Insofar as we suggest a typical coarse fraction gamma distribution which tends rapidly to zero, this seems reasonable.

The fines fraction is then given by
*K* when *x*<*x*_{f}, which represents the production of fine *x*<*x*_{f} from coarse *y*>*x*_{f}. In keeping with our earlier discussion, we wish to retain the power law of the Brown kernel, but this needs to be modified so that *K*→0 as *x*→*x*_{f}−. The simplest choice with this behaviour is the modified Brown kernel
*ν*>−1, and we have introduced the logarithmic factor so that *K*_{f}=0 at *x*=*x*_{f}, assuming as we do that *μ*>0. As for the Brown kernel, the exponent *ν* as well as *μ* will be assessed by comparison with the experimental data, but the pre-factor is determined by a mass conservation argument. This is detailed in appendix A, where it is shown that we can write (4.17) in the form
*B** is an *O*(1) number which depends on the precise amount of fines produced in the fracturing process.

Note that this fines production kernel is *not* self-similar. Putting (4.18) into (4.16), we obtain
*W*(*ϕ*,*t*), we have, with an obvious notation,
*ϕ*>*ϕ*_{f}. If we assume that *ϕ*_{f}−*ϕ*_{0} is large, and that *W*_{n} is given by the gamma distribution in (4.14), then (4.20) simplifies at time step *n*>0 to
*p*_{0} (thus *t* and thus *n*), and that they decay weakly with *ϕ*. Indeed, comparison with (A 7) indicates that *D* is fractal dimension, as described in appendix A. Typical observed values of *D* are in excess of 1.8 [35,34], p. 44, so that *ν*<−0.6, which would suggest that (*ν*+1)*A*<0.8, which seems consistent with the slow decay in figure 2. In principle, the value of *ν* should depend on the distribution of vesicle spacing, as shown in figure 5. Hartmann [35] notes that the lower range of *D* is associated with ‘simple fragmentation’ with little abrasion, as may be appropriate here. It should be pointed out his value for ‘ash and pumice’ is *D*=3.54, and commonly *D*>3 for such deposits, but he associates such high values with extreme grinding and abrasion, which is less applicable here. It is an interesting question how such exponents are obtained, but not one of relevance here. Note that (A 7) requires that *D*<3 for sufficiently large *ϕ*, so that values *D*>3 must represent the rising limb of the distribution; it is natural to suppose that the falling limb is lost as the airborne ash cloud.

In summary, we have used assumed power-law type kernels to produce predicted forms of GSD for experiments with different initial confining pressures. These pressures serve as proxy measurements of fragmentation time, as explained in §4a. We have shown qualitatively that the theory has the possibility of explaining the data, but we have as yet no quantitative comparison. This will be done in the following section, where we will choose values for the wealth of parameters which are consistent with the data.

## 5. Discussion

We define the gamma distribution function as
*x*>0, and *x*<0. Note that its integral is one. The fitting curve in figure 4 is thus *W*(*ϕ*,*p*) (measured as weight per cent) are thus of the form
*p*, *α*, *ϕ**, *a** and *n* corresponding to the fitting approximations in figure 4 are given in table 1. Noting that

We now wish to relate these fitted curves to the predictions given by (4.10) and (4.14) for the coarse fraction, and to then adjust the prediction by including the fine fraction prediction (4.21). Equations (4.10) and (4.14) suggest a coarse approximation
*n* representing time is a proxy for initial confining pressure *p*, and in (4.14) we would have *αa**=2*A*, or for the more general power law *g*∝*η*^{λ} (see after (4.14)), *αa**=*A*(2+*λ*), thus *λ*<0 would indicate a weighting towards fines production if *A** should be 50, since the sum of the weight fractions *ϕ*=0.5. Since the integral of *g* is one, this implies *A**=50. Bearing in mind (5.3), (5.5) is equivalent for any *α* to

We now add the fines production formula (4.21). This can be written in the form
*B*<1 is a number associated with the fracturing process, defined in appendix A before (A 12).

We add the fines fraction to the coarse fraction, subtracting a corresponding part of the coarse distribution to conserve mass. Writing *b*=*C***β*^{n}, this suggests a combined approximation

We have used this formula to fit the data, and the results are shown in figure 7. In fitting the data, we used initially the values in table 1 with some adjustment to fit the coarse part, and then the parameters of the fine part were adjusted (principally the decay parameter *γ* and the amplitude parameter *b*) to fit the heavy tail. The fits are obviously very good, though we have not endeavoured either to find the best fit or to choose parameters which best suit our suggestion that the data represent increasing values of *n*, although it can be seen from table 2 (columns for *b*, *n*Δ*χ* and *n*) that this is largely the case. We think it more significant that the analytic form of the prediction can be made to fit the data relatively easily; as anyone will know who has tried, this is not always the case.

One comment on the parameters is the value of *A**, which ought to be 50. This is approximately the case, although the value for *p*=10 is a bit high. This is due to the fact that while our fitted curve has a component with non-zero area in *ϕ*>4, the data are truncated by the finest mesh size at *ϕ*=4; thus, the final bin of the data includes all the fines with *ϕ*>4, although the fact that the final bin has a quantity commensurate with the trend of the other bins actually suggests that there is a cut-off size in the fine particle production; this might explain why the fitted curves do so well up to *ϕ*=4.

A second comment concerns the choice of *γ* given in (5.9). As indicated in (A 7), the corresponding fractal dimension is
*γ*<*A*≈2.08. This it is, and the fitted range 0.4–0.8 corresponds to a value of *D*=1.85–2.42. These are well within the range quoted by Turcotte [1], table 1.

## 6. Conclusion

Grain size distributions from experiments on the fragmentation of volcanic rocks have a characteristic shape which depends on the initial confining pressure, and which can be well represented by the sum of two gamma distributions. We suggest that the variation with confining pressure represents an equivalent evolution of the distribution with time following pressure release, based on our earlier theory [13], and we then pose a stochastic process theory to describe the evolution of the distribution with time. This theory introduces the concept of self-similar fragmentation kernels, and leads to a predicted gamma distribution for coarse fragments assuming a general power-law kernel. We further associate the heavy tail of the distributions with a secondary process, whereby the primary fracturing produces fine particles. A generalized power law for this (non-similar) process predicts a second gamma distribution whose amplitude grows geometrically in time. We then show that the combined distributions can easily fit the observations, and that the fitting parameters are broadly consistent with the hypothesis that the data represent an effective time evolution of the grain size distribution.

## Data accessibility

The data plotted in figures 1 and 2 are contained in the electronic supplementary material.

## Authors' contributions

B.S. performed the experiments and found the grain size distribution, and produced the experimental figures. A.C.F. devised and solved the model, and wrote the paper. Both authors worked on the final submitted text.

## Competing interests

We have no competing interests.

## Funding

A.C.F. acknowledges the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland grant 12/1A/1683. B.S. acknowledges the support of the European Commission (FP7-MC-ITN, grant no. 289976: NEMOH).

## Appendix A

We relate the number density *n* to the volume or weight fraction density *W* measured in terms of *ϕ*, as discussed in §3a. *nx* *dx* is the volume of fragments in the range (*x*,*x*+*dx*), and so we have
*ϕ* increases with decreasing *x*), and it follows that
*G*=0 for *χ*<0. Then the master equation (3.2) takes the form
*P*=1.

**(a) Fractal dimension**

Another commonly measured quantity is the number of fragments of radius larger than *R*, and we will denote this as *M*(*R*). We have
*W*, this is
*D* of a power-law distribution is given by the relation *M*∼*R*^{−D} for small *R* [34], and in terms of *W*, this implies
*ϕ*. Convergence at large *ϕ* requires *D*<3, as is commonly found [34].

**(b) Reductive fragmentation theory**

We consider the equation (4.9) for *W*
*χ*≈Δ*η*/*A*. We define
*G**(*ψ*) satisfies the integral constraint (4.12), and is zero for *ψ*<0. The solution of (A 11) then proceeds as in the main text.

**(c) Fines production**

We consider the value of the pre-factor *C* in (4.17). The volume of fines produced from a fracture of a block of size *y* is *d* is *B*<1, equating the two implies

- Received December 11, 2015.
- Accepted May 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.