A multivariate statistical approach for mineral geographic provenance determination using laser-induced breakdown spectroscopy and electron microprobe chemical data: A case study of copper-bearing tourmalines

LIBS analyses

LIBS is a relatively recent analytical technique that is finding utility in the geosciences [e.g., see reviews by Fabre (2020) and Harmon and Senesi (2021)]. The information-rich spectra contain signatures of all elements in concentrations above detection limits (e.g., Cremers and Radziemski 2013), molecular emissions, select isotopic ratios (e.g., Smith et al. 2002; Doucet et al. 2011: Russo et al. 2011), and some structural information (Serrano et al. 2015) resulting in a detailed chemical fingerprint of the material analyzed. To take advantage of the rich chemical data set embedded in tourmaline, this LIBS study uses the spectrum of relative peak intensities of each tourmaline rather than absolute quantities of individual elements within the tourmaline.

Minimal sample preparation is required for LIBS [see, e.g., McMillan et al. (2018, 2019) for additional information]. Rough samples were cleaned with isopropyl alcohol to remove oils and surface residue and air-dried. Most tourmalines are individual grains or clusters of grains. Originally, samples were mounted on a plexiglass sheet with BlueTac to secure the grains; later, the BlueTac was eliminated. The sheet was placed into the sample holder in the LIBS instrument chamber. LIBS data were acquired prior to EMP data analyses to avoid any possible contamination from EMP sample preparation, such as polishing and carbon coating of the grains.

Tourmalines were analyzed with an Applied Spectra J200 LIBS instrument at Materialytics, Inc., fitted with a Q-switched Quantel ULTRA 100 Big Sky Nd:YAG laser operated at a fundamental wavelength of 266 nm and <6 ns pulse width. The instrument utilized an Andor Mechelle ME 5000 spectrograph (λ/Δλ = 5000) and an Andor iStar ICCD (intensified charge-coupled device) camera, model DH334T-18F-03. Analytical conditions were a laser power of 150 mJ, with a delay of 0.5 microseconds (μs) between the time of the laser shot and light collection, a gate width (time of light collection) of 10 μs, and a nominal spot size of 50 μm (subsequent analyses demonstrated a larger ablation pit of nearly 80 μm). Spectra were obtained at 1 atm at room temperature in an argon atmosphere to confine the LIBS plasma and thus enhance emission intensity. Where grain size allowed, 64 shots were obtained per sample in an 8 × 8 grid with a spacing of 100 μm between shots—an area covering about 1 × 1 mm. An ancillary study suggested that 64 shots were optimal for characterizing the samples (McMillan et al. 2019). At each analytical location, a cleaning shot was done prior to the analytical shot. The spectral emission was collected over the 26 000+ channels of the detector/spectrometer system to assemble the spectrum in the wavelength range from 200–1000 nm for each analytical shot. Spectra were truncated at 771 nm which preserves the potassium peaks at 766.5 and 769.9 nm but masks the primary argon peaks at higher wavelengths. Multiple shots per sample and their corresponding spectra are averaged and normalized to the mean peak intensity to produce a single spectrum per sample. Averaging LIBS spectra helps mitigate variations caused by inherent shot-to-shot variability (McMillan and Dutrow 2024). Background correction was not applied. Intensities were converted to log values for modeling purposes. Where necessary, identification of LIBS peak positions utilized the online NIST database of optical emission lines (Kramida et al. 2022).

Acquisition of such a large data set requires statistical methods and/or machine-learning techniques for data analyses and interpretation. This study employs the multivariate statistical techniques principal component analysis (PCA) (Esbensen 2004) and partial least-squares regression (PLSR) (Wold et al. 2001; Esbensen 2004) to quantitatively classify spectra with reference to the geographical source of the tourmaline. The strong emission response of some major elements required the masking of select peaks from the spectra to allow subtler chemical variations to be enhanced. For these tourmalines, masking of peaks for the elements Si, Al, Li, Na as well as the Ca peaks at 393.3, 396.8, and 422.7 nm resulted in improved models. While other multivariate statistical techniques may be advantageous, for this test case, methods used previously were followed (e.g., McMillan et al. 2018).

