Mapping global kimberlite potential from reconstructions of mantle flow over the past billion years.

Kimberlites are the primary source of economic grade diamonds. Their geologically rapid eruptions preferentially occur near or through thick and ancient continental lithosphere. Studies combining tomographic models with tectonic reconstructions and kimberlite emplacement ages and locations have revealed spatial correlations between large low shear velocity provinces in the lowermost mantle and reconstructed global kimberlite eruption locations over the last 320 Myr. These spatial correlations assume that the lowermost mantle structure has not changed over time, which is at odds with mantle flow models that show basal thermochemical structures to be mobile features shaped by cold sinking oceanic lithosphere. Here we investigate the match to the global kimberlite record of stationary seismically slow basal mantle structures (as imaged through tomographic modelling) and mobile hot basal structures (as predicted by reconstructions of mantle flow over the past billion years). We refer to these structures as "basal mantle structures" and consider their intersection with reconstructed thick or ancient lithosphere to represent areas with a high potential for past eruptions of kimberlites, and therefore areas of potential interest for diamond exploration. We use the distance between reconstructed kimberlite eruption locations and kimberlite potential maps as an indicator of model success, and we find that mobile lowermost mantle structures are as close to reconstructed kimberlites as stationary ones. Additionally, we find that mobile lowermost mantle structures better fit major kimberlitic events, such as the South African kimberlite bloom around 100 Ma. Mobile basal structures are therefore consistent with both solid Earth dynamics and with the kimberlite record.

1. Introduction

The majority of diamonds, including economic grade diamonds, were formed at the base of thick continental lithosphere [1]. Using xenolith data from kimberlitic eruptions, ref. [2] determined the temperature-pressure window that allows for diamond formation and stability as between 1200°C and 1570°C and above 5 GPa. “Clifford’s rule” postulates that cratonic diamonds are formed in cratons with thick lithospheric keels, providing environments with sufficiently high pressures and low temperatures for diamond stability [3, 4].

Cratons of Archean age are particularly thick [possibly up to ~300 km; e.g., 5] and are an ideal environment for diamond formation [6-8]. Kimberlite eruptions are rapid magmatic events with magma ascending at speeds up to 20 m s-1 [9, 10], which is geologically instantaneous compared to average plate motion at ~1.26 × 10−9 m s-1 (4 cm yr-1) over the last 200 Myr [11]. Kimberlites originate from depths in excess of 120–150 km [12], forming from melts that may pool from up to ~300 km depth [13]. Some rare superdeep (sub-lithospheric) kimberlites may have ascended from as deep as 800 km [14-17], although these diamonds could have been transported to the base of the lithosphere before being entrained by kimberlite magmas [18]. The kimberlite record extends deep in geological times [the oldest known kimberlite pipe is 2.85 Gyr old; 19] and is considered to be episodic from 2 Ga [20-22]. Over 60% of known kimberlite pipes erupted during Mesozoic and Cenozoic times [22], which could be either the result of the preferential preservation of younger kimberlites or of increased eruption rate during the last 200 Myr [23]. Southern Africa underwent intensive kimberlite activity during this period, and this area preserves the most known kimberlite occurrences [12, 21, 22].

Combining databases of volcanic products with tectonic reconstructions and SMEAN [an average of three S-wave tomographic models; 24], ref. [25] showed that the reconstructed locations of the majority of large igneous provinces (LIPs) and kimberlites from 320 Ma (to the exception of Canadian kimberlites) are within ~1,000 km of the edges of two large seismically slow provinces above the core-mantle boundary (CMB) commonly known as large low-shear velocity provinces [LLSVPs; 26]. The LLSVPs, presently positioned under the Pacific Ocean (the Pacific LLSVP) and under the African continent (the African LLSVP), are ~7,500 km in radius and 500–1,000 km in height test–together they cover up to 50% of the CMB [26, 27]. The low shear-wave velocity of LLSVPs indicates that these structures are hotter than surrounding ambient mantle, and gradients in seismic waves at the edges of LLSVPs [28] suggest that the structures are intrinsically denser than surrounding mantle, although the density contrast remains debated [e.g., 29].

Due to their spatial relationship with volcanic surface features, the regions along the edge of LLSVPs are sometimes called plume generation zones (PGZs) [30]. The specific relationship between the edges of LLSVPs and reconstructed locations of volcanic products is debated, with studies showing that the relationship is instead between the interior of LLSVPs and reconstructed locations of volcanic products [31, 32]. It is generally accepted that major LIPs are linked to deep mantle plumes [33, 34], and the synchronicity of Mesozoic-Cenozoic kimberlite blooms with periods of LIP eruption suggests a potentially common heat source for kimberlite and LIP eruptions [25, 35]. However, the physical process connecting kimberlite magmatism [from 120–300 km depth; 12, 13] to LLSVPs [between 2000 km and 2500 km depth; 26] remains to be established. The spatial link between LLSVPs and reconstructed kimberlite eruption locations implies that hot basal mantle structures could have been stationary and rigid over time [25]. In contrast, global geodynamic models have revealed that similar hot basal mantle structures can form in response to the history of surface plate motions and plate subduction [36, 37].

Here we assume that kimberlite eruptions are linked to mantle upwelling occurring above basal mantle structures, following ref. [25]. We use the global tectonic reconstruction of ref. [38] to derive boundary conditions for forward mantle convection models that predict the evolution of the structure of the mantle over the past billion years. We map the intersection of basal mantle structures (either stationary from tomography or mobile from flow models) and reconstructed lithosphere, which represent areas with a high potential for past eruptions of kimberlites and are therefore of potential interest to diamond exploration. We use the distance between the reconstructed locations of known kimberlite eruptions and model high-potential areas as an indicator of model success and evaluate the fit of different models to the kimberlite record for selected periods and regions of high kimberlite activity.

2. Datasets and methods

Our approach uses tectonic plate reconstructions and kimberlite emplacement ages to link past kimberlite eruption locations, thick or ancient lithosphere, and basal mantle structures either imaged by tomographic models or predicted by mantle flow models. Key model results and scripts are available at https://doi.org/10.5281/zenodo.5760115.

2.1. Tectonic reconstructions

We use the global tectonic reconstruction by ref. [38], which is the first continuous full-plate reconstruction for the last billion years. This tectonic reconstruction links full-plate reconstructions by ref. [39] from 1,000–520 Ma, by ref. [40, 41] from 500 Ma to 410 Ma and by ref. [42] from 410 to 0 Ma. We consider the tectonic reconstruction in its original paleomagnetic reference frame (reconstruction hereafter referred to as M21), as well as in a reference frame in which the net rotation of the lithosphere with respect to the mantle was removed, which is more appropriate for mantle flow models [e.g., 43, 44]. This latter reconstruction is hereafter referred to as M21NNR (“no-net-rotation”).

2.2. Kimberlite locations and emplacement ages, and thick or ancient lithosphere

We use the kimberlite database of ref. [22] that contains 1,133 kimberlites with a peak in occurrences at around 100 Ma (Fig 1B)–notably in southern Africa (Fig 1A). Kimberlite occurrences were sampled at 20 Myr intervals to match the temporal resolution used for the output of global mantle flow models (see below). For each age of interest a, kimberlite occurrences were selected if they erupted within the period a + 10 Myr < a ≤ a—10 Myr.

Fig 1

Distribution of thick or ancient lithosphere and kimberlite eruptions.

A/ 957 kimberlite eruptions from ref. [22] from 640 Ma to the present, sampled in 20 Myr increments and coloured by age, and outlines of thick or ancient lithosphere considered in this study: magenta polygons show lithosphere thicker than 150 km in all four lithospheric thickness models considered in ref. [5], beige polygons are tectonic blocks from ref. [38] and grey polygons are blocks of tectonothermal age greater than 2.5 Ga from ref. [46]. Coastlines are shown in grey. This figure was created with GMT6 [47] and the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastlines are republished from [48] under a CC BY license, with permission from Paul Wessel, original copyright 1996. B/ Histogram showing the temporal distributions of kimberlites occurrences from ref. [22] from 640 Ma. None of the data in this figure are proprietary.

The age of the lithosphere can be established locally based on xenoliths carried to the surface by kimberlites [e.g., 45], and continental scale maps require interpolation. We therefore consider three end-member sets of thick or ancient lithosphere: 1/ lithosphere thicker than 150 km (a proposed proxy for cratons, ref. [4]) in all four lithospheric thickness models from ref. [5]; 2/ regions with tectonothermal ages greater than 2.5 Ga in the TC1 model of ref. [46] and 3/ tectonic blocks inferred to have existed for at least one billion years from ref. [38]. Almost all present-day locations of kimberlites fall within the tectonic blocks from ref. [38], however, this is not the case for lithosphere thicker than 150 km and for Archean lithosphere from ref. [46] that cover smaller areas (Fig 1A).

