Two billion years of episodic and simultaneous websteritic and eclogitic diamond formation beneath the Orapa kimberlite cluster, Botswana

The Sm–Nd isotope systematics and geochemistry of eclogitic, websteritic and peridotitic garnet and clinopyroxene inclusions together with characteristics of their corresponding diamond hosts are presented for the Letlhakane mine, Botswana. These data are supplemented with new inclusion data from the nearby (20–30 km) Orapa and Damtshaa mines to evaluate the nature and scale of diamond-forming processes beneath the NW part of the Kalahari Craton and to provide insight into the evolution of the deep carbon cycle. The Sm–Nd isotope compositions of the diamond inclusions indicate five well-defined, discrete eclogitic and websteritic diamond-forming events in the Orapa kimberlite cluster at 220 ± 80 Ma, 746 ± 100 Ma, 1110 ± 64 Ma, 1698 ± 280 Ma and 2341 ± 21 Ma. In addition, two poorly constrained events suggest ancient eclogitic (> 2700 Ma) and recent eclogitic and websteritic diamond formation (< 140 Ma). Together with sub-calcic garnets from two harzburgitic diamonds that have Archaean Nd mantle model ages (TCHUR) between 2.86 and 3.38 Ga, the diamonds studied here span almost the entire temporal evolution of the SCLM of the Kalahari Craton. The new data demonstrate, for the first time, that diamond formation occurs simultaneously and episodically in different parageneses, reflecting metasomatism of the compositionally heterogeneous SCLM beneath the area (~ 200 km2). Diamond formation can be directly related to major tectono-magmatic events that impacted the Kalahari Craton such as crustal accretion, continental breakup and large igneous provinces. Compositions of dated inclusions, in combination with marked variations in the carbon and nitrogen isotope compositions of the host diamonds, record mixing arrays between a minimum of three components (A: peridotitic mantle; B: eclogites dominated by mafic material; C: eclogites that include recycled sedimentary material). Diamond formation appears dominated by local fluid–rock interactions involving different protoliths in the SCLM. Redistribution of carbon during fluid–rock interactions generally masks any potential temporal changes of the deep carbon cycle.


Introduction
Diamonds and incorporated inclusions allow researchers to investigate large-scale processes such as plate tectonics and crustal recycling (Farquhar et al. 2002;Gurney et al. 2010;Schulze et al. 2003;Shirey and Richardson 2011). They provide samples of Earth's otherwise inaccessible sub-continental lithospheric mantle (SCLM) and a record of deep cycling of volatile elements (Cartigny et al. 1999;Deines 1980;Howell et al. 2020;Shirey et al. 2019). Together, these observations have led to the current paradigm that diamond formation involves metasomatic processes in the mantle comprising volatile-rich, supercritical high-density fluids, possibly triggered by tectono-magmatic events (Deines Communicated by Dante Canil. 1 3 54 Page 2 of 25 1980; Schrauder et al. 1996;Sobolev et al. 1998;Stachel et al. 2005; Thomassot et al. 2007). Insights from diamond inclusion ages, spectroscopic observations of internal growth patterns and from the heterogeneous stable isotope compositions of diamonds worldwide have established evidence for episodic diamond growth in discrete peridotitic (P-type), eclogitic (E-type) or websteritic (W-type) parageneses (Aulbach et al. 2018;Cartigny et al. 2014;Koornneef et al. 2017;Pearson et al. 1998;Richardson et al. 2004;Smit et al. 2016;Timmerman et al. 2017;Wiggers de Vries et al. 2013).
P-type diamonds formed in highly melt-depleted, harzburgitic, residues and lherzolitic substrates Richardson et al. 1984;Stachel et al. 2004a). Based on mineralogy and stable isotope record, several hypotheses have been proposed for the formation of E-type diamonds. In particular, they have been ascribed to subduction of oceanic lithosphere and its incorporation into the peridotitic SCLM (Deines et al. 1984;Li et al. 2019;Taylor et al. 1990), high-temperature volatile fractionation (Cartigny et al. 1998;Galimov 1991;Javoy et al. 1986) or primordial mantle heterogeneities (Deines et al. 1997). Websteritic diamonds have a transitional character with inclusion compositions similar to P-type diamonds but carbon isotopic compositions resembling E-type diamonds (Deines et al. 1997). Their formation appears to be related to re-fertilisation reactions (Bodinier et al. 2008;Smit et al. 2014;Taylor et al. 2003) through mixing of mafic, possibly slab-derived, fluids with peridotitic SCLM (Aulbach et al. 2002) or reaction with eclogites in a halo around an ascending mantle plume (Sleep 2006;Viljoen et al. 2018).
Unprecedented access to inclusion-bearing diamonds from the 93.1 Ma Orapa kimberlite cluster in north-east Botswana (Allsopp et al. 1989;Davis 1977) allows us to further examine diamond formation and the extent and timescale over which these events occurred, to place better constraints on the geological processes that lead to their formation. Recently, Timmerman et al. (2017) presented Sm-Nd isochron ages for eclogitic silicate inclusions from Letlhakane (245 ± 38 Ma; 998 ± 140 Ma; 2334 ± 22 Ma) and Orapa (140 ± 93 Ma; 1096 ± 230 Ma; 1699 ± 340 Ma) mines that confirmed multi-stage diamond growth. Their work supported observations of Deines et al. (1993) for mixing of distinct sources comprising organic carbon (δ 13 C < − 35‰), carbonate sediments (δ 13 C from − 10‰ to − 5‰) and depleted mantle material (δ 13 C ~ − 6‰) to account for the observed compositional variability in E-type diamonds.
In combination with the data from Timmerman et al. (2017), the present study examines diamonds of E-, W-and P-type paragenesis from Letlhakane (LK; n = 78 inclusions), Orapa (OR; n = 2) and Damtshaa (DM; n = 2) mines. Overall, the Orapa cluster is known for its unusual predominance of eclogitic (~ 85%) over peridotitic diamonds (Deines et al. 1993;Gurney et al. 1984) and it is one of the few locations on the Kalahari Craton where a websteritic diamond suite (5%) is present (see Viljoen et al. (2018) and references therein). Thus, in addition to new eclogitic garnet and clinopyroxene inclusion ages from the three mines, we provide the first websteritic and harzburgitic Sm-Nd ages for Letlhakane along with their major and trace element compositions. These data are set in context of the host diamond growth structure, the nitrogen content and aggregation data, as well as carbon and nitrogen isotope compositions of corresponding diamond growth zones. This approach allows assessment of whether diamond-forming processes occur simultaneously or at different times in different parageneses beneath individual mines or entire regions. These constraints are crucial for understanding the nature and style of diamond formation and the broader interpretation of the evolution of the SCLM. Details on regional geology and the major tectono-magmatic events that potentially influenced diamond formation are discussed in Online Resource 1.
After morphological characterisation central diamond plates were cut and polished in {110} orientation from 87 LK diamonds containing multiple large (> 50 μm) garnet, clinopyroxene and sulphide inclusions and these plates were photographed and imaged by cathodoluminescence (CL). Identification of kyanite and coesite was confirmed by Raman spectroscopy. Nine of the plates (LK09, LK10, LK33, LK40, LK42, LK47, LK50, LK75, LK79) were described in Timmerman et al. (2017). The inclusion-bearing diamonds from Orapa (n = 2) and Damtshaa (n = 2), for which inclusion ages and general characteristics are presented here, were selected by the same criteria from a selection of ~ 160,000 diamonds (Gress et al. 2021).