Multivariate statistical modeling

PCA is a dimension-reducing multivariate technique that calculates linear regressions, or Principal Components (PCs), through the data set in multivariate space (24 350 variables). A PCA score plot (sample analyses in n-dimensional space projected onto the plane of two principal components, e.g., PC 1 and PC 2) displays the spectral/compositional relationships of the data set in the two directions of the principal components. This comparison is used to determine the order in which the geographic localities (SJdB, RGdN, Moz, Nig) are modeled beginning with the compositionally most distinct group modeled first (Multari et al. 2010; Kochelek et al. 2015; McMillan et al. 2018).

PLSR models were used to quantitatively discriminate between the samples of the locality of interest and all other localities. PLSR is similar to PCA but includes the value of an independent variable, in this case, the Provenance Variable (PV), in the regression. Spectra of samples from the locality of interest were assigned a Provenance Variable value of 1; spectra of samples from all other localities were assigned a PV value of 0. To calibrate the model, 50% of the spectra from the geographic localities were selected; spectra from the 50% remaining samples were used for test-set validation in a later step. Because the database contained one spectrum per sample, no individual sample was present in both the calibration and validation sets, although samples from a given geographic locality were present in both sets. Statistical modeling was accomplished using the Unscrambler software by Camo. The nonlinear iterative partial least squares (NIPALS) algorithm was applied with 15 PLSR components; no weighting was applied to variables. All models are mean-centered [see also McMillan et al. (2018) for further discussion].

To quantitatively assign a spectrum to a locality group, a numerical value that separates calculated Provenance Variable values for the two groups in the calibration set is defined: the value of apparent distinction (VAD) (Kochelek et al. 2015). The VAD is calculated as the value that gives the highest number of correctly assigned samples during calibration. Any sample with a calculated PV value greater than or equal to the VAD is classified as a tourmaline within the group of interest; those with calculated locality variables less than the VAD are classified as belonging to the group of the remaining localities. Once a VAD is assigned, it does not change during validation.

PLSR models were validated using test-set validation. PV values are calculated for tourmaline spectra not used to calibrate the PLSR model. The VAD determined during calibration is used to predict whether each spectrum in the validation set belongs to the locality of interest or the group of the remaining localities. The prediction accuracy is calculated as the percent of correctly assigned test-set spectra for which locality information is known. For example, Model 1 evaluates São José da Batalha (SJdB) samples. Applying the VAD of 0.45 to the spectra not used in the calibration set, all of the São José da Batalha samples are predicted to be from this locality, as well as one African sample, and the other samples are predicted to belong to the group of remaining samples (Fig. 3). Thus, Model 1 is 96% successful, one sample is miscategorized. Once validated, the decision tree of PLSR plots is developed for each remaining group of samples (RGdN, Moz, Nig).

Figure 3

Graphs showing the calibration (left column) and validation (right column) results for the LIBS decision tree based on geographic locality. Samples used in the calibration set were not used in the validation set, leading to a different data distribution. The value of apparent distinction (VAD) determined in the calibration, shown as a dashed line, remains the same in the validation set. Note the change in scale. Each model indicates the number of samples correctly identified of the total number of samples from that locality and is given as a percent success. Because Mozambique is the final locality distinguished, it has the same success as the final model—Nigeria. Overall, 95% of validation spectra were correctly classified. See text for details.

Each PLSR model identifies spectra that belong to one group (i.e., the geographic locality). After a group is distinguished, those samples are removed from the data set and all subsequent models. In this case, São José da Batalha samples in Model 1 are removed. The order of the models may be critical to obtaining sufficient separation of samples. Each model is determined by choosing the compositionally most distinct group at each step, as defined by the relationships on a PCA score plot. Because the most distinct group is always eliminated, the samples near the final decision tree are those with the most compositional similarities. Typically, samples in those groups are indistinguishable from each other when modeled in the presence of the other samples, but the small differences between them can be extracted and used to separate these groups when they are modeled in isolation after the other groups are removed.

