Diamond ascent by rift-driven disruption of cratonic mantle keels

Diamonds are erupted at Earth’s surface in volatile-rich magmas called kimberlites 1,2,3 . These enigmatic magmas, origi-nating from depths exceeding 150 kilometres in Earth’s mantle 1 , occur in stable cratons and in pulses broadly synchronous with supercontinent cyclicity 4 . Whether their mobilization is driven by mantle plumes 5 or mechanical weakening of cratonic lithosphere 4,6 remains unclear. Here we show that most kimberlites spanning the past billion years erupted approximately 25 million years after the onset of continental fragmentation, suggesting an association with rifting processes. Our dynamic models show that physically steep lithosphere-asthenosphere boundaries formed during terminal rifting (necking) generate convective instabilities in the asthenosphere that slowly migrate many hundreds of kilometres inboard of the rift, causing destabilization of cratonic mantle keel tens of kilometres thick. Displaced lithosphere is replaced by hot, upwelling asthenosphere in the return ﬂow, causing partial melting of carbonated mantle and variable assimilation of lithospheric material. The resulting small-volume kimberlite magmas ascend rapidly and adiabatically, exsolving amounts of carbon dioxide (CO 2 ) that are consistent with independent constraints 7 . Our model reconciles diagnostic kimberlite features including association with cratons and geochemical characteristics that implicate a common asthenospheric mantle source contaminated by cratonic lithosphere 8 . Together, these results provide a quantitative and mechanistic link between kimberlite episodicity and supercontinent cycles via progressive disruption of cratonic keels.

temporal dependencies in the global tectonic cycle. 23 [2] We assessed the link between kimberlites and global tec- 24 tonics through geologic time using a database of well-dated 25 kimberlites 6 and a measure of the degree of fragmentation of 26 continental plates from tectonic reconstructions 16 . We calcu- 27 lated the rate of change of fragmentation (∆F) as a proxy for 28 dynamic plate reorganization (Methods; Extended Data Fig. 29 1), and then calculated multi-million-year time lags 17 between 30 ∆F and the kimberlite count, K, at 1 million year (Myr) in-31 tervals (Methods; Fig. 1b). Our cross-correlation analysis re-32 veals a statistically significant association between fragmenta- 33 tion and kimberlites over the past 500 Myrs (Fig. 1b), which 34 persists when we account for auto-correlation in the time series 35 (Extended Data Fig. 2a). The strongest correlation coefficient 36 (ρ=0.41) prevails at lags of -26 ± 4 Myr, indicating that conti- 37 nental fragmentation typically leads kimberlite magmatism by 38 22-30 Myr (Fig. 1b). When we extend our data compilation 39 to 1 Ga (Extended Data Fig. 1), using more uncertain data but 40 cycle.  . Temporal relationships between tectonics and kimberlites | a, Frequency distribution of kimberlites through geologic time (data from ref. 6 ), showing peaks coinciding with the breakup phase of supercontinent cycles (from ref. 18 ). b, Cross-correlations between ∆F and global kimberlites 6 spanning 500-0 Ma (Methods), showing dominant time lags at -26 ± 4 Myr (i.e., fragmentation leading kimberlites); positive correlations show that kimberlites are linked to continental fragmentation, not assembly; dashed blue lines show 95% confidence intervals. c, Time lags between emplacement of African kimberlite clusters and onset of local rift separation. d, Box and whisker plot showing time lags between rift separation and kimberlite eruption versus distance from Atlantic rifted margins (n=45; IQR: interquartile range; med: median). Shaded pink field shows a steady increase in distance between rifts and kimberlites in the ∼30 Myrs following breakup. [3] We scrutinized this linkage further by analyzing spatial 49 and temporal lags between kimberlites and continental plate 50 margins, targeting the Mesozoic kimberlite fields of southern 51 Africa (Fig. 1c). This region is perfectly suited for our purposes 52 because the rift kinematics along adjacent plate boundaries are 53 well constrained at this time 19 , and the abundant kimberlites 54 are well understood in terms of their age distributions 6 and 55 the structure of lithosphere they sample 20,21,22 . Using plate re- 56 constructions, we measured distances between kimberlites and 57 adjacent rift boundaries (Extended Data Fig. 3) and calculated 58 time lags between rift separation and eruption (Methods). This 59 analysis confirms that a peak in kimberlite emplacement occurs 60 25-27.5 Myr after the onset of regional breakup (Fig. 1c), cor- 61 roborating the global result (Fig. 1b). Importantly, we find that 62 kimberlites tend to erupt closer to the rift margins earlier in the 63 rifting cycle-as qualitatively noted before in the São Francisco 64 and Kaapvaal Cratons 11 -and migrate toward the cratonic in- 65 terior as fragmentation progresses, initially at a rate of about 66 18-23 km Myr −1 (Fig. 1d; Methods). 67 [4] Our results indicate that kimberlite magmatism is 68 strongly associated with continental breakup. However, both 69 breakup 23,24 and kimberlite magmatism 5 are commonly at- 70 tributed to mantle plumes. A crucial question, therefore, is 71 whether rifting itself is the primary driver of kimberlites, or 72 whether mantle plumes drive both rifting and kimberlite mag-73 matism. To address this question, we quantified the relation- 74 ship between rifting, plumes and kimberlites globally, using the 75 techniques described previously (Methods). We find that the 76 strongest statistically significant correlation between Large Ig-77 neous Provinces, LIPs (the accepted main surface expression of 78 mantle plumes) 24