Nitrogen content and aggregation
Infrared absorbance spectra were collected on whole diamonds during sampling at DTCB and along traverses on central plates on several FTIR instruments with different experimental settings (for more details see Gress et al. (2018) and Online Resource 1). Errors for N content and N aggregation (expressed as %B = 100B/[A + B]) are typically ± 10% and ± 5%, respectively, but uncertainties increase with plate Page 3 of 25 54 thickness and decreasing N abundance. The significance of the N aggregation data with low N abundance is limited because N aggregation effectively only takes place with contents > 200 at. ppm (Taylor et al. 1990). Therefore, the data were filtered for spectra with < 200 at. ppm N, and N aggregation data from these samples were excluded from the figures (see Online Resource 1). The terminology used to describe diamond growth structures distinguishes major growth zones where the N abundance showed a variation > 100 at. ppm (i.e., core, intermediate, rim).

Carbon and nitrogen isotopes
Carbon isotope ratios were determined from selected diamond growth zones using a large-geometry secondary ion mass spectrometer (LG-SIMS) at CRPG, Nancy, and gas source mass spectrometry (GSMS) at VU Amsterdam. For GSMS, the major growth zones were sampled from diamond fragments of about 0.05-0.1 mg from known location within a diamond (n = 193 growth zones) and combusted in an element analyser, as described in Timmerman et al. (2017). These analyses are referred to as 'combustion' analyses and are expressed in delta notation relative to Vienna PeeDee Belemnite (VPDB), δ 13 C VPDB [‰] = ( 13 C/ 12 C sample / 13 C/ 12 C VPDB -1) × 10 3 . Reported uncertainties (Online Resource 1) were obtained from the reference materials analysed in the corresponding analytical session and yield ± 0.12‰ to ± 0.20‰ (SD). Overall, duplicates of 14 diamond growth zones were analysed in different sessions with indistinguishable data (Online Resource 2: Table ESM_2.3). For LG-SIMS, a Cameca IMS-1280-h multi-collection ion microprobe equipped with high sensitivity Faradays cups (see Online Resource 1 and Bouden (2021)) measured carbon isotopes on central plates (n = 9 diamonds) in proximity to FTIR analyses and on polished fragments (n = 7 diamonds) obtained from major growth zones. Instrumental mass fractionation (IMF) was monitored using reference material Dp-418 (δ 13 C = − 5.32‰), Nam-56 (δ 13 C = − 29.3‰), Nam-114 (δ 13 C = − 26.93‰), see Cartigny et al. (2004), and give an absolute standard error of 0.14‰ for the entire analytical session. Long-term error propagation (i.e. for the entire session) gives uncertainties ranging from ± 0.15‰ to ± 0.17‰ (SD).
In-situ nitrogen content and isotopic composition were measured on polished fragments (n = 5 diamonds) adjacent to C isotope analysis locations, on the same instrument with similar settings for the primary 133 Cs + beam (Online Resource 1). Nitrogen isotope data are expressed in delta notation relative to air, δ 15 N AIR [‰] = ( 15 N/ 14 N sample / 15 N / 14 N AIR − 1) × 10 3 . Uncertainties associated with counting statistics mainly depend on N content (decreasing uncertainties with increasing concentration), with a typical 2 SD uncertainty of 1‰ in diamonds containing 500 at. ppm.
Instrumental mass fractionation was monitored after every 10 unknowns using reference material Nam-22 and Nam-114 (δ 15 N = 1.4‰ and 4‰, respectively, see Cartigny et al. (2004)) and the long-term standard deviation by IMF is better than 0.32‰. Uncertainties for individual measurements range from ± 0.46‰ to ± 1.73‰ (SD). More details for the methods can be found in Online Resource 1 and Gress et al. (2020).

Inclusion geochemistry
Silicate inclusions were analysed for their major element composition using Jeol JXA-8800 M and Jeol JXA-8350F Electron probe microanalyzers following the method of Timmerman et al. (2015). Analyses on both instruments were conducted with 15 kV accelerating voltage, 25 nA beam current, and 1 µm beam size. After EPMA analyses (Online Resource 2: Table ESM_2.4), the inclusions were washed in ethanol, deionised and Milli-Q water to remove the carbon coating prior to spiking and digesting the inclusions for trace element and isotopic analysis. The detailed procedures are described in Online Resource 1; all corresponding data are provided in Online Resource 2: Tables ESM_2.5 and ESM_2.6.

Results
On average, 3.4% of the 77,900 diamonds at Letlhakane are inclusion-bearing (Table 1). Sulphide inclusions are the most abundant (~ 81%), followed by eclogitic silicate (6.6%), unidentified silicate (5.1%) and peridotitic silicate (3.1%) inclusions. Diamonds with combined eclogitic silicate and sulphide inclusions account for 4.2% of the inclusion-bearing diamonds. Overall, the eclogitic inclusions dominate the silicate population, similar to reports in Shirey et al. (2002). In contrast, Deines and Harris (2004) reported an estimated distribution of diamond inclusions from Letlhakane production of 29.4% eclogitic, 31.7% peridotitic, and 38.9% sulphide inclusions, but it remains unclear according to which sampling parameters the inclusion abundance counting was performed.

Diamond characteristics
All diamonds fall within the resorption spectrum from octahedra to dodecahedra (Online Resource 2: Table ESM_2.1). Four octahedra are twinned and 29 diamonds were broken or chipped. The majority (n = 64) were colourless, 17 were colourless to pale yellow and 10 showed uneven colouration with yellow to brown graining or irregular brown cores with colourless rims.

