Geochemical Signatures of Mafic Volcanic Rocks in Modern Oceanic Settings and Implications for Archean Mafic Magmatism

For over 40 years, geochemical signatures of volcanic rocks have been used to infer possible geodynamic settings of ancient greenstone belts and their contained mineral deposits (e.g., Lesher et al., 1986; Paradis et al., 1988; Kerrich and Wyman, 1997; Hart et al., 2004; Piercey, 2011). Comparisons with the modern oceans, however, have been limited by a lack of samples. This is changing with increased ocean exploration and the availability of high-precision analytical data on submarine volcanic rocks from diverse locations. In a companion paper, we compiled data on more than 2,200 samples of felsic volcanic rocks from 70 different locations in the oceans and compared them to Archean rhyolites and dacites that host some of the world’s most important volcanogenic massive sulfide (VMS) deposits. We found significant geochemical diversity in the modern samples, spanning the full range of compositions of felsic volcanic rocks found in greenstone belts (Fassbender et al., 2023). The present study takes a broader look at the geochemistry of mafic volcanic rocks (<60 wt % SiO₂), which are found in nearly every niche of the modern oceans and are the dominant lithology of greenstone belts. We compiled high-quality geochemical analyses of more than 2,800 unique samples of submarine mafic volcanic rocks from mid-ocean ridges (MORs), ridge-hotspot intersections, intraoceanic arc and back-arc spreading centers, and ocean islands (Fig. 1) and compare these data to well-studied assemblages of the Abitibi greenstone belt of Canada, the largest and best-preserved Neoarchean greenstone belt in the world. Precise analogs of ancient greenstone belts in modern settings seem unlikely due to the hotter mantle in the Archean, thicker oceanic crust, a potential lack of continental basement, and subduction regimes that were in their infancy (e.g., Bédard et al., 2013; Wyman, 2013). However, many of the fundamental petrogenetic indicators are very similar between modern and ancient systems. Major and trace element geochemistry have been widely used to infer conditions, such as high heat flow associated with rifting, shallow depths of melting, and increased fractionation, that are favorable for the formation of VMS deposits (e.g., Galley et al., 1995; Gibson et al., 2007; Piercey, 2011). However, the most common lithotectonic discrimination diagrams for mafic volcanic rocks have been developed by inspection of relatively small data sets. Here, we use a much larger compilation to identify the geochemical characteristics of mafic volcanic rocks in a full range of modern settings. We also used agglomerative hierarchical clustering and principal component analysis (PCA) to identify geochemical differences and similarities between the data sets. The classes identified in the data reflect highly variable melt sources and processes at different stages of basin evolution and reveal unexpected complexity in rapidly evolving microplate domains. We suggest these findings have implications for understanding similar magmatic processes in older terranes such as the Abitibi greenstone belt.

Mid-ocean ridges account for 75% of global mafic magmatism, with the balance in back-arc basins, volcanic arcs, and ocean islands (Perfit and Davidson, 2000). Mafic melts in these different settings are produced by partial melting of the mantle due to (1) increased temperatures, as in hotspot-related regimes, (2) decreasing pressure, as in mid-ocean ridge spreading centers, and (3) addition of volatiles in subduction zones. These major differences in melting process and different sources involved result in very different trace element geochemical signatures of the basalts (e.g., Arevalo and McDonough, 2010; Hofmann, 2014; Ueki et al., 2018). In the following, we review the systematics of mafic melt production in modern oceanic settings (Fig. 1) and summarize the key geochemical characteristics (Tables 1 and 2). As shown in our previous paper (Fassbender et al., 2023), the petrogenesis of the mafic rocks has an important impact on felsic rock compositions, as fractionation of mafic melts or melting of the associated basaltic crust are the main sources of felsic magma.

Mid-ocean ridges are dominated by tholeiitic basalt (mid-ocean ridge basalt [MORB]) produced by decompression partial melting of the mantle (McKenzie and O’Nions, 1991). High degrees of partial melting result in low incompatible element concentrations, including light rare earth elements (LREEs) and large ion lithophile elements (LILEs), relative to the more compatible heavy rare earth elements (HREEs) and high field strength elements (HFSEs) (Hofmann, 2014). Complexity is introduced by varying magmatic heat input and hydrothermal cooling, which change with spreading rate (Michael and Cornell, 1998; Haase et al., 2005). End member examples include the fast-spreading ridges such as the East Pacific Rise at 20-21°S and the Pacific-Antarctic Rise at 37°S, which have stable axial melt lenses developed in a warm lithosphere and thickened oceanic crust, and slow-spreading ridges such as the Mid-Atlantic Ridge at 45°N, where magma chambers are dispersed in a colder lithosphere and magmatism is episodic. These differences result in distinct geochemical signatures, especially at slow-spreading ridges (White and Klein, 2014 and references therein) where greater tectonic extension can lead to strong hydrothermal cooling and fractionation. Mafic melts at ridge-hotspot intersections are produced by lower degrees of partial melting in thicker crustal settings, resulting in intermediate incompatible element concentrations (Jicha et al., 2013). Isotopic systems such as ⁸⁷Sr/⁸⁶Sr and ¹⁴³Nd/¹⁴⁴Nd show that mantle heterogeneity is also important.

Hotspot volcanism, such as at Hawaii and the Marquesas Archipelago, is characterized by deeply sourced alkaline melts produced by low degrees of partial melting of the mantle at high pressure and high temperature (Hofmann, 2014; Haase et al., 2019). This results in the enrichment of incompatible elements (LREEs, LILEs) relative to compatible elements (HREEs, HFSEs), producing the typical signatures of ocean island basalt (OIB). Additional complexity is introduced by mantle reservoirs with higher ⁸⁷Sr/⁸⁶Sr and lower ¹⁴³Nd/¹⁴⁴Nd than at MORs (Jackson et al., 2010; Price et al., 2014).

