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

HAL
open science

# A naive Bayesian method to chase mantle plumes in global tomography models

M van Camp, O de Viron, A Ferreira, O Verhoeven

► To cite this version:

M van Camp, O de Viron, A Ferreira, O Verhoeven. A naive Bayesian method to chase mantle plumes in global tomography models. Geophysical Journal International, 2022, 10.1093/gji/ggac415. hal-03831519

HAL Id: hal-03831519

https://hal.science/hal-03831519v1

Submitted on 4 Apr 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# A naive Bayesian method to chase mantle plumes in global tomography models

M. Van Camp<sup>1</sup>, O. de Viron<sup>1,2</sup>, A.M.G. Ferreira<sup>3</sup> and O. Verhoeven<sup>4</sup>

Accepted 2022 October 17. Received 2022 October 3; in original form 2022 February 23

## SUMMARY

This study provides a quantitative approach to search for mantle plumes in global seismic tomography models without any preconceived notions about the associated mantle velocity anomalies, other than the assumption that the plumes are not significantly deflected horizontally by more than 6°, anywhere in the mantle. We design identification tests with a reasonable detection threshold while keeping false alarms at a level lower than 5 per cent. This is based on a naive Bayesian clustering analysis, which is possible thanks to the varimax principal component analysis that provides components of the tomography models that are much more independent than the original number of depth slices in the models. We find that using such independent components greatly reduces detection errors compared to using an arbitrary number of depth slices due to correlations between the different slices. We detect a wide range of behaviour of the seismic velocity profiles underneath the hotspots investigated in this study. Moreover, we retrieve locations away from hotspots that have similar seismic velocity profile signatures to those underneath some hotspots. Hence, it is not possible to obtain a unique definition of seismic velocity anomalies that are associated with hotspots and thus care needs to be taken when searching for mantle plumes beneath hotspots using prior assumptions about the velocity anomalies that might be associated with them. On the other hand, we establish a probability distribution of the seismic velocity profiles that is specific to a sublist of hotspots. Overall, the mantle plume zones identified in our analysis do not appear to surround the Africa and Pacific large low shear velocity provinces (LLSVPs), but are rather within them. This rules out the idea that LLSVPs correspond to compact, dense piles with mantle plumes rising from their edges. Instead, our analysis suggests two possible options that either the LLSVPs: (1) correspond to bundles of thermochemical mantle plumes or (2) are made up of compact piles topped by a bundle of plumes.

**Key words:** Composition and structure of the mantle; Statistical methods; Seismic tomography; Hotspots.

## 1 INTRODUCTION

The hotspot paradigm has been extensively debated over the past 60 yr, since it was proposed that volcanic chains like the Hawaiian Islands are the surface expression of hot material—mantle plumes—that originate at the core–mantle boundary (CMB; Wilson 1965; Morgan 1971; Courtillot et al. 2003; Anderson 2005; French & Romanowicz 2014; Duncombe 2019). Alternatively, it has been suggested that hotspots result from melting due to shallow convective processes, with their shapes being controlled by stresses and cracks in the lithosphere (e.g. Foulger & Natland 2003). The coexistence of these processes with mantle plumes is actually possible and is the current subject of investigation.

Mantle plumes have been suggested to start as an instability generated by heat flux across the CMB and then vertically extending to the surface (Morgan 1971). More recently, several studies showed that plumes can be thermochemical and have a wide range of shapes and internal structures (Lin & van Keken 2006; Kumagai et al. 2008; Flament et al. 2017). Moreover, many hotspots lie on top or near the so-called large low shear velocity provinces (LLSVPs), which are continent-sized structures sitting on top of the CMB beneath Africa and the Pacific (Garnero & McNamara 2008; Austermann et al. 2014; French & Romanowicz 2015) whose nature and origin are still debated.

While some studies have suggested that LLSVPs are warmer and lighter than the average, and directly related to convective

<footnote n="1"><sup>1</sup>Royal Observatory of Belgium, Avenue Circulaire, 3, 1180 Uccle, Belgium. E-mail: m.vancamp@oma.be</footnote>

<footnote n="2"><sup>2</sup>Littoral, Environnement et Sociétés (LIENSs UMR7266 La Rochelle University - CNRS), 17000 La Rochelle, Charente Maritime, France</footnote>

<footnote n="3"><sup>3</sup>Department of Earth Sciences, University College London, London WC1E 6BT, United Kingdom</footnote>

<footnote n="4"><sup>4</sup>Nantes Université, Université d'Angers, Le Mans Université, Laboratoire de Planétologie et Géosciences, CNRS UMR 6112, France</footnote>
upwellings (e.g. Austermann et al. 2014; Koelemeijer et al. 2017), other studies indicated that they may be dense, long-lived chemically heterogeneous regions (e.g. Trampert 2004; Lau et al. 2017). Alternatively, a recent study has proposed that LLSVPs consist of bundles of thermochemical, dense plumes rooted in large ultra low-velocity zones (ULVZs) in the lowermost mantle (Davaille & Romanowicz 2020). Seismic imaging is arguably the most powerful way to discriminate mantle plumes in the Earth's interior and sheds light on their relations to the lowermost mantle. By considering tomography models together with other information such as, for example, the presence of long-lived volcanic chains, large igneous provinces, high buoyancy fluxes, ³He/⁴He ratios and variations in the thickness of the transition zone, one can investigate the existence and structure of deep mantle plumes (Montelli 2004; Zhao 2004; Suetsugu et al. 2009; French & Romanowicz 2015; Chang et al. 2020). Although tomography has been widely used to chase mantle plumes, other recent efforts are being developed, such as modelling of S waveforms (Cottaar & Romanowicz 2012) or the use of array seismology stacking techniques (Stockmann et al. 2019).

In the past decades, the global seismic imaging of the Earth's interior—global seismic tomography—has greatly progressed, so that recent tomographic models have started to robustly image mantle plumes (French & Romanowicz 2015; Chang et al. 2015; Bozdağ et al. 2016). For example, French & Romanowicz (2015) have identified the presence of fat mantle plumes beneath many hotspots, which may have a thermochemical origin. They classified hotspots as 'primary', 'learly resolve' and 'somewhat resolved' depending on the strength and on the vertical continuity of the negative shear wave speed anomalies associated with them in the mantle. Some regions exhibit substantial complexity, such as, for example a continuous, strong low shear wave speed anomaly in the mantle beneath Africa, from Afar to South Africa, with many parts of the anomaly not being associated with hotspots. Some studies (e.g. Davaille & Romanowicz 2020; Chang et al. 2020) suggested that this apparently continuous anomaly may actually result from the superimposition of several instabilities associated with different plumes at different stages of their evolution. For example, it is possible that beneath South Africa there is a plume that has not reached the upper mantle yet (and so there is no hotspot yet), while the Kenya and Afar hotspots are possibly underlain by currently active and dying plumes, respectively (e.g. Chang et al. 2020).

