Felsic volcanic rocks are abundant in ancient greenstone belts and important host rocks for volcanogenic massive sulfide (VMS) deposits. About half of all VMS deposits are hosted by dacite or rhyolite, an association that reflects anomalous heat flow during rifting, partial melting of basaltic crust, and fractional crystallization in high-level magma chambers. For over 30 years, geochemical signatures of these rocks (e.g., F classification of Archean rhyolites) have been widely used to identify possible hosts for VMS deposits in ancient greenstone belts. However, comparisons with modern oceanic settings have been limited, owing to a lack of samples of felsic volcanic rocks from the sea floor. This is changing with increasing exploration of the oceans. In this study, we have compiled high-quality geochemical analyses of more than 2,200 unique samples of submarine felsic volcanic rocks (>60 wt % SiO2) from a wide range of settings, including mid-ocean ridges, ridge-hot-spot intersections, intraoceanic arc and back-arc spreading centers, and ocean islands. The compiled data show significant geochemical diversity spanning the full range of compositions of rhyolites found in ancient greenstone belts. This diversity is interpreted to reflect variations in crustal thickness, the presence or absence of slab-derived fluids (dry melting versus wet melting), and mantle anomalies. Highly variable melting conditions are thought to be related to short-lived microplate domains, such as those caused by diffuse spreading and multiple overlapping spreading centers. Systematic differences in the compositions of felsic volcanic rocks in the modern oceanic settings are revealed by a combination of principal components analysis, unsupervised hierarchical clustering, and supervised random forest classification of the compiled data. Dacites and rhyolites from midocean ridge settings have moderately depleted mantle signatures, whereas rocks from ridge-hot-spot intersections and ocean islands reflect enriched mantle sources. Felsic volcanic rocks from arc-back-arc systems have strongly depleted mantle signatures and well-known subduction-related chemistry (strong large ion lithophile element enrichment in combination with strong negative Nb-Ta anomalies and low heavy rare earth elements [HREEs]). This contrasts with felsic volcanic rocks in Archean greenstone belts, which show high field strength element and HREE enrichment (so-called FIIIb-type) due to a less depleted mantle, a lack of wet melting, and variable crustal contamination. The differences between modern and ancient volcanic rocks are interpreted to reflect the lower mantle temperatures, thinner crust, and subduction-related processes in present-day settings. We suggest that the abundance of FIIIb-type felsic volcanic rocks in Archean greenstone belts is related to buoyant microplate domains with thickened oceanic crust that were better preserved on emerging Archean cratons, whereas in post-Archean tectonic settings most of these rocks are subducted.

In Archean greenstone belts, bimodal volcanic successions host the majority of volcanogenic massive sulfide (VMS) deposits (Barrie and Hannington, 1999; Franklin et al., 2005). A close association with felsic volcanic rocks is widely viewed as a consequence of rifting, during which large-scale hydrothermal systems were developed. Geochemical characterization of the felsic volcanic rocks has focused on their rare earth element (REE) and high field strength element (HFSE) signatures, which have been interpreted to reflect different conditions of melting (e.g., temperature and depth) and led to the well-known F classification of Archean rhyolite (Lesher et al., 1986; Hart et al., 2004), with felsic lavas as expressions of the high heat-flow environment. In developing the F classification, a causal link was established between the cooling history of the subvolcanic magma chambers (at least partly driven by hydrothermal convection) and the formation of VMS deposits. However, attempts to assign the geochemical signatures of the felsic volcanic rocks to different geodynamic settings have had mixed success (Prior et al., 1999; Yang and Scott, 2003; Huston et al., 2010; Piercey, 2011; Yeats et al., 2017). It has been shown that REE signatures of felsic volcanic rocks differ significantly between Archean and post-Archean systems, complicating any comparisons (Lentz, 1998; Hart et al., 2004; Piercey, 2011). Early studies that compared Archean rhyolites with samples from modern oceanic settings were limited by the small numbers of samples available from the present-day ocean floor. However, in the last decades ocean exploration and seabed sampling have significantly improved. An important finding has been that felsic volcanic rocks are common in modern oceanic assemblages, and they have a remarkably large compositional range that reflects highly variable processes, including assimilation of hydrated oceanic crust and fractional crystallization in a range of geodynamic settings (Brophy, 2008; Tamura et al., 2009; van der Zander et al., 2010; Wanless et al., 2010; Piercey, 2011; Freund et al., 2013; Embley and Rubin, 2018). With increased sampling, the different conditions of melt generation can now be comprehensively examined and compared to ancient analogues.

In this paper, we compiled data for more than 2,200 unique samples of submarine felsic and intermediate lavas from 71 different locations at mid-ocean ridges, ridge-hot-spot intersections, oceanic back-arc spreading centers, rifted island arcs, and intracontinental back-arc rifts (Fig. 1). Using a combination of principal components analysis, unsupervised hierarchical clustering, and supervised random forest classification, we show consistent geochemical diversity between the different settings. These differences are compared to the compositions of felsic volcanic rocks from the Archean Abitibi greenstone belt in the Superior province of Canada. Although it is clear that the geochemical characteristics of rhyolites have changed because of the evolution of tectonic styles (e.g., Huston et al., 2019), the diversity of modern oceanic felsic volcanic rocks highlights a number of different processes that could explain the signatures of similar rocks in ancient greenstone belts.

In modern oceanic settings, felsic volcanic rocks have been sampled from (1) propagating tips of mid-ocean ridge spreading centers (Wanless et al., 2010; Freund et al., 2013), (2) central volcanic complexes above mantle hot spots (Zellmer et al., 2008; Flude et al., 2010; van der Zander et al., 2010), (3) volcanic fronts of active island arcs (Hunter, 1998; Leat et al., 2007), and (4) arc-related rifts and back-arc spreading centers (Fretzdorff et al., 2006; Tamura et al., 2009) (Fig. 1; Table 1). The complex chemistry and petrogenesis of the felsic magmas closely match the diversity of tectonic settings in which these rocks occur (Tables 1, 2).

Mid-ocean ridge assemblages

Felsic lavas are now widely reported from mid-ocean ridges, including the East Pacific Rise, Galapagos spreading center, Juan de Fuca Ridge, and the Pacific-Antarctic Rise (Wanless et al., 2010; Freund et al., 2013). The silicic lavas are associated with the propagating tips of overlapping spreading centers (OSCs: Perfit and Fornari, 1983), with ridge-hot-spot intersections (Haase et al., 2005), and with ridge-transform intersections (Regelous et al., 1999). The melts are interpreted to be delivered as episodic injections caused by rift propagation or ridge-transform interactions, creating isolated magma chambers that undergo fractional crystallization. Such melts are found on the eastern limb of the overlapping spreading center at 9° N on the East Pacific Rise and at the intersection of Juan de Fuca Ridge with the Blanco transform zone (Wanless et al., 2010). Although the felsic lavas are erupted from small, isolated magma chambers, they are commonly found at or adjacent to axial highs where the largest volumes of melt reside in the crust (Haase et al., 2005). Assimilation of hydrothermally altered oceanic crust is also common at the propagating tips of the spreading centers (Wanless et al., 2010; Freund et al., 2013), further enhancing melting. Where influenced by a nearby hot spot or mantle plume, as at the Foundation hot spot near the Pacific-Antarctic Rise, the hotter mantle leads to increased melting, episodic magma replenishment, and a thicker crust. Under these conditions, more evolved melts develop in multiple magma reservoirs that are influenced by greater wall-rock assimilation and efficiently cooled by hydrothermal convection.

Ocean-island assemblages

Ocean-island dacites and rhyolites, such as those in Iceland, Ascension Island, and the Azores, are products of hot-spot or plume-influenced magmatism (van der Zander et al., 2010; Jicha et al., 2013; Jeffery et al., 2016). The felsic magmas are produced by deep partial melting (van der Zander et al., 2010) or shallow fractional crystallization and differentiation (Jicha et al., 2013). The volcanic products belong to well-developed bimodal successions, such as the large central volcanoes and fissure-dike swarms of Iceland (Jónasson, 2007; Sigmarsson and Steinthórsson, 2007). The Krafla central volcano, situated in the North Iceland rift zone, for example, hosts numerous fissure-dike swarms occupied by doleritic and felsic intrusions (Thordarson and Larsen, 2007). Rhyolite lenses in the rift zone are interpreted to be products of partial melting of hydrothermally altered basalt (Zierenberg et al., 2013) or fractional crystallization with assimilation of crustal material (Brophy, 2009). Although Icelandic magmatism is largely tholeiitic, having originated from melting at low fO2, calculated fH2O values can be as high as those of oxidized island-arc melts (e.g., derived from 90% fractional crystallization: Portnyagin et al., 2012). Similarly, on Ascension Island (Chamberlain et al., 2016) and at the Marquesas Archipelago (Chauvel et al., 2012), closed-system fractionation in shallow magma chambers leads to evolved melts with high H2O concentrations. Cyclic eruptions of evolved melts at the Furnas volcano in the Azores have been interpreted to reflect a stratified magma chamber (Jeffery et al., 2016). Shallow differentiation and deep hydrous melting produce tholeiitic and calc-alkaline rhyolites in the late, shield-building stages of the Waianae volcano, Hawaii (van der Zander et al., 2010), whereas felsic volcanic rocks associated with the shield-building stage in Samoa are alkaline (Konter and Jackson, 2012). The latter are similar to partial melts of deep-seated basaltic dike swarms at Huahine in the Society Islands (Harnois and Stevenson, 2006).

