# Assembly of the basal mantle structure beneath Africa

Nicolas Flament, Ömer F. Bodur, Simon E. Williams, Andrew S. Merdith

► To cite this version:

Nicolas Flament, Ömer F. Bodur, Simon E. Williams, Andrew S. Merdith. Assembly of the basal mantle structure beneath Africa. Nature, 2022, 603, pp.846-851. 10.1038/s41586-022-04538-y. insu-03710170

HAL Id: insu-03710170

https://insu.hal.science/insu-03710170v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assembly of the basal mantle structure beneath Africa

https://doi.org/10.1038/s41586-022-04538-y

Received: 22 April 2021

Accepted: 9 February 2022

Published online: 30 March 2022

<figure><img src="image_2.png" /></figure>

Check for updates

Nicolas Flament<fnref n="1" />,²,³, Ömer F. Bodur<fnref n="1" />, Simon E. Williams² & Andrew S. Merdith³,⁴

Plate tectonics shapes Earth's surface, and is linked to motions within its deep interior<fnref n="1" />,². Cold oceanic lithosphere sinks into the mantle, and hot mantle plumes rise from the deep Earth, leading to volcanism³,⁴. Volcanic eruptions over the past 320 million years have been linked to two large structures at the base of the mantle presently under Africa and the Pacific Ocean⁵,⁶. This has led to the hypothesis that these basal mantle structures have been stationary over geological time⁷,⁸, in contrast to observations and models suggesting that tectonic plates⁹,¹⁰, subduction zones¹¹⁻¹⁴ and mantle plumes¹⁵,¹⁶ have been mobile, and that basal mantle structures are presently deforming¹⁷,¹⁸. Here we reconstruct mantle flow from one billion years ago to the present day to show that the history of volcanism is statistically as consistent with mobile basal mantle structures as with fixed ones. In our reconstructions, cold lithosphere sank deep into the African hemisphere between 740 and 500 million years ago, and from 400 million years ago the structure beneath Africa progressively assembled, pushed by peri-Gondwana slabs, to become a coherent structure as recently as 60 million years ago. Our mantle flow models suggest that basal mantle structures are mobile, and aggregate and disperse over time, similarly to continents at Earth's surface⁹. Our models also predict the presence of continental material in the mantle beneath Africa, consistent with geochemical data¹⁹,²⁰.

Two large antipodal structures have been imaged by tomographic models at the base of the mantle under Africa and the Pacific Ocean²¹⁻²⁷. These structures, called large low-shear-velocity provinces (LLSVPs), are characterized by seismic velocities about one per cent slower than ambient mantle, indicating that they are hotter than ambient mantle<fnref n="1" />,²⁸ and possibly compositionally distinct²⁹,³⁰. At long wavelengths, these structures are about 15,000 km in diameter and cover a fraction of the area of the core–mantle boundary that is similar to that covered by continents at the surface²⁹. Analyses of tomographic models at shorter wavelengths suggest that lowermost mantle heterogeneity occurs at smaller scales³¹, and laboratory experiments suggest that thermal plumes could rise from a dense layer 300–400 km thick³². LLSVPs have been linked to volcanic eruptions at Earth's surface: the locations of large igneous provinces⁵,⁶ and kimberlites⁶ reconstructed at the time of their eruption from 320 Myr ago (Ma) are spatially close to the edge of LLSVPs. This relationship has led to the concept of LLSVP margins as plume-generation zones⁵, and to the hypothesis that LLSVPs could have been stationary over hundreds of millions of years⁷,⁸. Whether past volcanic eruption locations are linked to the margins⁵,⁶,³³ or interiors³⁴,³⁵ of LLSVPs is debated³³⁻³⁵. Regardless of this debate, mantle flow models have confirmed a link between predicted time-dependent basal mantle structures (BMSs; the model equivalent of LLSVPs) and mantle plumes¹⁶,³⁶, and shown that sinking slabs of cold oceanic lithosphere shape BMSs<fnref n="1" />,²,³⁷. This process makes the proposed fixity of BMSs difficult to reconcile with the mobility of subduction zones apparent in tectonic reconstructions¹⁰,¹³,¹⁴.

Previous work<fnref n="1" />,²,³⁷⁻³⁹ has shown that reconstructions of mantle flow back to approximately 450 Ma predict a present-day mantle structure that matches seismic tomography³⁸,³⁹. However, comparisons to tomographic models are limited to the present day, and evaluating the success of modelled mantle structure throughout geological time requires using different data and approaches. Here we extend previous statistical analyses³³⁻³⁵ to reconstructions of the structure of Earth's mantle to investigate the time dependence of the link between BMSs and surface volcanism⁴⁰,⁴¹. Specifically, we address the question of whether stationary BMSs are uniquely required to explain the distribution of volcanism since late Neoproterozoic times, or whether mobile BMSs provide an equally likely explanation.

## Modelling the deep Earth through time

To investigate the link between volcanism and BMSs (Fig. 1), we considered seven S-wave tomographic models²¹⁻²⁷, 22 mantle flow model cases (Extended Data Table 1) across which input parameters were varied, and volcanic eruptions (Fig. 1 and Extended Data Fig. 1) including kimberlites (T18, ref. ⁴¹), large igneous provinces associated with plume heads and volcanism associated with plume tails (J18, ref. ⁴⁰). As in previous work³⁸,⁴², lower mantle structure was represented by cluster analysis of seismic wave velocities and mantle temperature in tomographic models and flow models, respectively. The reconstructed locations of volcanic eruptions back to 320 Ma appear visually close to

<footnote n="1">¹GeoQUEST Research Centre, School of Earth, Atmospheric and Life Sciences, University of Wollongong, Wollongong, New South Wales, Australia. ²Department of Geology, Northwest University, Xi'an, China. ³UnivLyon, Université Lyon 1, ENS de Lyon, CNRS, UMR 5276 LGL-TPE, Villeurbanne, France. ⁴Present address: School of Earth and Environment, University of Leeds, Leeds, UK. ⁵e-mail: nflament@uow.edu.au</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Structure of the lower mantle and volcanic eruption locations from 320 Ma to the present day. **a**, High-velocity (white) and low-velocity (grey) regions revealed by cluster analysis between 1,000 km and 2,800 km depth for seismic tomographic model Savani²⁶, and location of volcanic eruptions (diamonds, J18, ref. ⁴⁰) and kimberlites (circles, T18, ref. ⁴¹) reconstructed at their time of eruption and shown at present day. **b, c, e** As in **a**, but with eruption locations reconstructed at 160 Ma (**b**) and 300 Ma (**c**). **d–f**, As in **a**, but for low-temperature (white) and high-temperature (grey) regions in a reconstruction of mantle flow (case 7), and with eruption locations</figcaption></figure>

reconstructed between 0 and 320 Ma and shown for the present day (**d**), reconstructed at 160 Ma (**e**) and reconstructed at 300 Ma (**f**). In **b, c, e, f**, locations are shown for eruptions within 10 Myr of the considered age. In **a–f**, the black lines indicate a value of five (solid) and a value of one (dotted) in a vote map for low-velocity regions in S-wave tomographic models⁴² and the brown transparent polygons are the reconstructed locations of ancient continental blocks¹⁰. Present-day coastlines are shown in brown in **a** and **b**. Symbols are coloured by age. Robinson projection at Earth's surface.

or above LLSVPs (Fig. 1a), and the spatial relationship depends on the considered reconstruction (Extended Data Fig. 2), as expected. The present-day mantle structure predicted by our flow model is consistent with tomographic models (Fig. 1d). Volcanic eruption locations from 320 Ma to the present day are spatially close to stationary LLSVPs imaged by tomographic models (Fig. 1b, c) and to mobile BMSs predicted by flow models (Fig. 1d–f).

To quantify model success, we adopted the statistical comparisons used to propose long-term fixity of LLSVPs⁶,³³. We computed the angular distance $\theta$ between each volcanic eruption location and the closest

BMS for a given tomographic or flow model (Supplementary Video 1). Distances were measured with respect to the edge of BMSs³³,³⁵, with distances positive inside BMSs, and negative outside³³ (Extended Data Fig. 3c). ‘Sample’ empirical density functions (EDFs) showed the cumulative probability that these distances, computed in 20-Myr intervals, were less than or equal to a given value (Fig. 2 and Extended Data Fig. 4). To assess whether sample EDFs were significantly better than chance, we compared them to 1,000 ‘random’ EDFs (Fig. 2 and Extended Data Fig. 4) for distances between BMSs and a set of points located anywhere on Earth’s surface with the same temporal distribution as volcanic
<figure><img src="image_4.png" /><figcaption>**Fig. 2 |** Distributions of distances between basal mantle structures and reconstructed eruption locations. **a–f**, Sample empirical distribution functions (EDFs; blue lines) showing the cumulative probability of minimum angular distances between 931 volcanic eruptions and the nearest basal mantle structure (BMS) over the last 320 Myr for three tomographic models (**a–c**) and</figcaption></figure>

