The Glencoe caldera is a well-studied example of a caldera system exposed to intermediate depths along the glacially excavated glen. We present a first quantitative assessment of clast-size population and matrix chemistry from the flinty crush rock that occurs on the main ring faults. Size–shape metrics of clasts differ from those of a ‘normal’ pseudotachylite from the Outer Hebrides. Both samples display good power-law clast-size populations, once allowance is made for dissolution of a portion of clasts into the melt that contained them, with fractal dimensions of 2.7 and 4.0 for the Outer Hebrides and Glencoe samples respectively. Mass-balance calculations of flinty crush rock matrix chemistry imply an origin by mixing between rhyolite and pseudotachylite that was derived from semipelitic host rock. This means that the flinty crush rock was transported some distance from the point of frictional heating, as previously proposed, as semipelitic rocks are not present at the surface at Stob Mhic Mhartuin but are likely to be present at depth. This transport, and mixing with rhyolite magma, would have provided the time and thermal energy for clast dissolution beyond that possible during normal pseudotachylite generation and quenching on fault wall rocks.

Supplementary material: Results of simulations, major element compositions of the matrix in the flinty crush rock sample and a supplementary figure are available at

The study of pseudotachylites has developed substantially since the term was introduced by Shand (1916), with clear evidence from experiments and natural samples that frictional melting can produce pseudotachylite during seismic rupture. Evidence for partial melting includes the following: (1) high temperatures in rupture models (Cardwell et al. 1978; Lachenbruch 1986) and high-velocity friction experiments, combined with direct generation of frictional melt in the experiments (Ohtomo and Shimamoto 1994; Di Toro et al. 2010); (2) quench textures indicative of crystallization from a melt (Di Toro and Pennacchioni 2004); (3) high-temperature polymorphs of plagioclase (Nestola et al. 2010); (4) preferential preservation of refractory minerals as clasts (Maddock 1983; Spray 1992; Lin and Shimamoto 1998; Jiang et al. 2015); (5) presence of diffusion gradients around clasts, consistent with dissolution into a melt (Dobson et al. 2018); (6) chemistry of glassy matrix, which requires peritectic melting of micas (Dobson et al. 2021). In addition to well-documented tectonic and impact-related examples, pseudotachylites have been described from large landslides (e.g. Masch et al. 1985; Legros et al. 2000, and references therein) and from a range of volcanic settings. Volcanic occurrences include slip surfaces around volcanic plugs in acidic arc settings (Tuffen and Dingwell 2005; Kendrick et al. 2014) and bounding faults of caldera collapse events (Han et al. 2019). Indeed, the first suggested example of a frictional melting origin in the geological literature was for the ‘flinty-crush rock’ occurrence on the bounding faults of the Glencoe caldera (Clough et al. 1909; Bailey 1960). The origin of the Glencoe flinty crush rock is, however, a matter of some debate, with researchers attributing it to shear fluidization, extreme explosion breccia or gas fluidization (e.g. Reynolds 1954; Hardie 1963; Roberts 1966) or citing shape and size characteristics of clasts within the flinty crush rock (and surrounding highly damaged quartzite host rock) as evidence of an origin by extreme comminution and melting (Kokelaar 2007) or gas attrition (Roberts 1966); both Kokelaar and Roberts still describe it as a pseudotachylite. Despite these divergent claims about the implications of size and shape variations of clasts in the debate on the origins of the Glencoe flinty crush rock a systematic analysis of their metrics is yet to be published. Studies of the Glencoe caldera have been seminal in developing a conceptual understanding of caldera formation processes, and a more complete understanding of the flinty crush rock that occurs along several sections of the main caldera ring faults might help constrain the dynamics of caldera eruptions.

The unequivocal identification of natural pseudotachylites can be difficult, and earthquake pseudotachylites are much rarer than the energetics of seismic rupture would suggest (e.g. Spray 1995). This might be due to subsequent modification of glassy matrix, which can destroy textural evidence of co-seismic melting (Kirkpatrick and Rowe 2013; Fondriest et al. 2020). There has therefore been recent interest in developing diagnostic metrics for pseudotachylites, based on two lines of investigation: size–shape metrics of survivor clasts and chemistry of the matrix material, thought to represent the melted portion of the fault material, between the clasts.

Dobson et al. (2021) compiled a large dataset of published pseudotachylite matrix compositions and showed that in many cases it is not possible to generate matrix compositions from linear mixtures of the minerals in the wall rock of the fault. In this case, rapid non-equilibrium melting of low-melting mica minerals requires crystallization of peritectic potassium feldspar and other peritectic phases (most commonly Al-rich phases such as corundum, mullite or an Al2SiO5 polymorph). The Al–(Fe + Mg)–K ternary plot was found to be diagnostic for pseudotachylites derived from igneous protoliths.

Clast-size metrics and clast-shape metrics are predicated on the fact that deformation of rocks causes grain-size reduction and modifies grain shapes in ways that are distinct from other geological processes such as grain growth or melting. In the brittle field, deformational energy is dissipated through friction and through fracturing grains. Theoretical work shows that fracturing of grains produces self-similar distributions of grain size (here we use the effective radius; r=(A/π)1/2, where A is the measured area of a grain) and number of grains of a given size or larger, n, following a power-law distribution:
with a fractal dimension, D, around 2.6 (Sammis et al. 1987; Marone and Scholz 1989); J is a fitting constant associated with the size of the dataset. On the basis of these pioneering studies, it is common to plot clast-size–number distributions of natural fault gouges on log–log graphs with linear regions (the fractal range) in these plots (often towards the larger clast sizes) taken to be representative of the fractal dimension of the system. Natural gouges can have fractal dimensions substantially above 2.6, which is commonly attributed to processes preferentially acting on smaller grains, such as abrasion and extreme comminution (Morgan et al. 1996), and has been taken to be an indicator of seismic slip (Balasmo and Storti 2011). For a discussion of fractal dimension in rocks deformed by brittle processes the reader is referred to Glazner and Mills (2012). Most studies of pseudotachylite clast grain size analyse 2D images of sectioned rocks, which modifies equation (1a) to
(see, for example, Sammis et al. (1987) or Glazner and Mills (2012) for a discussion of random 2D sampling of 3D grain populations).