2.3. Reconstructions of past global mantle flow

We reconstruct the evolution of mantle flow from one billion years ago to the present using CitcomS [49] in which the mantle is a shell represented with finite-elements in spherical geometry. We use 129 × 129 × 65 × 12 ≈ 13 million elements to obtain an average resolution of ~50 km × 50 km × 15 km at the surface, ~40 km × 40 km × 100 km in the mid-mantle, and ~28 km × 28 km × 27 km at the core-mantle boundary (CMB). We use a version of CitcomS [50] modified to progressively assimilate the thermal structure of the lithosphere and of subducting slabs (to ~350 km depth) determined from the synthetic age of the ocean floor [51] as well as surface velocities obtained from tectonic reconstructions with continuously closing plate polygons [52]. These boundary conditions are assimilated in one-million-year intervals. This approach makes it possible to obtain one-sided subduction in time-dependent global mantle convection models with computationally affordable resolution and viscosity variations.

We use the extended-Boussinesq approximation, which accounts for viscous dissipation, an adiabatic temperature gradient, internal heating and a decrease in the coefficient of thermal expansion by a factor of two over the thickness of the mantle [53]. Convection vigour was determined by the Rayleigh number , where the subscript “0” indicates reference values, with α0 the coefficient of thermal expansion, ρ0 = 4,000 kg m-3 is the density, g0 = 9.81 m s-2 the gravitational acceleration, ΔT = 3,100 K the temperature change across the mantle, h = 2,867 km the thickness of the mantle, κ0 = 1 × 10−6 m2 s-1 the thermal diffusivity and η0 = 1.1 × 1021 Pa s is the viscosity. With the values listed above Ra = 7.8 × 107. The dissipation number that controls shear heating is , and with R0 = 6,371 km being the radius of Earth and = 1200 J kg-1 K-1 the reference heat capacity we obtain Di = 1.56. Viscosity varies with depth, composition, temperature and pressure following , with η(r) = 0.02 above 160 km depth and between 310–660 km depth, η(r) = 0.002 between 160–310 km depth (asthenosphere), and η(r) = 0.02 below 670 km depth (lower mantle). The compositional viscosity pre-factor η was varied across model cases for the basal layer (see Table 1). r is the radius, R = 3,504 km is the radius of the core, E = 275 kJ mol-1 is the activation energy, Z = 2.1 × 10−6 m3 mol-1 is the activation volume, R = 8.31 J mol-1 K-1 is the universal gas constant, T is the dimensional temperature, Toff = 452 K is a temperature offset and TCMB = 3380 K is the temperature at the core-mantle boundary. The viscosity pre-factor, activation energy, activation volume and temperature offset were selected to obtain variations in viscosity over three orders of magnitude (viscosity variations were limited to the range 1.1 × 1020 Pa s—2.2 × 1023 Pa s) across the range of temperatures and pressures (Fig 2).

Table 1

Parameters varied across mantle flow models.

Case number	Reconstruction	Basal layer viscosity pre-factor	Initial slab depth (km)	Basal layer density δρb (%)
C1	M21NNR	10	1,000	+1.02
C2	M21NNR	10	1,000	+1.43
C3	M21NNR	1	550	+2.05
C4	M21NNR	1	550	+1.02
C5	M21NNR	10	550	+1.02
C6	M21	10	1,000	+1.02

Values differing from Case 1 are in bold.

Fig 2

Temperature and viscosity with depth.

Horizontally averaged present-day mantle temperature (left-hand side) and viscosity (right-hand side) for Case 2. Solid lines show the average, and dashed lines show the minimum and maximum. The grey line on the right-hand side panel is a viscosity profile adjusted to fit the geoid and post-glacial rebound [54].

The initial condition at 1 Ga consists of an adiabatic temperature profile between two thermal boundary layers, with slabs inserted down to either 550 km or 1,000 km depth (Table 1) with a dip of 45° down to 425 km and a dip of 90° below 425 km depth. For cases with slabs inserted down to 1,000 km depth, slabs are twice as thick in the lower mantle than in the upper mantle to account for advective thickening in the more viscous lower mantle [55]. The basal thermal boundary layer is initially 225 km thick and includes a 113-km thick layer [2% of the volume of the mantle; 56] that is denser than ambient mantle. The density contrast between the basal material and ambient mantle, obtained using tracers, was varied between +1.02% and +2.05% across model cases (Table 1).

As slabs sink into the mantle, the initially homogenous basal thermochemical layer is deformed [36, 37] and forms basal mantle structures of changing area. To assess this change, we quantify the evolution of the fractional area (f) of high temperature clusters in the models; we also compare it to the area of slow-velocity clusters in tomographic models.

2.4. Imaging lower mantle structures through cluster analysis

We consider six global S-wave tomographic models that use variations in seismic wave velocities to map mantle structure with respect to a reference model such as PREM [57]: SAW24B16 [58], S362ANI [59], GyPSuM-S [60], S40RTS [61], Savani [62] and SEMUCB-WM1 [63]. As in ref. [64], the structure of the lower mantle is represented by cluster analysis of seismic wave velocities in tomographic models (only for the present-day) and mantle temperature in mantle flow models (present-day and back in time in 20 Myr increments), respectively. We use k-means clustering [65], a procedure that minimizes the variance in squared Euclidean distance between vectors to separate ~200,000 points equally-spaced on Earth’s surface into two clusters with similar variations between 1,000–2,800 km depths in either seismic velocity (for tomographic models, Fig 3A and 3C) or mantle temperature (for mantle flow models, Fig 3B and 3D).

Fig 3

Cluster analysis of tomographic and mantle flow models.

A, B/ spatial distribution of two lower mantle regions revealed by k-means cluster analysis between 1,000 km and 2,800 km depths for A/ tomographic model S40RTS [61] and B/ mantle flow model Case 2. Coastlines are shown in black, and the grey lines indicate a value of five (solid) and a value of one (dashed) in a vote map of low-velocity regions from five S-wave tomographic models [66]. These figure panels were created with GMT6 [47] and the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastlines are republished from [48] under a CC BY license, with permission from Paul Wessel, original copyright 1996. C, D/ depth profiles of clustered properties in each cluster: C/ velocity profiles in high-velocity (fast, in blue) and low-velocity (slow, in red) regions for S40RTS and D/ temperature profiles in low-temperature (cold, in blue) and high-temperature (hot, in red) regions for Case 2. The solid curves are the mean, and the transparent envelopes are the associated standard deviation. None of the data in this figure are proprietary.

2.5. Quantifying model success

2.5.1. Success for the present-day: Match with tomographic models

We assess the success of mantle flow models in reproducing the present-day structure of the mantle by quantifying the match between cluster maps for tomographic models and flow models as in ref. [64]. The accuracy Acc = (TP+TN)/A involves the area of true positive (TP, orange in Fig 4) match (where high-temperature clusters from a given flow model intersect slow-velocity clusters from a given tomographic model), and the area of true negative (TN, grey in Fig 4) match (where low-temperature clusters from a given flow model intersect fast-velocity clusters from a given tomographic model) over the total area A. The sensitivity (sometimes called recall or true positive rate) S = TP/(TP+FN) involves the true positive area (TP, orange in Fig 4) and the false negative (FN, blue in Fig 4) area (where low-temperature clusters from a given flow model intersect slow-velocity clusters from a given tomographic model).

Fig 4

Spatial match of lower mantle structure between mantle flow and tomographic models with respect to tomographic model S40RTS.

A/ Comparison of present-day mantle structure for S40RTS and for mantle flow model Cases 1–6 as indicated. Orange (true positive) indicates low-velocity tomography and high-temperature regions, grey (true negative) indicates high-velocity tomography and low-temperature regions, green (false positive) indicates high-velocity tomography and high-temperature regions and blue (false negative) indicates low-velocity tomography and low-temperature regions. B/ Comparison of mantle structure imaged by different tomographic models, using S40RTS as a reference. Orange (true positive) indicates low-velocity regions for both models, grey (true negative) indicates high-velocity regions for both models, green (false positive) indicates high-velocity regions for S40RTS and low-velocity regions for other models and blue (false negative) indicates low-velocity regions for S40RTS and high-velocity regions for other models. This figure was created with GMT6 [47].

2.5.2. Success and predictions in the past: Distance between reconstructed kimberlite eruptions and model high kimberlite potential areas