three mantle flow models (**d–f**). Grey lines are a series of 1,000 EDFs showing the same quantity for points randomly distributed around the globe and BMSs. $\tilde{\theta}$ is the median distance for the sample EDF and $f_s$ is the fraction of random EDFs compared to which the sample EDF passes a statistical test.

eruptions (Extended Data Fig. 1). We analysed the time-averaged median angular distance $\tilde{\theta}$ (with cumulative probability equal to 0.5) of each sample EDF and the fraction $f_s$ of random EDFs for which (i) the mean of absolute values of distances ($|\tilde{\theta}|$) was greater than in the sample EDF and (ii) a Kolmogorov–Smirnov statistical test$^{43}$ indicated a significant difference from the sample EDF. Reconstructed volcanic eruption locations are close to BMSs ($|\tilde{\theta}| < 5.5^\circ$ in Fig. 2) for both tomographic and mantle flow models (Figs. 1, 2, Extended Data Fig. 2 and Supplementary Video 1), and $f_s \approx 1$ supports a link between the sample locations of volcanic eruptions and BMSs (Fig. 2b–f), whereas $f_s \approx 0$ indicates that the sample EDF is not significantly different from the random EDF (Fig. 2a).

## Factors controlling model success

To evaluate model success and trade-offs between parameters, we analysed four quantities for each mantle structure (Figs. 2, 3 and Extended Data Figs. 4–7): (i) the fractional area $f_a$ of deep mantle covered by BMSs in cluster space (time-averaged for flow models), (ii) the fit $\bar{\text{Acc}}$ to present-day mantle structure averaged across independent tomographic models from $\bar{\text{Acc}} = (\text{TP} + \text{TN})/A$, where TP is the area of true positives (orange in Extended Data Fig. 3), TN the area of true negatives (grey in Extended Data Fig. 3) and A the total area, (iii) $\tilde{\theta}$ and (iv) $f_s$. We identified $f_a$ as a key factor affecting $\bar{\text{Acc}}$ and $\tilde{\theta}$. For tomographic models, $0.32 < f_a < 0.52$ and $0.74 < \bar{\text{Acc}} < 0.87$; $\bar{\text{Acc}}$ was smallest and $\tilde{\theta}$ was positive for reconstruction M21-NNR (NNR, no net rotation) when $f_a$ was largest (GyPSuM-S). For SEMUCB-WM1, a relatively large $f_a$ and more numerous BMSs (Extended Data Fig. 3) resulted in smaller $|\tilde{\theta}|$ in random EDFs than in the sample EDF, leading to $f_s = 0$ (Figs. 2, 3); otherwise, $f_s \ge 0.9$ (Fig. 3d).

We assessed the effect of the tectonic reconstruction on results (Extended Data Figs. 2, 5), by considering reconstructions M16 (ref. 44) and Y19 (ref. 45), both extending back to 410 Ma, and M21 (ref. 10), extending back to 1,000 Ma. We considered two reference frames for reconstruction of M21: the original palaeomagnetic reference frame$^{10}$

(reconstruction M21), and a reference frame in which lithospheric net rotation—the motion of the entire lithosphere with respect to the mantle—was removed because it induces wholesale motion of the lower mantle in palaeogeographically constrained mantle flow models$^{46}$ (reconstruction M21-NNR). We obtained $|\tilde{\theta}| < 5^\circ$ across all tomographic models and reconstructions (Extended Data Fig. 5), even though M16 was the only reconstruction with inherent circularity in matching large igneous province locations to LLSVPs assumed to have been stationary since Palaeozoic times$^{12,44,47}$, which should lead to low $|\tilde{\theta}|$. Reconstructed volcanic eruption locations were furthest within BMSs for greater $f_a$ (GyPSuM-S); however, large $f_a$ may lead to $f_s = 0$ for some reconstructions (SEMUCB-WM1 and GyPSuM-S; Extended Data Fig. 5). We focus on mantle flow models based on reconstruction M21-NNR in the following because $\bar{\text{Acc}}$ was significantly larger for case 2 based on reconstruction M21-NNR than for case 12 based on reconstruction M21 and otherwise identical to case 2 (Extended Data Fig. 5b).

In the flow models, $f_a$ depends on the intrinsic density $\delta\rho_b$ and viscosity $\eta_b$ of the basal layer. We first determined that a basal layer ten times more viscous than ambient mantle ($\eta_b = 10 \times \eta_0$) led to greater $\bar{\text{Acc}}$ than $\eta_b = \eta_0$ or $\eta_b = 0.1 \times \eta_0$ (Fig. 3b). We recognized a trade-off between $f_a$, $\tilde{\theta}$ and $\bar{\text{Acc}}$ by varying $\delta\rho_b$ in models with $\eta_b = 10 \times \eta_0 \cdot f_a$ was consistent with tomographic models if $\delta\rho_b \ge +1.3\%$ (Fig. 3a). $\bar{\text{Acc}}$ reached up to 0.74 (best fit to present-day lower mantle structure) for $\delta\rho_b = +1.0\%$, which led to the smallest $f_a$ (case 6, Fig. 3a, b, Extended Data Fig. 3a and Supplementary Video 2), and $\bar{\text{Acc}}$ was smallest (worst fit to present-day lower mantle structure) for $\delta\rho_b = +0\%$ (case 4). $|\tilde{\theta}|$ was largest (poor fit to volcanic history) for $\delta\rho_b = +0.5\%$ (case 5) and decreased towards $\delta\rho_b = +2.0\%$ (case 9), owing to an increase in $f_a$ (Fig. 3a, c), and towards $\delta\rho_b = +0\%$ (case 4), which reflects an increase in the number of BMSs (Extended Data Fig. 3a). $f_s = 1$ for cases 1–9 (Fig. 3d). No single case produced metrics that simultaneously fall within the observed range. Our preferred case is case 7 ($\delta\rho_b = +1.3\%$), for which $\bar{\text{Acc}}$ is larger than for most other flow model cases, and $|\tilde{\theta}|$ and $f_a$ are within the range of results for tomographic models (Fig. 3).
<figure><img src="image_5.png" /><figcaption>**Fig. 3 | Match of model basal mantle structures to volcanic eruption locations and tomographic models.** a, Fractional area *f*<sub>s</sub> of deep mantle cluster maps covered by slow (in tomography) or hot (in flow models, averaged over 320 Myr) BMSs. b, Spatial match *Acc* between present-day BMSs for a given case and LLSVPs imaged by tomographic models. c, Time-averaged median of minimum angular distances *θ̂* between BMSs and volcanic eruptions from 320 Ma to the present day. d, Fraction *f*<sub>s</sub> of random EDFs compared to which the</figcaption></figure>

sample EDF passes a statistical test. In a–d, the first row shows results for a series of tomographic models (stationary LLSVPs), and the last two rows show results for two series of flow models (mobile BMSs) across which parameters are varied as indicated. The horizontal lines denote reference cases. The grey shadings in a–c highlight the range of results for tomographic models. C1–C9, mantle flow model cases 1–9; T1–T7, tomographic models 1–7.

In a second set of experiments, BMSs were solid, with distances to volcanic eruptions equal to zero within BMSs and positive outside BMSs (Extended Data Fig. 3c). These experiments confirmed the trends described above, and that case 7 was most consistent with tomographic models (*f*<sub>a</sub> and *Acc*) and reconstructed locations of volcanic eruptions (*θ̂* and *f*<sub>s</sub>; Extended Data Fig. 6a–d). A third series of experiments back to 640 Ma was possible because of the long period covered by T18 (Extended Data Fig. 1b) and by our mantle flow models. In these experiments, *f*<sub>s</sub> = 1 for flow models based on reconstruction M21-NNR, and |*θ̂*| was similar for mobile BMSs predicted by cases 4–9 and for stationary LLSVPs imaged by tomographic models, whereas |*θ̂*| was larger and *f*<sub>s</sub> was smaller for case 12 based on reconstruction M21 (Extended Data Fig. 6e–h). We investigated the sensitivity of results to the initial conditions and presence of mineral phase changes. In model cases with δρ<sub>b</sub> = +2.0%, *f*<sub>s</sub> > 0.85 for cases starting at 1,000 Ma, 900 Ma, 500 Ma or 400 Ma (Extended Data Fig. 7a). The initial slab depth and phase changes at 410 km and 670 km depth had a minor influence on the location of BMSs (Supplementary Video 1), with initially deeper slabs (compare case 2 to case 6 in Fig. 3) and no phase change giving favourable results across *f*<sub>a</sub>, *Acc*, |*θ̂*| and *f*<sub>s</sub> (Extended Data Fig. 7a–d).

