Crystal Size Distribution (CSD) Analysis of Volcanic Samples: Advances and Challenges

Introduction

The crystal population of volcanic samples represents the integration of processes that occur throughout magmatic systems. Reconstructing these processes from the crystal record is challenging, however, and requires integration of compositional and physical analysis of both individual crystals and crystal ensembles. A wide array of tools is now available for such analysis; the challenge lies in integrating analysis of crystal populations with evolving understanding of the physical processes that control magma evolution (e.g., Jerram et al., 2018) and modulate volcanic activity (e.g., Martel et al., 2019). Addressing the former requires unraveling complex textures and compositional zonation preserved in the phenocryst (macrocryst) population; addressing the latter involves interpreting groundmass measurements using kinetic constraints provided by laboratory experiments.

The past decades have seen accumulating evidence that many of the crystals present in volcanic samples did not grow exclusively from the melt that transported them to the Earth’s surface (e.g., Nakada and Motomura, 1999; Davidson et al., 2005, 2007; Humphreys et al., 2006, 2008; Smith et al., 2009; Viccaro et al., 2010; Cashman and Blundy, 2013; Neave et al., 2013). Instead, crystal cores commonly preserve records of formation in diverse parts of the magmatic system and complex histories of growth, reaction, resorption and agglomeration. For this reason, the phenocryst population (sensu lato) of a volcanic sample is often separated into indigenous phenocrysts that have grown entirely from the transporting melt and entrained exogenous antecrysts and xenocrysts. Many exogenous crystals, however, have mantles and/or rims that are in equilibrium with the melt and thus have a mixed exogenous-indigenous origin. Signatures of crystal growth histories are easily mapped via stunning advances in compositional imaging (e.g., QEMSCAN, Neave et al., 2014, 2017), tomographic imaging (e.g., Gualda et al., 2004; Mock and Jerram, 2005; Pamukcu and Gualda, 2010; Jerram et al., 2018; Polacci et al., 2018; Tripoli et al., 2019), microstructural analysis (e.g., Kahl M. et al., 2017; Holness et al., 2019) and isotopic analysis (e.g., Davidson et al., 2007; Morgan et al., 2007). Textural (size and shape) analysis of phenocryst populations, however, has yet to capitalize on these advances, particularly from the perspective of analyzing and interpreting the crystallization histories of crystals with mixed origins.

The past few decades have also seen a dramatic increase in both analysis of groundmass textures and experimental studies of crystallization relevant to syn-eruptive volcanic processes. Measurements of groundmass crystal populations can be linked to volcanic processes where sample times are well constrained (e.g., Cashman et al., 1999; Hammer et al., 1999, 2000; Cashman and McConnell, 2005; Wright et al., 2012; Preece et al., 2016; Harris et al., 2020) but textural interpretations must rely on experimental constraints when temporal data are not available. Experimental studies, in turn, underline the sensitivity of groundmass crystal populations to initial conditions (super- or sub-liquidus) and melt composition as well as cooling and decompression paths (e.g., Brugger and Hammer, 2010a, b; Martel, 2012; Riker et al., 2015a; Befus and Andrews, 2018).

Here I illustrate ways in which analysis of crystal textures can be linked to dynamic views of crystallization processes in volcanic systems. I focus primarily on studies of plagioclase crystallization, as this phase is both the most common and the most commonly analyzed; an excellent review of crystallization of mafic phases is provided by Hammer (2008). Context is provided by a brief review of the basic tenets of textural analysis and the relation between sample textures and the crystallization pathways that produced them. This is followed by an overview of cooling-driven crystallization of basaltic lava on the Earth’s surface, where experimental results can be linked to textural and compositional measurements of well-constrained samples. Syn-eruptive decompression-driven crystallization is more challenging to study in both experimentally and active volcanic systems; studies in these systems show broad agreement but highlight the need for experimental constraints on very shallow subvolcanic processes. Unraveling complex histories of phenocryst populations poses more profound challenges that will require, ultimately, combining techniques of isotopic analysis, diffusion chronometry and compositionally based crystal size distribution (CSD) and shape analysis.

Background

The basic concepts of crystallization and crystal size distribution (CSD) analysis are well reviewed by Marsh (1998) and Hammer (2008) and are only briefly reviewed here. Crystallization requires both crystal nucleation and growth of those nuclei. Rates of crystal nucleation and growth are functions of magma supersaturation, which is commonly described by the effective undercooling (ΔT_eff), or the temperature difference between the magma and the saturation temperature of the phase in question.