Intraoceanic arc and back-arc assemblages

Felsic volcanic rocks are common in modern intraoceanic arc and back-arc assemblages. Typically, the arc front is intruded by fluid-fluxed melt, whereas the back-arc basin is intruded by mantle rising into the back-arc spreading center in response to basin opening. At the volcanic front, devolatilization of the subducted slab enriches the mantle wedge in H2O, resulting in abundant calc-alkaline melts (Grove et al., 2003; Tatsumi and Suzuki, 2009). These commonly are stalled in the upper arc crust and fractionate (Tamura et al., 2009; Deering et al., 2011). The erupted melts have variable chemistry, reflecting input of different slab-derived components and different melting and crystallization processes (Tamura et al., 2009; Haraguchi et al., 2017). For example, on the Izu-Bonin arc, intra-arc volcanism is due to remelting of the basaltic or andesitic crust and fractionation in high-level magma chambers, producing a range of bimodal rift-related volcanic rocks (e.g., in the Sumisu rift) and rhyolite-dominated submarine volcanoes (e.g., Myojin Knoll: Shukuno et al., 2006; Tani et al., 2008; Tamura et al., 2009). Seismic profiles of the Izu-Bonin arc have confirmed that rhyolite-dominated volcanoes at the arc front are situated above thinned middle crust where partial melting occurred (Tamura et al., 2009). Bimodal volcanism at the front of the Tonga-Kermadec arc is related mostly to basalt fractionation with intrusions of mafic melts into the middle crust where they stagnate, fractionate, and eventually rise to the surface (e.g., Raoul, Macauley, and Healy volcanoes: Turner et al., 2012; Barker et al., 2013). The sources of fluids that induce melting are variable. Volcanic rocks of the southern Kermadec arc show the influence of fluids derived from subducted sediments (e.g., Brothers volcano), whereas in the northern Kermadec arc, the melts have been influenced by fluids from altered oceanic crust (e.g., Healy: Haase et al., 2006). Rhyolites of the northern Tonga arc volcanoes (e.g., Fonualei volcano, Niuatoputapu volcano) and nearby rear-arc complexes (Niuatahi volcano) are mainly related to intra-arc rifting and basalt fractionation (Turner et al., 2012; Park et al., 2015). Felsic rocks of the Proctor shoal in the South Sandwich arc are products of partial melting of basaltic crust within the arc (Leat et al., 2007). Felsic melts of the Anatahan volcano in the Mariana arc are derived solely from fractional crystallization without any crustal contamination (Wade et al., 2005). In contrast, felsic rocks from the Soufriere Hills in the Lesser Antilles are products of complex multistage magmatism from partial melting of the mantle wedge, deep fractional crystallization, and magma mixing and assimilation of crustal material (Zellmer et al., 2003).

Silicic magmatism in back-arc basins is focused at propagating rifts, overlapping spreading centers, and off-axis volcanoes (Fretzdorff et al., 2006; Beier et al., 2015). Active spreading centers in mature back-arc basin settings closely resemble mid-ocean ridges, and the felsic lavas are mostly related to fractional crystallization in isolated magma chambers, in some cases with local assimilation of hydrated crust. Enhanced cooling and fractionation of the melt lenses is promoted by stagnation in older and colder crust behind the propagating rift tips (Fretzdorff et al., 2006). This may be particularly common in areas of greater crustal thickness, as in the Manus basin (e.g., Beier et al., 2015). Back-arc spreading centers that are located farthest from the arc are most similar to mid-ocean ridges, with the least influence of subduction-derived components (e.g., no large ion lithophile element [LILE] enrichment or negative Nb-Ta anomalies: see below). Felsic volcanic rocks in these settings have been sampled at the Central Lau spreading center and Northwestern Lau spreading center in the Lau basin and at the Manus spreading center in the Manus basin (Jenner et al., 2012; Lytle et al., 2012; Beier et al., 2015). In contrast, back-arc spreading centers located immediately behind the arc have significant subduction-derived components, as at the Valu Fa Ridge behind the Tofua arc and at Pual Ridge and SuSu Knolls in the Eastern Manus rifts (Sinton et al., 2003; Zhang et al., 2018).

Intracontinental arc and back-arc assemblages

Felsic volcanic rocks are abundant in the early stages of rifting of continental crust, such as in the Taupo Volcanic Zone. Magmatism in the Taupo Volcanic Zone, both onshore and offshore, is linked to faulting, diking, uplift, and localized subsidence where the crust is most thinned (Parson and Wright, 1996; Rowland et al., 2010). The transition from andesitic to dominantly rhyolitic volcanism is associated with the onset of strong extension, allowing the felsic melts to rise into the upper crust (Deering et al., 2011). Felsic volcanic rocks in the Okinawa trough back-arc basin are similarly associated with melting of continental crust that is passively rifting with anomalously high heat flow (Shinjo and Kato, 2000). In both settings the melts show evidence of assimilation of crustal material and mixing of basaltic and more evolved magma (Shinjo and Kato, 2000; Deering et al., 2011). In contrast, felsic melts of the Bransfield Strait, an extending marginal basin of the West Antarctic Peninsula, are produced by fractional crystallization with minor crustal assimilation (Fretzdorff et al., 2004).

Geochemical analyses of silicic lavas (>60 wt % SiO2) were compiled for each of the locations described above. From an original database of more than 2,200 published analyses, we removed any samples analyzed by inappropriate methods (e.g., X-ray fluorescence for trace elements), duplicate samples, results of duplicate analyses, and samples with location errors or incomplete descriptions. We also rejected samples with incomplete major and trace element data, so all results contain at least Si, Ti, Al, K, Y, Zr, Nb, Ba, La, Sm, Eu, Yb, and Th. These elements were chosen based on their availability in the published data sets and because of their common use in lithotectonic discrimination diagrams. Every sample in the final data set was analyzed for each of these elements. No samples were included where the concentration of any element was equal to or below detection or exceeded the upper limits of detection or where the limits of detection were not specified. We further inspected the data by comparison with fresh-rock compositions of dacite and rhyolite, including inspection of REE profiles and immobile elements. We chose not to exclude potentially mobile elements from the analyzed data, such as K, Si, Ba and La, because they are important discriminators of rock type in unaltered samples. Pearce element ratio (PER) plots of these data (App. Fig. A1) show that the modern samples are not significantly affected by alteration. Figure 2 shows the range of felsic volcanic rock compositions in several common major and trace element discrimination diagrams that have been used for both basalt and rhyolite: a total alkali-silica plot (after Le Maitre, 1989), a Winchester and Floyd-type immobile element plot (after Pearce, 1996), and a magmatic affinity plot (after Ross and Bédard, 2009). The final database includes 404 high-quality analyses from modern oceanic settings (App. Table A1) plus a comprehensive data set of 239 samples from the well-studied Taupo Volcanic Zone.

The mid-ocean ridge felsic volcanic rocks from the East Pacific Rise, Galapagos spreading center, Pacific-Antarctic Rise, and Juan de Fuca Ridge, include andesite, dacite, and rhyolite (Fig. 2A, B). Samples from ridge-hot-spot intersections and ocean islands, such as Iceland, the Azores, and Ascension Island, are mostly trachytes and low-alkali dacites and rhyolites. Samples from intraoceanic back-arc assemblages, such as the Lau basin, are mainly andesites and dacites and plot in a strong linear array from low to high SiO2 with increasing Na2O + K2O. Samples from the Manus basin have a similarly large range of compositions from andesite to rhyolite. Samples from island-arc assemblages, such as Izu-Bonin, range from dacite to rhyolite and have a narrower range of Na2O + K2O. Felsic volcanic rocks from intracontinental arc-back-arc assemblages, such as the Okinawa trough, Western Antarctic Peninsula, and Taupo Volcanic Zone, are mostly rhyolite and dacite. A number of the arc and back-arc samples plot outside the fields of andesite in Pearce (1996) (Fig. 2B). The original Winchester and Floyd (1977) discrimination diagram is based on a limited number of samples, and therefore the field boundaries for andesite and rhyolite do not reflect the full variation of data available today. Pearce (1996) addressed this using a more comprehensive data set, but increasing numbers of samples dredged from the modern ocean sea floor show that some boundaries may still need to be adjusted. Data for a selection of felsic rocks from the Abitibi greenstone belt (Hart et al., 2004; Gaboury and Pearson, 2008) are shown for comparison; as discussed further below, they are not meant to be classified by these plots.