Growth structure
The complexity of the diamond growth histories is illustrated by the CL images (Fig. 1). Blue CL fluorescence colours are visible in 67% of the diamonds and are related to both the N3 centre (V 1 N 3 ), as well as B and A emission caused by nitrogen impurities in diamond (Collins 1992;van Wyk 1982;Welbourn et al. 1996). In addition to blue fluorescence, 22% of the diamonds have growth zones with yellow-green CL (Fig. 1a) colours either related to singlesubstitutional nitrogen (N 0 ; Boyd et al. (1994)), H3 (NVN; Clark and Davey (1984)) or NiN centres (V 2 NiN + and V 2 NiN; Dischler (2012)). Another 28% of the diamonds have slip lines with green fluorescence (Fig. 1b,c) formed during plastic deformation in the mantle (Varma 1970). The dislocations seem to be associated with the non-luminescent Type IIa layers, which are more vulnerable to plastic deformation (Brookes et al. 2000;Yu et al. 2012).
The growth patterns range from a single well-defined growth zone (Fig. 1d) to samples with multiple (2 to 4+) growth zones (e.g., Fig. 1b, e) that are recognised by sharp boundaries that may border genetically unrelated growth layers. In contrast, resorption boundaries (e.g., Fig. 1e) separating distinct growth events may be more difficult to distinguish, but can be recognised by truncation of growth layering often followed by limited non-luminescent (Type II) growth (15-120 µm). Approximately, 75% of the diamonds show at least one ( Fig. 1f) distinct resorption event and hence record a complex growth history (Fig. 1g). The remaining 25% of samples had either no internal resorption boundary or no clear growth pattern.