ΔT_eff is a straightforward concept for cooling-driven crystallization; it is less intuitive for decompression-driven crystallization, where supersaturation is more logically formulated as Δφ, the difference between the reference crystallinity (e.g., the crystallinity at the starting pressure) and the crystallinity at some lower pressure (Brugger and Hammer, 2010a, b; Riker et al., 2015a, b; Befus and Andrews, 2018). Importantly, in multicomponent systems supersaturation is not a fixed value, but varies with temperature and/or pressure. This is illustrated in Figure 1A, which shows equilibrium Δφ as a function of pressure for decompression of different H₂O-saturated bulk compositions. Experimentally, Δφ can be imposed by rapid decompression from liquidus conditions to a final pressure (P_f). Under these conditions, rates of crystal nucleation (J) peak at larger Δφ than rates of crystal growth (G), as illustrated in Figure 1B, where the ratio J/G is shown as a function of P_f. At a given P_f, the rate of crystallization is time-dependent, decreasing as the system approaches equilibrium. For this reason, measured nucleation and growth rates will vary depending on the time interval of measurement. For example, time-sequential measurements of crystal size (filled circles and solid line in Figure 1C) yield instantaneous growth rates (dashed lines) that are both faster and slower than the time-averaged value derived from the final experiment alone (line labeled apparent G).

The relative rates of crystal nucleation and growth determine the final texture of a sample, where here I use the term texture to mean the abundance and size distribution of the constituent crystals. For the same total crystallinity, high nucleation rates produce numerous small crystals, while high rates of growth on a limited number of sites produce fewer but larger crystals. Theoretical studies of crystallization provide expressions for crystal nucleation and growth that can be parameterized to link crystal size (L) and number density (N_v, per volume) to time-averaged rates of J and G (e.g., Brandeis and Jaupart, 1987):

In both expressions, the constant of proportionality is of order 1 except where nucleation is an exponential function of time (Marsh, 1998). This parameterization illustrates the strong control of J/G on crystal number and can be linked directly to experimental data (Figure 1B).

Quantitative Textural Analysis

The simplest (and often most expedient) form of textural analysis relates measured crystal number density (N_a, number per area) to the total area fraction occupied by the crystals (φ_a) through the average crystal area (d², where d is a measure of crystal size) as

Where progressive crystallization is achieved primarily by adding new crystals (nucleation), φ_a increases linearly with N_a and the slope of the line yields the average crystal area (1/slope = d²; Hammer et al., 1999). In contrast, where an initial period of nucleation is followed by crystal growth, crystallinity increases at a near-constant (or even decreasing) N_a. A plot of N_a vs. φ_a for sequentially crystallized samples can thus be used to determining the relative importance of nucleation and growth in driving crystallization (Hammer et al., 1999, 2000).

Measurements of N_a and φ_a in two dimensions can be related to three-dimensional measurements through use of simple stereological assumptions. When crystals are randomly oriented, the measured area fraction is directly equivalent to the volume fraction (φ_v). The number per volume (N_v) can be related to the measured number per area (N_a) as

if the average size (d_av) is known. Since the average crystal size can be determined from Eq. (3; d_av = √d²), N_v for each sample may also be determined (N_v = N_a^1.5/ φ_a^0.5).

Crystal size distributions analysis provides a more complete view of a sample’s crystallization history. CSD analysis was developed for steady state crystallizers (e.g., Randolph and Larson, 1988). Here the desired output CSD is controlled by (1) the seed crystals in the input solution, (2) the time spent in the crystallizer, and (3) the rate at which crystals nucleate and grow within the crystallizer (controlled by the supersaturation). Under steady conditions, the crystal number density n (the slope of the cumulative distribution, with units of number per size class per volume) is related to crystal size L as

where n°, the nucleation density, is the intercept at L = 0; L_d, the dominant size, is measured by the slope of the distribution (L_d = −1/slope = −1/Gt). Mathematically, L_d is the mode of the size-based distribution (the first moment of the integral form of Eq. 5) and is analogous to the half-life described by the radioactive decay equation; the equivalent half-size would be L_d/ln 2. In an industrial crystallizer, t is known and G can be calculated. The nucleation rate J is related to the nucleation density as J = n°G, the total number of crystals is n°Gt, and the total volume fraction (V_T) of the measured phase is

where k_v is the shape factor required to convert L_d³ to the true three-dimensional volume of a single crystal.

The application of CSD analysis to volcanic samples presents numerous challenges. The first relates to definition of crystal size. In thin section, crystal size can be measured variously as the short axis, long axis, √area or average intersection length; the choice of size metric affects both the slope and the intercept of the resulting CSD (Cashman, 1988; Muir et al., 2012). A second challenge relates to conversion of data collected on 2D thin sections to the 3D formulation of CSD analysis. A solution to this problem is provided by the widely used program CSDCorrections, where a constant average crystal shape (determined from CSDSlice) is used to correct for both the cut effect and the intersection probability (Higgins, 2006). An alternative approach is to make measurements directly in 3D. Techniques to obtain direct 3D measurements have included deconstruction of vesicular volcanic material to extract individual crystals (Dunbar et al., 1994; Gualda et al., 2004), serial sectioning of samples (Castro et al., 2003; Duchene et al., 2008), and, more recently, tomographic techniques that allow non-destructive 3D measurements (e.g., Gualda et al., 2004; Mock and Jerram, 2005; Gualda and Rivers, 2006; Jerram et al., 2010; Pamukcu and Gualda, 2010), and even 4D (e.g., Polacci et al., 2018; Tripoli et al., 2019).

