Many plagioclase phenocrysts from volcanic and plutonic rocks display quite complex chemical and textural zoning patterns. Understanding the zoning patterns and variety of crystal populations holds clues to the processes and timescales that lead to the formation of the igneous rocks. However, in addition to a “true” natural complexity of the crystal population, the large variety of plagioclase types can be partly artifacts of the use of two-dimensional (2D) petrographic thin sections and random cuts of three-dimensional (3D) plagioclase crystals. Thus, the identification of the true number of plagioclase populations, and the decision of which are “representative” crystal sections to be used for detailed trace element and isotope analysis is not obvious and tends to be subjective.
Here we approach this problem with a series of numerical simulations and statistical analyses of a variety of plagioclase crystals zoned in 3D. We analyze the effect of increasing complexity of zoning based on 2D chemical maps (e.g., backscattered electron images, BSE). We first analyze the random sections of single crystals, and then study the effect of mixing of different crystal populations in the samples. By quantifying the similarity of the compositional histogram of about a hundred 2D plagioclase sections it is possible to identify the so-called reference and ideal sections that are representative of the real 3D crystal populations. These section types allow filtering out the random-cut effects and explain more than 90% of the plagioclase compositional data of a given sample. Our method allows the identification of the main crystal populations and representative crystals that can then be used for a more robust interpretation of magmatic processes and timescales.
The compositional and textural features of plagioclase from igneous rocks have been investigated for a long time, first using its optical properties and the petrographic microscope (e.g., Hibbard 1995; Shelley 1993) and more recently using a combined approach of scanning electron microscope (e.g., Ginibre et al. 2002), electron and ion microprobes (e.g., Blundy and Shimizu 1991; Singer et al. 1995), rim-to-core analysis (e.g., Bouvet De Maisonneuve et al. 2012; Neill et al. 2015), cathodoluminescence studies (e.g., Higgins et al. 2015), and in situ microdrilling for isotopes (Davidson et al. 2001). Feldspar studies have provided critical clues about the magmatic processes (e.g., magma mixing, assimilation, fractionation, magma ascent, crystal recycling; e.g., Anderson 1984; Feeley and Dungan 1996; Landi et al. 2004; Streck 2008) and their associated timescales (e.g., Costa et al. 2003; Druitt et al. 2012; Stelten et al. 2015; Zellmer et al. 1999). However, a common observation of all of these studies is the extreme variety and the chemical and textural complexity of the plagioclase crystals, including those of open degassing volcanoes such as Stromboli, Etna, Llaima, or Mayon (Bouvet De Maisonneuve et al. 2012; Landi et al. 2004; Nicotra and Viccaro 2012), e.g., the example from Mayon shown in Figure 1 that we will discuss in detail below. The variety of plagioclase phenocrysts makes it very difficult to establish how many crystal populations are present in the deposit, and to objectively decide whether there are any phenocryst section that can be considered as “representative” of a population and thus used to derive conclusions about the processes from electron or ion microprobe data (e.g., Singer et al. 1995).
The variety and complexity of plagioclase textures and zoning can reflect the large number of processes and variables that affect its stability (e.g., Yoder et al. 1957), but may also partly be the result that we typically study them using 2D petrographic thin sections that are derived from a set of randomly cut 3D crystals. The 3D to 2D conversion of the crystals can significantly affect their shape (e.g., Higgins 1994) but also produce an artificial variety of zoning patterns and increase the apparent number of crystal populations (Pearce 1984; Wallace and Bergantz 2004). Wallace and Bergantz (2002) proposed a methodology to correlate between different crystals using one-dimensional (1D) traverses by doing a wavelet analysis of the anorthite (An = 100·Ca/[Ca+Na]) content. They designed a new crystal phylogeny analysis that showed how different crystals in the same deposits could be related, and shared a larger portion of their history going from core to rim. Later Wallace and Bergantz (2005) recognized the effects of random cuts and proposed a methodology to correct for them, although also acknowledged that loss of information by random cuts missing the inner parts of the crystals is unavoidable.
In this paper, we take a complementary approach by analyzing the 2D compositional features of plagioclase based on BSE images (e.g., An histograms of 2D plagioclase sections; e.g., Cashman and Blundy 2013). Given the ease with which the BSE images of a large number of crystals can be gathered, it is possible to obtain a large statistical data set that can be used to characterize the deposits and the crystals. Numerical simulation of 2D crystal sections in thin section may allow understanding the 3D crystals and magma processes in a similar manner to the quantification of crystal size distribution studies (Higgins 1994; Morgan and Jerram 2006), olivine and pyroxene zoning patterns, and timescales (Pearce 1984; Shea et al. 2015). However, such analysis for complexly zoned crystals like plagioclase is basically not available.
We first present numerical simulations of 3D plagioclase zoned crystals, which we use to quantify the effect of the 3D to 2D conversion using compositional maps. We show how it is possible to define some special sections that are fully representative of the 3D crystals by quantifying the similarity between the different 2D sections. We apply these numerical models to a variety of zoning patterns and mixed crystal populations and demonstrate how our approach can retrieve the original true populations and explain >90% of the compositional data for the samples.
Model setup and strategy for characterizing and comparing compositional zoning
When studying natural samples we do not know a priori how many crystal populations there are in a given sample, and it is difficult to recognize the artifacts of random sampling, in particular if crystals are geometrically and compositionally complex. Our approach is to first construct numerical crystals and perform forward models of their zoning patterns to find measurable variables and statistical procedures to identify the 2D sections that can be confidently taken as representative of the 3D crystal populations. We first study various relationships between populations using simple crystals, which we then make progressively more complex. Later we use these findings to study mixed crystal populations, first in the numerical experiments and then in natural samples.
Strategy: 1D vs. 2D data
The methods proposed in the literature to group crystals, such as wavelet-based correlations (Wallace and Bergantz 2002) and shared characteristic diagrams (Wallace and Bergantz 2005) are based on comparison of 1D profiles. However, given the complexity seen in 2D BSE crystal images, it is not straightforward to decide which 1D profile is representative of a given 2D section. Moreover, identification of crystal populations becomes increasingly difficult with increasing degrees of geometric complexity. Thus, we have used another way to classify plagioclase population based on the area frequency compositional distribution of 2D plagioclase sections (Cashman and Blundy 2013). The advantage is that we have a better overall characterization of the compositional data of the 2D section, but it has the disadvantage that we lose the spatial information, in particular the core-rim relationship. This hampers a detailed understanding of the processes that created the zoning patterns, but this is not the goal of our contribution. We aim at first identifying the crystal populations, and later we can turn to detailed core to rim traverses only for representative sections. Otherwise, given the textural and compositional variability of many natural plagioclase crystals it is impossible to unambiguously choose the crystal sections that are representative to do detailed 1D compositional traverses. A crystal population is made of crystals that have experienced very similar magmatic processes, although in the context of our analysis a crystal population is defined by those crystals that have them same 2D An (or grayscale) distribution (within error; see below). Moreover, the fact that the 2D sections are the result of random cuts already provides a certain degree of information about the core to rim zoning, because for example, the cores can be under sampled, but the rims are not. A large number of compositional information from 2D sections of crystals can be obtained from BSE images calibrated for their grayscale values using a few quantitative analyses and the electron microprobe (Ginibre et al. 2002). It is nowadays possible to produce BSE images of a full thin section using SEM under the same analytical conditions so that all plagioclase crystals can be compared (e.g., Fig. 1).
Numerical 3D simulation of plagioclase crystals and generation of random 2D sections
We first constructed a numerical model of a 3D plagioclase crystal with the geometry and number of faces according to theory (e.g., Deer et al. 1992; Higgins 1991, 2006; Fig. 2) and with five compositional zones using Matlab 2014b software environment (Mathworks 2014). Note that although plagioclase belongs to the triclinic system (Deer et al. 1992) we used a reference frame of three perpendicular axes that are parallel to the zoning patterns for our numerical simulations. We also performed a numerical simulation of a crystal with an angle of 115° between X and Y axes (close to the theoretical value for anorthite according to Deer et al. 1992) and we later show that this does not affect our results. We varied the ratios between the three dimensions of the crystal (S = shortest dimension, I = intermediate dimension, and L = longest dimension; (Deer et al. 1992; Higgins 1991, 2006) (Fig. 2). Previous studies showed that the S:I:L of most plagioclase phenocrysts may range from 1:2.8:4 and 1:5:8 to 1:5:5 (e.g., Cheng et al. 2014; Higgins 1994; Morgan and Jerram 2006) and we tested an extreme range of shapes from 1:1:2 and 1:2:5 to 1:5:5. To simplify the comparisons, we normalized the lengths using an arbitrary “unit.” The longest dimension of the plagioclase was fixed to 200 units along the X-axis; then, the intermediate and short dimensions were computed according to the different shapes along the Y and Z axes. For example, if the model uses the shape 1:1:2, then the intermediate and shortest axes (Y and Z) were both set to 100 units (Fig. 2).
Characterization of 2D plagioclase zoning patterns and compositional maps
The random cuts produce a large variety of zoning patterns (Fig. 3), ranging from compositionally homogenous sections (e.g., section 4 in Fig. 3) to up to a five-zoned pattern (e.g., section 9 in Fig. 3). To characterize the 2D sections we used the normalized area compositional histograms of each section. This means that for each section we calculated the area with the same composition (i.e., the same An content within a given tolerance, see below) and normalized it to the total area of the crystal section (Fig. 4). We call the random individual sections RI where I = 1, 2, 3…Z and Z is the total number of sections. To characterize the compositional distribution of all 2D random sections from a 3D crystal we added all the areas of the zones with the same composition from all sections and normalized them to the total area of the sections. We call this distribution the “all random population,” abbreviated as RA. This is an important distribution because it can be used to characterize the various crystal populations in natural samples for which we do not know a priori the zoning of the 3D crystal. In the same way, we can calculate the area compositional distribution for the ID sections (Fig. 4). We found that the compositional distributions of the ID and all RA for the simple five zoned crystal are similar but not identical; typically the compositions of the rims are over-represented and those of the crystal centers under-represented in the RA distribution (Fig. 5). This reflects that although all random sections pass through the crystal rims, some sections miss completely the cores (e.g., off-core sections). This gives a similarity of about 75% between the RA and ID distributions (see below for the precise definition of similarity used) and this difference varies depending on the details of the crystal zoning.
Comparing between 2D sections using the similarity of compositional maps
For example, the similarity between the RA and the ID distributions (S_ID*RA) for the simple crystal is about 75% (Fig. 6). The similarity depends on the number of bins, which in our case depends in turn on the An range we choose for each bin. Many authors advise that for real data sets histograms based on 5–20 bins usually suffice, noted by Scott (1979). Calibration of grayscale of BSE images of plagioclase typically produces An values with a precision of 1 to 2 mol% (Ginibre et al. 2002). We tested bin sizes of 1, 2, and 3 mol% An and found that the 3 mol% (the number of bins is more than 30) is quite similar to using 1 mol% (see Supplemental1 Appendix 1), so we used a bin size of 3 mol% An unless otherwise noted. We will use the similarity between different sections or groups of sections in this manner for the rest of the manuscript (e.g., Fig. 6). The uncertainty of the similarity was calculated to be up to 1% by assuming a relative analytical uncertainty of 0.5% on the An content in this study.
Reference and ideal sections
Using the method we described above, we can calculate the similarity between any section types and/or between groups of sections. This is the first step to be able to identify ideal or representative 2D sections and thus different crystal populations in natural samples. For example, the compositional similarity between the ID and RA distributions can be calculated to be about 75% (Fig. 6). The simplest approach to determine the number of crystal populations and representative sections would be to use directly the ID and RA distribution, but this is not possible for several reasons. The RA distribution actually does not correspond to any individual crystal section, so although it characterizes the overall plagioclase population we still need to identify representative sections. Moreover, although we know what are the ID sections in our numerical experiments, this is not the case in natural samples. These problems appear when we want to filter out the random cuts for a single-crystal population, but they become even more apparent when there are multiple crystal populations rather than a single one. Thus, we have designed a strategy where we first try to filter out the effects of random cuts using proxies to characterize the RA distribution and the number of populations, and then we can identify the best 2D sections using more strict criteria and the ID sections.
The first step is to calculate the similarity distribution between the overall random population (RA) and each random section (R1, R2, R3…RZ) (named S_RA*R1-RZ; Fig. 7). We find that it shows one main peak about 70–80% similarity that includes 25% of the sections (a total of 180 sections out of 745), and another at 30% similarity that includes about 14% of the sections. Moreover, more than 50% of the random sections have a similarity (S_RA*R1-RZ) ≥70% with that of the overall population, and thus they could be considered as a reasonably representative sample of the 3D crystal. A much smaller number of sections (about 5%) have a similarity ≥90% with the overall random population and we call them the reference sections (REF). Thus, we can use the RA distribution to identify individual reference sections from the overall population. However, this method works if there is only a single 3D crystal population, but not if there are more, since then the similarity values are much lower and it is not possible to identify the REF sections although they must exist in data sample (see detail in the population mixing sections).
Another way of identifying the reference sections is by calculating the similarity between a given individual section and all the rest of the random sections (e.g., not using the RA directly; S_R1*R1-RZ; Fig. 8). The similarity distribution pattern of each individual section and the rest is quite variable, but it is apparent that reference sections are more abundant for similarities ≥70% (Fig. 8a). This finding is akin to the effect of random cuts of crystal shapes, where some characteristic and important shapes, areas, and dimensions are more frequent because of the geometrical symmetry of the object (e.g., Higgins 1994). In our case, it is apparent that even if the cuts were done at random, the compositional histograms of sections parallel or perpendicular to the main geometry axes have the same compositional distribution and thus are more frequent to begin with than any other section (Fig. 8b). See Supplemental1 Appendix 2 for more details about the method involved in identification of these sections. We also discuss later how this method can be improved by using thresholds of similarities.
Another important group of sections that need to be characterized are the IDs, and we have calculated the similarity distribution between ID and each individual section (RI) (here called S_ID*R1-RZ; Fig. 7). The distribution has a peak at about 20% similarity, which means that most random sections are quite different from the ID, which reflects the fact that the random 2D cut effect significantly changes the compositional maps, with many sections that are off-center. However, there are more sections with a similarity (S_ID*R1-RZ; Fig. 7) ≥90% to the ID than to the overall random distribution (S_RA*R1-RZ; Fig. 7). Thus, although the random cut effects are important, there are still many 2D sections that are very similar (S_ID*R1-RZ ≥ 90%) to the ideal sections. These random sections record the most complete information and we will also call them ideal sections because they are ≥90% similar to the theoretical ID sections.
The REF are different from the ideal sections in that they typically do not record the inner parts of crystals. However, since we do not know a priori how many crystal populations there are in a given natural data set, we have to use the method of calculating the similarity between each random section (Fig. 8). It is apparent that the most frequent sections for similarities ≥90% are the ideal ones (see Supplemental1 Appendix 2 for more details). Later we quantify better these relationships using threshold values; here we wanted to illustrate that it is possible to identify the ID and REF by comparing each individual section to the rest in a systematic and statistical manner. In typical studies of plagioclase zoning, petrologists tend to use sections with shapes close to the most classical shapes of plagioclase, i.e., rectangular or square. If we do a subsample of the rectangular and square sections (sides ≥ 5) from our random samples we find that about 36% of those have a similarity with the ID of 90% or higher, which means that, in the case of a single-crystal population, the likelihood of choosing an ID section if the shape is rectangular is much higher.
So far, we have used a crystal with three perpendicular axes, but plagioclase belongs to the triclinic system and thus one axis at an angle that deviates significantly from 90° from the rest (Deer et al. 1992). To test whether this has effect on our method we have also done a simulation with a crystal that has different angles between the axes. Figure 9 shows the wire-structure of a zoned crystal with an angle of 115° between two axes that is close to triclinic system of anorthite plagioclase (according to Deer et al. 1992). In total, we produce 900 cuts to get 714 2D patterns. We still could find that it shows one main peak about 70–80% similarity that includes 25% of the sections, and another at 35% similarity that includes about 20% of the sections. The ID distribution has a peak at about 20% similarity, and there are more sections with a similarity ≥90% (Fig. 9). The two distributions are similarity to that of three perpendicular axes model. Thus, the approach described above for the identification of the different ID and REF sections is still valid. This is because the effects of random cuts are much larger than the small deviations deriving from the “incorrect” use of three perpendicular axes.
Models with complex crystal zoning patterns
Natural plagioclase crystals are however often more complex that the five-zone crystal simulation (Fig. 2) we have used so far (e.g., Ginibre et al. 2002). Thus, we also generated three other types of 3D crystals (Fig. 10): Type A with patterns of An values that generally increase from rim to core with oscillatory zoning and fine variations of less than 5% An, Type B with large and abrupt An changes between core and rim, and Type C that combines the zoning characteristics of the previous two types. Calculation of the similarity properties for the three types of crystals shows that they can be treated in the same manner as the simpler crystal, although in detail the similarity thresholds for the identification of the reference and ideal sections are somewhat different (Fig. 10). The similarity histogram of RA and each individual random section has main peaks at 70–90% similarity for crystals of Type A, at 70–80% similarity for those of Type B, and at 80–90% similarity for those of Type C, whereas the similarities of ID and each random section have also high frequency for all crystal zoning types at similarity ≥90% (Fig. 10).
To further refine our approach to find the reference and ideal sections we need to do an additional step that is critical when we are dealing with multiple crystal populations (e.g., for cases of magma mingling). We first calculate the similarity between each random section, and then we calculate the similarity between the section that has the largest number of section pairs with a minimum similarity threshold, starting from 50 up to 100%, with that of the RA distribution (please see Supplemental1 Appendix 2 for a step by step description of this approach). We determine where this similarity is at a maximum (Fig. 11). By doing this we find that for example, for Type A crystals we have a similarity ≥90% between some random sections and the RA for a similarity threshold between 70–80% (Fig. 11a). This means that these random sections can be taken as representatives of the Reference sections and thus used to characterize one crystal population. Note that this threshold is also higher than the mean threshold, which provides additional constraints to identify the appropriate threshold value. Different zoning types give slightly different thresholds, but at about 80% threshold, the sections are all ≥90% similar to the RA and thus can be used to identify the reference section. Similar relations can be used to identify the ideal sections (Fig. 11b) but these require higher thresholds at 90–100%. Finally, we also did a Monte Carlo simulation where we tested up to 1200 3D crystals with different An zoning patterns. We found that our inferences of similarity thresholds at 70–80% and 90–100% are robust and able to identify the Reference and Ideal sections, respectively, independent of the type of zoning (Fig. 12).
Multiple crystal populations
The next level of complexity is to be able to identify the reference and ideal sections for mixed crystal populations. We designed a numerical experiment were we first generated 50:50 mixtures of two crystal populations involving three different scenarios (Fig. 13):
Two crystal populations with different An core compositions but the same rim compositions, as representative of mafic-silicic magma interactions (e.g., Feeley and Dungan 1996); the similarity between these two populations is relatively low, about 15% (Figs. 13a, 13d, and 13g).
We find that the similarity histograms for the three mixing scenarios (Figs. 13g, 13h, and 13i) are quite different from those of single-crystal populations only affected by the random cut effect (e.g., Fig. 10). Most notably the maximum of similarity tends to reach much lower values (Figs. 13g, 13h, and 13i). In these examples it is apparent that the similarity information of the 2D crystal sections and the RA reflects the combination of both the cut effect and the mixing of the two populations. In effect, the overall RA that we might measure in a natural sample is a “weighted average” that is the result of mixing the RA of each individual population (Figs. 13g, 13h, and 13i).
To identify the two crystal populations we calculated the similarities between the different sections and we focused only on those with a similarity threshold around 80% (as suggested by Fig. 11). We then calculated the similarity of these sections with the overall RA for each of the scenarios and found low similarities ranging from 46 to 77% (Figs. 14a, 14c, and 14e). This means that there has to be other reference sections in the population that would explain the full data set. Thus, we removed from the pool the sections with a given similarity threshold (e.g., >80%) to the first reference section, and chose another section with a threshold around 80% and calculated whether the similarity with the overall population increased by mixing them in different proportions. We keep doing this using least-squares minimization between the RA distributions and that of mixing different sections in different proportions. Once we found that the overall similarity increased to >90% (Figs. 14a, 14c, and 14e) we stopped, and we considered that we were able to explain >90% of the overall compositional distribution of the plagioclase.
For the three scenarios above we were able to obtain similarities of the two reference sections and the RA of each scenario between 92 and 95% (Figs. 14b, 14d, and 14f). This approach worked for the reference sections as well as for the ideal sections. We also tested the effect of varying proportions of mixing ratios between two crystal populations (e.g., 10:90, 70:30) using scenario two (Fig. 15). We found that using the same approach of reference/ideal sections and least-squares minimization we are still able to recover the two populations. However, when the proportion of one population gets close to 90% or more it becomes more difficult to identify any other, until it becomes eventually undetectable. An example of application of our approach to mixed crystal populations from a lava of Mayon volcanoes can be found in Supplemental1 Appendix 3.
How many crystals are needed to represent a given sample
A related important question when studying natural crystals is how many random sections are needed to characterize the population(s) and build a reliable RA distribution. We used the Type C plagioclase (see Fig. 10) as the sample test. We made 400 random cuts and calculated their RA (RA400). We then compared their similarity to an increasing amount of random sections, randomly extracting 10 of them (RA10) until 300 (RA300). Because of the large variability of the RA when we sample a small number of sections at random (e.g., RA10; Fig. 16a), we calculated their RA multiple times, and thus their similarity also varies (Fig. 14b). We found that within about 60–100 random picks of random sections we characterize the full variability (Fig. 16b). Moreover, we find that the similarity between the distributions obtained from 80 (RA80) to 100 (RA100) random sections with that of the 400 sections (RA400) is consistently higher than 95%, and thus we conservatively suggest that 100 random 2D sections are enough to represent RA of the entire crystal population.
Implications for identification of plagioclase populations and representative crystal sections
The large variety of compositional and textural features of many plagioclase crystals from igneous rocks can be party explained by the effects of 2D random cuts of 3D plagioclase crystals. Using the compositional histogram of 2D sections of many plagioclase crystals and statistical analyses based on the compositional similarity between the different sections it is possible to account for the effects of the random cuts and separate them from the effects of the presence of different crystal populations. With our method, we can identify the different crystal populations that are at least 90% similar and the proportions of the different populations. The identification of the different populations by using reference and ideal sections removes the problem of a subjective choice of sections and allows studying the crystal sections that are more representative of the samples. These can be further studied using detailed electron microprobe, ion microprobe, and/or microdrilling of isotopes, which should lead to a much more robust understanding of magmatic processes based on plagioclase crystal records.
We acknowledge discussion with D. Ruth on this topic, and C. Newhall for the 1947 sample of Mayon. G. Fabbro, J. Herrin, D. Ruth, Sri Budhi, W. Li, D. Schonwalder, C. Widijiwayanti, L. Zeng, T. Flaherty, O. Bergal, and D. Nurfani are thanked for their participation in the expert elicitation exercise on the number of plagioclase types. Reviews by G. Bergantz, P. Izbekov, M. Higgins, and editorial handling and reviews of J. Hammer significantly improved the presentation of ideas in this manuscript and greatly appreciated. This research is supported by the National Research Foundation of Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative, under the “Crystal pattern” project.