Modified power-law behaviour of pseudotachylite survivor clasts

Survivor clasts in pseudotachylites often show significant (convex-up) curvature away from power-law behaviour, with a fall-off of clast numbers at small grain sizes. Although resolution issues can produce similar curvature at small grain sizes, multiple sampling at a range of resolutions can be used to test whether this curvature is a genuine population property in pseudotachylites. Curvature towards smaller numbers of small grains in pseudotachylites has been attributed to preferential melting of small clasts (e.g. Deb et al. 2015). The study of Montheil et al. (2020) on experimentally generated frictional gouges and pseudotachylites showed that the pseudotachylite survivor clasts had a high fractal dimension of 3.0, compared with 2.4 for the gouge where deformation was halted before melting, but a reduced range of grain sizes that could be reasonably fitted by a power-law distribution (the fitted range) in the pseudotachylite. Most previous studies have fitted power-law distributions to pseudotachylite survivor clast populations using an arbitrary choice of fitted range (for example, to have a coefficient of determination, R2, >0.99 in the study by Montheil et al. 2020) to exclude data at small grain sizes that show strong deviations from power-law behaviour. Here we use a simple modification of the power-law equation:
which models the effect of reducing every clast radius in a log-normal distribution by some constant amount, r*. J and K in equations (1a) and (2a) are constants, and D is the 3D fractal dimension of the clast population.

To demonstrate the effect of removing a constant radius from clasts that initially follow a power-law distribution we plot equation (2a) for a fractal dimension D = 3 and for three cases of removed radius, r*, of 0, 2.5 and 20 in Figure 1. Values of n are calculated for r values between 0.1 and 1000. The fall-off of n at small grain sizes is clear. The form of the curves for r*= 2.5 and r*= 20 is identical but shifted in log(r) space by log (2.5/20). Simple power-law fits to the three curves, using data with radii between 10 and 1000, return apparent fractal dimensions of 3.0 (R2 = 1), 2.9 (R2 = 0.9996) and 2.55 (R2 = 0.99) respectively for r*= 0, 2.5 and 20. Even relatively small reductions in diameter can significantly alter the apparent fractal dimension. In the case of D = 3, a grain-size reduction of 2% of the maximum grain size reduces the fitted D by 0.45 in the ‘well-fitted’ (R2 = 0.99) region; the discrepancy becomes larger for larger values of D (see the Supplementary material for further results of simulations).

Equation (2a) is equivalent to the modified power-law equation adopted by Nagahama et al. (1992) and subsequently by Tsutumi (1999) and Ray (1999, 2004):
But we chose the form in equation (2a), as it clarifies the physical meaning; that is, removal of a constant radius from every grain. Ray (1999) explained the left-hand tail-off in clast numbers, and equation (2b), in terms of isochemical melting of clasts as they are heated by conduction from the surrounding superheated melt. In this case, the melting front in a clast is controlled by diffusion of heat into the clast, which, being more efficient for clasts with smaller radii of curvature, should produce higher melting rates for smaller clasts (Bizzarri 2014). This produces a left-hand tail-off the removed radius, which should be should be larger for smaller clast sizes, not constant across all clasts, and hence a more strongly curved distribution. Isochemical melting might be appropriate for initial melting on a shearing fault plane where heating is so rapid that equilibrium between grains cannot be maintained but it is unlikely to persist once melt is lubricating the fault and heat production rates are not as high. The alternative to isochemical melting is dissolution of refractory survivor clasts into the surrounding melt as the system seeks to re-establish chemical equilibrium. Dobson et al. (2018) observed that Si diffusion profiles are preserved around quartz survivor clasts in response to dissolution into a melt that is undersaturated in silica with respect to quartz. They successfully used these diffusion profiles to model the thermal history of pseudotachylite veins on the assumption that the diffusion profiles only develop once fault motion has ceased. This implies that while the fault is in motion any excess silica added to the melt around dissolving clasts is advectively mixed into the flowing melt. Under this ‘well-mixed’ scenario the rate of dissolution is controlled by detachment kinetics, rather than diffusion, resulting in a rate of dissolution that is independent of clast size (e.g. Shaw 2004), and the equations (2a) and (2b) are valid.

The log-normal grain-size distribution is common in geological materials that have experienced significant grain growth (i.e. metamorphic and igneous rocks) and provides an alternative base-line grain-size distribution for comparison with those produced by equation (2a). In the absence of a differential stress grain size will increase, with the rate of grain growth being dependent on the temperature, dominant grain size and grain growth mechanism (Marchant and Gordon 1972; German 2010, and references therein). It produces log-normal grain-size population distributions that are also upward-concave in log-size, log-number space, but with curvature of the distribution at all grain sizes. It might therefore also be possible to distinguish pseudotachylite and phenocryst populations from their grain-size distribution.

