A Mesoarchean Large Igneous Province on the Eastern Kaapvaal Craton (Southern Africa) Confirmed by Metavolcanic Rocks from Kubuta, Eswatini

The Archean eon covers much of Earth’s history and was a period of major crust forming events (e.g., [1–3]) that culminated in the stabilization of the ancient cratons. One of these major crust-forming events occurred in the Mesoarchean (ca. 3.0 to 2.9 Ga) on the eastern Kaapvaal Craton in southern Africa producing large volumes of continental flood basalts in an intracratonic rift basin (e.g., [4–7]). These Mesoarchean basalts have well-preserved magmatic textures and preserve their primary geochemical composition, which can be used to place constraints on the genetic relationships to other units of the Kaapvaal Craton.

The correlation of lavas of the eastern Pongola basin and western Witwatersrand basin in terms of age and trace element composition led Gumsley et al. [6] to the conclusion of the existence of a so far unknown Mesoarchean large igneous province, which might extend further West to the Dominion Group [8]. Mainly based on precise age determinations and field observations, these authors proposed that the large igneous province was formed by two magma pulses derived from a short-lived mantle plume. In their petrogenetic model, the first magma pulse formed the basaltic lavas of the ca. 3.0 Ga Nsuze Group (Pongola Supergroup) and the gabbros and cumulate rocks of the Usushwana Igneous Complex (Figure 1; Gumsley et al. [4]). After some time of subsidence and sedimentation, the second pulse of magma formed the ca. 2.9 Ga lavas of the Mozaan Group (Pongola Supergroup; Figure 1) and the mafic sills of the Hlagothi Complex [6].

However, there is some disagreement on the nature of the parental magmas that formed the Nsuze volcanic rocks. While Hegner et al. [5] proposed an ultramafic komatiitic parental magma for the Nsuze volcanic rocks, Armstrong et al. [9] suggested a primitive basaltic parental magma. Likewise, the proposed crustal contamination processes that affected these parental magmas vary from bulk assimilation [5] and assimilation-fractional crystallization [9] to magma mixing of basaltic and felsic melts [10]. However, all aforementioned studies agree that a contamination by felsic crustal rocks from the Ancient Gneiss Complex (AGC) of Eswatini occurred. Yet, these studies do not specify, if the crustal contaminant to the Nsuze volcanic rocks belongs to the 3.2 Ga or 3.5 Ga crust-forming event within the AGC.

In this study, we sampled amphibolites from the Kubuta area of southwestern Eswatini. The affiliation of the Kubuta samples is unknown, but from a combination of major and trace element, as well as Lu-Hf and Sm-Nd isotope data, we will show that they could either be related to volcanic rocks of the 2.98-2.84 Ga Pongola Supergroup (Figure 1) [10–14] or the ca. 2.99 to 2.87 Ga Usushwana Igneous Complex (Figure 1) [4, 5] and the ca. 2.86 Ga Hlagothi Complex [6]. Alternatively, the Kubuta volcanic rocks could be related to the ca. 2.66 Ga White Mfolozi Dike Swarm (WMDS) [15] further south in the KwaZulu-Natal province.

Using our extensive data set, we aim to better constrain the composition of the mantle sources and the contamination history of the aforementioned sequences to test the reliability of the petrogenetic model of Gumsley et al. [4] that was mostly based on geochronology. In combination with literature data, we propose a refined petrogenetic model for the Nsuze, Mozaan, Usushwana, Hlagothi, and Kubuta volcanic rocks that also takes geochemical evidence into account.

The eastern Kaapvaal Craton, located in Eswatini and eastern South Africa, records a continuous geological history from the early Archean to the Neoarchean recording the stabilization of continental lithosphere in this region. The first stabilized continental crust is represented by the 3.66-3.20 Ga Ancient Gneiss Complex (AGC) [19–22] that predominantly contains complexly deformed Paleo- to Mesoarchean grey gneisses that are rarely interlayered with amphibolites (e.g., [23–25]). The oldest component of the AGC is the ca. 3.66-3.45 Ga Ngwane Gneiss (Figure 1) that underwent ductile metamorphism under amphibolite to granulite conditions [20, 26]. The Ngwane Gneiss is intruded by the genetically linked 3.48-3.43 Ga Tsawela Gneiss [26–28] as well as by several younger gneisses and granitoids like the 3.23-3.22 Ga Usutu Intrusive Suite in northern and central Eswatini [29] and the 3.28-3.24 Ga Nhlangano and Mahamba Gneisses in south-central Eswatini (Figure 1) [30]. The ca. 3.1 Ga Pigg’s Peak and Mpuluzi Batholiths [31, 32] represent the last major intrusions exposed in the eastern Kaapvaal Craton [22, 30, 33]. These intrusions separate the AGC in the Southeast from the Barberton Greenstone Belt in the Northwest. In its most southerly occurrence, the Mpuluzi Batholith forms the basement to the ca. 2.98-2.84 Ga Pongola Supergroup (Figure 1) [10–14].

The Pongola Supergroup was deposited in an area of about 24000 km² extending from south-central Eswatini into northern and central KwaZulu-Natal in eastern South Africa and belongs to one of the oldest cratonic cover sequences in the world [34]. Recently, new age determinations lead to the finding that the Pongola Supergroup also extends further West into the 2.96 Ga Dominion Group revealing that a craton-wide igneous province was active between ca. 2.99 and 2.96 Ga [8]. The succession was deposited in a volcanic rift setting [35] between contrasting crustal blocks of the AGC with generally older granitoids ranging in age from ca. 3.6 to 3.2 Ga to the Northeast [11] and ca. 3.3 to 3.2 Ga to the Southwest [36] of the Pongola basin. In the northern and central areas of the basin, the metamorphic grade of the volcanic and sedimentary rocks of the Pongola Supergroup ranges from largely greenschist up to granulite facies in some areas [11].

The Pongola Supergroup is divided into the lower, mostly volcanic, ca. 2.98 Ga Nsuze Group and the upper, predominantly sedimentary, ca. 2.86 Ga Mozaan Group (Figure 1) [37, 38]. The volcanic rocks of the Nsuze Group range from metabasalts to metarhyolites [38] and were deposited in an intracontinental environment [39] under subaerial to locally subaqueous conditions [9, 40]. The Mozaan Group rests unconformably on the Nsuze Group and mostly contains sandstone, shale, mudstone, and locally banded iron formations and stromatolites indicating deposition in a shallow marine environment [37, 41–44]. Towards the top of the succession, volcanic rocks become more dominant again [44].