We found that mobile BMSs predicted by palaeogeographically constrained mantle flow models are consistent with the history of volcanism. In the preferred model case 7, the basal layer is ten times more viscous and 1.3% denser than ambient mantle. This intrinsic density is consistent with that inferred from sharp gradients in seismic velocity along the eastern edge of the African LLSVP³⁰, although it is 2.6 times as large as that inferred from tidal tomography for LLSVPs globally⁴⁸. Increasing the bulk modulus of the basal layer could stabilize BMSs⁴⁹

and decrease the BMS density required to best fit observations. Converting the mantle structure from temperature and pressure to seismic velocities, and applying a seismic filter²⁵ could also reduce the BMS density that best fits tomographic models²⁸. Seismic filtering, which does not affect *Acc* (ref. ³⁸), is limited to the present day and cannot be considered in our time-dependent analysis.

## Deep slabs link BMSs to palaeogeography

Interestingly, our reconstructions of past mantle flow starting at 800 Ma or earlier all predict that the African BMS was not yet a coherent structure at 160 Ma (Fig. 1e and Supplementary Video 1). This result contrasts with previous palaeogeographically constrained mantle flow models that suggested that the African BMS could have been formed by approximately 200 Ma³⁷,⁴⁵. A key difference is that the tectonic reconstructions used here extend much further back in time and capture key geological events and evolving plate boundary configurations that have shaped the structure of the African mantle domain. In flow models based on reconstruction M21-NNR, cold slabs sank into the mantle from approximately 740 Ma in the African hemisphere along a network of subduction zones around western Gondwanan tectonic blocks (Extended Data Fig. 8a). These slabs progressively sank into the mantle and accumulated above the core–mantle boundary by 500 Ma (Fig. 4a). The African BMS then aggregated as slabs sinking from 400 Ma along the northern margin of the Palaeo-Tethys Ocean and the western margin of Gondwana (Extended Data Fig. 8c) pushed three distinct BMSs together (Fig. 4b, Supplementary Videos 1, 3–5). The aggregation of these BMSs warmed up the cold material subducted from 740 Ma that
<figure><img src="image_6.png" /></figure>

**Fig. 4 | Dynamic assembly of the African basal mantle structure and predicted continental crust in the deep Earth.** Mantle structures 310 K hotter than ambient mantle below 170 km depth (coloured by depth with warm colours), mantle structures 155 K colder than ambient mantle below 1,500 km depth (coloured by depth with cold colours), and motion paths of sinking oceanic lithosphere below 250 km depth (coloured by depth with cold colours). **a**, From 740 Ma around western Gondwanan tectonic blocks (shown at 500 Ma and centred on 30° E). **b**, From 400 Ma from the northern margin of the Palaeo-Tethys Ocean and the western margin of Gondwana (shown at 80 Ma

and centred on 60° E). **c**, **d**, Deep mantle structures (coloured by depth with warm colours) and mantle predicted to presently consist of at least 5% continental crust (below 170 km depth, coloured by depth with cold colours) for the African mantle domain (**c**) and Pacific mantle domain (**d**). Ancient continental blocks are shown in pink, and coastlines are shown in white. In **a–d**, reconstructed plate boundaries are shown in green, the graticule spacing is 30°, and the white arrow indicates north. WGB, west Gondwanan blocks; WGS, west Gondwanan slabs.

separated them, leading to the formation of a coherent African BMS by approximately 60 Ma (Supplementary Videos 1, 3–5). The evolution of the mobile African BMS is consistent with the history of volcanism in the African hemisphere: $f_s = 1$ for cases 1–9, even though $f_s < 0.8$ for four tomographic models (Extended Data Fig. 7e–h).

As previously observed in tomographic and flow models², at long wavelengths the model Pacific BMS is more rounded than the model African BMS, which is ridge-like (Fig. 4c, d). Both BMSs present smaller-scale heterogeneity³¹ with approximately 15 long-lived lower mantle thermal mantle plumes rising from topographic highs¹⁶ (Fig. 4 and Supplementary Videos 3–7) above a relatively thin dense layer³². By tracking the evolution of numerical tracers representing compositionally distinct continental crust in the models, we found that model sinking slabs entrained continental crust down to the African BMS, but not to the Pacific BMS (Supplementary Videos 6, 7), because subduction zones have mostly been adjacent to continents in the African hemisphere over the last billion years (Extended Data Fig. 8). This first-order result is consistent with a recent reinterpretation of deep mantle domains based on the isotopic geochemical signatures (Pb, Nd and Sr) of plume-related basalts¹⁹. However, such analyses depend on plume selection²⁰, and continental crust could have been entrained to the Pacific hemisphere earlier in Earth's history²⁰, before one billion

years ago. Our reconstructions of past mantle flow show that the structure of the deep Earth is linked to the evolution of its surface, and that BMSs assemble and disperse over time in a way that is reminiscent of the past motion of continents.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04538-y.

1. Bunge, H.-P. et al. Time scales and heterogeneous structure in geodynamic Earth models. *Science* **280**, 91–95 (1998).

2. McNamara, A. K. & Zhong, S. Thermochemical structures beneath Africa and the Pacific Ocean. *Nature* **437**, 1136–1139 (2005).

3. Bercovici, D., Schubert, G. & Glatzmaier, G. A. Three-dimensional spherical models of convection in the Earth's mantle. *Science* **244**, 950–955 (1989).

4. Campbell, I. H. Large igneous provinces and the mantle plume hypothesis. *Elements* **1**, 265–269 (2005).

5. Burke, K., Steinberger, B., Torsvik, T. H. & Smethurst, M. A. Plume generation zones at the margins of large low shear velocity provinces on the core–mantle boundary. *Earth Planet. Sci. Lett.* **265**, 49–60 (2008).
6. Torsvik, T. H., Burke, K., Steinberger, B., Webb, S. J. & Ashwal, L. D. Diamonds sampled by plumes from the core-mantle boundary. *Nature* **466**, 352–355 (2010).

7. Conrad, C. P., Steinberger, B. & Torsvik, T. H. Stability of active mantle upwelling revealed by net characteristics of plate tectonics. *Nature* **498**, 479–482 (2013).

8. Dziewonski, A. M., Lekic, V. & Romanowicz, B. A. Mantle anchor structure: an argument for bottom-up tectonics. *Earth Planet. Sci. Lett.* **299**, 69–79 (2010).

9. Irving, E. Drift of the major continental blocks since the Devonian. *Nature* **270**, 304–309 (1977).

10. Merdith, A. S. et al. Extending full-plate tectonic models into deep time: linking the neoproterozoic and the phanerozoic. *Earth-Sci. Rev.* **214**, 103477 (2021).

11. Moresi, L., Betts, P. G., Miller, M. S. & Cayley, R. A. Dynamics of continental accretion. *Nature* **508**, 245–248 (2014).

12. Domeier, M. & Torsvik, T. H. Plate tectonics in the late Paleozoic. *Geosci. Front.* **5**, 303–350 (2014).

13. Flament, N., Williams, S., Müller, R. D., Gurnis, M. & Bower, D. J. Correspondence: Reply to 'Numerical modelling of the PERM anomaly and the Emeishan large igneous province'. *Nat. Commun.* **8**, 822 (2017).

14. Li, Z. X. et al. Decoding Earth's rhythms: modulation of supercontinent cycles by longer superocean episodes. *Precambrian Res.* **323**, 1–5 (2019).

15. Tarduno, J., Bunge, H.-P., Sleep, N. & Hansen, U. The bent Hawaiian–Emperor hotspot track: inheriting the mantle wind. *Science* **324**, 50–53 (2009).

16. Hassan, R., Müller, R. D., Gurnis, M., Williams, S. E. & Flament, N. A rapid burst in hotspot motion through the interaction of tectonics and deep mantle flow. *Nature* **533**, 239–242 (2016).

17. Frost, D. A. & Rost, S. The P-wave boundary of the large-low shear velocity province beneath the Pacific. *Earth Planet. Sci. Lett.* **403**, 380–392 (2014).

18. Lynner, C. & Long, M. D. Lowermost mantle anisotropy and deformation along the boundary of the African LLSPV. *Geophys. Res. Lett.* **41**, 3447–3454 (2014).

19. Doucet, L. S. et al. Distinct formation history for deep-mantle domains reflected in geochemical differences. *Nat. Geosci.* **13**, 511–515 (2020).

20. Jackson, M., Becker, T. & Steinberger, B. Spatial characteristics of recycled and primordial reservoirs in the deep mantle. *Geochem. Geophys. Geosyst.* **22**, e2020GC009525 (2021).