We reconstruct the three considered sets of thick or ancient lithosphere back in time and compute the intersections of these polygons with either present-day LLSVPs (assumed to be stationary and rigid back in time) mapped from cluster-analysis of the six considered S-wave tomographic models, or time-dependent basal mantle structures mapped with cluster-analysis in the seven considered flow models. These intersections are considered regions where kimberlites are likely to occur (high kimberlite potential regions). We compute time-dependent potential maps as the distance from these intersections in the mantle reference frame, as well as relative potential maps that summarise potential over a given period in the plate frame of reference; in the latter, potential is obtained by summing intersections between thick or ancient lithosphere and basal mantle structures over time, with values of one within intersections (prospective area at a given time) and zero outside. “Relative” potential is obtained by normalising the result with respect to the highest.

To evaluate the success of models in predicting past kimberlite eruptions, we measure the distance between each reconstructed kimberlite and the nearest intersection between thick or ancient lithosphere and basal mantle structure. These minimum distances, cumulated over time, are summarised in one empirical distribution function for each combination of thick or ancient lithosphere and basal mantle structures, and the median and average minimum distances from these EDFs are used as indicators of model success.

3. Results

3.1. Area of basal mantle structures in tomographic and mantle flow models

The fractional area f of basal mantle structures affects the match between mantle flow models and tomographic models, as well as the minimum distances between kimberlites and potential areas because models with larger f result in larger potential areas, leading to smaller minimum distances to kimberlites.

For LLSVPs mapped as slow-velocity clusters from tomographic models, f ranges between 0.31 and 0.39 to the exception of GyPSuM-S for which the area is 0.52 (Fig 5). In the mantle flow models, the fractional area covered by the high-temperature cluster changes from f = 1.0 in the initial condition until it stabilizes somewhere between 0.2 and 0.5 depending on the density of the basal layer (Fig 5A and Table 1). Over time, sinking slabs reduce the fractional area of the high-temperature cluster. This effect is most pronounced in cases with slabs initially inserted down to 1,000 km depth (Cases 1, 2 and 6 in Fig 5A). In all cases, the area of basal mantle structures stabilises by about 640 Ma (360 Myr into the model run), which is about twice the time it takes for slabs to sink to the deep mantle in these models [67]. For this reason, we analyse models for the last 640 Myr.

Fig 5

Fractional area covered by hot/slow basal mantle structures in mantle flow and tomographic models.

A/ Evolution of the fractional area of basal mantle structures f in flow models (blue circles and lines) and for tomographic models (red diamonds at present-day). The dashed vertical line is at 640 Ma. B/ Fractional area of basal mantle structures averaged from 640 Ma (, blue diamonds with standard deviation shown as error bars on the left-hand-side panel, sorted by decreasing ), compared to present-day f for tomographic models (open green diamonds on the right-hand-side panel, sorted by increasing f).

In the flow models, the fractional area covered by basal mantle structures f also depends on the buoyancy of the basal layer: it is larger for a denser basal layer ( ≈ 0.47 for Case 3 with δρ = +2.05%) and smaller for a less dense basal layer ( ≈ 0.26 for Cases 1, 4 and 5 with δρ = +1.02%). Case 2 with δρ = +1.43% results in ≈ 0.37 which is consistent with tomographic models (0.31 < f < 0.39). For a given density of the basal layer, f is larger for M21 ( ≈ 0.31 for Case 6) than for M21NNR ( ≈ 0.26 for Cases 1, 4 and 5), due to the different location of subduction zones with different reference frames (Fig 4A).

3.2. Match between predicted mantle structure and tomographic models

The shape and size of basal thermochemical structures differ between model cases. Flow models that better match present-day tomographic models can be extrapolated back in time with more confidence, because the structure of the lower mantle is shaped by the history of subduction. To evaluate how the mantle structure from flow models matches that from tomographic models, we compute the accuracy Acc = (TP+TN)/A and sensitivity S = TP/(TP+FN) between tomographic models and between each mantle flow model and tomographic model (Fig 6). The accuracy is a global match, whereas the sensitivity is a match between imaged and predicted and slow/hot basal mantle structure.

Fig 6

Quantitative match between predicted and imaged basal mantle structures.

Match between hot basal mantle structures predicted by flow models and LLSVPs in seismic tomography models as mapped between 1,000 km and 2,800 km depths using cluster analysis. Left-hand-side panels show the sensitivity and right-hand-side panels the accuracy. In the top panels, quantities are derived between six flow models and six tomographic models, and in the bottom panels they are derived between six tomographic models.

Overall, the match between tomographic models (0.58 < S < 0.86 and 0.76 < Acc < 0.88, excluding auto-correlation for which both scores are equal to one, Fig 6) is better than the match between mantle flow and tomographic models (0.33 < S < 0.74 and 0.54 < Acc < 0.78). This result is expected since tomographic models are generally similar to one another in the way they were built, whereas the forward mantle flow models that cover one billion years of tectonic history are completely independent from tomographic models. Within tomographic models, the largest average sensitivities and accuracies are obtained for tomographic models S40RTS and SEMUCB-WM1. In contrast, GyPSuM-S appears to be an outlier for which the accuracy and sensitivity are systematically lower than for other tomographic models, because f is greater in that tomographic model (Figs 4–6). Amongst mantle flow models, the sensitivity is largest between Case 1, 4 and 5 and S40RTS, Savani and S362ANI; it is also relatively large for Case 2. The accuracy is largest for the same mantle flow models but with S40RTS and SEMUCB-WM1. The accuracy and sensitivity are lowest for Case 6 that is based on reconstruction M21, in which net lithospheric rotation present in the surface boundary conditions induces wholesale motion of the lower mantle and basal mantle structures [44]. For this reason, we focus on M21NNR in the following.

There is a trade-off between f and S and Acc: model cases in which the basal layer is less dense (Cases 1, 4 and 5) result in larger S and Acc (match to tomographic models), however in these cases f is smaller than in tomographic models. Case 2 (δρ = +1.43% and reconstruction M21NNR), for which f is comparable to tomographic models and Acc and S are reasonably large, is our preferred mantle flow model case. In the following, we first illustrate results for Case 2 and tomographic models S40RTS (which is the best tomographic model based on Acc and S) before summarising results for all models.

3.3. Reconstructing kimberlites, thick or ancient lithosphere and basal mantle structures back to 640 Ma

As an example, we reconstruct kimberlites occurrences from ref. [22] and lithosphere thicker than 150 km, and consider basal mantle structures in S40RTS and Case 2 back to 640 Ma (Fig 7 and S1 Video). S40RTS better captures reconstructed southern African kimberlites than Case 2 back to 80 Ma (S1 Video). However, before 80 Ma, Case 2 better captures the peak of South African kimberlite activity at 100 Ma. All flow models based on reconstruction M21NNR capture the southern African kimberlite blooms back to around 200 Ma, as illustrated by comparing S40RTS and Case 2 at 120 Ma, 160 Ma and 240 Ma (Fig 7). S40RTS and Case 2 are both generally consistent with kimberlite eruptions on the Siberian block back to 360 Ma. S40RTS does not match the kimberlite records in Antarctica at 120 Ma and Australia at 240 Ma, whereas Case 2 does (Fig 7). Case 2 almost captures the kimberlite occurrence in North China at 240 Ma (Fig 7G), although the North China craton is thinner than 150 km following lithospheric delamination [68] and the model basal mantle structure is offset to the East of the North China block from ref. [38] (Fig 7G).

Fig 7

Evolution of the location of slow/hot basal mantle structures and kimberlite eruptions.

Kimberlite eruption locations (stars) sampled in 20 Myr increments from ref. [22] continent outlines from ref. [38] (white transparent polygons) and lithosphere thicker than 150 km derived from ref. [5] (dark grey transparent polygons), all reconstructed at 120 Ma (A, E), 160 Ma (B, F), 240 Ma (C, G) and 360 Ma (D, H) with reconstruction M21NNR. The background shows cluster maps with basal mantle structures in pink for tomographic model S40RTS (A-D) and for mantle flow model Case 2 (E-H). The intersections of white and pink areas are high-potential regions for kimberlites. Labels: A/ “Si”: Siberia; B/ “An”: Antarctica; C/ “Afr”: Africa, “Au”: Australia, “Ca”: Canada; D/ “N.Am”: North America. This figure was created with GMT6 [47] and the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastlines are republished from [48] under a CC BY license, with permission from Paul Wessel, original copyright 1996.