Nitrogen systematics
The N content of 246 representative whole stones ranges from 12 to 1153 at. ppm with 9 to 100%B aggregation (Online Resource 2: Table ESM_2.2). Less than 1% of these diamonds are pure Type IaA, 56% are Type IaAB (i.e. 10-90%B), 6% are Type IaB, while the remaining 37% have < 200 at. ppm N and hence, no efforts were made to quantify corresponding time/temperature constraints based on the analysed N aggregation state. The whole stone data are used for reference in plots of [N] vs %B. FTIR traverses were determined on 62 representative Letlhakane plates and plotted together with their C-(N) isotope composition in Online Resource 3, with data provided in Online Resource 2. Of these plates, ~ 20% show a constant C-N distribution close to analytical error (± 10% N and ± 2%B; ± 1‰ δ 13 C), whereas the majority have different N contents, aggregation and C isotope ratios in individual growth zones. Spot analyses (n = 1034) from traverses along the plates of all parageneses range from N below detection limit to 1944 at. ppm N and 0 to 100%B aggregation (Online Resource 2). Of these, data from the five known W-type diamonds are restricted to N contents between 6 and 286 at. ppm with a mean of 86 ± 64 at. ppm (n = 65 spot analyses). The largest N intra-diamond variation (1934 at. ppm) is recorded in eclogitic plate LK03, which ranges from below detection limit to 1934 at. ppm. The largest change in aggregation state is recorded in LK311, 95%B. About 2% of the spot analyses traversing the plates are pure Type IaA, 48% are Type IaAB, 9% are Type IaB, 29% of the measurements have < 200 at. ppm N and an additional 9% have N below detection limit.
Mantle residence temperatures (T MR ), calculated after Leahy and Taylor (1997), were modelled for all plates with sufficient N (i.e. > 200 at. ppm, n = 380 spot analyses) based on diamond formation ages presented below corrected for the 93 Ma kimberlite eruption age. The calculated T MR range from 1088 °C to 1373 °C assuming young diamond formation at kimberlite eruption, to 1021-1271 °C assuming . a LK244 consists of two intergrown diamonds of which one has solid CO 2 in its brown core; b LK241 shows oscillatory growth with a gradual shift in its composition; c LK250 is a low N-bearing diamond; d LK61has high N in combination with very low aggregation and implies young diamond formation; e the successive growth layers in LK272 show a similar trend to LK61, however slight changes in fluid composition between core, intermediate zone and rim are visible; f LK231 has a clear distinction between core and rim in its C-N data, the core that formed at 0.75 Ga was partially resorbed at a later stage potentially in combination with precipitation of the rim growth zone; g LK338 shows very complex growth patterns and different episodes of plastic deformation, dated inclusions were recovered from the core (~ 3 Ga), intermediate zone (~ 1 Ga) and rim (~ 0.3 Ga) zones, the light δ 13 C indicates the involvement of organic carbon in the diamond-forming fluid; h the core of LK47 formed at ~ 2.3 Ga whereas the intermediate and rim zones indicate formation at a later stage diamond formation at 2.3 Ga. Temperature differences between individual growth zones within a single diamond are generally < 50 °C i.e., within the uncertainty (Kohn et al. 2016) but can exceed 150 °C. For example, assuming diamond formation at 1 Ga results in T MR of 1044 °C in LK311 core and 1208 °C in its intermediate growth zone. These results illustrate the limits of one-stage thermal evolution models when applied to different growth zones within a diamond. Unfortunately, there are few diamonds with multiple dated inclusions from different growth zones that have sufficient N (> 200 at. ppm) to model the time-temperature characteristics within a diamond (LK75, LK338; Fig. 1g).
Diamond growth: constant C-N A minority (~ 20%) of the studied diamond plates consist of a single or multiple growth zones that display oscillatory growth with a constant N distribution (± 10%) and minor internal variation in δ 13 C (± 1‰) e.g., LK48 with N < 20 at. ppm and δ 13 C of − 22.3 ± 0.4‰; LK61 with 705 ± 30 at. ppm N and 1 ± 1%B ( Diamond growth: variable C-N The majority of the studied diamonds (~ 80%) display internal variations in N content, N aggregation, δ 13 C (and δ 15 N) between growth zones. Several diamonds record abrupt changes between major growth zones (i.e., Δδ 13 C > 4‰) associated with differences in CL-characteristics and FTIR data (Fig. 1f) e.g., LK01 with δ 13 C from − 14.3‰ (core) to − 5.6‰ (rim). The largest intra-diamond variation of 14.4‰ in δ 13 C is recorded in LK305 between core (− 22.0‰) and rim (− 7.5‰), while LK60 and LK231 have differences of up to 12.8‰ and 12.2‰, respectively (Fig. 1f). The variation within individual major growth zones is normally less than 2‰ e.g., in LK18 (T MR : 1088-1134 °C) with δ 13 C from − 7.2 to − 5.1‰, in LK79 (T MR : 1102-1128 °C) with δ 13 C from − 19.1 to − 17.8‰ and an inclusion age > 2.3 Ga for the intermediate zone (Timmerman et al. 2017), in LK95 (T MR : 1144-1156 °C) with δ 13 C from − 4.9‰ to − 4.5‰. Internal differences of 4.2‰ were found in the core of LK362, of 2.7‰ and 4.4‰ in the intermediate zones of LK50 and LK60, respectively, and of 2.6‰ in the rim of LK231. These data potentially illustrate the difficulties in three-dimensional sampling of diamond fragments of a specific growth zone for GSMS analyses of C isotopes. A Rayleigh fractionation process (Cartigny et al. 2001;Deines 1980) could explain the variations observed within some major growth zones, but few diamonds (e.g., LK241, LK272; Fig. 1b, e) have detailed traverses to fully investigate this possibility.

Trace element patterns
The harzburgitic, eclogitic and websteritic garnets have C-1 chondrite-normalised trace element patterns (Online Resource 2: Table ESM_2.5) generally within the range previously reported in garnet inclusions in diamonds worldwide (Stachel et al. 2004a). The harzburgitic garnets (Fig. 4a) have sinusoidal REE N patterns with LREE N increasing from La to a maximum at Nd, decreasing MREE N and flat to increasing HREE N . In contrast, some eclogitic garnets have a steep positive slope in the LREE N and fairly flat MREE N and HREE N , while others have humped REE N patterns with a positive slope in the LREE N peaking at Nd and decreasing in the MREE N to flat or decreasing HREE N (Fig. 4b). Some garnets have Eu anomalies (Eu/Eu*: 0.25-1.7). Websteritic garnets (Fig. 4c) have flatter REE N patterns than eclogitic garnets but also show a positive slope in LREE N with relatively flat MREE N and HREE N and negative Eu anomalies (Eu/Eu*: 0.34-0.88).
The dated eclogitic and websteritic clinopyroxenes typically show trace element patterns (Online Resource 2: Table  ESM_2.5) with a positive slope in the LREE peaking at Nd and decreasing in the MREE-HREE (Fig. 4d). The slopes from the peak at Nd to Lu are steeper in the eclogitic (13-15 times chondrite) and websteritic (27-77 times chondrite) clinopyroxenes in this study than worldwide sources (Nd/ Lu of ± 5) reported in Stachel et al. (2004a). This may be caused by the different analytical techniques, as this study did not analyse bulk inclusions by laser ablation, but matrix residua corrected to reference materials for > 90% loss of Nd and Sm during column chemistry. Additionally, the overall abundance of REE in cpx inclusions is lower compared to garnets and in cases of small inclusion size (< 5 µg), a strong decrease from MREE to HREE in chondrite-normalised patterns is seen in the garnets resulting in elevated (Ce/Yb) N ratios. Overall, inclusions recovered from the same growth zone of a diamond (e.g., LK338 cpx C, D, E) have comparable REE N patterns and Eu anomalies (Fig. 4). Validity of the trace element compositions was demonstrated in Gress et al. (2021) that confirmed a general agreement between such matrix residua collected from chromatography columns and bulk analysis of inclusions.
Parental fluid compositions Fluid compositions in equilibrium with the silicate inclusions were calculated from partition coefficients for supercritical liquids at 1000 °C and 4 GPa (Kessel et al. 2005). The metasomatic nature of diamond formation involves volatile-rich liquids which are commonly referred as high-density C-H-O fluids or low degree melts at supra-solidus conditions (see Shirey et al. (2019)). Henceforth, the term fluid is used in context with diamondforming processes, because distinction between the role of high-density fluids and melts is difficult to resolve. Overall, the calculated fluids (Fig. 5) show a wide range of LREE enrichment (La N = 12-7013, (La/Gd) N = 2-993 and (Ce/ Yb) N = 20-19174) comparable to the range of compositions recorded in kimberlites (Mitchell 1986) and high-density fluids trapped in fibrous diamonds from other locations in Southern Africa (Schrauder et al. 1996;Weiss et al. 2013). Fig. 3 (right) a Peridotitic, eclogitic and websteritic garnet inclusions in diamonds from the Orapa kimberlite cluster classified by molar Mg# versus Cr 2 O 3 contents with compositional fields from Deines et al. (1993). b Classification scheme of Grütter et al. (2004) with distinction between individual compositional garnet groups i to vi. c Composition of clinopyroxene inclusions in Letlhakane diamonds by molar Mg# versus Cr# for individual compositional groups according to garnets. Literature data (empty markers) from Deines and Harris (2004) Kessel et al. (2005). The fluids show a wide compositional range comparable to high-density fluids trapped in fibrous diamonds from other southern African mines (Weiss et al. 2013); worldwide average kimberlite composition from Schrauder et al. (1996) The inferred fluid compositions for individually dated inclusions recovered from the same growth zone are reproducible with (La/Gd) N for example in E-type DM012 (8-11), P-type LK101 (295-313) and W-type LK113 (13-25).

Sm-Nd isotope age arrays
The broad compositional range of garnet and cpx inclusions define five well-(Group 1-5) and three less well-constrained (Group a-c) arrays in a Sm-Nd isochron diagram (Fig. 6). Published data from Letlhakane are included and combined with the new data from the Orapa and Damtshaa mines. From old to young we distinguish: (Group a) three sub-calcic garnets with T CHUR model ages ≥ 2.9 Ga (2.87-3.38 Ga); (Group b) although with a significant blank contribution (~ 30%), one eclogitic garnet has a T CHUR model age > 2.7 Ga; (Group 1) five cpx and one garnet define a 2341 ± 21 Ma age array with an εNd i of + 5.9; (Group 2) four cpx describe a 1698 ± 280 Ma age array with εNd i of + 3.3; (Group 3) eleven garnets define a 1100 ± 64 Ma age with negative εNd i of − 1.7; (Group 4) four garnets plot along a 746 ± 100 Ma array with εNd i of − 2.7; and (Group 5) seven garnets and three cpx define a line with a slope equivalent to an age of 220 ± 80 Ma with εNd i of 1.3; (Group c) eight garnets with unradiogenic 143 Nd/ 144 Nd (εNd i from − 42.5 to − 2.2) plot below a 140 Ma reference line and are outside the uncertainty of Group 5. Notably, inclusions on the age arrays (Group 1-5) are of both eclogitic and websteritic paragenesis suggesting diamond formation occurs contemporaneously in both eclogitic and websteritic protoliths.

Discussion
The recognition of multiple growth zones and traces of plastic deformation within the diamonds (Fig. 1a, c, f, g) is in line with their isotopic heterogeneity and the marked compositional range of the inclusions. Together, diamonds and inclusions provide incontrovertible evidence for diamond formation in a variety of distinct eclogitic/websteritic/harzburgitic/lherzolitic environments in the SCLM beneath the Orapa kimberlite cluster at different times. Similar observations were made elsewhere at various mines [see compilation in Stachel and Harris (2008)].

Sm-Nd isotope constraints on diamond formation ages
Garnet within a harzburgitic assemblage contains the vast majority of the REE (> 90%) and hence, the model ages of sub-calcic garnets have direct age significance. Initially, we will review the peridotitic inclusions in an attempt to constrain the formation age of the oldest diamond populations. Inclusions from eclogitic assemblages, the majority of inclusions in this study and in Timmerman et al. (2017), will be considered thereafter, before examining inclusions of websteritic paragenesis from Letlhakane and subsequently eclogitic inclusions from Orapa and Damtshaa.
Defining which inclusions are co-genetic is a key aspect of the approach. Most growth zones surrounding analysed inclusions have low nitrogen and thus, the nitrogen systematics cannot be used as a main discriminator for determining which groups of minerals are co-genetic and hence potentially form isochron relationships. Since depleted mantle model ages (T DM ) yield unrealistically old ages due to intersection of the depleted mantle model curve with CHUR at ~ 2.9 Ga, Nd isotope T CHUR model ages are reported. Overall, the model ages are only considered indicative for E-/W-type diamonds due to the evidence of recycled crustal material in many studied eclogitic diamond populations (Deines et al. 1984;Smart et al. 2011;Taylor et al. 1990) and due to variable trace element partitioning into clinopyroxene and garnet. Such partitioning depends on the presence of accessory phases and, for example, the CaO contents of garnets fractionating the Sm/Nd ratios in the mineralfluid system (Aulbach and Jacob 2016; Barth et al. 2001;Harte and Kirkley 1997;O'Reilly and Griffin 1995). This complicates interpretation of the time-integrated Sm-Nd  (2017) and Koornneef et al. (2017) i.e., inclusions were divided into groups with comparable T CHUR ages of < 0.4 Ga, < 0.75 G a, < 1.15 Ga, < 2.7 Ga and > 2.7 Ga and further subdivided based on mineral composition, 87 Sr/ 86 Sr and the calculated enrichment of LREE relative to MREE using (La/Gd) N of the parental fluids (Fig. 7).

Peridotitic diamonds at Letlhakane
The T CHUR model ages (≥ 2.9 Ga) of the three harzburgitic garnets (Group a) provide the first age constraints for harzburgitic diamond formation at Letlhakane (and the entire Orapa cluster). Diamond formation apparently resulted in metasomatic re-enrichment of an originally highly meltdepleted protolith.  ) and Lu-Hf, Sm-Nd and Re-Os peridotite xenolith data in the Kalahari Craton (Branchetti et al. 2021;Carlson et al. 1999;Simon et al. 2007). The variable but high Cr 2 O 3 contents (12-17 wt%) of the LK peridotitic inclusions suggest different degrees of original source depletion and/or varying pressure in diamond formation (Deines et al. 2009;Griffin et al. 1999;Stiefenhofer et al. 1997). Application of the Cr/Ca-in-garnet barometer of Grütter et al. (2006) using a 40 mW/m 2 geotherm yields 5.9-6.1 GPa (~ 200 km depth) for LK058 gnt A, B and 7.8-8.0 GPa for the high-Cr LK101 gnt A, B (~ 260 km depth). There are potentially increased uncertainties from microprobe measurements when using unpolished surfaces of the inclusions; however, these pressure estimates are comparable to those recorded in a harzburgitic garnet from the nearby Karowe mine (7 GPa; Motsamai et al. (2018)). Deep formation/storage near the base of the lithosphere is supported by the inferred T MR of 1195 °C in LK058 (66 at. ppm N with 71%B) that overlaps with equilibration temperatures calculated from olivine-garnet diamond inclusion thermometry from Orapa (1230 ± 80 °C; Stachel et al. (2004b)), Karowe (1200 °C, 5.2 GPa; Motsamai et al. (2018)) and the inferred geophysical structure of the SCLM beneath Letlhakane (Griffin et al. 2003;Shirey et al. 2002).
The LK inclusions overlap with the radiogenic 87 Sr/ 86 Sr ratios of garnet harzburgite xenoliths (0.7091) and the Sm-Nd model age (T CHUR : ~ 3.2 Ga) that were interpreted to demonstrate the ancient origin of the lithospheric root beneath Letlhakane (Stiefenhofer et al. 1997). The trace element contents of the harzburgitic garnet inclusions (Online Resource 2: Table ESM_2.5; Fig. 4a) indicate metasomatic enrichment of a depleted protolith, as the garnets have low Zr (24-48 ppm) and Y contents (1.8-3.5 ppm). They are characterised by high LREE and low HREE ((La/Gd) N 295-467) as has been observed in garnet inclusions from Orapa and Karowe (Motsamai et al. 2018;Stachel et al. 2004b). Overall, harzburgitic diamonds (86% of peridotitic diamonds) are widespread in the SCLM (Stachel and Harris   Fig. 6 Sm-Nd isochron compilation of individual harzburgitic, eclogitic and websteritic garnet (square) and clinopyroxene (circle) inclusions from Letlhakane, Orapa (stippled markers) and Damtshaa (dotted markers) in this study and literature data for Letlhakane and Orapa (white diamonds) from Timmerman et al. (2017). Note: εNd i of the isochron in reference to CHUR; individual colours represent cogenetic inclusions Fig. 7 The evolution of eclogitic (E-type), websteritic (W-type) and harzburgitic (H-type) diamonds from the Orapa kimberlite cluster in context with the geodynamic setting (see Sect. 4.4.1) in southern Africa. Isochron ages (coloured font) refer to diamond formation in Letlhakane (LK), Orapa (OR) and Damtshaa (DM) with the number of dated inclusions (n) in brackets. Isochron ages from Orapa and Letlhakane (grey font) are from Timmerman et al. (2017), from Venetia (grey font) from Koornneef et al (2017). Box and whisker plots for carbon isotopes and nitrogen content together with general properties of 87 Sr/ 86 Sr composition, T CHUR model ages, mineral equilibration temperature (T K ) and integrated mantle residence time (T MR ) relate to the corresponding age population. Note: the figure is not to scale 1 3 54 Page 16 of 25 2008) but appear rare in Botswana. Two small (< 2 µg) lherzolitic cpx were recovered from the core of LK200 that has relatively low N content (~ 50 at. ppm). The N aggregation is above 80% but no firmer age constraints can be made.

Age group
Page 17 of 25 54 1698 ± 280 Ma). An additional age of 751 ± 150 Ma was found for eclogitic garnet inclusions at Letlhakane. No new eclogitic inclusions were found associated with the reported 2322 ± 22 Ma age (Fig. 6). Initial ratios of all the isochron ages are close to CHUR with ɛNd i ranging from − 2.9 to 3.3 (Table 3). Additionally, two less well-constrained eclogitic events (Groups b and c) are indicated with notable un-/ radiogenic 143 Nd/ 144 Nd. Combining eclogitic ages with data of websteritic inclusions from Letlhakane (and subsequently with eclogitic inclusions from Orapa and Damtshaa) suggests simultaneous and episodic diamond formation in both diamond parageneses. Hence, the following section will critically assess the diamond and inclusion properties of the different petrogenetic diamond suites to evaluate a genetic relationship and constrain the scale and nature of diamondforming processes (also see Online Resource 1.3: Tectonomagmatic implications).