21. Mégnin, C. & Romanowicz, B. A. The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms. *Geophys. J. Int.* **143**, 709–728 (2000).

22. Houser, C., Masters, G., Shearer, P. & Laske, G. Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms. *Geophys. J. Int.* **174**, 195–212 (2008).

23. Kustowski, B., Ekström, G. & Dziewoński, A. M. Anisotropic shear-wave velocity structure of the Earth's mantle: a global model. *J. Geophys. Res. Solid Earth* **113**, B06306 (2008).

24. Simmons, N. A., Forte, A. M., Boschi, L. & Grand, S. P. GyPSuM: a joint tomographic model of mantle density and seismic wave speeds. *J. Geophys. Res. Solid Earth* **115**, B12310 (2010).

25. Ritsema, J., Deuss, A., van Heijst, H. J. & Woodhouse, J. H. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. *Geophys. J. Int.* **184**, 1223–1236 (2011).

26. Auer, L., Boschi, L., Becker, T. W., Nissen-Meyer, T. & Giardini, D. Savani: a variable resolution whole-mantle model of anisotropic shear velocity variations based on multiple data sets. *J. Geophys. Res. Solid Earth* **119**, 3006–3034 (2014).

27. French, S. W. & Romanowicz, B. A. Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography. *Geophys. J. Int.* **199**, 1303–1327 (2014).

28. Davies, D., Goes, S. & Lau, H. C. P. In *The Earth's Heterogeneous Mantle* (eds Khan, A. & Deschamps, F.) 441–477 (Springer, 2015).

29. Garnero, E. J., McNamara, A. K. & Shim, S.-H. Continent-sized anomalous zones with low seismic velocity at the base of Earth's mantle. *Nat. Geosci.* **9**, 481–489 (2016).

30. Ni, S., Tan, E., Gurnis, M. & Helmberger, D. V. Sharp sides to the African superplume. *Science* **296**, 1850–1852 (2002).

31. Tkalčić, H., Young, M., Muir, J. B., Davies, D. R. & Mattesini, M. Strong, multi-scale heterogeneity in Earth's lowermost mantle. *Sci. Rep.* **5**, 18416 (2015).

32. Davaille, A. & Romanowicz, B. Deflating the LLSPVs: bundles of mantle thermochemical plumes rather than thick stagnant “piles”. *Tectonics* **39**, e2020TC006265 (2020).

33. Doubrovine, P. V., Steinberger, B. & Torsvik, T. H. A failure to reject: testing the correlation between large igneous provinces and deep mantle structures with EDF statistics. *Geochem. Geophys. Geosyst.* **17**, 1130–1163 (2016).

34. Austermann, J., Kaye, B. T., Mitrovica, J. X. & Huybers, P. A statistical analysis of the correlation between large igneous provinces and lower mantle seismic structure. *Geophys. J. Int.* **197**, 1–9 (2014).

35. Davies, D., Goes, S. & Sambridge, M. On the relationship between volcanic hotspot locations, the reconstructed eruption sites of large igneous provinces and deep mantle seismic structure. *Earth Planet. Sci. Lett.* **411**, 121–130 (2015).

36. Garnero, E. J. & McNamara, A. K. Structure and dynamics of Earth's lower mantle. *Science* **320**, 626–628 (2008).

37. Zhong, S. & Rudolph, M. L. On the temporal evolution of long-wavelength mantle structure of the Earth since the early Paleozoic. *Geochem. Geophys. Geosyst.* **16**, 1599–1615 (2015).

38. Flament, N., Williams, S., Müller, R., Gurnis, M. & Bower, D. J. Origin and evolution of the deep thermochemical structure beneath Eurasia. *Nat. Commun.* **8**, 14164 (2017).

39. Flament, N. Present-day dynamic topography and lower-mantle structure from palaeogeographically constrained mantle flow models. *Geophys. J. Int.* **216**, 2158–2182 (2019).

40. Johansson, L., Zahirovic, S. & Müller, R. D. The interplay between the eruption and weathering of large igneous provinces and the deep-time carbon cycle. *Geophys. Res. Lett.* **45**, 5380–5389 (2018).

41. Tappe, S., Smart, K., Torsvik, T., Massuyeau, M. & de Wit, M. Geodynamics of kimberlites on a cooling Earth: clues to plate tectonic evolution and deep volatile cycles. *Earth Planet. Sci. Lett.* **484**, 1–14 (2018).

42. Lekic, V., Cottaar, S., Dziewonski, A. & Romanowicz, B. A. Cluster analysis of global lower mantle tomography: a new class of structure and implications for chemical heterogeneity. *Earth Planet. Sci. Lett.* **357**, 68–77 (2012).

43. Kolmogorov, A. Sulla determinazione empirica di una Igge di distribuzione. *Gior. Inst. Ital. Attuari* **4**, 83–91 (1933).

44. Matthews, K. J. et al. Global plate boundary evolution and kinematics since the late Paleozoic. *Global Planet. Change* **146**, 226–250 (2016).

45. Young, A. et al. Global kinematics of tectonic plates and subduction zones since the late Paleozoic Era. *Geosci. Front.* **10**, 989–1013 (2019).

46. Rudolph, M. L. & Zhong, S. History and dynamics of net rotation of the mantle and lithosphere. *Geochem. Geophys. Geosyst.* **15**, 3645–3657 (2014).

47. Torsvik, T. H. et al. Deep mantle structure as a reference frame for movements in and on the Earth. *Proc. Natl Acad. Sci.* **111**, 8735–8740 (2014).

48. Lau, H. C. et al. Tidal tomography constrains Earth's deep-mantle buoyancy. *Nature* **551**, 321–326 (2017).

49. Bower, D. J., Gurnis, M. & Seton, M. Lower mantle structure from paleogeographically constrained dynamic Earth models. *Geochem. Geophys. Geosyst.* **14**, 44–63 (2013).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2022, corrected publication 2022
Mantle flow models

The evolution of mantle flow was reconstructed from one billion years ago using the code CitcomS⁵⁰, in which the mantle was represented with finite elements in a spherical domain, with velocities of tectonic plates from tectonic reconstructions with continuously closing plate polygons⁵¹ imposed as boundary conditions. We represented the mantle with 129 × 129 × 65 × 12 ≈ 13 million nodes which, with radial mesh refinement, resulted in average resolutions of roughly 50 km × 50 km × 15 km at the surface, roughly 40 km × 40 km × 100 km in the mid-mantle, and roughly 28 km × 28 km × 27 km at the core–mantle boundary. We used a version of CitcomS⁵² modified to progressively assimilate the age of the ocean floor and the location and polarity of subduction zones determined in one-million-year intervals. This semi-empirical approach makes it possible to investigate the time dependence of the plate–mantle system with one-sided subduction. This problem remains intractable in fully dynamic models that are either limited to the present day⁵³ or result in two-sided intra-oceanic subduction⁵⁴. Here we give an overview of the parameters and set-up of the model, which are further described elsewhere³⁹,⁵⁵.

We approximated Earth's mantle as an incompressible fluid under the extended Boussinesq approximation, which accounts for viscous dissipation, an adiabatic temperature gradient, internal heating and the depth dependence of the coefficient of thermal expansion⁵⁶. In the models, the vigour of convection was controlled by the Rayleigh number $Ra = (\alpha_0 \rho_0 g_0 \Delta T h_M^3) / (\kappa_0 \eta_0)$, where $\alpha_0$ is the coefficient of thermal expansion, $\rho_0 = 4,000$ kg m⁻³ is the density, $g_0 = 9.81$ m s⁻² is the gravitational acceleration, $\Delta T = 3,100$ K is the temperature change across the mantle, $h_M = 2,867$ km is the thickness of the mantle, $\kappa_0 = 1 \times 10^{-6}$ m² s⁻¹ is the thermal diffusivity, $\eta_0 = 1.1 \times 10^{21}$ Pa s is the viscosity, and the subscript '0' indicates reference values. The coefficient of thermal expansion decreased with increasing depth by a factor of approximately two over the thickness of the mantle (see refs. ⁵⁷,⁵⁸) from $\alpha_0 = 3 \times 10^{-5}$ K⁻¹ at the surface. With the values listed above, $Ra = 7.8 \times 10^7$. Heat dissipation by shear heating was controlled by the dissipation number $Di = (\alpha_0 g_0 R_0) / C_{p0}$, where $R_0 = 6,371$ km is the radius of Earth and $C_{p0} = 1,200$ J kg⁻¹ K⁻¹ is the reference heat capacity. The rate of internal heating was $H = 33.6$ TW, which represents heat produced by radioactive elements (between 9–17 TW)⁵⁹ and primordial sources of heat dissipated by secular cooling of the solid Earth (between 8–29 TW)⁵⁹. The predicted present-day surface heat flux (excluding radioactive sources in the crust) ranged between -36.5 TW and -38.5 TW, and the predicted present-day core–mantle boundary heat flux ranged between -8 TW and -13.5 TW, both of which are consistent with constraints⁵⁹. For the preferred case 7, the present-day surface heat flux was equal to 38.6 TW, and the predicted present-day core–mantle boundary heat flux was equal to 10.1 TW.