132
To investigate the essential physics of the process, we carried where λ * d and q * d are analytically determined scales for the char-152 acteristic horizontal wavelength and growth rate respectively, 153 g ′ = g∆ρ/ρ is reduced gravity, ∆ρ is the density difference that   Table 1a). Models 1, 4 and 5 (black lines) are most applicable, and describe a dense, viscous layer representing the lithospheric keel that overlies a less dense, viscous half-space representing asthenosphere. A kinematic viscosity of 7 × 10 15 m 2 s −1 gives good agreement with constraints from kimberlite migration (Fig. 1d), numerical simulations of basal lithospheric instabilities (Fig. 2), and lithospheric root thicknesses inferred from modelling xenolith P-T data. b, Phase diagram showing carbonate-bearing lherzolite solidi (solid orange line) for proposed kimberlite starting compositions, JADSCM-7 and JADSCM-14 27 ; shaded area shows the CaO-MgO-Al 2 O 3 -SiO 2 -CO 2 (CMAS-CO 2 ) field of carbonatite and kimberlite melt generation from refs. 27,28 ; hatched box shows P-T conditions in the thermal boundary layer, derived from analysis of xenolith geotherms (Extended Data Fig. 5); Dol=dolomite, Mgs=magnesite, Sp=spinel, Grt=garnet, An=anorthite. The graphite-diamond phase boundary is from ref. 29 . Xenolith P-T estimates for various southern African 30   As we are interested in the correlation of (lagged) ∆F with 331 K t , and wish to remove the effects of shorter lags at each step 332 (e.g., see ref. 17 ), we define the Bayesian network node hierar- elling and observations (Fig. 1). We repeated the above proce-345 dure to analyze the relationship between ∆F and large igneous 346 provinces (LIPs) using the well-established ages of LIP mag-347 matism from ref. 24 (Extended Data Fig. 4). Here, the input is a in agreement with our models (Fig. 2). In regions such as the is located at least 400 km inboard of the rift (ref. 52  to our modelled propagation rates of Rayleigh-Taylor instabili-401 ties (Fig. 3a).  Table 1a).

427
To apply these analytical solutions, note that length is scaled Now we consider a single instability of dominant wavelength.

440
Since τ d is a e-folding time, soon after t = τ d the instability  gives a good match between this scaling law, numerical mod-477 els described below, and observed migration rates of kimberlite 478 magmatism (Fig. 3a).  25 . In further agreement with previous mod-492 els, breakup is delayed due to rift migration 57 .

493
We will now briefly describe the geometrical, mechanical   To plot estimated P-T conditions of llherzolite/peridotite 657 nodules (Fig. 3b), it was necessary to convert some depths from 658 published compilations to pressures, which was done using where ρ is the density of the lithosphere, g is acceleration due 662 to gravity and h is the depth. In our calculations, we took ρ = 135 km and density of 3,300 kg m −3 (i.e., using a middle value 667 from ref. 68 ).

668
Statistical analysis of isotope variations 669 We compiled existing data to assess whether any temporal

717
Additional references are cited in Extended Data 75,76 .

Data availability
All data generated or analyzed during this study will be provided in the online version of this article and in Extended Data Table 1.

Code availability
More details on the computational methods and tools used for this study are available from the corresponding author (Thomas.Gernon@noc.soton.ac.uk) upon request.