Archaean eclogitic diamond formation Older than 2.7 Ga (Group b)
The Sm-Nd isotope systematics of the eclogitic garnet LK224 are extreme but the age poorly constrained (T CHUR = 3.9 ± 1.2 Ga). The inclusion has supra-chondritic REE contents with LREE depletion, flat MREE/HREE and a negative Eu anomaly (Eu/Eu* = 0.25, Fig. 4b). The existence of Archaean eclogitic diamonds (> 2.5 Ga) has previously been indicated in the Kaapvaal Craton (Pearson and Harris 2004;Richardson et al. 2001Richardson et al. , 2004Shirey et al. 2008;Westerlund et al. 2004). The diamond is low in N (< 20 at. ppm) and does not allow firmer constraints on its age. The moderately light δ 13 C (− 9.2‰) in combination with the negative Eu anomaly suggests the involvement of recycled crustal component(s) in the source of the metasomatic fluid or eclogitic protolith (Aulbach and Jacob 2016).

Paleoproterozoic diamond-forming events At 2341 ± 21 Ma (Group 1)
Websteritic LK338 records a complex growth history (Fig. 1g) and its inclusions were recovered from distinct growth zones. Three cpx from its core have T CHUR ages between 1.4 and 2.0 Ga with moderately radiogenic 87 Sr/ 86 Sr (0.706). In combination with the 2.3 Ga eclogitic inclusions (0.7045 and 0.7097) of Timmerman et al. (2017) that have similar Sm/Nd ratios (0.2 and 0.3), these data imply diamond formation at 2341 ± 21 Ma with an εNd i of + 5.9 (Table 3). The websteritic LK338 core has light δ 13 C (− 25.6‰) as do the eclogitic inclusions (δ 13 C − 33.6‰ to − 9.2‰) and suggest the involvement of subducted and recycled (Archaean) organic material (Cartigny et al. 1999;Timmerman et al. 2017). This conclusion is re-enforced by heavy δ 15 N in the cores (Fig. 2a) of LK42 (+ 2.3‰ ± 0.9‰) and LK47 (~ + 9.3‰). Because no LREE data are available for these eclogitic inclusions, only the MREE and HREE of the inferred diamond-forming metasomatic fluid compositions can be compared. Calculated (Dy/Yb) N ratios for the eclogitic (5-9) and websteritic (8-29) inclusions overlap and potentially suggest simultaneous interaction of a fluid or series of fluids of similar composition with eclogitic and websteritic portions of the SCLM. The similar T MR of websteritic LK338 (core: 259 at. ppm N, 63%B; T MR : 1155 °C) and these eclogitic diamonds (1119-1165 °C) imply comparable storage conditions and thus, suggest that the Eand W-type diamonds precipitated at a similar level in the SCLM.
At 1698 ± 280 Ma in Orapa (Group 2) Four eclogitic Orapa cpx with unradiogenic 87 Sr/ 86 Sr (0.7032-0.7046) define an Sm-Nd age of 1698 ± 280 Ma and indicate the involvement of isotopically depleted source material with εNd i of + 3.3 (Table 3). The diamonds are low in N (< 20 at. ppm) and have mantle-like δ 13 C (− 7.5‰ to − 5.6‰). Only cpx inclusions define this event in comparison to the predominance of garnet inclusions that define the other events, but notably, the cpx have comparable T CHUR model ages (~ 2 Ga). As yet, there is no evidence for this event at Letlhakane or Damtshaa, suggesting localised diamond formation. This finding could, however, represent sampling bias towards diamonds carrying multiple, large silicate and sulphide inclusions.
Combining data from the entire Orapa kimberlite cluster constrains the time that eclogitic and websteritic diamonds formed to 1100 ± 64 Ma (n = 11) with εNd i − 1.7. Notably, the age does not change significantly when different parageneses are combined but the uncertainty improves. The near Bulk Earth initial ratio of the isochron suggests the presence of recycled crustal material in the source of the diamond-forming fluids, which is supported by the light δ 13 C signatures observed in the eclogitic (− 20.5‰ to − 7.0‰) and websteritic (− 25.0‰ to − 19.8‰) diamonds (Online Resource 2: Table ESM_2.4) and variable Eu anomalies (Eu/ Eu* = 0.7-0.9). The integrated residence temperatures of corresponding N-bearing diamond growth zones are comparable (T MR = 1109-1206 °C) and indicate similar conditions for subsequent storage in the SCLM. These consistent temperatures suggest that subsequent diamond-forming events that clearly affected some samples (e.g., LK338) have not significantly increased the average T MR , implying formation of the rims of these diamonds was not associated with major regional heating events.