Arc volcanism in intraoceanic settings, such as the Tofua and Izu Bonin arcs, is mainly calc-alkaline. Melting in the sub-arc mantle is triggered by addition of volatiles from a dehydrating slab, which lowers the solidus temperature of the mantle (e.g., Hawkesworth et al., 1993; Langmuir et al., 2006). The volatiles are sourced from sediments, pore fluids, and hydrous minerals (from earlier alteration of the oceanic crust) and are released into the mantle wedge via pervasive or channeled flow during dehydration of the slab. The melts have high water contents, strong LILE enrichment, and high oxidation states. Certain trace elements, such as Ba, are mobile in the slab-derived fluids and become enriched in fluid-fluxed melts, while Nb and Yb are immobile, resulting in high Ba/Nb and Ba/Yb ratios (e.g., Pearce et al., 2005).

Back-arc rifts and spreading centers, such as in the Lau basin and Manus basin, produce both calc-alkaline and tholeiitic melts, depending on the volatile input from the slab (Langmuir et al., 2006). This mainly varies with proximity of the spreading center to the adjacent arc and trench, but also with the age of the subducting lithosphere, the angle of subduction, the rate of slab rollback, the trajectory of the overriding plate, and different mantle flow regimes (Sdrolias and Müller, 2006; Castillo, 2012; Schmidt and Poli, 2014; Schellart, 2020). These, in turn, vary with the progress of basin opening, which begins with initial rifting and eventually leads to back-arc spreading (Taylor and Martinez, 2003; Pearce et al., 2005; Langmuir et al., 2006). With basin opening, melting regimes that may initially overlap become separated into zones of hydrous melting proximal to the arc and dry (MORB-like) melting of depleted mantle in the deeper back-arc basin. At the Valu Fa Ridge in the Lau basin and the Eastern Rifts of the Manus basin, calc-alkaline melts are produced where the rifts are propagating into arc crust. At more distal spreading centers, such as the Northwest Lau spreading center and the Manus spreading center in the Manus basin, the melts are tholeiitic.

Boninitic melts have been found in a wide range of modern oceanic settings, including back-arc spreading centers and mantle plumes (Cooper et al., 2010; Resing et al., 2011; Golowin et al., 2017; Pearce and Reagan, 2019; Pearce and Arculus, 2020). They are produced by melting of previously depleted mantle material. In subduction zones, early rollback of the subducting slab triggers melting that leaves behind a depleted mantle (Arculus et al., 2019). Melting of this more refractory material is promoted by eventual addition of volatiles and increased heat. The resulting melts are widely interpreted to be indicative of subduction initiation (Arculus et al., 2019; Pearce and Arculus, 2020). However, similar melts are known from back-arc spreading centers and plume-related regimes (Pearce and Reagan, 2019). In the Fonualei rift and spreading center of the northern Lau basin, boninites are variably interpreted to be products of melting of subarc mantle that was depleted by earlier episodes of back-arc spreading or by melting of refractory material related to the Samoan plume (Falloon et al., 2007; Cooper et al., 2010; Resing et al., 2011; Escrig et al., 2012). In mantle-plume settings, the melting is thought to be initiated by a pulse of hot material from the plume (Pearce and Arculus, 2020). The origins of these melts are of particular interest, as they have been widely recognized in greenstone belts (Kerrich et al., 1998).

Geochemical analyses of mafic lavas (<60 wt % SiO₂) from 433 locations in the oceans were compiled from the literature (Fig. 1). From an original database of more than 3,800 analyses, a subset of 2,857 unique samples was selected for this study. The rejected samples included duplicate samples, repeat analyses of the same sample, samples analyzed by inappropriate methods, and samples with location errors or incomplete descriptions. Only samples with complete major and trace element data were retained; every sample in the final data set was analyzed for at least Ti, Al, Ba, Th, Zr, Nb, La, Y, and Yb. No samples were included where the concentration of any element was equal to or below its detection limit or exceeded the upper limits of detection or where the limits of detection were not specified. To check for altered samples, we compared the data to fresh rock compositions of basalts and basaltic andesites, including immobile element ratios and REE profiles. In particular, we excluded samples with unusually high or low mobile element concentrations, such as Si, Na, K, and Eu. The final database is provided in Appendix Table A1.

The MORBs in the database have low SiO₂ (48.7–52.1 wt %), low K₂O (0.01–0.47 wt %), and intermediate TiO₂ (0.55–2.57 wt %) and are depleted in LREEs. Mafic volcanic rocks from ridge-hotspot intersections have variable SiO₂ concentrations (43.7–54.8 wt %) but high K₂O (0.2–2.3 wt %), high TiO₂ (1.6–3.9 wt %), and LREE enrichment. Ocean island basalts also have high K₂O (0.2–2.0 wt %), high TiO₂ (2.4–4.5 wt %), and LREE enrichment but low SiO₂ concentrations (42.4–51.2 wt %). Samples from intraoceanic back-arc and island arc assemblages have high SiO₂ concentrations (44.9–58.9 wt %), intermediate K₂O (0.01–1.1 wt %), and low TiO₂ (0.1–2.4 wt %) and are only moderately LREE enriched. Samples from intracontinental arc-back arc assemblages have high SiO₂ concentrations (45.9–58.6 wt %), intermediate K₂O (0.2–1.9 wt %), intermediate TiO₂ (0.6–2.8 wt %), and moderate LREE enrichment. Figure 2 (see App. Fig. A1 for details) shows the range of mafic volcanic rock compositions in the database on a total alkali-silica plot (after Le Maitre, 1989) and an immobile element plot (after Winchester and Floyd, 1977). The original Winchester and Floyd (1977) diagram is based on a limited number of samples, and so the full range of modern oceanic basalts considered here is not well classified by the depicted field boundaries. Pearce (1996) used a more comprehensive data set; however, some of the modern oceanic basalts are still not well classified (App. Fig. A2). Data for mafic volcanic rocks from the Abitibi greenstone belt are shown for comparison and discussed further below.