Electron microprobe analysis (EMP)

To test the applicability of the multivariate statistical approach on widely available tourmaline compositional data from EMP, a subset of 15 tourmaline samples for which LIBS data were obtained (Figs. 4a–4c; 5 grains Brazil; 6 Mozambique; 4 Nigeria; and two additional samples), were analyzed by wavelength-dispersive spectrometry using the JEOL 8230 electron microprobe at LSU. Quantitative compositional analyses for major and minor elements were obtained at an accelerating potential of 15 kV and a 10 nA beam current using a 2 μm spot size, with Na analyzed first. Natural minerals and synthetic materials were used as standards including andalusite (Al), diopside (Ca, Mg, Si), fayalite (Fe), chromite (Cr), kaersutite (Ti), rhodonite (Mn), willemite (Zn), chalcopyrite (Cu), galena (Pb), albite (Na), sanidine (K), fluorite or fluor-phlogopite (F), tugtupite (Cl) with synthetic Bi₂Te₃ (Bi), V-diopside glass (V), and GaAs (Ga). EMP detection limits are given in the Online Materials^[1]. Li, H, or B cannot be effectively analyzed by the EMP and were not included in the data modeled. Two well-characterized elbaite tourmalines served as secondary standards. Count times for major elements were 10 s on the peak, 20 s on the background, and for minor and trace elements 60 s peak, 30 s background. Analytical precision is estimated to be ±1% relative for the major elements and ±5% for the minor elements. Where color zoning is apparent, analytical traverses were made across the samples; in other cases, 10–30 analytical spots per grain were randomly selected.

Figure 4

Selection of polished Paraíba tourmaline samples, in epoxy, used for EMP data collection with sample numbers. (a) Colorzoned Brazilian Paraíba tourmaline grain from the São José de Batalha (SJdB) and Rio Grande do Norte (BZ) localities. White arrow shows the location of a detailed EMP traverse. (b) Variety of colored Paraíba tourmaline grains from Mozambique. (c) Paraíba tourmaline grains from Nigeria. (see Table 1, Online Materials^[1] Appendix 1, and Online Materials^[1] for more information).

Mineral formulas were normalized following the recommended procedures of Henry et al. (2011) permitting B, H, and Fe³⁺ to be calculated based on stoichiometry and charge balance and Li estimated by the procedures of Pesquera et al. (2016). Calculating atoms per formula unit (apfu) served as an additional quality check for EMP data but the normalized data are not used for the statistical analysis. To avoid calculation artifacts, oxide weight percentages of measured elements were used for multivariate statistical modeling and are given in the Online Materials^[1].

Evaluating the efficacy of multivariate statistical models for separating the provenance of Paraíba tourmaline using EMP data followed the same methodology as for separating the LIBS data. However, only 18 variables per chemical analysis are available for modeling. Although the data set comprised 295 analyses, only 15 samples were analyzed. All analyses for each sample were restricted to either the calibration or the validation set to ensure that the models focused on fundamental characteristics of the tourmalines rather than simply identifying analyses from the same sample. Because of the low number of samples, calibrations were based on analyses from 2–4 samples per country, and models were validated with 2 samples from each country. As a result, the calibration set comprised analyses from 4 (Mozambique), 3 (Brazil), or 2 (Nigeria) samples, and the validation set comprised analyses from 2 samples from each country.

Multivariate statistics using LIBS data

