Kimberlite eruptions driven by slab flux and subduction angle

Jul 12, 2023

Scientific Reports volume 13, Article number: 9216 (2023) Cite this article

55 Accesses

16 Altmetric

Metrics details

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.

Kimberlites are mafic volcanic rocks erupted from the Earth's mantle and are the host rocks of most diamonds1. Kimberlites occur on every craton and have sporadically been emplaced since 3 Ga2, yet the greatest number of kimberlite eruptions preserved on Earth today formed during the last 250 to 50 million years, primarily in Africa and North America3. While the distribution of kimberlites has been associated with the edges of large low-shear-wave-velocity provinces (LLSVPs)4, and changes in angular plate velocity3, these do not explain the frequency of kimberlite eruptions nor enriched radiogenic isotope signatures indicating a subducted slab component in some Cretaceous kimberlite populations1,5. The steep subduction of oceanic lithosphere into the mantle has been proposed to drive strong mantle return flow and pulses of magmatism6. Yet, despite the theoretical connection between volcanic eruptions and high rates of slab flux7, difficulties in estimating the volume and subduction angle of oceanic lithosphere being recycled at ancient subduction zones has thwarted any correlation with kimberlite eruptions. Previous attempts to characterise the dip angle of subducting slabs applied multivariate analysis of subduction zone characteristics to search for correlations among key parameters8,9,10,11,12,13,14. However, these approaches are mostly useful for reproducing the present-day slab dip and have limited application to subduction zones reconstructed through deep geological time. Here, using a recent tectonic plate reconstruction model15 and models of plate cooling16,17,18, we revisit the estimation of slab dip from simple plate kinematic parameters that characterise most subduction zones across the globe to explore the potential role of steeply subducting slabs in controlling kimberlite eruptions in Africa and North America.

Depths of subducting slabs obtained from the Slab2 model19 overlain with kimberlites erupted within the last 250 million years20. Trench segments have the following abbreviations in (i) Oceania: Ton, Tonga; Ker, Kermadec; NH, New Hebrides; Sol, Solomon; (ii) Southeast Asia: PNG, Papua New Guinea; Sum, Sumatra; Mar, Marianas; IZB, Izu-Bonin; Ryu, Ryukyu; Man, Manila; Ph, Philippine; (iii) Asia: Mak, Makran; SJ, South Japan; NJ, North Japan; Kur, Kurile; (iv) Europe: Hel, Helenic; Cal, Calabria; (v) North America: Al, Aleutians; Cas, Cascades; (vi) Central America: Mex, Mexico; MAM, Middle America; LAT, Lesser Antilles; (vii) South America: EC, Ecuador; SA, South America; SC, South Chile; SSW, South Sandwich. Approximate trench locations, boundary types, names, and corresponding abbreviations are listed in Table S1. White areas indicates regions of non-oceanic crust; thick red lines indicate mid-ocean ridges; thin red lines indicate transform boundaries. Map was generated using Cartopy21.

The downward pull of slabs into the Earth's mantle is the largest driving force in plate tectonics22, where the bending of the subducting plate plays an important role in modulating the amount of slab pull which is transferred to the surface to drive plate motion23. Since the 1980s8, multiple studies have explored correlations between subduction angle and a number of parameters, including the duration of subduction, convergence rate, nature of the over-riding plate, trench length and others8,9,10,11,12,13,14. A recent study returned to a multivariate regression and reinforced the long-held view that slab dip is influenced by subduction duration and, to a lesser extent, the age of the down-going plate and whether the overriding plate is continental or oceanic in nature13. While these parameters are useful for approximating the present-day slab dip, the demarcation of subduction zone segments are subjective and time-dependent, which poses a difficulty when attempting to cast these slab dip relationships back through deep geological time with constantly evolving subduction boundaries. For example, data on the duration or 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 rotation24. 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 model19, 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 model25 using pyGPlates26. The dip angle, \(\theta\), of the down-going slab is taken at multiple depth intervals inboard from the trench (Fig. 2). The average dip angle, \(\theta _{\textrm{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/brmather/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.