Figure 3 shows the range of magmatic affinities of modern oceanic rock (after Ross and Bédard, 2009). The mafic volcanic rocks from Juan de Fuca Ridge and the East Pacific Rise, Galapagos spreading center, Pacific Antarctic Rise, Mid-Atlantic Ridge, and South Indian Ridge plot mainly in the tholeiitic field. Samples from ridge-hotspot intersections and ocean islands, such as Ascension Island, Iceland, Hawaii, and volcanoes of the Marquesas Archipelago, are mostly transitional to alkaline basalts, basanites, and basaltic trachytes. Samples from intraoceanic back-arc assemblages, such as in the Manus basin, Lau basin, Mariana Trough, East Scotia Ridge, and New Hebrides, are similar to MORBs (i.e., dominantly tholeiitic basalts). Samples from island-arc assemblages, such as the Mariana, Izu-Bonin, Fiji, Tonga-Kermadec, and Lesser Antilles volcanic arcs, span the range from tholeiitic to calc-alkaline basalts and andesites. Mafic volcanic rocks from intracontinental arc-back arc settings, such as the Okinawa Trough and Western Antarctic Peninsula, are mostly calc-alkaline andesitic basalts. The corresponding tectonic affiliations are shown in the discrimination plot of Pearce (2014).

Extended trace element plots of the mafic volcanic rocks in different settings are shown in Figure 4. They are grouped according to typical MOR-like signatures (Fig. 4A, B), intraoceanic arc-back arc rift and spreading center signatures (Fig. 4C, D), and arc-back arc signatures (Fig. 4E, F). The MOR mafic volcanic rocks, represented by the East Pacific Rise and Mid-Atlantic Ridge, have mostly flat REE patterns and no negative Nb-Ta anomaly and plot in the N-MORB array of Pearce (2014). Samples from ridge-hotspot intersections and ocean islands have more steeply dipping REE patterns, no negative Nb-Ta anomaly, and plot in the E-MORB array. Mafic volcanic rocks from intraoceanic back-arc assemblages, represented by the East Scotia Ridge, have mostly flat REE patterns with strong negative Nb-Ta anomalies and low HREE concentrations and plot in the oceanic arc array. Samples from the terminations of back-arc spreading centers, represented by the Northeast Lau spreading center, the Southern East Scotia Ridge, and the northern volcanic tectonic zone of the Mariana Trough, all have weak negative Nb-Ta anomalies and dipping REE patterns and plot between the E-MORB and continental arc arrays. Samples from back-arc rifts that are propagating into arc crust, represented by the Fonualei rift and spreading center, have even stronger negative Nb-Ta anomalies and low HREE concentrations. Mafic volcanic rocks from intracontinental arc-back arc settings, represented by the Okinawa Trough, have dipping REE patterns and intermediate negative Nb-Ta anomalies and overlap with the continental arc array.

Lithotectonic discrimination diagrams, such as those of Pearce (2014) and Ross and Bédard (2009), rely on a limited number of trace elements and trace element ratios that can be interpreted in nonunique ways. For example, the Zr/Y ratio is controlled by very slight differences in compatibility (Pearce and Norry, 1979), and the Nb/Yb ratio is strongly dependent on highly variable source depletion (Pearce, 2008). Th/Yb ratios reflect not only volatile input but also assimilation of subduction-related crust (Pearce, 2014). Increased Th may also be related to higher temperatures that are required to mobilize this element. These difficulties can be overcome by considering multiple trace element concentrations simultaneously, as we did in our study of oceanic rhyolites and dacites (Fassbender et al., 2023). We applied multivariate statistical techniques (principal component analysis [PCA] and agglomerative hierarchical clustering) to identify and verify the geochemical differences between sample suites using multielement data rather than individual element ratios. In this study we chose nine elements: Ti, Al, Ba, Th, Zr, Nb, La, Y, and Yb. All were recalculated to elemental concentrations in ppm, and centered-log ratio (CLR) transformations were applied to overcome closure effects (Aitchison, 1982; see App. Table A2). The PCA and clustering were performed with the R statistical software package (R Core Team, 2020).

We used PCA to illustrate the minimum number of “components” that account for the largest variance in the data set (Bartholomew, 2010). Each principal component (PC) is affected by all the variables; the first and second components, PC1 and PC2, account for most of the variance in the data set (66 and 22%, respectively; see App. Table A3). The contribution of the different elements is described in terms of loadings of the element on each component. In Figure 5, the red labels represent the loadings of the elements for the entire data set (see App. Table A4). The plotted points are individual sample scores in PC1-PC2 space, given in Appendix Table A5. Yb and Y have strong positive loadings on PC1, whereas Th and La have strong negative loadings, similar to the felsic volcanic rocks (Fassbender et al., 2023). Zr and Nb have strong positive loadings on PC2, whereas Ba and Al have strong negative loadings. The PCA shows very consistent behavior among the mobile elements (e.g., Ba and La) and less mobile elements (Th and Al).

To identify geochemically distinct groups of samples in the data set, we used unsupervised agglomerative hierarchical clustering of the CLR-transformed data (App. Table A6). This type of cluster analysis can identify subtle differences in geochemical signatures between sample types, particularly for volcanic rocks that show small but systematic variations between lithotectonic settings. In Fassbender et al. (2023), we performed cluster analysis on the PCA results for the felsic volcanic rocks rather than the raw geochemical data. Here, we performed clustering on the CLR-transformed trace elements of the mafic volcanics, which yielded much better separation of individual clusters. The analysis was performed in R using the hierarchical clustering function, “hclust.” Agglomerative hierarchical clustering uses a set of queries to classify samples in a way that maximizes the differences between the groups. In this type of analysis, no preassigned groups or clusters are identified (i.e., the analysis is unsupervised). The analysis starts by assuming that every individual sample forms its own cluster, and then the statistical differences (or “distances”) from all other samples are calculated. The two samples with the closest distance to each other form the next level of clusters. The distances are calculated using Ward’s criterion, which determines the center of each cluster and then the sum of the squared distances of individuals from the center. Each sample was assigned to one cluster (App. Table A6). The smallest number of clusters representing all of the samples was determined using an iterative approach discussed below. The workflow is presented in Appendix Figure A3 and can be replicated with the code and instructions provided in the Appendix. Seven different clusters (C1-C7) were selected to train our classifier and are outlined in Figures 5 and 6.

