Probing the Southern African Lithosphere With Magnetotellurics: 2. Linking Electrical Conductivity, Composition, and Tectonomagmatic Evolution

The tectonic history of southern Africa includes Archean craton formation, multiple episodes of subduction and rifting and some of the world's most significant magmatic events. These processes left behind a compositional trail that can be observed in xenoliths and measured by geophysical methods. The abundance of kimberlites in southern Africa makes it an ideal place to test and calibrate mantle geophysical interpretations that can then be applied to less well‐constrained regions. Magnetotellurics (MT) is a particularly useful tool for understanding tectonic history because electrical conductivity is sensitive to temperature, bulk composition, accessory minerals, and rock fabric. We produced three‐dimensional MT models of the southern African mantle taken from the SAMTEX MT data set, mapped the properties of ∼36,000 garnet xenocrysts from Group I kimberlites, and compared the results. We found that depleted regions of the mantle are uniformly associated with high electrical resistivities. The conductivity of fertile regions is more complex and depends on the specific tectonic and metasomatic history of the region, including the compositions of metasomatic fluids or melts and the emplacement of metasomatic minerals. The mantle beneath the ∼2.05 Ga Bushveld Complex is highly conductive, probably caused by magmas flowing along a lithospheric weakness zone and precipitating interconnected, conductive accessory minerals such as graphite and sulfides. Kimberlites tend to be emplaced near the edges of the cratons where the mantle below 100 km depth is not highly resistive. Kimberlites avoid strong mantle conductors, suggesting a systematic relationship between their emplacement and mantle composition.

the presence of extensive outcropping rocks and mantle xenolith-bearing kimberlites (e.g., Griffin et al., 2003;Jelsma et al., 2004), the Kalahari Craton is a great natural laboratory to understand craton formation and survival as well as plate tectonic and magmatic processes throughout geological time. These rocks and xenoliths strongly indicate that the mantle beneath Kalahari Craton composition is highly variable (e.g., Griffin et al., 2003;Humbert et al., 2019). Variations in mantle composition are either a result of the age-dependent thermal state and composition of the mantle as it initially formed (e.g., Griffin et al., 1999;Pearson et al., 2004) or subsequent alterations imposed by melts and fluids that infiltrated the lithospheric mantle (Alard et al., 2000;Griffin et al., 2003). The Kalahari Craton and its surroundings also hold crucial economic deposits such as the PGE-rich layers of the Bushveld Complex (VanTongeren, 2017) and diamond-bearing kimberlites (Jelsma et al., 2004;A. G. Jones et al., 2009), which formed as the result of lithospheric-scale processes (Begg et al., 2010;Griffin et al., 2013). Therefore, studying lithospheric composition may improve models for economic geology as well as continental evolution.
The magnetotelluric (MT) method is a geophysical technique that images the subsurface electrical conductivity structure of the Earth to upper mantle depths. Electrical conductivity can provide knowledge on bulk composition and temperature as well as the presence of interconnected accessory materials (e.g., fluids, melt, hydrous minerals, and sulfides). Previous studies have shown that the cratonic mantle has highly variable electrical conductivity that cannot be accounted for by temperature differences alone (e.g., Evans et al., 2011;Selway, 2015). These conductivity variations can partly be ascribed to hydrogen species structurally bound to the nominally anhydrous minerals (NAMs) that constitute the bulk of the upper mantle, such as olivine, pyroxenes, garnet, and spinel (Demouchy & Bolfan-Casanova, 2016).
In this article, we investigate the relationships between the Kalahari Craton's tectonic setting, thermal structure, magmatic events, and metasomatic signatures with its geoelectric structure. We use the continental-scale Southern African Magnetotelluric Experiment (SAMTEX) data set (A. G. Jones et al., 2009) to produce 3D MT models of the mantle, and we estimate the composition, metasomatic signatures, and thermal structure from analyses made on 36,066 garnet xenocrysts taken from southern African kimberlites. This article represents the second part of a two-part study. In the first part (Moorkamp et al., 2021), we investigated 3D MT modeling of the SAMTEX data set using different inversion algorithms (Kelbert et al., 2014;Moorkamp et al., 2011) and aimed to understand the effects of strategies used to model the mantle conductivities. In this article, we rely mainly on the inversion of the "selected data" data set run with ModEM from a median starting model described in that work, while being mindful of the robustness of these features implied by different modeling attempts. We refer any reader that is interested in the MT modeling aspect of this study to the first part (Moorkamp et al., 2021).