Despite the diversity of sample sites, the major and trace element compositions of felsic volcanic rocks from different settings are remarkably uniform. The andesites and dacites from mid-ocean ridges have SiO2 concentrations up to 68.0 wt %, low K2O (0.9–1.4 wt %), high TiO2 (0.6–1.9 wt%), flat REE profiles, and high values of Yb(CN) (54.5–99.9; CN = chondritenormalized: Table 1). Felsic volcanic rocks from ocean islands have similar SiO2 (60.3–68.1 wt %), intermediate TiO2 (0.1–1.8 wt %), and high K2O (2.8–6.0 wt %), with slightly to strongly elevated light rare earth elements (LREEs) compared to midocean ridge samples and intermediate Yb(CN) (6.3–37.7). The felsic volcanic rocks from ridge-hot-spot intersections have a large range of SiO2 concentrations (61.1–76.5 wt %), high K2O (1.9–5.9 wt %), and REE concentrations that are slightly enriched with intermediate Yb(CN) (12.9–52.5). Felsic volcanic rocks from intraoceanic back-arc rifts and spreading centers have SiO2 concentrations up to 74.6 wt % and higher K2O (0.3–2.0 wt %) than the felsic rocks from mid-ocean ridges, as well as low TiO2 (0.3–2.1 wt %), flat REE patterns, and generally low Yb(CN) (8.6–79.7). Samples from intraoceanic island arcs have higher SiO2 concentrations (up to 75.8 wt %) than samples from back-arc spreading centers and consistently low Yb(CN) (9.80–24.05). Felsic volcanic rocks from intracontinental back-arc rifts and spreading centers have the highest SiO2 (71.1–80.4 wt %) and K2O (2.1–4.9 wt %), low TiO2 (0.06–0.45 wt %), slightly elevated La/Yb(CN) (4.22–8.64), and mostly low Yb(CN) (6.9–15.2). The trace and REE patterns of specific assemblages are summarized in Figure 3.

The origins of the silicic magmas in these different settings are debated (Haase et al., 2006; Brophy, 2008; Tamura et al., 2009; Wanless et al., 2010; Piercey, 2011). Although the dominant processes are thought to be partial melting of hydrated mafic crust and fractional crystallization in mafic magma chambers, there are important variations. Melting in the mantle to produce the precursor basalt can be influenced by geochemical heterogeneities caused by earlier partial melting, input of plume-derived material, or mantle mixing. Thus, variations in the composition of the mafic crust can have a major influence on the derived silicic magmas. This is monitored using trace elements that track enrichment and depletion of incompatible components (Hofmann, 2014). In subduction settings, melting in the mantle is enhanced by the presence of aqueous fluids released from the dehydrating slab (so-called “wet melting”) (Pearce et al., 2005; Langmuir et al., 2006; Escrig et al., 2012). Fluids from the slab are highly enriched in LILEs, whereas the less mobile HFSEs remain in the slab. Thus, arc basalts derived from wet melting are enriched in LILEs and significantly depleted in moderately incompatible elements such as Ti and heavy rare earth elements (HREEs) (Kelemen et al., 2014). A diagnostic feature is a relative depletion of Nb and Ta, which are insoluble in aqueous fluids and are retained in the mantle wedge above the dehydrating slab (Briqueu et al., 1984; Baier et al., 2008). Others have argued that this is occurring in the slab where rutile and other phases retain Nb, Ta, and Ti (e.g., Ryerson and Watson, 1987; Rudnick et al., 2000). This is reflected in a negative Nb-Ta anomaly in normalized trace element plots of the basalts and any derived felsic volcanic rocks (see below). The high water contents of arc-related melts also suppress plagioclase fractionation, which results in less pronounced negative Eu anomalies, lower Ti and Fe concentrations, and higher Al concentrations compared to anhydrous melts (Danyushevsky, 2001; Tamura et al., 2005; Langmuir et al., 2006). The high water content also results in high oxygen fugacity (fO2), suppressing sulfide saturation in the melts (Jenner et al., 2010).

Dry melting of the mantle occurs at mid-ocean ridges and at mature back-arc spreading centers in response to adiabatic decompression (Langmuir et al., 2006). The felsic melts generally lack the negative Nb-Ta anomaly and show no LILE enrichment typical of wet melting, and they have pronounced negative Eu anomalies. When mafic melts rise into axial magma reservoirs and become isolated and fractionate, any differentiated products will similarly lack the Nb-Ta anomaly and LILE enrichment. Crustal thickness determines the temperature and pressure and therefore the stability of different phases (e.g., garnet or amphibole) during melting. High pressures associated with thicker crust stabilize garnet, which concentrates HREEs (Hart et al., 2004). At lower pressures, amphibole is the residual phase resulting in the enrichment of La over Yb in the melt and lesser HREE (Yb) depletion. Melting at shallow crustal levels in the absence of garnet or amphibole fractionation results in flat REE patterns in the associated rocks (Hart et al., 2004).

A third process—assimilation of hydrothermally altered crust—can result in significant enrichment of HREEs, owing to the presence of HREE-rich hydrothermal amphibole in the melted material (Wanless et al., 2010). This is important at mid-ocean ridges but appears to be less so in back-arc assemblages (Beier et al., 2015).

From the above, the characteristic features of silicic lavas at mid-ocean ridges (i.e., negative Eu anomalies, significant HREE enrichment, and flat REE patterns) may be related to (1) a lack of garnet fractionation during shallow/low-pressure melting causing flat REE patterns, (2) plagioclase fractionation under dry melting conditions leading to depletion of Eu, and (3) partial melting and assimilation of altered crust containing HREE-rich hydrothermal amphibole. In contrast, in subduction settings, the felsic rocks have weak negative Eu anomalies, LILE enrichment, and strong Nb-Ta anomalies related to (1) the suppression of plagioclase fractionation under wet melting conditions, (2) enrichment of LILE by fluids from the dehydrating slab, and (3) the low solubility of Nb and Ta in the subduction fluids. In continental margin arcs, such as the Okinawa trough and Taupo Volcanic Zone, the melt compositions additionally reflect variable assimilation of crustal material with high LILEs, high Th, low HFSEs, and low HREE concentrations (Shinjo and Kato, 2000).

Principal component analysis (PCA) and unsupervised hierarchical clustering

Multivariate statistical analysis was undertaken to determine whether the modern felsic volcanic rocks from different tectonic assemblages can be clearly distinguished in the whole-rock geochemical data. The analysis was performed with the R statistical software package (R Core Team, 2020). The workflow is presented in Appendix Figure A2 and can be replicated with the code and instructions provided in the Appendix Text and Workflow. All data were recalculated to elemental concentrations in parts per million, and centeredlog ratio transformations were applied to overcome closure effects (Aitchison, 1982; see App. Table A2). PCA (Bartholomew, 2010) was conducted with all 13 elements (Si, Ti, Al, K, Y, Zr, Nb, Ba, La, Sm, Eu, Yb, Th) for all samples in the data set. PCA reduces the variables into a limited number of components that account for the largest variance in the data set. Each principal component (PC) is affected by all the variables and represents a fraction of the total variance; the first and second principal components, PC1 and PC2, account for most of the variance (54 and 31%, respectively; see App. Table A3). The loadings of each element on each component describe how the different elements contribute individually to the variance (Fig. 4; App. Table A4). In Figure 4, the red element labels correspond to loadings for the entire data set. The plotted points for each sample are a linear combination of the loadings for Si, Ti, Al, K, Y, Zr, Nb, Ba, La, Sm, Eu, Yb, and Th plotted in PC1-PC2 space. The loadings for each sample are given in Appendix Table A5. The analysis shows Yb and Y have strong positive loadings on PC1, whereas Th, La, and K have strong negative loadings. Zirconium, Nb, and Sm have strong positive loadings on PC2, whereas Si, Al, and Ba have strong negative loadings (Fig. 4; Table 3). The PCA shows very consistent behavior between elements that are commonly more mobile (e.g., K, Ba, and La) and the immobile elements (Th and Al), suggesting that the results are reflecting primary melt compositions and not alteration. To determine if alteration of the samples has affected the results, we repeated the PCA using only Al, Nb, Ti, Th, Y, Yb, and Zr, which yielded nearly identical PC1 and PC2 scores (see App. Fig. A3).

We used an unsupervised agglomerative hierarchical clustering of the PCA data to identify geochemically distinct groups of samples in the data set. This type of cluster analysis applies a set of queries or algorithms to classify samples into different groups and maximize the differences between the groups. This is done using Ward’s criterion, which is the sum of the squared distances of individuals from the center of the cluster. In our study, we performed a cluster analysis on the first three principal components from the PCA. They account for ~90% of the variance and the most important geochemical features in the data set. The cluster analysis was performed with the R statistical software using the unsupervised hierarchical clustering function HCPC (Husson et al., 2009). In this type of analysis, no preassigned groups or clusters are identified (i.e., the analysis is unsupervised), and the number of clusters is determined by a hierarchical procedure that builds the data tree (i.e., a dendrogram of decision steps that identify the different clusters). Agglomerative hierarchical clustering starts by assuming that every individual sample forms its own cluster, and then the statistical differences (or “distances”) from all other samples are calculated (e.g., using a proximity matrix). The two samples with the closest distance to each other form the next level of clusters. Ward’s criterion calculates the center of each cluster and the distance between clusters. The smallest number of clusters representing all of the samples was determined using an iterative approach. Figure 4A shows nine clusters (C1-C9) that were identified with clear separation of the samples based on the loadings of the different elements on PC1 and PC2. The elements with the greatest influence in identifying the individual clusters are Nb, La, and Zr (C1); Nb, La, and Th (C2); Y, Yb, and Sm (C3); Ti, Eu, and Al (C4); Si, Al, and Th (C5); Si, Al, and Ba (C6); Yb, Zr, and Sm (C7); Th, La, and Ba (C8); and K, Ba, and Th (C9). The clustering results and the cluster metadata are given in Appendix Tables A5 and A6.