The Nsuze and Mozaan Groups were intruded by the Usushwana Igneous Complex (UIC; Figure 1), the Thole Complex, and the Hlagothi Complex. The UIC has an emplacement age of $2871±30$ Ma [5]. However, recent U-Pb baddeleyite age dating by Gumsley et al. [4] on gabbroic samples from the Piet Retief Suite (now re-named to Mkhondo Suite) that is part of the UIC yielded significantly older ages up to 2989 Ma falling in the age range of Nsuze volcanic rocks. Based on their new ages, Gumsley et al. [4] interpreted the Usushwana Igneous Complex as being a feeder dike system for Nsuze lavas. The $2866±2$ Ma Hlagothi Complex and the related $2874±2$ Ma Thole Complex consist of several gabbroic sills that are found in the KwaZulu-Natal province and are interpreted to be the feeder dikes for the Mozaan flood basalts [6]. Based on the similarity of ages and composition, Gumsley et al. [6] proposed that the UIC and the Thole and Hlagothi Complexes, together with the Mozaan volcanic rocks form part of a Mesoarchean large igneous province located in the southeastern part of the Kaapvaal Craton and possibly extending into the eastern Witwatersrand block.

Later intrusions of ca. 2.85-2.70 Ga granitoids [30] were triggered by extensional deformation that also caused local deformation of the Pongola Supergroup [45]. Following this event, ca. 2.70 to 2.66 Ga extensional dolerite dike swarms, sometimes referred to White Mfolozi Dike Swarm, crosscut the eastern Kaapvaal Craton terminating its Archean crust formation history (e.g., [15, 46–48]).

Seven samples have been collected from the area of the Kubuta Ranch in south-central Eswatini (Figure 2). These samples are referred to as ‘Kubuta volcanic rocks’ in the following. Samples KS-SW 31 to 35 have been collected east of the Sibowe shear zone at the MR 25 (Figure 2) in a geologic succession that was mapped by the Geological Survey of Eswatini as “amphibolite of uncertain correlation” [49]. In contrast, samples KS-SW 45 and 47 have been collected west of the Sibowe shear zone as boulders in the river bed of the Matsanjeni River (Figure 2). Hence, one of the aims of this study is to constrain the relationship of these rocks in terms of mantle source and contamination history to other volcanic sequences in this region, i.e., the 2.98-2.84 Ga Pongola Supergroup (Figure 1) [10–14], the ca. 2.99 to 2.87 Ga Usushwana Igneous Complex (Figure 1) [4, 5], the ca. 2.86 Ga Hlagothi Complex [6], or the ca. 2.66 Ga White Mfolozi Dike Swarm (WMDS) [15] further south in the KwaZulu-Natal province.

All samples appear relatively fresh with preserved primary magmatic mineral textures. However, when taking a closer look at the samples in thin section, the rocks can be identified to have undergone greenschist-facies metamorphism expressed by variably chloritized mineral assemblages and formation of epidote. Sample KS-SW 31 has a microgabbroic texture indicating a moderate cooling rate and containing minor relict olivine and abundant relict clinopyroxene and plagioclase (Figures 3(a) and 3(b)). Samples KS-SW 33 to 35 exhibit an intergranular texture with usually bigger needle-shaped amphibole crystals and a medium-grained matrix of clinopyroxene, orthopyroxene, plagioclase, and epidote. Sample KS-SW-32 is a cumulate containing orthopyroxene of up to 5 mm in a diameter. Samples KS-SW 33 and 35 have radiating and branching fibrous plagioclase crystals preserved (Figures 3(c) and 3(d)). Carbonatization along mineral cracks is rare in most samples, except in KS-SW 32, which shows relatively strong carbonatization and chloritization. However, the original magmatic texture of this sample is still preserved.

In contrast, samples KS-SW 45 and 47 show a porphyric texture having either 1-2 mm amphibole prisms and clinopyroxene crystals floating in a fine-grained plagioclase-bearing matrix (KS-SW 45; Figure 3(e)) or plagioclase and clinopyroxene crystals floating in a fine-grained orthopyroxene-bearing matrix (KS-SW 47). These two samples show almost no carbonatization. Furthermore, apatite can be found in plagioclase phenocrysts and different opaque phases occur as accessory components within the matrix.

The rock samples were cut using a diamond blade saw to remove alteration rims, and thin sections were prepared. The remaining sample material was crushed in a steel jaw crusher and ground in an agate mill. Whole rock major and some trace element compositions including the contents of H₂O and CO₂ were analyzed by X-ray fluorescence (XRF) on Li₂B₄O₇-flux fusion discs using a PANalytical Axios X-ray spectrometer at the ElMiE Lab at the German Research Centre for Geosciences (GFZ) in Potsdam, Germany (Table 1). Loss on ignition (LOI) was determined by weight difference after fusion. Reproducibility was determined on three certified reference materials (crm) and is within the analytical precision, which is better than 2% for main elements and better than 10% for trace elements.

The trace element compositions of all samples were analyzed using a Thermo Finnigan Element XR SF-ICP-MS at Freie Universität Berlin following the protocols of Schneider et al. [50]. Sample powders were digested in Parr® pressure vessels in concentrated HF-HNO₃ for 72 hours and were subsequently dried down. The evaporation step was repeated with 2 ml concentrated HNO₃. After evaporation to dryness, samples were diluted in 3% HNO₃ for trace element measurements. The analytical procedure for trace element analysis was carried out following Schneider et al. [50]. Relative total uncertainties (RSD) for all measured elements usually range from 4% to 18%. The accuracy of our measurements can be inferred from multiple analyses of the BHVO-2 reference material and was better than 5% for most elements and between 6% and 12% for Ga, Sr, Y, and Nb (Table 2). Additionally, the total procedural uncertainty can be inferred from sample replicates that were analyzed in four analytical sessions and were better than 5% for the REEs, Ba, Th, U, and Co, better than 11% for Ga, Y, Ta, Sc, V, and Cu, and 25% for Mo and 40% for Pb (Table 2).

Hafnium and Nd isotope compositions and Lu, Hf, Sm, and Nd concentrations by isotope dilution were analyzed (Table 3) using a Thermo Scientific Neptune multicollector ICP-MS at the University of Cologne (Germany), operated in static mode. The analyses were obtained by isotope dilution using mixed ¹⁷⁶Lu-¹⁸⁰Hf and ¹⁴⁹Sm-¹⁵⁰Nd tracers following the protocols of Münker et al. [51] and Weyer et al. [52]. Chemical separation of Lu, Hf, Sm, and Nd followed the procedures outlined in Pin and Zalduegui [53] and Münker et al. [51].