Here we determine the size metrics, determined by a semi-automated procedure, of clasts from a purely fault-generated pseudotachylite from the Outer Hebrides and flinty crush rock from Stob Mhic Mhartuin in the Glencoe caldera, which might contain an additional magmatic contribution. Both samples show upward-convex deviations from simple power-law distributions that are well reproduced by the modified power-law function. The Outer Hebrides pseudotachylite sample is best fitted by the modified power-law function. The Glencoe flinty crush rock is equally well fitted by the modified power-law or a log-normal distribution, which might be indicative of very large amounts of dissolution of clasts. Major element chemical analyses of the matrix of Glencoe flinty crush rock are presented, which support it being a mixture between a pelite-derived pseudotachylite and rhyolite magma, consistent with the clast-size metrics for this sample.

Note on terminology

We use the term ‘clast’ here to denote any individual crystalline object that is embedded in glassy matrix that was, or appears to have been, originally molten. This definition encompasses both the products of brittle failure on shears (clasts sensu stricto) and crystallization from magma (often termed grains or -crysts). We reserve the term ‘grain’ for crystalline components of rocks (i.e. multi-grain aggregates) and rock fragments, which themselves might form clasts. No genetic origin is implied through the use of these terms unless specifically stated.

Outer Hebrides pseudotachylite

A low-angle fault zone, the Outer Hebrides Fault Zone, cuts Lewisian gneiss, cropping out along the eastern shore of the Outer Hebrides for their entire length (Fig. 2). The 1–6 km wide fault zone has a long and complex deformation history with active periods of deformation from the Neoproterozoic to the Carboniferous (Imber et al. 2001). Within the shear zone, pseudotachylites show complex relationships with ductile shear zones, both cutting and being overprinted by mylonites (Sibson 1975). To the west, structurally below the base of the Outer Hebrides Thrust, pseudotachylites cut largely undeformed Lewisian rocks, and their occurrence increases in frequency and size southwards through the Outer Hebrides, culminating in brecciated pseudoctachylite bodies that reach a metre or more in thickness in Barra. Individual pseudotachylites in the Outer Hebrides have been 40Ar/39Ar dated, yielding ages of 1900, 1200 and 700 Ma (Sherlock et al. 2009) or 430 Ma (Kelley et al. 1994), highlighting the long and complex deformation history recorded in these rocks.

Samples of pseudotachylite for the present study were collected from Rubha Áird a’ Mhuile [NF715297] on South Uist, approximately 10 km west of the Outer Hebrides Thrust zone. This locality contains numerous cobbles and boulders containing pseudotachylite as developments on fault planes, injection veins and breccias within gneiss. The source of beach boulders is likely to be the headland immediately west of the beach, where pseudotachylite fault-veins of significant thicknesses up to 10 cm occur (Fereé et al. 2016). The sample analysed here is from a pseudotachylite-hosted breccia that contained pseudotachylite volumes of several cubic centimetres between rounded breccia clasts of unfoliated granodioritic gneiss.

Glencoe flinty crush rock

The valley and surrounding mountains of Glencoe incise an andesitic–rhyolitic caldera volcano system, 420 Ma in age, exposing the upper magma chamber and deep sections of the volcanic system (Fig. 2). The bounding faults of the caldera system are exposed and show complex faulting and intrusion relationships. Displacement on the caldera faults is estimated at 0.5–1.7 km, with evidence of rapid, ‘superfault’, movement associated with eruptive phases and caldera subsidence (Kokelaar 2007). The surrounding country rocks are quartzites, pelites and psammites of late Precambrian age, with well-developed shear fabrics where the caldera faults cut them (Fig. 3). A dark, aphanitic rock reminiscent of pseudotachylite decorates some fault surfaces of the main bounding faults, forming flow-banded sheets and commonly showing mixing relationships with the phenocrystic rhyolite that was intruded along the fault system. This is the ‘flinty crush rock’ of early workers. It is particularly well developed around Stob Mhic Mhartuin [NN208574] where it forms two sets of parallel-sided sheets between 2 and 20 cm in thickness, forming the edges of the ring-fault intrusion. These flinty crush rock sheets are bounded on one side by country rock (quartzites and psammites of the Eilde Quartzite and Eilde Flags Formations) and on the other by porphyritic monzonites and rhyolites – the ring-fault intrusion. Where the flinty crush rock contacts the country rock around this locality it strongly resembles fault-generated pseudotachylites, including geometric relations reminiscent of injection veins, but the contact with the ring-fault intrusion shows strong textural evidence for magma mixing, leading Kokelaar (2007) to conclude that the pseudotachylite had been transported for some distance, along with rhyolite magma during an early phase of a caldera eruption. A full description of the flinty crush rock and its field relationships has been given by Kokelaar (2007) and in the references therein. A small sample of flinty crush rock was collected from the talus immediately below the main exposure of the inner (southern) wall of the fault intrusion on Stob Mhic Mhartuin, which showed typical flow banding and had a contact with the neighbouring quartzite (Fig. 3). The flinty crush rock portion of the sample was ∼2 cm thick but the original thickness of the layer it was derived from is uncertain as it was collected as a loose fragment. In addition, the images presented by Kokelaar (2007) were analysed for comparison with a literature example of Glencoe flinty crush rock.

Synthetic pseudotachylite

The granite-hosted experimentally generated pseudotachylite studied by Montheil et al (2020) was re-analysed here as an example of a pseudotachylite sample for which the process of pseudotachylite formation is well characterized, and which has not undergone subsequent modification (post freezing of the melt). This sample was produced in a high-velocity rotary shear experiment, at a shearing velocity of 0.08 m s−1 and a normal stress of 40 MPa. The clast-size data for ‘granite’ from figure 3 of Montheil et al. (2020) are refitted here by equations (2a) and (2b). (See Montheil et al. (2020) for further details of the sample and experiment.)

Image collection and analysis