Tectonic Evolution of the Kalahari Craton and Surrounding Terranes
The Kalahari Craton consists of terranes with Archean to Neoproterozoic basement ages, namely the Zimbabwe Craton, Kaapvaal Craton, Limpopo Belt, Kheis Belt, Magondi Belt, and Rehoboth Terrane (Figure 1). Many of the older basement terranes are hidden under the cover of Neogene Kalahari Group sediments and thick marine-sedimentary sequences of the Jurassic Karoo Group. Due to this cover, most geological understanding of the Rehoboth Terrane and northwestern Kalahari Craton is derived from geophysical and borehole data (e.g., Chisenga et al., 2020;Corner & Durrheim, 2018). In contrast, the basement rocks of the Kaapvaal and Zimbabwe cratonic nuclei and their immediate surrounding blocks in the eastern Kalahari Craton are better exposed ( Figure 1) and have been the subject of rigorous geological and geophysical study for many decades (Corner & Durrheim, 2018;De Beer, 2016;Oriolo & Becker, 2018).
The basement of the Kaapvaal Craton is dominated by 3.6-2.9 Ga gneiss, granitoid, and greenstone belts (Poujol et al., 2003). The orientations of the greenstone belts are markedly different in the Kimberley Block in the western Kaapvaal Craton and the Witwatersrand Block in the eastern Kaapvaal Craton, separated by the Colesburg Magnetic Lineament (Figure 1), which suggests the independent formation of these Archean terranes before amalgamation into a single craton (Jacobs et al., 2008). In the northern Kaapvaal Craton, the Pietersburg Terrane is thought to have accreted to the Witwatersrand Block at 2.73-2.65 Ga (Laurent et al., 2019). This collision created the Thabazimbi-Murchison Lineament (TML), a trans-lithospheric structure observed in aeromagnetic data (Good & De Wit, 1997). Assembly of the Kaapvaal Craton was accompanied by the formation of the foreland basin (3.0-2.78 Ga), the rift-related Ventersdorp Volcanics (2.79-2.65 Ga, Gumsley et al., 2020) and deposition 3 of 28 of Transvaal Supergroup (2.6-2.058 Ga, Zeh et al., 2020). The final amalgamation of the Pietersburg Terrane with the Witwatersrand Block occurred during the collision between the Kaapvaal Craton and the Zimbabwe Craton to the north. The Limpopo microcontinent was wedged between the two cratons during this collision (2.71-2.67 Ga, Laurent et al., 2019), giving rise to the Limpopo Orogeny and the formation of the Limpopo Belt.
Accretion of the Paleoproterozoic belts and blocks started around 2.06 Ga, corresponding to the timing of the Bushveld and syn-Bushveld magmatism along the TML (Molopo Farms Complex, Okwa Complex, Zeh et al., 2015) and amalgamation of Archean terranes (Laurent et al., 2019;Oriolo & Becker, 2018). The Magondi Belt accreted onto the northwestern end of the Zimbabwe Craton through subduction processes (Jacobs et al., 2008;Master et al., 2010) and then experienced coeval transpressional deformation and metamorphism with the Limpopo Block in response to a collision with an unknown terrane from the northwest (Oriolo & Becker, 2018). Contemporaneously, Kheis Belt sediments were deposited onto the western end of the Kaapvaal Craton and the region subsequently collided with the Rehoboth Terrane (Jacobs et al., 2008) along the prominent Kalahari Magnetic Lineament (Corner & Durrheim, 2018). The tectonic history and basement geology of the Rehoboth Terrane are somewhat enigmatic due to the thick sedimentary sequences covering the surface, especially in the central region. Geochemical and geochronological data from the westernmost inliers suggest that the Rehoboth Terrane  (Chisenga et al., 2020;Corner & Durrheim, 2018;Hanson, 2003;McCourt et al., 2013 Following Paleoproterozoic evolution and cratonization, the assembled Kaapvaal-Limpopo-Zimbabwe-Rehoboth Block was behaving as a rigid block by Mesoproterozoic (Jacobs et al., 2008). After stabilization, the region was impacted by 1.4-1.35 Ga intraplate alkaline magmatism (Pilanesberg Complex and Premier kimberlite, Hanson et al., 2006). This was followed by more abundant and widespread magmatism concentrated near the northwestern border of the craton associated with the ∼1.1 Ga Umkondo Large Igneous Province, an event thought to be related to rifting (De Kock et al., 2014;Hanson et al., 2006) that used the lithospheric-scale weakness zones formed during assembly of the proto-Kalahari Craton (e.g., Xade Complex, Hanson et al., 2006). The Umkondo-aged rhyolitic units (Kgwebe Formation and southwestern correlatives) mostly outcrop along a ridge in the Ghanzi-Chobe Belt stretching from northeastern Botswana to central Namibia. This ridge was uplifted in response to the Neoproterozoic Pan-African Orogeny, which also formed the Damara Belt (Modie, 2000). To the south, the Kalahari Craton is bounded by the Proterozoic Namaqua-Natal Belt, which formed as an assemblage of numerous microcontinental terranes in a convergent setting. The Namaqua-Natal Belt collided with and accreted onto the Kalahari Craton during series of tectonic events between 1.2 and 1.0 Ga (Jacobs et al., 2008).
The Jurassic breakup of Gondwana had a significant impact on the Kalahari Craton, including the emplacement of widespread rift-related Karoo units (e.g., Drakensberg Lavas, Karoo diorite sills, Okavango and Save-Limpopo Dyke Swarms, Svensen et al., 2012). Following the Karoo event, extensive Group II (∼110-127 Ma) and Group I (∼110−72 Ma) kimberlites were emplaced in the Kalahari Craton. Differences in mantle xenolith and xenocryst compositions between kimberlites from these two time windows indicate an intervening major metasomatic event in the mantle (Kobussen et al., 2009).

Magnetotelluric Data and Modeling
In the first part of this two-part study (Moorkamp et al., 2021), we explored the effects on 3D MT modeling of the SAMTEX data set imposed by data selection, regularization, starting model, and the inversion algorithms ModEM (Kelbert et al., 2014) and jif3D (Moorkamp et al., 2011). Results demonstrated that modeled mantle conductivities could be affected by the regularization schemes. For instance, ModEM tends to converge toward the initial model, whereas the jif3D model remains relatively constant as the data sensitivity becomes poorer. The model presented here was produced using the ModEM algorithm (Kelbert et al., 2014) and the sparsely selected, good-quality data. We choose this over the jif3D model since the inversion code was more accessible to the primary authors of this study and to the broader MT community but the features on which we base our interpretations are common to all models. The model setup and robustness of the conductivity features are described in the companion study (Moorkamp et al., 2021). Given our finding that the ModEM regularization biases the result toward the initial model, we base this analysis on the ModEM inversion run from an initial model constructed from the median station apparent resistivities, as also described in the companion article. Depth and vertical cross-sections from the inversion result are shown in Figures 2 and 3, respectively. A 3D contour plot of the models is shown in Figure 4.
Compositional interpretation of mantle conductivity requires sensitivity tests to be carried out to estimate absolute resistivity values from MT models. We did this by selecting the areas of interest for calculating water contents and other compositional parameters, replacing the modeled resistivities in these areas with blocks of different resistivity values, and testing the impact on data misfit (Figures S1-S5 in Supporting Information S1).

Garnet Xenocryst Analyses
Analyses were carried out on existing data collected from 36,066 garnet xenocryst samples taken from Group I kimberlites around southern Africa (Table S1 in Supporting Information S1). Garnet xenocrysts are a valuable resource used in exploring mantle composition since they are more abundant than much rarer mantle xenoliths (Ryan et al., 1996). Although they do not provide the range of compositional and textural information available from xenoliths, the ability to analyze a large number of samples allows extensive investigation of the lithospheric mantle beneath. The data include major-element analyses by electron microprobe and trace-element data collected using with laser-ablation ICPMS techniques, in both the GEMOC ARC National Key Centre, Macquarie University, and the DeBeers Group Services Laboratory, Johannesburg, South Africa (Kobussen et al., 2009(Kobussen et al., , 2008. To make depth-dependent classifications, we first performed thermobarometry on these samples. We estimated pressure and temperature using the following steps: (a) max − conditions were calculated with nickel-in-garnet thermometry and the chromium solubility barometer of Ryan et al. (1996). (b) A generalized cratonic geotherm of Hasterok and Chapman (2011) was fitted based on the locus of maximum P Cr at each T Ni . (c) Temperature at the base of the depleted lithosphere (BDL; T BDL ) was defined as the temperature at which the proportion of garnets with ≤10 ppm (wt) yttrium decreases sharply (e.g., Figure S8 in Supporting Information S1). This decrease in garnets with low Y contents indicates a sudden decrease of depleted garnets since Y-in-garnet is highly correlated with melt-related metasomatism (O'Reilly & Griffin, 2006). The BDL does not therefore equate to the base of the thermal, seismic or rheological lithosphere, but is the base of the geochemically depleted lithosphere that has not been substantially melt metasomatized. (d) The inflection point of the geotherm was defined as T BDL , the temperature at the BDL. (e) At depths greater than BDL temperatures estimated from the max − thermometer are less reliable. At these temperatures, a kinked geotherm was defined as a line parallel to the trend of the diamond-graphite transition (Day, 2012), based on the commonly observed distribution of max − parallel to this trend (Griffin et al., 2003). While the temperature estimates may not be entirely correct after this point, we interpret our results on depths shallower than the BDL. (f) Finally, calculated equilibrium temperatures of the garnet grains were projected onto the defined geotherm to determine their depth of origin. The calculated geotherm parameters for each pipe are given in Tables S1 and S2 in Supporting Information S1, and fitted generalized cratonic geotherm surface heat flow (SHF) values are mapped in Figure 5b.
For the garnet xenocryst classifications depicted in Figures 6 and 7, we used the method of cluster analysis by regressive partitioning (CARP, Griffin et al., 2002). The CARP method was designed to exploit the abundance of garnet mantle xenocryst samples by classifying Cr-pyrope garnets into statistically significant populations of Information S1), which Griffin et al. (1999) suggested to be a feature related to specific Archean metasomatic processes. Depleted lherzolites with phlogopite metasomatism (blue) have relatively depleted major element compositions but trace-element (Ti, Zr) signatures characteristic of phlogopite crystallization. The most populous metasomatic class, fertile lherzolites (green), represents rocks with compositions enriched in major and trace elements with more diverse characteristics. These could be rocks that never experienced depletion, or depleted Archean material that was later refertilized. While depleted and refertilized signatures are usually associated with the Archean mantle, signatures of the mantle that has never been depleted are usually associated with Proterozoic or younger mantle (Griffin et al., 2002). Samples from the "melt-metasomatized" class (red) are associated with very enriched, high-temperature lherzolites. They commonly show sheared microstructures which indicate melt infiltration into the rock (Griffin et al., 2003) and are largely located below the BDL. Whole-rock Al 2 O 3 contents used in Figure 6a are calculated from regression analyses made on yttrium-in-garnet (O'Reilly & Griffin, 2006), while Mg ol # is calculated from the garnet data by the method described in Gaul et al. (2000). REE N patterns of the CARP classes mostly demonstrate significantly different characteristics. The sinuosity of the REE patterns is associated with depleted material in the mantle (Griffin et al., 1999) and can be quantified by using the log-Nd N /Dy N ratio, in which values above zero indicate sinuous patterns and values below zero indicate less sinuous patterns. Yb N , on the other hand, is used as a proxy of overall HREE enrichment. Metasomatic CARP classes exhibit less sinuous patterns, heavy REE (HREE) enrichment, and light REE (LREE) depletion. During refertilization, garnet REE N trends become less sinuous (lower log-Nd N /Dy N ), with higher HREE N /LREE N ratios, as the metasomatic fluids percolate more extensively and equilibrate with the environment (Griffin et al., 2003). When the mean values of log-Nd N /Dy N and Yb N are plotted against depth (Figures 9b and 9c), they are correlated (negatively and positively, respectively) with the population of total metasomatic classes within the fertile layer (usually ∼100-140 km). Therefore, we use these parameters as a proxy for the lateral extent of metasomatic fluid percolation in the lithospheric column, or in other words, for the intensity of metasomatism at the corresponding depth.

Methods for Determining the Links Between Garnet Xenocryst and MT Data
Electrical conductivity can be used to infer the composition of the mantle (e.g., Karato & Wang, 2012;Selway et al., 2019). Compositional interpretation can be achieved by integrating experimental conductivity and petrology studies, thermal structure, and phase-mixing models (Özaydın & Selway, 2020). Primarily, the electrical conductivity of mantle minerals depends on temperature through semi-conduction processes. Some additional materials such as structurally bound hydrogen (expressed below as water) can change the thermal energy required to enhance conductivity (e.g., Dai & Karato, 2014;Wang et al., 2006). The extent to which water can be incorporated into mantle minerals is limited (e.g., Férot & Bolfan-Casanova, 2012;Padrón-Navarta & Hermann, 2017) and in some cases, the measured mantle conductivities may require an additional conductive phase, either a mineral with a higher activation enthalpy (e.g., phlogopite) or a very conductive mineral with low-temperature dependence (e.g., graphite, sulfides, Watson et al., 2010;Zhang & Yoshino, 2017).
We used the program MATE (Özaydın & Selway, 2020) to investigate such relationships. For the water calculations, we used the olivine water-partitioning coefficients of Demouchy et al. (2017) and Novella et al. (2014) for pyroxenes and garnet, respectively. We chose to seek solutions of water content up to limits determined by the olivine solubility model of Padrón-Navarta and Hermann (2017) since it reflects the near-pure H 2 O state of the cratonic mantle in subsolidus conditions. Water contents were modified for both water-solubility and electrical conductivity models to reflect the calibrations of Withers et al. (2012) for olivine and Bell et al. (1995) for pyroxenes and garnet. Our figures show the water content results using the olivine conductivity model of Gardés et al. (2014). Different selections of pyroxene and garnet electrical conductivity models do not make considerable differences (Özaydın et al., 2021). We chose to use the conductivity models of Dai and Karato (2009b), Y. Li   Zhang et al. (2012) for orthopyroxene, clinopyroxene, garnet, and phlogopite, respectively. For the mixing model, we used the Modified Archie's model (Glover, 2010). Olivine was set to be a perfectly connected matrix (m << 1), while interconnectivity of orthopyroxene was set to m = 2.5. Clinopyroxenes and garnet were set to be not connected with a value m = 4.
A recent study made detailed comparisons of xenolith water measurements and MT-derived water calculations for the Kimberley-Jagersfontein region (Özaydın et al., 2021) and showed that water contents measured from mantle xenoliths broadly match those interpreted from MT models with similar configurations on experimental parameters. Since there are no water content measurements made outside this region in southern Africa, in this work we focus on the trends of the modeled water contents rather than the specific water contents interpreted from the data. We show the resistivity-depth sections and relevant water content calculations at seven locations near the kimberlite clusters at Figure 8.

Lithospheric Architecture From Magnetotelluric Models
MT models produced in this study demonstrate highly variable mantle conductivities across Precambrian terranes with different ages (Figure 4). Archean cratons (Kaapvaal and Zimbabwe cratons) are depicted as complex regions of mostly resistive lithosphere carved by conductive features, reflecting their metasomatic history. In the Kaapvaal Craton (Figure 3a), this is exemplified by the contrast between the conductive mantle beneath the Bushveld Complex and the resistive center of the Kaapvaal Craton (C 1 -C 2 , Figure 3a), which appears to have a convex shape very similar to the mantle keels modeled (e.g., Afonso et al., 2008) and imaged by seismic tomography (e.g., Fouch et al., 2004;A. Li & Burke, 2006;Ortiz et al., 2019). North of the Bushveld Complex, at the northern end of the Kaapvaal Craton, the Archean Pietersburg Terrane also appears as a deep resistive feature. Still further north, the Archean Limpopo Belt consists of conductive mantle starting around ∼50 km depth (LC) and a complex crustal assemblage . The distribution of electrical conductivity in the MT models in this region generally resembles previous 2D models  and 3D models made in Limpopo Belt , while being modestly different in terms of absolute resistivities and small-scale features in the crust.
One of the MT profiles (Z 1 -Z 2 , Figure 3h) crosses the Archean Zimbabwe Craton on its proposed southwestern edge. The same profile previously modeled with 2D methods observed a complex crustal assemblage with a thick resistive root beneath the Zimbabwe Craton (Miensopust et al., 2011). These results are similar to ours at the mantle scale, while our model does not exhibit the same complexity in the crust due to larger mesh sizes and restricted frequency selection. Going toward the north, the Magondi Belt is modeled with a moderately conductive lower crust and mantle, bounded to the north by the thick resistive mantle of the Mesoproterozoic Ghanzi-Chobe Belt.
Due to the Neogene Kalahari sediments that cover its surface, the structure of the Archean-Proterozoic Congo Craton is not well known in most places. The craton is a well-defined resistive feature in the westernmost profile near the outcropping Kamanjab Inlier (K 1 -K 2 , Figure 3b). However, its high resistivity seems to be limited to specific sections within the proposed craton boundaries. Further east (Figure 3c), mantle resistivities become considerably lower in the vicinity of profile B 1 -B 2 and then increase to become again more keel-like in the vicinity of profile D 1 -D 2 (Figure 3e). A north-south striking conductive feature at lower lithospheric mantle depths beneath the central Congo Craton is hinted at in this model but is more prominent in inversions run with more of the MT stations included (Moorkamp et al., 2021) and is similar to a structure modeled in S and P wave tomography studies. Those studies also show that the region surrounding the Kamanjab Inlier has the highest velocities; the lowest velocities are in the central craton and more moderate velocities are modeled near the eastern margin of the craton (White-Gaynor et al., 2020). These results suggest that the Congo Craton, as it is often mapped  ( Figures 3 and 4), might be a fragmented tectonic unit and may consist of either a complex tectonic arrangements of blocks of different lithospheric thicknesses or may contain relics of past magmatism beneath its Neogene cover.
The Proterozoic Reheboth Terrane is imaged as a fragmented feature in which the central parts (K 2 -K 3 , Figure 3b) have a more conductive lower lithosphere (>100 km), while the northern end nearing the Ghanzi-Chobe region and the Southern Gibeon Fields is modeled with a resistive lithosphere. Similarly, other Proterozoic regions between the Kaapvaal-Zimbabwe-Limpopo Craton and Congo Craton have similar attributes. The youngest of the mobile belts, the Damara Belt, is a Pan-African orogenic zone formed during the collision of the Kalahari and Congo cratons in the late Neoproterozoic and early Cambrian (Goscombe et al., 2017). As imaged by MT models, the Damara Belt consists of a complex crustal assemblage  and a prominent lower crustal and upper mantle conductor (Figures 3c, 3e, and 4). Variations in the mantle conductivity of Proterozoic regions suggests that different processes had different effects on the compositional evolution of the mantle. Comparison of CARP sections depicted as 10 km interval histograms and olivine Mg ol # with modeled mantle conductance (50-150 km). The relative proportions of metasomatic styles are indicated with different colored bars in the histograms, where red is melt metasomatized lithosphere, green is lithosphere that is fertile in major and trace elements, blue is lithosphere that is depleted in major elements but shows signatures of phlogopite precipitation, and bright and dark yellow are lithosphere that is depleted in major and trace elements with lherzolitic and harzburgitic compositions, respectively. BDL refers to the depth below which depleted garnet populations are greatly diminished. Locations of the kimberlite clusters are denoted on the map (a-p). SHF: surface heat flow value of fitted generalized cratonic geotherm (Hasterok & Chapman, 2011). Horizontal dotted line indicates the BDL.

Southern African Mantle Imaged by Garnet Xenocrysts
A selection of the garnet data from 100 to 150 km depth is depicted in Figures 5 and 6. The ages of the kimberlite eruptions ( Figure 5a) are taken from multiple sources (Tables S1 and S2 in Supporting Information S1). In situations where the age of the pipe was unknown and the geochemical stratification was very similar to the other pipes within the same cluster, we assumed the kimberlite to be roughly the same age as that cluster. The rest of the unknown ages are represented with black squares. The garnet xenocryst analysis illustrates the complexity of the style of metasomatism in southern Africa and emphasizes that, from region to region, metasomatic fluids had distinct compositions and/or that the initial formation processes of the mantle rocks differed. For instance, iron enrichment (lower Mg ol #) and Al 2 O 3 enrichment are not always correlated even though they can both regarded as proxies for the fertility of the mantle (Griffin et al., 2002).
Lithospheric CARP sections from selected clusters are depicted as histograms in Figure 7 alongside the mantle conductance map ( ∑ = ; i.e., integrated conductivity (σ) times thickness (d)) at 50-150 km depth derived from the MT inversion model. Histograms here show relative proportions of the rocks where garnet xenocrysts are derived from.
Some generalizations can be made from the existing garnet xenocryst database for Group I kimberlites: 1. Lower whole-rock Al 2 O 3 , as calculated from Y-in-garnet contents, is a feature of a depleted mantle and mostly appears in Archean regions. On the other hand, higher values point to either a formerly depleted and then refertilized mantle (e.g., Bushveld region) or a mantle that may never have been depleted (e.g., Gibeon, Figure 6a) 2. Similarly, higher Mg ol # values signify higher proportions of depleted material, while lower values of Mg ol # (<92) are more indicative of widespread melt metasomatism (Figure 6b). In contrast, higher whole-rock Al 2 O 3 values are more likely to be related to infiltration of lighter fluids rather than dense and iron-rich melts (Griffin et al., 2003). Therefore, lower values of Mg ol # usually correlate with a shallower BDL, and are more likely to be a feature of the thinner Proterozoic terranes with lithospheric thicknesses between 100 and 150 km (e.g., Okwa and Uintjiesberg clusters, Figure 7) 3. As expected, the total proportion of metasomatized xenocryst classes is higher in non-Archean terranes ( Figure 6c). One of the main differences between Archean and Proterozoic terranes is the distinct stratification trends observed in lithospheric CARP sections shown in Figure 7. Geochemical tomography sections of the Archean mantle usually show a "depleted layer" between roughly 120-170 km. In contrast, in Proterozoic terranes, this depleted material is both less abundant and/or spread throughout the whole depth range in lower proportions 4. The CaO − Cr 2 O 3 classification scheme (Grütter et al., 2004) gives a similar story; Archean units are marked by lower proportions of lherzolitic material compared to their Proterozoic counterparts (Figure 6d) 5. Ti/Zr can be used as a proxy for phlogopite metasomatism whereby medium to low Ti/Zr ratios generally indicate phlogopite-related metasomatism (Griffin et al., 2002, Figure 6e). Very high values of Ti/Zr in some Proterozoic areas (Gibeon, Okwa, Sikereti) might indicate that the fertile material at these locations was not affected by phlogopite-related metasomatism. The pipes with high Ti/Zr ratios mostly plot in the fertile fields on Ti-Zr plots. Most of the xenocryst material from which these fields were originally classified came from garnet lherzolites xenoliths in basalts from off-craton areas and are thought to reflect fertile mantle which never experienced a depletion event. The garnets plotted in fertile fields on Ti/Zr plots are typically a feature of younger mantle (Griffin et al., 2002) 6. "Unclassified" samples denote garnet analyses that cannot be grouped into the statistically significant types of metasomatism (Figure 6f). These garnets most likely reflect the effects of complex metasomatic overprinting, either through multiple episodes by compositionally different fluids. The very high proportions of unclassified garnets observed in the mantle beneath the Bushveld region may indicate such complex metasomatic processes and are associated with a prominent mantle conductor 7. Figures 6g and 6h show the mean chondrite-normalized Nd/Dy and Yb ratios of xenocrysts from the different "Fertile" classes. The meaning of the garnet REE patterns are explained in Section 3.2. Garnets with more sinuous REE and HREE-poor characteristics are more abundant near the core of the Kaapvaal Craton (Kimberley and Kuruman clusters), mirroring the amount of depletion observed. The least sinuous and HREErich areas appear on the edges of the cratons or areas immediately surrounding them (Uintjiesberg, Jwaneng) indicating the metasomatism in these areas was more extensive

Discussion
If we assume that the connection between composition and electrical conductivity of the mantle can be adequately established through experimental petrology studies, information from garnet xenocrysts and MT models should be compatible. Therefore, we developed interpretations of the southern African lithospheric mantle considering both the broad lithospheric architecture and case by case evaluations of mantle composition. This improves not only our interpretations of MT models, but also of the three-dimensional composition of the southern African lithosphere and our ability to make such interpretations in regions with poorer xenolith constraints. First, we will discuss the general relationships between lithospheric architecture inferred from MT models (Figures 2-4) and broad metasomatic signatures derived from garnet xenocryst data (Figures 6 and 7). Then, we will interpret magnetotelluric models near kimberlite clusters with petrological information from garnet xenocrysts sections on a case-by-case basis.

Resistivity of Archean, Proterozoic, Depleted, and Metasomatized Domains
Electrical conductivity distribution of the Archean terranes varies significantly. This suggests that metasomatic events left a compositional mark on the initially depleted lithosphere. Similarly, Proterozoic terranes assembled near the Archean cratons do not show consistent characteristics and architectures that would provide a simple relationship between lithospheric thickness, age, and composition from the electrical conductivity distribution. The south-central Rehoboth Terrane, for instance, demonstrates lower resistivities (K 2 -K 3 , Figure 3b) and higher velocity bodies (White-Gaynor et al., 2020). In contrast, the Gibeon area, which includes one of the clusters of Cretaceous kimberlitic volcanism, is imaged as a highly resistive region (Figure 3f), with relatively lower velocities (White-Gaynor et al., 2020). Similarly, the northern end of the Rehoboth Terrane and adjacent Ghanzi-Chobe belt possibly shows the greatest depths at which 1,000 Ωm resistivity is observed (Figure 12b). Such high resistivities have been previously imaged in Proterozoic regions surrounding older cratons (e.g., Selway, 2015;Selway et al., 2011) and their high resistivities have been interpreted as reflecting magmatic and plate tectonic processes that induced mantle melting and depletion. In southern Africa, the resistors around the Ghanzi-Chobe Belt line up with the 1.1 Ga Umkondo magmatic units (Modie, 2000, Figures 2 and 3g). This area was the focus of rifting and extensive melting during the breakup of the Rodinia Supercontinent (De Kock et al., 2014) which might have created the depleted resistive mantle with low volatile contents that we observe today; similar features are seen in models of the Mozambique Belt, a region of Pan-African aged magmatism and high-grade metamorphism and a currently active rift zone (Selway, 2015). If the relationship between rifting events and volatile extraction from the melting of the mantle can be more confidently linked, the MT models of the Ghanzi-Chobe belt may provide evidence for a rifting origin for the Umkondo magmatism (De Kock et al., 2014).

Kimberley
The Kimberley domain is rich in xenolith-bearing diamondiferous kimberlites from which we have 5,521 garnet xenocrysts to map the lithospheric mantle. Previous studies of garnet xenocrysts have demonstrated that the area was metasomatized by an event that occurred between the eruptions of Group II (120-180 Ma) and Group I (90-120 Ma) kimberlites (Griffin et al., 2003;Kobussen et al., 2009). Group II kimberlites are populous in this area, reflecting modal metasomatic modification of the mantle prior to kimberlite emplacement (Giuliani et al., 2015).
The Kimberley-Jagersfontein area is geochemically stratified into three layers: (a) A fertile layer between 90 and 130 km that is marked by lower Mg ol # values, (b) a depleted layer between 125 and 175 km, and (c) an underlying melt-metasomatized layer. Our MT results demonstrate that the Kimberley domain overlies a lithosphere with relatively low resistivities (the greatest depth at which 1,000 Ω m is modeled is ∼100-120 km, Figure 12b) compared to the main Kaapvaal resistor (>150 km, Figure 12b). The conductive nature of the lithosphere beneath 100 km might reflect the chemical alteration observed in garnet xenocrysts.
In a recent study, Özaydın et al. (2021) compared garnet xenocrysts, xenolith water contents, and MT models in the Kimberley-Jagersfontein region. The 3D model presented in this study agrees well with the previous 3D modeling and probabilistic inversion ( Figure S6 in Supporting Information S1). These analyses demonstrated that water contents calculated from MT models decreased with depth between 100 and 160 km, which roughly agrees with the measured xenolith water contents. However, a larger misfit was observed in the depleted layer, suggesting a local metasomatic control around kimberlite conduits. Since the fertile layer corresponds to a water-saturated zone and calculated water contents can account for observed xenolith water contents, the authors considered that the metasomatism in the fertile layer has to be laterally extensive enough to be sensed by electromagnetic fields. In support of this argument, we found that log-Nd/Dy N and Yb N values correlate with the number of garnets classified as metasomatized in the fertile layer ( Figure 9). Given that increasing metasomatism leads to less sinuous log-Nd/Dy N patterns and more HREE enrichment in garnet xenocrysts (Griffin et al., 2002), we suggest that the metasomatic fluids pervasively percolated through the fertile layer, peaking around 110-120 km depth. As a result, the kimberlites sampled more fertile material at this depth.
A conductor bounds the Kimberley region to the east (Figure 3a). Like the Molopo Farms and the Bushveld complexes, this conductor coincides with a large Paleoproterozoic mafic intrusive complex (1.9 Ga, Maier et al., 2003), the Trompsburg Complex (TC, Figure 1). Since it does not have surface exposure, existing knowledge about the complex is limited to gravitational modeling and few borehole measurements (Rezaie et al., 2017). We think the conductor here may be related to the emplacement of the mafic magmas and a metasomatized residue left from the event.

Bushveld Magmatic Event and Its Effects on Composition of the Lithosphere
The Bushveld Complex (BC) is a large igneous body emplaced between 2.055 and 2.06 Ga, intruding Transvaal Group sediments (Scoates et al., 2021). The complex outcrops as mafic-ultramafic eastern, western, and northern lobes (Rustenberg Layered Suite) and slightly younger felsic rocks toward the center (Figure 10a). A strong stratigraphic and geochemical concordance between the mafic intrusives in the different lobes suggests connectivity between them or a co-genetic relationship (Kruger, 2005). Gravity models of the Bushveld Complex demonstrate that surface exposures of the complex are continuous at depth and define a lopolith structure with a maximum thickness exceeding 8 km (Cole et al., 2014). The Palmietgat kimberlite pipe, near the center of the Bushveld Complex, contains xenoliths of plagioclase-bearing lower-crustal pyroxenites petrologically akin to the pyroxenites of the Rustenburg Layered Suite, indicating connectivity between the underlying western and eastern lobes (Webb et al., 2011). Zircon U-Pb dating constrains the emplacement of the Rustenberg Layered Suite to a one million year time window, indicating extremely high magma fluxes (Zeh et al., 2015). The lower units of the Rustenberg Layered Suite hold economically valuable, PGE-hosting chromitite layers (e.g., Merensky Reef) that contain more than 80% of global platinum reserves (VanTongeren, 2017). Both the source of the parental magma and the magma chamber differentiation processes necessary to produce such PGE-rich chromitite layers are still debated (Cawthorn, 2015;VanTongeren, 2017). Considering the targeting depths of our MT models, the most relevant ideas to investigate are speculations about the mantle component required in the parental magmas to generate the specific properties observed in the intrusives (Richardson & Shirey, 2008;Zirakparvar et al., 2014).
More than 10 adjacent magmatic bodies, each with geochemical affinities to the BC, were emplaced during the same narrow time window (Figure 10a, Rajesh et al., 2013), and generally are considered to be satellite intrusions. One of these satellite magmatic bodies, the Molopo Farms Complex (MFC, Figure 1), is a layered mafic intrusion situated in southern Botswana roughly 400 km west of BC (Beukes et al., 2019). The stratigraphic layering of the MFC is similar to the lower layers of the BC, but available borehole data indicate it lacks the PGE-rich chromitite layers (e.g., Kaavera et al., 2018). The BC and these satellite intrusions may have formed as part of a larger Bushveld Large Igneous Province, which used a structural weakness corridor along the TML (Rajesh et al., 2013).

Garnet Xenocrysts and Composition of the Bushveld Mantle
Garnet xenocrysts that sampled the mantle in the vicinity of the BC show more variable geochemical signatures from pipe to pipe, in contrast to those from Proterozoic terranes or those from the Kimberley region. They are concentrated around the BC or close to the limbs of the Rustenberg Layered Suite. Only the Palmietfontein and Palmietgat pipes fall within the boundaries of the complex. Compared to those from the southwestern pipes, the garnet xenocrysts from the pipes near the BC demonstrate higher whole-rock Al 2 O 3 , lower Mg ol #, a general increase in metasomatic classes and a lack of depleted material (Figures 6c-6e).
The Palmietgat pipe is one of the very few kimberlite pipes emplaced near the center of the BC (Figure 7i). Quite a high proportion of the garnets here (41%) do not fall into statistically significant CARP classifications associated with a consistent metasomatic signature. The high proportion of unclassified garnets at the Palmietgat pipe may indicate complex metasomatic overprinting, with multiple episodes of melts and fluids using the same conduits. A general increase in the proportion of unclassified samples near the BC can be seen in Figure 6f, suggesting the metasomatic signature in the area is highly complex.
The kimberlite cluster immediately west of or penetrating the western limb of the Bushveld Complex (Western Bushveld, Figure 7i) exhibits a thick layer of depleted lherzolites with phlogopite metasomatism between 130 and 160 km. Like other sections near the Bushveld Complex, many unclassified samples come from this middle section (130-170 km). These unclassified samples fall into wehrlite and low-Cr CaO − Cr 2 O 3 classifications at the fertile layer and beneath the BDL, whereas they appear to be lherzolitic in the middle section. The Palmietfontein pipe ( Figure S22c in Supporting Information S1), like the Premier sections ( Figure S23c in Supporting Information S1), completely lacks samples from this middle section where the unclassified garnets are most populous.
Just south of the BC, the Premier kimberlite cluster (Figure 7j) is considerably older (928-1,150 Ma) and yields higher equilibration temperatures (40 mW/m 2 ) compared to the Bushveld West cluster and the Palmietgat pipe (90-153 Ma, 38 mW/m 2 ). The amount of secular cooling expected between 928 Ma and 153 Ma is around (82°C, Shu et al., 2014), a value similar to the difference in geotherms between the two clusters (∼50°C−100°C). Therefore, the area is not likely to have been impacted by a large thermal event between these times. This is consistent with the similarity of the geochemical "tomography" sections showing a large portion of unclassified samples and depleted lherzolites with phlogopite metasomatism in the middle of the section. Like the Palmietfontein pipe of the Bushveld West cluster, the Premier and Franspoort pipes have entrained a relatively low number of garnets from this middle section (Figures S11c, S23c in Supporting Information S1); it is probable that the strong metasomatism has destroyed most pre-existing garnet in the rocks at these levels.
The mantle beneath the MFC, on the other hand, is sampled most closely by the Triassic-aged Jwaneng and Thankane kimberlite pipes (Figure 7d). The equilibration temperatures of the sampled xenocrysts are lower than in the BC (35 vs. 38 mW/m 2 ). The lower geotherm might indicate that the event that increased the geotherm beneath the Bushveld Complex, whether it be the main Bushveld event or not, did not affect the mantle near MFC at a similar scale. Geochemical stratification in these pipes is also more akin to the mantle beneath the Kimberley domain, where an increased population of depleted material is evident at mid-lithosphere depths. REE N patterns display trends of metasomatic intensity at depths similar to those observed in Kimberley. A peak metasomatic intensity (lowest log-Nd N Dy N , highest Yb N ) is observed at 125 km ( Figure S44 in Supporting Information S1). However, base-level rates of log-Nd N Dy N are much lower than the Kimberley (Figure 6g), suggesting metasomatic effects might be more extensive or that the area originally underwent less depletion.

Bushveld Geoelectric Structure and Other Geophysical Studies
In our MT models, the Bushveld region shows marked lateral heterogeneity with some strongly conductive vertical sections extending toward the crust (e.g., BCC, Figure 10). The TML crosses three contrasting mantle electrical domains: a conductive lower lithospheric mantle and resistive upper lithospheric mantle in the vicinity of the MFC, a resistive keel-like area beneath the Kanye Basin and a very conductive lithospheric mantle associated with the BC (Figure 10a). Broadly, such resistivity differences can either have thermal or compositional causes. S wave seismic receiver function studies give no clear indication of a thinned mantle lithosphere beneath this region (Ravenna et al., 2018;Sodoudi et al., 2013). However, the Bushveld region is associated with lower S and P wave velocities (e.g., Fouch et al., 2004;Ortiz et al., 2019;Youssof et al., 2013) which suggests a chemical/ metasomatic origin of the conductive anomalies rather than a thermal one .
Seismic receiver function analyses indicate that the western-central area of the BC is marked by a relatively high crustal thickness (∼49 km) with a more gradual Moho transition compared to sharper waveforms observed in the surrounding crust (Delph & Porter, 2015;Youssof et al., 2013). Kgaswane et al. (2012) also suggest structural complexity at crustal depths to account for observed receiver functions. The depressed Moho can be a sign of underplating by mafic intrusions . This anomalous zone overlies the Bushveld complex conductor (BCC), extending from >200 km to near-surface where the MT transect crosses the TML. This region is a good candidate for a zone of mechanical weakness in the lithosphere where the Bushveld magmas could ascend and precipitate interconnected conductive minerals. The TML has been invoked many times before as the feeder location for the Bushveld magmatic units (e.g., Clarke et al., 2009). Such Moho complexity and localized increased Moho depths are not observed beneath the MFC, which instead overlies gradual decreases in both Moho complexity and crustal thickness (Delph & Porter, 2015;Youssof et al., 2013).
The mantle beneath the MFC has been modeled to be conductive at depths greater than ∼110 km but the resolution of this conductor is not robust and this feature is not as strong in the other inversion models (Moorkamp et al., 2021). At shallower depths, resistivities are comparable to the mantle beneath the Kanye Basin. If we assume similar parental magmas for the BC and MFC, the differences in resistivity of the mantle beneath them may suggest: (a) The magma emplacement that caused the MFC was not as focused as that which caused the BC, or the volume of magma intruding the lithosphere was smaller, (b) the MFC was not emplaced from deep feeder dykes but instead magmas flowed horizontally from the primary Bushveld magma source, most likely along the TML (Prendergast, 2012). Our models can not differentiate between these emplacement models but provide evidence for conductive metasomatic signatures beneath. However, if one expects similar mantle compositional sections from the BC and the MFC and a geotherm indicating lack of thermal impact beneath the MFC, the second option may seem more probable. A denser MT data collection around both areas may help to address this situation.

Compositional Causes of Conductivities
We have tested different compositional scenarios to fit the observed ranges of conductivity of the mantle beneath the Bushveld region. The thermal structure was constrained using the 38 mW/m 2 geotherm derived from garnet xenocryst data for the Cretaceous Group I kimberlites around the western lobe of the Bushveld Complex ( Figure 7k) and the Jwaneng-Thankane pipes (35 mW/m 2 ). The Bushveld geotherm agrees well with other calculations derived from measured heat flow and crustal heat generation data (M. Jones, 2017). Similarly, thermobarometry data calculated with the method of Nimis and Taylor (2000) from Jwaneng samples fits a geotherm of 36.3 mW/m 2 ( Figure S63 in Supporting Information S1, Preston & Sweeney, 2003). These results imply that with the maximum water contents allowed by the water solubility model (Padrón-Navarta & Hermann, 2017) and water-partitioning coefficients (Figure 11a, Demouchy et al., 2017), the water in NAMs in a lherzolitic mantle cannot explain the conductivities observed.
Another possible explanation for such conductive mantle is the existence of well-connected fluorine-bearing phlogopite (Y. Li et al., 2017). Abundant phlogopite style metasomatism is observed in almost all compositional sections near the Bushveld region (Figure 7). However, these signatures do not necessarily indicate the presence of perfectly connected phlogopite grains since they are the byproducts of reactions in which phlogopite replaces sparsely distributed garnets (Van Achterbergh et al., 2001). However, phlogopite could still precipitate on boundaries between NAM grains during percolation of metasomatic fluid or melt. Such a model with conductive phlogopites can also be envisioned as peridotites veined with MARID (Mica-Amphibole-Rutile-Ilmenite-Diopside) and PIC (Phlogopite-Ilmenite-Clinopyroxene) assemblages (Foley, 1992). This model of hydrous mineral emplacement is more likely to produce a mantle with highly interconnected phlogopite. We used the median F (fluorine) contents of glimmerite (MARID) xenoliths (∼0.425 w.t.) as a realistic estimate and the maximum F content of all peridotite mantle xenoliths (∼1.5 w.t.) as a bound on maximum conductivity that can be observed ( Figure S62 in Supporting Information S1). For modal compositions, we derive a lherzolitic matrix with 12% and 6% of perfectly interconnected phlogopite. With these compositions, the conductivity of the mantle beneath Bushveld Complex can be explained below 100-120 km with 6% and 12% phlogopite with 1.5 wt% and 0.425 wt% F, respectively. The same compositions would explain conductivities at relatively greater depths (∼120-150 km) in the colder mantle (35 mW/m 2 ) beneath the MFC.
It is impossible to model the conductor beneath the Bushveld Complex above 110 km (BCC) using experimental electrical conductivity data for the major rock-forming minerals existing in the mantle (olivine, pyroxenes, garnet, phlogopite, and amphibole; Figure 11). Therefore, the conductivity requires minor accessory minerals such as graphite or sulfides, either as films on the edges of grains or as crystalline material precipitated from a metasomatic fluid (Pearson et al., 1994). However, the behavior of these materials in a peridotitic matrix is not well understood, and some existing studies are in apparent conflict with each other (Wang et al., 2013;Zhang & Yoshino, 2017). Graphite, for instance, may establish an effective conductivity with higher than usual amounts of carbon in the mantle (>1.6% vol., Wang et al., 2013). However, other authors suggest that graphite films are not stable with more realistic concentrations of carbon in the sample (≪0.8% vol. Watson et al., 2010;Zhang & Yoshino, 2017) and the previous results of Wang et al. (2013) may have only been the result of temporally limited crystallization, which would not retain its interconnection given enough time (Zhang & Yoshino, 2017). Furthermore, the nature of carbon speciation in the mantle is not well understood, and the processes that bring xenoliths to the surface may result in compositionally biased xenolith samples (Stagno et al., 2019). For instance, some rare xenoliths from the Premier and Jagersfontein kimberlites are very rich in large veins of crystalline graphite (Pearson et al., 1994), suggesting that crystalline graphite might be underrepresented in xenoliths since they are likely to be destroyed by the carrier magma. On the other hand, the Palmietgat kimberlite is diamondiferous (Webb et al., 2011), which suggests the possible existence of graphite above the graphite-diamond transition.
Recent detailed studies on peridotite mixtures with magnetite (Dai et al., 2019) and chromite (W. Sun et al., 2021) showed that these types of minor constituents might have very different percolation thresholds before they dominate electrical conduction in an assemblage, probably reflecting their precipitation behavior in the matrix. The percolation threshold represents the volume needed for a mineral to dominate the electrical conduction in an assemblage. While it is hard to reconcile any of the required volumetric abundances used in these experiments to reflect the actual state of the mantle (16% for chromite, 1.5% for magnetite), they illustrate the possibility that such phases could have a significant effect on metasomatized portions of the mantle.
Nevertheless, all Bushveld aged magmatism in the Kaapvaal Craton is strongly associated with a conductive signature in the MT models, which supports the idea that the generation and emplacement of Bushveld magmas left a metasomatized signature in the lithospheric mantle (e.g., Richardson & Shirey, 2008;Zeh et al., 2015).

Namaqua-Natal Belt and Uintjiesberg
Southwest of the Kimberley area and over the proposed craton boundary, a different geoelectrical structure is apparent. Here, MT models indicate a lithosphere composed of a body (Eastern Namaqua-Belt Resistor, ENBR) that is highly resistive above 100-120 km and much less resistive below that depth (Figure 3a). Due to the thick Karoo sedimentary sequence, it is not entirely clear whether the southern limit of this resistive area is the border of the Kaapvaal Craton or not (Kobussen et al., 2008;Weckmann, 2012). However, the continuation of potential field lineaments alongside the Kheis Belt units suggests that the existing craton boundary coincides with our so-called suture zone (Corner & Durrheim, 2018). The significantly deeper Moho depths south of this boundary (Delph & Porter, 2015;Youssof et al., 2013) and the Proterozoic ages of the lower crustal granulite xenoliths from kimberlite clusters (Schmitz & Bowring, 2004) also suggest that the existing craton boundary runs through this location. We suggest that this resistor represents a Paleoproterozoic microcontinental block at the northernmost front of the Namaqua-Natal Belt where it collided with the Kaapvaal Craton.
Around the eastern Namaqua-Belt Resistor, the Uintjiesberg cluster consists of non-or weakly diamondiferous Group I kimberlites (76-103 Ma, Figure 5a). Garnet xenocrysts indicate a mantle section with higher equilibration temperatures (40 mW/m 2 ) than the nearby Kimberley region (37 mW/m 2 ). The Group I thermal structure (Figure 5b), geoelectric architecture (Figure 3a), and calculated water contents (Figure 8) are consistent with a shallower base of the depleted lithosphere (135 km, Figure S41 in Supporting Information S1) and a higher proportion of melt-metasomatized classes, suggesting the region experienced extensive metasomatism during a thermal event (Kobussen et al., 2008). Differences in temperature and chemistry between Group I and Group II kimberlites in the region also suggest that extensive melt-infiltration and metasomatism occurred between these times, raising the base of the depleted lithosphere by ∼40 km (Kobussen et al., 2008).
From just above the base of the depleted lithosphere, a sudden drop in mean Ti/Zr ratios is evident and continues up to 110 km, suggesting phlogopite-related metasomatism at these depths. Considering that this layer (110-140 km) coincides with a layer of low resistivity in the MT models, it might be rich in interconnected phlogopites. However, electrical conductivity calculated with a perfectly connected phlogopite of 5%-10% with 0.425 wt% F is roughly an order of magnitude more conductive than modeled resistivities ( Figure S7 in Supporting Information S1). This calculation indicates that electrical conductivities observed in the region are likely due to water in NAMs and not due to minor interconnected phlogopite. Fertile garnets in Uintjiesberg cluster also demonstrate the least sinuous patterns with strong HREE enrichment (Figures 6g and 6h). However, most of these garnets are between 100 and 110 km depth, where lower Ti/Zr values are not observed, suggesting that the metasomatic signature might not be as extensive at these depths.

Limpopo Belt
The Limpopo Belt is a continental block wedged between the Zimbabwe and Kaapvaal cratonic assemblages. Seismic tomography models demonstrate the existence of a thick cratonic lithosphere with fast velocities beneath the Limpopo Belt (Fouch et al., 2004;A. Li & Burke, 2006;Ortiz et al., 2019;Youssof et al., 2015). However, these models also show a seemingly fragmented structure, with relatively slower velocities in the eastern Limpopo Belt. These slower velocities may be attributed to metasomatism related to intense and pervasive Jurassic Karoo magmatism. Many authors have envisioned that the Karoo plume head was centered roughly on the eastern Limpopo Belt, based on the distribution of dyke swarms and volcanism radiating from the area (Figure 1, KTJ, Jourdan et al., 2006). Whether or not the Karoo event was initiated by a mantle plume related to the Gondwana breakup, the eastern Limpopo Belt seems to be the center of intense Jurassic volcanism.
Our MT models show that the Eastern Limpopo Belt comprises a highly conductive lower crust (≤10 Ωm) underlain by a highly conductive mantle. Sensitivity tests show the mantle between 50 and 100 km depth can be fit by a block with resistivities ∼60-100 Ωm. Like the Bushveld Complex, the region also exhibits locally increased Moho depths with gradual velocity gradients (Delph & Porter, 2015;Youssof et al., 2013), which may suggest mafic underplating.
Garnet xenocrysts from kimberlites in the Venetia cluster, which was emplaced around the time of Gondwana assembly (∼572-583 Ma), indicate a considerable proportion of depleted harzburgites and lherzolites compared to the Kimberley area (Figure 7h). The high conductivity (Figures 7 and 3a) and lower seismic velocities modeled in the area do not reflect these garnet xenocryst compositions. The discrepancy between xenocrysts and geophysical models point to metasomatism of the lithospheric mantle that post-dated eruption of the Venetia kimberlites, likely to be effects of the ∼180 Ma Karoo event, as previously suggested by Griffin et al. (2003).

Gibeon Fields
The Gibeon region stands out on the MT models with its resistive root and young, non-diamondiferous kimberlite cluster (Figures 7 and 12) around a relatively resistive region. The region is roughly situated within the Rehoboth Terrane along the boundary with the Namaqua-Natal Belt and is marked by a sudden drop in magnetic intensity (Corner & Durrheim, 2018). The kimberlites here are among the youngest in southern Africa (60)(61)(62)(63)(64)(65)(66)(67)(68)(69)(70)(71)(72) and are considered to be related to the same plume responsible for the emplacement of Group I kimberlites in Kimberley (Davies et al., 2001). Garnet xenocrysts show a mantle section mostly made up of fertile lherzolites with relatively constant Mg ol # (Figure 7a). Unlike what is observed in the Archean regions and the Uintjiesberg cluster, almost all of the fertile samples have high-Ti and low-Zr signatures likely to indicate low degrees of phlogopite metasomatism (Figure 6e). The geotherm derived from garnet xenocrysts can only be estimated from low-temperature garnets while high-temperature estimates seem to be scattered, similar to observations in other thermobarometric calculations  and likely to be an effect of melt-metasomatism beneath the depleted lithosphere. Two-pyroxene thermobarometry (Brey & Köhler, 1990) carried out by Bell et al. (2003) fits best to a conductive 41.5 mW/m 2 geotherm if the high-temperature samples between 4 and 5 GPa are excluded. While these high-temperature samples may indicate melt-related metasomatism, there is no indication of large-scale melt infiltration in the geochemical tomography section and stable (Mg ol # = 92) values are observed.
Since the 3D MT model shows a highly resistive region with low estimated water contents, we propose that the mantle here might not be extensively metasomatized. The high resistivities suggest that there is no pervasive phlogopite metasomatism, while the abundant fertile lherzolites may indicate that the mantle never really experienced strong depletion. Due to a lack of data, we do not know whether this is a feature along the entire boundary between the Rehoboth Terrane and Namaqua-Natal Belt or if it is limited to the area around the Gibeon kimberlite field.

Distribution of Kimberlites
Kimberlites are hydrous and carbonated volcanic rocks that originate from low-degree melting of the mantle below ∼150 km (Giuliani et al., 2020). Because they ascend very rapidly (>4-20 m/s, Sparks et al., 2006) and entrain xenoliths and xenocrysts from surrounding wall-rock, they provide an invaluable window to the deeper mantle. The current most accepted mechanism for the faster ascent rates is crack-tip propagation via CO 2 exsolution, in which carbonated melts cleave the mantle by carbon degassing at the top of the percolation front while dissolving and assimilating orthopyroxene from the wall rock (Giuliani et al., 2020;Russell et al., 2012).
Clifford's rule (Clifford, 1966) suggests that diamondiferous kimberlites occur in terranes with Archean crust. Today we know that diamondiferous kimberlites frequently sample the cratonic mantle because diamonds are only stable below 130-150 km at colder cratonic geotherms (Day, 2012). Temporally, the genesis of kimberlites is linked to major tectonic events associated with global plate reorganizations such as the assembly and disruption of Gondwana (Jelsma et al., 2009). Spatially, kimberlites are mostly restricted to the edges of depleted cratonic keels comprised of chemically refertilized lithospheric mantle (Faure et al., 2011;Griffin et al., 2009). This specific spatiotemporal distribution of the kimberlites is correlated with fossil and active tectonic boundaries (Jelsma et al., 2004). These are commonly observed as lineaments in potential field geophysical studies (Corner & Durrheim, 2018), possibly indicating the tendency of kimberlites to use pre-existing zones of weakness in the lithosphere.
Through interpolated maximum resistivity maps of the lithospheric mantle, A. G. Jones et al. (2009) showed that kimberlites in southern Africa tend to concentrate near the edges of deep resistors that may correspond to depleted and dry cratonic keels. To compare the distribution of kimberlites and our magnetotelluric models, we made maps of mantle conductance ( ∑ = ; i.e., integrated conductivity (σ) times thickness (d)) between 50 and 150 km (Figure 12a). Another map (Figure 12b) shows where the greatest depth to the 1,000 Ωm value is observed. While containing similar information to the conductance map, this map emphasizes the spatial distribution of thick resistive roots. These maps demonstrate that kimberlites tend to be emplaced at the edges of the deepest and most resistive parts of the cratonic lithospheric mantle while also avoiding the most conductive parts.
One reason for the distribution of kimberlites around the deep resistors might be the inability of kimberlitic melts to penetrate thick and depleted cratonic roots. In such depleted regions, the kimberlite magmas may react with the refractory wall-rock materials to the extent that they cannot ascend to the surface (Giuliani et al., 2016). If these highly resistive keels have an orthopyroxene-poor (<15%) dunitic composition as suggested by some authors (e.g., Griffin et al., 2009), the rapid-ascent mechanism of crack-tip propagation via CO 2 degassing through assimilation of orthopyroxene from the wall rock might not operate efficiently (Russell et al., 2012). Generation of kimberlitic melt requires either the lithospheric mantle or the magma source beneath the lithosphere to be enriched in volatiles. This could be possible via metasomatism and carbonate freezing near the base of the lithosphere, including by Phanerozoic subduction processes (C. Sun & Dasgupta, 2020). A recent study suggests that the enriched signatures observed in kimberlite-hosted xenoliths globally reflect assimilated lithospheric mantle material (Giuliani et al., 2020). These factors also make penetration of kimberlites through deep resistors less likely, because higher water and carbon contents would be required to melt mantle peridotite in keels that extend to greater depths (Foley & Pintér, 2018) and thick Archean cratonic keels are unlikely to have been affected by subduction-related metasomatism. All of these things considered, it is more likely for kimberlites to episodically erupt at places that have previously been enriched in orthopyroxene, metasomatized and carbonated, and therefore have a moderate to low resistivity. Such events are more likely to concentrate around fossil continent/ micro-continent collision zones with metasomatized compositions, as supported by correlations of kimberlite occurrences with magnetic and geological lineaments (Jelsma et al., 2004).
In the model, there are regions at the edges of the resistive roots (>1,000 Ωm) where the lithospheric mantle has lower resistivities (100-1,000 Ωm) at depths greater than 100 km. There are clear associations between these areas and kimberlite clusters ( Figure 12) and almost all diamondiferous kimberlites occur within these relatively conductive zones. These conductive zones are indicated in Figure 12b as lighter blue areas usually appearing near darker blue, thick resistive zones. In contrast, non-diamondiferous kimberlites can be observed in these lower resistivity areas as well as in the more resistive portions of the cratons. These lower-resistivity regions are likely to have gone through broad-scale metasomatism, as also shown by their garnet xenocryst geochemical "tomography" (Figure 7) and the compositions of whole-rock xenoliths, which include modal metasomatic minerals (e.g., Grégoire et al., 2003). Accordingly, a general recipe for diamond exploration with magnetotelluric data would be to search for regions located between the resistors of the cratonic nuclei and the conductors of metasomatized trans-lithospheric weakness zones.
On the other hand, it is not clear why kimberlite clusters avoid the prominent mantle conductors around the resistive cratonic nuclei. Similar behavior has also been observed seismically, whereby kimberlites avoid regions of low shear wave velocity . To examine this question further, detailed descriptions of these conductors have to be considered. One example of a prominent mantle conductor is the one at the southwest of the Kimberley area (SKC, Figures 3a and 3b). This conductor extends sub-vertically toward the surface at the edge of the craton where the Kheis Orogenic Belt and related fault systems are situated. This feature may be related to the collision between the proto-Namaqua-Natal Belt and the Kaapvaal Craton, and such conductors are often observed along fossil and active suture zones (e.g., Kelbert et al., 2019;Kirkby et al., 2020). Using the structural weakness corridor, Group I kimberlites may have further mineralized this pathway with rising metasomatic fluids. If kimberlites preferentially ascend along pre-existing zones of lithospheric weakness (e.g., Jelsma et al., 2004), this would be a perfect candidate. However, the area around this conductor contains only a few kimberlites, while most of the kimberlites in this region instead overlie relatively resistive mantle.
The absence of kimberlites associated with the Molopo Farms and BCCs may provide insight into the reasons kimberlites appear to avoid conductors more generally. Since both mantle conductors are associated with large mafic intrusions, one could argue that low-degree kimberlitic melts might have a rheological preference to move laterally along the contact with the Transvaal Basin when they meet the mafic material. This model could explain the near absence of kimberlites erupted over these conductors. Another possibility might be that a very depleted mantle exists near the base of the lithosphere (>200 km) due to melt pooling in large volumes at this depth until it finally penetrates the lithospheric mantle. This model of Silver et al. (2006) is suggested to explain the very short duration (<1 m.y., Zeh et al., 2015) of emplacement of the Bushveld Complex. In such a scenario, kimberlites might not be able to penetrate this depleted layer or be generated within in due to low amounts of volatiles left in the layer after this massive melt extraction (Figure 13), similar to emplacement behavior observed at resistive cratonic nuclei.

Concluding Remarks and Summary
3D MT models of southern Africa and geochemical information from garnet xenocrysts in Group I kimberlites are compared in the light of their tectonic and magmatic history. Both qualitative and quantitative interpretations were made for MT-derived composition analysis. Some of the most critical discussion points are: 1. No general relationships can be made between the age of the terranes and conductivity distribution in the mantle since magmatic processes can either deplete or refertilize the lithospheric mantle 2. Most of the Archean cratonic regions are imaged as relatively depleted regions with lower whole-rock Al 2 O 3 values, higher mean Mg ol #, and a lower proportion of metasomatized CARP classes between 100 and 150 km. Some Archean areas around the Bushveld and MFC es appear to have a more fertile signature and metasomatized characteristics, matching well with electrically conductive mantle and the large-scale magmatic event that occurred at 2.05 Ga 3. The area around the Bushveld Complex contains an increased proportion of unclassified garnets in the CARP classification scheme, suggesting complex overprinting of multiple magmatic episodes. Conductivities below 100 km cannot be explained by water but can be explained by perfectly connected phlogopites. Conductivities above 100 km suggest that well-connected accessory minerals (e.g., sulfides, graphite, or chromite) are likely to be present, precipitated extensively using the lithospheric weakness zone along the TML 4. The electrical conductivity of the cratonic mantle reflects the style of metasomatism, that is, the composition of metasomatic fluids, how they were emplaced and the extent of the metasomatism, rather than reflecting simply the fertility of the mantle peridotite. This is best exemplified by the lithospheric mantle in the vicinity of Gibeon Fields, which is geochemically fertile yet highly resistive 5. Kimberlites in southern Africa are more likely to be observed around the edges of highly resistive mantle regions but also avoid the most conductive areas. They are most populous where the mantle below 100 km depth is relatively conductive. In resistive, depleted, orthopyroxene-poor, and volatile-poor mantle, kimberlites are likely unable to penetrate through the mantle or cannot be generated 6. Conductive regions with no co-located kimberlites are usually associated with a mafic intrusion (e.g., Bushveld Complex, MFC). We suggest this is due to: (a) A rheological preference of kimberlitic melts to intrude the contact with the sedimentary units when they meet the mafic layer; and/or (b) during an event like Bushveld magmatism, in which significant volumes of magma are rapidly emplaced, volatiles near the lithosphere-asthenosphere boundary are also intensively extracted, created a zone where kimberlites either cannot penetrate or be generated Figure 13. Conceptual sketch of our interpretation of mantle beneath the Bushveld Region. Scenarios 1 and 2 represent the possible reasons why kimberlites might avoid the center of the conductors. In Scenario 1, kimberlites avoid the conductors because they cannot penetrate or be generated within the depleted mantle formed by intense melt extraction during the Bushveld event. Scenario two posits the possibility that the kimberlites might be deflected by the mafic layer.