Viscosity varied with depth, composition, temperature, and pressure following

$$ \eta = \eta(r) \eta_0 \eta_c \exp \left\{ \frac{[E_{\eta} + \rho_0 g_0 Z_{\eta} (R_0 - r)]}{[R(T + T_{\text{off}})]} - \frac{[E_{\eta} + \rho_0 g_0 Z_{\eta} (R_0 - R_c)]}{[R(T_{\text{CMB}} + T_{\text{off}})]} \right\}, \quad (1) $$

where $\eta(r) = 0.02$ above 160 km depth and between 310–660 km depth, $\eta(r) = 0.002$ between 160–310 km depth (asthenosphere), and $\eta(r) = 0.02$ below 670 km depth (lower mantle). The compositional viscosity pre-factor was $\eta_c = 100$ for the continental lithosphere (modelled as in ref. ⁵⁵) and $\eta_c = \eta_b$ for the basal layer was varied across model cases (see Fig. 3 and Extended Data Figs 6, 7 and Extended Data Table 1). $r$ is the radius, $R_c = 3,504$ km is the radius of the core, $E_{\eta} = 284$ kJ mol⁻¹ is the activation energy, $Z_{\eta} = 2.1 \times 10^{-6}$ m³ mol⁻¹ is the activation volume, $R = 8.31$ J mol⁻¹ K⁻¹ is the universal gas constant, $T$ is the dimensional temperature, $T_{\text{off}} = 452$ K is a temperature offset

and $T_{\text{CMB}} = 3,380$ 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 \times 10^{20}$ Pa s to $2.2 \times 10^{23}$ Pa s) across the range of temperatures and pressures (Extended Data Fig. 9c). Such lateral viscosity contrasts are smaller than expected to occur within the solid Earth⁵³, however they can be computed with a resolution that makes it possible to model time-dependent mantle flow. The parameters used in equation (1) result in horizontally averaged present-day mantle viscosity profiles that are consistent with a viscosity profile adjusted to fit the geoid and post-glacial rebound⁶⁰ (Extended Data Fig. 9d). Phase changes at 410 km and 670 km depth were considered in model cases 22 and 23 (Extended Data Table 1), assuming a phase change to occur over a depth range of 40 km. At 410 km depth, the Clapeyron slope was taken equal to 4 MPa K⁻¹, and the phase change resulted in a 3% increase in density (ref. ⁶¹ and references therein). At 670 km depth, the Clapeyron slope was taken equal to −2 MPa K⁻¹, and the phase change resulted in a 7% increase in density (ref. ⁶¹ and references therein).

In the initial condition and for progressive data assimilation, the thickness and temperature of the lithosphere were derived using a half-space cooling model taking 80 Myr as the maximum age of ocean floor (there was no maximum age for case 10). Synthetic past ages of the ocean floor (Extended Data Fig. 8) were computed using the method of ref. ⁶², and simplified tectonothermal ages for the continental lithosphere as in ref. ⁵⁵. The global thermal structure of slabs was constructed from the location of subduction zones and from the age of the ocean floor⁵². The global thermal structure of the lithosphere and of subducting slabs was assimilated in the dynamic models in 1-Myr increments, down to 350 km depth at subduction zones⁵². In the initial condition (Extended Data Fig. 9a, b), subduction zones were inserted into the mantle assuming a descent rate of 3 cm yr⁻¹ in the upper mantle and 1.2 cm yr⁻¹ in the lower mantle⁶³. Subduction zones that appeared during the model were progressively inserted to 350 km depth based on the age of subduction initiation and on the plate convergence rate.

The initial condition at age $a_0$ (which was varied across cases 14–20, Extended Data Table 1 and Extended Data Fig. 7) consisted of an adiabatic temperature profile between two thermal boundary layers (Extended Data Fig. 9a, b). The basal thermal boundary layer was initially laterally uniform and 225-km thick and included a 113-km-thick layer (2% volume of Earth’s mantle)⁶⁴ of material that was denser than ambient mantle (Extended Data Fig. 9a). The composition of the continental lithosphere and of the basal layer was modelled using tracers as in ref. ⁵⁵. The density contrast between the basal material and ambient mantle $\delta\rho_b$ was deduced from the buoyancy ratio $B = \delta\rho_b / (\rho a \Delta T)$, where $\rho = 5,546$ kg m⁻³ and $a = 1.32 \times 10^{-5}$ K⁻¹ are the average density⁶⁵ and thermal expansivity⁵⁷ in the 100-km layer above the core–mantle boundary. Varying $B$ between 0 and 0.5 resulted in $0 < \delta\rho_b < +2.0\%$ across cases 4–9 (Fig. 3, Extended Data Table 1 and Extended Data Figs. 6, 7). The continental crust was 7.7% lighter than ambient mantle ($B = -0.83$ with $\rho = 3,380$ kg m⁻³ as in ref. ⁶⁵ and $\alpha_0 = 3 \times 10^{-5}$ K⁻¹ close to the surface). Slabs were initially inserted down to a depth of either $z_0$ (varied between 550 km and 1,200 km across cases 2, 6 and 21; Extended Data Table 1 and Extended Data Fig. 7), or to the depth derived from their initiation age and sinking rates if that depth was shallower than $z_0$, with a dip of 45° down to 425 km and a dip of 90° below 425 km depth. Slabs were initially twice as thick in the lower mantle compared to their thickness in the upper mantle, to account for advective thickening in the more viscous lower mantle⁶⁶. Motion paths of sinking slabs were computed using the Particle Trace filter in ParaView⁶⁷ to advect seed locations within slabs at -250 km depth at the initial age of interest (either 740 Ma in Fig. 4a or 400 Ma in Fig. 4b) according to predicted velocities.

We considered four global tectonic reconstructions ('TR' in Extended Data Table 1). Reconstruction M16, used for case 10 (Extended Data Table 1 and Extended Data Figs. 2, 5), is the reconstruction of ref. ⁴⁴, which
is based on the reconstructions of ref. ¹² for the period 410–250 Ma, with a palaeomagnetic absolute reference frame corrected for true polar wander⁴⁷, and of ref. ⁶⁸ for the period 230–0 Ma, although with absolute plate motions based on a different true polar wander-corrected reference frame⁶⁹ for the period 230–70 Ma. The reference frame is that of ref. ⁷⁰ for the last 70 Myr. Reconstruction Y19, used for case 11 (Extended Data Table 1 and Extended Data Figs. 2, 5), is the reconstruction of ref. ⁴⁵, which adopts different scenarios than reconstruction M16 for the closure of the Rheic Ocean and motion of circum-Paleaeotethys blocks, and the palaeomagnetic reference frame of ref. ⁷¹ for the period 410–250 Ma, to decrease global tectonic speeds and trench migration rates. Reconstructions M16 and Y19 are essentially identical for the period 250–0 Ma. Reconstruction M21 (ref. ¹⁰) links the reconstruction by ref. ⁷² for the period 1,000–500 Ma, to reconstructions from refs. ⁷³,⁷⁴ for the period 500–410 Ma and to reconstruction Y19 for the last 410 Myr. We considered two endmembers of that reconstruction. The first endmember (M21, used for case C11, Extended Data Table 1 and Extended Data Figs. 2, 5) is based on a palaeomagnetic reference frame, including a new apparent polar wander path for Gondwana from 540 to 320 Ma and a global apparent polar wander path from 320 Ma to 0 Ma (ref. ¹⁰). In the second endmember (M21-NNR, used in cases C1–9 and C12–23) lithospheric net rotation, which is the motion of the entire lithosphere with respect to the mantle, was removed because net rotation is known to induce rotation of the lower mantle⁴⁶, and because removing net rotation from the tectonic reconstruction greatly improves the fit to present-day mantle structure as imaged by tomographic models (Extended Data Fig. 5b).

Model results, including the predicted temperature of the mantle and the composition field, were analysed in 20-million-year increments.

## Cluster analysis of lower mantle structure

As in ref. ³⁸, we used cluster analysis to objectively classify a set of 196,596 points equally spaced on Earth's surface (average distance -0.45°) into groups of points with similar variations in a given property with depth. For each mantle flow model case (or tomographic model), time-dependent temperature (or present-day seismic velocity) profiles were treated as independent vectors of 31 coordinates specifying the temperatures (or seismic velocities) sampled at depths between 1,000–2,800 km (as in refs. ³⁸,⁴²), with an average spacing of 58 km. The vectors were grouped into two clusters using *k*-means clustering⁷⁵, a procedure that minimizes the variance in squared Euclidean distance between vectors within each cluster. We used the scientific Python implementation of the *k*-means algorithm (http://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.vq.kmeans2.html).