In order to identify mantle plumes in tomography models, previous studies typically searched for continuous low seismic wave velocity anomalies spanning both the Earth's upper and lower mantle (e.g. Boschi et al. 2007, 2008; Marignier et al. 2020). Moreover, Boschi et al. (2007) also included in their analysis comparisons with geodynamic models, which may add further constraints, but that are also substantially affected by uncertainties. Marignier et al. (2020) recently proposed a method to assess the probability of mantle plumes in global tomography models based on the assumption that a plume shows a negative shear wave speed anomaly throughout the whole mantle not tilting more than 10°. A limitation of this approach is that this plume definition does not account for well-known complexities in plume behaviour, such as plume heads, deflections and gaps. Hence, general and objective criteria for mantle plume detection in tomography models are still lacking.

In this study, we propose to chase mantle plumes in global tomography models following a different approach from previous studies, whereby we do not rely on prior assumptions about the low-velocity anomalies that may be associated with mantle plumes. We use a novel, alternative approach based on naive Bayesian clustering methods whereby we investigate to what extent the mantle

structure beneath hotspots is substantially distinct from other different regions. As explained below, these methods allow the natural separation of distinct velocity profiles in tomography models beneath hotspots, oceans and continents. Other statistical methods that have been used in seismology to classify structures in tomography models include, for example the k-means clustering method (e.g. Lekic et al. 2012; Cottaar & Lekic 2016). Unlike the k-means clustering, which classifies the data points based only on their values, the naive Bayesian clustering used in this study is based on a pre-classification used to define the clustering criteria.

Using the hotspot list given by Marignier et al. (2020), which is based on French & Romanowicz (2015), Sleep (1990) and Courtillot et al. (2003), we show that our Bayesian method allows a fairly robust test on the presence of mantle plumes in global tomography models. It recovers shear wave velocity anomaly patterns associated with hotspots while limiting false detections (type I error). We also show that false alarm locations, that is geographical locations with no hotspots but identified as corresponding to mantle velocity anomaly profiles associated with hotspots, make geophysical sense. Finally, we show that a sublist of seismic velocity anomaly profiles associated with hotspots does not behave similarly as the majority of hotspots. Based on these results, we conclude that (1) the hotspots considered in this study are not the surface expression of a unique type of mantle shear wave velocity anomaly, which would require further classification (i.e. mantle plumes beneath hotspots are quite diverse) and (2) it is possible to establish a probability distribution of the seismic velocity profiles that is specific to a sublist of hotspots and which does not occur significantly in locations away from hotspots.

# 2 DATA AND METHOD

## 2.1 Data

In this study we search for mantle plumes in eight 3-D global seismic tomography models: (i) S20RTS (Ritsema 1999); (ii) S40RTS (Ritsema et al. 2011); (iii) SEISGLOB2 (Durand et al. 2017); (iv) SEMUCB-WM1 (French & Romanowicz 2014); (v) SGLOBE-rani (Chang et al. 2015); (vi) S362WMANI+M (Moulik & Ekström 2014); (vii) SAVANI (Auer et al. 2014) and (viii) GLAD-M25 (Lei et al. 2020). These models were built with different data sets and modelling approaches, as summarized in Table 1. The models used show the main features in current global isotropic and radially anisotropic tomography models. Thus, they are representative of the current state of global seismic tomography. Given the disagreement in the anisotropic structures in global anisotropy models, in this study we focus only on the isotropic structure and always refer to it in the remainder of this paper.