Size and shape metrics of clasts in the Outer Hebrides and Glencoe samples were collected from photomicrographs of petrographic thin sections in plane- and cross-polarized light. Images were collected using a Jenoptik ProgRes digital camera mounted on a Leitz Orthoplan polarizing microscope. Clasts were distinguished from dark matrix material by their difference in optical absorption, which was sufficiently large to allow a semi-automatic fitting procedure using the ImageJ suite of programs (Schneider et al. 2012). Clasts in photomicrographs were highlighted by choosing a threshold intensity with pixels with higher intensity than the threshold belonging to clast material and pixels with lower intensity belonging to matrix material. The threshold was chosen manually for each image such that clast shapes were defined but very little, or no, matrix material contained pixels with intensity above the threshold. The image was then made binary based on this threshold and any holes in clasts (i.e. regions of pixels with below-threshold intensity) were filled in, using the Fill Holes algorithm. There is a trade-off in threshold intensity value between defining clasts such that they produced a single object at this stage and having touching clasts defined as single objects; this was dealt with in part by manually tracing around larger or darker clasts (such as poly-grain clasts containing a significant amphibole content) so that, after holes were filled, they produced a single object rather than a constellation of smaller objects. Touching clasts were separated using the watershed algorithm in ImageJ with any inappropriate splitting of larger clasts manually removed. Figures 3 and 4 show typical examples of photomicrographs and identified clasts for the two samples. The Analyse Particles algorithm was then used to measure the area, A, of all clasts excluding ones touching the edges of the image. This process was repeated for between 10 and 15 photomicrographs for each sample at objective lens magnifications of ×2.5 and ×10. At ×2.5 magnification the entire thin section was sampled by the grid of photomicrographs. Photomicrograph locations at ×10 magnification were chosen at random within areas of pseudotachylite to avoid any sampling bias, with the proviso that images were not overlapping. There was some variability in clast concentration between images, but the population distributions were consistent between each image; this was confirmed by visual comparison of the clast number–size population distributions as plotted on log–log plots (see next section). Because each photomicrograph appeared to be sampling from the same clast-size population is it is reasonable to combine data from each photomicrograph at a particular magnification into a single dataset. This has the possible effect that the largest clasts might be underrepresented, as they have a disproportionately large probability of touching the edge of an image and being disregarded (as seen in Figs 4 and 5); however, this will affect only the very largest clasts, which do not contribute significantly to the results of this study.

Particle size metrics

Data are reported here as a function of clast effective radius, defined as r=A/π, where A is the measured area of a given clast. This single metric is assumed to be representative of the size of clasts even though it is a mean radius of a 2D slice through a complex non-spherical 3D object. Data are ranked in order of decreasing size, with the rank, n, starting at unity for the largest value of r.

We have collected very large datasets, with total numbers of clasts, N, in many thousands for each sample, meaning that statistical inferences are likely to be robust. For comparison between different magnifications for the two natural samples the numbers of grains were normalized to the area analysed across each combined dataset: n10=n10(AA2.5/AA10), where AA2.5 and AA10 are the total area of photomicrographs analysed at ×2.5 or ×10 magnification, and n2.5=n2.5.

Three functional forms were fitted to the data: equations (1a) and (2a) represent a simple fractal distribution across all the data, and include a constant correction to the grain size respectively, and
represents the data as a cumulative a log-normal distribution, appropriate for the case of normal grain growth. The constants μ and σ in equation (3) are the mean and standard deviation of the log-normal distribution. The scalar constants J, K and L are applicable across data at all magnifications if the fits are applied to n*, but would otherwise vary between datasets at different magnifications. Each of the three equations provides values of n* as a continuous function of r, and discrete values of n* were calculated at each observed r for comparison with the data collected from ×10 magnification images. The ×10 data were used as they have resolution to smaller clast sizes than the ×2.5 data, but if a given fit is truly a good estimator of the clast-size population it should also fit the ×2.5 data to larger clast sizes than the ×10 data can resolve. The data were fitted using a generalized reduced gradient algorithm in Excel, minimizing the log-differences between observed and expected n* values: RSS = ∑{[ln(n*obs) – ln(n*exp)]2}. Log-differences were used as the underlying power-law function that we expect the data to follow is linear in log–log space and this procedure avoids biassing the fit towards the large n* values at small size.

Chemical analyses

The thin section used for optical microscopy was polished for electron microscopy, finishing with 0.3 μm alumina powder and coated with 30 nm of carbon. It was imaged using a Jeol JSM-648OLV Scanning Electron Microprobe operating at 15 kV accelerating voltage with a spot size of 50, corresponding to a beam current of approximately 10 nA. The source was a traditional tungsten filament giving a resolution for imaging of ∼0.1 μm. Characteristic elemental X-ray fluorescences were measured for compositional analysis using a 30 mm Xplore solid-state energy-dispersive X-ray detector running under Oxford Aztec software. Regions of matrix (without survivor clasts) of between 100 and 1000 μm2 were analysed by rastering the beam over the selected area to obtain representative average compositions. The analytical uncertainty is typically 0.2% for these analyses; however, there is always the possibility of unintentional inclusion of small survivor clasts into the analysed area, which would affect the composition to an unknown extent. Because the most common survivor clast is quartz the most likely component to be affected by this is SiO2, and elevated silica content was observed in three analyses.

Mass-balance calculations