Schematic representation of slab dips and subduction parameters. The convergence rate, \(v_c\), is the summation of the velocity of the subducting plate, \(v_s\), and the overriding plate velocity, \(v_o\). We assume the velocity of the overriding plate and the trench are equal and opposite (\(v_o = -v_t\)). \(v_s\) and \(v_o\) are positive towards the trench and each vector is orthogonal to the trench. \(v_{\textrm{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 \(\theta\) is the dip angle of the subducting slabs calculated at different intervals.

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 lithosphere8,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_{\textrm{plate}}\) (\(P=0.59\)) (Fig. 3). Plate thickness was predicted from the thermal evolution of oceanic lithosphere,

where \(\sqrt{\kappa t}\) derives from models of plate cooling16,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 \(\kappa = 1\) mm\(^2\)/s. We employ a plate model of oceanic lithosphere cooling, 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 Earth16. The product of the plate thickness, \(h_{\textrm{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.

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.

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 boundary12,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 zones10. 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.

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 plate30. The source of volatile enrichment occurs at mid-ocean ridges, where volatile-bearing melt circulates through channels in newly-forming oceanic lithosphere31. The sequestration rate of volatiles within oceanic lithosphere, \(q_v\), is proportional to the seafloor spreading rate, \(v_{\textrm{hsp}}\)31. We sample seafloor spreading rate grids, generated using workflows in32, at trench boundaries in the same way that we interpolate age grids to calculate plate thickness, \(h_{\textrm{plate}}\). We find that \(v_{\textrm{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 crust33, sepertinisation during ultraslow spreading34, and during bending and cracking of the plate before entering the subduction zone35, however, the spreading rate relationship we outline above constitutes the most general source of volatiles subducted at most trenches globally.

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 depletion36. The crust is the main source of positive buoyancy in oceanic lithosphere due to its relatively low density (\(\sim\) 2900 kg/m\(^3\))37. Assuming modern mantle temperature conditions, 10–20 Myr old lithosphere is negatively buoyant36 and becomes more negatively buoyant with age. A 60 Myr segment of oceanic lithosphere is approximately 79.4 km thick16, 7 km of which is crust, which equates to a mean density of \(\rho _{\textrm{av}} = 3278\) kg/m\(^3\). While the buoyancy of oceanic plates may not be sufficient to initiate subduction alone38, 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 zones39 or lead to flat slab subduction40.

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, \(\theta _{\textrm{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_{\textrm{train}}\) to the test dataset \(X_{\textrm{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, \(d_k = \Vert X_{\textrm{train}} - X_{\textrm{test}} \Vert _k\). 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.

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.

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/brmather/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 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 stagnation41,42 and potentially slab anchoring43. 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.

Kimberlites are volcanic rocks which ascend rapidly from the mantle and are emplaced in cratons worldwide2. Kimberlite eruptions have been associated with mantle upwellings from the LLSVP4 and extensional tectonics associated with lithospheric unloading44 or changes in plate velocity3. 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 ratios1,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 Orogeny45. Here, magma may have been generated from water-fluxed decompression melting of the mantle transition zone46 and transported upwards by subduction-induced mantle return flow47. 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 mantle48. Subduction-driven mantle upwellings have been linked to the formation of the 260 Ma Emeishan large igneous province in SW China49 and Cenozoic volcanism in NE China50. 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 pyGPlates26. We used a 170 Myr plate model, modified from15, with significantly improved subduction zone boundaries along the western margin of North America, and improved resolution of the Caribbean plate51 (refer to Methods). The locations of kimberlites \(\le\) 170 Ma were reconstructed back to their time of eruption, using a compilation of kimberlites with eruption ages52 that we resampled to a 3-times refined icosohedral mesh53 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\(^\circ\)) predicted from our slab dip analysis on reconstructed subduction 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 flow35, 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.

Spatio-temporal association between kimberlite eruptions and LLSVP boundaries. (a) Shear wave velocity anomaly from the SMEAN2 tomography model54 overlain with the spatial distribution of kimberlites (diamonds) reconstructed to their eruption age. We define LLSVP boundaries as the 1% slow contour of the 2800 km depth slice from the velocity variation which is represented by the red dashed line. Present-day coastlines are added for reference. Map was generated using Cartopy21. (b) Cumulative density function of the distance between LLSVP boundaries and kimberlite populations in Africa and North America with random continental locations at the present day and reconstructed to 170 Ma when the Pangea supercontinent was assembled. African kimberlites hold a statistically significant relationship to LLSVP boundaries while North American kimberlites do not.

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 \(^\circ\) 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 Africa55, and may be associated with some kimberlite eruptions4, 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 released56. 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 LLSVP7, 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 depths49. 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 frequency3,4 (Fig. 6b).

Relationship between kimberlite eruption density, slab flux, and plate velocity over 170 million years. (a) Pulses in combined kimberlite eruptions correlate with periods of high slab flux of the east-dipping Farallon plate along the western margin of North America. The low dip angle of subduction (30–35\(^\circ\)) predicted from our analysis is indicated by a strong sideways component of slab flux. (b) the first peak in African kimberlite eruptions is correlated with high slab flux between 120 and 130 Ma as Pangea dispersed; the second peak at 80–90 Ma is more likely to be explained by an increase in African plate velocity3 than slab flux due to the the increased distance between the Americas and Africa from the opening of the South Atlantic Ocean. (c) Pulses in North American kimberlite eruptions is closely correlated to pulses in slab flux integrated along east-dipping subduction zones within 2500 km radius of the nearest kimberlite eruption in North America with no correlation observed with plate velocity.

Evolution of slab dip angle at global subduction zones from 160 Ma to the present day overlain with the spatial distribution of kimberlite eruptions. A large peak in slab flux and kimberlite eruption density occurs at 120 Ma during the breakup of the Gondwana supercontinent, and a second (smaller) peak at 80 Ma associated with the removal of the flat Farallon slab from the base of overriding continental lithosphere in North America. White regions indicate non-oceanic crust, grey regions indicate present-day coastlines. Maps were generated using Cartopy21. A full timeseries from 170 Ma to 0 Ma is available as an animation in Movie S1.

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 trench45, 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 wedge57. From our reconstructions of slab dip, the average dip angle along the western margin of North America varies between 30 and 36\(^\circ\) 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 plate40. 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 USA58 (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 belt57. 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 flow47. Removal of the flat Farallon slab from the base of overriding continental lithosphere at 50 Ma2.3.CO;2 (1995)." href="/articles/s41598-023-36250-w#ref-CR59" id="ref-link-section-d64025640e2584">59 would further stimulate mantle return flow, triggering more widespread kimberlite eruptions which occur within the formerly amagmatic zone of central USA (Fig. 6c).

Evolution of flat slab subduction along North America. (a) Flat slab subduction of the Farallon Plate is caused by the Shatsky Rise which enters the trench at 95 Ma. (b) Much of the Shatsky conjugate is consumed by 70 Ma by which time the Farallon underplates much of central USA, producing tectonic uplift associated with the Laramide Orogeny. (c) At 50 Ma the Farallon plate is removed from the base of overriding continental lithosphere resulting in widespread kimberlite eruptions. Black polygons indicate reconstructed large igneous provinces, grey polygons indicate inferred LIP conjugates, arrows indicate absolute plate velocity. Maps were generated using Cartopy21.

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 LLSVP3. 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 generation45, 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.

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 slab. Applying this framework to plate reconstructions provides new insights into the dynamics of past subduction zones, the spatial distribution of arc volcanism through deep geological time, and the fate of subducted slabs. Using this predictive framework, we reconstruct the slab dip angle of subduction zone segments over the last 170 million years to help explain pulses in kimberlite eruptions. High subduction rates stimulate mantle return flow which promotes partial melting and magma generation. Kimberlite eruptions in Africa and North America are linked to the subduction of the Farallon plate beneath the Americas. In Africa, peaks in kimberlite eruptions exhibit a strong correlation with high slab flux during the initial stages of supercontinent breakup and high plate velocities (up to 6 cm/yr) as Africa migrates over the LLSVP. In North America, the subduction of the Shatsky Rise conjugate from 95–50 Ma results in flat subduction beneath central USA, associated with the Laramide Orogeny, which limits magmatism to the edges of the flat slab until its removal from the base of the lithosphere after 50 Ma. Our results highlight the important role of subduction angle in modulating horizontal slab flux and thereby the distribution and timing of volcanism. This helps to explain the dichotomy of kimberlite populations in Africa and North America and has important implications for the formation of ancient kimberlites in Australia, India, and South America which may be explained by reconstructions of slab dip through deep geological time.

The plate model used in this paper was modified from a recently-published model15 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 ridge51. Finally, the absolute plate motion model was constrained using an iterative optimisation workflow24. A ZIP archive containing plate reconstruction files for use in GPlates is available on Zenodo (doi.org/10.5281/zenodo.5769002).

The datasets generated and/or analysed during the current study are available in the Zenodo repository, https://doi.org/10.5281/zenodo.5831990. A series of Jupyter notebooks containing a Python workflow to calculate slab dip using a plate reconstruction model is also available at Zenodo via https://doi.org/10.5281/zenodo.5831990 and GitHub (https://github.com/brmather/Slab-Dip). The plate reconstruction software, GPlates and pyGPlates, are freely available from www.gplates.org/download.

Woodhead, J. et al. Kimberlites reveal 2.5-billion-year evolution of a deep, isolated mantle reservoir. Nature 573, 578–581. https://doi.org/10.1038/s41586-019-1574-8 (2019).

Article ADS CAS PubMed Google Scholar

Shirey, S. B. & Richardson, S. H. Start of the Wilson cycle at 3 Ga shown by diamonds from subcontinental mantle. Science 333, 434–436. https://doi.org/10.1126/science.1206275 (2011).

Article ADS CAS PubMed Google Scholar

Tappe, S., Smart, K., Torsvik, T., Massuyeau, M. & de Wit, M. Geodynamics of kimberlites on a cooling Earth: Clues to plate tectonic evolution and deep volatile cycles. Earth Planet. Sci. Lett. 484, 1–14. https://doi.org/10.1016/j.epsl.2017.12.013 (2018).

Article ADS CAS Google Scholar

Torsvik, T. H., Burke, K., Steinberger, B., Webb, S. J. & Ashwal, L. D. Diamonds sampled by plumes from the core-mantle boundary. Nature 466, 352–355. https://doi.org/10.1038/nature09216 (2010).

Article ADS CAS PubMed Google Scholar

Tappe, S., Graham Pearson, D., Kjarsgaard, B. A., Nowell, G. & Dowall, D. Mantle transition zone input to kimberlite magmatism near a subduction zone: origin of anomalous Nd-Hf isotope systematics at Lac de Gras, Canada. Earth Planet. Sci. Lett. 371–372, 235–251. https://doi.org/10.1016/j.epsl.2013.03.039 (2013).

Article ADS CAS Google Scholar

Faccenna, C. et al. Subduction-triggered magmatic pulses: a new class of plumes?. Earth Planet. Sci. Lett. 299, 54–68. https://doi.org/10.1016/j.epsl.2010.08.012 (2010).

Article ADS CAS Google Scholar

Cao, X., Flament, N., Bodur, Ö. F. & Müller, R. D. The evolution of basal mantle structure in response to supercontinent aggregation and dispersal. Sci. Rep. 11, 1–16. https://doi.org/10.1038/s41598-021-02359-z (2021).

Article CAS Google Scholar

Jarrard, R. D. Relations among subduction parameters. Rev. Geophys. 24, 217–284. https://doi.org/10.1029/RG024i002p00217 (1986).

Article ADS Google Scholar

Lallemand, S., Heuret, A. & Boutelier, D. On the relationships between slab dip, back-arc stress, upper plate absolute motion, and crustal nature in subduction zones. Geochem. Geophys. Geosyst.https://doi.org/10.1029/2005GC000917 (2005).

Article Google Scholar

Syracuse, E. M. & Abers, G. A. Global compilation of variations in slab depth beneath arc volcanoes and implications. Geochem. Geophys. Geosyst.https://doi.org/10.1029/2005GC001045 (2006).

Article Google Scholar

Schellart, W. P., Freeman, J., Stegman, D. R., Moresi, L. & May, D. Evolution and diversity of subduction zones controlled by slab width. Nature (London) 446, 308–311 (2007).

Article ADS CAS PubMed Google Scholar

MacDougall, J. G., Kincaid, C., Szwaja, S. & Fischer, K. M. The impact of slab dip variations, gaps and rollback on mantle wedge flow: Insights from fluids experiments. Geophys. J. Int. 197, 705–730. https://doi.org/10.1093/gji/ggu053 (2014).

Article ADS Google Scholar

Hu, J. & Gurnis, M. Subduction duration and slab dip. Geochem. Geophys. Geosyst. 21, 1–34. https://doi.org/10.1029/2019GC008862 (2020).

Article Google Scholar

Schellart, W. P. Control of subduction zone age and size on flat slab subduction. Front. Earth Sci. 8, 1–18. https://doi.org/10.3389/feart.2020.00026 (2020).

Article Google Scholar

Clennett, E. J. et al. A quantitative tomotectonic plate reconstruction of Western North America and the Eastern Pacific basin. Geochem. Geophys. Geosyst. 21, 1–25. https://doi.org/10.1029/2020GC009117 (2020).

Article Google Scholar

Grose, C. J. Properties of oceanic lithosphere: Revised plate cooling model predictions. Earth Planet. Sci. Lett. 333–334, 250–264. https://doi.org/10.1016/j.epsl.2012.03.037 (2012).

Article ADS CAS Google Scholar

Richards, F. D., Hoggard, M. J., Cowton, L. R. & White, N. J. Reassessing the thermal structure of oceanic lithosphere with revised global inventories of basement depths and heat flow measurements. J. Geophys. Res. Solid Earth 123, 9136–9161. https://doi.org/10.1029/2018JB015998 (2018).

Article ADS Google Scholar

Richards, F., Hoggard, M., Crosby, A., Ghelichkhan, S. & White, N. Structure and dynamics of the oceanic lithosphere-asthenosphere system. Phys. Earth Planet. Inter.https://doi.org/10.1016/j.pepi.2020.106559 (2020).

Article Google Scholar

Hayes, G. P. et al. Slab2, a comprehensive subduction zone geometry model. Science 362, 58–61. https://doi.org/10.1126/science.aat4723 (2018).

Article ADS CAS PubMed Google Scholar

Giuliani, A. & Pearson, D. G. Kimberlites: From deep earth to diamond mines. Elements 15, 377–380. https://doi.org/10.2138/GSELEMENTS.15.6.377 (2019).

Article CAS Google Scholar

Met Office. Cartopy: a cartographic python library with a Matplotlib interface. Exeter, Devon, https://doi.org/10.5281/zenodo.7430317 (2023).

Conrad, C. P. & Lithgow-Bertelloni, C. The temporal evolution of plate driving forces: Importance of slab suction versus slab pull during the Cenozoic. J. Geophys. Res. Solid Earth 109, 1–14. https://doi.org/10.1029/2004JB002991 (2004).

Article Google Scholar

Capitanio, F. A., Morra, G. & Goes, S. Dynamics of plate bending at the trench and slab-plate coupling. Geochem. Geophys. Geosyst.https://doi.org/10.1029/2008GC002348 (2009).

Article Google Scholar

Tetley, M. G., Williams, S. E., Gurnis, M., Flament, N. & Müller, R. D. Constraining absolute plate motions since the Triassic. J. Geophys. Res. Solid Earth 124, 7231–7258. https://doi.org/10.1029/2019JB017442 (2019).

Article Google Scholar

Müller, R. D. et al. A global plate model including lithospheric deformation along major rifts and orogens since the Triassic. Tectonics 38, 1884–1907. https://doi.org/10.1029/2018TC005462 (2019).

Article ADS Google Scholar

Müller, R. D. et al. GPlates: Building a virtual earth through deep time. Geochem. Geophys. Geosyst. 19, 2243–2261. https://doi.org/10.1029/2018GC007584 (2018).

Article ADS Google Scholar

Parsons, B. & Sclater, J. G. An analysis of the variation of ocean floor bathymetry and heat flow with age. J. Geophys. Res. 82, 803–827. https://doi.org/10.1029/JB082i005p00803 (1977).

Article ADS Google Scholar

Capitanio, F. A., Stegman, D. R., Moresi, L. N. & Sharples, W. Upper plate controls on deep subduction, trench migrations and deformations at convergent margins. Tectonophysics 483, 80–92. https://doi.org/10.1016/j.tecto.2009.08.020 (2010).

Article ADS Google Scholar

Schepers, G. et al. South-American plate advance and forced Andean trench retreat as drivers for transient flat subduction episodes. Nat. Commun. 8, 1–9. https://doi.org/10.1038/ncomms15249 (2017).

Article CAS Google Scholar

Gonzalez, C. M., Gorczyk, W. & Gerya, T. V. Decarbonation of subducting slabs: Insight from petrological–thermomechanical modeling. Gondwana Res. 36, 314–332. https://doi.org/10.1016/j.gr.2015.07.011 (2016).

Article ADS CAS Google Scholar

Keller, T., Katz, R. F. & Hirschmann, M. M. Volatiles beneath mid-ocean ridges: Deep melting, channelised transport, focusing, and metasomatism. Earth Planet. Sci. Lett. 464, 55–68. https://doi.org/10.1016/j.epsl.2017.02.006 (2017).

Article ADS CAS Google Scholar

Seton, M. et al. A global data set of present-day oceanic crustal age and seafloor spreading parameters. Geochem. Geophys. Geosyst. 21, e2020GC009214. https://doi.org/10.1029/2020GC009214 (2020).

Article ADS CAS Google Scholar

Dutkiewicz, A., Müller, R. D., Cannon, J., Vaughan, S. & Zahirovic, S. Sequestration and subduction of deep-sea carbonate in the global ocean since the early cretaceous. Geology 47, 91–94. https://doi.org/10.1130/G45424.1 (2019).

Article ADS CAS Google Scholar

Merdith, A. S., Atkins, S. E. & Tetley, M. G. Tectonic controls on carbon and serpentinite storage in subducted upper oceanic lithosphere for the past 320 Ma. Front. Earth Sci. 7, 1–23. https://doi.org/10.3389/feart.2019.00332 (2019).

Article ADS Google Scholar

Faccenda, M., Gerya, T. V. & Burlini, L. Deep slab hydration induced by bending-related variations in tectonic pressure. Nat. Geosci. 2, 790–793. https://doi.org/10.1038/ngeo656 (2009).

Article ADS CAS Google Scholar

Korenaga, J. Initiation and evolution of plate tectonics on earth: Theories and observations. Ann. Rev. Earth Planet. Sci. 41, 117–151. https://doi.org/10.1146/annurev-earth-050212-124208 (2013).

Article ADS CAS Google Scholar

Afonso, J. C., Ranalli, G. & Fernàndez, M. Density structure and buoyancy of the oceanic lithosphere revisited. Geophys. Res. Lett.https://doi.org/10.1029/2007GL029515 (2007).

Article Google Scholar

Gurnis, M., Hall, C. & Lavier, L. Evolving force balance during incipient subduction. Geochem. Geophys. Geosyst.https://doi.org/10.1029/2003GC000681 (2004).

Article Google Scholar

Knesel, K. M., Cohen, B. E., Vasconcelos, P. M. & Thiede, D. S. Rapid change in drift of the Australian plate records collision with Ontong Java plateau. Nature 454, 754–757. https://doi.org/10.1038/nature07138 (2008).

Article ADS CAS PubMed Google Scholar

Liu, L. et al. The role of oceanic plateau subduction in the Laramide orogeny. Nat. Geosci. 3, 353–357. https://doi.org/10.1038/ngeo829 (2010).

Article ADS CAS Google Scholar

Mao, W. & Zhong, S. Slab stagnation due to a reduced viscosity layer beneath the mantle transition zone. Nat. Geosci. 11, 876–881. https://doi.org/10.1038/s41561-018-0225-2 (2018).

Article ADS CAS Google Scholar

King, S. D., Frost, D. J. & Rubie, D. C. Why cold slabs stagnate in the transition zone. Geology 43, 231–234. https://doi.org/10.1130/G36320.1 (2015).

Article ADS Google Scholar

Chen, Y. W., Wu, J. & Suppe, J. Southward propagation of Nazca subduction along the Andes. Nature 565, 441–447. https://doi.org/10.1038/s41586-018-0860-1 (2019).

Article ADS CAS PubMed Google Scholar

Jelsma, H., Barnett, W., Richards, S. & Lister, G. Tectonic setting of kimberlites. Lithos 112, 155–165. https://doi.org/10.1016/j.lithos.2009.06.030 (2009).

Article ADS CAS Google Scholar

Currie, C. A. & Beaumont, C. Are diamond-bearing Cretaceous kimberlites related to low-angle subduction beneath western North America?. Earth Planet. Sci. Lett. 303, 59–70. https://doi.org/10.1016/j.epsl.2010.12.036 (2011).

Article ADS CAS Google Scholar

Kjarsgaard, B. A., Heaman, L. M., Sarkar, C. & Pearson, D. G. The North America mid-Cretaceous kimberlite corridor: Wet, edge-driven decompression melting of an OIB-type deep mantle source. Geochem. Geophys. Geosyst. 18, 2727–2747. https://doi.org/10.1002/2016GC006761 (2017).

Article ADS CAS Google Scholar

Chen, Y. et al. Reconciling seismic structures and Late Cretaceous kimberlite magmatism in northern Alberta, Canada. Geology 48, 872–876. https://doi.org/10.1130/G47163.1 (2020).

Article ADS CAS Google Scholar

Mather, B. R. et al. Intraplate volcanism triggered by bursts in slab flux. Sci. Adv. 6, 1–8. https://doi.org/10.1126/sciadv.abd0953 (2020).

Article Google Scholar

He, C. Imprints of subducted Palaeo-Tethys oceanic lithosphere on upper-mantle discontinuities and the formation of the Emeishan large igneous province. Geophys. J. Int. 231, 1298–1308. https://doi.org/10.1093/gji/ggac251 (2022).

Article ADS Google Scholar

Yang, J. & Faccenda, M. Intraplate volcanism originating from upwelling hydrous mantle transition zone. Nature 579, 88–91. https://doi.org/10.1038/s41586-020-2045-y (2020).

Article ADS CAS PubMed Google Scholar

García-Reyes, A. & Dyment, J. Structure, age, and origin of the Caribbean Plate unraveled. Earth Planet. Sci. Lett. 571, 117100. https://doi.org/10.1016/j.epsl.2021.117100 (2021).

Article CAS Google Scholar

Giuliani, A. et al. Kimberlite genesis from a common carbonate-rich primary melt modified by lithospheric mantle assimilation. Sci. Adv. 6, 1–10. https://doi.org/10.1126/sciadv.aaz0424 (2020).

Article CAS Google Scholar

Moresi, L. & Mather, B. Stripy. A Python module for (constrained) triangulation in Cartesian coordinates and on a sphere. Journal of Open Source Software 4, 1410. https://doi.org/10.21105/joss.01410 (2019).

Article ADS Google Scholar

Jackson, M. G., Konter, J. G. & Becker, T. W. Primordial helium entrained by the hottest mantle plumes. Nature 542, 340–343. https://doi.org/10.1038/nature21023 (2017).

Article ADS CAS PubMed Google Scholar

Celli, N. L., Lebedev, S., Schaeffer, A. J. & Gaina, C. African cratonic lithosphere carved by mantle plumes. Nat. Commun.https://doi.org/10.1038/s41467-019-13871-2 (2020).

Article PubMed PubMed Central Google Scholar

Regier, M. E. et al. The lithospheric-to-lower-mantle carbon cycle recorded in superdeep diamonds. Nature 585, 234–238. https://doi.org/10.1038/s41586-020-2676-z (2020).

Article CAS PubMed Google Scholar

Axen, G. J., van Wijk, J. W. & Currie, C. A. Basal continental mantle lithosphere displaced by flat-slab subduction. Nat. Geosci. 11, 961–964. https://doi.org/10.1038/s41561-018-0263-9 (2018).

Article ADS CAS Google Scholar

Currie, C. A. & Copeland, P. Numerical models of Farallon plate subduction: Creating and removing a flat slab. Geosphere 18, 476–502. https://doi.org/10.1130/GES02393.1 (2022).

Article Google Scholar

Humphreys, E. D. Post-Laramide removal of the Farallon slab, western United States. Geology 23, 987. 2.3.CO;2">https://doi.org/10.1130/0091-7613(1995)023<0987:PLROTF>2.3.CO;2 (1995).

2.3.CO;2" data-track-action="article reference" href="https://doi.org/10.1130%2F0091-7613%281995%29023%3C0987%3APLROTF%3E2.3.CO%3B2" aria-label="Article reference 59" data-doi="10.1130/0091-7613(1995)0232.3.CO;2">Article ADS Google Scholar

Download references

This study was supported by AuScope Simulation, Analysis & Modelling node funded by the Australian Government through the National Collaborative Research Infrastructure Strategy, NCRIS. We gratefully acknowledge funding from the Australian Research Council through grant DP200100966 (M.S.).

EarthByte Group, School of Geosciences, The University of Sydney, Sydney, 2006, Australia

Ben R. Mather, R. Dietmar Müller, Christopher P. Alfonso, Maria Seton & Nicky M. Wright

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

B.M. conceived the idea, developed the Python tools, generated the figures and wrote the manuscript. R.D.M. conceived the idea and applications of the slab dip formulation to explain kimberlite eruptions. C.P.A. improved the tectonic reconstruction in North America and provided age grids and seafloor spreading rate grids. M.S. conceived the experiments and improved the manuscript. N.M.W. provided the reconstructions of large igneous provinces.

Correspondence to Ben R. Mather.

The authors declare no competing interests

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information 1.

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://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Mather, B.R., Müller, R.D., Alfonso, C.P. et al. Kimberlite eruptions driven by slab flux and subduction angle. Sci Rep 13, 9216 (2023). https://doi.org/10.1038/s41598-023-36250-w

Download citation

Received: 25 October 2022

Accepted: 30 May 2023

Published: 06 June 2023

DOI: https://doi.org/10.1038/s41598-023-36250-w

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.