The results of the unsupervised clustering can be used to train a classifier. We did this by performing a supervised classification on the newly defined clusters identified above. We then conducted a blind test to determine if the classification scheme can accurately predict geodynamic influences in the mafic volcanic rock geochemistry. Supervised classification is a machine-learning task used for large data sets that sorts samples into known classes based on a “training” data set. The training is an iterative process by which a machine-learning classifier, such as random forest (RF), learns what elements achieve the best classification outcome and then builds a model that can be used as a classification scheme. Random forest has proven to be particularly useful for supervised classification based on multivariate data sets (Breiman, 2001; Fernandez-Delgado et al., 2014). The steps involve initial processing of the data to create the training set, training of the classifier, and then evaluation of its prediction success. The input is a selection of parameters (i.e., the analyzed elements) and the predetermined clusters (i.e., target classes). The starting database is randomly split into the training data set, used to build and validate the RF classifier, and the blind test data used to evaluate the prediction rate.

RF functions using multiple decision trees with true or false outcomes (e.g., “are SiO₂ concentrations <60 wt % ?”; if yes, the sample is a mafic rock; if not, the sample is a felsic rock). The decision points are referred to as “nodes,” and the cutoffs are determined by the classifier’s ability to separate the training data. For example, the correct cutoff for SiO₂ is identified by testing several SiO₂ concentrations until the best separation of the clusters is achieved. “Pure” separation is achieved if the cut-off value of 60 wt % SiO₂ correctly identifies all basalts. At every node, the “impurity” of the separation can be measured in terms of the probability of a parameter falsely classifying a sample. This is commonly determined using the “Gini index”:

where “t” is the node in the decision tree (e.g., Ti), “c” is the target class, “j” refers to the possible classes, “g_c” is the relative frequency of classification, “n_c” is the number of samples belonging to target class “c,” and “n” is the total number of samples evaluated at this node. Nodes, or parameters (e.g., Ti), with a low Gini index have a high probability of correctly classifying a sample and therefore greater importance in making the decision.

Machine-learning classifiers may work very well with the training data but perform poorly in blind tests. Random forest overcomes this problem by using a large number of uncor-related decision trees operating as a “committee” that casts votes for each decision. Many decision trees are created by randomly selecting elements from a training data set as the nodes. The least number of decision trees (usually 500 or more) is selected for the final random forest model. Typically, after this many decision trees, the performance of the classifier does not improve, and the smallest number is chosen to lower computational cost and minimize any potential artifacts in the data analysis. Each decision tree is then “trained” with randomly chosen samples equal to the total number of samples in the training data base, a method commonly known as “bootstrapping.” During this process samples can be picked multiple times, and a pool of samples that were not picked and therefore not involved in training (so-called “out of bag” samples) is left behind. This is typically about one third of all samples. These “out of bag” samples are used to validate the classification scheme (see below). To classify an unknown sample, each of the randomized decision trees casts a single “vote” for a single target class or cluster from the training set. The votes cast by all of the decision trees for each cluster are given in percentages (summing to 100%), and RF uses the majority vote to finally classify a sample (Breiman, 2001). A disadvantage of the RF classifier is that it will force a decision regarding the classification of a sample even if the sample does not match any of the target classes. Forced decisions commonly return a randomly distributed set of votes.

We selected 80% of our database at random (see App. Table A7) to train and validate our RF classifier; the remaining 20% was used for the blind test (see App. Table A8), which is a standard approach (Gholamy et al., 2018). We used the target classes identified by the unsupervised hierarchical clustering to train the classifier. Different numbers of samples assigned to the target classes in the training data set can create a class imbalance and lead to poor performance of the classifier. To overcome this, we equalized the number of samples in each target class in the training set using the synthetic minority oversampling technique (SMOTE: Chawla et al., 2002). SMOTE randomly “removes” samples from classes with too many samples, and “creates” new samples for classes that have too few by slightly altering the geochemical composition of individual samples until each class contains the same number of samples. In our study, this resulted in a modified training set with 2,288 samples (App. Table A9).

To validate the model, we used the subset of samples remaining from the bootstrapping procedure (i.e., “out of bag” samples). Of the “out of bag” samples from all target classes, 97.2% were correctly classified (App. Table A10). The probability of misclassification by a particular node, averaged over the entire RF classification scheme, is referred to as the Mean Decrease Gini. This parameter identified Nb, Th, Y, Yb, and Al as the most important variables in determining whether a sample is classified correctly (Fig. 7). Nodes corresponding to La and Zr have low Mean Decrease Gini values and therefore less influence in the correct classification of a sample. We then conducted a blind test on the 20% of the original samples in the database that were not used in the training set (see App. Table A8; none were created by SMOTE) to determine how well the classifier performs on new data. The test samples were correctly classified in almost every case, identifying 96% or 547 of 569 samples (App. Table A11). Because RF does not use any element independently, this significantly reduces the effect of potentially anomalous behavior of any individual element for individual samples. To exclude any potential control of Ba mobility on RF classification, we created a second classifier excluding Ba yielding almost identical results (App. Table A12).

We identified seven target classes of mafic volcanic rocks from the modern oceans, corresponding to distinct geodynamic settings and different melting conditions. Based on these target classes, we constructed an RF classifier that accurately identifies rocks from those settings. Here, we discuss the tectonic and petrogenetic significance of the target classes and then assess the application of the RF classifier to ancient volcanic rocks.

In the PCA, we interpreted the loadings on PC1 to represent different mantle sources, from depleted MORs to enriched ocean island settings, whereas the loadings on PC2 mainly reflect different melting regimes, from dry melting at MORs and distal back-arc basin spreading centers to wet melting in near-arc environments. Mafic volcanic rocks from MORs and mature back-arc spreading centers have positive loadings on PC1 and PC2 (App. Table A5) and belong to target class C1 (Table 3; App. Table A6). These rocks have low Nb and La, typical of MOR magmas and an incompatible element-depleted mantle source. Mafic rocks from less-depleted MORs and back-arc spreading centers have intermediate positive loadings on PC1 and positive loadings on PC2 and belong to target class C2. Mafic rocks from enriched MORs and back-arc spreading centers have slightly negative loadings on PC1 and positive loadings on PC2 and belong to target class C3, with low to intermediate Nb and La. Mafic rocks from ocean island-like settings have negative loadings on PC1 and positive loadings on PC2 and belong to target class C4, with mantle sources that are highly enriched in incompatible elements (i.e., high Nb and La). Mafic rocks from immature back-arc rifts have small positive loadings on PC1 and negative loadings on PC2 and belong to target class C5, influenced by at least some fluid-fluxed melt components and having intermediate Al and Ba. Mafic rocks from arc volcanoes and arc-related rifts have intermediate loadings on PC1 and negative loadings on PC2 and belong to target class C6, reflecting mostly fluid-fluxed melts. These rocks have high Al and Ba, typical of island arcs. Finally, mafic volcanic rocks from nascent back-arc rifts, especially in continental margin settings, have small loadings on PC1 and negative loadings on PC2 and belong to target class C7, including both fluid-fluxed melt components and enriched mantle or crustal sources. These rocks have high Ba, Nb, Th, and La, at least partly due to crustal contamination.