Mass-balance calculations were performed, following Dobson et al. (2021), to assess the origin of the matrix material. Measured matrix compositions were fitted by a linear mixture of mineral phases that make up the source rocks. For a given mixing model there are six mineral phases, whose phase proportions can vary. These minerals, which comprise the rhyolite and monzonite melts found in the fault-intrusion at Stob Mhic Mhartuin are (from the monzonite and granite analyses of Garnham 1988) plagioclase, muscovite, biotite, ulvöspinel, quartz and orthoclase. Chemical compositions for each mineral phase were those determined by Garnham (1988), where given (for biotite and plagioclase feldspar) or were assumed to be pure end-member compositions (for quartz, muscovite, orthoclase and ulvöspinel). The eight chemical components SiO2, TiO2, Al2O3, ‘FeO’, CaO, MgO, Na2O and K2O were included in the calculation. These eight components were recalculated as atomic percentages and normalized to 100%. The model concentration of the ith chemical component, Xim, is expressed as a sum of the concentration of that component (Xi) over all phases weighted by the phase fraction, φj:
where Xij is the concentration of the ith chemical component in the jth phase. Under the constraint that phase fractions must sum to unity, q=6, so there are five independent variables in the system (five of the six phase fractions) and the model composition will sum to 100% for normalized chemical compositions. The model composition is fitted to the composition of the analysed matrix material by minimizing the sum, over all components, of squared differences between the analysed fault concentration, XiMatrix, and the model concentration, SSm=i=18(XiMatrixXim)2 and the residual sum of squares (RSS) is used as a measure of fit quality. A fit was considered to be good if RSS is <0.2.

A second set of (‘magma mixing’) mass-balance calculations were performed to study the possibility of mixing between pseudotachylite and fault intrusion as an origin for the flinty crush rock. Equation (4) is used, with q=3 ‘phases’ corresponding to the bulk compositions of pseudotachylite, ring-fault intrusion magma and quartz. In this case three possible pseudotachylite compositions were considered, corresponding to average values of pseudotachylites hosted by pelites, semipelites or psammites, from the compilation of pseudotachylite matrix compositions of Dobson et al. (2021). The fault intrusion composition used was the porphyritic granite of Garnham (1988), which corresponds to the rhyolite that is in intimate contact with flinty crush rock at Stob Mhic Mhartuin (Kokelaar 2007). In practice, all best-fitting models needed only the semipelite pseudotachylite composition, but also extra silica, corresponding to small quartz clasts, which could not be avoided in the matrix analyses. For these models, therefore, q=3, corresponding to the bulk compositions of semipelitic pseudotachylite and porphyritic granite, and quartz; there are two independent variables. A fit was considered to be good in this case if RSS is <1.

The locality where flinty crush rock was collected, along with a thin section photograph and photomicrographs of the sample analysed here are presented in Figure 3. Flow banding and mixing between darker and lighter regions in the thin section is clear, being the pseudotachylite and magma components, respectively, of Kokelaar's (2007) description. Some of the larger clasts are polycrystalline fragments of quartzite but others, particularly in the lighter regions, appear to be magmatic phenocrysts of plagioclase (Fig. 3b and f). However, the colour variation marking out the flow banding at hand specimen scale appears to be largely due to variations in clast density, rather than differences in matrix colour, when viewed under the microscope (Fig. 3c). Smaller clasts are sub-angular or rounded and dominated by quartz, and to a lesser extent feldspar, and the semi-automated clast selection process identifies these fairly well (Fig. 3c and d). The exception is where mafic grains are included within a clast, as in the clast marked by a cross in Figure 3c. The mafic material is misidentified as pseudotachylite matrix, resulting in this clast being subdivided into smaller objects (Fig. 3d). This misidentification is relatively rare and is dealt with by manually tracing around larger clasts containing mafic components to ensure that they are identified as single objects.

Figure 4 shows the hand specimen and several photomicrographs of the Outer Hebrides sample. The gneiss host rock for this pseudotachylite contains appreciable amounts of amphibole and biotite but the semi-automatic procedure with manual tracing of larger clasts reproduces the clast contents of the photomicrographs well (Fig. 4b and c). The largest clasts in the hand specimen are much too large to be captured by the present analysis, being breccia clasts and not survivor clasts from a shearing and melting fault rupture. Clasts are angular to subrounded and many of the larger clasts are polycrystalline. The clast density in the Outer Hebrides sample is lower than in the flinty crush rock with a lower proportion of smaller clasts in particular.

Grain-size populations

Figure 5 presents the number of clasts, nr, plotted against effective radius for samples of pseudotachylite from the Outer Hebrides, flinty crush rock from Glencoe and the synthetic pseudotachylite. Fitting parameters for the power-law (equation 1b), modified power-law (equation 2a) and log-normal (equation 3) equations are given in Table 1.

In the case of the Other Hebrides pseudotachylite (Fig. 5a) there is a clear region of the clast-size population that is well fitted by a power-law distribution (black points). The power-law fit to the linear region in the data collected from photomicrographs at ×2.5 magnification extends over one order of magnitude, from 400 to 30 μm (black data in Fig. 5a), and provides a good fit to the ×10 magnification data to 6 μm, suggesting that clast size follows a fractal distribution in this sample. The fractal dimension derived from the fit to the ×2.5 data between 400 μm > r > 30 μm is D3d=2.6, within with the range of previous measurements for pseudotachylites. At clast sizes below 6 μm the data fall below the power-law fit, but the modified power-law equation (2a) fits the data well. All bar the largest six clasts in the ×10 data are very well fitted by this equation (bold black line in Fig. 5a), with a modified fractal dimension of 2.7 and a dissolved radius, r, of 2.3 μm. The log-normal function (equation 3; blue dashed line in Fig. 5a) fits the ×10 pseudotachylite data equally well at small clast sizes but shows downward curvature at large clast sizes, falling below the data for clast sizes greater than 20 μm.