## Spatial match between present-day basal structures for a given case and basal structures imaged by tomographic models

We quantified how well clusters obtained for each case (mantle flow model or tomographic model) matched the global geographic distribution of tomographic models by computing the global accuracy $Acc = (TP + TN)/A$, where TP is the area of true positives (orange in Extended Data Fig. 3 and Supplementary Video 2), TN is the area of true negatives (grey in Extended Data Fig. 3 and Supplementary Video 2) and A is the total area. The accuracy was computed for each model case and for each tomographic model, against seven global S-wave tomographic models: SEMUCB-WM1²⁷, Savani²⁶, S4ORTS²⁵, GyPSuM-S²⁴, S362ANI²³, HMSL-S²² and SAW24B16²¹. For mantle flow model cases, the accuracy Acc was averaged across all seven tomographic models, and for tomographic models, the accuracy Acc was averaged across the six independent tomographic models. Even though tomographic models are only available for the present day, we quantified Acc back to 640 Ma as an indication of the change in lower mantle structure over time (Supplementary Video 2).

## Fractional area covered by hot (or slow) basal mantle structures

We computed the fraction of the global geographic distribution of cluster maps covered by hot (for mantle flow models) or slow

(for tomographic models) BMSs. The fractional area $f_a$ was reported for the present day for tomographic models and averaged over the period of interest (either from 320 Ma or from 640 Ma) for mantle flow models.

## Distance between volcanic eruption locations and basal mantle structures

We considered volcanic eruptions from ref. ⁴⁰, herein referred to as J18, and ref. ⁴¹, herein referred to as T18. J18 contains large igneous provinces (LIPs) from refs. ⁷⁶⁻⁷⁸ back to ~400 Ma. These LIPs are typically associated with plume heads⁷⁹, covering >100,000 km² and erupted as igneous pulses 1–5 Myr long over less than 50 Myr (ref. ⁷⁶). J18 also contains smaller oceanic islands and seamounts associated with plume tails⁷⁷, which are more persistent than transient plume heads⁷⁹. T18 contains kimberlites from all major known kimberlite clusters, without over-representing diamondiferous clusters. We used pyGPlates (www.gplates.org/pygplates), an application programming interface (API) for the plate tectonic reconstruction software GPlates⁸⁰, to identify the centroid of LIP polygons in database J18, only considering LIPs covering an area greater than 20,000 km². We adjusted the temporal resolution of databases J18 and T18 to that of the model output (20 Myr) as follows: volcanic eruptions that occurred within $(a - a_w) < a \le (a + a_w)$, where $a_w = 10$ Myr is an age window, were assigned the age $a$. This temporal resolution is appropriate given that LIPs may be emplaced over up to 50 Myr (ref. ⁷⁶). We combined centroids from J18 with kimberlites from T18 to carry out the presented analyses.

We created maps of angular distances from either the edge of open BMSs (with distances positive inside and negative outside BMSs: Extended Data Fig. 3c; refs. ³³,³⁵) or the edge of solid BMSs (with distances equal to 0 above BMSs, and positive outside BMSs; Extended Data Fig. 3c). For a given tectonic reconstruction, minimum angular distances to BMSs were queried from these maps at reconstructed volcanic eruption locations, or at random locations for the statistical analysis of results. These minimum distances were cumulated over the period of interest and represented as empirical distribution functions (EDFs) for volcanic eruption locations (sample EDFs) or for random points (random EDFs). The median $\tilde{\theta}$ of sample EDFs was reported (Fig. 3 and Extended Data Figs. 5–7).

## Statistical significance of results

For each mantle flow or tomographic model case, we created 1,000 distributions of spatially random points. Each distribution consisted of exactly as many points as are present in the database of volcanic eruption locations over the period of interest, and random points were distributed over time as in the database (Extended Data Fig. 1a, b). For the case limited to the African hemisphere (Extended Data Fig. 7e–h), random points were selected from a global population large enough to obtain as many points as were present within the African hemisphere in the database of volcanic eruptions. For each random EDF, we carried out a Kolmogorov–Smirnov statistics test⁴³ based on the distance between two EDFs. We used the scientific Python implementation of the Kolmogorov–Smirnov test for two samples (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats ks_2 samp.html). In these tests, the null hypothesis that the sample EDF and the mean of random EDFs were accidentally selected from the same underlying population could be rejected if the Kolmogorov–Smirnov statistic *D* was greater than a critical value, a relationship which we expressed as⁸¹: $D_n > \sqrt{-\ln(\alpha/2)/n}$, where *n* was the number of LIPs in both the sample and random EDFs, and $\alpha = 0.05$ was the confidence level. In addition to this Kolmogorov–Smirnov test, we computed the mean of absolute values of the minimum angular distance ($|\tilde{\theta}|$) between volcanic eruption locations and BMSs for each distribution and compared $|\tilde{\theta}|$ for each random EDF to $|\tilde{\theta}|$ for the sample EDF. We reported the fraction $f_s$ of random EDFs that passed the Kolmogorov–Smirnov test and for which $|\tilde{\theta}|$ was greater than in the sample EDF. Values of $f_s$ close to 1 indicate that the considered volcanic eruption locations are closer to the
considered BMSs than random points (in a mean sense) and that the sample distribution is significantly different from the random distribution. By contrast, values of $f_s$ close to 0 indicate that random points are closer to the considered BMSs than the considered volcanic eruption locations (in a mean sense), and/or that the sample distribution is not significantly different from the random distribution.

## Data availability

Data generated for this study are available at https://doi.org/10.5281/zenodo.6031641. Source data are provided with this paper.

## Code availability

The code used to compute the mantle flow models is available at https://github.com/EarthByte/citcoms. Figure 1 was created with the Generic Mapping Tools⁸² (GMT6) which is open-source software licensed under the GNU Lesser General Public License. Key Python scripts to compute the results shown in Figs. 2, 3 are available at https://doi.org/10.5281/zenodo.6031641. Figures 2, 3 were produced with the open-source Python library Matplotlib⁸³, and Fig. 4 was created with ParaView⁶⁷, which is shared openly under the 3-Clause BSD License.

50. Zhong, S., McNamara, A., Tan, E., Moresi, L. & Gurnis, M. A benchmark study on mantle convection in a 3-D spherical shell using CitcomS. *Geochem. Geophys. Geosyst.* **9**, Q10017 (2008).

51. Gurnis, M. et al. Plate tectonic reconstructions with continuously closing plates. *Comput. Geosci.* **38**, 35–42 (2012).

52. Bower, D. J., Gurnis, M. & Flament, N. Assimilating lithosphere and slab history in 4-D Earth models. *Phys. Earth Planet. Inter.* **238**, 8–22 (2015).

53. Stadler, G. et al. The dynamics of plate tectonics and mantle flow: from local to global scales. *Science* **329**, 1033–1038 (2010).

54. Coltice, N., Husson, L., Faccenna, C. & Arnould, M. What drives tectonic plates? *Sci. Adv.* **5**, eaax4295 (2019).

55. Flament, N. et al. Topographic asymmetry of the South Atlantic from global models of mantle flow and lithospheric stretching. *Earth Planet. Sci. Lett.* **387**, 107–119 (2014).

56. Chopelas, A. & Boehler, R. Thermal expansivity in the lower mantle. *Geophys. Res. Lett.* **19**, 1983–1986 (1992).

57. Tosi, N., Yuen, D. A., de Koker, N. & Wentzcovitch, R. M. Mantle dynamics with pressure- and temperature-dependent thermal expansivity and conductivity. *Phys. Earth Planet. Inter.* **217**, 48–58 (2013).

58. Hassan, R., Flament, N., Gurnis, M., Bower, D. J. & Müller, R. D. Provenance of plumes in global convection models. *Geochem. Geophys. Geosyst.* **16**, 1465–1489 (2015).

59. Jaupart, C., Labrosse, S. & Mareschal, J. In Treatise on Geophysics. Volume 7: Mantle Dynamics 1st edn (ed. Bercovici, D.) 253–303 (Elsevier, 2007).

60. Steinberger, B. & Calderwood, A. R. Models of large-scale viscous flow in the Earth's mantle with constraints from mineral physics and surface observations. *Geophys. J. Int.* **167**, 1461–1481 (2006).

61. Billen, M. I. Modeling the dynamics of subducting slabs. *Annu. Rev. Earth Planet. Sci.* **36**, 325–356 (2008).

62. Williams, S., Wright, N. M., Cannon, J., Flament, N. & Müller, R. D. Reconstructing seafloor age distributions in lost ocean basins. *Geosci. Front.* **12**, 769–780 (2021).