We retrieved the global tomography models either from the IRIS earth model collaboration repository (http://ds.iris.edu/ds/product/emc-earthmodels/; Trabant et al. 2012) or directly upon request from their authors. Since some of the models have distinct reference 1-D models and use different parametrizations, for consistency they have been all converted into variations in shear wave speed and in radial anisotropy with respect to the 1-D reference model PREM (Dziewonski & Anderson 1981) on a common 1° × 1° horizontal grid and for 29 depth slices starting at 50 km and with a 100 km spacing from 100 to 2900 km depth. This conversion into the same 1-D reference model simplifies the comparison of the various tomography models. Moreover, since both the vertical structure and the horizontal patterns in the tomography models are normalized,
<table><caption>**Table 1.** Global tomography models used in this study, including a short description of the data set, parametrization and the modelling approach used in their construction. The models were built using least-squares inversions with different regularization choices or non-linear gradient-based optimization techniques including smoothing of the misfit function's gradient (Table 1). Except for the very recent GLAD-M25 model, the same models were used by de Viron *et al.* (2021) for varimax model analyses and comparisons.</caption><thead><tr><td>Model</td><td>Data set</td><td>Parametrization</td><td>Modelling approach</td></tr></thead><tbody><tr><td>SGLOBE-rani (Chang et al. 2015)</td><td>Fundamental mode group-velocity data (T ~ 16–150 s), fundamental and overtone phase-velocity data (25–374 s, up to the 4th overtone), body-wave traveltimes.</td><td>Spherical harmonics up to degree 35 (laterally) and 21 spline functions (radially). 1-D reference model: PREM.</td><td>Ray theory. Regularization: norm damping.</td></tr><tr><td>SAVANI (Auer et al. 2014)</td><td>Fundamental and overtone phase velocities (T ~ 25–370 s, up to the 4th overtone), body-wave traveltimes.</td><td>Variable size blocks on a 5° × 5° base grid (laterally). 28 variable thickness depth layers. 1-D reference model: PREM.</td><td>Ray theory. Regularization: vertical and horizontal smoothing.</td></tr><tr><td>S20RTS (Ritsema 1999)</td><td>Fundamental mode and overtone phase velocities (T ~ 40–275 s, up to the 4th overtone), body-wave traveltimes, even-degree self-coupling normal mode splitting functions.</td><td>Spherical harmonics up to degree 20 (laterally) and 21 spline functions (radially). 1-D reference model: PREM.</td><td>Ray theory. Regularization: norm damping.</td></tr><tr><td>S40RTS (Ritsema et al. 2011)</td><td>Expanded data set of S20RTS (increased amount of measurements).</td><td>Spherical harmonics up to degree 40 (laterally) and 21 spline functions (radially). 1-D reference model: PREM.</td><td>Ray theory. Regularization: norm damping.</td></tr><tr><td>SEMUCB-WM1 (French &amp; Romanowicz 2014)</td><td>Body waveforms (T > 32 s), surface waveforms (T > 60 s).</td><td>Spherical splines with spacing <2° (laterally) and 20 cubic b-splines (radially). 1-D reference model: own model.</td><td>Spectral element method for forward modelling, non-linear asymptotic coupling theory for inverse modelling. Regularization: vertical and horizontal smoothing.</td></tr><tr><td>S362WMANI+M (Moulik &amp; Ekström 2014)</td><td>Fundamental mode phase velocities (T ~ 35–150 s), body wave traveltimes and waveforms, normal mode splitting functions.</td><td>362 spherical splines (laterally) and 16 cubic splines (radially). 1-D reference model: STW105 (Kustowski et al. 2008).</td><td>Ray theory. Regularization: vertical and horizontal smoothing.</td></tr><tr><td>SEISGLOB2 (Durand et al. 2017)</td><td>Fundamental and overtone phase velocities (T ~ 40–360 s, up to the 5th overtone), body-wave traveltimes, normal mode self- and cross-coupling coefficients.</td><td>Spherical harmonics up to degree 40 (laterally) and 21 spline functions (radially)</td><td>Ray theory. Regularization: lateral smoothing controlled by a horizontal correlation length.</td></tr><tr><td>GLAD-M25 (Lei et al. 2020)</td><td>Body and surface waveforms in three different period ranges (17–40 s, 40–100 s and 90–250 s)</td><td>Cubic sphere parametrization (same as the mesh of the SPECFEM3D_GLOBE package)</td><td>Spectral element method for forward modelling. Adjoint method for inverse modelling using a quasi-Newton technique including smoothing of the misfit function's gradient. Source and receiver geographical weighting to compensate for uneven spatial sampling.</td></tr></tbody></table>

this also implies that the reference model used does not influence the analysis. In addition to the gridded representation of the models, we also interpolate the model parameters from the horizontal grids into a regular polyhedron of 9002 equiareal faces of which the vertices were generated through the icosahedron tool of Zechmann (2019). This transformation leads to a uniform sampling of the sphere, with each vertex having a surface corresponding to that of a 250 × 250 km² square. By using regularly gridded data in our analysis we would overweight the contribution of velocity structure at high latitudes. Moreover, as the shallowest depth slices of the tomography models contain the majority of the variability in the velocity anomalies, most of the analysis performed in this study would get focused on those slices, which is not suitable for our purpose. Hence, we normalize the shear wave speed perturbations by depth slice, that is the mean value of the slice is subtracted from

each value, and each value is divided by the standard deviation of the slice (de Viron et al. 2021).

The wave speed data are then compressed using a varimax principal component (PC) algorithm, as described in de Viron et al. (2021). The varimax PC analysis provides components that are much more independent than the original 29 slices of the model, which is better suited to the naive Bayesian clustering method that we use in this study, as explained in the next section. We computed two sets of varimax PCs, corresponding to 95 and 80 per cent of the total model variance, respectively. In this paper, we refer to experiment E-100 when we use the full 29 slices, and we refer to experiments E-95 and E-80 when we perform the two sets of varimax PC analyses. Six (SAVANI) to thirteen (SGLOBE-rani) components are necessary to reach the variance level of 95 per cent, while this reduces to between three (SAVANI) and eight (SGLOBE-rani) for
the 80 per cent variance level (Table 2). As discussed by de Viron et al. (2021), the number of PCs required to reach a given variance level, given in Table 2, depends on the data types used to build the models and on model regularization.

We use the list of 50 hotspots compiled by Marignier et al. (2020) to perform our analysis. Then, we create four sets of geographical locations (Fig. S1):

(A) The seven gridpoints closest to each hotspot location, which represents a total of slightly less than $7 \times 50 = 350$ points on Earth's surface—there is some overlap associated with close hotspots; the choice of seven points is due to the hexagonal geometry of the icosahedral grid. The width of a whole A region is 6° as each seven-point cell is about 6° wide, which corresponds to the nominal resolution of the most detailed global tomography models.

(B) The gridpoints in a 10° range from each hotspot. The construction of this separate location set guarantees that the profiles beneath the C and D points introduced below are not affected by those beneath the locations defined in A.

(C) The ocean gridpoints that are not included in A nor B.

(D) The continent gridpoints that are not included in A nor B.

C and D complement A+B for oceans and continents, respectively. Our classification is based on the strong assumption that the plumes are not significantly deflected horizontally by more than 6°.

## 2.2 Naive Bayesian clustering

The naive Bayesian algorithm (see for example Koller & Friedman 2009) is a fast classification method. It is based on the Bayes theorem, which expresses how the probability of event K considering the observation set V can be estimated from the probability of the observation set V under the condition of event K. The naive part of the classification comes from the hypothesis that the observations from the set V are supposed to be independent from each other. In the case of our mantle plume detection application, it translates into the fact that the probability $P(K|V)$ for a point to belong to a given category K—for example, hotspot—considering its velocity profile can be computed from the probability $P(V|K)$ to have the observed velocity profile if it belongs to the category K, assuming that the probability laws from the different slices are independent of each other. This hypothesis of independence is incorrect for the raw data, that is for the 29-depth slice version of the models, on the other hand, this is a fair approximation for the variance compression of the varimax representation and the associated decoupling of the geographic patterns.

As the velocity profile probability density function for each category is unknown, we first fit the observed distribution of each category and model on 50 continuous probability distributions (Table S1). The distribution parameters were fit using the Maximum Likelihood Estimation (MLE) criterion (Kay 1993; Virtanen et al. 2020). From our tests, we found that the best-fitting probability law, for any of the locations and depths considered, was always one of the three following distributions: $S_U$ and $S_B$ Johnson distributions, and the continuous generalized extreme value distribution. This provides us with the best-fitting distribution shown in Fig. 1 with the hotspots in red, the oceans in blue and the continents in green. The points corresponding to region B (gridpoints in a 10° range from each hotspot) are not considered in this study in order to sample only the regions A (seven points closest to each hotspot location), C (oceans) and D (continents), which are expected to present clearer signatures than those in region B. This means that we would rule

out plumes that would be deflected by more than 6° underneath the hotspots. Then, we estimate the probability of obtaining the velocity profile observed at a given location if the location belongs to the continent, ocean or hotspot category, by a multiplication of the probability for each slice—this is the naive part of the algorithm. Using the Bayes theorem, we can estimate the probability that a given location belongs to each category given the velocity profiles, and the location is then attributed to the most likely category. Explicitly, the probability $P(K_j|V_i)$ that a location i belongs to each of the three categories $K_j$ given the velocity profiles $V_i$ is computed as

$$P(K_j|V_i) = \frac{P(V_i|K_j)P(K_j)}{P(V_i)}$$

with $P(V_i|K_j)$ being the probability of obtaining the velocity profile $V_i$ observed at a given location i if the location belongs to the continent, ocean or hotspot category $K_j$ and $P(K_j)$ being the uniform probability that any point belongs to category $K_j$. The denominator is computed as

$$P(V_i) = \sum_{j=1}^{3} P(V_i|K_j)P(K_j).$$

## 3 RESULTS

Fig. 2 shows the locations that our approach clustered as hotspots, for the three different experiments considered (E-100, E-95 and E-80). These locations are more likely to belong to the hotspot category than to any of the other two categories of geographical locations used in this study (C and D), according to the maximum likelihood criterion. They are mostly located in the southern Pacific and Indian oceans, around Africa up to Iceland in the northern Atlantic and within the East African rifting zone, roughly corresponding to the LLSVPs [Large Low Shear Velocity Provinces, see e.g. Lekic et al. (2012) or French & Romanowicz (2015)]. The areas classified as hotspots, which extend around hotspots, are consistent with the hypothesis that plumes spread laterally at the base of the lithosphere and volcanoes appear where the lithosphere is sufficiently young and thin (Marzocchi & Mulargia 1993; O'Connor et al. 2012; Wang & Li 2021).

Recovering between 1/5 and 1/4 of the Earth's surface, the E-100 experiment suffers from a too high rate of false alarms. Fig. 3 summarizes, for each tomography model and experiment, the false alarm rate as a function of the probability threshold for a point to be a hotspot. For example, if we set the threshold to 0.50 in the E-95 case, for all models but S20RTS, the false alarm rate is below 8 per cent, and this decreases under 5 per cent for all models if we set the threshold to 0.67. With the 29 slices E-100 experiment, a threshold of 0.99 is not enough to limit the false alarms to 5 per cent or equivalently, the significance level to 95 per cent, for 6 of the 8 models.

The white rectangles inside the coloured rectangles show the false alarm rate for continents, and it is found to be always lower than for the oceans. This is expected, given that the continents show depth distributions of velocity anomalies that are much more different than those from the hotspots and the oceans, as illustrated by Fig. 1.

We decide to consider a zone as a hotspot if the probability to belong to that category is equal to or larger than 50 per cent while ensuring a false alarm rate smaller than 5 per cent. We will use the notation P(50,95) in the rest of this article to describe this criterion. As shown in Table 2, applying this P(50,95) choice, the E-100 and E-95 experiments recover comparable surfaces for all models but SAVANI, while the E-80 case recovers from
<table><caption>**Table 2.** Percentage of Earth's surface ('Surface per cent') identified as hotspots with the P(50,95) criterion, for the global tomography models used in this study, considering the whole 29 slices (E-100) or the varimax analysis (E-95 and E-80). The number of hotspots found in the recovered zones ('\#HS') is also given, as well as the number of principal components obtained from the varimax analysis ('\#PCs').</caption><thead><tr><th>Model</th><th>Surface per cent (E-100)</th><th>#HS (E-100)</th><th>#PCs (E-95)</th><th>Surface per cent (E-95)</th><th>#HS (E-95)</th><th>#PCs (E-80)</th><th>Surface per cent (E-80)</th><th>#HS (E-80)</th></tr></thead><tbody><tr><td>SGLOBE-rani</td><td>8.5</td><td>20</td><td>13</td><td>9.0</td><td>26</td><td>8</td><td>6.2</td><td>20</td></tr><tr><td>SAVANI</td><td>10.0</td><td>22</td><td>6</td><td>3.3</td><td>12</td><td>3</td><td>1.4</td><td>7</td></tr><tr><td>S20RTS</td><td>9.1</td><td>23</td><td>11</td><td>10.0</td><td>23</td><td>5</td><td>4.3</td><td>11</td></tr><tr><td>S40RTS</td><td>8.9</td><td>20</td><td>11</td><td>8.9</td><td>21</td><td>6</td><td>3.9</td><td>16</td></tr><tr><td>SEMUCB-WM1</td><td>10.0</td><td>23</td><td>8</td><td>8.4</td><td>25</td><td>5</td><td>6.6</td><td>21</td></tr><tr><td>S362WMANI+M</td><td>10.9</td><td>27</td><td>7</td><td>8.1</td><td>20</td><td>4</td><td>1.0</td><td>6</td></tr><tr><td>SEISGlob2</td><td>9.5</td><td>23</td><td>11</td><td>6.9</td><td>23</td><td>7</td><td>2.5</td><td>12</td></tr><tr><td>GLAD-M25</td><td>10.6</td><td>25</td><td>8</td><td>7.3</td><td>22</td><td>5</td><td>3.4</td><td>19</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distribution of the velocity anomalies (x-axis) as a function of depth (y-axis) for eight different tomography models and the three different sets of geographical locations considered in this study (see main text for details): the 50 hotspots (locations A, red), the oceans (locations C, blue) and the continents (locations D, green). This is done for the 29 depth slices (E-100, top row) directly obtained from the models, as well as for the components recovered by the varimax analysis, keeping the number of components necessary to recover 95 per cent (E-95, middle row) and 80 per cent (E-80, bottom row) of the total variance of the signal. We better recover the signal's details in the upper mantle and transition zone in the E-95 case than in the E-80 one. For each distribution, a probability law is adjusted. The vertical lines connect the maximums values of the different distributions, for the three categories of geographical locations considered.</figcaption></figure>

30 per cent (S362WMANI+M) to 84 per cent ((SEMUCB-WM1) of the hotspots found in the E-95 case. The E-80 experiment is too selective. Hence, it is a fair compromise to select the E-95 experiment considering its reasonable false alarm rate and the criterion of independence of the slices required as a hypothesis of the Bayesian analysis.

There are also differences between the areas identified as hotspots by the different models, whatever the experiment (E-100, E-95 and E-80). In the E-95 case, the model S20RTS, with its high false alarm rate (Fig. 3), classifies 10.0 per cent of the Earth's surface as hotspots, which is 309 per cent larger than for the SAVANI model, which classifies as hotspots only 3.3 per cent of the Earth's surface.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Locations identified as most probable for belonging to the hotspot category (set A of geographical locations; see main text for details), according to the maximum likelihood criterion (MLE), for the eight models used in this study, for the E-100 (29 slices, top), E-95 (95 per cent of the total variance, middle, also shown for the P(50,95) criterion) and E-80 (80 per cent of the total variance, bottom) cases. The 50 hotspots from the list of Marignier *et al.* (2020) are plotted as red dots.</figcaption></figure>

For the six other models, with comparable surfaces, between 7.3 and 9.0 per cent, the recovered zones can differ to a notable extent; see for example the differences in the south Atlantic and the central Pacific in Fig. 2.

For the experiment E-95 and the P(50,95) criterion, the Bayesian test recovers at least one of the 7 points closest to the hotspot locations, ranging from 12 (for SAVANI) to 26 (for SGLOBE-rani) hotspots out of the 50 included in the list (Fig. 4 and Table S2). At least seven models recover 9 of the 12 ‘highly likely’ hotspots of

Marignier *et al.* (2020), while the three remaining ones are found by six models (MacDonald and Marquesas) and three models (Marion). On the other hand, at least seven models recover 8 of the 11 hotspots classified as category 1 (‘primary plume’) by French & Romanowicz (2015) and one (Caroline) as category 2 (‘clearly resolved’; Fig. 4). Out of the three other hotspots belonging to category 1 of French & Romanowicz (2015), Mac Donald and Marquesas are recovered by six models and the Comoros one is not evidenced by our E-95 analysis. This is consistent with the data used, as the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Percentage of false alarms as a function of the probability threshold considered as acceptable to consider a point on Earth as a hotspot, for the cases E-100 (top), E-95 (middle) and E-80 (bottom), for the eight different global tomography models considered in this study. For example, for E-95, if we consider that a location is a hotspot if the depth distribution of the velocity anomalies has a probability of being larger than 0.5 (thresh. = 0.5) to obey the hotspot-derived velocity distribution, then all models experience less than 7 per cent of false alarms, except S20RTS, which is above 8 per cent. 'Max. likelihood' is the threshold above which the hotspots are the most probable. The white rectangles show the results of the same analysis but for the continents (set D of geographical locations; see the main text for details). The dotted horizontal line shows the 5 per cent false alarm level.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Hotspots recovered by our analysis using the eight global tomography models considered in this study, applying the P(50,95) criterion, for the principal component analysis experiment E-95. This corresponds to at least one reddish slice of the middle pie charts in Figs 6 and S2–S8. A detection (one coloured bin, with different colours being used for different tomography models; see key at the bottom of the figure) means recovery in at least one of the seven locations that are the closest to the hotspot. The hotspots that are underlined correspond to the category identified as 'Primary Plume' by French & Romanowicz (2015) while the hotspots in bold correspond to the 'Highly likely' mantle plume category identified by Marignier *et al.* (2020).</figcaption></figure>

Comoros profile is similar to the oceanic one for S20RTS, S40RTS, SEMUCB-WM1, GLAD-M25 and SAVANI (Fig. 5), or oscillates around the oceanic average profile for SGLOBE-rani, SEISGLOB2 and S3622WMANI+M. The other retrieved hotspots are model-dependent and it is difficult to assess whether the 12 hotspots (Fig. 4) recovered by only one or two models result from false alarms or not.

Fig. 6 (for model SGLOBE-rani) and Figs S2–S8 (for other models) display the significance of the probability of a location to correspond to a hotspot for each of the seven points closest to the hotspot, considering the 29 slices (E-100, outer ring) and the varimax principal component analysis E-95 (middle ring) and E-80 (inner ring). The probabilities above 0.5 are enhanced in reddish colours on the figure, whereas the pie charts are dashed when the significance level
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Velocity profiles for the SAVANI model for the ocean off the hotspots points in a radius of 10° (locations B). All ocean profiles are in light blue, the average is in dark blue and the average over the seven points around the hotspot (locations A) is in red.</figcaption></figure>

is below 95 per cent or, equivalently, when the false alarm rate is above 5 per cent. This false alarm rate is estimated from the number of ocean locations from category C classified as hotspots over the total number of category C locations.

Unfortunately, keeping the false alarm rate smaller than 5 per cent is always at the expense of event detection. For example, considering the E-95 case for the mantle model SGLOBE-rani, 48 per cent of the hotspots are not detected as we only find a hotspot signature for 26 locations: Afar, Amsterdam, Ascension, Azores, Bowie, Canary, Cape Verde, Caroline, Crozet, Darfur, Discovery, East Africa, Easter, Galapagos, Great Meteor, Hawaii, Iceland, Jan Mayen, Juan de Fuca, Louisville, Mac Donald, Pitcairn, Reunion, Samoa, Tahiti and Yellowstone. This high level of missed detections (type II error), which rises to 76 per cent for SAVANI, means that based on our list, there is no common definition of the shear wave velocity profiles underneath hotspots. It also highlights that some hotspots are probably not the surface expression of quasi-radial mantle plumes or correspond to plumes difficult to image in tomography models (e.g. very thin plumes beyond the models' limits of resolution), as other studies have also pointed out (French & Romanowicz 2015; Marignier et al. 2020). This also implies that discrimination of the hotspots based on the velocity profiles alone would be specific to a given type of hotspot.

To further demonstrate the specificity of the hotspot velocity profiles, that is to show that the hotspot list corresponds to locations with specific velocity profiles shared by a large proportion of them, we made 250 Monte Carlo simulations. For each Monte Carlo run, we randomly picked 50 points at the Earth's surface that we use as surrogate hotspots. Then, we applied exactly the same process as

above: take the seven closest points from each as hotspot points, remove the points in the 10° range, separate what is left between ocean and continent, fit the probability distributions, and apply the naive Bayesian algorithm to recover the surrogate hotspot locations. From this test, we found that for 67–78 per cent of the cases, depending on the model, no hotspots are recovered, while one to three hotspots are recovered in 20–26 per cent (1 false alarm hotspot), 2–4 per cent (2 false alarm hotspots) and 0–2 per cent (3 false alarm hotspots) of the cases (Table S3). Hence, we can conclude that the probability distributions of the seismic velocity profiles associated with hotspots (red lines, Fig. 1, E-95 test) are specific to a sublist of hotspots.

## 4 DISCUSSION

Anderson (2005) proposed twelve criteria to diagnose a potentially deep mantle origin of hotspots. Among other criteria, Anderson (2005) suggested considering seismic tomography information at various depths, and pleaded for appropriate statistical tests to evaluate the significance of the correlation between hotspots and low-velocity zones in the mantle. Here, we propose a new approach, which is different from previous methods based on a priori assumptions about the character of low velocity anomalies associated with mantle plumes (Boschi et al. 2007, 2008; Marignier et al. 2020; French & Romanowicz 2015). In this study, rather than directly searching for low-velocity regions in the mantle, we propose an alternative statistical approach to chase mantle plumes in global tomography models. We build an objective definition based on a naive Bayesian classification of the mantle seismic structure
SGLOBE

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Probabilities that the velocity profiles beneath a given candidate hotspot location obey the hotspot signature, given the distribution for the 29 slices (E-100, outer ring), and the PCs (E-95, middle ring and E-80, inner ring) for the isotropic part of the model SGLOBE-rani. This is done for the seven locations that are the closest to the hotspot (see main text for details). We show detections with a probability lower (grey) and higher (reddish colours) than 50 per cent, as well as detections with probabilities with significance lower than 95 per cent (dashed colours). 1: Afar; 2: Amsterdam-St Paul; 3: Ascension; 4: Azores; 5: Baja/Guadalupe; 6: Balleny; 7: Bermuda; 8: Bouvet; 9: Bowie; 10: Cameroon; 11: Canary; 12: Cape Verde; 13: Caroline; 14: Comores; 15: Crozet/Pr.Edward; 16: Darfur; 17: Discovery; 18: East Africa; 19: East Australia; 20: Easter; 21: Eifel; 22: Fernando; 23: Galapagos; 24: Great Meteor; 25: Hawaii; 26: Hoggar; 27: Iceland; 28: Jan Mayen; 29: Juan de Fuca/Cobb; 30: Juan Fernandez; 31: Kerguelen; 32: Lord Howe; 33: Louisville; 34: McDonald; 35: Marion; 36: Marquesas; 37: Martin/Trindade; 38: Pitcairn; 39: Raton; 40: Reunion; 41: Samoa; 42: San Felix; 43: St Helena; 44: Socorro; 45: Tahiti/Society; 46: Tasmanid; 47: Tibesti; 48: Tristan; 49: Vema and 50: Yellowstone.</figcaption></figure>

beneath oceans, continents and hotspots. For that purpose, we identify a velocity profile probability density function based on the list of 50 hotspots of Marignier *et al.* (2020). This assumes that the list of hotspots is representative enough, that at least part of the hotspots shares similar seismic velocity profiles, with a low horizontal deflection. This is done on the 6 (SAVANI) to 13 (SGLOBE-rani) varimax components that suffice to capture 95 per cent of the variance in the models, for eight different global seismic tomography models. Based on this classification, we measure a high level of type II error (i.e. the missed detections), but a low level of type I errors (i.e. the false alarm rate).

### 4.1.1 Missed detection or type II error

The naive Bayesian analysis could recover a sublist of hotspots, from 24 per cent (SAVANI) to 52 per cent (SGLOBE-rani). We could not

recover any hotspots for 12 geographical locations: Bermuda, Comores, East Australia, Eifel, Hoggar, Lord Howe, Martin, Raton, San Felix, Socorro, Tibesti and Vema. In French & Romanowicz (2015), San Felix is ‘somewhat resolved’ while none but Comores are considered as ‘primary plume’, and none as ‘clearly resolved’. Vema is identified as ‘likely’ by Marignier *et al.* (2020), but all other ones were found to be ‘unclear’ or ‘unlikely’ in their study. Moreover, the missed detections are at least partly model-dependent: 12 hotspots were recovered in one or two models while being absent in the other ones (Fig. 4 and Table S2). Among them, Cameroon, recovered by S20RTS, and Galapagos, recovered by SGLOBE-rani, are mentioned as ‘clearly resolved’ by French & Romanowicz (2015), while Bouvet, identified by S362WMANI+M, Discovery, identified by SGLOBE-rani and S362WMANI+M, and Fernando, identified by S20RTS, are considered as ‘likely’ by Marignier *et al.* (2020).

As discussed briefly above, this partial recovery of hotspots suggests that at least some of the missed detections are possibly not the
surface expression of mantle plumes or are underlain by very thin
plumes beyond the resolution of current global tomography models.
These results are also influenced by the choices made when building
the tomography models, such as, for example the different data types
used and, importantly, the distinct choices of model regularization
(Table 1). Let us also remind that Marignier et al. (2020) and French
& Romanowicz (2015) classified the plumes in very distinct ways
to each other and to this study. Based on the analysis of the noise
from the highest harmonic degrees, Marignier et al. (2020) deter-
mined the probability that negative velocities are statistically robust,
while French & Romanowicz (2015) established their classification
based on the continuity and amplitude of observed anomalies. For
example, French & Romanowicz (2015) consider plumes as ‘pri-
mary’ if they identify a low-velocity conduit where $\frac{dV_s}{V_s}$ is lower
than 1.5 per cent for most of the depth interval 1000–2800 km. We
must also consider that our method could not detect plumes that
are deflected on a significant portion of their ascension. Hence, it is
not surprising that there are some differences between those studies
and our results. Nevertheless, overall the main hotspots recovered
in this study were also classified as highly likely mantle plumes
in those previous studies; this includes the deflected Pitcairn and
Hawaii plumes, and thus there is a good level of agreement.

4.1.2 False alarm rate or type I error

In our study, we found that while keeping the 29 slices (E-100) in the analysis, the false alarm rate cannot be brought below 5 per cent if the detection threshold for a point is not larger than 0.99. In the E-95 case, a threshold not larger than 0.66 ensures the false alarm rate remains below 5 per cent. This strongly highlights the advantages of using the varimax approach when analysing global tomography models, as also shown by de Viron et al. (2021). It is important to emphasize that the null hypothesis test performed on random sets of surrogate hotspots further demonstrated the specificity of the velocity profiles from the hotspot list.

The false alarms correspond to a larger set of locations mostly
around the hotspots, but occasionally also away from them (e.g.
west of Hawaii, detected by S20RTS, southwest of Cameroon for
GLAD-M25 and west of Juan Fernandes for SGLOBE-rani). This
may support the hypothesis that these detected hotspots are the sur-
face expression of large-scale mantle phenomena and that many
mantle plumes are fat and thus with a thermochemical origin, as
suggested by some previous studies (e.g. French & Romanowicz
2015; Davaille & Romanowicz 2020). This is also consistent with
a possible substantial horizontal deflection of some plumes (Stein-
berger 2000; Nolet et al. 2007; Stockmann et al. 2019). Alterna-
tively, the anomalies may result from the superimposition of several
thinner plumes that cannot be distinguished due to the limited res-
olution of the tomography models as proposed, for example for
the East African region (Chang et al. 2020). Beneath Africa, some
recovered zones are further away from hotspots (e.g. in the south-
ern Africa region) and may correspond to mantle plumes that have
not yet reached the upper mantle, as suggested, for example by
Davaille & Vatteville (2005). Interestingly, the mantle plume zones
identified in our analysis do not appear to surround Africa and Pa-
cific LLSVPs, but are rather within the LLSVPs. This supports the
hypothesis that the LLSVPs may correspond to bundles of ther-
mochemical, fat plumes as suggested by Davaille & Romanowicz
(2020) or that the LLSVPs are made up of compact piles, topped
by a bundle of plumes, rather than to compact, dense piles with
plumes rising from their edges (e.g. Burke et al. 2008). Davaille

& Romanowicz (2020) used a combination of geodynamic con-
straints from fluid mechanics laboratory experiments and seismic
tomography information to support the thermochemical plume bundle idea. Their results were partly based on the analysis of minima in shear wave low velocity in the lower mantle in the tomography model SEMUCB-WM1 along with synthetic resolution tests of various possible structures of LLSVP-like anomalies. Since our statistical analysis of the tomography models considered in this study is based on a completely independent method to the analysis of Davaille & Romanowicz (2020), it provides further alternative support to the hypothesis of LLSVPs consisting of bundles of thermochemical plumes, possibly enriched in denser than average material.

5 CONCLUSION

We provide a quantitative way to look for hotspots using eight
different global seismic tomography models by looking for co-
herency between the list of 50 hotspots established by Marignier
et al. (2020) and the depth distribution of shear velocity anoma-
lies provided by eight global tomography models. The test is based
on a naïve Bayesian clustering method, which is possible thanks
to the varimax PC analysis, providing components that are much
more independent than the original 29 slices of the tomography
models used in our study. The proposed method also assumes that
the plumes are not significantly deflected horizontally. Based on
the 6 (SAVANI) to 13 (SGLOBE-rani) varimax components that
suffice to capture 95 per cent of the model variance, the hotspot
list considered did not allow to establish a generic characteristic,
given the high rate of missed detection, or type II error. Hence, the
list of 50 hotspots and the associated distribution of shear velocity
anomalies do not allow a unique definition of a characteristic veloc-
ity profile that could be associated with hotspots. This supports
the view that different geodynamic processes underlie different
hotspots.

Nevertheless, we establish a probability distribution of the seis-
mic velocity profiles underneath the hotspots and showed that this
distribution does not significantly exist elsewhere, that is it is spe-
cific to a sublist of hotspots. Our Bayesian analysis could recover
from 24 per cent (SAVANI) to 52 per cent (S20RTS) of the hotspot
list. In particular, it recovered eight of the eleven hotspots classified
as ‘primary plume’ by French & Romanowicz (2015) for at least
seven global tomography models, and ten by at least six models.
A noticeable exception is the East African Rift, which obeys the
same profiles as the hotspots for all the tomography models con-
sidered, supporting the paradigm that the origin of the Rift may
come from the deep mantle. Moreover, we found that the mantle
plume regions detected in our study do not surround the African
and Pacific LLSVPs, but are rather within them. This rules out the
idea that LLSVPs correspond to compact, dense piles with mantle
plumes rising from their edges. Instead, our analysis suggests two
possible options that either the LLSVPs: (1) correspond to bundles
of thermochemical mantle plumes or (2) are made up of compact
piles topped by a bundle of plumes.

Our open, freely available codes are modular and easy to use, and
the method can be easily adapted for further analysis, for example
of regional seismic tomography models. Moreover, other tectonic
settings can also be investigated in future work, such as, for example
subduction zones and cratons. In the future, it will be worth applying
the proposed tests when new tomographic models are available, and
when other criteria allow establishing new lists of hotspots.
## ACKNOWLEDGMENTS

AMGF is grateful for support from NERC grant NE/N011791/1 and from the ERC (UPFLOW grant number 101001601). We greatly thank Lapo Boschi, Barbara Romanowicz and Raj Moulik for providing us with files containing their model parametrizations. This paper benefited from valuable comments and suggestions from two anonymous reviewers and editor Ebru Bozdag.

*Author contributions:* OdV and MVC performed the formal analysis, OdV and MVC developed and adapted the methodologies used, OdV made the visualization, and all the authors contributed to writing the paper and to the investigation and validation of the results.

## CONFLICT OF INTEREST

No potential conflict of interest was reported by the authors.

## DATA AVAILABILITY

No new data were generated in support of this research. The tomographic models S20RTS, S40RTS and SEMUCB-WM1 are available directly upon request from their authors; the other five models are publicly available from the IRIS earth model collaboration repository: http://ds.iris.edu/ds/products/emc-earthmodels/, Trabant *et al.* 2012.

## REFERENCES

Anderson, D.L., 2005. Scoring hotspots: the plume and plate paradigms, in *Special Paper 388: Plates, Plumes and Paradigms*, Vol. 388, pp. 31–54, Geological Society of America.

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

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

Boschi, L., Becker, T.W. & Steinberger, B., 2007. Mantle plumes: dynamic models and seismic images: mantle plume tomography, *Geochem. Geo-phys. Geosyst.*, **8**(10), doi:10.1029/2007GC001733.

Boschi, L., Becker, T.W. & Steinberger, B., 2008. On the statistical significance of correlations between synthetic mantle plumes and tomographic models, *Phys. Earth planet. Inter.*, **167**(3-4), 230–238.

Bozdağ, E., Peter, D., Lefebvre, M., Komatitsch, D., Tromp, J., Hill, J., Podhorszki, N. & Pugmire, D., 2016. Global adjoint tomography: first-generation model, *Geophys. J. Int.*, **207**(3), 1739–1766.

Burke, K., Steinberger, B., Torsvik, T.H. & Smethurst, M.A., 2008. Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary, *Earth planet. Sci. Lett.*, **265**(1-2), 49–60.

Chang, S.-J., Ferreira, A. M.G., Ritsema, J., van Heijst, H.J. & Woodhouse, J.H., 2015. Joint inversion for global isotropic and radially anisotropic mantle structure including crustal thickness perturbations: radially anisotropic mantle structure, *J. geophys. Res.*, **120**(6), 4278–4300.

Chang, S., Kendall, E., Davaille, A. & Ferreira, A. M.G., 2020. The evolution of mantle plumes in East Africa, *J. geophys. Res.*, **125**(12), doi:10.1029/2020JB019929.

Cottaar, S. & Lekic, V., 2016. Morphology of seismically slow lower-mantle structures, *Geophys. J. Int.*, **207**(2), 1122–1136.

Cottaar, S. & Romanowicz, B., 2012. An unusually large ULVZ at the base of the mantle near Hawaii, *Earth planet. Sci. Lett.*, **355-356**, 213–222.

Courtillot, V., Davaille, A., Besse, J. & Stock, J., 2003. Three distinct types of hotspots in the Earth’s mantle, *Earth planet. Sci. Lett.*, **205**(3-4), 295–308.

Davaille, A. & Romanowicz, B., 2020. Deflating the LLSVPs: bundles of mantle thermochemical plumes rather than thick stagnant ‘piles’, *Tectonics*, **39**(10), doi:10.1029/2020TC006265.

Davaille, A. & Vatteville, J., 2005. On the transient nature of mantle plumes: mantle plumes, *Geophys. Res. Lett.*, **32**(14), doi:10.1029/2005GL023029.

de Viron, O., Van Camp, M., Grabkowiak, A. & Ferreira, A. M.G., 2021. Comparing global seismic tomography models using varimax principal component analysis, *Solid Earth*, **12**(7), 1601–1634.

Duncombe, J., 2019. The unsolved mystery of the Earth blobs, *EOS, Trans. Am. geophys. Un.*, **100**, doi:10.1029/2019EO117193.

Durand, S., Debayle, E., Ricard, Y., Zaroli, C. & Lambotte, S., 2017. Confirmation of a change in the global shear velocity pattern at around 1000 km depth, *Geophys. J. Int.*, **211**(3), 1628–1639.

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

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

Foulger, G.R. & Natland, J.H., 2003. Is ‘hotspot’ volcanism a consequence of plate tectonics?, *Science*, **300**(5621), 921–922.

French, S.W. & Romanowicz, B., 2015. Broad plumes rooted at the base of the Earth’s mantle beneath major hotspots, *Nature*, **525**(7567), 95–99.

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

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

Kay, S.M., 1993. *Fundamentals of Statistical Signal Processing*, Prentice Hall Signal Processing Series, Prentice-Hall PTR.

Koelemeijer, P., Deuss, A. & Ritsema, J., 2017. Density structure of Earth’s lowermost mantle from Stoneley mode splitting observations, *Nat. Commun.*, **8**(1), doi:10.1038/ncomms15241.

Koller, D. & Friedman, N., 2009. *Probabilistic Graphical Models: Principles and Techniques, Adaptive Computation and Machine Learning*, MIT Press.

Kumagai, I., Davaille, A., Kurita, K. & Stutzmann, E., 2008. Mantle plumes: Thin, fat, successful, or failing? Constraints to explain hot spot volcanism through time and space, *Geophys. Res. Lett.*, **35**(16), doi:10.1029/2008GL035079.

Kustowski, B., Ekström, G. & Dziewoński, A.M., 2008. Anisotropic shear-wave velocity structure of the Earth’s mantle: a global model, *J. geophys. Res.*, **113**(B6), doi:10.1029/2007JB005169.

Lau, H. C.P., Mitrovica, J.X., Davis, J.L., Tromp, J., Yang, H.-Y. & Al-Attar, D., 2017. Tidal tomography constrains Earth’s deep-mantle buoyancy, *Nature*, **551**(7680), 321–326.

Lei, W. *et al.*, 2020. Global adjoint tomography-model GLAD-M25, *Geophys. J. Int.*, **223**(1), 1–21.

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

Lin, S.-C. & van Keken, P.E., 2006. Dynamics of thermochemical plumes: 2. Complexity of plume structures and its implications for mapping mantle plumes, *Geochem. Geophys. Geosyst.*, **7**(3), doi:10.1029/2005GC001072.

Marignier, A., Ferreira, A. M.G. & Kitching, T., 2020. The probability of mantle plumes in global tomographic models, *Geochem. Geophys. Geosyst.*, **21**(9), doi:10.1029/2020GC009276.

Marzocchi, W. & Mulargia, F., 1993. Patterns of hot spot volcanism, *J. geophys. Res.*, **98**(B8), 14 029–14 039.

Montelli, R., 2004. Finite-frequency tomography reveals a variety of plumes in the mantle, *Science*, **303**(5656), 338–343.

Morgan, W.J., 1971. Convection plumes in the lower mantle, *Nature*, **230**(5288), 42–43.

Moulik, P. & Ekström, G., 2014. An anisotropic shear velocity model of the Earth’s mantle using normal modes, body waves, surface waves and long-period waveforms, *Geophys. J. Int.*, **199**(3), 1713–1738.

Nolet, G., Allen, R. & Zhao, D., 2007. Mantle plume tomography, *Chem. Geol.*, **241**(3-4), 248–263.
O'Connor, J.M., Jokat, W., le Roex, A.P., Class, C., Wijbrans, J.R., Keßling, S., Kuiper, K.F. & Nebel, O., 2012. Hotspot trails in the South Atlantic controlled by plume and plate tectonic processes, *Nat. Geosci.*, **5**(10), 735–738.

Ritsema, J., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland, *Science*, **286**(5446), 1925–1928.

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

Sleep, N.H., 1990. Hotspots and mantle plumes: some phenomenology, *J. geophys. Res.*, **95**(B5), 6715, doi:10.1029/JB095iB05p06715.

Steinberger, B., 2000. Plumes in a convecting mantle: models and observations for individual hotspots, *J. geophys. Res.*, **105**(B5), 11 127–1152.

Stockmann, F., Cobden, L., Deschamps, F., Fichtner, A. & Thomas, C., 2019. Investigating the seismic structure and visibility of dynamic plume models with seismic array methods, *Geophys. J. Int.*, **219**(Supplement.1), S167–S194.

Suetsugu, D. et al., 2009. South Pacific mantle plumes imaged by seismic observation on islands and seafloor: imaging of south pacific mantle plumes, *Geochem. Geophys. Geosyst.*, **10**(11), doi:10.1029/2009GC002533.

Trabant, C., Hutko, A.R., Bahavar, M., Karstens, R., Ahern, T. & Aster, R., 2012. Data products at the IRIS DMC: stepping stones for research and other applications, *Seismol. Res. Lett.*, **83**(5), 846–854.

Trampert, J., 2004. Probabilistic tomography maps chemical heterogeneities throughout the lower mantle, *Science*, **306**(5697), 853–856.

Virtanen, P. et al., 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python, *Nat. Methods*, **17**(3), 261–272.

Wang, Y. & Li, M., 2021. The interaction between mantle plumes and lithosphere and its surface expressions: 3-D numerical modelling, *Geophys. J. Int.*, **225**(2), 906–925.

Wilson, J., 1965. Evidence from ocean islands suggesting movement in the Earth, *Phil. Trans. R. Soc. Lond.*, A., **258**(1088), 145–167.

Zechmann, E., 2019. Make Icosahedron (https://www.mathworks.com/matlabcentral/fileexchange/19169-make-icosahedron), MATLAB Central File Exchange. Retrieved October 30, 2022.

Zhao, D., 2004. Global tomographic images of mantle plumes and subducting slabs: insight into deep Earth dynamics, *Phys. Earth planet. Inter.*, **146**(1-2), 3–34.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** The four sets of geographical locations considered in this study: red: hotspots (Locations A); white: gridpoints in a 10° range from each hotspot (Locations B); blue: oceans (Locations C); green: continents (Locations D).

**Figure S2.** Same as in Fig. 6, but for the tomography model SA-VANI.

**Figure S3.** Same as in Fig. 6, but for the tomography model S20RTS.

**Figure S4.** Same as in Fig. 6, but for the tomography model S40RTS.

**Figure S5.** Same as in Fig. 6, but for the tomography model SEMUCB-WM1.

**Figure S6.** Same as in Fig. 6, but for the tomography model S362WMANI+M.

**Figure S7.** Same as in Fig. 6, but for the tomography model SEIS-GLOB2.

**Figure S8.** Same as in Fig. 6, but for the tomography model GLAD-M25.

**Table S1.** The 50 continuous probability distributions tested in this study. The distributions are proposed in the SciPy.stats library (Virtanen et al. 2020).

**Table S2.** List of the hotspots recovered by the eight S-wave isotropic global tomography models used in this study, applying the P(50,95) criterion for the principal component analysis experiment E-95. This corresponds to at least one reddish slice of the middle pie charts on Figs 6 and S2–S8. The hotspots that are underlined indicate the zones classified as ‘Primary Plume’ by French & Romanowicz (2015) while the bold ones correspond to the ‘Highly likely’ mantle plumes of Marignier et al. (2020).

**Table S3.** Percentage of cases for which we obtain 0, 1, 2 or 3 false alarms (FA) for the eight global tomography models considered in this study, using 50 surrogate hotspots to perform the statistical analysis (see the main text for further details).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.