The flinty crush rock (Fig. 5b) does not show a clear region of power-law behaviour but is convex-up for all clast sizes below 150 μm. The best-fitting power-law region gives a fractal dimension of D3d=3.1, but data collected at ×2.5 and ×10 magnification both clearly deviate from this below r = 20 μm. The two datasets collected at different magnifications are indistinguishable down to r = 10 μm, implying that the deviation from power-law behaviour is not related to resolution problems but is intrinsic to the population. The modified fractal dimension, D=4.0, is substantially larger than the simple power-law fit, consistent with the large dissolved radius, r=8.3μm, required by the fit.

For both samples the modified power-law and log-normal equations both fit the x10 magnification data well in the regions where there is resolution, with RSS of the log-normal fit approximately two times larger than for the modified power-law fit. However, if the fits are truly representative of the population then they should be applicable to the ×2.5 data with no modification as the two datasets were normalized by their analysed area. It is clear from Figure 5a that the modified power-law fit, but not the log-normal fit, extrapolates very well to larger grain size along the ×2.5 data in the region where there is resolution. In Figure 5b, however, the log-normal fit gives a slightly better fit, as seen in the reduced χ2 metric for the ×2.5 data (Table 1).

The synthetic pseudotachylite only includes data collected from a single image at a single magnification so no population normalization is necessary (Fig. 5c). As with the two natural samples there is significant upward convexity at small clast sizes, which is very well fitted by the modified power-law distribution at all clast sizes above 1 μm. In this case 1 μm is very close to the instrument resolution. The modified fractal dimension is identical within error to that of the flinty crush rock (D=4.0) but there is a smaller dissolved radius of r=4.5μm, as expected for this sample, which was molten for just the duration of shearing. The log-normal distribution is almost indistinguishable from the modified power-law distribution across the size range sampled by the data, deviating below the sample population and the modified power-law distribution only at radii greater than 77μm. This is seen in the RSS of the two fits, which are 4.1 and 7.3 respectively for the modified power-law and log-normal fits.

The only feature not explained by the fits in any of the three samples is the higher than expected number of large clasts in the flinty crush rock sample. This feature is also seen in our clast-size analysis of the photomicrographs of flinty crush rock presented by Kokelaar (2007) (Supplementary material, Fig. S1) suggesting that it is not just due to random fluctuations at low number density but rather is a property of the population distribution. This requires there to be a second population of objects (such as phenocrysts from a mixed-in magma component), which has a large mean grain size but relatively small standard deviation. This is consistent with the observation in thin section that several of the largest ‘clasts’ in the flinty crush rock do not appear to be survivor clasts but rather are primary magmatic feldspar phenocrysts (Fig. 3b and f).

Chemistry of the flinty crush rock matrix

Major element compositions of the matrix in the flinty crush rock sample are given in Supplementary material, Spreadsheet S2, along with best-fitting mass-balance calculations. Of the nine areas analysed, six were well fitted by the major minerals in the fault intrusion ‘granite’ from Garnham (1988) (with RSS ≤ 0.1), with very similar modal proportions (Ksp 0.41±0.04; Qtz 0.36±0.04; Plag 0.20±0.05; Bt 0.017±0.0003; Ox 0.01±0.0002). The main difference is that the analyses here contain slightly higher modal quartz (0.36 compared with 0.22 in the granite of Garnham 1988), possibly owing to inadvertent inclusion of quartz clasts in the analysed regions. The remaining three matrix analyses required significant muscovite to bring the residuals below 0.25, despite there being very little muscovite (<1%) in the granitic fault intrusion. The variation in fitted potassium feldspar phase fraction v. mica phase fraction is presented in Figure 6a. The data fall into two distinct groups, with five of the six ‘granite-like’ analyses clustering tightly around 0.43 Ksp and 0.01 Mica, whereas the remaining four analyses show a clear negative correlation between potassium feldspar and mica content, with φKsp=0.340.78(φBt+φMsc). This relationship is very similar to that seen in the global compilation of pseudotachylite compositions by Dobson et al. (2021), where constants of proportionality between potassium feldspar and total mica varied from −0.64 to −0.77. Dobson et al. (2021) took the relationship to be strongly indicative of peritectic melting of micas during fluid-absent pseudotachylite formation, and the present data support Kokelaar's (2007) suggestion that the flinty crush rock found around Stob Mhic Mhartuin is a product of mixing between pseudotachylite and rhyolitic fault intrusion magma.

It is reasonable to ask what was the protolith of the pseudotachylite, as the muscovite and biotite contents (0.08<φMsc<0.15;0.01<φBt<0.05) would require extreme enrichments of muscovite, but only little or no enrichment of biotite compared with the rocks currently hosting them (φMsc0.01;0.03<φBt<0.07). As Kokelaar (2007) observed, the pseudotachylite melt might have been transported for some distance before being deposited on the wall rock of the caldera fault system and so it is possible, in principle, for the protolith to differ from the wall rocks currently in contact with the quenched pseudotachylite. The A–K–F plot ((Al2O3–Na2O–K2O–CaO)–K2O–(FeO + MnO + MgO)) shown in Figure 6b is informative here. The granite and monzonite compositions of Garnham (1988), corresponding to the rhyolite and monzonite of Kokelaar (2007), which form the fault intrusion at Stob Mhic Mhartuin, plot close to the K–F vertex whereas the matrix compositions of the flinty crush rock are significantly enriched in the A-component. Pseudotachylites generated from granitic, tonalitic and basaltic protoliths are strongly enriched in the F-component, following the black arrow in Figure 6b (Dobson et al. 2021), meaning that remelting of an igneous protolith is not a viable way to produce the observed matrix compositions. The range of pseudotachylite compositions, from the compilation of Dobson et al. (2021), produced in pelitic and semipelitic host rocks is plotted on the A–K–F diagram and falls close to the biotite–muscovite mixing line, consistent with the dominance of mica melting in this type of pseudotachylite. The matrix compositions of the flinty crush rock analysed here all fall on mixing lines between this composition of pseudotachyite and the rhyolitic fault intrusion.