The cluster analysis illustrates the significant geodynamic control on the trace element signatures (Figs. 5, 6, and 8): C1, C2, and C3 are MOR-like assemblages; C4 corresponds to ocean island-like assemblages; and C5-C7 correspond to arc and back-arc assemblages (Table 3). The type localities for clusters C1-C3 are Juan de Fuca Ridge, the East Pacific Rise, the Galapagos spreading center, the Pacific Antarctic Rise, the Mid-Atlantic Ridge, and the South Indian Ridge as well as mature back-arc spreading centers such as the Central Lau spreading center (Lau basin) and the Manus spreading center (Manus basin). Samples in C1 (e.g., East Pacific Rise and Central Lau spreading center) have very low La/Sm, La/Yb, Th/Yb, Nb/Yb, and Ba/Nb, and no negative Nb-Ta anomaly. Samples in C2 (e.g., Mid-Atlantic Ridge and Central Scotia back-arc spreading center) have low La/Sm, La/Yb, Th/Yb, and Nb/Yb and very low Ba/Nb, with no negative Nb-Ta anomaly. Samples in C3 (e.g., Galapagos spreading center and Northwest Lau spreading center) have intermediate La/Sm, La/Yb, Th/Yb, and Nb/Yb, very low Ba/Nb, and no negative Nb-Ta anomaly. The type localities for cluster C4 are the Samoa hotspot, Hawaii and the Marquesas Archipelago (ocean island settings), Ascension Island (a ridge-hotspot intersection), and other settings with ocean island-like signatures, such as the Mid-Atlantic Ridge north of Iceland. Samples in C4 have high La/Sm, La/Yb, Th/Yb and Nb/Yb, low Ba/Nb, and no negative Nb-Ta anomaly. Samples in C5 include mafic volcanic rocks from the Myojin rift (Izu-Bonin) and Ngatoroirangi rift (Havre Trough), which have intermediate La/Sm, La/Yb, and Th/Yb, low Nb/Yb, intermediate Ba/Nb, and moderate negative Nb-Ta anomalies. Samples from the Fonualei rift and spreading center (Lau basin) and the Tofua volcanic in C6 have intermediate La/Sm, high Th/Yb, low Nb/Yb, high Ba/Nb, and strong negative Nb-Ta anomalies. Finally, samples from the nascent back-arc rift of the Northeast Lau spreading center (Lau basin) and continental margin arc settings, such as the Bransfield Strait, fall in C7 and have high La/Sm and Th/Yb, intermediate Nb/Yb and Ba/Nb, and moderate negative Nb-Ta anomalies. The results of training the RF classifier using these different clusters is discussed below.

Highly variable melt sources and conditions, ranging from the thickness of the crust to the composition of the sub-arc mantle wedge and mantle mixing (e.g., corner flow and incursion through slab tears and windows: Taylor and Martinez, 2003; Haase et al., 2009; Cooper et al., 2010; Price et al., 2017), all contribute to the overall complexity observed in modern oceanic basalts. It can be difficult to unravel all of the possible influences from just a few trace-element ratios. By using nine elements simultaneously, even samples with a weak influence from subduction fluids (target class C5) can be distinguished from samples with a stronger influence (target class C6). We can also distinguish truly arc-related assemblages (e.g., target classes C5 and C6) from mafic volcanic suites that are only affected by crustal contamination, such as the Icelandic basalts in target class C3 (Caracciolo et al., 2022; Fig. 4A). Such distinctions often require difficult isotopic studies to unravel.

Although no prior knowledge of the geodynamic settings of the samples influenced the clustering, the classification outcome presented here corresponds closely to established geochemical characteristics associated with different melt sources. The most important are the absence of negative Nb-Ta anomalies in target classes C1-C4, indicating no subduction-derived fluids; increasing La/Sm and Nb/Yb ratios from target class C1 to C4, indicating the transition from depleted mantle sources to enriched sources; negative Nb-Ta anomalies and LILE enrichments with low Th in target classes C5 and C6, related to subduction zone processes; and weak negative Nb-Ta anomalies and LILE enrichments but high Th in target class C7, with enriched-mantle components and crustal contamination. The increase in La/Sm and Nb/Yb ratios from C1 to C3 reflects undepleted mantle in C2 and more enriched mantle in C3, but they all lack high Ba/Nb and LILE (Pearce, 2008; Hofmann, 2014; White and Klein, 2014). Although the high La/Sm and Nb/Yb ratios are related mainly to enriched mantle sources (e.g., Haase et al., 2019), they may also be produced by different degrees of partial melting and variations in the depth of melting (Engel et al., 1965; Wanless and Shaw, 2012).