Neither S40RTS nor Case 2 capture kimberlite magmatism in the Americas at 120 Ma. However, Case 2 better captures kimberlite occurrences than S40RTS in southern Africa and Antarctica at that time. S40RTS encompasses the kimberlite occurrence in Western Australia at 120 Ma whereas Case 2 does not (Fig 7A and 7E). In contrast, Case 2 captures the kimberlite occurrence in the Northern Territory of Australia at 240 Ma and 360 Ma, whereas S40RTS does not (Fig 7C and 7G).

At 160 Ma (Fig 7B and 7F), both tomographic and flow models capture some north-eastern Canadian kimberlites although neither models matches the North American kimberlites until 170 Ma; as in previous work [25]. At 240 Ma, S40RTS captures Siberian and north-eastern Canadian kimberlites (Fig 7C), and Case 2 captures kimberlites in Siberia, Africa, and Australia (Fig 7G). Both models capture some Siberian kimberlites at 360 Ma (Fig 7D and 7H).

3.4. Global kimberlite potential maps

3.4.1. Global kimberlite potential maps in the mantle frame of reference

We next present the time-dependent distance to polygons defined by the intersection of reconstructed lithosphere thicker than 150 km and basal mantle structures as mapped using cluster analysis. In these kimberlite potential maps, the distance is zero within the intersection and increases away from the intersections (Fig 8 and S2 Video).

Fig 8

Global kimberlite potential maps and kimberlite eruptions through time.

Distance to the intersection between basal mantle structures and reconstructed lithosphere thicker than 150 km [5] at 120 Ma (A, E), 160 Ma (B, F), 240 Ma (C, G) and 360 Ma (D, F) for tomographic model S40RTS (A-D) and mantle flow model Case 2 (E-H). Continent outlines from ref. [38] are shown as transparent white polygons, and kimberlite eruptions sampled in 20 Myr increments from ref. [22] are shown as stars; both are reconstructed using M21NNR. Labels: A/ “Gr”: Greenland, “Si”: Siberia; B/ “An”: Antarctica, “N.Eu”: Northern Europe; C/ “Afr”: Africa, “Au”: Australia, “Ca”: Canada, “In”: India; D/ “N.Am”: North America, “N.Ch”: North China, “S.Am”: South America. This figure was created with GMT6 [47] and the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastlines are republished from [48] under a CC BY license, with permission from Paul Wessel, original copyright 1996.

Kimberlite potential maps for S40RTS and Case 2 both indicate southern African, eastern Canada, Greenland, northern Europe, the eastern African coast, and the west coasts of India and Madagascar as high-potential areas for kimberlite magmatism over the last 240 Myr. Again, the predicted kimberlite potential before 80 Ma in southernmost Africa is smaller for S40RTS than for Case 2 (S2 Video). S40RTS predicts lower kimberlite potential than Case 2 for South America, Australia and Antarctica at 360 Ma (Fig 8D and 8H). Since reconstructed kimberlites are generally within or close to high kimberlite potential areas (Fig 8), we next quantify the distance of reconstructed kimberlites to high kimberlite potential areas to evaluate the success of different models.

3.4.2. Global kimberlite potential maps in the plate frame of reference

We next present relative kimberlite potential over a given period in the plate frame of reference and within presently emerged continents (Fig 9A–9D). This shows that for reconstruction M21NNR, the largest predicted kimberlite potential is in central and western Africa for S40RTS (Fig 9A and 9C), whereas it is in southern Africa for Case 2 (Fig 9B and 9D). This result is similar over 320 Myr (Fig 9A and 9B) and 640 Myr (Fig 9C and 9D). Because more regions are at the intersection between thick lithosphere and basal mantle structures when a longer model duration is considered, the area of high-potential regions is larger over 640 Myr (Fig 9C and 9D) than over 320 Myr (Fig 9A and 9B). Overall, the relative kimberlite potential at kimberlite locations is greater for Case 2 than for S40RTS (Fig 9E and 9F): the peak in the number of kimberlites is for a relative kimberlite potential of about 0.3 for S40RTS (median 0.21 and mean 0.19) and about 0.55 for Case 2 (median 0.24 and mean 0.35) over 320 Ma (Fig 9E); over 640 Ma (Fig 9F) the peak is at ~0.15 for S40RTS (median 0.16 and mean 0.18) and ~0.55 and ~0.75 for Case 2 (median 0.22 and mean 0.33). We note that there are also significant peaks at low-potential (<0.15), and that these peaks are higher for Case 2 than for S40RTS (Fig 9E and 9F).

Fig 9

Predicted relative kimberlite potential.

A-D/ Global relative kimberlite potential maps for lithosphere thicker than 150 km from 320 Ma (A-B) and from 640 Ma (C-D) for tomographic model S40RTS (A, C) and Case 2 (B, D). Kimberlites from ref. [22] are shown as black circles for the last 320 Ma (A-B) and 640 Ma (C-D). These figure panels were created with GMT6 [47] and the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastlines are republished from [48] under a CC BY license, with permission from Paul Wessel, original copyright 1996. E-F/ Distribution of relative kimberlite potential at kimberlite locations shown in A-D for S40RTS and Case 2 from 320 Ma (E) and from 640 Ma (F).

3.5. Distance between kimberlites and predicted high kimberlite potential areas

3.5.1. Empirical distribution functions

We sample the distance to the closest high kimberlite potential region (referred to as minimum distance) at each reconstructed kimberlite locations (Fig 8) and report results cumulated over time as empirical distribution functions (EDFs) for each combination of basal mantle structures and thick or ancient lithosphere and for the last 640 Ma (Fig 10). These EDFs show the likelihood that a given proportion of kimberlites falls within a given minimum distance from a high kimberlite potential region. Thus, EDFs with smaller minimum distances reflect cases for which reconstructed kimberlites are closer to predicted high kimberlite potential areas, which indicates more successful models that better account for the kimberlite record. Models with a greater likelihood that distances are equal to zero are also more successful. As expected, distances are smallest for tomographic model GyPSuM-S and flow model Case 3 (Fig 10) because these models result in larger basal mantle structures (Figs 4 and 5).

Fig 10

Distances between kimberlite eruption locations and high kimberlite potential areas.

Empirical Distribution Functions (EDFs) of cumulative minimum distances from 640 Ma between reconstructed kimberlites (sampled in 20 Myr increments) and predicted high kimberlite potential areas for lithosphere thicker than 150 km [5] reconstructed using M21NNR, for all six tomographic models, and mantle flow models Cases 1–5 (that are based on M21NNR).

3.5.2. Metrics from empirical distribution functions

To summarise EDFs, we report the mean minimum distance and associated standard deviation for each distribution, as well as the median minimum distance (distance for which the likelihood is equal to 0.5; Fig 11). The mean and median minimum distances confirm a link between reconstructed kimberlite eruption locations, reconstructed thick or ancient lithosphere, and basal mantle structure locations: median minimum distances indicate that 50% of kimberlites are within ~800 km of high kimberlite potential areas for all models (Fig 11A and 11C) and mean minimum distances are less than 1500 km (Fig 11B and 11D). Mean and median minimum distances are similar for mantle flow models and tomographic models (Fig 11), suggesting that mobile and deforming basal mantle structures are as consistent with the kimberlite record as stationary and rigid ones. This first order result, as well as trends between models, are similar over 320 Ma (Fig 11A and 11B) and over 640 Ma (Fig 11C and 11D). Distances are largest (least favourable) for Cases 4 and 5, in which slabs are initially inserted down to 550 km depth and the density of the basal layer is intermediate (+1.02%; Table 1). Standard deviations are large for tomographic model S362ANI and flow model Cases 1, 4 and 5 (Fig 11), reflecting more spread distributions (Fig 10) due to smaller f (Fig 4B). In contrast, mean and median distances are smallest for tomographic model GyPSuM-S and mantle flow Case 3, however these two cases do not match other tomographic models (Fig 6) due to large basal mantle structures (Fig 5B). Interestingly, mantle flow model Case 2 results in relatively small minimum distances (Fig 11) as well as being compatible with tomographic models (Figs 5 and 6). For Case 2, 50% of kimberlites are within ~250 km of high kimberlite potential areas over 320 Ma (Fig 11A), and 50% of kimberlites are within ~350 km of high kimberlite potential areas over 640 Ma (Fig 11C). When different lithospheric extents are used, trends are similar, and mean and median distances are larger when using the smaller blocks of Archean tectonothermal age from ref. [46], and smaller when using the larger tectonic blocks are used from ref. [38] (Fig 12). These trends are expected, and we note that the difference varies between cases, because the distance depends on the location of basal mantle structures as well as the location of reconstructed thick or ancient lithosphere (Figs 7 and 8, S3 and S4 Videos). When reconstruction M21 (in which the net rotation of the lithosphere is not removed) is used, mean and median distances are approximately twice as large than with reconstruction M21NNR for tomographic models (S1 Fig).