LIBS spectra (unmasked) for the Cu-bearing elbaites display prominent Na, Al, Si, Li, and B peaks, in addition to Cu and Mn peaks as expected (Fig. 5). In several samples, LIBS detected minor and trace elements such as K, Mg, Bi, Zn, Ga, and Sr. The presence of these elements was confirmed by previous LA-ICP-MS analyses of Paraíba tourmaline (Z. Sun, personal communication). Although Ca and Mg are minor components, the high intensity of these emission lines reflects the relatively low ionization energy of the alkaline earth elements (Cremers and Radziemski 2013).

Figure 5

Representative LIBS spectrum from each of the four different localities for Paraíba tourmaline, stacked to show alignment of peaks. Brazilian localities are separated into: Brazil, SJdB for São José de Batalha; and Brazil, RGdN for Rio Grande do Norte. Selected major and minor elements are labeled. The black vertical line combines several features, two Ca and two Al emission lines within the labeling line.

The decision tree for these sample suites consists of three models (Fig. 6; Dutrow et al. 2019). In an initial PCA that includes all the tourmaline spectra from the four localities (São José da Batalha, Brazil, SJdB; Rio Grande do Norte, Brazil, RGdN; Mozambique; and Nigeria), no single group clustered tightly and the groups overlapped in PC1-PC2 space (Fig. 7). The SJdB spectra were chosen as the first group to model because the São José da Batalha, Brazil PLSR model had the highest success rate of all possible first models. Model 1, which classifies spectra as either belonging to the SJdB group or to the group of all other tourmalines, is excellent (Fig. 6), despite the overlap of groups in PCA space (Fig. 7). The calibration shows separation between the groups with a VAD of 0.45 (Fig. 3). The validation is 96% successful, correctly classifying 25 of 26 samples. The one false positive is a sample of Nigerian tourmaline classified as SJdB.

Figure 6

Decision tree for PSLR modeling of LIBS spectra for Paraíba tourmalines shown with percent correctly predicted (success) for each locality. After locality samples are modeled, they are removed from all subsequent models. Numerical value of 1 refers to samples belonging to that model data set, 0 indicates all others. See Figure 3 and text for details.

Figure 7

PCA score plot calculated using LIBS spectra of Paraíba tourmaline samples from four localities, shown by different symbols. Samples from different localities overlap and lack distinct data clustering per locality. PC1 accounts for 36% of the variance in the data set; PC2 accounts for 13%. PCA plots are used to determine the sequence of PLSR models.

The spectra of SJdB tourmalines were removed from all subsequent models. Model 2 classifies spectra as belonging to RGdN or to the group of all other tourmalines (Mozambique and Nigeria). There is a clear separation between the two groups in the calibration of Model 2 (Fig. 3), which used a VAD = 0.50 value. The validation is 96% successful, correctly predicting the provenance of 22 of 23 samples (Figs. 3 and 6). Again, one Nigerian sample yielded false positive results. This sample is the same as that which was incorrectly classified as SJdB in Model 1.

Finally, Model 3 discriminates between tourmaline spectra from Nigeria and Mozambique (Fig. 6). Spectra are well separated in the calibration with a VAD = 0.52 (Model 3; Fig. 3). The calibration is 94% successful, correctly classifying 16 of 17 samples. One Nigerian sample was misclassified as belonging to the Mozambique group; however, it is a different sample than the false positive sample in Models 1 and 2. The consistent misclassification of Nigerian samples suggests that the sample set is too small to be representative of the actual dispersion of compositions. Alternatively, on visual examination, this sample has a saw mark, which may have left a surface contamination or varied the surface texture of the sample, affecting plasma properties. Overall, the decision tree correctly classified 63 of 66 spectra (one spectrum per sample), resulting in a cumulative prediction accuracy of 95%. The overall true positive rate (only considering the location assigned to PV 1) is 94% (16 of 17; Fig. 3).

Based on the success of the previous geographic modeling, the geographic origin of two unknown Brazilian samples was predicted. Using the LIBS decision tree developed, both unknown samples are classified as being from the Rio Grande do Norte, Brazil, locality.

Multivariate statistics using EMP data