The magma mixing model, including semipelitic pseudotachylite, rhyolite and quartz clasts, produced good fits to all the measured matrix compositions (0.1<RSS<0.9;meanRSS=0.35), with pseudotachylite contents of the matrix (100×φPST/(φRHY+φPST)) of 5–41%. The four matrix regions that define the mica–potassium feldspar trend require the highest pseudotachylite contents of 16, 27, 32 and 41%. None of the fits were improved by addition of pelitic or psammitic pseudotachylite components, despite the increased degrees of freedom that addition of extra components provides. It seems likely therefore that the flinty crush rock at Stob Mhic Mhartuin is a product of mixing between pseudotachylite and rhyolitic fault intrusion magma, as suggested previously by Kokelaar (2007), and that the pseudotachylite component was generated in a region where the caldera fault passed through semipelitic host rocks.

Both the clast-size–shape metrics and the major element chemistry of flinty crush rock from Stob Mhic Mhartuin support Kokelaar's (2007) suggestion of a mixing origin between pseudotachylite and rhyolitic magma, with clear differences between this sample and pseudotachylite sensu stricto. The chemistry of all the analysed regions is well fitted, with RSS < 1, by a mixture of rhyolite magma (the porphyritic granite of Garnham 1988) and pseudotachylite generated from a semipelitic host rock (the mean of semipelite pseudotachylite compositions in the compilation of Dobson et al. 2021), with 0.05<(φPST/φRHY)<0.68. We reiterate that this is consistent with mixing between molten rhyolite magma and semipelitic pseudotachylite, rather than two pseudotachylite compositions derived from rhyolite and semipelite. A pseudotachylite derived from a rhyolite or monzonite wall rock would plot to the right of the biotite–muscovite line in Figure 6b (black arrow), owing to closed-system-like peritectic melting of micas, well away from the observed mixing trend.

The modified power-law equation, which takes account of removal of a constant rim thickness by dissolution into a melt, fits data from Stob Mhic Mhartuin flinty crush rock and from an Outer Hebrides pseudotachylite very well, but with some significant differences in the fitting parameters. Whereas the Outer Hebrides pseudotachylite showed a typical fractal dimension of 2.7 and a dissolved radius of 2.3 μm, the flinty crush rock sample had a much larger fractal dimension of 4.0 and a large dissolved radius of 8.3 μm. This large value of dissolved radius might be expected in the case where pseudotachylite is emplaced in intimate contact with magma, as this will inhibit the rapid quenching against cold wall rock that is normal for pseudotachylites, allowing more time for dissolution and re-equilibration of clasts with pseudotachylite matrix. The fractal dimension is larger than found by previous workers for pseudotachylites but, as previously discussed, the dissolution of clast material can have a profound effect on the apparent fractal dimension obtained by simple straight-line fitting of the log–log plot. The fractal dimension of clasts in the synthetic granite pseudotachylite of Montheil et al. (2020) increases from 2.9 to 4.0, in remarkable agreement with the values for flinty crush rock presented here when equation (1a) is replaced by equation (2) (Fig. 5; Table 1). In the pseudotachylite and flinty crush rock samples the modelled dissolved radius is large compared with the maximum clast size (4–5%), compared with the Outer Hebrides sample (0.7%), and hence the apparent fractal dimension obtained from straight-line fitting significantly underestimates the true value of the clast population prior to modification by partial melting. In the case of highly modified fractal distributions (i.e. where the dissolved radius is a significant fraction of the maximum clast size) the modified fractal distribution becomes difficult to distinguish from a log-normal distribution.

The mixed sources that contribute to the flinty crush rock further complicate the situation, as the rhyolite component introduced magmatic phenocrysts into the pseudotachylite mixture (e.g. Roberts 1966; Kokelaar 2007) and the semi-automated clast analysis procedure adopted here does not distinguish between rhyolite phenocrysts and pseudoctachylite clasts. Kokelaar (2007) presented four photomicrographs of regions of intermixed flinty crush rock and rhyolite from Stob Mhic Mhartuin and we have performed the same clast-size analysis on these images, treating regions of rhyolite and flinty crush rock separately (Supplementary material, Fig. S1a). Although the resolution of the images analysed is poor, as reflected in the very large r* value for the dissolved radius, the fractal dimensions for rhyolite-dominated and flinty crush rock-dominated regions are rather similar and agree with the large value determined for the sample studied here. We also note that Kokelaar's (2007) images and the flinty crush rock sample analysed here both show significant deviation to larger values for the largest few clasts. Although the numbers of these large clasts are too small for statistically robust fitting, this finding is consistent with there being a population of large phenocrysts mixed in with the survivor clast population. Supplementary material, Figure S1b shows how a two-population fit reproduces the data (modified power-law distribution similar to the fits shown in Fig. S1a plus a log-normal distribution centred at large grain size). We suggest that this is also consistent with the photomicrographs of flinty crush rock, where the largest clasts appear to be primary magmatic feldspar phenocrysts (squat-tabular shape with some resorption rounding of corners; no undulose extinction; simple magmatic twinning; Fig. 3b and f). Given the complexity of the mixed magma system, we conclude that there is no evidence from the clast-size analysis of flinty crush rock that would suggest an origin other than as fault-generated pseudotachylite. Combined with the clear chemical evidence for mixing between pseudotachylite and rhyolite we conclude that a frictional melting origin, combined with magma mixing during emplacement, is most likely for the the Glencoe flinty crush rock at Stob Mhic Mhartuin.