Felsic volcanic rocks from ridge-hot-spot intersections cluster in the field with negative loadings on PC1 and positive loadings on PC2 (C1 and C2), reflecting the incompatible element-rich mantle source and high Nb and La of oceanic hot-spot assemblages. Samples from mid-ocean ridges cluster in the field with positive loadings on PC1 and positive loadings on PC2 (C3), reflecting the significant enrichment in HREEs due to assimilation of hydrated oceanic crust containing HREE-rich hydrothermal amphibole. Felsic volcanic rocks from intraoceanic arc-back-arc environments cluster in three distinct fields with small to large positive loadings on PC1 and slightly negative loadings on PC2 (C4-C6), reflecting strong enrichment of elements that are mobile in slab-derived fluids, such as Ba, and depletion of insoluble elements, such as Zr. A few samples from intraoceanic back-arc environments, such as the Manus spreading center in the Manus basin and the Northwest Lau and Central Lau spreading centers in the Lau basin, plot in the same cluster as mid-ocean ridge samples (C3). This reflects significant HREE enrichment due to assimilation of mid-ocean ridge-like crust. Felsic volcanic rocks from intracontinental arc-back-arc systems, including the Taupo Volcanic Zone, mostly plot with a narrow range of positive to negative loadings on PC2 and a large range of positive to negative loadings on PC1 (C7-C9), reflecting strong enrichment of elements in slab-derived fluids (Ba, K) and enrichment of Th due to crustal assimilation. Felsic volcanic rocks from the marginal back-arc basin of the West Antarctic Peninsula cluster in a field with positive loadings on PC1 and PC2 (C7), reflecting less enrichment of fluid-mobile elements (Ba, K) due to stalled subduction (see below).

The characteristics of each cluster and their relationship to specific lithotectonic assemblages are summarized in Table 3 and Figure 4A (see also App. Fig. A4). Cluster C1 contains mainly samples from the Society Island group, Marquesas Archipelago, and the Azores, whereas cluster C2 contains mostly samples from Iceland and Ascension Island. Samples in C1 have higher Nb, Th, and Zr concentrations and stronger pronounced negative Eu anomalies. The samples from the Society Islands are related to partial melting of deep-seated basalt with significant crustal contamination (Harnois and Stevenson, 2006), whereas the samples from the Azores are related to fractional crystallization of cyclically rejuvenated magma chambers with trachytic caps (e.g., at the Furnas volcano located at the Terceira rift: Jeffery et al., 2016). The samples from Iceland in C2 have a transitional affinity, high La/Yb ratios, and negative Eu anomalies related to the mantle plume and increased crustal thickness, which causes more alkaline geochemical affinity. These rocks are products of near-solidus differentiation (Jónasson, 2007), partial melting of hydrothermally altered and metasomatized basalt, and further fractional crystallization (Zellmer et al., 2008). Although Ascension Island is located 90 km west of the Mid-Atlantic Ridge, the felsic lavas erupted when Ascension was closer to, or at, the ridge axis. These samples occur in both C1 and C2, reflecting two different fractionation trends: one with increased Zr/Nb ratios due to ilmenite fractionation (i.e., partitioning Nb) and the other related to zircon/titanite fractionation (i.e., partitioning Zr: Jicha et al., 2013). Cluster C3 contains all of the mid-ocean ridge dacites from the Pacific-Antarctic Rise, East Pacific Rise, Galapagos spreading center, and Juan de Fuca Ridge. These samples have anomalously high HREE concentrations, typical for assimilation of hydrothermally altered crust at mid-ocean ridge spreading centers. Cluster C3 also contains a number of samples from intraoceanic back-arc spreading centers, notably the Manus spreading center and the Northwestern Lau and Central Lau spreading centers (Table 3). These samples have flat REE patterns, mostly lack negative Nb-Ta anomalies, and show no LILE enrichment normally associated with subduction-related assemblages. Instead, they are located far from a subduction zone where the melts are influenced mainly by the mantle source and spreading regime. One sample from the West Antarctic Peninsula (Bransfield Strait, which is a failed intracontinental back-arc spreading center behind the South Shetland island arc), plots in the same cluster as mid-ocean ridge samples (C3), reflecting no enrichment of fluid-mobile elements (Ba, K) (see App. Fig. A4), despite its proximity to the continental margin. The subduction zone at Bransfield has stalled (Fretzdorf et al., 2004), resulting in a decrease in the contribution from a dehydrating slab and a rhyolite signature closer to that of the East Pacific Rise and Pacific-Antarctic Rise. Cluster C4 contains felsic volcanic rocks from the Raoul volcano of the Tonga-Kermadec arc and Sumisu Knoll at the volcanic front of the Izu-Bonin arc. These samples have strong arc signatures (e.g., deep negative Nb-Ta anomalies, LILE enrichment, and weak negative Eu anomalies) and flat REE patterns, and they are less enriched in HREEs than rhyolites from mid-ocean ridge or most back-arc spreading centers. An exception is the Valu Fa rift and back-arc spreading center of the Southern Lau basin, which is propagating into arc crust of the Tofua volcanic arc. Samples from the Valu Fa rift fall in cluster C4 (see App. Fig. A4) and share a strong arc signature resulting from fractional crystallization of mafic magma chambers with the addition of components from the nearby arc (Fretzdorff et al., 2006; Barker et al., 2013). The volcanic rocks at Sumisu Knoll are thought to be products of partial melting of the middle arc crust (Tamura et al., 2009; Haraguchi et al., 2017). Cluster C5 contains felsic volcanic rocks from the volcanic front of the Mariana arc (Anatahan), the Lesser Antilles arc (Montserrat), the Kermadec arc (Brothers and Healy volcanoes), and the Niuatahi rear-arc volcanic complex behind the active Tofua arc in the northern Lau basin. Samples from arc-related rifts are also included in this cluster, namely from the Sumisu rift and Myojin rift behind the Izu arc, two samples from the Bransfield Strait, and several samples from the Pual Ridge in the Manus basin. All samples have strong arc signatures, including LILE enrichment, slight LREE enrichment, and weak negative Eu anomalies related to fractional crystallization and/or partial melting of the middle arc crust. They exhibit a negative Nb-Ta anomaly, but it is less pronounced than the samples belonging to cluster C4. Two samples from the Okinawa trough also fall into cluster C5. Cluster C6 contains felsic volcanic rocks from the Fonualei and Niuatoputapu volcanic islands on the Tofua arc, Macauley volcano on the Kermadec arc, and Proctor shoal on the South Sandwich arc. These samples have strong arc-like signatures similar to the samples in cluster C5 but with a slightly weaker negative Eu anomalies, stronger pronounced negative Nb-Ta anomalies, and generally lower HFSE and REE concentrations. They are related to mafic melt lenses that were emplaced into the cold arc crust and underwent fractional crystallization (Turner et al., 2012; Barker et al., 2013). Cluster C6 also contains samples from the Bayonnaise Knoll, located in a back-arc depression behind the Izu volcanic arc, from rift-related assemblages at SuSu Knolls in the Manus basin, and from the Torishima rift and Aogashima rift behind the Izu volcanic arc. These rocks, including samples from Pual Ridge in the Manus basin, are also products of fractional crystallization of mafic melt lenses with minor contributions from assimilated crustal material or partially melted middle arc crust (Tamura et al., 2009; Beier et al., 2015). Samples from some locations (e.g., Macauley volcano, Torishima rift, and Pual Ridge) plot both in cluster C5 and C6, reflecting a transitional character. Cluster C7 is mainly defined by the felsic volcanic rocks from Hook Ridge in Bransfield Strait. However, it also contains four samples from backarc rifts behind the Izu-Bonin arc (Sumisu rift and Myojin rift) and two samples from the Mangatolu triple junction within the northern Lau back-arc basin. These rocks all have slight negative Nb-Ta anomalies, slight LILE enrichment, intermediate negative Eu anomalies, mostly flat REE patterns, and slightly less HREE enrichment than mid-ocean ridge samples, reflecting limited input from the subduction zone compared to other arc-related samples (C4-C6). Clusters C8 and C9 contain mostly felsic rocks from the Okinawa trough back-arc rift and the Taupo Volcanic Zone, which are enriched in LREEs, alkali elements, and Th and have negative Nb-Ta anomalies. However, samples from the Taupo Volcanic Zone fall into two clusters, capturing both dry (i.e., cluster C8) and wet (i.e., cluster C9) melting conditions, as previously documented by Deering et al. (2008) (see discussion). A few samples from hotspot volcanoes (Hawaii, Samoa, and Ascension Island) plot in C8 (Fig. 4A) because of their low HREE concentrations and stronger alkaline signatures compared to other hot-spot volcanoes. These samples are related to late-stage shield volcano development and deep hydrous crustal melting (van der Zander et al., 2010).