Direct comparison of CSDs calculated from 2D and 3D measurements shows that the cumulative area distribution in 2D closely approximates the cumulative volume distribution in 3D (Jerram et al., 2009). This observation lends support to the definition of average crystal size as d = √area (Eq. 3). CSD slopes and number densities are generally in good agreement except at the smallest sizes, where CSDs calculated using CSDCorrections may show high number densities that contrast with downturns in the number of small crystals measured in 3D (Castro et al., 2003). One explanation for this discrepancy is that crystals of different sizes have different shapes. For example, prismatic pyroxene microlites measured by Castro et al. (2003) become more anisotropic with increasing size, while the tabular plagioclase crystals imaged by Duchene et al. (2008) are more equant with increasing size. More generally, the response of specific morphologies to sectioning affects the accuracy of shape correction. Tabular crystals, for example, are more easily recognized and quantified than prismatic crystals because of statistical problems associated with intersecting the longest axis in very acicular forms (Morgan and Jerram, 2006); this problem is further complicated in mixed crystal populations. The sensitivity of stereological corrections to crystal shape suggests caution when calculating CSDs over a large size range using a single shape factor (see also Mock and Jerram, 2005; Brugger and Hammer, 2010a).

More fundamental challenges to CSD analysis lie in the extent to which conditions of crystallization in natural systems deviate from those in steady state industrial crystallizers. Even in simple volcanic systems, crystallization does not occur at fixed undercoolings (Figure 1), or with constant crystal and melt compositions, and the final crystal size distribution preserved in any sample necessarily represents an integration of nucleation and growth processes over time. This means that extraction of kinetic parameters is possible only for samples with independent time constraints. Other complications relate to the effects of simultaneous crystallization of multiple phases, which may increase the potential for heterogeneous nucleation and as well as the likelihood of crystal-crystal impingement (e.g., Zieg and Lofgren, 2006; Holness et al., 2007). The presence of pre-existing crystals also affects the conditions under which a new batch of nuclei will form (Fokin et al., 1999). Below I explore these challenges through the lens of recent experiments and their applications to the interpretation of crystal size characteristics of volcanic samples.

Cooling-Driven Crystallization of Basalt at 1 Atm

Experimental investigation of cooling-induced crystallization at 1 atm has a long history, reaching back to the Neptunist-Plutonist arguments of the late 18th and early 19th centuries. By the early 20th century, establishment of the Hawaiian Volcano Observatory had prompted volcanologists to examine the chemical and physical properties of basaltic lava and controls on flow emplacement. Included in this work was investigation of the contrasting surface morphologies of pāhoehoe and ‘a‘ā flows and the observed along-channel transformation of lava from pāhoehoe to ‘a‘ā with increasing transport distance. Emerson (1926) addressed this problem by melting lava in a forge and stirring it as it cooled, experiments that demonstrated the importance of shear in generating the numerous small crystals that characterize ‘a‘ā lava.

Cooling Rate Experiments

Kinetic experiments on low viscosity basaltic melts performed at 1 atm (e.g., Kouchi et al., 1986; Sato, 1995; Pupier et al., 2008; Vona and Romano, 2013) provide quantitative insight into the crystallization kinetics of lava flows on the Earth’s surface. First, experiments that start at super-liquidus conditions show small changes in initial temperatures can have profound impacts on the final sample textures (Sato, 1995; Vetere et al., 2013). Second, for similar initial conditions, cooling rate exerts a strong control on plagioclase number density (and size) but does not affect the crystallinity, which rapidly approaches equilibrium values (Pupier et al., 2008). Deformation experiments (Kouchi et al., 1986; Vona and Romano, 2013; Tripoli et al., 2019) have validated Emerson’s (1926) conclusions about the role of stirring in generating high crystal nucleation rates, interpreted to reflect the role of melt advection during stirring, although mechanical breakage to form secondary nuclei may also be important.

Textural analysis of cooling rate experiments highlights additional controls on final sample textures. Experiments initiated at superliquidus temperatures (Pupier et al., 2008) produce N_a-φ trends positive (low cooling rates) to negative (high cooling rates) slopes and CSDs that “fan” with decreasing quench temperature (increasing crystallinity; Figure 2). The reduction in crystal number accompanying increasing crystallinity requires a substantial increase in crystal size as well as loss of some crystals; both signatures can be explained by crystal agglomeration (synneusis; Schwindinger and Anderson, 1989). Alternative explanations come from recent innovations in 1 atm experiments (e.g., Schiavi et al., 2009) that allow in situ observations of olivine growth (Ni et al., 2014). In these experiments, a decrease in the number of early-formed nuclei is attributed to Ostwald ripening, consistent with results of Cabane et al. (2005) that show ripening is most effective for very small crystals. Also important is size-dependent growth, such that large crystals grow faster than small crystals (e.g., Eberl et al., 2002; Kile and Eberl, 2003). In the experiments of Ni et al. (2014), size-dependent growth is attributed to melt advection, as measured growth rates are faster than those expected for diffusion alone. Size-dependent growth also produces fanning CSDs, as required to conserve the total number of crystals (=n°Gt).

Natural Cooling Experiments

Natural examples of cooling-driven crystallization are provided by studies of crystallization in Hawaiian lava flows (e.g., Cashman et al., 1999; Soule et al., 2004; Riker et al., 2009; Cashman and Mangan, 2014). Here time constraints are provided by samples collected along active lava channels on the same day (e.g., Crisp et al., 1994; Cashman et al., 1999), flow cooling is measured by glass geothermometry and eruption of lava at near-liquidus temperatures means that crystallization is driven by cooling accompanying flow. Rates of cooling (∼18°C/hr) and crystallization (∼1.8 φ/hr) are obtained by combining distance with known flow velocities along proximal to medial flow reaches. Temporal rates of both cooling and crystallization appear to be independent of eruption rate, which suggests that textural analysis of older flows can be used to infer eruption rates (e.g., Riker et al., 2009).

Crystallization of Hawaiian lavas is dominated by cotectic precipitation of plagioclase and pyroxene, which typically form crystal clusters that indicate the prevalence of heterogeneous, rather than homogeneous, nucleation (Figure 3A). Textural analysis of along-channel samples shows positive correlations between N_v and φ for both phases (Figure 3B), consistent with rapid and nucleation-dominated crystallization (e.g., Cashman et al., 1999; Riker et al., 2009). CSDs are linear, as expected for continuous cooling from near-liquidus conditions, and, when sample locations and flow rates are well constrained, yield nucleation rates of ∼50–150/mm³s and growth rates ∼10^–7 mm/s (Cashman and Mangan, 2014). These high rates likely reflect strong advection driven by efficient heat loss at the channel margins (e.g., Cashman et al., 2006).

The nucleation-controlled crystallization observed in open channel samples contrasts with the growth-dominated crystallization observed in insulated pahoehoe flows and lava lakes (Cashman and Marsh, 1988; Cashman and Mangan, 2014). The latter show both decreasing N_v with decreasing φ (Figure 3B) and fanning CSD trends similar to those produced by super-liquidus experiments (Pupier et al., 2008; Ni et al., 2014; Figure 2). It seems likely that crystal agglomeration, loss of the smallest crystals via Ostwald ripening and size-dependent growth rates all contribute to these textural trends. Mafic enclaves may also show linear CSDs with characteristics (n°, Ld; e.g., Martin et al., 2006a) that are similar to the most crystalline samples from lava lakes (e.g., Cashman and Mangan, 2014). The implied slow cooling suggests that crystallization of these enclaves occurred before, not after, entrainment in the melt. Pre-entrainment crystallization is further suggested by the lack of correlation between measured textures and enclave size (Coombs et al., 2003) and suggests that existing plagioclase frameworks may help to preserve enclaves during the entrainment process (Martin et al., 2006b; Andrews and Manga, 2014).

Decompression-Driven Crystallization of Hydrous Magmas

Questions about conditions leading to explosive vs. effusive eruptions (e.g., Cassidy et al., 2018) have driven an explosion of experiments designed to characterize decompression-driven crystallization of hydrous magma (e.g., Hammer, 2008). These experiments span a wide range of compositions and initial conditions (P, T, and X_H2O), although most have been conducted on H₂O-saturated silicic melts. Plagioclase is the most abundant crystallizing phase, and the most universally quantified for crystal size and shape. Experimental aims are to determine rates of crystallization, often with the goal of constraining magma ascent paths for specific eruptions. The past two decades have also seen an impressive increase in the number of textural studies of volcanic samples, many of which have been done in conjunction with experimental studies. Although very rapid decompression prevents crystallization, syn-eruption crystallization is observed in products of (1) short explosions that mark the buildup to large Plinian eruptions, (2) cycles of Vulcanian activity, and (3) episodic or continuous lava effusion (including ash venting and dome collapse).

Decompression Experiments

Decompression experiments are conducted using (1) single-step decompressions (SSD), which involve an initial rapid pressure decrease followed by time at the final pressure, (2) multi-step decompressions (MSD), where steps and dwell times are varied to approximate different decompression rates, and (3) continuous decompression (CD). SSD experiments typically produce high N_v (high initial nucleation rates) but require annealing at final pressures to achieve equilibrium crystallinity values; nucleation and growth rates reported for these experiments vary depending on anneal duration (e.g., Befus and Andrews, 2018; Figure 1C). MSD experiments are difficult to generalize as they vary in both step size and hold time at each step. CD experiments (Brugger and Hammer, 2010a, b; Riker et al., 2015a, b; Waters et al., 2015; Befus and Andrews, 2018) are not annealed at the final pressure (P_f), so that the crystallization time is controlled by the decompression rate and pressure interval of interest.

Decompression experiments have been performed over time periods of hours to weeks; this means that they are most relevant to syn- and short intra-eruptive periods of crystallization. Resulting plagioclase textures span a wide range that reflects variations in melt composition, initial conditions and rates and modes of decompression. Experimental data are shown in Figure 4; here N_v and φ are used as plotting parameters as they are the easiest textural parameters to measure and correlate well with the results of more detailed CSD analysis (Brugger and Hammer, 2010a). Experimental crystallinities are typically less than 0.3 and number densities range from 10² to 10⁷/mm³; they vary most strongly with P_f and melt composition such that low P_f, low temperatures and silicic compositions produce the highest N_v.

As observed in cooling experiments, the initial conditions matter. Super-liquidus starting conditions inhibit crystal nucleation and promote crystal growth relative to sub-liquidus samples run at the same conditions (Martel, 2012; Waters et al., 2015). Pre-existing crystals, in contrast, provide sites for crystal growth and therefore allow crystallization even when Δφ is sufficiently small to prohibit nucleation of new crystals. Growth on pre-existing crystals can be significant. For example, if growth rates are linear and similar on all major growth faces, adding rims to large crystals will produce more efficient crystallization, volumetrically, than adding new small crystals (Befus and Andrews, 2018). One consequence is development of those pre-existing crystals as a separate population that adds pronounced curvature to CSD plots (Brugger and Hammer, 2010b; Riker et al., 2015b; Figure 5).

Decompression-Crystallization During Eruptions

A compilation of groundmass plagioclase data from pyroclasts and lavas produced by recent eruptions is shown in the N_v-φ plot of Figure 6. The data are from a range of bulk compositions and eruption styles, span wide range of N_v (10² to > 10⁸ mm^–3) and include variations of φ from near 0 to 0.8. Although all samples were quenched at the surface, they are interpreted to reflect a range of P and T conditions that record different conduit processes. In detail, rhyolitic samples from Pinatubo, Philippines (PIN), have the lowest crystallinities while basaltic andesite samples from Merapi, Indonesia (MER), and Spurr, Alaska have the highest. Pyroclasts from an individual eruption sequence can show a wide range of crystallinities, but less common is a wide range in N_v.

Comparison of Figures 4, 6 shows good agreement between experimental data and many of the sample suites, thus demonstrating the power of combining these two approaches to constrain eruption conditions. Where experiments attempt to replicate textures of natural samples (e.g., Mount St. Helens – Muir et al., 2012 and Riker et al., 2015b; Pinatubo – Hammer et al., 1999 and Befus and Andrews, 2018) the CSDs are similar in form but experiments fail to reproduce the large number of small microlites observed in some natural samples. This mismatch can be seen in Figure 6, where number densities of natural samples can exceed experimental values by 1–2 orders of magnitude.

Eruption Conditions That Produce Very High N_v

Very high groundmass crystal number densities were first measured in pyroclasts from eruptions that preceded the 1991 climactic eruption of Pinatubo, Philippines, the largest well-documented volcanic eruption in the 20th century. The climactic eruption was preceded first by a dome and then 3 days of explosive activity that decreased in duration and increased in frequency approaching the climactic eruption (Hoblitt et al., 1996); this pattern of activity has been interpreted to record construction of a connected conduit from the magma reservoir to the surface (Scandone et al., 2007). The very high plagioclase number densities occur in the most crystalline samples and increase in maximum crystal size with increasing inter-eruptive repose time (Hammer et al., 1999; Figure 6). Crystals in individual samples form linear CSDs with very high n° and steep slopes (Figure 7A). A limited range in matrix glass H₂O records magma arrest in the conduit at effective pressures of c. 8–16 MPa. Together the textural data and glass volatile contents suggest crystallization by rapid decompression during an explosion, followed by annealing at shallow levels in the conduit during the short inter-explosion repose periods (28–262 min) and expulsion during the next explosion.

The high number densities of Pinatubo pyroclasts are not unique to this eruption but have also been observed in pyroclasts associated with the P1 Plinian eruption of Mont Pelee (650 ybp; Martel and Poussineau, 2007) and in some samples from the 1980 eruptions of Mount St. Helens (Klug and Cashman, 1994; Cashman and McConnell, 2005; Figure 6). As at Pinatubo, the high N_v pyroclasts likely record rapid decompression followed by short (hours) anneal times. Interestingly, very high crystal numbers are also reported in cryptodome samples from Mount St. Helens (MSH; Cashman and Hoblitt, 2004), in surge deposits from Mont Pelee (Martel and Poussineau, 2007) and in block-and-ash flow deposits formed by collapse of the Soufriere Hills dome in 2010 (Murch and Cole, 2019; Figure 8). Although individual pyroclasts from MSH reach higher overall crystallinities than those from Pinatubo (Figure 8), crystal populations form linear CSDs that have a similar range of intercept and slope (Figure 7B). The very high crystal number densities in dome and cryptodome samples are puzzling. Both observational and experimental constraints show that effusive eruption and emplacement of silicic lava domes (and cryptodomes) requires very slow magma ascent to very low pressures. Under these conditions, we might expect extensive decompression-driven crystallization at relatively low Δφ. How is it possible, then, to generate a uniform, finely crystalline but still melt-rich groundmass?

Complex Crystallization in Volcanic Systems

Overwhelming documentation of the prevalence of antecrysts in magmatic systems (e.g., Bacon and Lowenstern, 2005; Humphreys et al., 2006; Smith et al., 2009; Viccaro et al., 2010; Neave et al., 2013, 2014, 2017; van der Zwan et al., 2013; Witter et al., 2016) has been critical for the development of mush-dominated paradigms of magmatic systems (e.g., Edmonds et al., 2019; Sparks et al., 2019). Entrained crystals can be identified by major element and/or isotopic compositions that are not in equilibrium with the transporting magma (e.g., Morgan et al., 2006, 2007; Berlo et al., 2007; Martin et al., 2010; Sakyi et al., 2012; Cashman and Blundy, 2013; Sides et al., 2014; Brenna et al., 2018; Shane et al., 2019) and by evidence of deformation in crystal-rich mush zones (e.g., Kahl W. et al., 2017; Befus et al., 2019; Holness et al., 2019; Wieser et al., 2020). These crystals are often mantled by crystal rims that are compositionally and/or texturally distinct from the entrained core and typically in equilibrium with the carrier liquid.

Complex crystallization histories are particularly well recorded by plagioclase crystals. Examples from Mount St. Helens include inherited cores with a range of compositions and dissolution/reaction textures (Figure 9), mantles with complex oscillatory zones (Figure 9A), rims that are distinct (Figure 9B) to almost non-existant (Figure 9C) and evidence of agglomeration both before (Figure 9B) and after (Figure 9D) final rim growth. Importantly, patterns of rim growth may be mirrored by adjacent microphenocrysts (Figures 10A,B), demonstrating that both formed in response to the same decompression conditions, in this case multi-stage decompression in the volcanic conduit (e.g., Cashman, 1992; Geschwind and Rutherford, 1995; Blundy and Cashman, 2005; Cashman and McConnell, 2005).

Although crystal rims associated with decompression are often thin, they may contribute significantly to the total volume of both individual crystals (e.g., Befus and Andrews, 2018) and the total amount of decompression-driven crystallization (Riker et al., 2015b). For example, Figure 10C shows the relation between crystal volume and rim volume for the same rim thickness but variable crystal shapes and Figure 10D shows the relation between area % and volume % of rims for variable shapes and rim thicknesses. Specifically, measurement of core and rim areas for the phenocryst (Ph) and microphenocryst (MP) in Figure 10A show that the rim comprises 16% of the phenocryst area and 46% of the microphenocryst area; this illustrates the relative importance of adding approximately the same linear thickness of rim growth to a larger and smaller crystal. The absolute addition of crystal area via rim growth, however, is ∼3.6 times larger for the phenocryst (2577 μm²) than for the microphenocryst (710 μm²); the volumetric contribution will be substantially larger. From the perspective of crystallization history, therefore, the growth of phenocryst rims can contribute substantially to the amount of decompression-driven crystallization, which is typically assessed via groundmass crystallization alone.

Alternative Approaches to CSD Analysis of Complex Crystal Populations

The diverse and complex crystal populations described above raise important questions about current approaches to CSD analysis. It is clear that the phenocrysts (sensu lato) often record more than one crystallization episode and may have time breaks, for example at resorption boundaries. Crystals in explosively erupted pyroclasts may also be extensively fragmented (e.g., Bindeman, 2005; Jerram et al., 2009; van Zalinge et al., 2018). Crystal cores may be drawn from different populations and therefore do not always represent a single crystallization event. Finally, crystal rims represent a late stage of crystallization which is often correlative with formation of the groundmass population. Therefore, if the intent of CSD analysis is to decipher crystallization histories, interpretation of CSDs by absolute crystal size may not be adequate.

Several alternative approaches to textural analysis have been suggested. Early stages of crystallization may be preserved within oikocrysts or megacrysts; measurement of these crystals is particularly useful for determining the initial crystallization conditions of plutonic rocks (e.g., Higgins, 2017). Patterns of crystal zoning can be used to identify different crystal populations and magma recharge events, which in turn can be combined with diffusion analysis to constrain the timing of magma inputs (e.g., Morgan et al., 2004, 2006; Kahl et al., 2013). Isotopic microdrilling can be used to identify individual crystallization episodes and to synchronize individual crystals for use as chronometers (e.g., Martin et al., 2010); when combined with CSD analysis (iCSDs; Morgan et al., 2007), identified crystallization events can be tied to episodic magma recharge. These methods provide important constraints on the relation between crystallization and magma input; none, however, isolates the late-stage crystallization responsible for growth of both phenocryst rims and groundmass crystals.

Here I suggest an additional approach that complements the applications of isotope tagging and diffusion chronometry described above. When crystal cores and rims/groundmass differ in composition, the contribution of different crystal components to the overall sample crystallinity can be assessed using compositional maps made directly from BSE images or element maps (e.g., Muir et al., 2012), including automated analysis such as QEMSCAN (e.g., Neave et al., 2014, 2017). To illustrate the use of compositional information in BSE images, I use an example from Mount St. Helens, where the diversity of the plagioclase population has been well documented (e.g., Berlo et al., 2007; Cashman and Blundy, 2013).