Values of ¹⁷⁶Hf/¹⁷⁷Hf and ¹⁴³Nd/¹⁴⁴Nd were corrected for mass bias using the exponential law and assuming $146Nd/144Nd=0.7219$ and $179Hf/177Hf=0.7325$⁠, respectively. Repeated analyses of the Cologne AMES Hf solution standard yielded an external reproducibility of ±41 ppm (2 s.d., $n=4$⁠) for ¹⁷⁶Hf/¹⁷⁷Hf, whereas repeated analysis of the AMES Nd and JNdi-1 standards yielded external reproducibilities of 21 ppm (2 s.d., $n=4$⁠) and 70 ppm (2 s.d., $n=10$⁠) for ¹⁴³Nd/¹⁴⁴Nd, respectively, over the course of the analytical session. Reported values are given relative to $143Nd/144Nd=0.512115$ for JNdi (i.e., $143Nd/144Nd=0.511858$ for La Jolla Nd [54]) and $176Hf/177Hf=0.282160$ for the Münster AMES Hf standard, which is isotopically indistinguishable from the JMC-475 standard for Hf [55]. All procedural blanks handled along with the sample batches during the course of this study were better than <3 pg for Lu, <37 pg for Hf, <22 pg for Sm, and<132 pg for Nd and therefore were insignificant. CHUR parameters used to calculate epsilon values are from Bouvier et al. [56] and an age of 2980 Ma was estimated (see discussion in Section 6.2).