It remains to consider the source of the semipelitic pseudotachylite component. The rocks hosting the ring faults and fault rock intrusion in the vicinity of Stob Mhic Mhartuin are Eilde Flags (micaceous psammites) to the north and quartzite to the south in the downthrown block. Neither of these contains sufficient pelitic material to produce the pseudotachylite compositions seen here and, given the likely emplacement mechanism of the fault intrusion and associated flinty crush rock mixture, it seems likely that the pseudotachylite component was generated at depths below the present land surface. The stratigraphy underlying the Eilde Flags is unclear, as these are the oldest of the Grampian rocks exposed locally; however, the general trend is of decreasing pelite content with increasing age, suggesting that the semipelite source rocks for the pseudotachylite are not stratigraphically below the Eilde Flags. The structure of the region is complex, with several conflicting structural interpretations (e.g. Stephenson et al. 2013); however, all studies agree that Eilde Flags in the area immediately north of here overlie rocks in the overturned limb of the recumbent Kinlochleven Antiform–Treig Syncline system (e.g. Treagus 1974; Stephenson and Gould 1995). The next younger rocks in the Grampian stratigraphic succession here are the Eilde Schist group, which would provide ample semipelitic source material for the pseudotachylite component. The presence of semipelite-sourced pseudotachylite might therefore imply that, at depth, the caldera-bounding ring faults at Stob Mhic Mhartuin cut and displace younger rocks of the overturned limb of the Kinlochleven Antiform–Treig Syncline system.

Conceptual model for evolution of clast size in a molten pseudotachylite

The present study has demonstrated rational relationships between the curvature of clast-size distributions at small radius and the amount of melting in the pseudotachylite. The good fits produced by assuming that a constant radius is removed from all clasts implies that clast melting occurred by dissolution into the pseudotachylite melt, rather than by diffusion of heat into, and isochemical melting of, clasts. It further implies that diffusion gradients did not form around clasts while this was occurring, implying that most of the clast dissolution occurred while seismic slip was still active (see Dobson et al. 2018, for further discussion of diffusive dissolution of survivor clasts). We expect the clast-size population on a fault to evolve during seismic slip as shown schematically in Figure 7. Prior to melting, ultracataclasis results in a high fractal dimension and power-law distribution to small clast sizes (stage a, Fig. 7). On first melting of the least refractory minerals, survivor clasts start to dissolve into the melt (stage b, Fig. 7). This modifies the clast-size population, resulting in upward convexity at smallest clast sizes. It is thought that injection veins form during this stage of pseudotachylite formation in response to the volume of melting, causing increased overpressure of the system. Any injection veins would rapidly quench in the cold country rock (e.g. Dobson et al. 2018), preserving the clast-size population at this stage. Continued shear maintains melt on the fault plane, resulting in more dissolution of survivor clasts and more extreme modification of the size population until seismic shear ends and the system quenches (stage c, Fig. 7). Quenching rates on the fault plane will depend on the size of the pseudotachylite body and on the duration of the seismogenic rupture (Dobson et al. 2018). In tectonic fault systems one might expect to encounter clast populations from stages (a) to (c) and deviations of clast-size populations from power-law behaviour might provide another indicator of rupture duration (or at least time that the pseudotachylite was molten). The cooling history of pseudotachylite-hosted breccias, such as the Outer Hebrides sample studied here, might be complex, given that fault-generated melt is isolated from the heat source but also has the potential to collect into large volumes. In volcanic fault systems interactions with magma become possible. In stage (d) (Fig. 7), the magma system will have its own phenocryst population with a different size distribution. In stage (e) (Fig. 7), mixing between the magmatic and pseudotachylite melt will delay freezing of the pseudotachylite, further modifying the survivor clast population. Additionally, the overall ‘clast’-size distribution on the fault plane will be a combination of the magmatic and pseudotachylite populations, as seen in the flinty crush rock at Glencoe.

Clast-size metrics and matrix compositions of the flinty crush rock at Stob Mhic Mhartuin suggest a formation history that is more complex than simple frictional melting and quenching on the fault plane(s) of an erupting caldera system. The combined facts of intense mixing between pseudotachylite and rhyolite components and the likely semipelitic wall-rock composition for the pseudotachylite component imply that the flinty crush rock was transported some distance prior to emplacement. This is also consistent with the large melted-clast radius required by the modified power-law fit, which implies that the flinty crush rock was not rapidly quenched against cold wall rocks within seconds to minutes of its formation (e.g. Dobson et al. 2018). The current findings are consistent with the model of Kokelaar (2007), which involves mobilization of a (frictional melting) pseudotachylite by, and mixing with, an expanding foam of rhyolitic magma during early stages of caldera collapse eruptions. The demonstration of a rational relationship between dissolved radius and amount of heat available in the molten system raises the possibility that systematic studies of dissolved clast radius might shed light on the time that a pseudotachylite body remained molten, and hence the duration of seismogenic rupture. This will be investigated in a future study.

We thank T. Mitchell for helpful discussions, and A. D. Hollinsworth and an anonymous reviewer for constructive reviews, which significantly improved the paper.

DPD: conceptualization (lead), data curation (lead), formal analysis (lead), funding acquisition (lead), investigation (equal), methodology (lead), project administration (lead), supervision (lead), writing – original draft (lead), writing – review & editing (lead); VM: formal analysis (supporting), investigation (equal), methodology (supporting), writing – original draft (supporting)

This work was funded in part by the Natural Environment Research Council (NE/X009807/1).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All data generated or analysed during this study are included in this published article and its supplementary information files.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (