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

Studies of magmatic systems have long used the textures of erupted samples to infer processes that control the location and duration of magma storage and drive volcanic eruptions from these storage regions. Models of volcanic processes and magmatic systems have evolved substantially over the past decades, in large part because of advances in analytical and experimental techniques. Cooling- and decompression-experiments have greatly enhanced our understanding of crystal textures produced by crystallization associated with volcanic eruptions, while advances in compositional mapping, isotopic analysis and diffusion chronometry provide the tools to unravel complex histories of individual crystals. Experiments, however, have failed to replicate the full range of groundmass textures observed in volcanic samples and the recognition that magma commonly includes both indigenous (grown from the transporting liquid) and exogenous (incorporated from elsewhere in the system) crystals complicates interpretation of crystal populations in volcanic samples. Analysis and interpretation of crystal size distributions (CSDs) and other physical measures of crystal populations, in particular, have yet to fully account for crystal populations with diverse origins and growth histories. Here I assess the extent to which experiments replicate observed crystal populations and thus can be used to improve understanding of volcanic processes. I then review conditions under which the size characteristics of crystal populations can be reasonably interpreted, examine possible reasons for experimental failure to achieve the very high crystal number densities that characterize some eruptive samples, and suggest ways to link CSD analysis to other techniques that seek to constrain the origin of the complex crystal populations. Finally, I show that compositionally based crystal size measurements are critical for interpreting different stages of crystal growth and can be yield well constrained growth histories if linked to diffusion time scales and phase constraints on crystallization conditions.


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., 2005Davidson et al., , 2007Humphreys et al., 2006Humphreys et al., , 2008Smith 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., 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., 1999Hammer et al., , 2000Cashman 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 (superor 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 2 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 timedependent, 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 FIGURE 1 | Experimental constraints on phase relations and kinetics. (A) Equilibrium crystallinity as a function of pressure. Data from Riker et al. (2015a;R2015), Couch et al. (2003;C2003), Martel (2012M2012), and Brugger and Hammer (2010b;BH2010). (B) Average ratio of nucleation rate (J) to growth rate (G); data from experiments of Hammer and Rutherford (2002). (C) Crystal size as a function of time at P f = 75MPa. Dashed lines show growth rates between different experiment times; dotted line shows the time-averaged, or apparent, growth for the entire duration of the experiment.
crystal nucleation and growth that can be parameterized to link crystal size (L) and number density (N v , per volume) to timeaveraged 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 2 , 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 2 ; 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(Hammer et al., , 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 2 ), 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 3 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 . 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  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 Frontiers in Earth Science | www.frontiersin.org (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 3 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 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 2 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) singlestep 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 2 to 10 7 /mm 3 ; 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 2 to > 10 8 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 andRiker et al., 2015b;Pinatubo -Hammer et al., 1999 andAndrews, 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 welldocumented 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 2 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),   Hammer et al. (1999). (B) MSH cryptodome. 541 samples were erupted early, and 505 samples were erupted late; D and dense and V are vesicular. Samples described in 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 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   Cashman and McConnell, 2005). BSE images show compositional variations; brighter zones are An-rich. (A) Crystal with resorbed An-rich core and thick oscillatory zoned rim. (B) Crystal with extensive patchy resorption in the core and overgrowth rim with two zones; note the agglomerated oscillatory zoned crystal enclosed within the outer rim. (C) Highly resorbed and anhedral An-rich crystal with no rim growth. (D) Resorbed Ab-rich core with thick normally zoned rim; note late-adhering crystal with its own late-stage rim. 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 2 ) than for the microphenocryst (710 µm 2 ); 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., 2004Kahl 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., 2014Neave et al., , 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 3 ( 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 temperaturesequential 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. Nucleationdominated 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 3 . This intrusion rate estimate overlaps with extrusion rates measured during continuous dome growth in -2008Dzurisin, 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 2 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 intereruption 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 populationsgroundmass 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 syneruptive 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.