As illustrated in Figure 9, plagioclase crystals from Mount St. Helens have diverse core populations and mantles that vary in thickness, most likely as a consequence of different durations of pre-eruptive storage in the upper crust. Figure 11A shows the distribution of resorbed (yellow), high-An (turquoise) and low-An (blue) core components as well as mantle zones (green). For measurement simplicity, the crystal populations were separated into only core and mantle components; the former includes cores of variable compositions and textures, the latter includes thin rims that were difficult to isolate at the full thin section scale. Cumulative length (diamonds) and volume (triangles) distributions of cores, rims and whole crystals were calculated assuming size d = √area and volume = d³ (Figure 11B). Resulting CSDs (with N_v calculated as in Eq. 4) are linear; individual core and mantle populations have steeper slopes than measured for the whole crystal population (Figure 11C), consistent with the smaller average sizes of individual core and mantle components. The diversity of the core population prohibits interpretation of core size data as a singe crystallization event. Size information related to the mantles, however, yield information on crystal growth during pre-eruptive magma storage at ∼100–125 MPa (Cashman and Blundy, 2013). Timescales of pre-eruptive magma storage derived from diffusion studies of orthopyroxene rim growth are ≤ ∼2 years (Saunders et al., 2012); if these time scales also apply to plagioclase mantles, then the measured Gt of 0.1 mm yields growth rates ≥1.5 × 10^–9 mm/s. Interestingly, this is similar to rates inferred for crystallization in mafic magmas (e.g., Patwardham and Marsh, 2011; Fornaciai et al., 2015) but faster than rates of 10^–10 to 10^–11 mm/s typically assumed for more silicic melts (e.g., Witter et al., 2016).

Discussion and Conclusion

The examples provided above illustrate both the power and the challenges of using the textures of volcanic rocks to reverse engineer the processes that formed them. Fundamental questions to ask when assessing such samples include:

(1) Does the crystal population of interest represent a single crystallization event or a mixed population?

(2) How can the measured textures be linked to the physical processes responsible for transporting magma to, and across, the Earth’s surface?

Eruption Conditions That Produce Linear CSDs

Examination of both experiments and natural samples shows that cooling-driven crystallization of basaltic lava at the Earth’s surface produces linear CSDs when the erupted magma is at near-liquidus temperatures (few to no pre-existing crystals). Studies of well-constrained lava samples have shown that the kinetics of cooling-driven crystallization are tied closely to the conditions of flow emplacement. Flow emplacement conditions, in turn, are controlled by a balance between rates of advection and rates of cooling and crust formation, coupled with the internal flow dynamics (Griffiths et al., 2003). Time- and temperature-sequential samples may show either nucleation-dominated crystallization, where plots of N_v-φ are positively and linearly correlated (Figure 3B), or growth-dominated crystallization, where CSDs fan around one or more pivot points. Nucleation-dominated crystallization characterizes open-channel basaltic lava flows, where nucleation is enhanced by shearing and melt advection in response to rapid head loss from crust-free marginal shear zones. Crystal growth dominates when cooling occurs after emplacement, for example in well-insulated lava flows and lava lakes. Sequential samples have crystal numbers that often decrease with increasing crystallinity (Figures 2A,3B), a pattern attributed to crystal agglomeration (synneusis), Ostwald ripening and/or size-dependent growth rates.

Decompression of hydrous melts can also produce N_v-φ trends that are positively and linearly correlated and associated linear CSDs; these trends are common for groundmass plagioclase crystals preserved in pyroclasts produced by pulsatory explosions (e.g., Hammer et al., 1999; Cashman and McConnell, 2005; Martel and Poussineau, 2007; Martel, 2012; Figure 7). Co-eruption of crystal-free pyroclasts that require very rapid decompression from depth and crystal-bearing pyroclasts that annealed at shallow levels provide a snapshot of the conduit stratigraphy (Cashman and McConnell, 2005); crystallization rates at different depths (pressures) can be inferred if the time of magma residence in the conduit is known. Textural analysis of these samples provides important constraints on models of magma ascent, arrest, and pressure buildup associated with cyclic eruptive activity (e.g., Diller et al., 2006; Clarke et al., 2007; Burgisser et al., 2011). Inferring time scales from measured CSDs is more difficult because growth rates vary substantially with decompression path (e.g., Brugger and Hammer, 2010a). Importantly, decompression experiments have failed to replicate the very high N_v values observed in some pyroclasts (Figures 4, 6).

The Paradox of Very High N_v

The paradox of very high N_v can be illustrated by the cryptodome that intruded into the Mount St. Helens edifice in the months prior to the 1980 eruption. Here magma intrusion into the edifice was slow and therefore the expectation would be that decompression-driven crystallization would proceed under conditions of low nucleation rates and high growth rates (small Δφ) to produce high crystallinities and relatively low crystal number densities. Instead, cryptodome samples have very high plagioclase number densities and variable plagioclase crystallinities (Figure 8). The very high groundmass number densities therefore raise two questions: (1) what conditions of magma ascent prohibited extensive decompression-driven crystallization during slow magma ascent and (2) what produced the high N_v values of the cryptodome samples?

In the 2 months prior to May 18, the volcano’s north flank was moving outward at ∼2 meters per day (0.08 m/h), which equates to a decompression rate of ∼0.002 MPa/h for an overburden thickness of density of 2500 kg/m³. This intrusion rate estimate overlaps with extrusion rates measured during continuous dome growth in 2004–2008 (0.07–0.36 m/h, or ∼0.002–0.01 MPa/h; Dzurisin, 2018) but is two orders of magnitude slower than the minimum decompression rate used in experiments. Experiments show that crystallization efficiency is reduced when decompression is slow, particularly when growth is on pre-existing crystals. This raises the interesting question of whether very slow ascent might allow magma to rise to shallow levels without (significant) groundmass crystallization, perhaps explaining the stealth nature of cryptodome emplacement.

Reconstruction of the cryptodome geometry (Donnadieu and Merle, 1998) suggests that the shallowest part of the cryptodome occupied a pressure range of ∼7–18 MPa consistent with a maximum pressure of ∼20 MPa suggested by measured H₂O of 0.07–1.5 wt% in cryptodome glass (Hoblitt and Harmon, 1993; Newman and Lowenstern, 2002; Neill et al., 2010). Phreatomagmatic eruptions in the months prior to the 18 May eruption ejected material from the cryptodome (Cashman and Hoblitt, 2004); there is also isotopic evidence of open system degassing from the cryptodome during this time period (Neill et al., 2010). This raises the question of whether rapid decompression of shallow-stored magma (during a precursory eruption) followed by annealing during an inter-eruption repose period could generate the characteristic high N_v of the cryptodome samples. Similarly, the high N_v pyroclasts at Pinatubo occur in magma that stalled at pressures of 8–16 MPa between precursory eruptions (Hammer et al., 1999); here the low crystallinities can be explained by repose intervals that were much shorter than at Mount St. Helens (hours instead of days to weeks). Importantly, most decompression experiments have used initial pressures P_i ≥ 125 MPa. This raises the questions of whether very high crystal number densities could be replicated experimentally by single step decompression from low pressures followed by annealing times of hours to weeks.

Re-thinking CSD Analysis of Volcanic Material

There is abundant evidence that most volcanic materials contain a complex mixture of crystals from multiple sources. Experiments show that the presence of pre-existing crystals creates CSDs that are either kinked or curved because new groundmass crystals are added as rims grow on pre-existing crystals (Figure 5). Kinked or curved CSDs measured in volcanic samples, however, are typically treated as two populations – groundmass and phenocrysts - that experienced two separate (and unrelated) crystallization histories (e.g., Higgins, 1996; Neave et al., 2013; Witter et al., 2016). There are several problems with this interpretation. First, heterogeneous crystal cores record crystallization over different time scales (e.g., Cooper and Kent, 2014), and under different P_H2O-T-X conditions (e.g., Berlo et al., 2007; Humphreys et al., 2009), and therefore do not represent a single and unique crystallization episode. Second, most crystals also have mantles (often oscillatory zoned) that record growth in pre-eruptive upper crustal magma storage regions. Variations in the thickness and composition of these mantles requires individual crystals erupted together to have spent different amounts of time, and even in different parts of, magma reservoirs prior to eruption (Cashman and Blundy, 2013). Finally, crystal rims form during magma ascent (decompression) and therefore represent the same crystallization event responsible for nucleation and growth of groundmass crystals. Treating large crystals as a single crystallizing population is therefore misleading. This problem is compounded when resulting CSD data are used to infer a crystal residence time, which requires assumption of a single average growth rate for the entire phenocryst population.

To move forward, it is important to examine the goals of textural analysis. The reason to study crystallization during eruptive episodes is clear: groundmass crystallization affects magma rheology and permeability and thus understanding syn-eruptive crystallization histories provides critical constraints for estimating magma rheology (e.g., Klein et al., 2018), modeling lava flow advance (e.g., Dietterich et al., 2017) and anticipating transitions in eruptive styles (e.g., Cassidy et al., 2018). To examine the evolution of magmatic systems and assembly of eruptive magma bodies (e.g., Flaherty et al., 2018) further requires not only identification of different crystal populations (e.g., Morgan et al., 2007; Kahl et al., 2013; Wieser et al., 2020) and sources of those populations (e.g., Morgan et al., 2007) but also the extent and conditions of their shared ascent histories (e.g., Cashman and Blundy, 2013). Measuring the amount of crystallization that can be attributed to each population would allow textural analysis to be linked, for example, to diffusion chronometry (e.g., Jerram et al., 2018) and to models of magmatic systems (e.g., Blundy and Cashman, 2008). Unraveling the crystallization histories of diverse crystal populations in thin section requires, ideally, large area compositional mapping at high resolution (e.g., QEMSCAN) coupled with detailed isotopic measurements and diffusion chronometry of individual crystals to identify the origin of, and time of entrainment, different components of each crystal population. Advances in tomographic imaging could further allow compositional textural analysis in 3D. The primary challenges of this approach are balancing spatial coverage with resolution, as well as handling the computational expense.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Funding

KC acknowledges the AXA Research Fund (Volcanology Chair) and a Royal Society Wolfson Research Merit Award for supporting this research.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

I thank reviewer DJ for a very helpful review, Jenny Riker for help with interpreting the experimental literature (although mis-interpretations are all mine!) and Jon Blundy for years of collaborative work on the Mount St. Helens magmatic system. I also acknowledge the AXA Research Fund and a Royal Society Wolfson Research Merit Award for support for this work.

References