The cluster analysis illustrates the significant geodynamic control on the geochemical signatures of modern oceanic felsic volcanic rocks. Clusters C1 and C2 correspond mainly to hot-spot-related assemblages, C3 corresponds to mid-ocean ridge-like assemblages, C4 to C6 correspond to intraoceanic arc assemblages, C7 corresponds to arc-back-arc assemblages with weak subduction zone influence, and C8 and C9 correspond to intracontinental arc assemblages. In the following section we perform a supervised classification based on the identified clusters and conduct a blind test to determine if the classification scheme can accurately predict geodynamic influences in the felsic volcanic rock geochemistry.

Training a new classifier

Supervised classification is a machine learning task that sorts samples into known classes identified using a training data set. The training is an iterative process by which a machine learning classifier, such as random forest, learns what values are meaningful to achieve the best classification outcome and builds a model that can be used as a classification scheme. We used our database of 643 samples of modern felsic volcanic rocks to build the model. 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 initial input is a selection of parameters (i.e., analyzed elements) and predetermined clusters (i.e., target classes) for each sample. These are randomly divided into a training data set, which is used to build and validate the random forest classifier, and a blind test data set, which is used to evaluate the prediction rate. Random forest has proven to be particularly useful for supervised classification based on multivariate input (Fernandez-Delgado et al., 2014). It functions using multiple decision trees with true or false outcomes (Breiman, 2001) (e.g., “Are SiO2 concentrations >73 wt % ?”; if yes, the sample is a rhyolite; if not, the sample is a dacite). These decision points are commonly referred to as “nodes”; follow-on decision points are referred to as “child nodes.” The cutoffs for each cluster are determined by the classifier’s ability to separate the training data most efficiently. For example, if dacites and rhyolites are to be divided at a node based on their SiO2 concentrations, a value of 73 wt % correctly identifies all dacites and all rhyolites. The correct cutoff is identified through an iterative approach testing several SiO2 concentrations until the best separation of the clusters is achieved, commonly referred to as “pure” separation. At every node, the “impurity” can be measured, corresponding to the probability of a parameter falsely classifying a sample, and therefore a measure of the importance of the parameter in the decision making. 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 (e.g., C1), j refers to the possible classes (e.g., C1-C9), gc is the relative frequency of classification, nc 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.

A single decision tree alone may be very sensitive to overfitting; it 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 uncorrelated decision trees operating as a “committee” that casts votes for each decision. Multiple decision trees (usually a thousand or more) are created by randomly selecting elements (i.e., the nodes) from a training data set. Each decision tree is then evaluated or “trained” with randomly chosen samples (equal to the total number of samples in the data base), a method commonly known as “bootstrapping.” Importantly, during this process samples can be picked multiple times, which leaves behind a pool of samples that were not picked and therefore not involved in some decision trees (so-called “out of bag” samples), typically about one-third of all samples. These out of bag samples are used to validate the classification (see below). To classify an unknown sample, each of the randomized decision trees casts a single vote for a single target class or cluster (C1-C9). The votes cast by all of the decision trees are given in percentages (summing to 100%), and random forest uses the majority decision of these votes to finally classify a sample (Breiman, 2001). A disadvantage of the random forest 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 do not receive a majority vote but rather a randomly distributed set of votes that can be investigated in the metadata.

We randomly selected 80% of our database to train and validate our random forest classifier; the remaining 20% was used for a blind test (see below), which is a standard approach to validate the random forest classifier (Gholamy et al., 2018). The target classes are the nine clusters shown in Figure 4. The differences in the number of samples in the target clusters can create a class imbalance, which can lead to poor performance of the random forest classifier. To overcome this, we equalized the number of samples for each cluster prior to the training process using the synthetic minority oversampling technique (SMOTE: Chawla et al., 2002). SMOTE creates new samples by slightly altering the geochemical composition of original samples until each class contains the same number of samples as the largest class, resulting in an artificial training set with 1,440 samples. The data created by this process are given in Appendix Table A7.

Building and testing our classifier was done using the R statistical software, in which a defined number of decision trees and a minimum number of nodes are selected, based on classification performance. Gini indexes were calculated to evaluate the performance of each node, which were ordered from their highest to their lowest predictive power. One thousand decision trees, each trained with three randomly selected elements yielded the highest predictive power in our classifier. To validate the model, we used the subset of samples remaining from the bootstrapping procedure (i.e., out of bag samples). These samples were separated from the original 643 samples, and none were created via SMOTE. The method correctly classified 100% of the out of bag samples of all target classes C1-C9 (App. Table A8). The impurity for each node averaged over the complete random forest classification scheme is referred to as the mean decrease Gini. It corresponds to the probability of misclassification and identifies Al, Ba, Nb, Th, Yb, and Y as the most important variables in determining whether a sample is classified correctly (Fig. 5). In contrast, nodes corresponding to Ti and Eu have low mean decrease Gini values and therefore less influence in the correct identification of a sample. We then conducted a blind test on the randomly selected samples that were not used in the training of the model (20% of our original database) to determine how well the classifier performs on new data. At least three samples were available for each target class in the blind test (see App. Table A9). The test samples were correctly classified in almost every case for all target classes, reaching an overall accuracy of 96% (App. Table A10).

Using automated PCA and hierarchical clustering, we have identified nine target classes of felsic volcanic rock compositions in the modern oceans corresponding to distinct geodynamic settings and reflecting different melting conditions. Based on these target classes, we then constructed a random forest classifier that accurately identifies rocks from those settings. Here, we discuss the petrogenetic significance of the target classes and then assess the application of the random forest classifier to ancient volcanic rocks.

Petrogenetic significance of the hierarchical clusters

The most important geochemical characteristics of the target classes in our model include high La/Sm in target class C1 and C2, indicating enriched mantle contributions to the source melts; flat and HREE-rich REE patterns in target class C3, indicating assimilation of hydrothermally altered crust; negative Nb-Ta anomalies and LILE enrichments with low Th in target classes C4, C5, and C6, related to subduction-zone processes in intraoceanic arc-back-arc systems; and similar negative Nb-Ta anomalies and LILE enrichments but with high Th in target classes C7, C8, and C9, reflecting contributions from continental crust. These characteristics correspond closely to the established petrogenesis of felsic volcanic rocks in different geodynamic settings described above and summarized in Table 3. However, no prior knowledge of the geodynamic settings of the samples has influenced the clustering described here.

The characteristics of felsic rocks in target class C1 are most consistent with an ocean island setting, with increased heat flow and enriched mantle sources typical of the hot-spot environment (van der Zander et al., 2010; Chauvel et al., 2012; Konter and Jackson, 2012). In addition to high La/Sm, they are characterized by intermediate to strongly negative Eu anomalies, no negative Nb-Ta anomaly, and an alkaline affinity. Low Ti, Si, and Al and high Nb and La result in negative loadings on PC1 and positive loadings on PC2. Felsic rocks in target class C2 also have enriched ocean island basalt (OIB)-like mantle signatures, high La/Sm, intermediate to strongly negative Eu anomalies, and no negative Nb-Ta anomalies, typical of ridge-hot-spot intersections, but they have lower Nb, Th, and Al and more variable and higher Si and are less alkaline than rocks from target class C1. This reflects shallower melting and a shift to tholeiitic affinity in rift regimes. Felsic rocks in target class C3 have significant enrichment in HFSEs with high total REEs and flat REE patterns but no negative Nb-Ta anomalies and low potassium, typical of rocks from mid-ocean ridges. The pronounced negative Eu anomaly and lack of LILE enrichment reflect dry melting conditions. The characteristics of felsic volcanic rocks in this class overlap with felsic rocks from intraoceanic back-arc spreading centers far from any subduction zone or where subduction has stalled, resulting in no dehydration of the downgoing slab. Felsic volcanic rocks in target classes C4, C5, and C6 have geochemical signatures reflecting a range of melting regimes in intraoceanic arc and back-arc systems. Felsic volcanic rocks in C4 and C6 have strong arc signatures with significant input from a dehydrating slab, as in arc-front volcanoes and arcrelated rifts. These rocks are characterized by large negative Nb-Ta anomalies, slightly less HREE enrichment than target class C7, and small to intermediate negative Eu anomalies. Strong negative Eu anomalies, flat REE patterns, and high HREEs in target class C4 indicate drier melting conditions and a shallower melting regime than samples in target class C6. Felsic volcanic rocks in target class C5 also have strong arc signatures characterized by negative Nb-Ta anomalies, LREE enrichment, and small to intermediate negative Eu anomalies observed at near-arc volcanoes and in related backarc rifts. They are similar to rocks in target class C6 but with higher REEs, Th, and HFSEs and less pronounced Nb-Ta anomalies. Felsic volcanic rocks in target class C7 have weak arc signatures with minor input from a dehydrating slab indicated by small negative Nb-Ta anomalies, less HREE enrichment than at mid-ocean ridges, and intermediate negative Eu anomalies, typical of felsic rocks in extended continental arc or back-arc crust. Felsic volcanic rocks in target classes C8 and C9 have geochemical signatures mostly reflecting the role of continental crust with variable input from a dehydrating slab, as in intracontinental back-arc rifts and spreading centers. They have a calc-alkaline affinity, negative Nb-Ta anomalies, intermediate negative Eu anomalies, LILE enrichment, and high Th. Target class C8 has higher HREEs, indicating drier melting conditions than target class C9.

