Kimberlite eruptions driven by slab flux and subduction angle

Kimberlites are sourced from thermochemical upwellings which can transport diamonds to the surface of the crust. The majority of kimberlites preserved at the Earth’s surface erupted between 250 and 50 million years ago, and have been attributed to changes in plate velocity or mantle plumes. However, these mechanisms fail to explain the presence of strong subduction signatures observed in some Cretaceous kimberlites. This raises the question whether there is a subduction process that unifies our understanding of the timing of kimberlite eruptions. We develop a novel formulation for calculating subduction angle based on trench migration, convergence rate, slab thickness and density to connect the influx of slab material into the mantle with the timing of kimberlite eruptions. We find that subduction angles combined with peaks in slab flux predict pulses of kimberlite eruptions. High rates of subducting slab material trigger mantle return flow that stimulates fertile reservoirs in the mantle. These convective instabilities transport slab-influenced melt to the surface at a distance inbound from the trench corresponding to the subduction angle. Our deep-time slab dip formulation has numerous potential applications including modelling the deep carbon and water cycles, and an improved understanding of subduction-related mineral deposits.


Scientific Reports
| (2023) 13:9216 | https://doi.org/10.1038/s41598-023-36250-w www.nature.com/scientificreports/ re-initiation of subduction is increasingly sparse with age and plate reorganisations render subduction longevity difficult to quantify. Similarly, the subdivision of slab segments through deep time is highly subjective and a break in subduction zone topology may result in vastly different slab dip estimates. To overcome these challenges, we explore slab dip correlations with a multivariate analysis of plate rheology and kinematic parameters that are less ambiguous through deep geological time due to the optimisation of plate models for no-net rotation 24 . By taking this approach, the estimation of slab dip can be applied to evolving subduction zones stretching back as far as a given plate reconstruction will allow. We extract present-day slab dip data from the Slab2 model 19 , which estimates the depth and geometry of subducting slabs across the globe from earthquake depths and tomography models (Fig. 1), and combine those with close to present-day plate kinematic properties obtained from a recent plate reconstruction model 25 using pyGPlates 26 . The dip angle, θ , of the down-going slab is taken at multiple depth intervals inboard from the trench (Fig. 2). The average dip angle, θ av is simply the arithmetic mean of all depth intervals resolved by the Slab2 model orthogonal to the trench, where d is the depth interval and n is the total number of depth intervals that may be resolved by the Slab2 model. The entire workflow to calculate the dip angle from the Slab2 model as well as plate rheology and kinematic parameters of the subducting plate are openly available on GitHub (https:// github. com/ brmat her/ Slab-Dip). The most statistically significant combinations of parameters that are sensitive to the present-day dip angle of subducting slabs are described in following sections.
Slab flux. Previous studies on the multivariate analysis of subduction coefficients did not find a statistically significant relationship between slab dip and the age of subducting oceanic lithosphere 8,9,13 . However, we find a significant relationship where the slab dip is proportional to the convergence velocity, v c ( P = 0.39 ), and the thickness of the down-going plate, h plate ( P = 0.59 ) (Fig. 3). Plate thickness was predicted from the thermal evolution of oceanic lithosphere, (1) θ av = 1 n n d=0 θ d Figure 1. Depths of subducting slabs obtained from the Slab2 model 19 16,27 . These models describe the thickening of the thermal boundary layer as a function of the age of seafloor, which approaches a maximum thickness around 80 Myr for a constant thermal diffusivity coefficient of κ = 1 mm 2 /s. We employ a plate model of oceanic lithosphere cooling, . v s and v o are positive towards the trench and each vector is orthogonal to the trench. v hsp is the half-spreading rate at mid-ocean ridges, which is proportional to the rate of volatile influx ( q v ) into the plate, and θ is the dip angle of the subducting slabs calculated at different intervals. Correlations between slab dip and convergence velocity, slab thickness, and spreading rate. 2D histograms illustrate the probability density between each parameter with slab dip for all subduction zone segments tessellated at 0.5 degree increments (a-c); subduction zone segments are grouped for each subduction zone with whiskers indicating one standard deviation from the mean (d-f). P r is the Pearson correlation coefficient value and P is the p-value. www.nature.com/scientificreports/ which has been demonstrated to produce an optimal fit to depth and heat flow data, to calculate the thickness of oceanic lithosphere being recycled at trenches across the Earth 16 . The product of the plate thickness, h plate , and the convergence velocity, v c , integrated along a trench segment gives the volumetric rate of lithospheric recycling, or "slab flux". We sample seafloor age grids at trench boundaries to calculate the thickness of subducting oceanic lithosphere.
Rollback. The second significant parameter to consider is the rate of trench advance or retreat. Slab rollback is widely regarded to lead to low-angle subduction and slab stagnation at the upper and lower mantle boundary 12,28,29 , which partly explains why lower subduction angles are more widely observed at subduction zones between two oceanic plates, where trench retreat is more common, than at ocean-continent subduction zones 10 . The rate of trench migration, v t , is calculated relative to the mantle reference frame and can be compared to the convergence velocity, v c , to characterise different subduction dynamics: • If v c = −v t the entire convergence rate is partitioned to rollback; • If v t = 0 the trench is stationary and the velocity of the subducting plate is equal to the convergence rate ( v s = v c ); • If v t > 0 the trench is advancing in the direction of subduction.
Volatile enrichment. A third parameter we considered is the volatile enrichment of the subducting plate.
An increased abundance of volatile components, such as water and carbon, enhances the coupling between the subducting plate to the overriding plate 30 . The source of volatile enrichment occurs at mid-ocean ridges, where volatile-bearing melt circulates through channels in newly-forming oceanic lithosphere 31 . The sequestration rate of volatiles within oceanic lithosphere, q v , is proportional to the seafloor spreading rate, v hsp 31 . We sample seafloor spreading rate grids, generated using workflows in 32 , at trench boundaries in the same way that we interpolate age grids to calculate plate thickness, h plate . We find that v hsp exhibits a strong negative correlation with slab dip ( P = 0.34 , Fig. 3c,f). Volatiles are also added to the plate by other processes, such as hydrothermal alteration of the upper crust 33 , sepertinisation during ultraslow spreading 34 , and during bending and cracking of the plate before entering the subduction zone 35 , however, the spreading rate relationship we outline above constitutes the most general source of volatiles subducted at most trenches globally.