Mafic volcanic rocks in MOR-like settings and mature back-arc spreading centers, captured in target classes C1-C3, are products of decompression melting due to adiabatic rise of mantle material into MOR spreading centers (McKenzie and O’Nions, 1991; Hofmann, 2014). Although they all have depleted mantle sources with very low La/Sm, Nb/Yb, and no negative Nb-Ta anomaly, the geochemical signatures are influenced by a range of melting conditions. Increased depth of melting can leave garnet in the melt residue, which results in depletion of Yb (Hofmann, 2014). The transition from C1-C3 to C4 reflects deep, low-degree partial melting of anhydrous depleted mantle, characteristic of hotspots and ocean island settings, which results in high La/Sm, high Nb/Yb, and no negative Nb-Ta anomaly (Chauvel et al., 2012; Jicha et al., 2013; Haase et al., 2020). Melting that is influenced to different degrees by fluids from a dehydrating slab is captured in target classes C5, C6, and C7. In this case, the mafic volcanic rocks are products of mixing between hydrous melts derived from devolatilization of the downgoing slab and a water-poor MOR-like melt with negative Nb-Ta anomalies and LILE enrichment compared to MOR basalts (Taylor and Martinez, 2003; Langmuir et al., 2006). Less fluid from the dehydrating slab at back-arc spreading centers distant from the volcanic arc results in weak negative Nb-Ta anomalies, intermediate Ba/Nb, and low Th/Yb (Fretzdorff et al., 2006; Bezos et al., 2009), which is captured in target class C5. The large volatile fluxes above a dehydrating slab that result in strong negative Nb-Ta anomalies, high Ba/Nb, and intermediate Th/Yb, characteristic of arc volcanoes and arc-related rifts (Keller et al., 2008; Escrig et al., 2012), are captured in target class C6. Arc rifting and initial back-arc spreading may be related to processes entirely within the upper plate (e.g., Caratori Tontini et al., 2019) or to complex mantle flow regimes related to the subducting slab (Heuret and Lallemand, 2005; Schellart, 2008). Finally, in nascent back-arc rifts and continental margin rifts, crustal assimilation produces the high Th/Yb and La/Sm (Fretzdorff et al., 2004; Haase et al., 2020) captured in target class C7.

The Lau basin is the type example of an intraoceanic back-arc basin and a key location to explore the geochemical variations in mafic volcanic rocks in response to subduction, microplate interactions, and mantle heterogeneity (Taylor et al., 1996; Keller et al., 2008; Yan et al., 2012; Sleeper and Martinez, 2016). It shows remarkable variations in arc- and non-arc petrogenesis that may be encountered in a single back-arc basin at scales of tens to just a few hundred kilometers. The V-shaped basin is ~450 km wide and ~1,000 km in length (Fig. 9). It is flanked by the active Tofua arc in the east, the Lau Ridge remnant arc in the west, and the Vitiaz lineament or fracture zone in the north (Hawkins, 1995). The opening of the basin, which began approximately 6 m.y. ago and has been propagating southward (Hawkins, 1995; Taylor et al., 1996), is accommodated by the Eastern Lau spreading center in the south, Valu Fa Ridge and Central Lau spreading center in the middle, and the Fonualei rift and spreading center, Mangatolu Triple Junction, Northeast Lau spreading center, and Northwest Lau spreading center in the north. The southern part is tectonically simple, but the north is a complex mosaic of microplates that are variably influenced by the subduction of the Pacific Plate and large-scale transcurrent faulting (Stewart et al., 2022). In addition to the structural complexity, a number of researchers have suggested that the mantle below the back-arc region is characterized by complex flow regimes including contributions from the nearby Samoan mantle plume (Pearce et al., 2007; Escrig et al., 2009, 2012; Price et al., 2014; Yan et al., 2020; Haase et al., 2022). Isotopic studies of Hf, Nd, Sr, and Pb indicate significant mantle mixing and crustal assimilation in the northeastern Lau basin (Regelous et al., 2008; Tian et al., 2011; Price et al., 2014; Nebel and Arculus, 2015).

Mafic volcanic rocks from spreading centers closest to the Tofua arc (VFR and FRSC) have trace element signatures that reflect the subduction, including strong negative Nb-Ta anomalies, low HREEs, high LILEs, and high Ba/Nb ratios (Figs. 9 and 10). Rocks from the Valu Fa Ridge show the greatest influence of volatiles from the dehydrating slab (Haase et al., 2009), indicated by the high Ba/Nb and Ba/Yb (Fig. 10). With increasing distance from the slab, the amount of fluid that enters the melting regime decreases, resulting in decreasing Ba/Nb, Ba/Yb, Th/Yb, and Th/Yb (Fig. 10). Similarly, mafic rocks from the Fonualei rift and spreading center show decreasing influence of the arc (e.g., decreasing Ba/Nb) from south to north (Fig. 10), while Th/Yb increases due to enriched mantle sources and potential assimilation of preexisting arc crust. In contrast, the dominant source for melts along the Central Lau spreading center and Northwest Lau spreading center is the Indian Ocean mantle, with only a minor subduction component (Tian et al., 2008). The Mangatolu Triple Junction, Fonualei rift and spreading center, Eastern Lau spreading center, and Northeast Lau spreading center show a range of contributions from the slab, the modified mantle wedge, as well as Indian mid-ocean ridge basalt (MORB) (Yan et al., 2012; Zhang et al., 2018). Sampling has revealed complex melt sources in a very small area, with mixing of depleted mantle, enriched mantle, and boninitic sources (Haase et al., 2022). The enriched mantle component, indicated by the high Nb/Yb in Figure 10, has been attributed to Samoan plume material crossing into the northern Lau basin (Price et al., 2014), although this is debated (e.g., Falloon et al., 2008; Lupton et al., 2015; Nebel et al., 2018; Haase et al., 2022). Another possibility is that the enriched melts are related to the subduction of the Louisville Seamount Chain; another is assimilation of crustal material noted above (Beier et al., 2017; Price et al., 2017; Schönhofen, 2021). A contribution from the Samoan plume is also considered to be the source of the high Nb/Yb at the Northwest Lau spreading center and Rochambeau rifts (Zhang et al., 2018). High He isotope ratios, which are a diagnostic feature of the Samoan plume, are well known in mafic melts in the Northwest Lau spreading center and Rochambeau rifts (Lupton et al., 2009, 2015). However, volcanic rocks from the Fonualei rift and spreading center and Northeast Lau spreading center lack the He isotope anomaly. Boninitic melts from the Fonualei rift and spreading center and Northeast Lau spreading center, indicated by elevated SiO₂, low TiO₂, and low REE concentrations, reflect enhanced melting caused by a high volatile flux from the subduction zone (Falloon et al., 2007; Regelous et al., 2008; Resing et al., 2011, Escrig et al., 2012).