63. van der Meer, D. G., Spakman, W., van Hinsbergen, D. J., Amaru, M. L. & Torsvik, T. H. Towards absolute plate motions constrained by lower-mantle slab remnants. *Nat. Geosci.* **3**, 36–40 (2010).

64. Hernlund, J. W. & Houser, C. On the statistical distribution of seismic velocities in Earth's deep mantle. *Earth Planet. Sci. Lett.* **265**, 423–437 (2008).

65. Dziewonski, A. M. & Anderson, D. L. Preliminary reference Earth model. *Phys. Earth Planet. Inter.* **25**, 297–356 (1981).

66. Ricard, Y., Richards, M., Lithgow-Bertelloni, C. & Le Stunff, Y. A geodynamic model of mantle density heterogeneity. *J. Geophys. Res.* **98**, 21895–21909 (1993).

67. Ahrens, J., Geveci, B. & Law, C. Paraview: an end-user tool for large data visualization. In *The Visualization Handbook* (eds Hansen, C. D. & Johnson, C. R.) 717–731 (Academic Press, 2005).

68. Müller, R. D. et al. Ocean basin evolution and global-scale plate reorganization events since Pangea breakup. *Annu. Rev. Earth Planet. Sci.* **44**, 107–138 (2016).

69. Torsvik, T. H. et al. Phanerozoic polar wander, palaeogeography and dynamics. *Earth Sci. Rev.* **114**, 325–368 (2012).

70. Torsvik, T. H., Müller, R. D., van der Voo, R., Steinberger, B. & Gaina, C. Global plate motion frames: toward a unified model. *Rev. Geophys.* **46**, RG3004 (2008).

71. Torsvik, T. H. & Voo, R. V. D. Refining Gondwana and Pangea palaeogeography: estimates of Phanerozoic non-dipole (octupole) fields. *Geophys. J. Int.* **151**, 771–794 (2002).

72. Merdith, A. S. et al. A full-plate global reconstruction of the Neoproterozoic. *Gondwana Res.* **50**, 84–134 (2017).

73. Domeier, M. A plate tectonic scenario for the Iapetus and Rheic oceans. *Gondwana Res.* **36**, 275–295 (2016).

74. Domeier, M. Early Paleozoic tectonics of Asia: towards a full-plate model. *Geosci. Front.* **9**, 789–862 (2018).

75. MacQueen, J. Some methods for classification and analysis of multivariate observations. In *Proc. Fifth Berkeley Symp. Mathematical Statistics and Probability: Volume 1* (eds Le Cam, L. M. & Neyman, J.) 281–297 (Univ. California Press, 1967).

76. Bryan, S. E. & Ernst, R. E. Revised definition of large igneous provinces (LIPs). *Earth Sci. Rev.* **86**, 175–202 (2008).

77. Coffin, M. F. et al. Large igneous provinces and scientific ocean drilling: Status quo and a look ahead. *Oceanography* **19**, 150–160 (2006).

78. Ernst, R. E. Large Igneous Provinces (Cambridge Univ. Press, 2014).

79. Richards, M. A., Duncan, R. A. & Courtillot, V. E. Flood basalts and hot-spot tracks: plume heads and tails. *Science* **246**, 103–107 (1989).

80. Müller, R. D. et al. GPlates: building a virtual Earth through deep time. *Geochem. Geophys. Geosyst.* **19**, 2243–2261 (2018).

81. Knuth, D. E. Art of Computer Programming. Volume 2: Seminumerical Algorithms (Addison-Wesley, 2014).

82. Wessel, P., Smith, W. H., Scharroo, R., Luis, J. & Wobbe, F. Generic mapping tools: improved version released. *Eos* **94**, 409–410 (2013).

83. Hunter, J. D. Matplotlib: a 2D graphics environment. *IEEE Ann. Hist. Comput.* **9**, 90–95 (2007).

**Acknowledgements** N.F. and Ö.F.B. were supported by Australian Research Council grant LP170100863 (industry partner: De Beers). S.E.W. was supported by National Natural Science Foundation of China grant 41972237 and by Australian Research Council grants DP180102280 and DP200100966. A.S.M. was supported by the Deep Carbon Observatory and the Richard Lounsbery Foundation. 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 and through the University of Wollongong (UOW).

**Author contributions** N.F.: conceptualization, methodology, software, validation, formal analysis, investigation, writing of the original draft, and visualization. Ö.F.B.: methodology, software, formal analysis, review and editing of the paper, and visualization. S.E.W.: methodology, software, and review and editing of the paper. A.S.M.: resources and review and editing of the paper.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-04538-y.

**Correspondence and requests for materials** should be addressed to Nicolas Flament.

**Peer review information** Nature thanks Allen McNamara and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 1** | Temporal distribution of considered volcanic eruptions. Histograms showing the number of volcanic eruptions as a function of age. a, Back to 320 Ma. b, Back to 640 Ma. The number of volcanic

<figure><img src="image_8.png" /></figure>

eruptions over the considered period for each database is given in brackets.
Databases: J18 (ref. ⁴⁰) and T18 (ref. ⁴¹).
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 2** | Location of volcanic eruptions from 320 Ma in different reconstructions. High-velocity (white) and low-velocity (grey) regions revealed by k-means cluster analysis between 1,000 km and 2,800 km depth for seismic tomographic model Savani²⁶, and location of volcanic eruptions (diamonds, J18, ref. ⁴⁰) and kimberlites (circles, T18, ref. ⁴¹) reconstructed at their time of eruption and shown at present day using

tectonic reconstruction M16 (ref. ⁴⁴; a), tectonic reconstruction Y19 (ref. ⁴⁵; b), tectonic reconstruction M21 (ref. ¹⁰; c), and tectonic reconstruction M21-NNR without net rotation (no net rotation or 'NNR') (d). In a–d, the brown lines are present-day coastlines. Symbols are coloured by age. Robinson projection at Earth's surface.
<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 3** | Spatial match of lower mantle structure between mantle flow and tomographic models with respect to tomographic model Savani, and distance to hot basal mantle structures from Savani. a,</figcaption></figure>

Intersection of mantle structure revealed by *k*-means cluster analysis between 1,000 km and 2,800 km for tomographic model Savani²⁶ and mantle flow model cases 1–9 as indicated. Orange (true positive) indicates high-temperature model regions and low-velocity tomographic regions, grey (true negative) indicates low-temperature model regions and high-velocity tomographic regions, green (false positive) indicates high-temperature model regions and high-velocity tomographic regions and blue (false negative) indicates low-temperature regions and low-velocity tomographic regions.

b. As in *a* but for tomographic model Savani and other tomographic models as

indicated. 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 Savani and low-velocity regions for other models and blue (false negative) indicates low-velocity regions for Savani and high-velocity regions for other models.

c. Angular distance (AD) to the edge or interior (as indicated) of LLSVPs (delineated by a red contour) as derived by *k*-means cluster analysis of tomographic model Savani between 1,000 km and 2,800 km depth. In *a* and *b*, the white lines indicate a value of five (solid) and a value of one (dotted) in a vote map for low-velocity regions in S-wave tomographic models⁴². Present-day coastlines are shown in black in *a*–*c*. Robinson projection at Earth's surface.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 4** | Distributions of distances between basal mantle structures and volcanic eruption locations. Sample empirical distribution functions (EDFs; blue lines) showing the cumulative probability of minimum angular distances between volcanic eruption locations (J18 and T18) and the

closest BMS for the last 320 Myr. Grey lines are a series of 1,000 random EDFs each consisting of points in random locations with the same temporal distribution as in the sample EDF. **a**, For S-wave tomographic models as indicated. **b**, For mantle flow models cases 1–22 as indicated.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 5** | Match between volcanic eruption locations, tomographic models and mantle flow models for different tectonic reconstructions. a, Fractional area *f*<sub>a</sub> of the surface of cluster maps covered by slow (in tomography) or hot (in flow models, averaged over 320 Myr) BMSs. b, Spatial match *Acc* between present-day BMSs for a given case and LLSVPs imaged by tomographic models. c, Time-averaged median of minimum angular distances θ̃ between BMSs and volcanic eruption locations from 320 Ma. d, Fraction *f*<sub>s</sub> of random EDFs compared to which the sample EDF passes a