Phanerozoic diamond-forming events At 220 ± 80 Ma (Group 5)
Eclogitic garnets from Letlhakane (LK42 rim, LK75 rim, LK362 rim) define an age of 293 ± 250 Ma but show an extended range in (La/Gd) N fluid ratios (20-336) and light δ 13 C (− 21.3‰ to − 14.3‰). Including the websteritic garnet LK338 ((La/Gd) N 9; δ 13 C − 25.7‰) and cpx LK40-1 (δ 13 C − 20.5‰) recovered from rims of the diamonds, improves the isochron uncertainty to 294 ± 110 Ma with εNd i + 5.7. Notably, diamond LK338 is a rare example providing evidence for episodic websteritic diamond growth, potentially spanning > 2 Gyr with inclusions lying on the 2.34 Ga, 1.10 Ga and 0.22 Ga isochrons. Such protracted diamond growth has only been established in a few diamonds Timmerman et al. 2017;Wiggers de Vries et al. 2013). The websteritic mineral equilibration temperature, T K (Krogh 1988), from LK40 rim (N < 200 at. ppm) yields a diamond formation temperature of T K = 1113 °C that overlaps with integrated mantle residence temperatures of 1119 °C to 1165 °C from the eclogitic diamonds (N > 200 at. ppm) suggesting that diamond formation was not associated with a major thermal pulse and subsequent cooling.
Younger than 140 Ma (Group c) Two eclogitic garnets from Letlhakane (LK170) have unradiogenic 143 Nd/ 144 Nd (0.5119-0.5123) and 147 Sm/ 144 Nd above 1 (Fig. 6). They have negative T CHUR ages and are outside the uncertainty of the 220 ± 80 Ma isochron implying extremely recent formation (< 140 Ma). A websteritic inclusion (LK346) has similar characteristics and the fluids inferred to have formed these inclusions are highly enriched in LREE ((La/Gd) N 24-108). Further evidence for very recent diamond formation is given by three garnets from Orapa ( 87 Sr/ 86 Sr 0.7054-0.7081) and two from Damtshaa ( 87 Sr/ 86 Sr 0.7055-0.7101). One of these samples (DM012 B) has an extreme 147 Sm/ 144 Nd of 2.98. The negative εNd i (− 42.5 to − 2.2) of the inclusions (Table 2) indicates diamond formation involved a protolith and/or metasomatic fluids with time-integrated LREE enrichment such as Archaean to Proterozoic eclogites or orangeites/lamproites (εNd i − 28 to 0; see Fraser et al. (1985); Pearson et al. (2019)). These data are also consistent with the involvement of highly trace element enriched sources in the SCLM, such as stalled kimberlitic melts . Nitrogen systematics of the host diamonds appear consistent with young ages (< 140 Ma) and yield high T MR (1226-1246 °C) suggesting deep formation or involvement of recent hot magmas. Stable isotope compositions (Fig. 2a, b) record a wide range in δ 13 C (− 21.3‰ to − 4.9‰) and δ 15 N (7.5‰; LK170 rim) suggesting the involvement of a wide variety of fluids and protoliths of crustal and mantle origin.