Buoyancy.
A fourth parameter we considered was the mean density of the subducting plate. Lithospheric buoyancy can be estimated from plate age, thermal structure, crustal thickness, and depletion 36 . The crust is the main source of positive buoyancy in oceanic lithosphere due to its relatively low density ( ∼ 2900 kg/m 3 ) 37 . Assuming modern mantle temperature conditions, 10-20 Myr old lithosphere is negatively buoyant 36 and becomes more negatively buoyant with age. A 60 Myr segment of oceanic lithosphere is approximately 79.4 km thick 16 , 7 km of which is crust, which equates to a mean density of ρ av = 3278 kg/m 3 . While the buoyancy of oceanic plates may not be sufficient to initiate subduction alone 38 , lateral changes in the density structure of the down-going plate change the buoyancy force in established subduction zones which may affect slab dip. Such positive buoyancy anomalies are associated with oceanic plateaus which often congest subduction zones 39 or lead to flat slab subduction 40 .

Estimating slab dip
Taking all of these rheological and kinematic relationships into consideration for present-day subduction zones, we applied a nearest-neighbour regression to predict the dip angle, θ av , of subducting slabs. This regression implements a k-d tree to efficiently lookup k number of neighbours with the shortest euclidean distance from the training dataset X train to the test dataset X test , and takes the weighted mean of corresponding slab dips to predict the test slab dip.
where d k is the euclidean distance between the training and test data for k nearest neighbours, Using a subset of the present-day configuration of subduction boundaries and slab dips obtained from the Slab2 model as the test dataset, we compare the performance of the slab dip prediction against the training dataset (Fig. 4). The training score ( R 2 value) measures the closeness of fit between the training and test dataset which is a maximum of 1 for k = 1 . This is because there is an exact match for the test dataset from the training dataset at the present day. As k gets larger, more neighbours are incorporated within the average which reduces the R 2 value and the cross-validation score. The combination of these two metrics evaluate the estimator performance such that the problem is not over-fit or under-fit. We opt for k = 5 where there is an optimal trade-off between the cross validation and training scores.
We have developed a flexible object-oriented Python package to estimate slab dip based on present-day plate rheology and kinematic parameters (https:// github. com/ brmat her/ Slab-Dip). The default regression algorithm and training dataset have been described above, however, the user can also define their own training dataset and any regression algorithm bundled within the scikit-learn Python package to create a bespoke slab dip estimator. The GitHub repository has several examples provided as Jupyter notebooks for the user as well as installation instructions. The nearest neighbours regression we have chosen is general and establishes a robust link between present-day properties of subducting plates with those through deep time. While it may not be applicable in specific geodynamic contexts that differ from those operating at the present day, it encapsulates some of the primary www.nature.com/scientificreports/ drivers of subduction that are predicted from plate tectonic theory such as the relationship between seafloor age, slab thickness, trench migration, and convergence velocity. This recognises that the slab suction force, associated with older dense slabs which sink sharply into the mantle, drives the velocity of the tectonic plates at the surface of the Earth, and an increase in the amount of rollback leads to a flattening of the down-going slab.
It is important to note that the relationship we have formulated predicts the slab dip down to the maximum depth resolved by the Slab2 model. The trajectory of the down-going slab into the mantle may deviate from the predicted slab depths as it encounters viscosity contrasts in the mantle which sometimes leads to slab stagnation 41,42 and potentially slab anchoring 43 . These dynamics are not captured by our regression which may preclude its application to specific subduction zones over a nominal time range. Nevertheless, our slab dip formulation holds for most subduction zones globally and poses an advantage that it may be applied back through deep geological time using tectonic plate reconstructions to predict the trajectory of subducting slabs into the mantle. This can be used to estimate the distance between trenches and volcanic arcs, the incidence of flat slab subduction, and the recycling of oceanic lithosphere into deeper portions of the mantle.

Subduction controls on kimberlite eruption
Kimberlites are volcanic rocks which ascend rapidly from the mantle and are emplaced in cratons worldwide 2 . Kimberlite eruptions have been associated with mantle upwellings from the LLSVP 4 and extensional tectonics associated with lithospheric unloading 44 or changes in plate velocity 3 . However, these mechanisms do not explain the strong subducted slab signatures observed in Cretaceous kimberlites in Africa, Brazil, and North America from increased strontium-isotope ratios 1,5 . While African kimberlites exhibit a statistically significant relationship with the distance to LLSVPs, kimberlites in North America hold no such relation (Fig. 5). In contrast, kimberlite eruptions in North America have been linked to flat subduction of the Farallon plate during the Laramide Orogeny 45 . Here, magma may have been generated from water-fluxed decompression melting of the mantle transition zone 46 and transported upwards by subduction-induced mantle return flow 47 . High slab flux has been previously connected to volcanic eruption frequency, where mantle upwellings are driven by large volumes of oceanic lithosphere being subducted into the mantle 48 . Subduction-driven mantle upwellings have been linked to the formation of the 260 Ma Emeishan large igneous province in SW China 49 and Cenozoic volcanism in NE China 50 . To reconcile the role of subduction in generating distinct kimberlite populations in Africa and North America, we separated the vertical and horizontal components of slab flux using our formulation of slab dip in the previous section, and reconstructed global subduction boundaries using pyGPlates 26 . We used a 170 Myr plate model, modified from 15 , with significantly improved subduction zone boundaries along the western margin of North America, and improved resolution of the Caribbean plate 51 (refer to Methods). The locations of kimberlites ≤ 170 Ma were reconstructed back to their time of eruption, using a compilation of kimberlites with eruption ages 52 that we resampled to a 3-times refined icosohedral mesh 53 to avoid duplication and geographic sampling bias. We then separated the kimberlites into (i) North America and (ii) Africa populations. Combined, these kimberlite populations constitute 91% of the global kimberlite dataset which have erupted within the last 250-50 million years termed the "kimberlite bloom" 3 . In comparing the timing of kimberlite eruptions with rates of slab flux, we only consider trenches where subduction is in the direction of kimberlite populations. We find a strong correlation between high slab flux along the western margin of North and Central America, associated with the subduction of the Farallon plate, with both African and North American kimberlite populations (Fig. 6a). The low subduction angle (30-35 • ) predicted from our slab dip analysis on reconstructed subduction . Cross-validation score and training score ( R 2 value) for the nearest-neighbours regression algorithm to predict the dip angle of subducting slabs. The algorithm is tested using increasing numbers of nearest neighbours (k) used in the calculation from Eq. 3. Shaded regions indicate one standard deviation from the mean (solid lines). Cross validation score is shown as negative number for visual clarity to compare the optimal trade-off with the training score. www.nature.com/scientificreports/ zones indicates that slabs would extend more than 1,000 km from the trench before intersecting the 660 km mantle transition zone, and may penetrate deeper into the lower mantle. Rapid subduction of slab material can produce mantle return flow 35 , from which mantle upwellings may drive kimberlite eruptions. In the following sections we explore how subduction-induced mantle return flow from high rates of slab flux may be connected with kimberlite eruptions within Africa and North America.  www.nature.com/scientificreports/ African kimberlites. A strong correlation exists between slab flux and the frequency of kimberlite eruptions during the peak in African kimberlite eruptions between 120 and 130 Ma (Fig. 6b). Subduction persisted along the western margin of the Pangea supercontinent during its assembly until rifting commenced between Africa and South America at approximately 120 Ma (Fig. 7). From 160 to 120 Ma, two peaks in kimberlite eruptions correlate with pulses of high slab flux from the subducting Farallon plate beneath the Americas which plunged 30-35 • into the mantle. The larger peak in kimberlite eruptions at 120 Ma corresponds to a slab flux of 60 km 3 /yr. A second peak in African kimberlite eruptions occurred between 80 and 90 Ma, which correlates with a second pulse in slab flux (up to 80 km 3 /yr) and a maxima in plate velocity (6 cm/yr) as the rate of seafloor spreading increases between Africa and South America (Fig. 6b). While it has been shown that mantle plumes associated with the LLSVP have eroded a significant proportion of cratonic lithosphere in Africa 55 , and may be associated with some kimberlite eruptions 4 , this does not explain the subduction signatures in kimberlites or the timing of their eruption. We propose that a reservoir of recycled slabs occupy the mantle from pervasive subduction during the assembly of Pangea, which results in dehydration melting of the overlying mantle as water entrained with cold slabs is released 56 . Then, as Pangea begins to break-up, the rapid subduction of slab material at a low angle drives mantle return flow from this fertile mantle reservoir to drive kimberlite eruptions. Since subducting slabs influence the deep mantle structure, potentially triggering enhanced plume flux at the edges of the African LLSVP 7 , this could accelerate the delivery of slab-influenced melt from the uppermost lower mantle to the surface corresponding to an enhanced frequency of kimberlite eruptions at 120 Ma. This may be similar to a process contributing to the formation of the 260 Ma Emeishan large igneous province, where recycled Palaeo-Tethys Ocean beneath SW China has been proposed to induce large-scale mantle upwellings from 410-660 km depths 49 . The second pulse in kimberlite eruptions at 80-90 Ma cannot readily be linked to slab flux due to the vast distance from the nearest subduction zone following the opening of the South Atlantic Ocean (Fig. 7. Instead, an increase in African plate velocity would expose more cratonic lithosphere to mantle upwellings connected to the LLSVP, thereby increasing the kimberlite eruption frequency 3,4 (Fig. 6b).

North American kimberlites. A second population of kimberlite eruptions occurred between 110 and
40 Ma while North America migrated westward during the opening of the North Atlantic Ocean. It has been proposed that the dehydration of hydrous minerals stored within the flat-subducting Farallon plate promoted magmatism and kimberlite generation approximately 1500 km from the nearest trench 45 , however, geodynamic models suggest that flat subduction inhibits arc magmatism as the release and convection of fluids from the slab are obstructed by the asthenospheric wedge 57 . From our reconstructions of slab dip, the average dip angle along the western margin of North America varies between 30 and 36 • and the slab flux predicts the peaks and troughs in kimberlite eruption frequency between 110 and 40 Ma (Fig. 6c). Slab dip is spatially and temporally variable along North American subduction boundaries during the Laramide period, which has been attributed to the flat subduction of the Shatsky Rise conjugate on the northernmost section of the Farallon plate 40 . Its subduction predicts the distribution of magmatic and amagmatic zones in North America. From 95 to 60 Ma, the subduction of relatively young seafloor (5-50 Ma) combined with subduction of the buoyant conjugate Shatsky Rise leads to flat slab subduction beneath central USA 58 (Fig. 7). The distribution of kimberlite eruptions during this period are focused in Canada and the south of North America on either side of the conjugate (Fig. 8). Abrupt changes in subduction angles could be accommodated by slab tears adjacent to the Arizona-New-Mexico magmatic belt 57 . It is likely that melts associated with the dehydration of recycled slab material in the mantle transition zone were delivered to the surface through subduction-induced return flow 47 . Removal of the flat Farallon slab from the base of overriding continental lithosphere at 50 Ma 59 would further stimulate mantle return flow, triggering more widespread kimberlite eruptions which occur within the formerly amagmatic zone of central USA (Fig. 6c).
A unifying mechanism for kimberlite eruptions? The dichotomy of kimberlite eruption between African and North American populations is important because, while the mantle return flow mechanism is consistent, the stimulation of source regions in the mantle is different. In Africa, the mantle return flow associated with subducting slabs had likely stimulated upwellings along the edges of the LLSVP and invigorated a fertile mantle reservoir derived from sinking slab remnants left over from the assembly of Pangea, which explains the subduction signatures observed in this Cretaceous kimberlite population. The second pulse of African kimberlite eruptions at 80-90 Ma is likely connected to an increase in plate velocity as Africa migrates over upwellings associated with the LLSVP 3 . Meanwhile, kimberlite eruptions in North America are driven by upper mantle return flow in regions adjacent to the flat subduction of the Shatsky Rise and within the region affected by the Laramide Orogeny following the removal of the flat slab beneath North America. Importantly, both kimberlite populations are linked to the rapid subduction of the Farallon plate at low angle beneath the Americas, which suggests that this slab plays an important role in driving upwelling-induced volcanism. Whether the the Farallon slab is imbued with a high enrichment of volatiles, such as H 2 O which lowers the solidus temperature promoting partial melting and magma generation 45 , or if it stimulates preexisting fertile mantle reservoirs is unclear. Nevertheless, this study highlights the importance of subduction on the generation of kimberlites, and challenges previous conceptions that kimberlites are primarily generated by mantle plumes.

Conclusions
The dip angle of subducting oceanic lithosphere is a key parameter which characterises mantle and continental dynamics at subduction zones. We propose a simple framework for predicting slab dip from the thickness of the down-going slab, the convergence rate, the rate of trench migration, the density and volatile enrichment of the

Methods
The plate model used in this paper was modified from a recently-published model 15 as follows. Subduction zones to which no motions assigned were initially stationary through time; these were assigned new motions relative to the global plate circuit to model moderate trench retreat while still remaining consistent with tomographic constraints. The Orcas plate was split into two separate plates at 170-130 Ma, consistent with its configuration after 130 Ma, in order to accommodate divergence at the plate boundaries. The Caribbean plate was also split by a new back-arc spreading centre at 140-120 Ma, to allow for the formation of the Caribbean large igneous province at a spreading ridge 51 . Finally, the absolute plate motion model was constrained using an iterative optimisation workflow 24 . A ZIP archive containing plate reconstruction files for use in GPlates is available on Zenodo (doi.org/10.5281/zenodo.5769002).

Data availibility
The datasets generated and/or analysed during the current study are available in the Zenodo repository