A more widely used analytical technique for characterizing tourmaline mineral chemistry is by electron probe microanalyses (EMP). As such, this multivariate statistical approach was developed using an EMP analytical data set obtained for a subset of the tourmalines for which LIBS data had been acquired (see Online Materials^[1] for all oxide weight percentages used for multivariate statistics). Importantly, in addition to the major elements, Cu, Mn, and F are present in the tourmalines in amounts readily analyzed by the EMP. F is not easily detected by LIBS but is with EMP. V, Cr and Pb are at, or below, EMP detection limits (Online Materials^[1]).

Modeling EMP data with multivariate statistics followed similar procedures as the modeling for the LIBS data. Because of the smaller sample set size, both Brazilian localities were combined. The character of the EMP data set is different than the LIBS data set, in which each sample is represented by a single spectrum. For EMP data, 10–30 points were analyzed for each of the 15 tourmaline samples (Brazil: 5; Mozambique: 6; Nigeria: 4), resulting in a total of 295 analyses. This data set captures the variability within each sample well, but there are too few samples to be representative of the variability within each country of origin.

A PCA score plot for the calibration EMP analyses in the models shows good clustering for analyses of each tourmaline sample but lacks distinct clustering of samples from each country (Fig. 8). Some relationships are consistent with those found via LIBS. For example, the Brazilian samples plot at negative values of PC2 over a large range of PC1 values. The Nigerian and Mozambican samples cover broad areas that intersect near the origin of the score plot. Analysis of more samples could help the PCA discern different relationships that might provide better separation of the groups.

Figure 8

PCA score plot for all EMP analyses used in this study. Specific samples and localities are given in different colors and symbols by sample number. Analyses from each sample plot in discrete clusters, but clear distinctions between tourmalines from different countries are not apparent.

PLSR is a supervised method where the variables (EMP analyses) are correlated with known provenance variables (PV). Because of this, PLSR models can be successful, regardless of messy relationships in PCA. Model 1 in the EMP decision tree (Fig. 9) separates Brazilian tourmaline analyses from the group of Mozambican and Nigerian samples. The calibration is 96% successful, correctly predicting the origin of 173 of 180 calibration samples with a VAD of 0.51 (Fig. 10). The validation is also 96% successful, correctly predicting the origin 110 of 115 analyses. Five Brazilian analyses are predicted to belong to the group of all others; there are no false positives.

Figure 9

PLSR decision tree based on EMP data with percent of the samples correctly identified as belonging to the known locality shown as success. Both localities in Brazil were grouped together because of the small sample set and represented by “Brazil”. See Figure 6 for details.

Figure 10

Calibration (left) and validation (right) results for EMP decision tree. The dashed line indicates the selected value of apparent distinction (VAD). Each colored symbol represents a different locality as given. Overall, 87% of validation analyses were correctly classified.

Model 2 is more complex. The calibration (Fig. 10) establishes relatively consistent Provenance Variable values for Mozambican calibration analyses with an average near 1 (average = 0.91; range = 0.36–1.36; standard deviation = 0.18). In contrast, the PV values calculated for Nigerian samples, while less than 1, are different from each other. One sample clusters at an average of 0.45 and the other with an average of 0.03 (Fig. 10). Because one Nigerian sample has relatively high calculated PV values, the VAD that results in the best model success is 0.62. This VAD value results in a calibration accuracy of 97% (112 correct predictions of 115), with two false negatives and one false positive. However, this VAD value is not the best choice for the validation (Fig. 10). A higher VAD would have yielded a higher success, as all of the Mozambican-validated analyses have comparable high calculated PV values, as do 19 of the 40 Nigerian validation analyses. This results in a prediction accuracy of 75% for this model (Figs. 9 and 10; 56 of 75 analyses). More samples with analytical data are needed to calculate more successful models. Overall, the EMPA decision tree correctly predicts the country of origin of 87% of the analyses.