Diamond formation beneath the Orapa cluster
The Sm-Nd isotope data demonstrate that the timing of websteritic diamond formation is simultaneous with eclogitic diamond formation at 2341 ± 21, 1110 ± 64, 220 ± 80 Ma and < 140 Ma extending over an area of ~ 200 km 2 beneath the Orapa cluster. The synchronous formation of E-and W-type diamonds on this scale places constraints on general diamond-forming processes. A fundamental question is whether a single metasomatic fluid was responsible for diamond formation on a regional scale (10's of km) or whether low-volume fluids were mobilised locally (on a scale of metres to kilometres) but contemporaneously in several locations (i.e., dispersed) across the Orapa cluster. Similar considerations were recently made by the authors (Gress et al. 2020) for the Jwaneng mine, 360 km southwest of the Orapa cluster. Metasomatic fluids derived from the asthenosphere or within the SCLM could have originated from peridotitic, websteritic or eclogitic sources. Such metasomatic sources may account for fluids of varying compositions when infiltrating the overlying lithosphere.
Websterites are commonly considered to be transitional between peridotites and eclogites, representing the result of fluid-rock interaction (Aulbach et al. 2002;Rapp et al. 1999). The geochemical systems involved in diamond petrogenesis can therefore be considered to have involved a minimum of three components (Fig. 8); peridotitic mantle (A); eclogites dominated by mafic material (B); and eclogites that included recycled sedimentary material (C). Simultaneous W-and E-type diamond formation would have produced compositional arrays as the diamond-forming metasomatic fluid(s) would have undergone fluid-rock reactions with different protoliths. Based on C-N-O-S isotope heterogeneity in eclogitic xenoliths and E-type inclusions, previous workers have also argued for the evolution of metasomatic fluids (Deines and Harris 2004;Ickert et al. 2013;Smit et al. 2019;Thomassot et al. 2009). The inclusion chronology enables us to examine the inferred mixing processes and to assess any temporal variations (Figs. 7 and 8).
The dated W-type inclusions have restricted Al/Cr (< 50) and light δ 13 C − 25.7‰ to − 19.8‰ in the diamonds compared to the wide range in E-type diamonds. Despite the episodic growth history (Fig. 1g), adjacent growth zones in the websteritic diamonds record little carbon isotope variation (< 1.5‰; Online Resource 3) and inclusions have similar parental fluid compositions (Fig. 5). Comparable T MR and T K demonstrate insignificant thermal impact from more recent diamond-forming events and are consistent with W-type diamond formation at depth of 140-180 km as suggested by Deines et al. (1993). In this setting, the carbon isotope compositions of different growth zones seem buffered by the local W-type protolith, suggesting carbon was (re-)mobilised under conditions of relatively low fluid/rock ratios.
The generation of the websteritic protolith must have occurred before 2.34 Ga i.e., the age of the oldest websteritic diamond (LK338), and potentially involved organic-rich subducted material (component C) mixing with mantle material (component A) because eclogitic 2.34 Ga diamonds with low to moderate Al/Cr (< 600) and variable δ 13 C (− 33.6 to − 9.0‰) also fall along this hypothetical mixing array, as does the > 2.7 Ga eclogitic garnet (Fig. 8). The REE enrichment of the inferred parental fluid compositions of these E-and W-type inclusions with the same age are consistent with a genetic link (Fig. 5).
The bi-modal δ 13 C distribution in eclogitic diamonds (Online Resource 1: Figure ESM_1.1) may be related to the involvement of different fluids (e.g., mantle vs recycled) and eclogitic protoliths which they infiltrate. Diamonds with lighter δ 13 C (mode at − 16.2‰) potentially have a common metasomatic fluid with the W-type diamonds whose protoliths formed in the transition from the Archaean to the Paleoproterozoic. E-type diamonds with a mode at − 6.6‰ represent mixtures of mantle component A with material predominantly derived from altered oceanic lithosphere (component B) that seem to represent the dominant source for E-type diamonds since the Mesoproterozoic. Simultaneous E-and W-type diamond formation was caused by fluids that migrated through different protoliths and were modified by local fluid/rock interaction. Depending on the prevailing buffering capacity of the protoliths (eclogite > > websterite > > peridotite) at specific times and depth (e.g., fO 2 ; Burness et al. (2020); Woodland and Koch (2003)) the diamonds and inclusions would preferentially record component A mixed with either with component B or C.
Due to insufficient data on mineral equilibration temperatures in our sample set, we cannot yet verify that Letlhakane diamonds with δ 13 C < − 8‰ originate from a specific part Fig. 8 Time-dependant changes in Al 2 O 3 /Cr 2 O 3 of garnet (square) and cpx (circle) inclusions versus the carbon isotope (δ 13 C in ‰) composition of its corresponding diamond growth zone after Deines et al. (1993) and Deines et al. (2009). The figure illustrates estimated mixing trends between mantle (component A) with no crustal signature (low Al 2 O 3 /Cr 2 O 3 and heavy δ 13 C (~ − 5.0 ‰ ± 1 ‰), hydrothermally altered subducted oceanic lithosphere (component B) with variable amounts of Al-and Cr-oxide and δ 13 C (> − 24 ‰) and organic-rich material (component C) represented by LK33 with variable amounts of Aland Cr-oxide and δ 13 C (< − 24 ‰). Note: individual colours represent co-genetic inclusions, not dated samples from this study (grey fill), literature (white fill). Data sources as in Fig. 3 of the mantle (1125-1240 °C and depth of ~ 150-170 km) as suggested by Deines et al. (2009) but we can confirm W-and E-type diamond formation included a significant subducted organic-rich (Archaean) component. These observations support models that propose diamond growth under various redox conditions e.g., from oxidised CO 2 -/carbonate-bearing or reduced CH 4 -bearing fluids (Shirey et al. 2019), which complicates assessment of potential temporal changes in carbon sources.

Implications for the deep carbon cycle
Less than 20% of the studied diamonds have homogeneous C isotopes and N systematics (Online Resource 3) associated with oscillatory growth that may indicate formation within a geologically short period of time (millions of years; Bulanova et al. (2014)). The difference in T MR (> 100 °C) and the range in carbon isotopic composition and nitrogen content of these diamonds suggests that they are not all cogenetic and resided at different depths and conditions in the SCLM.
Some diamonds (e.g., LK241, LK272; Fig. 1b, e) record systematic variations in C isotopes and N systematics that could be ascribed to small changes in the conditions (P-T-fO 2 ) under which the diamonds precipitated (Boyd et al. 1992;Deines et al. 1989) while diamonds with abrupt compositional changes between major growth zones (Fig. 1f) clearly demonstrate formation of adjacent growth zones from different metasomatic fluids e.g., in LK09 δ 13 C varies from − 10.0‰ (core: 1110 ± 64 Ma; T MR : 1109 °C) to − 5.0‰ (rim). Review of the global dataset establishes that limited intra-diamond C isotope variation is typical (96% have < 3‰) suggesting that the oxidation state and carbon speciation of diamond-forming fluids have been similar throughout history (Howell et al. 2020;Jablon and Navon 2016;Weiss et al. 2014). Local variations do exist, however, and, potentially in combination, changes in the scale and mechanism of diamond formation, may be depth or time dependent.
Based on his previous work, Deines et al. (2009) concluded that light δ 13 C values of bulk diamonds from the Orapa cluster were more prevalent with increasing mantle depth (> 200 km). Depletion in LREE of the E-and W-type garnets (Fig. 4b, d) is similar to that found in eclogitic xenoliths from the base (170-220 km) of the Kalahari Craton while minerals from shallower levels (140-180 km) were found to be more LREE enriched (Burness et al. 2020). A δ 13 C depth dependence was also proposed at Argyle (Australia), where eclogitic diamonds with lighter carbon and lower 3 He/ 4 He ratios are mostly formed at the base of the lithosphere where subduction input is largest (Jaques et al. 1989;Timmerman et al. 2019).
A correlation of δ 13 C values with relative time has been proposed based on rim growth zones of diamonds that record heavier δ 13 C than their cores Chinn et al. 2018;Thomson et al. 2014). Moreover, the internal range in δ 13 C values was found to become smaller (towards heavier values) in dated diamonds from 3.0 to 2.1 to 1.0 Ga (Howell et al. 2020), indicating some homogenization of carbon isotope compositions in the SCLM. Globally, significant changes in δ 13 C values within individual diamonds (4% has > 3‰) are rare (Howell et al. 2020). In contrast to this general trend, diamonds from Letlhakane have large internal shifts in δ 13 C values (33% have > 3‰) and they predominantly (80%) evolve towards heavier δ 13 C values from core to rim (Online Resource 2: Table ESM_2.3). This observation either supports a relative time-related dependence for the carbon isotope composition (Timmerman et al. 2017) and may indicate fluid input from a different part of the subducting slab and/or increasing input of carbonates with time or some homogenization of carbon in the SCLM (potentially restricted to a local scale) resulting in a higher relative sampling of heavier δ 13 C during diamond rim growth. In the latter case, the change of δ 13 C then is only dependent on the prevailing carbon composition of the local SCLM.
Overall, δ 13 C values of diamond hosts show no systematic variation with absolute inclusion age (Figs. 7 and 8). Largescale events in the region, such as the Bushveld Large Igneous Province (LIP) at 2.06 Ga or Proterozoic crustal accretion such as the Rehoboth Terrain between 2.0 and 1.75 Ga are capable of modifying major parts of the SCLM and may lead to a redistribution of volatile elements . In this scenario, it appears that the observed variability in δ 13 C and inclusion composition in the studied diamonds has little systematic control and the deep carbon cycle seems controlled by a combination of time and protolith variability, which in turn may be partly depth controlled.

Conclusions
Eclogitic and websteritic silicate inclusion populations from Letlhakane record multiple and contemporary diamondforming events. In combination with inclusions from the nearby Orapa and Damtshaa mines, it is established that eclogitic and websteritic diamond formation occurred beneath the Orapa kimberlite cluster at 220 ± 80 Ma, at 746 ± 100 Ma, at 1110 ± 64 Ma, at 1698 ± 280 Ma and at 2341 ± 21 Ma in events that are related to major tectonomagmatic perturbations in the SCLM (see Online Resource 1 for greater discussion). Additionally, two less-well-defined events close to kimberlite eruption (< 140 Ma) and during craton stabilisation (> 2.7 Ga) can be established. Harzburgitic diamond formation in the roots of the SCLM also tope measurement by LG-SIMS greatly benefit from technical discussion with Johan Villeneuve. ET thanks Dorrit Jacob and Daniel Howell for generous supply of the reference material to aid development of high-resolution C and N isotope analyses at CRPG-Nancy. MUG also gratefully acknowledges Karen V. Smit and Wuyi Wang, Gemological Institute of America, for access to perform additional FTIR analyses during an internship realised through the American Immigration Council. Richard Smeets, Suzan Verdegaal-Warmerdam, Bouke Lacet, Ciaran Kelly, Peter van Krieken, Sergei Matveev and Tilly Bouten are thanked for assistance with analyses. Pieter Ouwerkerk and Gassan Diamond B.V. is thanked for help with polishing. The Dr. Schürmannfonds foundation provided support for field work to Botswana during which sample selection was made with support of Quint van den Heuvel, Ellen Schulten and Cas Nooitgedacht. Thomas Stachel is thanked for providing access to the De Beers diamond laboratory at the University of Alberta, Edmonton, for FTIR analyses. The manuscript benefitted from constructive comments of Stephen Richardson, an anonymous reviewer and editor Dante Canil.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.