The interpretations of the different classes are supported by isotopic data that confirm different conditions and sources of melting (Fretzdorff et al., 2006; Haase et al., 2006; Barker et al., 2013; Beier et al., 2015). In the case of ocean islands (target class C1) the felsic volcanic rocks are products of deep, low-degree mantle melting of anhydrous depleted mantle, whereas at ridge-hot-spot intersections (target class C2), they result from shallow fractionation of melt lenses. The latter have transitional alkaline to tholeiitic geochemical affinity but with plume-like signatures that are less pronounced than in ocean islands because of the additional influence of the spreading regime, including stagnation of melt lenses, assimilation of crustal material, and partial melting of crustal rocks (Brophy, 2009; Zierenberg et al., 2013). At mid-ocean ridge spreading centers (target class C3) the volcanic rocks are products of upwelling asthenospheric mantle, with felsic rocks derived from fractional crystallization of isolated magma chambers and assimilation of hydrated oceanic crust at overlapping spreading centers, propagating rift tips, and transform boundaries. Felsic rocks from intraoceanic arc and back-arc systems (target classes C4-C7) all show some influence of the arc and subducting slab, including high degrees of partial melting of the basaltic crust due to the strong volatile input (C4), volcanoes built on older arc or forearc crust with fractional crystallization occurring in stalled magma chambers (C5), or highly variable influence of slab-derived fluids (C6, C7). Felsic rock in continental margin arcs (target classes C8, C9) all show evidence of crustal melting and the input of metasedimentary material.

At back-arc spreading centers there is commonly a transition from wet to dry melting with distance from the volcanic arc (Fretzdorff et al., 2004; Langmuir et al., 2006; Escrig et al., 2012). Back-arc basins that are characterized by diffuse spreading and multiple spreading centers, such as the Northeast Lau basin (Stewart et al., 2022), experience a wide range of wet and dry melting conditions (Zhang et al., 2018), whereas volcanic rocks from back arcs dominated by a single spreading center are geochemically less diverse (e.g., Mariana back arc: Taylor and Martinez, 2003). These differences are illustrated schematically in Figure 6. Thermal convection below spreading centers causes mantle flow orthogonal to the ridge axis (Sparks et al., 1993). Because back-arc spreading centers are commonly asymmetric, rocks from across the basin commonly reflect a range of different melting conditions (Langmuir et al., 2006; Bézos et al., 2009)—the arc-proximal side influenced by wet melts, and the arc-distal side containing drier mantle melts. Differences in the degree of partial melting are indicated by flat HREE-rich and flat HREE-poor end members in close proximity (Tamura et al., 2005; Escrig et al., 2012), as observed in the northern Lau basin (Fig. 7). Felsic volcanic rocks from mature back-arc spreading centers, especially in the west of the Lau basin, lack a negative Nb-Ta anomaly and have overall higher HFSE and REE concentrations and more pronounced negative Eu anomalies than observed in back-arc rifts closer to the arc, such as the Valu Fa rift (Fig. 3D). Identical relationships are seen in the Manus basin, with felsic volcanic rocks from the Eastern rift, located close to the arc, exhibiting strong negative Nb-Ta anomalies, low HSFE and HREE concentrations, and weak negative Eu anomalies compared to spreading centers farther from the arc (e.g., at Pual Ridge; Fig. 3F). Felsic volcanic rocks from the mature Manus spreading center in the center of the basin are nearly identical to those from mid-ocean ridges. These differences correspond closely to the FIV-type arc-like signatures and FIIIb-type intraoceanic back-arc spreading signatures inferred for some Archean rhyolites (see below).

Greater crustal thickness, which causes deeper melting, is reflected in increasing La/Yb(CN), Th, and K, as observed, for example, in the Eastern rifts in the Manus basin (Beier et al., 2015), at the Healy and Brothers volcanoes of the Tonga-Kermadec arc (Barker et al., 2013), and some near-arc volcanoes in thickened crust of the Northeast Lau basin (Niuatahi: Park et al., 2015). These variations are similar to what has been observed in ancient greenstone belts where thickened crust has been proposed (McKenzie and Bickle, 1988).

Comparison with Archean felsic volcanic rocks

On average, Archean crust was more mafic than oceanic crust today (e.g., Palin et al., 2020), and the mafic melts were formed from mantle sources that required melting at very high temperatures (Prior et al., 1999; Piercey, 2011; Thurston, 2015). Felsic volcanic rocks in Archean greenstone belts occur in bimodal mafic-felsic volcanic successions in locally thickened crust. They were derived in part from differentiation of high-level magma chambers and by partial melting of the hydrated oceanic crust, as in the modern oceans (Hart et al., 2004). The volumes of felsic volcanic rocks are significantly greater than in modern intraoceanic settings (e.g., Franklin et al., 2005), probably reflecting the higher mantle temperatures, longer residence times of melt in the crust, and more extensive fractionation. The source of the felsic melts is recorded in part by Zr/Y ratios, HFSE concentrations, and, in particular, REE patterns that underpin the F classification of Archean rhyolites (Table 2; Fig. 8; Lesher et al., 1986; Barrie et al., 1993; Hart et al., 2004).

In the F classification, different types of Archean rhyolites were grouped according to the compatible-incompatible element behavior of Zr, Y, LREEs, and HREEs—the latter reflecting the strong partitioning of HREEs into garnet and amphibole in the source regions of the rhyolitic melts. FI rhyolites have alkaline to calc-alkaline affinities and formed by low-degree partial melting at high pressure (i.e., in the stability field of garnet) with HREE depletion due to garnet fractionation at the source. They are interpreted to have formed from arc-related magmas at deep crustal levels. FII rhyolites have a calc-alkaline affinity and formed from high-degree partial melts with slightly depleted HREEs. They are thought to have formed in continental arc and back-arc settings with melts originating at intermediate crustal depths. FIII rhyolites have tholeiitic affinity and formed by low-pressure partial melting of basalt with variable HREE enrichment (no residual garnet or amphibole). FIIIa rhyolites, which are thought to have formed in shallow melt lenses in rifted island arcs, have flat REE patterns and intermediate HREE concentrations. FIIIb rhyolites, which are interpreted to have formed in oceanic rifts, have flat REE patterns, significant HREE enrichment, and strong negative Eu anomalies caused by partial melting of hydrated mafic crust and plagioclase fractionation. FIV rhyolites have a tholeiitic affinity and are thought to have formed in juvenile volcanic arcs by low-pressure, moderatedegree partial melting of a LREE-poor mantle source. Hart et al. (2004) considered FIV-type rhyolites to be restricted to post-Archean intra-arc assemblages associated with volatilefluxed melting. Figure 9 compares the trace element signatures of modern oceanic felsic volcanic rocks to the fields originally chosen for the F classification of Archean rhyolites. The major differences and similarities are discussed below, with reference to the felsic volcanic rocks of the Abitibi greenstone belt.

The Abitibi greenstone belt is the largest intact Neoarchean greenstone belt and economically one of the most important volcanic terranes in the world. Submarine volcanism produced six temporally distinct assemblages between 2795 and 2695 Ma, traditionally divided into the Northern and Southern volcanic zones (Monecke et al., 2017). They contain a number of mafic and felsic synvolcanic intrusive complexes, such as at Matagami and Noranda, and mafic-ultramafic successions, such as in the Kidd-Munro assemblage. Significant volcanogenic base metal sulfide deposits are associated with the felsic volcanic rocks at Matagami, Noranda, and Kidd Creek. The volcanic and sedimentary rocks of the Abitibi greenstone belt are widely interpreted to be part of an intraoceanic arc-backarc system with similarities to modern oceanic microplate mosaics (Jackson et al., 1994; Percival et al., 2012; Wyman, 2013). However, a number of authors have suggested significant differences between Archean and post-Archean magmatism and tectonics, with some Archean terranes possibly having no modern analogue (see below). We have used the random forest classifier developed in this paper to test whether Archean felsic volcanic rocks of the Abitibi greenstone belt can be classified using similar geochemical features as modern samples and thereby possibly interpreted in terms of processes like those observed in the modern oceans. We make no assumptions about the nature of Archean tectonics or the applicability of Phanerozoic systems to Archean geodynamics, which is a larger discussion beyond the scope of this paper (Bédard et al., 2013; Wyman, 2013).

Felsic volcanic rocks from 20 different locations in the Abitibi greenstone belt (see App. Table A11) were examined with the random forest classifier developed for the modern samples. We used data compiled by Gaboury and Pearson (2008) and the MRD085 data set from the Ontario Geological Survey (2001), which was also used by Hart et al. (2004). These data were selected because the samples were collected specifically to avoid altered rocks, which was confirmed by petrography, and any suspect samples were flagged for removal based on differences between mobile and immobile element behavior. To mitigate any effects of alteration that were not already recognized and filtered out by Gaboury and Pearson (2008), we tested the data using pearce element ratio (PER) plots [e.g., (Ca Na K)/Zr versus (Al Ca)/Zr: App. Fig. A1] and also ran the random forest classifier with and without the mobile elements (see App. Tables A12, A13). The random forest classifier addresses the problem of alteration by failing to classify any samples for which the major or trace elements fall outside identified ranges for unaltered samples in the training sets. The samples included in the analysis are tholeiitic to calc-alkaline felsic volcanic rocks from the Selbaie area, Gemini-Turgeon, Geant-Dormant, Joutel, Normetal, Lac Herbert, Quevillon, and Coniagas in the Northern volcanic zone, as well as samples from the Hunter mine group, Vald’Or, Noranda, Bousquet, the western Blake River Group, and the Timmins-Porcupine area and dominantly tholeiitic rocks from the Matagami area, Kidd Creek, and Kamiskotia. One hundred forty-five samples with >60 wt % SiO2 were selected for the final database, following the same quality assurance/quality control applied to the modern volcanic rocks (App. Table A11).