Fig 11

Results for lithosphere thicker than 150 km [5] and reconstruction M21NNR.

Metrics from EDFs for lithosphere thicker than 150 km [5] and reconstruction M21NNR for mantle flow and tomographic models from 320 Ma (A-B) and from 640 Ma (C-D). A, C/ Median and B, D/ average minimum distances with one standard deviation as error bar. In each panel, mantle flow model cases are to the left and tomographic models to the right of the dashed vertical line.

Fig 12

Results for different lithospheric extents and reconstruction M21NNR.

Same as Fig 11 but for tectonothermal ages greater than 2.5 Ga from ref. [46] (open downward pointing triangles, “Archean (A06)”), lithosphere thicker than 150 km [5] (diamonds, “150 km thick”) and tectonic blocks from ref. [38] (open upward pointing triangles, “Blocks (M21)”), using reconstruction M21NNR.

4. Discussion

4.1. Model predictions: Regions of potential interest for large-scale kimberlite exploration

The global kimberlite potential maps presented in this study (Figs 7 and 8, S1–S4 Videos) may be used for global-scale kimberlite exploration, placing them at the base of the hierarchy of scale-dependent targeting processes [69]. Thick or ancient lithosphere that have repeated intersections with model slow or hot clusters are more likely to have undergone kimberlite magmatism and are therefore of potential interest for large-scale kimberlite exploration. Kimberlite potential maps for S40RTS and Case 2 capture regions of prominent kimberlite magmatism such as Africa [70] and Siberia [71] (Figs 8 and 9D). Case 2 also predicts some kimberlite activity in northwest Australia, Greenland, northeast Canada and Northern Europe (Figs 8 and 9D). Interestingly, Case 2 indicates Antarctica as a location of kimberlite potential for the last 640 Ma (Fig 9D), although there is only one known instance of kimberlite magmatism on the continent [72].

4.2. Mobility of basal mantle structures

Previous work has been shown that the reconstructed location of ~1,100 kimberlites over the last 320 Myr are within ~1,650 km of the edge of the African LLSVP as imaged by SMEAN at 2,800 km depth [25]. Interestingly, ~240 kimberlite occurrences from the Slave Province of Canada were more distant from the edge of the African LLSVP. Nevertheless, reconstructed locations of 25 LIPs erupted in the last 297 Myr were also close to the edge of LLSVPs, leading to the hypotheses that 1/ mantle plumes arise from LLSVPs, and 2/ that LLSVPs could have been stationary over the last 300 Myr, and possibly longer [25, 73, 74].

Here, we assumed that a relationship between basal mantle structures and kimberlites exists, and we considered two main scenarios: one in which lower mantle structures are mobile (represented by mantle flow models) and one in which lower mantle structures are stationary over time (represented by tomographic models). Overall, we found that minimum distances between kimberlites and high-potential areas for kimberlites are similar for the two scenarios (Figs 11 and 12 and S1 Fig), and that mobile basal structures better predict the abundance of kimberlites in southern Africa than stationary basal mantle structures (Figs 7–9). This suggests that kimberlites are as close to mobile basal structures as they are to stationary basal mantle structures, and that there is no need to hypothesise that basal mantle structures have been stationary over time, which is consistent with the mobility of tectonic plates, subduction zones [38, 42] and mantle plumes [75, 76] over time.

4.3. Limitations

Our global kimberlite potential maps are based on the assumption that kimberlites are associated with deep mantle upwelling above basal mantle structures. We did not consider shallower controls such as fractures and faults cutting through the lithosphere [21]. In addition, not all kimberlites are necessarily linked to the deep mantle, and kimberlite eruptions may also occur due to partial melting at the base of the lithosphere or above subducting slabs [12, 77] and can be emplaced through weak zones within the lithosphere or close to plate boundaries [21]. Indeed, the lithospheric mantle might be the source of some South African orangeites (formerly known as group II kimberlites) [78], and kimberlites from the Slave Province that are not closely associated with the African LLSVP (Figs 7–9) could be associated to the history of subduction under North America [79]. Recent geochemical fingerprinting of kimberlites revealed that some kimberlite erupted above basal mantle structures could be derived from a primitive, depleted mantle source, whereas kimberlite erupted away from basal mantle structures could be derived from a source that contains deeply subducted crustal material [80]. Even though upwelling of hot mantle from the deep Earth may provide the source of heat leading to kimberlite eruption, the physical processes that could link basal mantle structures [>2000 km deep; 26] and kimberlite eruption from ~120–300 km depth [12, 13] are not yet firmly established.

Our global kimberlite potential maps also rely on a definition of relevant lithospheric extent. While diamondiferous kimberlites are generally associated with cratons [e.g., 12], the definition of craton may not be limited to the age of the overlying crust [46] and the geochemical composition of the subcontinental lithospheric mantle [45], and a recent definition suggests that lithosphere 150 km thick could be an appropriate proxy for craton extent [4]. We follow this definition here, which results in larger blocks than using tectonothermal ages [46] (Fig 1), leading to smaller distances between kimberlite eruptions and kimberlite potential maps (Fig 12). We derived lithosphere thicker than 150 km from published lithospheric thickness models [5] and we note that this simple proxy does not consider modified cratons [4], most notably the North China Craton that has lost its lithospheric root [68] and contains kimberlites (Fig 1).

The presented reconstructions of past mantle flow do not predict the location of kimberlite eruptions. Modelling kimberlite eruptions requires a horizontal resolution of a few hundred meters, modelling two-phase flow (melt and solid) and considering viscosity contrast over 26 orders of magnitude between the melt with viscosity 0.01 Pa s [81] and the lithosphere with viscosity about 1 x 1024 Pa s. The presented time-dependent global mantle flow models achieve a horizontal resolution of 50 km, viscosity contrast over three orders of magnitude, and do not consider melting. Achieving kilometre-scale resolution is possible in instantaneous global mantle flow models [82] but remains prohibitively expensive in time-dependent models. Viscosity contrasts of six orders of magnitude have been achieved in global mantle flow models [83], however these models were not paleogeographically-driven so could not be as readily compared to the geological record as the models presented here. Modelling magmatic processes in 3D global spherical geometry remains a frontier.

Our results are sensitive to the density of the basal mantle structures, which is a key control on their fractional area (f). In our preferred Case 2, basal structures were intrinsically 1.43% denser than surrounding mantle. This density is comparable to previous models [e.g., 37, 84], and with that inferred from Earth’s free oscillations [up to 1.7%; 85], however, it is larger than those estimated from tidal tomography 0.5% [29]. Further research is required to determine the density of basal mantle structures since seismically-filtered mantle flow models [86] and tomography including Stoneley modes [87] suggest that they could be purely thermal.

We modelled mantle flow forward in time, starting from an uncertain initial condition derived from tectonic reconstructions. A limitation of this approach is that the predicted present-day structure of the lower mantle is similar (up to 78% similarity, Fig 6) rather than identical to the structure of the lower mantle imaged by tomographic models. An alternative approach would consist of reconstructing mantle flow back in time from the present-day structure of the mantle inferred from tomographic models. However, such approaches have so far been limited to the last ~100 Myr [88], and seismically fast structures in tomographic images may constrain plate motions back to 250 Ma at most [67, 89]. We used the only available full-plate tectonic reconstruction extending back to 1 Ga, and considered the original reconstruction [M21; 38] and a version of the reconstruction in which net rotation of the lithosphere was removed (M21NNR). The latter is appropriate for mantle flow modelling since net rotation of the lithosphere is thought to arise from lateral viscosity variations within the mantle [43, 90]. Indeed, imposed lithospheric net rotation results in induced net rotation of the lower mantle and results in a poor match of lower mantle structures with tomographic models [44]. Nevertheless, lithospheric net rotation is unlikely to be equal to zero; it could be between 0.023°/Myr [91] and 0.26°/Myr [92]. Further work is required to develop new tectonic reference frames for paleogeographically-constrained global mantle flow models.

5. Conclusion