Mafic volcanic rocks with these highly variable geochemical signatures have been dredged in almost every back-arc basin investigated in this study, including the western Manus spreading center, the Mariana back-arc spreading center, the East Scotia Ridge, and the Bransfield Strait (Fig. 4). The locations of the anomalous samples are almost always at the terminations of the spreading centers, which are typically associated with complex mantle flow (Pearce and Stern, 2006). These locations are variably affected by localized upwelling of mantle, mantle flow around the edges of the slab, and mantle mixing (Fig. 11; Fretzdorff et al., 2002, 2004; Sinton et al., 2003; Pearce et al., 2005; Pearce and Stern, 2006; Beier et al., 2010; Haase et al., 2020). Strong toroidal flow has been proposed in the northern Lau basin (Schellart and Moresi, 2013) and elsewhere (northern East Scotia Ridge; Fretzdorff et al., 2002). This has been the subject of analog and numerical modeling (Chen et al., 2016; Magni, 2019) that suggest flow through a slab tear at the northern termination of the Tonga Ridge is introducing fertile mantle into the back-arc region.

The different mantle domains that are suggested by the hierarchical clustering model are most reasonably interpreted in terms of the complex microplate architecture of the basin (Fig. 12; Baxter et al., 2020; Schmid et al., 2020; Stewart et al., 2022). Enriched mantle material at the Northwest Lau spreading center and Rochambeau rifts is reflected in the predominance of samples in target classes C2 and C3. Mafic volcanic rocks of the Fonualei rift and spreading center are mainly in target class C6, with a shift to target class C7 among the samples from the Northeast Lau spreading center, where there is less fluid-fluxed melting, enriched mantle sources, and possible assimilation of crustal material. Low-degree partial melting of fertile mantle is also observed at several locations, producing OIB-like mafic volcanic rocks of target class C4.

The volcanic and sedimentary rocks of the AGB are widely interpreted to be part of an intraoceanic arc-back arc system with similarities to modern oceanic microplate mosaics (Jackson et al., 1994; Percival et al., 2012; Wyman, 2013). However, other authors have suggested significant differences between Archean and post-Archsean magmatism and tectonics, with some Archean terranes possibly having no modern analog (e.g., Bédard et al., 2013). Submarine volcanism in the Abitibi greenstone belt produced seven temporally distinct assemblages between 2795 and 2695 Ma, traditionally divided into the Northern and Southern Volcanic Zones (Fig. 13A). We used the RF classifier developed in this paper to test whether Archean mafic volcanic rocks of the Abitibi greenstone belt can be classified along the same geochemical lines as modern samples and thereby possibly interpreted in terms of processes like those observed in the modern oceans. The comparison complements our earlier study of modern felsic volcanic rocks, which showed significant geochemical diversity similar to Archean rhyolites from the Abitibi greenstone belt (Fassbender et al., 2023).

We used data from Mole et al. (2021) and several data sets of the Ontario Geological Survey (MRD 378, 292, 362, 377, 393, 355, 085), the Geological Survey of Canada (Open File 6623), and the SIGEOM (Système d’information géominière of Québec). The samples include tholeiitic to calc-alkaline volcanic rocks from the Selbaie area, Geant-Dormant, and Joutel in the Northern Volcanic Zone (NVZ), dominantly tholeiitic rocks from the northeast Matagami area, and tholeiitic to calc-alkaline volcanic rocks from from the Hunter Mine Group, Val-d’Or, Noranda, Bousquet, the western Blake River Group, and the Timmins-Porcupine area, including Kidd Creek and Kamiskotia, in the Southern Volcanic Zone (SVZ) (Fig. 13A). These locations host some of the world’s most important volcanogenic massive sulfide deposits (e.g., Franklin et al., 2005). Following the same QA/QC applied to the modern volcanic rocks, 5,746 samples were selected for the final database (see App. and App. Table A13). Effects of alteration are minimized, although the larger range of K₂O and Na₂O in the Abitibi greenstone belt samples (e.g., Fig. 2) suggests that some altered samples are present. In our analysis of the RF classifier, we showed that Ba mobility had no effect on classification performance for modern sea-floor basalts (App. Table A12). Barium mobility in the older Abitibi greenstone belt samples may lead to lower prediction rates for individual classes. To exclude samples that may have received lower prediction rates due to alteration, only classification results with a prediction rate >50% were included in the analysis. In the Abitibi greenstone belt data set, of the 5,746 samples in the database, 1,065 were not classified (i.e., no class received more than 50% of the votes for the sample).

The mafic volcanic rocks from the Abitibi greenstone belt are distinguished from modern ocean floor basalts by more variable K₂O and Na₂O, lower TiO₂, and flat to steeply dipping REE patterns. The larger range of Zr/TiO₂ could be related to higher temperatures, more contamination in the Abitibi greenstone belt, or a sampling bias. Extraction of melts from the Archean mantle required melting at very high temperatures under mostly dry melting conditions (Prior et al., 1999; Piercey, 2011; Thurston, 2015). However, examples of negative Nb-Ta anomalies and LILE enrichment, which are the dominant features of melting hydrated crust in modern submarine arcs, are found in every assemblage of the Abitibi greenstone belt (Ayer et al., 2002). Samples of mafic volcanic rocks from the Val-d’Or formation, for example, have strong negative Nb-Ta anomalies, low REE concentrations and LILE enrichment (Fig. 14), interpreted to indicate an arc-like regime (Scott et al., 2002). Samples from Noranda have overall weaker negative Nb-Ta anomalies with slightly dipping REE patterns that have been interpreted as originating in a back-arc regime (e.g., Yang and Scott, 2003). Samples from north east of Matagami lack negative Nb-Ta anomalies and LILE enrichment, and they have significantly higher REE concentrations that have been attributed to melt generation at an oceanic spreading center (Hart et al., 2004). The weak negative Nb-Ta anomalies and higher Th concentrations in Abitibi mafic volcanic rocks are commonly cited as important differences compared to basalts from modern oceanic subduction zones (Bédard et al., 2013). Arc-like geochemical signatures could have been created in a number of different ways, including assimilation of hydrated basalt (e.g., amphibole-containing altered crustal components) and/or ilmenite fractionation (e.g., Haase et al., 2005; Pearce, 2008; Mole et al., 2021). The example of Icelandic basalts perhaps best illustrates this. They have high Th/Yb and plot close to the continental arc field in the lithotectonic discrimination diagram of Pearce (2014) (Fig. 4B), although they lack the Ba enrichment that is typical for arc-like rocks. Mafic volcanic rocks nearly identical to those from the Eastern Blake River Group occur at the terminations of modern back-arc spreading centers (e.g., in the Northeast Lau spreading center, the southern East Scotia Ridge, the western Bransfield Strait, and the southern and northern Mariana Trough; Figs. 4 and 14). Similar rocks are clearly identified in the Abitibi greenstone belt by the RF classifier developed from modern oceanic basalts.