The SiO₂ and MgO contents in the samples range from 44.4 to 53.2 wt.% and from 7.17 to 22.9 wt.%, respectively (Table 1). Magnesium numbers (Mg#) range from 54 to 78 corresponding to the variability in MgO and FeO_total concentrations. Compared to literature data, our samples have much higher MgO contents than the Nsuze Group volcanic rocks (Pongola Supergroup; Figures 4(a)–4(d)) [9–11, 46] and volcanic rocks from the WMDS [15]. However, our Kubuta samples agree in MgO content with the volcanic rocks of the Mozaan Group (Pongola Supergroup) [10, 57], the Hlagothi Complex [6], and the Usushwana Igneous Complex [10]. The two samples of the latter unit from the study of Nhleko [10] were originally grouped as Mozaan volcanic rocks. However, Nhleko [10] supposed in his thesis that they might rather belong to the Usushwana Igneous Complex based on their distinct geochemical composition. Adopting this interpretation, we plot these two samples as UIC in the diagrams of this study. Unfortunately, there is currently no more literature data from the UIC to test a possible overlap in element abundances of our Kubuta samples with the UIC; only chondrite-normalized REE data are published for this unit [5].

In Fe₂O₃ vs. MgO space, a positive correlation of samples of the Nsuze and Mozaan suites below 5 wt.% MgO (Figure 4(a)) suggest fractionation of pyroxene and titanomagnetite. In addition, fractionation of plagioclase and pyroxene can be identified at MgO contents below 5 wt.% by the positive correlations of Al₂O₃ and CaO vs. MgO (Figures 4(b)–4(d)), respectively. Furthermore, the positive correlation of Eu anomalies with CaO contents indicates plagioclase fractionation. Fractionation of apatite is indicated by the negative correlation of P₂O₅ above a SiO₂ content of 60 to 65 wt.% (Figure 4(e)), and fractionation of zircon is indicated by the negative correlation of Zr and SiO₂ above a SiO₂ concentration of around 70 wt.% (Figure 4(f)). However, another trend is visible in Figure 4(f) which cannot be explained by zircon fractionation, where Zr increases with increasing SiO₂, possibly representing the involvement of Zr-rich partial melts of the continental crust. Compared to the fractionation trends observed from literature data, our samples mostly plot at the high-Mg side of these trends, resembling the high-Mg samples from the Usushwana Igneous Complex and Hlagothi Complex.

The trace element concentrations of the samples are summarized in Table 2. The primitive mantle-normalized trace element patterns for Kubuta samples show an enrichment of light rare earth elements (LREE) and flat to slightly depleted heavy REEs (Figures 5(a) and 5(b)). Consequently, most samples exhibit (La/Yb)_CN (i.e., chondrite-normalized) values greater than one, ranging from 1.90 to 6.15, except for sample KS-SW 47 that has (La/Yb)_CN of around one. Four samples have (Gd/Yb)_CN values around one and in three samples the (Gd/Yb)_CN values range from 1.28 to 2.21. Pronounced negative anomalies are observed for Nb $Nb/Nb∗PMN=0.44−0.92$⁠, Ta $Ta/Ta∗PMN=0.44−0.80$⁠, and Ti. Values of Th/U range from 3.3 to 6.0 in agreement with Archean upper crustal values, i.e., $Th/U>3$ (e.g., [58]). Furthermore, negative Eu anomalies are present for four samples ranging from $Eu/EuCN∗=0.54to0.95$⁠, indicating fractionation of plagioclase, whereas sample KS-SW 31 shows a small positive $Eu/EuCN∗$ of 1.15. All other samples lack Eu anomalies. Furthermore, all samples show elevated concentrations of fluid mobile elements such as Rb, Ba, and Pb.

When compared to the literature data for the Nsuze Group [9–11, 46] and Mozaan Group [10, 57] volcanic rocks, our samples are generally less enriched in incompatible trace elements, yet exhibit similar (La/Yb)_PMN and (Gd/Yb)_PMN values and elemental anomalies (Figure 5(a)). When compared to samples from the Usushwana Igneous Complex [10] and Hlagothi Complex [6], our samples match well their elemental abundances and exhibit similar (La/Yb)_PMN and (Gd/Yb)_PMN values and elemental anomalies (Figures 5(a) and 5(b)). However, compared to volcanic rocks from the White Mfolozi Dike Swarm [15], the Kubuta samples have higher (La/Yb)_PMN and (Gd/Yb)_PMN values and pronounced negative Nb and Ta anomalies that are not observed in the volcanic rocks of the WMDS.

The Lu-Hf and Sm-Nd isotopic compositions and Lu, Hf, Sm, and Nd concentrations of our sample suite are summarized in Table 3. For calculation of initial εHf and εNd values, an age of 2980 Ma was applied ([11]; see discussion in Section 6.2). Most of the Kubuta samples plot on the terrestrial array in εHf_(t) vs. εNd_(t) space after Vervoort et al. [60] that we use as a measure for a possible decoupling of initial Hf and Nd isotope compositions. However, sample KS-SW 31 has much higher εHf_(t) (+26.6) at its εNd_(t) value (+6.3) and plots off the array (see in Table 3). Thus, this sample was not included in the plot in Figure 6, considering it as disturbed. All other samples have εHf₍₂₉₈₀₎ and εNd₍₂₉₈₀₎ values ranging from -4.8 to +3.5 and from -2.5 to +1.1, respectively.

Literature data for rocks of the Nsuze Group and Usushwana Igneous Complex yielded εNd_(t) values ranging from -1.8 to -7.8 [5, 11] and -3.7 to -1.6, respectively. These values are in good agreement with three of our samples yielding negative εNd_(t) (and negative εHf_(t)) values. However, the other three Kubuta samples have chondritic to slightly elevated εNd_(t) (and εHf_(t)) values at comparable neodymium concentrations (Figure 6). Unfortunately, whole-rock Lu-Hf isotope compositions were not available in the literature for further comparison. Likewise, no Hf-Nd isotope data were available for the Mozaan volcanic rocks, the Hlagothi Complex, and the White Mfolozi Dike Swarm to compare with our Kubuta samples.

Great care must be taken when interpreting the geochemistry of Archean volcanic rocks, because these rocks underwent metamorphism that likely caused mobility of some elements which lead to alteration of the element inventory of these rocks. In Archean volcanic rocks, elements typically known as relatively mobile during metamorphism include Rb, K, Na, Sr, Ba, Ca, Fe, P, and Pb [62, 63]. Thus, these elements should be used with caution to constrain magmatic processes. In contrast, the effects of metamorphic alteration on high field strength elements (HFSE), REEs, Al, Cr, and Ni are usually minor, and therefore, these elements are mostly considered immobile on the whole-rock scale [63–66].

To distinguish the elemental disturbance of Archean greenstones, Polat et al. [67] developed five requirements that need to be fulfilled in order to consider the greenstones to be largely undisturbed in their trace element composition. Most of our investigated greenschist-facies metabasalts meet all five criteria of Polat et al. [67] in that they (1) generally have volatile (i.e., H₂O plus CO₂) contents lower than 6 wt.% (Table 1); (2) show coherent trace elements patterns of immobile elements normalized to PM (Figure 5); (3) display correlations of REEs, Th, Ti, and Nb abundances with Zr content; (4) show no correlation of their Th-Nb-La interelement ratios with volatile (i.e., H₂O plus CO₂) abundance; and (5) show no correlation of element abundance with the degree of chloritization or carbonatization. However, sample KS-SW 32 has a high volatile content (i.e., H₂O plus CO₂) of 7.9 wt.%, resulting from its relatively strong chloritization and carbonatization. This sample also has much higher MgO (22.9 wt.%), lower Na₂O (0.01 wt.%) and K₂O (0.03 wt.%) contents than less altered samples, likely indicating mobility of Mg, Na, and K during greenschist-facies metamorphism. Consequently, we consider the major element composition of sample KS-SW 32 to be disturbed. In contrast, chloritization and carbonatization appear not to have disturbed the REE, HFSE, and the Hf-Nd isotope systematics of this sample. All other samples do not seem to have disturbed major and trace elemental compositions.

As stated above, sample KS-SW 31 plots off the modern terrestrial array [60] in εHf_(t) vs. εNd_(t) space towards much higher εHf_(t) at a given εNd_(t) compared to the other Kubuta samples and compared to literature Nd isotope data for the Pongola and Usushwana suites [5, 11]. These high εHf_{(2980 Ma)} at εNd_{(2980 Ma)} could be a result of incorrect age assignment: This sample might belong to a different, much younger, volcanic suite, possibly to the Karoo flood basalts, leading to overestimated radiogenic εHf_(t) at εNd_(t) values. Other arguments in favor of this hypothesis are, for instance, the different petrographic textures (microgabbro) leading to the positive Eu anomaly of this rock compared to the other samples. Consequently, this sample is excluded from our further discussion.

Because of possible disturbance of the major element composition of our samples, we used the discrimination diagram of Winchester and Floyd [68] that was updated by Pearce [69] for all magma types applying the immobile trace element ratios Zr/Ti and Nb/Y, instead of the total alkalis vs. silica (TAS) diagram after Le Maitre et al. [70]. In the Zr/Ti vs. Nb/Y diagram, the former ratio reflects the degree of differentiation, and therefore is comparable to SiO₂, and the latter tracks enrichment, and, hence, is comparable to the alkali content. With this approach, the Kubuta lavas can be classified as subalkaline (⁠$Nb/Y<0.9$⁠) basalts and minor basaltic andesites overlapping with the Hlagothi and Usushwana Igneous Complexes, the White Mfolozi Dyke Swarm, and most of the Nszue and Mozaan lavas (Figure 7).

Based on our isotope analysis, we obtained a Sm-Nd whole-rock regression line corresponding to an errorchron of $2780±580$ Ma (⁠$MSWD=65$⁠; Figure 8(a)) and a Lu-Hf whole-rock regression line that corresponds to an errorchron of $2973±620$ Ma (⁠$MSWD=113$⁠; Figure 8(b)) for five samples of the Kubuta sample suite. However, the uncertainties on the calculated ages are quite large which may result from the relatively small spread in ¹⁴⁷Sm/¹⁴⁴Nd and ¹⁷⁶Lu/¹⁷⁷Hf ratios in our samples. Alternatively, the large uncertainties on the regression lines may indicate that our samples might not be cogenetic or be contaminated by crustal material. However, the data still give a first-order approximation on the geological unit to which the Kubuta metabasalts could be assigned to. Hence, the most likely units that could correlate with the Kubuta samples in terms of age and composition are the mafic/ultramafic Nsuze and Mozaan Groups of the Pongola Supergroup, as well as the mafic/ultramafic Usushwana Igneous Complex (UIC), the mafic/ultramafic Hlagothi Complex, and the mafic volcanic rocks of the White Mfolozi Dike Swarm.

The age of the Nsuze Group is well constrained by zircon U-Pb dating of evolved volcanic rocks which range from $2977±5$ Ma to $2985±1$ Ma in age [10–12]. In contrast, a Sm-Nd whole-rock errorchron ‘age’ obtained from basaltic to rhyolitic samples yielded a slightly younger, but within uncertainty overlapping, age of $2934±114$ Ma [5]. The age of the overlying Mozaan Group ranges from $2860±26$ Ma [13] to $2869±5$ Ma [4], respectively. The mostly intrusive Usushwana Igneous Complex (UIC) yielded Rb-Sr whole-rock age data of $2871±30$ Ma [5] and $2875±40$ Ma [71], postdating the Nsuze volcanic rocks. However, recent U-Pb baddelyite dating by Gumsley et al. [4] yielded slightly older ages ranging from $2978±2$ Ma to $2990±2$ Ma that indicate a synchronous evolution of the UIC and the Nsuze Group. Additionally, the ages for the Mozaan Group and UIC overlap within error with the ca. 2866 Ma mafic to ultramafic sills of the Hlagothi Complex [6] that intrude the basal units of the Nsuze Group. The WMDS was dated to $2662±2$ Ma [15], also overlapping with our Sm-Nd and Lu-Hf whole-rock errorchron dates.

In summary, based on our Sm-Nd and Lu-Hf errorchron ages, the Kubuta volcanic rocks could be assigned to all the aforementioned volcanic units. However, in terms of trace elemental patterns combined with the Nd isotope composition, the Kubuta volcanic rocks share more similarities with the lavas of the Nsuze and Mozaan Groups [9–11, 57] as well as with the Usushwana [c.f., 5] and Hlagothi Igneous Complexes [c.f., 6] rather than with the lavas of the White Mfolozi Dike Swarm [c.f., 15] (Figures 5(a), 5(b), and 6(a)).

Hence, we exclude the possibility that the Kubuta volcanic rocks are related to the White Mfolozi Dike Swarm. Consequently, based on the good match of major, trace element, and Nd isotope abundances, we assign the Kubuta samples to be part of the proposed large igneous province of Gumsley et al. [6], consisting of mafic and ultramafic rocks from the Mozaan Group, the Hlagothi Complex, and the Usushwana Igneous Complex. However, the role of these units in the proposed large igneous province and their relation to the Nsuze lavas needs further investigation.

Crustal contamination of Nsuze volcanic rocks has been suggested by several authors based on the elemental composition and negative εNd_(t) values of the lavas [5, 9, 10, 46]. As potential contaminants to the Nsuze volcanic rocks, granitoids from the Ancient Gneiss Complex (AGC) were favored [5, 9]. The proposed assimilation processes vary from bulk assimilation of ca. 3.5 Ga old granitoid crust [5] to assimilation-fractional crystallization (AFC) of more granitic AGC crust [9]. In addition, Nhleko [10] proposed magma mixing of basaltic melts and felsic melts derived from partial melting of the continental crust to explain the diverse compositions of Nsuze volcanic rocks.

The arguments in favor of an involvement of continental crust in the petrogenesis of the volcanic rocks from the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta sample suites are mostly their enrichment of incompatible elements, especially La and Th, in conjunction with negative anomalies for Nb and Ta (Figures 4 and 5). Additionally, consistent primitive mantle normalized trace element patterns and mostly negative εNd_(t) values for the aforementioned volcanic rocks suggest a similar contamination history for all suits. In Th/Yb vs. Nb/Yb space, all samples plot well above the modern MORB-OIB array towards higher Th/Yb and Nb/Yb indicating substantial contamination of the parental magmas resulting from Archean crustal input and/or partial melting and fractional crystallization (e.g., [72]) (Figure 9). However, one type of crustal contaminant alone cannot account for the diversity in compositions of the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta samples (e.g., Figure 9).

Possible crustal contaminants that could have modified the composition of the aforementioned ca. 2.9 to 2.8 Ga volcanic rocks are granitoid gneisses from the AGC either belonging to the ca. 3.5 Ga (e.g., [5, 9, 10]) or the ca. 3.2 Ga crust-forming events (for a review on the ca. 3.5 and 3.2 Ga crust-forming events in the eastern Kaapvaal Craton, please refer to [61]). Contamination most likely occurred at the upper crustal levels (e.g., [5, 9, 10]). Arguments in favor of AFC are fractionation of pyroxene and plagioclase as indicated by the fractionation trends in Figure 4 and by the depletions in Ti and Eu illustrated in Figure 5, as well as by the typical crustal contamination trend displayed in εNd_(t) vs. La/Sm space (Figure 6(b)).

In contrast, the increase in Zr concentration with increasing SiO₂ for Nsuze rhyolites, illustrated in Figure 4(f), is likely a result of the increased alkalinity of these rocks which hampers zircon fractionation (e.g., [73]), suggesting that these rocks formed from crustal melts (e.g., [10]). Hence, magma mixing of mafic melts with felsic partial melts of older crust could explain the evolved compositions of the volcanic rocks of the Nsuze and Mozaan Groups, too, as previously suggested by Nhleko [10].

In the following, we will discuss the mechanisms of the two different contamination processes, i.e., magma mixing and assimilation-fractional crystallization.

(1) The Role of Magma Mixing. Partial melting of AGC crust and subsequent magma mixing with mafic underplated magma was suggested by Nhleko [10] to have formed the Nsuze volcanic rocks. One argument for partial melting of AGC crust is based on one rhyolite from the study of Hegner et al. [5] which yields a strongly negative εNd_(t) value of -7.8. Partial melting of crustal material would shift the resulting magma composition towards higher La/Yb and Zr concentrations and towards slightly lower Th/Yb and Nb/Yb than the AGC crustal rocks. This is observed for rhyolites of the Nsuze and two andesites of the Mozaan and Hlagothi suites, respectively (Figure 9). Hence, it is possible that at least some andesites and basaltic andesites of the Nsuze and Mozaan Groups might have been formed by magma mixing with subsequent fractional crystallization.

To explore the role of magma mixing, we calculated mixing paths using the Microsoft Excel calculation spreadsheet of Ersoy and Helvacı [84]. As a mafic mixing end member, we chose basalt AG-F of the Hlagothi Complex [6], because it has the lowest Nb/Yb at low Th/Yb, low La/Yb, and medium MgO (8.36 wt.%) and therefore serves as the best approximation of the composition of the parental melt. As felsic mixing end member, we used rhyolite NZ-39 of the study of Nhleko [10], because it has high SiO₂ (72.2 wt.%), very low MgO (0.05 wt.%), and the highest Zr concentration (601 ppm), and consequently most likely represents a partial melt of the crust. The results of magma mixing are illustrated in Figure 10. With this approach, most of the Nsuze and Mozaan volcanic rocks having $Th/Yb>1.8$ can be reproduced by 10 to 50% mixing of a rhyolitic melt with a basaltic melt (Figures 10(a) and 10(b)). However, many basalts, basaltic andesites, and andesites of the Nsuze and Mozaan Groups plot towards lower Th/Yb at similar Nb/Yb (Figure 10(a)). These lavas cannot be reproduced by magma mixing. Furthermore, magma mixing would have caused more heterogeneity in Nd isotope composition as a function of magma composition, i.e., the more felsic the magma, the lower the εNd_(t) should be. However, this is not reflected by the relatively uniform negative εNd_(t) values in most Nsuze, Usushwana, and three Kubuta volcanic rocks ranging between -2 and -4 (Figure 6). Consequently, a process other than simple magma mixing must have played a role in the formation of the Nsuze and Mozaan, but also the Hlagothi and Kubuta volcanic rocks, i.e., assimilation-fractional crystallization.

(2) The Role of Assimilation-Fractional Crystallization. To constrain the possible crustal contaminants to the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta volcanic rocks, we calculated AFC paths using the Microsoft Excel calculation spreadsheet of Ersoy and Helvacı [84]. For modeling the contamination history of Nsuze and Mozaan volcanic rocks, a komatiite from the study of Schneider et al. [50] was used as analogue for the parental komatiitic magma, having 26.9 wt.% MgO (sample ZA 33). Our approach follows the conclusion of Hegner et al. [5] who proposed a komatiitic parental magma for Nsuze volcanic rocks. For all other volcanic suites including the Kubuta samples of this study, we tried different mafic endmembers and chose a modeled 30% batch melt of incompatible element-depleted mantle (after Salters and Stracke [74]; equations after Shaw [82]; partition coefficients for mantle melting are listed in Table S 1), because it gave the best fit for our modelling. As possible contaminants, different crustal rocks from the AGC were tested, ranging in age from ca. 3.5 to 3.2 Ga (Figure 11).

The crustal compositions which yield AFC paths that best describe the evolved compositions of the Nsuze and Mozaan volcanic rocks are two granitoid gneisses of the 3.2 Ga crust-forming event, of which one has a tonalitic composition and was sampled in the vicinity of the Kubuta samples of this study (AGC 473, 3241 Ma) [61]. The other one (AGC 449, 3229 Ma) [61] has a granodioritic to granitic composition and was sampled from the Piggs Peak inlier, the northernmost outcrop of AGC crust, close to the Phophonyane shear zone. Although this latter sample is not from the vicinity of exposed Nsuze or Mozaan volcanic rocks, we take its composition as an approximation of the composition of the crustal contaminant, because the two granitoid gneisses belong to the same crust-forming event at ca. 3.2 Ga, which is considered to have been widespread in the AGC (e.g., [61]). Based on Hf and Nd whole-rock isotope data, the ca. 3.2 Ga old AGC gneisses are interpreted to have formed by melting of mafic lower crust and variable assimilation of older ca. 3.5 Ga felsic crust, resulting in variable enrichments of incompatible elements in these rocks (e.g., [61]).

The chosen potential crustal contaminants are distinct in their concentrations of incompatible elements. Tonalite AGC 449 has high Th/Yb, high Nb/Yb, a low Zr concentration, and high La/Yb, whereas granodiorite AGC 473 has lower Th/Yb at elevated Nb/Yb, a higher Zr concentration, and high La/Yb [61]. These two granitoid gneisses sample the heterogeneity of ca. 3.2 Ga AGC crust and represent two distinct end members in our modeling (Figure 11(a)). With this approach, two different AFC trends for Nsuze and Mozaan andesites and basaltic andesites in Th/Yb vs. Nb/Yb space can be reproduced. The upper AFC trend in Figure 11(a) (red line) is described by a maximum of 18% AFC of AGC 449, whereas the lower trend in Figure 11(a) (black line) can be reproduced by maximum 12% AFC for AGC 473. Again, the rhyolites of the Nsuze Group likely represent partial melts of continental crust. In contrast, in Zr vs. La/Yb space (Figure 11(c)), the AFC modeling does not reproduce all Nsuze and Mozaan lavas indicating that the chosen gneisses do not fully match the composition of the actual crustal contaminant. Additionally, the calculated percentages of AFC for Zr vs. La/Yb are higher than for Th/Yb vs. Nb/Yb which could be due to an overestimation possibly caused by small amounts of magma mixing, because Zr is more sensitive to magma mixing than Th.

On the other hand, our modeling for the Hlagothi, Usushwana, and Kubuta units shows that parental magmas derived from a depleted mantle source likely underwent AFC processes involving a crustal contaminant that has a low concentration of incompatible over more compatible elements and a high Zr concentration. In our modeling, this crustal contaminant is represented by the 3534 Ma tonalitic gneiss sample AGC 200 from the Ngwane Gneiss Unit (Figure 11(b)) [20]. In Th/Yb vs. Nb/Yb space, our calculations show that the Hlagothi rocks and our Kubuta samples can be reproduced by maximum 10% AFC and subsequent fractional crystallization of the parental magma shifting the composition along a parallel line to the MORB-OIB array to higher Th/Yb and Nb/Yb (Figure 11(b)). In addition, in Zr vs. La/Yb space, the Hlagothi and Kubuta rocks can be reproduced by maximum 12% AFC and subsequent fractional crystallization as indicated by the arrow in Figure 11(d).

In summary, our calculations show that the variable compositions of the volcanic rocks from the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta suites are largely the result of assimilation-fractional crystallization of either a komatiitic parental magma (Nsuze and Mozaan) or a parental magma similar to depleted mantle melt (Hlagothi, Usushwana, and Kubuta). However, magma mixing between a fractionated melt that already underwent assimilation on the one hand and a felsic crustal melt on the other hand can also account for compositions of some of the Nsuze and Mozaan Groups lavas. Furthermore, our modeling shows that the AGC crustal rocks from both, the ca. 3.5 Ga and the ca. 3.2 Ga crust-forming events are suitable contaminants to the aforementioned volcanic sequences.

This observation, however, is inconsistent with the εNd_(t) values of Nsuze, Usushwana, and Kubuta rocks that generally are about -2 to -3 (Figure 6). We calculated that at the time of formation of these volcanic units, i.e., at ca. 3.0 Ga, the ca. 3.2 AGC crust developed only to εNd_(3.0) values of about -3 to -4, whereas the ca. 3.5 Ga AGC crust had already developed to εNd_(3.0) values of up to -8 (Figure 6). Considering the low calculated amounts of AFC ranging from 10% to 18%, crustal assimilation of the younger generation TTGs seems to be unlikely. Hence, taking the isotope compositions into account, it is more likely that ca. 3.5 Ga old AGC crust was assimilated, that had comparable elemental compositions to the ca. 3.2 Ga old crustal rocks which were used for modeling.

An alternative explanation for the enriched trace elemental composition in the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta lavas could be a lithospheric mantle source beneath the Eastern Kaapvaal Craton that had previously been overprinted metasomatically, similar to what has been proposed for the lithospheric mantle beneath the Pilbara Craton (e.g., [85]). As a consequence, partial melting of metasomatized lithospheric mantle or plume-interaction with such a metasomatized mantle would produce similar geochemical characteristics to crustal contamination, i.e., high Th/Nb and elevated LREE (e.g., [85, 86]), as observed in the aforementioned volcanic sequences.

The addition of fluids that caused metasomatism of the lithospheric mantle prior to magma formation could have occurred in two ways, either by ancient subduction processes during crust formation (e.g., [22, 29, 87]) or by a process analogous to sagduction such as proposed for the Pilbara Craton (e.g., [85, 86, 88, 89]). These two possibilities are controversially discussed and determining which of the two possibilities is more plausible for a possible metasomatic overprint of the lithospheric mantle beneath the Kaapvaal Craton is beyond the scope of this paper. Nevertheless, we cannot rule out the probability that partial melting of a metasomatized lithospheric mantle could have produced the observed enrichment in the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta volcanic rocks.

Another way to test the relationship of our Kubuta samples to the Nsuze Group or the supposed large igneous province consisting of Mozaan, Hlagothi, and Usushwana rocks (e.g., [6]) is to constrain the possible mantle sources for these units. The coherence in elemental anomalies of primitive mantle-normalized trace element patterns of our Kubuta samples compared to the Nsuze, Mozaan, Hlagothi, and Usushwana rocks (Figures 5(a) and 5(b)) suggests that all these volcanic rocks could be genetically related. We note that in Th/Yb vs. Nb/Yb space (Figure 9), the Kubuta samples, together with the samples from the other units, plot consistently above the modern MORB-OIB array towards higher Th/Yb. As discussed in the previous chapter, the high Th/Nb and LREE enrichment could either result from contamination of the parental magmas of these rocks with continental crust or, alternatively, from the presence of metasomatized lithospheric mantle beneath the underlying continental basement.

In previous studies on the volcanic rocks of the Nsuze Group, different types of parental magmas for the evolved rocks were suggested. For example, Armstrong et al. [9] proposed a basaltic parental magma for Nsuze volcanic rocks derived from a primitive mantle source comparable to a garnet-lherzolite, whereas Nhleko [10] suggested that the basaltic parental magma was derived from a mixture of depleted mantle and enriched mantle melts. In contrast, Hegner et al. [5] proposed that the Nsuze and Usushwana volcanic rocks were derived from an ultramafic komatiitic parental magma sourced from a depleted mantle reservoir.

In the previous section, we showed that the Kubuta, Nsuze, Mozaan, Hlagothi, and Usushwana volcanic rocks lie on assimilation-fractional crystallization trends involving different crustal components and at least two mantle source compositions. As a result of our modeling, the parental melts of the Nsuze and Mozaan volcanic rocks appear to be derived from komatiitic melts, similar to Barberton komatiites. These results support the hypotheses of Armstrong et al. [9] and Hegner et al. [5] for the origin of Nsuze volcanic rocks.

In contrast, all other samples from the Hlagothi, Usushwana, and Kubuta units do not appear to be derived from primitive mantle melts. Instead, these more mafic to ultramafic rocks seem to originate from a more depleted mantle source (Figure 11(a)). This feature indicates that despite the coherence of primitive mantle-normalized trace element patterns, the Kubuta samples together with the volcanic rocks from the Hlagothi and Usushwana units derived from a mantle reservoir distinct from that of the volcanic Nsuze and Mozaan Group rocks. This observation contradicts the conclusions of Hegner et al. [5] who proposed a komatiitic parental magma for Usushwana volcanic rocks, but supports the possible existence of a common mantle reservoir for the Hlagothi, Usushwana, and Kubuta units, as proposed by Gumsley et al. [6]. However, in our modeling we use a modern depleted mid-ocean ridge basalt (MORB) mantle source [74] which might not be representative of the depleted mantle in the Archean. However, taken as a rough approximation, our results require an Archean mantle comprising reservoirs that were at least as depleted as modern MORB mantle. These could either be represented by ambient depleted mantle or, alternatively, depleted subcontinental lithospheric mantle. Furthermore, different mantle reservoirs for Hlagothi, Usushwana, and Kubuta volcanic rocks on the one hand, and Nsuze and Mozaan Group rocks on the other hand, would argue for a more complex formation history of the large igneous province than previously proposed by Gumsley et al. [6] based on age dating.

Petrogenetic models for Nsuze, Mozaan, Hlagothi, and Usushwana volcanic rocks have been proposed earlier considering zircon and baddeleyite age dating, as well as common geochemical characteristics. These petrogenetic models invoke a plume origin for the aforementioned volcanic rocks with at least two pulses of magma generation (e.g., [4–6, 46]). In these models, the first pulse of magma produced the parental magmas of the Nsuze lavas and the Usushwana Igneous Complex, whereas a later mantle plume event produced the parental magmas of the Mozaan lavas and the Hlagothi Complex [4, 6]. Based on similar ages for the Usushwana and Nsuze suites and the more mafic composition of the Usushwana unit, Hegner et al. [5] and Gumsley et al. [4] interpreted the Usushwana dikes as being feeder dikes for the Nsuze volcanic rocks. Likewise, Gumsley et al. [6] suggested that the similar trace element patterns of the synchronous Hlagothi and Mozaan units together with the more mafic composition of Hlagothi rocks are indicative of the Hlagothi Complex being the feeding system for the Mozaan volcanism. Based on their observations, Gumsley et al. [6] and Gumsley et al. [4] proposed that the sills of the Hlagothi Complex and the grabbros of the Usushwana Igneous Complex together with the extrusive Mozaan lavas form part of a large igneous province in the southeastern part of the Kaapvaal Craton. This large igneous province is proposed to extend into the western Witwatersrand basin, based on synchronous volcanic sequences, dikes, sills, and layered complexes that are found on the eastern Witwatersrand block [4, 6]. A relation of synchronous volcanic rocks of the Mozaan Group and Witwatersrand Supergroup has also been suggested by Mukasa et al. [11] and recently confirmed by Paprika et al. [8] by a stratigraphic correlation based on zircon U-Pb ages with the Dominion group.

In the present study, we investigated the geochemistry and probable relation of Nsuze, Mozaan, Hlagothi, Usushwana, and the newly sampled Kubuta volcanic rocks in more detail. At first glance, similar trace element and Nd isotope compositions (where available) suggest a common mantle source and a common crustal contamination process for the aforementioned volcanic sequences. Our assimilation-fractional crystallization modeling of Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta volcanic rocks showed that the lavas from the Nsuze and Mozaan Groups were derived from a common parental melt similar to komatiite from a weakly depleted mantle source. In contrast, the volcanic rocks from the Hlagothi, Usushwana, and Kubuta suites were derived from a distinct, much more depleted mantle source.

Our findings challenge the hypothesis of Hegner et al. [5], Gumsley et al. [6], and Gumsley et al. [4] that the Usushwana Igneous Complex and Hlagothi Complex would represent feeder systems for Nsuze and Mozaan lavas, respectively. Instead, our results suggest that the Usushwana IC and the Hlagothi Complex were derived from a more incompatible element-depleted mantle source than the extrusive rocks of the Nsuze and Mozaan Groups. However, a synchronous formation of the Usushwana IC and Nsuze lavas on the one hand and the Hlagothi Complex and Mozaan lavas on the other hand is supported by precise age determinations [4–6]. Hence, we propose a refined petrogenetic model for Nsuze, Usushwana, Mozaan, Hlagothi, and the newly sampled Kubuta volcanic rocks (Figure 12) combining our findings and the conclusions of Hegner et al. [5], Gumsley et al. [6], and Gumsley et al. [4]

We propose that the komatiitic parental magma to the Nsuze lavas was sourced from an ascending mantle plume leading to the first pulse of magma generation at around 2980 to 2970 Ma. The plume likely stagnated beneath the lithospheric mantle and mafic underplating of the komatiitic magma occurred between the lithospheric mantle and the felsic TTG crust that ultimately sourced the komatiitic parental magmas to the Nsuze lavas. Contemporaneously, the heat of the mantle plume likely caused partial melting of incompatible element-depleted ambient mantle or depleted subcontinental lithospheric mantle possibly by conductive heating. These partial melts represent the parental melts of the Usushwana IC. Subsequently, both, the komatiitic and depleted mantle melts filled magma chambers in the felsic AGC crust and underwent AFC processes that reduced the density of the magmas and the evolved melts could erupt (Figure 12(a)). After this first magmatic event, thermal subsidence led to the deposition of Mozaan Group sediments in a shallow marine environment, likely an epicontinental sea [4, 18]. In a second event at around 2860 Ma, magmatism was reactivated and plume-sourced komatiitic parental melts produced the parental magmas of the Mozaan lavas. Analogous to the first magmatic event, the hot komatiitic plume caused partial melting of the incompatible element-depleted ambient mantle or, alternatively, the depleted subcontinental lithospheric mantle which resulted in the production of the parental magmas of the Hlagothi Complex. Again, these parental melts filled magma chambers in the felsic AGC crust and underwent AFC prior to eruption (Figure 12(b)).

In our refined petrogenetic model, the Kubuta rocks represent differentiation products of depleted mantle melts similar to the lavas of the Usushwana Igneous Complex or Hlagothi Complex. Due to the poorly constrained age of the Kubuta samples, it is currently not possible to distinguish if these rocks formed in association with the first or second pulse of magma generation.

The newly sampled metavolcanic rocks from the vicinity of the Kubuta Ranch in southern Eswatini resemble the major and trace element compositions of the Hlagothi and Usushwana Igneous Complexes, having similar trace element characteristics that indicate combined mineral fractionation and crustal assimilation. Also, our results show that the Nsuze and Mozaan volcanic rocks of the Pongola Supergroup were formed from a komatiitic parental magma. When compared to the Nsuze and Mozaan lavas, the Kubuta, Hlagothi, and Usushwana suites are more mafic to ultramafic in composition. In addition, assimilation and partial melting of continental crust is indicated in Nsuze, Usushwana, and Kubuta volcanic rocks by generally negative εNd_(t) values and negative Nb-Ta anomalies.

Assimilation-fractional crystallization modeling shows that the lavas of all volcanic suites (Nsuze, Mozaan, Usushwana, Hlagoti, and Kubuta) result from modification of their respective parental magmas by AFC processes involving Ancient Gneiss Complex crust of diverse trace elemental composition. The heterogeneity of AGC crust is represented in the crustal end members used for modeling that include samples of the 3.2 Ga and 3.5 Ga crust-forming events. However, based on the negative εNd_(t) values of about -3 in Nsuze, Usushwana, and Kubuta lavas and the calculated low amounts of AFC, we consider it more likely that AGC crust of the 3.5 Ga crust-forming event was assimilated. However, melting of a metasomatized lithospheric mantle or plume-lithospheric mantle interaction may also produce magmas with the observed geochemical characteristics.

Based on our results and previous studies, we propose a refined petrogenetic model for the formation of Nsuze, Mozaan, Hlagothi, Usushwana, and the newly sampled Kubuta volcanic rocks. In this model, we follow the interpretation of Gumsley et al. [6] that the aforementioned volcanic suits formed a Mesoarchean large igneous province that probably extended into the eastern Witwatersrand block, including the Dominion group [8]. However, unlike Gumsley et al. [6], we also assign the Nsuze lavas to belong to this large igneous province. Furthermore, our results show that the Nsuze, Mozaan, Hlagothi, Usushwana, and Kubuta lavas were not derived from a single mantle plume, as proposed by Gumsley et al. [6], but rather were derived from two distinct sources—a plume producing komatiitic magmas (Nsuze and Mozaan) on the one hand and partial melting of depleted mantle (Hlagothi, Usushwana, and Kubuta) triggered by the conductive heat of the mantle plume on the other hand.

To test our hypothesis of two different mantle sources and, additionally, to constrain the contamination history in more detail to be able to distinguish if 3.2 Ga or 3.5 Ga AGC crust was preferentially assimilated, more geochemical work on the Mozaan, Hlagothi, and Usushwana and Dominion volcanic rocks is needed in form of major and trace element as well as whole-rock Hf-Nd isotope analyses. Furthermore, precise age determinations of the Kubuta volcanic rocks are crucial to affiliate these rocks to either the first or second pulse of magma generation.

Nevertheless, our findings strengthen the geochemical fingerprint of the proposed large igneous province [8] and provide information to better constrain the petrogenesis of these Mesoarchean flood basalts that record an important crustal growth process in the eastern Kaapvaal Craton at ca. 2.9 Ga.

All data are available in the manuscript and supplementary material. Raw data can be provided on request from the corresponding author.

The authors declare that there is no conflict of interest regarding the publication of this paper.

JEH, KPS, CM, and AK thank the Deutsche Forschungsgemeinschaft for grants HO4794/1-1, HO4794/1-2, Mu1406/8, Mu1406/19-1, KR590/94-1, and KR590/94-2. Axel Hofmann is thanked for logistic support in the field and sample shipping. We thank Philipp Gleißner, Yogita Kadlag, and Chunhui Li (Berlin) for support with the ICP-MS as well as Frank Wombacher and Almut Katzemich for lab support at Cologne. Monika Feth is thanked for laboratory assistance at Berlin. This is a contribution within the framework of the SPP1833 Project ‘Building a habitable Earth’ funded by the DFG. We acknowledge support by the Open Access Publication Fund of the Freie Universität Berlin.