statistical test. In a–d, the first four rows show results for tomographic models (stationary LLSVPs) and different tectonic reconstructions as indicated, and the fifth and last row shows results for a series of mantle flow models (mobile BMSs) based on different tectonic reconstructions as indicated. NNR, no-net-rotation reference frame; MFM, mantle flow model. Tectonic reconstructions: M16, ref. ⁴⁴; Y19, ref. ⁴⁵; M21, ref. ¹⁰. The horizontal lines denote reference cases. The grey shadings in a–c highlight the range of results for tomographic models.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 6** | Match of model basal mantle structures to volcanic eruption locations and tomographic models: structure interiors and over 640 Myr. **a–d**, Match to BMS interiors. **a**, Fractional area *f*<sub>a</sub> of the surface of cluster maps covered by slow (in tomography) or hot (in flow models, averaged over 320 Myr) BMSs. **b**, Spatial match Acc between present-day BMSs for a given case and LLSVPs imaged by tomographic models. **c**, Time-averaged median of minimum angular distances θ̂ between BMSs and volcanic eruption locations from 320 Ma, considering distances to be zero within BMSs and positive outwards from their edges. **d**, Fraction *f*<sub>a</sub> of random EDFs compared to which the sample EDF passes a statistical test. **e–h**, Match to BMS edges over 640 Myr. **e**, Fractional area *f*<sub>a</sub> of the surface of cluster maps covered by hot (in flow

<figure><img src="image_15.png" /></figure>

models, averaged over 640 Myr) or slow (in tomography) BMSs. **f**, Same as **b**. **g**, Time-averaged median of minimum angular distances θ̂ between BMSs and volcanic eruption locations from 640 Ma. **h**, Same as **d**. **In g, h**, open symbols denote results for reconstruction M21 as opposed to M21-NNR. In **a–h**, the first row shows results for a series of tomographic models (stationary LLSVPs), and subsequent rows show results for mantle flow models (mobile BMSs) across which parameters are varied as indicated. BL, basal layer; T1–T7, tomographic models 1–7; C1–C12, mantle flow model cases 1–12. The horizontal lines denote reference cases. The grey shadings in **a–c, e–g** highlight the range of results for tomographic models.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 7** | Match of model basal mantle structures to volcanic eruption locations and tomographic models: varying model parameters and African hemisphere. **a–d**, Match for varied model parameters as indicated. **a**, Fractional area $f_a$ of the surface of cluster maps covered by slow (in tomography) or hot (in flow models, averaged over 320 Myr) BMSs. **b**, Spatial match $\bar{Acc}$ between present-day BMSs for a given case and LLSVPs imaged by tomographic models. **c**, Time-averaged median of minimum angular distances $\tilde{\theta}$ between BMSs and volcanic eruption locations from 320 Ma. **d**, Fraction $f_s$ of random EDFs compared to which the sample EDF passes a statistical test. **e–h**, Match in the African hemisphere. **e**, Same as **a** but for different model cases. **f**, Global spatial match $\bar{Acc}$ between present-day BMSs for a given case and LLSVPs imaged by tomographic models. **g**, Time-averaged median of

<figure><img src="image_17.png" /></figure>

minimum angular distances $\tilde{\theta}$ between BMSs and volcanic eruption locations in the African hemisphere (within 10,000 km of a point located at 0° N, 11° E; ref. 47) from 320 Ma. **h**, For points within the African hemisphere, fraction $f_s$ of random EDFs compared to which the sample EDF passes a statistical test. In **a–h**, the first row shows results for a series of tomographic models (stationary LLSVPs), and subsequent rows show results for mantle flow models (mobile BMSs) across which parameters are varied as indicated. C1–C22, mantle flow model cases 1–22; T1–T7, tomographic models 1–7; ISD, initial slab depth; PCD, phase change depth; BL, basal layer. The horizontal lines denote reference cases. The grey shadings in **a–c**, **e–g** highlight the range of results for tomographic models.
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 8** | Snapshots of tectonic reconstruction and synthetic seafloor ages. a–f, Location of tectonic blocks for tectonic reconstruction M21-NNR¹⁰ in the no-net-rotation frame with synthetic seafloor ages⁶². Snapshots are shown at: 740 Ma (a), 550 Ma (b), 400 Ma (c), 250 Ma (d), 100 Ma (e), 0 Ma (f). The reconstructed subduction zones are shown as black lines with triangles on the overriding plate, reconstructed mid-oceanic ridges and

transform faults are shown in light orange, reconstructed ancient continental interiors are shown as light grey polygons with a pattern fill, and present-day coastlines are shown in white in f. The colour scale indicates the palaeo-age of the ocean crust. WGB, west Gondwanan blocks; PTNM, Palaeo-Tethys northern margin; GWM, Gondwana western margin. Robinson projection at Earth's surface.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

**Extended Data Fig. 9** | Model initial temperature field and present-day temperature and viscosity. a, b, Initial temperature field for case 7 at 1,000 Ma. a, Temperature at 109 km depth. Reconstructed subduction zones are shown in red, mid-oceanic ridges and transform faults in yellow, and continental blocks in grey. Robinson projection at Earth's surface.

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

b, Temperature along an equatorial cross-section (green line in a). Numbers above the colour scale indicate non-dimensional temperature, and numbers below the colour scale indicate temperature in kelvin. c, d, Horizontally averaged present-day mantle temperature (c) and viscosity (d). The grey line in d is a viscosity profile adjusted to fit the geoid and post-glacial rebound⁶⁰.
Extended Data Table 1 | Boundary conditions and parameters for mantle flow model cases

<table><thead><tr><th>Model case</th><th>Tectonic reconstruction</th><th>a<sub>0</sub> (Myr ago)</th><th>z<sub>0</sub> (km)</th><th>η<sub>b</sub></th><th>δρ<sub>b</sub> (%)</th><th>z<sub>r</sub> (km)</th></tr></thead><tbody><tr><td>Case 1</td><td>M21-NNR</td><td>1,000</td><td>550</td><td>1</td><td>1.0</td><td>-</td></tr><tr><td>Case 2</td><td>M21-NNR</td><td>1,000</td><td>550</td><td>10</td><td>1.0</td><td>-</td></tr><tr><td>Case 3</td><td>M21-NNR</td><td>1,000</td><td>550</td><td>0.1</td><td>1.0</td><td>-</td></tr><tr><td>Case 4</td><td>M21-NNR</td><td>1,000</td><td>1,000</td><td>10</td><td>0</td><td>-</td></tr><tr><td>Case 5</td><td>M21-NNR</td><td>1,000</td><td>1,000</td><td>10</td><td>0.5</td><td>-</td></tr><tr><td>Case 6</td><td>M21-NNR</td><td>1,000</td><td>1,000</td><td>10</td><td>1.0</td><td>-</td></tr><tr><td>Case 7</td><td>M21-NNR</td><td>1,000</td><td>1,000</td><td>10</td><td>1.3</td><td>-</td></tr><tr><td>Case 8</td><td>M21-NNR</td><td>1,000</td><td>1,000</td><td>10</td><td>1.6</td><td>-</td></tr><tr><td>Case 9</td><td>M21-NNR</td><td>1,000</td><td>1,000</td><td>10</td><td>2.0</td><td>-</td></tr><tr><td>Case 10</td><td>M16</td><td>410</td><td>1,200</td><td>10</td><td>1.0</td><td>-</td></tr><tr><td>Case 11</td><td>Y19</td><td>410</td><td>1,200</td><td>10</td><td>1.0</td><td>-</td></tr><tr><td>Case 12</td><td>M21</td><td>1,000</td><td>550</td><td>10</td><td>1.0</td><td>-</td></tr><tr><td>Case 13</td><td>M21-NNR</td><td>400</td><td>1,200</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 14</td><td>M21-NNR</td><td>500</td><td>550</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 15</td><td>M21-NNR</td><td>600</td><td>550</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 16</td><td>M21-NNR</td><td>700</td><td>550</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 17</td><td>M21-NNR</td><td>800</td><td>550</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 18</td><td>M21-NNR</td><td>900</td><td>550</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 19</td><td>M21-NNR</td><td>1,000</td><td>550</td><td>1</td><td>2.0</td><td>-</td></tr><tr><td>Case 20</td><td>M21-NNR</td><td>1,000</td><td>1,200</td><td>10</td><td>1.0</td><td>-</td></tr><tr><td>Case 21</td><td>M21-NNR</td><td>1,000</td><td>1,200</td><td>10</td><td>1.0</td><td>670</td></tr><tr><td>Case 22</td><td>M21-NNR</td><td>1,000</td><td>1,200</td><td>10</td><td>1.0</td><td>410 &amp; 670</td></tr></tbody></table>

Tectonic reconstructions: M21 is Merdith et al. (2021)<sup>10</sup>; M16 is Matthews et al. (2016)<sup>43</sup> and Y19 is Young et al. (2019)<sup>44</sup>. NNR indicates a no-net-rotation reference frame.

a₀ is the model start age, z₀ is the initial slab depth, η₀ is the composition viscosity pre-factor for BMSs, δρ₀ is the density anomaly for BMSs, and zᵣ is the depth of phase changes.