We assumed that a link exists between upwelling above basal mantle structures and kimberlite magmatism to investigate the relationships between kimberlite eruptions, thick or ancient lithosphere, and lower mantle structure over the last 640 Myr. To investigate the hypothesis that basal mantle structures might have been stationary over time, we considered lower mantle structure to be either stationary and given by tomographic models, or mobile and predicted by reconstructions of past mantle flow over the past billion years. We showed that the present-day structure of the lower mantle predicted by flow models is compatible with that imaged by tomographic models. We created kimberlite potential maps by computing the time-dependent intersection between reconstructed thick or ancient lithosphere and lower mantle structures. We measured the distance of reconstructed kimberlite eruption locations to high kimberlite potential regions, which confirmed that kimberlites are generally close to high kimberlite potential regions, although this result depends on the reference frame being used and the lithosphere extent considered for reconstruction; for the preferred mantle flow model case 50% of kimberlites were within ~250 km of high kimberlite potential areas over 320 Ma and within ~350 km of high kimberlite potential areas over 640 Ma. We found that mobile basal mantle structures better predicted the peak in kimberlite eruption in southern Africa at ~100 Ma than stationary mantle structures. Overall, these results suggest that mobile basal mantle structures are as consistent with kimberlite eruptions as stationary ones. We therefore favour a model in which basal mantle structures have been mobile over time, which is consistent with observations and models suggesting that tectonic plates, subduction zones and mantle plumes have also been mobile over time.

(MP4)

Click here for additional data file.

(MP4)

Click here for additional data file.

(MP4)

Click here for additional data file.

(MP4)

Click here for additional data file.

(TIF)

Click here for additional data file.

7 Feb 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Shuan-Hong Zhang, PhD

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf.

2. Thank you for stating the following in the Acknowledgments Section of your manuscript:

[NF and OFB were supported by Australian Research Council grant LP170100863. This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government. This research was supported by the Australian Government's National Collaborative Research Infrastructure Strategy (NCRIS), with access to computational resources provided by the National Computational Infrastructure (NCI) through the National Computational Merit Allocation Scheme. Access to NCI was partly supported by resources and services from the University of Wollongong (UOW). Key model results and scripts are available at https://doi.org/10.5281/zenodo.5760115.]

We note that you have provided funding information that is not currently declared in your Funding Statement. However, funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.

Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:

[NF and OFB were supported by Australian Research Council (https://www.arc.gov.au/) grant LP170100863.]

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

3. We note that Figures 1, 3, 4, 7, 8, 9, S1 and S2 in your submission contain [map/satellite] images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

a) You may seek permission from the original copyright holder of Figures 1, 3, 4, 7, 8, 9, S1 and S2 to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

b) If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

The following resources for replacing copyrighted map figures may be helpful:

USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/

The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/

Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.html

NASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/

Landsat: http://landsat.visibleearth.nasa.gov/

USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#

Natural Earth (public domain): http://www.naturalearthdata.com/

Additional Editor Comments:

Dear Dr. Nicolas Flament,

Thank you very much for submitting your interesting work to PLOS One. I have received a very detail review on your manuscript. Both reviewer and I found that the topic is interesting and novel and are worthy of publication in PLOS One after a moderate revision. The reviewer has given a numerous comments and suggestions on your manuscript (including those in the attached PDF files), which I believe are very useful for improvement of quality of the manuscript. Please consider all these comments and suggestions carefully during revision.

I am looking forward to receiving your revision manuscript soon.

Yours sincerely,

Shuan-Hong Zhang

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This work builds up on the hypothesis of Torsvik et al. (2010 Nature) that kimberlites are mostly derived from plumes stemming from the margins of LLSVPs (i.e. seismically defined thermochemical piles above the core-mantle boundary), and assesses the likelihood that the position of mobile thermochemical piles better match the location of Phanerozoic kimberlites compared to fixed thermochemical piles. This is a good idea to test because it is increasingly recognised that LLVSPs are probably not fixed but must be swept around by tectonic processes and mantle convection. To achieve this goal the Authors employ a fluid dynamic model to simulate modification of the thermochemical piles position over the last 640 Myr and identify as zones of kimberlite generation those where the thermochemical piles overlay with the position of cratons (the latter monitored using state-of-the-art plate reconstruction models). They find a good match between kimberlite location and these possible zones of kimberlite generation, and conclude that using mobile rather than fixed thermochemical piles marginally improve the previously proposed hypothesis that most kimberlites are linked to upwellings from these deep mantle structures. They further conclude that this approach could be exploited to identify prospective regions for diamond exploration.

The manuscript is well written and easy to follow. The conclusions are justified by the modelling results. However, I must stress I am not a geodynamic modeller so I could not evaluate the numerical model and some of the boundary conditions. As my expertise is in kimberlite and mantle geochemistry, I have focused my review on the overall approach taken by the Authors while also assessing the soundness of the background information they provide.

I believe the manuscript has one major problem: It relies on craton outlines (those of Merdith et al., 2021) which, to me, have little geological meaning because they extend way beyond known cratonic mantle occurrences and, locally, even into oceanic domains. While the Authors also employ the craton outlines of Artemieva (2006), but hardly discuss this alternative except for one statement, there are more recent and robust definitions of craton (and especially their mantle roots) such as Pearson et al. (2021 Nature), which provide a much better fit for the theme under discussion here.

The problem with employing the craton outlines of Merdith et al becomes clear in figure 9. There is no thick lithosphere (say cratons) in most of Europe (except parts of Ukraine and NW Russia) or in most of US, including all eastern US (e.g., Pearson et al., 2021), which is the opposite of what figure 9 implies (as well as figure 1). A better definition of thick lithosphere will considerably improve the outcome of this work. This issue might seem trivial, but it is not because changing the definition of craton completely modifies the model results. As an example, at L475-476 the Authors state that using the craton outlines of Artemieva (2006) the average distance between predicted and actual position of kimberlites increases significantly, again consistent with the gross overestimation of craton sizes in Merdith et al.

In addition, or perhaps in alternative, if the ultimate aim is to assess diamond prospectivity, it might be better to run the model using as diamond prospective regions only those where the lithosphere is thicker than ~180 km (or ~150 km if the aim is to match the formation of kimberlites) for example based on global tomographic models. This definition of diamond prospective areas is not necessarily confined to cratons but also include thick pericratonic regions (e.g., Limpopo Belt between South Africa and Zimbabwe) where important diamond mines occur (e.g., Venetia). This empirical approach would also circumvent potential problems with craton definition taken from other publications such as that of Merdith and coworkers.

My second major comment will probably sound naïve and reflects my ignorance with fluid dynamic modelling. The Authors start their model at 1 Ga and the conditions are not well specified – i.e. what is the structure of deep-mantle thermochemical piles at this time? Could you please clarify it? They make sure that the model matches the present-day extension and configuration of LLSVPs above the core-mantle boundary to assess model validity, which seems to be logical. However, would not be more insightful to start with the present-day configuration of LLSVPs, which is well established (and not hypothetical) and reconstruct how the position and shape of LLSVP changed backward in time using existing geological constraints (e.g., position of paleo-subduction zones) and mantle flow models? Again, this might sound a naïve comment, but it would be nice to hear the Authors thoughts. It might perhaps be an interesting follow-up work.

Finally, referencing to the literature specific to the origin of kimberlites is often inappropriate and so are some statements – however, I do acknowledge that none of the Authors has a background in kimberlite petrology-geochemistry, so this is not a major issue from my perspective. I have suggested several amendments in the attached pdf and some notable examples follow.

L45. This statement does not make any sense (“the average depth of kimberlite eruptions is 200 km”) and no surprise it comes from obscure Russian literature. Kimberlites originate from depths in excess of 120 km and generally 150 km, this is all we know - e.g., Giuliani and Pearson (2019 Elements) for a concise review (cited later on).

L46. RE superdeep diamonds: also Pearson et al (2014 Nature) and Tschauner et al (2021 Science). However, it should be noted that the evidence supporting a deep kimberlite origin from superdeep diamonds is controversial because these diamonds could have been transported to the base of the lithosphere before being entrained by kimberlite magmas (e.g., Harte and Cayzer, 2007 Phys Chem Min).

L47. The record of kimberlite magmatism is actually episodic not continuous, e.g., Heaman et al., 2019 Elements

L74. I think the only option to link LLSVPs and kimberlites is provided by plumes, either their peripheries or weak ones (see Enrst, 2014 for an assessment of the relationship between LIPs and kimberlites)

L108-110. Unfortunately, the database of kimberlite geochronology by Tappe et al. (2018) has two biases, the first one being geographic areas of more intense exploration for diamonds and the second one being several ages for multiple kimberlite bodies from the same cluster or field (say area). In addition, there might be an inherent preservation bias, which cannot be addressed. I would just remove this statement.

L118-120. Not all kimberlites are located on craton so using Merdith et al. (2021) to obtain craton extensions is not a good idea. Many kimberlites are located in Proterozoic mobile belts and are commonly not diamondiferous. Also, cratons are confined to continental interiors unlike the assessment by Meredith et al. shown in figure 1a.

L364-365. I am not aware of any kimberlite emplaced in Western Australia at this time (~240 Ma). Not the Authors' fault, just one of the several problems with the database of Tappe et al 2018.

L541-544. Additional factors to consider when modelling potential kimberlite formation include lithospheric stress, pre-existing translithospheric structures, plume buoyancy and composition just to name a few.

Finally, it should be kept in mind that not all kimberlites (or related magmas such as those from North China; see Tompkins et al., 1999 Proc 7th Inter Kimb Conf) are necessarily linked to the deep mantle; some might be related to shallower sources including the mantle transition zone (e.g., Kjarsgaard et al., 2017 G-cubed; Chen et al., 2020 Geology) or also the lithospheric mantle (e.g., Giuliani et al., 2015 Nature Comms). This should be noted in the manuscript in my opinion.

Figures 7 and 8 would be easier to follow if the whole continental blocks were shown rather than cratons only. Similarly, the statement at L393-394 is not easy to picture because continental blocks are not shown but only craton outlines.

Additional minor edits and suggested references are included in the attached pdf.

I hope my comments will be helpful and my apologise for the late review (6-Feb-2022)

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: PONE-D-21-39227_reviewer_mod.pdf

Click here for additional data file.

11 Apr 2022

Our replies to comments are in blue; they include line numbers with reference to the manuscript with tracked changes.

Reviewer #1: This work builds up on the hypothesis of Torsvik et al. (2010 Nature) that kimberlites are mostly derived from plumes stemming from the margins of LLSVPs (i.e. seismically defined thermochemical piles above the core-mantle boundary), and assesses the likelihood that the position of mobile thermochemical piles better match the location of Phanerozoic kimberlites compared to fixed thermochemical piles. This is a good idea to test because it is increasingly recognised that LLVSPs are probably not fixed but must be swept around by tectonic processes and mantle convection. To achieve this goal the Authors employ a fluid dynamic model to simulate modification of the thermochemical piles position over the last 640 Myr and identify as zones of kimberlite generation those where the thermochemical piles overlay with the position of cratons (the latter monitored using state-of-the-art plate reconstruction models). They find a good match between kimberlite location and these possible zones of kimberlite generation, and conclude that using mobile rather than fixed thermochemical piles marginally improve the previously proposed hypothesis that most kimberlites are linked to upwellings from these deep mantle structures. They further conclude that this approach could be exploited to identify prospective regions for diamond exploration.

The manuscript is well written and easy to follow. The conclusions are justified by the modelling results. However, I must stress I am not a geodynamic modeller so I could not evaluate the numerical model and some of the boundary conditions. As my expertise is in kimberlite and mantle geochemistry, I have focused my review on the overall approach taken by the Authors while also assessing the soundness of the background information they provide.

Thank you for this supportive summary of our manuscript.

I believe the manuscript has one major problem: It relies on craton outlines (those of Merdith et al., 2021) which, to me, have little geological meaning because they extend way beyond known cratonic mantle occurrences and, locally, even into oceanic domains. While the Authors also employ the craton outlines of Artemieva (2006), but hardly discuss this alternative except for one statement, there are more recent and robust definitions of craton (and especially their mantle roots) such as Pearson et al. (2021 Nature), which provide a much better fit for the theme under discussion here.

The problem with employing the craton outlines of Merdith et al becomes clear in figure 9. There is no thick lithosphere (say cratons) in most of Europe (except parts of Ukraine and NW Russia) or in most of US, including all eastern US (e.g., Pearson et al., 2021), which is the opposite of what figure 9 implies (as well as figure 1). A better definition of thick lithosphere will considerably improve the outcome of this work. This issue might seem trivial, but it is not because changing the definition of craton completely modifies the model results. As an example, at L475-476 the Authors state that using the craton outlines of Artemieva (2006) the average distance between predicted and actual position of kimberlites increases significantly, again consistent with the gross overestimation of craton sizes in Merdith et al.

In addition, or perhaps in alternative, if the ultimate aim is to assess diamond prospectivity, it might be better to run the model using as diamond prospective regions only those where the lithosphere is thicker than ~180 km (or ~150 km if the aim is to match the formation of kimberlites) for example based on global tomographic models. This definition of diamond prospective areas is not necessarily confined to cratons but also include thick pericratonic regions (e.g., Limpopo Belt between South Africa and Zimbabwe) where important diamond mines occur (e.g., Venetia). This empirical approach would also circumvent potential problems with craton definition taken from other publications such as that of Merdith and coworkers.

We agree that the using the outlines of Merdith et al. (2021) and referring to them as “cratonic” was problematic since they are instead continental blocks inferred to have existed for one billion years. In the revised manuscript, we have implemented the suggestion of the reviewer to consider lithosphere thicker than 150 km, which has been proposed as a first-order proxy for cratonic lithosphere (Pearson et al., 2021). We show the sensitivity of the results to the extent of the continental lithosphere by considering the outlines of Artemieva (2006) for lithosphere older than 2.5 Ga and of Merdith et al. (2021). As a result of this change, we have modified Figs. 1 and 7-12 and the associated text. We also discussed the extent of lithosphere relevant to kimberlite potential mapping in the discussion (L. 938-948).

My second major comment will probably sound naïve and reflects my ignorance with fluid dynamic modelling. The Authors start their model at 1 Ga and the conditions are not well specified – i.e. what is the structure of deep-mantle thermochemical piles at this time? Could you please clarify it? They make sure that the model matches the present-day extension and configuration of LLSVPs above the core-mantle boundary to assess model validity, which seems to be logical. However, would not be more insightful to start with the present-day configuration of LLSVPs, which is well established (and not hypothetical) and reconstruct how the position and shape of LLSVP changed backward in time using existing geological constraints (e.g., position of paleo-subduction zones) and mantle flow models? Again, this might sound a naïve comment, but it would be nice to hear the Authors thoughts. It might perhaps be an interesting follow-up work.

We have clarified in the discussion that while the initial condition of forward models is uncertain, models starting at the present day and extending back in time are limited to the last ~100 million years (possibly 250 Myr from imaged subducted slabs) and therefore not best suited for the problem at hand (L. 972-979).

Finally, referencing to the literature specific to the origin of kimberlites is often inappropriate and so are some statements – however, I do acknowledge that none of the Authors has a background in kimberlite petrology-geochemistry, so this is not a major issue from my perspective. I have suggested several amendments in the attached pdf and some notable examples follow.

Thank you for these suggestions, which we have considered and addressed as outlined below. Comments in italics were originally made directly on the pdf and have been included here for completeness and transparency.

L14: “or over” �  “or through” (L. 18)

L15: “emplacement ages” �  “emplacement ages and locations” (L. 20)

L.17: actually the whole Phanerozoic if one follows the original hypothesis of Torsvik and coworkers �  we disagree. Torsvik et al. (2010) established a relationship between kimberlites and LLSVPs for the last 320 Myr, then assumed it for the rest of the Phanerozoic.

L. 40: Pearson et al., 2021 Nature is probably a better reference for this statement �  added reference to Pearson et al. (2021) (L. 58)

L42: Some of these references are either a bit obscure or not well suited for this statement. Better Stachel and Harris (2008 Ore Geol Rev), Helmstaedt and Gurney (1995 J Geochem Explor), Kjarsgaard et al. (2019, Elements) �  references updated; thank you for these suggestions (L. 60).

L.43: ascending �  change made (L. 61).

L.43: Inappropriate references. Better Sparks et al. (2006 JVGR) and Canil and Fedortchouk (1999 EPSL) who actually calculated ascent speeds �  references updated; thank you for these suggestions (L. 60).

L45. This statement does not make any sense (“the average depth of kimberlite eruptions is 200 km”) and no surprise it comes from obscure Russian literature. Kimberlites originate from depths in excess of 120 km and generally 150 km, this is all we know - e.g., Giuliani and Pearson (2019 Elements) for a concise review (cited later on).

We have changed this statement to: “Kimberlites originate from depths in excess of 120-150 km [12], forming from melts that may pool from up to ~300 km depth [13]” (L. 63-64).

L46. RE superdeep diamonds: also Pearson et al (2014 Nature) and Tschauner et al (2021 Science).

We added Pearson et al (2014 Nature) and Tschauner et al (2018 Science) (L. 65).

However, it should be noted that the evidence supporting a deep kimberlite origin from superdeep diamonds is controversial because these diamonds could have been transported to the base of the lithosphere before being entrained by kimberlite magmas (e.g., Harte and Cayzer, 2007 Phys Chem Min).

We stated “although these diamonds could have been transported to the base of the lithosphere before being entrained by kimberlite magmas” and cited Harte and Cayzer (2007 Phys Chem Min) (L. 65-67).

L47. The record of kimberlite magmatism is actually episodic not continuous, e.g., Heaman et al., 2019 Elements

“continuous” has been changed to “episodic” and references have been updated (L. 68).

L. 48: incomplete reference in the reference list �  reference updated (L. 68).

L.52: Another rather obscure reference. Perhaps better Giuliani and Pearson 2019 which include a map of kimberlite distribution �  references updated (L. 72-73).

L.55: (see also L. 17) the relationship was established back to 320 Ma, not back to 540 Ma

L. 61: If not mistaken, the best assessment of the size of LLSVPs is that of Cottaar and Lekic (2016 Geophys J Intern) �  we disagree. The cluster analysis method used by Cotaar and Lekic leads to much larger LLSVP volumes than derived from individual tomographic models.

L.64: comma added "e.g. ,” (L. 98 and elsewhere throughout the manuscript).

L. 70: also the various papers of Campbell and Griffith, e.g., 1990 EPSL, which come much earlier than Courtillot et al. �  we added a reference to Campbell & Griffiths 1990 (L. 104).

L.72: also Ernst (2014): Large Igneous Provinces �  we added a reference to Heaman & Kjarsgaard (2000) (L. 106).

L.73: [on average from ~ 200 km depth; 3] �  replaced with [from 120-300 km depth; 12, 13] (L. 107).

L74. I think the only option to link LLSVPs and kimberlites is provided by plumes, either their peripheries or weak ones (see Enrst, 2014 for an assessment of the relationship between LIPs and kimberlites)

We agree that heat must be advected from the deep mantle to the surface by mantle upwelling, however, we note that kimberlites are not systematically associated with large igneous provinces.

L.78: In this context the work of Conrad et al. (2013 Nature) is probably relevant, i.e. LLSVP might be fixed only over no more 100s Myr �  Conrad et al. (2013) did not present geodynamic models. They argued that LLSVPs are stationary and control mantle flow based on an analysis of surface plate motions over time and present-day mantle flow. This work is at odds with the statement: “global geodynamic models have revealed that similar hot basal mantle structures can form in response to the history of surface plate motions and plate subduction”

L.80: see also Giuliani et al (2021 PNAS) for an evaluation of this model in light of existing geochemical constraints for kimberlites worldwide. Not all kimberlites are related to deep-mantle structure and when they are not the isotopic features seem to be distinct from those of kimberlites that are linked to deep-mantle upwellings �  we agree, and we have made that point in the discussion (L. 924-933).

L.91: and kimberlite emplacement ages �  change made (L. 145).

L108-110. Unfortunately, the database of kimberlite geochronology by Tappe et al. (2018) has two biases, the first one being geographic areas of more intense exploration for diamonds and the second one being several ages for multiple kimberlite bodies from the same cluster or field (say area). In addition, there might be an inherent preservation bias, which cannot be addressed. I would just remove this statement.

We have deleted the statement.

L. 117: Please also check the recent Pearson et al. (2021 Nature), which include geophysical and petrological constraints �  this point has guided our revision and has resulted in significant changes to the text and to the figures

L118-120. Not all kimberlites are located on craton so using Merdith et al. (2021) to obtain craton extensions is not a good idea. Many kimberlites are located in Proterozoic mobile belts and are commonly not diamondiferous. Also, cratons are confined to continental interiors unlike the assessment by Meredith et al. shown in figure 1a.

We agree that using Merdith et al. (2021) to obtain craton extensions is not a good idea. We now consider lithosphere thicker than 150 km for the main analysis, and for a sensitivity study we also consider regions with a tectonothermal age greater than 2.5 Ga (Artemieva, 2006) and tectonic blocks inferred to have existed for at least one billion years (Merdith et al., 2021).

L.120-122: “this statement is out of context here, please move it to the previous paragraph” �  change made (L. 165-166).

L. 231: Or is this just a comparison of tomographic models? Please clarify for the uneducated reader (such as myself) �  change made (L. 372-373).

L. 334: please remind the reader the parameters employed in Case 2 - I had to go back to Table 1 to refresh my memory at this point. �  change made (L. 498-499).

L. 349-350: please specify the age of these kimberlites and/or the panels of figure 7 when these kimberlites are specifically shown �  change made (L. 516-519).

L364-365. I am not aware of any kimberlite emplaced in Western Australia at this time (~240 Ma). Not the Authors' fault, just one of the several problems with the database of Tappe et al 2018.

The Roper kimberlite field was emplaced in the Northern Territory (rather than in Western Australia). According to Hutchison (2012): “Age ranges of NT Roper field kimberlites from 65-250 Ma”. Tappe et al. (2018) list the Packsaddle-01 kimberlite as 238 Myr old (method: Rb-Sr on Phlogopite/biotite; no reference given).

L. 475-476: consistent with the likely overestimation of craton size in Meredith et al �  we have modified the text to reflect that this has been changed

L. 489: perhaps move Figure 13 to appendix �  we agree and have moved Figure 13 to the supplement (S1 Figure)

L.497-499: this is an excellent idea and something already under consideration in the diamond industry; however, a tighter definition of diamond-bearing lithospheric mantle is needed compared to the 'cratons' employed in the manuscript.

We now consider lithosphere thicker than 150 km.

it should also be considered the triggers of kimberlite magmatism because plumes are not continuously produced along LLSVPs

We agree that the mechanism linking kimberlite magmatism to LLSVPs should be investigated in the future. In this study, the link is assumed as in Torsvik et al. (2010).

L.516: also Conrad et al. (2013 Nature) �  the analysis by Conrad et al. (2013) is debated (see Rudolph and Zhong, 2013, Nature); we added a reference to Dziewonski et al. (2010) (L. 906).

L.533: see also Kjarsgaard et al. (2017 G-cubed) or Chen et al. (2020 Geology) �  we cited Chen et al. (2020 Geology) (L. 919) who explained why the model proposed by Kjarsgaard et al. (2017 G-cubed) does not work.

L. 536: Also Tappe et al. (2013 EPSL) and Tovey et al. (2021 Lithos) who both link recent kimberlite magmatism in the Slave craton to past subduction events �  we cited Tappe et al. (2013 EPSL) (L. 924)

L541-544. Additional factors to consider when modelling potential kimberlite formation include lithospheric stress, pre-existing translithospheric structures, plume buoyancy and composition just to name a few.

We mentioned that lithospheric controls on potential kimberlite formation were not considered (L. 917).

Finally, it should be kept in mind that not all kimberlites (or related magmas such as those from North China; see Tompkins et al., 1999 Proc 7th Inter Kimb Conf) are necessarily linked to the deep mantle; some might be related to shallower sources including the mantle transition zone (e.g., Kjarsgaard et al., 2017 G-cubed; Chen et al., 2020 Geology) or also the lithospheric mantle (e.g., Giuliani et al., 2015 Nature Comms). This should be noted in the manuscript in my opinion.

We have modified this paragraph, including citations to Chen et al. (2020 Geology) (L. 919) and to Giuliani et al. (2015 Nature Communications) (L. 922).

Figures 7 and 8 would be easier to follow if the whole continental blocks were shown rather than cratons only. Similarly, the statement at L393-394 is not easy to picture because continental blocks are not shown but only craton outlines.

We now show present-day coastlines as well as lithosphere thicker than 150 km on Figures 7 and 8.

Additional minor edits and suggested references are included in the attached pdf.

These have been addressed above (in italics).

I hope my comments will be helpful and my apologise for the late review (6-Feb-2022)

Thank you for these comments that have taken the manuscript to a new level in terms of referencing and have given us the opportunity to improve the definition of the extent of continental lithosphere relevant to kimberlite exploration.

Submitted filename: point-by-point-reply-to-reviewer-comments.pdf

Click here for additional data file.

22 Apr 2022

Mapping global kimberlite potential from reconstructions of mantle flow over the past billion years

PONE-D-21-39227R1

Dear Dr. Flament,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Shuan-Hong Zhang, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Dear Dr. Flament,

Thank you very much for revising the manuscript and addressed all the comments and suggestions clearly. After evaluating from the previous reviewers and my reading through of your manuscript, I am pleased to confirm that your paper has been accepted for publication in PLOS ONE.

Thank you for choosing PLOS ONE to plublish your interesting work.

Best wishes,

Shuan-Hong Zhang

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: (No Response)

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: (No Response)

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: (No Response)

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: (No Response)

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Andrea Giuliani

25 May 2022

PONE-D-21-39227R1

Mapping global kimberlite potential from reconstructions of mantle flow over the past billion years

Dear Dr. Flament:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Shuan-Hong Zhang

Academic Editor

PLOS ONE