The dacites and rhyolites from the Abitibi greenstone belt have generally higher SiO2 (60.2–85.9 wt %) than modern submarine felsic volcanic rocks and more variable K2O (0.02–4.43 wt %); Na2O + K2O concentrations are generally lower than in modern felsic volcanic rocks and span the entire total alkali-silica diagram range (Fig. 2A). They have low TiO2 (0.05–0.81wt %), flat to steeply dipping REE patterns, and variable Yb(CN) (9.10–50.10) (Fig. 10). Importantly, negative Nb-Ta anomalies and LILE enrichments, which are the dominant features of felsic volcanic rocks in modern submarine arcs, are found in every assemblage of the Abitibi greenstone belt. The random forest classifier identifies the Abitibi greenstone belt rhyolites as belonging to four main target classes: C2, C3, C5, and C8 (App. Table A12; Fig. 11). However, 87 of the 145 samples were not classified (i.e., no single class received more than 50% of the votes). This indicates either fundamental differences from modern oceanic felsic volcanic rocks or geochemical variation due to alteration and metamorphism among these samples. The majority of samples that were not classified have high and variable K2O that would suggest alteration (e.g., from Kidd Creek: App. Table A11). The remaining 58 samples were assigned to one of the target classes (Fig. 11; App. Table A12). Aluminum, Ba, Nb, Th, and Yb had the greatest influence on the identification of a target class for the Abitibi greenstone belt samples. None of the samples were classified as belonging to target class C1 (e.g., with strong plume influence). However, four rhyolite samples from Kidd Creek were classified as belonging to C2, indicating similarities to ridge-hot-spot intersections, such as Iceland, as previously suggested by Prior et al. (1999). Three samples from Matagami were assigned to target class C3, indicating similarities with modern mid-ocean ridges (FIIIbtype signatures: see below). The samples from Noranda show greater geochemical variability with one sample classified in C3, five samples in C5, one sample in C7, and one sample in C8, similar to rhyolites associated with rifting of thickened arc crust in the Eastern Manus basin, as previously observed by Yang and Scott (2003). A high proportion of felsic volcanic rocks in the Southern volcanic zone of the Abitibi greenstone belt are classified as belonging to intraoceanic back-arc and arc-related rifts (e.g., target class C5). None of the samples are fully classified into C4 or C6, which contain mostly felsic volcanic suites from nascent back-arc rifts and volcanoes at the volcanic fronts of active arcs. Four samples from the Hunter mine area and six samples from the Bousquet Formation were classified into C8 and C9, indicating geochemical similarities with rhyolites from the Taupo Volcanic Zone and Okinawa trough (with FI- to FII-type signatures: see below). Although there is no evidence of continental basement in the Hunter mine and Bousquet areas, the Bousquet Formation is interpreted to have been deposited on much older crust, as indicated by inherited zircon ages (Ayer et al., 2002; Mercier-Langevin et al., 2007). Samples from the Bousquet Formation have negative Nb-Ta anomalies, high Th concentrations, and negative Eu anomalies that are nearly identical to felsic volcanic rocks of the Taupo Volcanic Zone, although the La/Yb and Zr/Y ratios are slightly lower (Fig. 10B). In Appendix Table A13, we show the results for the Abitibi greenstone belt samples using the random forest classifier including only immobile elements, which resulted in more samples being classified into C7 instead of C5.

The comparison with the classification of modern samples suggests that felsic magmas generated above the hot Archean mantle were geochemically most similar to those of modern back-arc spreading centers with mid-ocean ridge-like signatures, such as the Manus, Central Lau, and Northwest Lau spreading centers (target class C3) and the large dacitic lava flows recently documented in the Northeast Lau basin (e.g., Niuatahi in target class C5). The voluminous eruption of dacite at Niuatahi is thought to be due to high upper mantle temperatures at this location (Embley and Rubin, 2018). An important difference is that the Abitibi felsic volcanic rocks have a weaker negative Nb-Ta anomaly than felsic volcanic rocks erupted above modern subduction zones, implying little input from a dehydrating slab or simply a less depleted mantle. In Figure 10A, the geochemical signatures of the felsic volcanic rocks of Matagami and Kidd Creek are compared to Iceland and mid-ocean ridge-like dacites from the northern Lau basin. In these plots, rhyolites from Kidd Creek and Matagami have the greatest similarity to felsic volcanic rocks from mid-ocean ridge-like back-arc spreading centers such as the Northwestern Lau spreading center and Mangatolu triple junction. In contrast, samples from Noranda bear a closer resemblance to felsic volcanic rocks from Niuatahi and Okinawa (Fig. 10D). The close juxtaposition of diverse geochemical signatures in centers such as Noranda (Gaboury and Pearson, 2008) may reflect complex spreading regimes similar to modern, high heat-flow domains in the northern Lau basin and Manus basin (Lagabrielle et al., 1997; Yang and Scott, 2003; Karson, 2017). In particular, the complex microplate tectonics of these settings (Sleeper and Martinez, 2016; Baxter et al., 2020) have been suggested as a possible modern analogue of Archean greenstone terranes, citing the high heat flow and more abundant spreading centers (e.g., Stewart et al., 2022). The melting conditions at Matagami may have been similar to dry melting conditions and/or thickened crust at the Pacific-Antarctic Rise and Northwestern Lau spreading center in the Lau basin today.

Arc-like geochemical signatures are common in the Abitibi (Ayer et al., 2002) but do not prove that volcanic arcs and subduction zones existed, as there could have been a number of different ways to produce the observed signatures (Pearce, 2008; Bédard et al., 2013; Mole et al., 2021). Because Archean rocks are not as depleted as modern volcanic rocks, many argue that modern plate tectonic processes could not have been the cause of negative Nb-Ta anomalies and high Th concentrations (Bédard et al., 2013; Thurston, 2015). A number of alternative models of Archean greenstone belt evolution have been suggested that could account for the differences (e.g., dome-and-keel and mantle wind models: Sawyer et al., 1994; Bédard et al., 2013; Thurston, 2015; Mole et al., 2021). Some have suggested that crustal contamination created the arc-like signatures, which is consistent with the Th enrichment identified in the random forest classification and with evidence of older crust beneath the Abitibi greenstone belt (e.g., Mole et al., 2021). Any process that introduced hydrated basalt to depth could have created conditions appropriate for generation of arc-like geochemical signatures (Pearce, 2008). In particular, negative Nb-Ta anomalies may be related to ilmenite fractionation and assimilation of amphibolite (seawater-altered crustal components: Haase et al., 2005), resulting in arclike rocks without subduction. We show that the only places in the modern oceans where arc-like geochemical signatures exist in the absence of subduction are in the few felsic rocks from the Pacific Antarctic Rise (Haase et al., 2005; Freund et al., 2013) and dacites with negative Nb-Ta anomalies in Iceland (Willbold et al., 2009). The comparison with the Abitibi greenstone belt shows that very similar rocks, divided along clearly defined geochemical lines, existed in the Archean and perhaps for similar reasons.

Implications for the F classification scheme

Attempts to identify modern analogues of FI-FIV rhyolites have had mixed results (e.g., Prior et al., 1999; Piercey, 2011; Huston et al., 2019). It is well established that geochemical characteristics of Archean rhyolites are not directly comparable to those in Phanerozoic belts, which are dominated by FII and FIV types (Lentz, 1998; Piercey et al., 2001; Hart et al., 2004). Piercey (2011) characterized felsic volcanic rocks from juvenile, mafic-dominated post-Archean settings as having tholeiitic FIV-like signatures with low HFSEs and HREEs, typical of modern fore-arc or intra-arc rifting. Felsic volcanic rocks from more mature arc-back-arc systems, and especially those derived from partial melting of continental crust and sedimentary rocks, were found to have intermediate FII to FIII signatures with high HFSEs and HREEs and withinplate (A-type) to peralkaline affinity (Piercey, 2011).