In the Abitibi greenstone belt data set, Nb, Ba, Yb, Y, and Th had the greatest influence on the classification of the rocks by the RF classifier. Of the 5,746 samples in the database, 1,065 were not classified (i.e., no class received more than 50% of the votes). This indicates either fundamental differences from modern oceanic mafic volcanic rocks or geochemical variation due to alteration and metamorphism. The remaining 4,681 samples were assigned to a target class comparable to that of modern oceanic basalts. Rocks belonging to C2, C3, C5, C6, and C7 were confidently identified in the Abitibi greenstone belt data set (Fig. 13B; App. Table A14). They can be broadly divided into two groups, which in the modern classification would correspond to dry melting (C2 and C3) and wet melting conditions (C5, C6, and C7). The latter could be interpreted as dehydration of hydrated basalts in a downgoing slab or, in other models, hydrated basaltic crust that sank into the mantle during delamination. None of the samples were fully classified (100% of the votes) into C1 (i.e., depleted MORS) or C4 (i.e., ocean island settings). Samples from Kidd Creek were classified as belonging to target class C2 and C3 (Figure 13C), indicating similarities to ridge-hotspot intersections, such as Iceland, as previously suggested by Prior et al. (1999). The samples from the Blake River Group (and Noranda specifically) show greater geochemical variability. Samples from the western part of the Blake River Group are mostly classified in C2 (i.e., undepleted MOR volcanic suites), whereas samples from the Eastern Blake River Group and Noranda are mostly classified in C7. The latter are most similar to mafic volcanic rocks from nascent back-arc rifts as in the Northeast Lau spreading center and the eastern Manus basin, as previously recognized by Yang and Scott (2003). A high proportion of mafic volcanic rocks in the Bousquet and Val d’Or districts were classified in target classes C5, C6, and C7 (i.e., intraoceanic back-arc and arc-related rifts). This is in agreement with previous studies characterizing the Val-d’Or formation as arc-like crust (Scott et al., 2002).

We compared the major and trace element geochemistry of modern submarine mafic volcanic rocks from different geodynamic settings using PCA, agglomerative hierarchical clustering, and RF. Automated hierarchical clustering reveals significant geochemical diversity related to mantle heterogeneity, highly variable melting regimes, and contributions of fluids from the subducting slab. Back-arc basalt and MOR-like mafic volcanic rocks, such as at the Northwest Lau spreading center and the Pacific Antarctic Rise, as well as ridge-hotspot intersections such as in Iceland, all share enriched mantle signatures typical of spreading centers far from any subduction. However, significant variation in the trace element geochemistry is found in complex microplate mosaics at very small scales (e.g., distances of less than 100 km along the Fonualei rift and spreading center and the Northeast Lau spreading center). The greatest diversity occurs at the terminations or propagating tips of back-arc spreading centers, such as the Northeast Lau spreading center, the northern and southern East Scotia Ridge, and the northern and southern Mariana Trough, owing to complex mantle flow regimes.

Comparisons with mafic volcanic rocks of the Abitibi greenstone belt indicate magmas generated above the hot Archean mantle were most similar geochemically to those of modern back-arc spreading centers with weak arc signatures and enriched mantle sources, similar to the Northwest Lau spreading center of the Lau basin and, in the open ocean, similar to the Pacific Antarctic Rise and Iceland. The mafic volcanic rocks from these settings bear a strong similarity to volcanic rocks from northeast of Matagami and from Kidd Creek. Significant geochemical diversity in mafic rock compositions in the Abitibi greenstone belt, including in the Eastern Blake River Group at Noranda and at Val d’Or, occur at about the same scale as in the northern Lau basin. We interpret these variations to reflect mantle heterogeneity and mantle flow regimes influenced by microplate formation. The ability to recognize possible ancient analogs in greenstone belts raises the possibility of improved area selection for mineral exploration, focusing on the enhanced magmatic productivity of back-arc spreading centers with enriched mantle sources, at the terminal or propagating tips of arc-related rifts, and at ridge-hotspot intersections where mantle upwelling likely provided the heat for melting.

We thank two anonymous referees for the helpful comments that improved the manuscript. This research was funded by the Canada First Research Excellence Fund (CFREF, Metal Earth at Laurentian University) and the Marine Mineral Resources Group at the GEOMAR Helmholtz Centre for Ocean Research Kiel. The Natural Sciences and Engineering Research Council of Canada (NSERC), the Helmholtz Association, and the German Ministry of Science and Education (BMBF, grant 03G0267) are acknowledged for the support of this work through research grants to the authors. The project was also supported by the NSERC Collaborative Research and Training Experience program (iMAGE-CREATE) on Marine Geodynamics and Georesources. This is contribution MERC-ME-2023-27.

Marc Fassbender is a freelance consulting geologist holding a Ph.D. from the University of Ottawa, Canada. He specializes in host rocks of actively forming VMS deposits on the modern sea floor and their ancient analogues, using petrology, geochemistry, and machine learning. Marc holds a B.Sc. and an M.Sc. in geology from the University of Erlangen-Nuremberg, Germany, where his research focused on HREE enrichment in olivine and petrogenetic evolution of the Vergenoeg F-Fe-REE deposit in South Africa. He has experience in the fields of economic geology, geochemistry, and integrating geochemical and other data sets using machine learning to solve geoscience problems.