Figure 8 shows the distribution of Abitibi greenstone belt rhyolites in FI-FIV space. Felsic volcanic rocks of the Bousquet Formation, with high Zr/Y and distinct negative Nb-Ta anomalies, have FI signatures similar to continental margin arc volcanic rocks of the Taupo Volcanic Zone but are less enriched than hot-spot volcanoes. Many of the Abitibi greenstone belt felsic volcanic rocks fall in the FII field, including samples from the Hunter mine area and Lac Herbert, although this type of rhyolite is mainly found in post-Archean felsic volcanic centers (Hart et al., 2004). Felsic volcanic rocks from Noranda, Joutel, Normetal, and Quevillon (App. Table A11; Fig. 10) overlap with the FII and FIIIa fields. The REE signatures of felsic volcanic rocks from Noranda are variable. They range from HFSE- and HREE-rich suites (almost FIIIb) to HFSE- and HREE-poor rocks (mostly FIIIa), similar to modern intraoceanic arc-back-arc systems. In contrast, FIIIb-type rhyolites from Matagami and Kidd Creek have significantly higher Yb(CN). The enriched sources were most likely hot mantle material that is thought to be directly linked to the formation of large VMS deposits in these areas. Prior et al. (1999) particularly compared FIIIb rhyolites from Kidd Creek to the rhyolites from plume-like sources in Iceland. Thus, Archean FIIIb rhyolites probably originated in extensional settings dominated by mid-ocean ridge-type spreading above an anomalous mantle, possibly in complex microplate mosaics like the northern Lau basin today. Several authors have pointed to the lack of post-Archean FIIIb signatures as an indication of the decrease in mantle temperatures and shift in the depth of melting through time (e.g., Lesher et al., 1986; Lentz, 1998; Hart et al., 2004). FIV-like felsic volcanic rocks are rare or absent in the Abitibi greenstone belt, probably reflecting the lack of volatile-fluxed melting.

Felsic volcanic rocks from modern settings span almost the entire range of the F classification (Fig. 9). Those associated with alkaline hot-spot magmatism and some calc-alkaline continental arcs overlap with the FI field. Felsic volcanic rocks from the Taupo Volcanic Zone, Myojin rift (Izu-Bonin), reararc volcanoes such as Niuatahi (Lau basin), and large arc volcanoes such as Montserrat (Lesser Antilles arc) overlap with FII. These rocks are all characterized by weak calc-alkaline affinity and steep REE patterns due to amphibole fractionation. Increasing crustal thickness, from Niuatahi to the Taupo Volcanic Zone, is reflected in increasing La/Yb ratios. Some samples from modern ridge-hot-spot intersections (Iceland, Azores, Ascension Island) also have FII signatures, but with higher Zr/Y, La/Yb and total HREE concentrations than Archean FII-rhyolites. Those from transitional or continental margin arc-back-arc systems, such as the Okinawa trough, Western Antarctic Peninsula, and the Eastern rifts of the Manus basin, overlap with FIIIa. Samples from intraoceanic arcback-arc systems, such as Izu-Bonin, the Mariana arc, and volcanoes of the Tonga-Kermadec arc, also plot in FIIIa. The arc influenced suite is distinguished from FIIIb by lower Yb(CN) and significant negative Nb-Ta anomalies compared to more distant back-arc spreading centers (e.g., Northwestern Lau, Central Lau, and Manus spreading centers: Fig. 3). Samples from distant back-arc spreading centers, together with felsic volcanic rocks from mid-ocean ridges and ridge-hot-spot intersections (e.g., Pacific-Antarctic Rise, East Pacific Rise, Galapagos spreading center, Juan de Fuca Ridge), overlap with FIIIb. They have flat REE patterns, slightly negative to absent Nb-Ta anomalies and significant HREE enrichment, indicating assimilation of hydrated oceanic crust (Haase et al., 2005; Wanless et al., 2010). Propagating rift tips and overlapping spreading centers, in particular, which are characterized by episodic magma injections and melt fractionation without continuous replenishment (Wanless et al., 2010), appear to be ideal settings for the generation of FIIIb-type rhyolites today. These complex spreading regimes may be the surface expression of small and hot mantle convection cells, leading to localized mantle upwelling and enhanced melting (Rychert et al., 2020). Felsic volcanic rocks from intraoceanic arc volcanoes and arc rifts, such as the Aogashima rift and Valu Fa rift, are FIV type. They are characterized by lower HREE concentrations, smaller negative Eu anomalies, and more pronounced negative Nb-Ta anomalies than the back-arc spreading centers (Fig. 3), consistent with the influence of the subducting slab.

Although there is general agreement with many of the inferred geodynamic settings of the F classification, it is clear that modern felsic volcanic rocks cross several different fields, and some fall completely outside the FI-IV classifications. Most Archean rhyolites show some degree of LREE enrichment due to the enriched mantle or greater crustal thickness (Lesher et al., 1986), whereas modern felsic volcanic suites have strongly depleted signatures. This, together with the numerous examples of alkali and peralkaline samples, accounts for many of the modern examples plotting outside the fields of Archean rhyolites (Hart et al., 2004; Piercey, 2011; Hollis et al., 2015). Although some have suggested that the F classification has outlived its usefulness with the widespread availability of immobile trace element data, it has yet to be replaced with a rhyolite classification scheme involving more elements, as is the case for mafic volcanic rocks (e.g., Agrawal et al., 2008).

The geochemical differences between Archean and post-Archean submarine felsic volcanic rocks have been interpreted mainly in terms of crustal and mantle conditions (Hart et al., 2004). An important observation of this study is that modern FIIIb-like dacites and rhyolites are found in a number of different tectonic settings, including at the tips of propagating rifts and overlapping spreading centers on mid-ocean ridges and at mature back-arc spreading centers such as the Northwestern Lau spreading center in the Lau basin and the Manus spreading center in the Manus basin. In the modern oceans, this type of crust has a low preservation potential owing to eventual subduction. Thus, the limited occurrences of FIIIb-type felsic rocks in post-Archean assemblages might simply reflect limited preservation. Studies of ophiolites indicate that the likelihood of preservation of these different types of oceanic crust is strongly dependent on the age of the crust, its thickness, and the thermal state and geometry of the plate boundaries (Shervais, 2001; Dilek and Furnes, 2014). Although the characterization of Archean greenstone belts as obducted oceanic crust is debated (Furnes et al., 2015; Palin et al., 2020), the abundance of FIIIb-type felsic volcanic rocks in the Archean may indicate more complicated microplate architectures and increased probability of preservation in buoyant microplate domains (e.g., Wyman, 2013).

We compared the major and trace element geochemistry of modern submarine felsic volcanic rocks using PCA, hierarchical clustering, and random forest to identify important similarities and differences between tectonic settings. Automated hierarchical clustering reveals significant geochemical diversity in felsic volcanic rocks from subduction zone settings related to fluids from the subducting slab but also considerable variation related to local geodynamic influences. Mid-ocean ridge-like dacites found in the complex microplate mosaics of the northern Lau basin and the Manus basin are related to anomalous heat flow, rising asthenospheric mantle, and rift propagation especially at overlapping spreading centers. They have flat REE patterns, significant HREE enrichment, and strong negative Eu anomalies, similar to Archean FIIIb rhyolites, such as felsic volcanic rocks from Matagami. Major differences in the abundance of Archean and post-Archean felsic volcanic rocks of different types may be at least partly explained by their preservation potential. Modern mid-ocean ridges and mature back-arc spreading centers have low preservation potential, owing to the eventual subduction at basin margins. In contrast, the complex microplate mosaics of Archean greenstone belts might have been important settings for preservation of this type of crust, promoted by the thickness of the crust and high heat flow in the Archean. While they are not exact tectonic analogues, the geochemical evolution of felsic volcanic rocks in Archean greenstone belts may share many of the same influences observed in modern microplate breakouts, such as in the northern Lau basin.

Felsic volcanic rocks from transitional arc-back-arc systems, such as the eastern Manus basin, are characterized by calcalkaline to tholeiitic affinities with negative Nb-Ta anomalies, positive La/Yb(CN) ratios, and variable HREE concentrations, similar to rhyolites of the Noranda Volcanic Complex. Felsic volcanic rocks at continental margin arcs, such as the Taupo Volcanic Zone and Okinawa trough, are characterized by calcalkaline affinities with negative Nb-Ta anomalies, positive La/Yb(CN) ratios and high Th concentrations, similar to the rhyolites of the Bousquet Formation. Dacites and rhyolites from high heat-flow environments, mid-ocean ridges, and distant spreading centers, characterized by tholeiitic affinity and flat REE patterns, are similar to rhyolites from Matagami. Rhyolites from ridge-hot-spot intersections, which have tholeiitic to weak alkaline affinity and flat REE patterns to slightly positive La/Yb(CN) ratios, are geochemically close relatives of rhyolites at Kidd Creek.

We thank Pierre-Simon Ross and Stephen Piercey for the helpful comments on the improvement of this manuscript and the comparisons with Abitibi greenstone belt felsic volcanic rocks. This research was funded through the Canada First Research Excellence Fund (CFREF, Metal Earth) and the Marine Mineral Resources Group at the GEOMAR Helmholtz Centre for Ocean Research Kiel. The Helmholtz Association, the Natural Sciences and Engineering Research Council of Canada (NSERC), 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 is also supported by the NSERC Collaborative Research and Training Experience program (iMAGE-CREATE) on Marine Geodynamics and Georesources. This is contribution MERC-ME-2022-12.

Marc Lorin Fassbender is a Ph.D. student at the University of Ottawa, Canada. His work is focused on host rocks of actively forming sea-floor massive sulfide deposits and their ancient analogues, volcanogenic massive sulfide deposits, using geochemistry, petrology, and machine learning. He completed his M.Sc. and B.Sc. degrees at the University of Erlangen-Nuremberg, Germany, working on the petrogenetic evolution and exceptional HREE enrichment in fayalites of the Vergenoeg F-Fe-REE iron oxide copper-gold deposit.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data