## Abstract

The thermochemical structure of lithospheric and asthenospheric mantle exert primary controls on surface topography and volcanic activity. Volcanic rock compositions and mantle seismic velocities provide indirect observations of this structure. Here, we compile and analyze a global database of the distribution and composition of Neogene-Quaternary intraplate volcanic rocks. By integrating this database with seismic tomographic models, we show that intraplate volcanism is concentrated in regions characterized by slow upper mantle shear-wave velocities and by thin lithosphere (i.e. <100 km). We observe a negative correlation between shear-wave velocities at depths of 125–175 km and melt fractions inferred from volcanic rock compositions. Furthermore, mantle temperature and lithospheric thickness estimates obtained by geochemical modeling broadly agree with values determined from tomographic models that have been converted into temperature. Intraplate volcanism often occurs in regions where uplifted (but undeformed) marine sedimentary rocks are exposed. Regional elevation of these rocks can be generated by a combination of hotter asthenosphere and lithospheric thinning. Therefore, the distribution and composition of intraplate volcanic rocks through geologic time will help to probe past mantle conditions and surface processes.

## Introduction

Mantle convection drives plate tectonics, redistributes heat, cycles chemical species, and generates dynamic topography at the Earth’s surface. Constraining spatio-temporal changes in the thermochemical structure of the mantle will help to refine our understanding of these interlinked phenomena. This significant challenge requires careful integration of diverse observations. Two well-established methodologies for probing the thermochemical state of the mantle are igneous geochemistry and seismic tomography. It is particularly useful to analyze the relationship between igneous geochemistry and upper mantle structure away from complications associated with active plate boundaries in order to understand the interaction between mantle composition, temperature and melting through geologic time.

An important difficulty is that composition of partial melts within the mantle and shear-wave velocity of mantle rocks are both sensitive to the combined effects of temperature and chemical composition. The uppermost mantle is subdivided into a mechanically strong lithosphere and a weak underlying asthenosphere that is ~100–200 km thick. For a given mantle source composition, the depth and degree of melting are principally controlled by a combination of asthenospheric temperature and lithospheric thickness^{1}. The extent of melting increases with increasing potential temperature, *T*_{p}, and decreasing pressure so that elevated asthenospheric temperature and/or thinner lithosphere produce greater volumes of melt. Given suitable assumptions, the location, volume and composition of volcanic rocks can be used to analyze the thermal structure of the uppermost mantle^{2,3,4}. Nevertheless, it is important to bear in mind that the distribution and composition of volcanic rocks are also significantly influenced by mantle composition, by the geometry of mantle flow, and by the interaction between melts and surrounding rocks as they ascend to the surface^{5,6,7,8,9,10,11}.

Global seismic tomographic models demonstrate that shear-wave velocity anomalies, Δ*V*_{s}, occur throughout the upper mantle. Although it is agreed that shear waves propagate slower through warm mantle and faster through cold mantle, it is also clear that *V*_{s} varies in a non-linear fashion with pressure and temperature^{12,13}. Furthermore, shear-wave speed is also sensitive to mantle composition and to the presence of interstitial melt^{14,15}. By analyzing the distribution and composition of volcanic rocks in conjunction with shear-wave velocity anomalies, it is possible to determine temperature variations within the upper mantle. A related approach has been successfully used along mid-oceanic ridges where lithospheric thickness can be assumed to be negligible^{16,17}.

Here, we investigate a more general problem by analyzing intraplate regions, where both asthenospheric temperature and lithospheric thickness vary. Our principal aim is to quantity the relationship between the composition of volcanic rocks and mantle seismic velocities with a view to constraining the size, extent and surface expression of putative thermal anomalies. We are less concerned with the more difficult problem of how these anomalies are generated in the first place. First, we examine correlations between the distribution and composition of Neogene-Quaternary intraplate volcanism, Δ*V*_{s}, and lithospheric thickness. Secondly, we estimate potential temperature and lithospheric thickness using a combination of geochemical and seismologic techniques. Finally, we scrutinize the link between intraplate volcanism and the distribution of emergent marine sedimentary deposits, which are a tangible manifestation of dynamic topographic uplift. Although our primary focus is on relatively youthful volcanic rocks, we are hopeful that our quantitative framework will provide helpful tools for interrogating the Phanerozoic record of these processes.

## Results

### Spatial distribution of intraplate magmatism

We have compiled a global database that consists of >20, 000 geochemical analyses of Neogene-Quaternary intraplate volcanic rocks (Database 1; Fig. 1; and Supplementary Tables 1 and 2). Since lithospheric plates translate across the globe, and since sub-plate mantle circulation evolves as a function of time and space, we deliberately restrict our study to samples that are <10 Ma and located <400 km whence they were erupted, based upon present-day plate speeds^{18}. Global surface-wave tomographic models have a horizontal resolution of 200–600 km, which means that these restrictions ensure co-location of intraplate volcanic rocks and pertinent observations of sub-surface structure^{19,20,21,22}. The majority of analyses are of mafic samples that are located far from active plate boundaries, although we have included analyses from locations adjacent to several extensional provinces and subduction zones that exhibit intraplate-like compositions (e.g., western North America, Anatolia, East African Rift). We have carried out a literature review to identify and remove samples pertaining to subduction zone processes (see Supplementary Information). Database 1 shows that most intraplate volcanism is concentrated within bands located on continental lithosphere away from cratonic regions where lithosphere >200 km thick can occur. The four most extensive bands reach along the length of western North America, along the eastern seaboard of Australia, throughout eastern China and southeast Asia, and across Europe through Anatolia and Arabia down to southern Africa. The smaller number of database entries from the oceanic realm compared with the continents reflects sampling bias toward sub-aerial locations.

If the distribution of volcanic rocks is compared with the pattern of upper mantle shear-wave velocities, it is clear that intraplate volcanism is concentrated within regions where negative shear-wave velocity anomalies occur at depths of 150 ± 25 km. Figure 1a shows a striking visual association for the SL2013sv tomographic model (see “Methods”^{21}). This global vertical shear-wave velocity model exploits body and surface waves that include both fundamental and higher modes with periods of 11–450 s. It is important to emphasize that similar results have been obtained for a range of other tomographic models (Supplementary Fig. 1). This global relationship is consistent with regional studies, which show that Δ*V*_{s} correlates with geochemical indicators of melt fraction within intraplate settings^{23,24,25}. Intraplate volcanism is almost completely absent in regions where Δ*V*_{s} is positive. There is a similarly compelling spatial relationship between the distribution of intraplate volcanic rocks and lithospheric thickness (Fig. 1b). To construct this map, shear-wave velocities from the SL2013sv model are converted into temperature by exploiting a global calibration between *V*_{s} and *T*, based upon the oceanic plate cooling model^{26}. Here, we utilize the revised and modified *V*_{s}-to-*T* parameterization described by ref. ^{27}, which is based upon an empirical anelastic parameterization (see “Methods”^{13}). Lithospheric thickness is estimated by tracking the depth to the 1175 °C isothermal surface^{27}. This isothermal surface is chosen since it provides a good fit to the peak change in the orientation of azimuthal anisotropy, which is characteristic of the transition between rigid lithosphere and convecting asthenosphere within oceanic regions^{28,29}. In the continental realm, it is clear that intraplate volcanic rocks are concentrated within regions where the lithosphere is <100 km thick.

These visual associations can be quantified by formally examining relationships between, Σ*A*_{I}, the cumulative areal distribution of binned intraplate volcanic samples, Δ*V*_{s}, and lithospheric thickness (see “Methods”). Eighty-nine percent of Σ*A*_{I} occurs in regions underlain by negative values of Δ*V*_{s} whose areal distribution only constitutes 61% of global area (Fig. 2a). Ninety-five percent of Σ*A*_{I} occurs in regions underlain by continental lithosphere <100 km thick whose areal distribution constitutes 62% of global area (Fig. 2b). These spatial associations are consistent for a range of tomographic and lithospheric thickness models (Supplementary Fig. 1). A Kolmogorov–Smirnov test is used to gauge the likelihood of these relationships being merely coincidental (“Methods”^{30}). The probability that co-location of intraplate samples, negative values of Δ*V*_{s}, and anomalously thin continental lithosphere arises from chance is <1 in 10^{100}.

### Relationships between tomography and geochemistry

The ratio of La and Sm concentrations in mafic igneous rocks is often used to track mantle melting^{31,32,33}. Since La is less compatible within the mantle source compared with Sm, the La/Sm content of a melt decreases as melt fraction increases. The relationship between La/Sm and Δ*V*_{s} is a useful way to explore the role that sub-plate thermal anomalies play in generating intraplate volcanism. In order to assess this relationship, we sub-divide our global database of La/Sm and Δ*V*_{s} measurements into 1° bins and determine the Pearson product-moment correlation coefficient, *R*, between these measurements (Fig. 2c; “Methods”). Since 96% of areal bins occur where lithospheric thickness is <100 km, we use the average value of Δ*V*_{s} at 150 ± 25 km depth to represent asthenospheric mantle. Mineral loss or accumulation can significantly alter the composition of a lava sample away from that of the original mantle melt. Removal or addition of mineral cargo will decrease or increase MgO concentration of the melt, respectively. Thus, MgO content is considered to be a useful indicator of either process. Here, we have excised all samples that have MgO contents <9 wt% and >14.5 wt% from Database 1 to ensure that we only retain samples that are close to the composition of primitive mantle melts. For this filtered version of Database 1, only bins with >5 samples that have been analyzed for La and Sm concentrations are exploited in order to mitigate the influence of local compositional heterogeneity.

There is a positive correlation between La/Sm ratios and Δ*V*_{s} values for the SL2013sv tomographic model at a depth of 150 ± 25 km (*R* = 0.65 ± 0.07; Fig. 2c). Thus, igneous rocks with lower La/Sm ratios, which are indicative of higher melt fractions, coincide with regions with lower values of Δ*V*_{s} at a depth of 150 ± 25 km, which are indicative of higher temperatures at this depth. The value of this correlation does not significantly change if different methods of binning and filtering are employed, if locations adjacent to mid-oceanic ridges are excluded, or if regions where plate subduction has recently occurred are excised (Supplementary Figs. 2–4). In Fig. 2d, correlations between La/Sm and Δ*V*_{s} are shown as a function of depth for the SL2013sv, CAM2016Vsv, SEMUCB-WM1 and S40RTS global tomographic models^{19,21,22,34}. In each case, the value of *R* is greatest between 100 and 200 km depth. For three of these models, the value of *R* sharply decreases at depths of >250 km. Since both La/Sm and Δ*V*_{s} are expected to decrease as mantle temperature increases, these observations suggest that intraplate volcanism is sensitive to temperature variations within the uppermost mantle.

Notwithstanding this correlation, melt chemistry and upper mantle shear-wave velocity can also be influenced by changes in lithospheric thickness, by the composition of the mantle source region, and by depth of the base of the melt column, which is controlled by a combination of mantle temperature and composition^{11,14,32,35}. Although we define the lithosphere-asthenosphere boundary to be the 1175 °C isothermal surface, the base of the actual thermal boundary layer probably extends to greater depths. Global surface-wave tomographic models have a vertical resolution of 25–50 km^{19,21,22,34}. Variations in thickness of the thermal boundary layer coupled with seismic resolution may, therefore, influence values of Δ*V*_{s} estimated within the asthenospheric mantle. Consequently, any correlation between La/Sm and Δ*V*_{s} could also be influenced by lithospheric thickness variations. At depths where Δ*V*_{s} values are affected by lithospheric thickness, a positive correlation between La/Sm and Δ*V*_{s} is expected. In Fig. 2e, correlations between lithospheric thickness and Δ*V*_{s} are shown for the suite of tomographic models as a function of depth. Beneath intraplate volcanic provinces, lithospheric thickness correlates significantly with Δ*V*_{s} at depths ≤100 km. At depths ≥150 km, this correlation rapidly becomes insignificant (i.e., where *R* < 0.28, which is the minimal correlation distinguishable from zero at a significance level of 0.001 for a sample size of 134; see “Methods”). Clear correlations between La/Sm and Δ*V*_{s} (i.e., *R* = 0.5–0.7) are observed at depths of 150–200 km where the influence of lithospheric thickness changes is deemed to be negligible (Fig. 2d). Thus, it is reasonable to assume that correlations between Δ*V*_{s} and La/Sm at these depths are principally controlled by temperature variations within the asthenospheric mantle.

Composition of the mantle source region affects the values of both La/Sm and Δ*V*_{s} since melting occurs to a lesser degree within a depleted source region, depleted mantle has lower initial La/Sm values than fertile mantle, and shear waves propagate faster through depleted mantle than through fertile mantle^{14,36}. To determine how the thermochemical structure of the mantle controls melt composition, we carried out principal component analysis on a suite of eight incompatible element concentrations taken from our filtered and binned version of Database 1 (see “Methods”). We then calculated correlations for each of these principal components with respect to geochemical and geophysical indicators of mantle temperature and depletion. The first principal component, *P*_{1}, accounts for 79 ± 1% of the variance within our filtered and binned version of Database 1. The weightings of *P*_{1} are ~0.4 for all elements with the exception of Yb (Fig. 3a). As melt fraction increases, incompatible element concentrations within the melt decrease. Since the elemental data within this analysis is both mean-centered and variance-scaled, the range of values for each element between the lowest and highest melt fractions should be approximately equal. Therefore, an equal weighting over a wide range of incompatible elements for *P*_{1} suggests that *P*_{1} is primarily sensitive to melt fraction variation.

As melt fraction increases, incompatible elemental concentrations, and thus *P*_{1}, will decrease within the melt. Alternatively, since these elemental concentrations can vary as a function of source depletion, *P*_{1} may instead be primarily sensitive to mantle composition. To investigate these contrasting hypotheses, we generate correlations for *P*_{1} with respect to La/Sm, Δ*V*_{s} and *ε*Nd. Note that the *ε*Nd value of igneous rocks is a commonly used proxy for mantle depletion since depleted mantle has a higher value of *ε*Nd than fertile mantle^{6,36}. If variations in *P*_{1}, La/Sm, Δ*V*_{s} and *ε*Nd were primarily controlled by mantle depletion, we would expect *P*_{1} to positively correlate with La/Sm, but to negatively correlate with Δ*V*_{s} and *ε*Nd. However, *P*_{1} positively correlates with both La/Sm and Δ*V*_{s}, and it negatively correlates with *ε*Nd (*R* = 0.87, 0.59 and −0.50, respectively; Supplementary Fig. 5a–c). If Δ*V*_{s} variations were primarily dependent upon mantle composition then the most depleted mantle would be associated with the fastest shear-wave speeds (i.e., *ε*Nd and Δ*V*_{s} would positively correlate). Since the opposite relationship is observed, it is less likely that the correlation between La/Sm and Δ*V*_{s} is controlled by source compositional variations. Instead, we conclude that *P*_{1} is sensitive to melt fraction variations. This result is anticipated since temperature changes are known to have a greater effect upon *V*_{s} than mantle composition^{12,14}. We infer that upper mantle temperature variations are a dominant control on the positive relationship between melt fraction (i.e., *P*_{1} and La/Sm values) and Δ*V*_{s}. Note that, while fractional crystallization can have a modest influence on the value of La/Sm, *P*_{1} does not correlate with MgO concentrations and the observed range of La/Sm values is too great to be controlled by fractional crystallization alone (Supplementary Fig. 6).

Melt fraction, and therefore both La/Sm and *P*_{1}, vary as a function of both asthenospheric temperature and lithospheric thickness. We can disentangle the relative contributions of these quantities by using the second principal component, *P*_{2}, which accounts for 14% of the variance and is dominated by variations in Yb. Compatibility of Yb within the mantle increases at depths greater than 60–80 km where garnet becomes stable at the expense of spinel^{32,37}. Since Yb is far more compatible in garnet-bearing mantle relative to spinel-bearing mantle, Yb concentrations within a melt are more sensitive to the depth at which melting occurs than to changes in melt fraction. Therefore, we can use *P*_{2} as a proxy for comparing the respective contributions of melting within the spinel and garnet stability fields. The contribution of melting within the garnet stability field relative to the spinel field is greater beneath thicker lithosphere or in regions where elevated asthenospheric temperature acts to deepen the onset of melting. Consequently, *P*_{2} can act as a proxy for lithospheric thickness, provided that asthenospheric *T*_{p} remains constant. Alternatively, *P*_{2} can act as a proxy for asthenospheric *T*_{p}, provided that lithospheric thickness is fixed. Δ*V*_{s} at depths of 150 ± 25 km is sensitive to asthenospheric *T*_{p} variations, but insensitive to lithospheric thickness (Fig. 2e). We do not observe a correlation between *P*_{2} and Δ*V*_{s} at 150 ± 25 km depths (*R* = −0.07; Supplementary Fig. 5e). This absence of correlation implies that asthenospheric temperature alone does not control the depth range over which melting occurs and that lithospheric thickness variation must play an important role. Consequently, these changes of lithospheric thickness will influence the observed values of La/Sm and, therefore, the inferred final melt fractions. However, if lithospheric thickness variations were the dominant control on melt fraction, we would not observe a correlation between La/Sm and Δ*V*_{s} at depths of 150 ± 25 km (Fig. 2d). We believe that this analysis demonstrates that asthenospheric temperature variations are the primary control on values of La/Sm but that lithospheric thickness changes exert an important secondary effect. Note that this conclusion partially contradicts previous studies, which identify a dominant lithospheric thickness signal by comparing melt chemistry with plate age in oceanic regions (e.g., refs. ^{11,38}). These studies were primarily concerned with major element measurements, which we do not consider here. Such opposing views are reconcilable since lithospheric thickness may have a more significant impact upon major element composition than *T*_{p} if melts stall and thermally re-equilibrate at or near the base of the lithosphere^{25}.

The third principal component, *P*_{3}, accounts for 7 ± 1% of the variance and is primarily sensitive to variations in Ba, Nb, and K. Subduction zone melting is commonly enriched in Ba and K but depleted in Nb. Partial melting of metasomatized lithospheric mantle can exhibit anomalously high, or anomalously low, normalized concentrations of K and Ba relative to Nb^{39,40}. The low variance of *P*_{3} suggests that there are minor contributions either from prior subduction zone magmatism, from melting of metasomatized lithospheric mantle, or from both. *P*_{3} has low variance and does not correlate either with La/Sm or with Δ*V*_{s} at depths of 150 ± 25 km (Supplementary Fig. 5). Thus, the degree of contamination by lithospheric melts is evidently not the primary control for incompatible element concentrations in intraplate mafic igneous rocks.

These results are corroborated by analysis of synthetic principal components that are calculated for a suite of geochemical models. One-hundred and fifty synthetic trace element distributions were calculated by randomly varying values of *T*_{p}, lithospheric thicknesses, and mantle compositions within fixed limits (i.e., 1250–1450 °C, 35–75 km, and primitive/depleted mantle, respectively; see “Methods”). This procedure is repeated 99 times. The first synthetic principal component, *M*_{1}, is similar to *P*_{1} both in terms of elemental weightings and variance (Fig. 3d). *M*_{1} strongly correlates with *T*_{p} and does not significantly correlate either with lithospheric thickness or with mantle depletion (*R* = − 0.85 to − 0.72, *R* = 0.05 to 0.42, *R* = −0.27 to 0.21, respectively; Supplementary Fig. 7). It is reasonable to conclude that both *M*_{1} and *P*_{1} are primarily sensitive to changes in *T*_{p}. *M*_{2} accounts for 8 ± 4% of the variance and, like *P*_{2}, it is dominated by variations in Yb (Fig. 3e). For the majority of model runs, *M*_{2} correlates with lithospheric thickness but does not correlate either with *T*_{p} or with mantle depletion (*R* = −0.81 to −0.26, *R* = −0.37 to 0.24, −0.22 to 0.18, respectively; Supplementary Fig. 7). Correlation between *M*_{2} and lithospheric thickness may indicate that *P*_{2} is also primarily sensitive to lithospheric thickness variations. *M*_{3} has a similar distribution of weightings as *P*_{3} but not as strong weighting for Ba, K, and Nb. As expected, *M*_{3} does not consistently correlate with *T*_{p}, lithospheric thickness or depletion. These observations add strength to the idea that *P*_{3} represents lithospheric contamination, which is not accounted for by this synthetic model. Note that conclusions based on these synthetic models depend upon the assumption that asthenospheric *T*_{p}, lithospheric thickness, and mantle depletion are uncorrelated with respect to each other on Earth.

### Calculating asthenospheric temperatures

The global distribution of intraplate magmatism together with a positive correlation between La/Sm and Δ*V*_{s} at a depth of 150 ± 25 km suggest that intraplate magmatism is principally associated with elevated asthenospheric temperatures beneath thin lithosphere. This inference can be quantitatively investigated by analyzing relationships between the chemical composition of volcanic rocks, shear-wave velocity anomalies and asthenospheric potential temperature, *T*_{p}. An important goal is to compare geochemical and geophysical estimates of *T*_{p}.

First, the chemical composition of intraplate volcanic rocks is used to constrain the depth and degree of mantle melting. Here, we exploit a geochemical inverse modeling strategy, which minimizes the misfit between observed and calculated concentrations of rare earth elements (REE) by systematically varying melt fraction as a function of depth (“Methods”^{2,41}). In this polybaric model, REE concentrations are calculated by integrating instantaneous melt compositions along isentropic melting paths. These melting paths are determined by a combination of potential temperature of the mantle source region, *T*_{p}, and lithospheric thickness, *a*. The top of the melting column is, by definition, the base of the lithosphere. Beneath the lithosphere, melt fraction as a function of depth is controlled by a combination of *T*_{p} and a chosen melting parameterization. Here, we exploit the hydrous lherzolitic melting model described by Katz et al.^{35}. Note that some model parameters have been updated to honor experimental constraints that have subsequently been obtained (“Methods”^{42}). This melting model requires a *T*_{p} value of 1312 °C in order to generate the average crustal thickness at a mid-oceanic ridge (i.e., 6.9 km^{29}). We assume that this value represents the potential temperature of ambient asthenospheric mantle. For each volcanic province, we systematically vary *T*_{p} and *a* from 1250 to 1550 °C and from 30 to 80 km, respectively. For each combination of *T*_{p} and *a*, the misfit between observed and calculated compositions is measured. In this way, the values of *T*_{p} and *a*, which yield the best fit, are identified. Note that the quality of fit between observed and calculated REE concentrations depends upon the relative contributions of melting within the spinel- and garnet-bearing peridotite stability fields. Combinations of values of *T*_{p} and *a* are deemed acceptable provided that the residual misfit between observed and calculated REE concentrations is not >1.5 times the smallest residual misfit value for any given model run. Volcanic provinces that exhibit a greater range of REE concentrations can generally be fitted within acceptable limits using a larger number of *T*_{p} and *a* combinations.

Figure 4 shows inverse modeling results for the Galápagos volcanic islands located upon the Nazca plate and for the Haruj volcanic province located within north Africa. In each case, an optimal distribution of melt fraction as a function of depth is calculated by a grid search in which only *T*_{p} and *a* are varied. Mantle depletion and water content are set by the value of *ε*Nd (“Methods”). This straightforward approach enables trade-off between values of *T*_{p} and *a* to be clearly identified. For the Galápagos archipelago, we obtain an excellent fit between observed and calculated REE concentrations for *T*_{p} = 1366 ± 15 °C and *a* = 47 ± 5 km, which is within ~30 °C and ~5 km of previous estimates (Fig. 4a, b^{24}). The misfit function shows that there is minimal trade-off between *T*_{p} and *a* and that the value of *a* could be smaller but not significantly larger (Fig. 4c). For the Haruj province, we obtain a satisfactory fit for *T*_{p} = 1367 ± 28 °C and *a* = 56 ± 2 km, in agreement with previous estimates (Fig. 4d, e^{43}). The misfit function shows that there is a negative trade-off between *T*_{p} and *a*, which means that a slightly thicker lithosphere with hotter asthenospheric mantle could fit the observations equally well (Fig. 4f). In both examples, a combination of elevated asthenospheric temperature and anomalously thin lithosphere is required, which is principally a consequence of the constraints imposed by low REE concentrations and by the depth of the spinel-garnet phase transition, respectively.

Inverse modeling has been carried out for intraplate volcanic analyses from our filtered and binned version of Database 1 (Supplementary Fig. 9). In this way, we calculate values of *T*_{p} and *a* for every 1° bin in locations where intraplate volcanism occurs. These values can be directly compared with asthenospheric temperatures and with plate thicknesses obtained from tomographic models that have been converted into temperature. Here, we compare geochemical estimates of *T*_{p} with values determined at 150 ± 25 km depths for the SL2013sv model (“Methods”^{21,27}). Note that the geochemical inverse model and the *V*_{s}-*T* conversion scheme assume hydrous and anhydrous mantle sources, respectively. Consequently, a slightly higher value of ambient potential temperature, *T*_{p} = 1333 °C, would be required to generate the average oceanic crustal thickness for the *V*_{s}-*T* scheme^{27,29}. Globally, we observe a positive correlation of *R* = 0.54 between geochemical and tomographic estimates of Δ*T*_{p}, which is the temperature anomaly with respect to ambient asthenospheric mantle (Fig. 5a).

We acknowledge that limitations of both geochemical and tomographic methodologies could give rise to Δ*T*_{p} discrepancies (e.g., refs. ^{14,26,43,44}). Global tomographic models are considerably damped and smoothed. Ray path coverage of the upper mantle is variable, spatial resolution is limited to ~200–600 km, and lateral smearing of fast velocities can sometimes occur adjacent to thick cratonic lithosphere. These spatial resolution issues could account for the small number of intraplate volcanic regions that fall on regions with positive values of Δ*V*_{s} at depths of 150 ± 25 km. Calibration of the *V*_{s}–*T* parametrization is predicated upon the validity of the plate cooling model since extrapolation of elastic and anelastic behavior determined from laboratory experiments might not be directly applicable to damped tomographic models^{26}. The presence of melt within the mantle may reduce *V*_{s} as a result of poroelastic effects that are not accounted for within this parameterization^{15}. However, the amount of melt retained within the mantle is probably very small (i.e., ~0.1%), especially at a depth of 150 km, which in most regions is beneath the peridotite solidus^{35,45}. Since the reduction of *V*_{s} caused by poroelastic effects should scale as a function of melt fraction, these effects can be considered to be negligible^{46}. The parameterization does not account for compositional variations, which means that temperature estimates for continental regions, especially depleted cratonic cores, may involve additional uncertainties. In continental regions with thin lithosphere (i.e., ≲75 km), additional uncertainties can arise as a result of inaccuracies in the parameterization of crustal structure. These inaccuracies can cause “bleeding” of slow crustal velocities into the uppermost mantle that are interpreted as hotter temperatures, which yield underestimates of lithospheric thickness^{26}. Nevertheless, temperature profiles calculated using this tomographic approach provide acceptable fits to continental geothermal profiles derived from xenolith thermobarometric observations for locations with a range of lithospheric thicknesses (50–200 km^{27}). The observed negative correlation between *ε*Nd and Δ*V*_{s} at depths of 150 ± 25 km also indicates that, in regions of thin lithosphere, shear-wave velocities are controlled by thermal anomalies rather than by compositional variations associated with depletion and enrichment (Supplementary Fig. 5). Within the oceanic realm, where lithospheric structure and composition are broadly uniform, the uncertainty of tomographically derived lithospheric thickness estimates is probably ±30 km^{47}.

Geochemical inverse modeling of intraplate volcanic rocks necessarily depends upon a simplified treatment of mantle source composition, upon the depth of the spinel-garnet phase transition, upon experimentally determined partition coefficients, and upon details of the melting model^{2}. Changes in values of input parameters inevitably affect absolute values of *T*_{p} and *a*. For example, increasing the depth to the spinel-garnet transition by 5 km, decreasing water content of the mantle by 0.01 wt%, or varying mantle composition from depleted to primitive mantle affects calculated values of *T*_{p} and *a* by +15 °C and +5 km, by +5–10 °C and +2 km, and by +20 °C and −3 km, respectively^{43}. Furthermore, the inverse model assumes a uniform lherzolitic source and does not consider harzburgitic and pyroxenitic lithologies. For example, the presence of harzburgite and pyroxenite within the mantle can lead to over- and underestimates of *T*_{p}, respectively, if a purely peridotitic model is implemented^{48}. These sources of uncertainty could be responsible for low temperature values recorded in some provinces and for the minor, but systematic, offset between geochemical and tomographic temperature estimates (Fig. 5a). The correlation coefficient between geochemical and seismic temperature estimates is less than the correlation coefficient between La/Sm and Δ*V*_{s} (*R* = 0.65 and *R* = 0.54, respectively). Although we have shown that asthenospheric temperature variations control the correlation between La/Sm and Δ*V*_{s}, lithospheric thickness changes represent an important secondary geochemical signal. Since our geochemical potential temperatures are calculated by accounting for lithospheric thickness changes, the decrease in correlation coefficient that we observe may result from removal of this secondary signal. Despite the inherent uncertainties associated with both modeling approaches, the positive correlation between these independent estimates affirms the underpinning assumptions and highlights the potential rewards of combining temperatures estimated from geochemical analyses of intraplate volcanic rocks with temperatures calculated by calibrating tomographic models.

In order to highlight the consistency between geochemical and tomographic estimates of mantle potential temperature and lithospheric thickness, we present a hemispherical transect that summarizes regional variations of La/Sm, Δ*V*_{s} at a depth of 150 ± 25 km, Δ*T*_{p} at a depth of 150 ± 25 km, and *a* across Europe and Africa (Fig. 5b). This transect starts at the Icelandic plume, crosses Central Europe and Arabia, traverses the East African Rift, and finishes on the South African craton (Fig. 5c). The lowest global values of La/Sm and Δ*V*_{s} are associated with the Icelandic plume (Fig. 5d). Further south, there is a gradual decline in both of these values between the Eifel region and Central Arabia. South of the Afar region, La/Sm and Δ*V*_{s} values increase again, although there is some scatter. As expected, there are concomitant increases and decreases in geochemical and tomographic estimates of Δ*T*_{p} (Fig. 5e). Geochemical and tomographic estimates of lithospheric thickness also broadly match along this transect (Fig. 5f). Although it is not guaranteed that a seismically defined lithosphere-asthenosphere boundary should coincide with the depth at which melting ceases, it is likely that the difference between these depths is minimal compared to the uncertainties inherent in both techniques. Calculated lithospheric thicknesses of ~50–70 km occur beneath European and Arabian intraplate volcanic provinces. Given that crustal thicknesses in these regions are ~30–40 km, our combined geochemical and tomographic values suggest that the lithospheric mantle is remarkably thin beneath these provinces^{49}. Note that in regions where the lithosphere is >125 km thick, a correlation between tomographically derived estimates of lithospheric thickness and the values of Δ*T*_{p} at a depth of 150 ± 25 km is expected since these values of *T*_{p} reside within the lithospheric mantle (e.g., North Sea, Tanzania).

## Discussion

Our global analysis of volcanic provinces and shear-wave velocity anomalies suggests that a combination of asthenospheric temperature anomalies and lithospheric thickness variations helps to determine the spatial distribution of intraplate magmatism. This inference is manifest by a range of geochemical and geophysical observations. First, the bulk of Neogene-Quaternary intraplate volcanism coincides with slow shear-wave velocity anomalies and thin lithosphere. Secondly, there is a positive correlation between La/Sm values of intraplate volcanic rocks and the average magnitude of these shear-wave velocity anomalies within an asthenospheric channel centered on 150 ± 25 km. Thirdly, potential temperatures and lithospheric thicknesses calculated by geochemical inverse modeling of REE concentrations generally match those values obtained from a *V*_{s}-*T* calibration of the SL2013sv tomographic model. Significantly, our results complement those obtained by a global study of the mid-oceanic ridge system, which shows that shear-wave velocity anomalies correlate with major element compositions of basaltic rocks and with axial ridge depths^{17}. Note that, in addition to the influence of temperature anomalies, these correlations could be partly modified by mantle compositional variation or by crystal fractionation beneath mid-oceanic ridges^{50,51}.

Intraplate volcanism is spatially associated with thin lithosphere and its melt composition can be used to estimate asthenospheric temperature. Both thinning of undepleted lithospheric mantle and elevated asthenospheric temperature are significant means for generating regional uplift on horizontal length scales of order 10^{3} km and on timescales of order 10 Ma. Therefore, our results have direct implications for the spatial and temporal evolution of dynamic topography, which is generated and maintained by mantle convective processes^{52,53}. We define dynamic topography to embrace long-wavelength topography generated by mantle flow, isostatic responses to thermochemical processes within the convecting mantle, as well as regional changes in thickness of the lithospheric mantle^{54,55}. If we assume an equilibrated lithospheric thickness of *a*_{0} = 120–150 km at mean sea level, the amount of regional uplift is given by

where *a*_{1} is the present-day lithospheric thickness, *α* = 3 × 10^{−5}^{ ∘}C^{−1} is thermal expansivity, Δ*T* is an asthenospheric temperature anomaly of thickness *h*, and *T*_{1} = 1390 °C is ambient asthenospheric temperature (i.e., the temperature at 150 km depth provided *T*_{p} = 1330 °C, mantle rock density = 3300 kg m^{−3}, and specific heat capacity = 1187 J kg^{−1 }K^{−1} (refs. ^{42,56,57}). If lithospheric thickness is reduced to, say, *a*_{1} = 60 km, if Δ*T* = 50 ^{ ∘}C, and if *h* = 100 km, *U* = 0.81–1.34 km for the disequilibrated case where the effects of pressure upon density are ignored. Variations in Δ*T* will change *U* by ~3 m °C^{−1}. Thus, a significant corollary of our proposed mechanism for generation of intraplate volcanism is rapid large-scale uplift. This prediction is easily tested by examining the relationship between topographic elevation and emergent (but undeformed) marine sedimentary rocks that crop out in regions where intraplate volcanism predominates, where asthenospheric temperatures are anomalously high, and where continental lithosphere is anomalously thin. Dramatic examples include western North America, southernmost South America, eastern Australia, and North Africa where marine deposits of negligible dip indicate that hundreds of meters of wholesale rapid uplift occurred during Cenozoic times (Fig. 1b^{55}). In western Arabia and in eastern Anatolia, Cenozoic marine rocks occur at Jabal Umm Himar and within the Sivas Basin at present-day elevations of 1.2 and 1.5 km, respectively (Fig. 5c, f^{[ 58,59}). Regional uplift of these undeformed marine deposits can be achieved by a combination of emplacement of an asthenospheric temperature anomaly and lithospheric mantle thinning^{60,61}.

Dynamic topography undoubtedly makes a profound contribution to relative sea-level changes, ancient oceanic circulation, fluvial drainage patterns, and sedimentary deposition^{62,63,64,65,66,67}. Here, we have shown that intraplate magmatism, positive asthenospheric temperature anomalies and thinned lithosphere are closely related during Neogene-Quaternary times. Removal or thinning of lithospheric mantle produces regional uplift but subsidence will eventually occur as the lithosphere conductively cools and re-thickens^{68,69}. These processes generate vertical motions of the order of hundreds of meters at the Earth’s surface that are superimposed on motions generated by mantle convection. There is considerable debate about the way in which mantle flow contributes to dynamic topography, and about the relative contributions of the upper and lower mantle^{70}. By combining the distribution and composition of intraplate volcanism with calibrated tomographic models, we have shown that asthenospheric temperature variation and lithospheric thickness changes can make a significant contribution to dynamic topography. Since intraplate magmatism occurs throughout the stratigraphic record, we suggest that a combined analysis of igneous rocks and stratigraphy could help to elucidate the history of mantle convective and dynamic topographic phenomena throughout the Phanerozoic Era.

## Methods

### Tomographic model

For this study, we primarily exploited the SL2013sv tomographic model^{21}. This global upper mantle model uses body and surface waves that include both fundamental and higher modes with periods of 11–450 s. A total of ~750,000 seismograms were analyzed by ref. ^{21} and misfits were calculated by comparing observed and calculated waveforms (≤18th overtone). Significantly, the SL2013sv model inverts for crustal structure during the optimization process whereas other global models mostly, but not always, use a defined a priori crustal architecture. The SL2013sv model has a vertical resolution of 25–50 km and chequerboard tests demonstrate that features with diameters of ~600 km can be clearly resolved at lithospheric depths. In areas of greater ray-path coverage (e.g., North America, Europe) finer scale features should be resolvable. Here, this model is shown relative to AK-135 reference model recomputed for a reference period of 50s and slightly modified following a preliminary optimization^{71}. This reference model also incorporates laterally varying crustal structure that is initially based upon CRUST2.0 and is continually updated during optimization^{21,49}. Descriptions of the CAM2016Vsv, SEMUCB-WM1 and S40RTS tomographic models are provided in Supplementary Information^{19,22,34}.

### Shear-wave velocity to temperature conversion

We exploit a *V*_{s}-to-*T* conversion scheme to estimate temperature variations of the upper mantle from shear-wave velocities. The temperature model exploited in our study is created by applying the empirical approach of ref. ^{27} to the SL2013sv tomographic model^{21}. This scheme was originally developed by ref. ^{26} and it has been updated to exploit an empirical formula, which describes *V*_{s} as a function of pressure and temperature^{13}. Several of the parameters used in this empirical formula were experimentally constrained^{13}. However, seven material constants are required to be independently calibrated for each different tomographic model. These constants are: \({\mu }_{U}^{0}\), \(\frac{\partial {\mu }_{U}}{\partial T}\), \(\frac{\partial {\mu }_{U}}{\partial P}\), *ν*_{r}, *E*_{a}, *V*_{a} and \(\frac{\partial {T}_{s}}{\partial z}\), which represent the unrelaxed shear modulus at surface pressures and temperatures, how the shear modulus changes as a function of temperature, how the shear modulus changes as a function of pressure, shear viscosity for a reference pressure, temperature and grain size (1200 °C, 1.5 GPa, 1 mm), activation energy, activation volume and solidus temperature as a function of depth, respectively.

All seven empirically determined parameters are permitted to vary within physically plausible limits (e.g., refs. ^{13,26,72}). The resultant *V*_{s}(*T*, *P*) parameterization is calibrated by minimizing four misfit functions between observed and calculated values in respect of four different sets of constraints (i.e., *H*_{1}–*H*_{4}). For oceanic lithosphere, *V*_{s} varies as a function of depth and plate age^{12}. Oceanic plate profiles parallel to plate spreading (i.e., flowlines) taken from the SL2013sv tomographic model are averaged and *V*_{s} as a function of depth and plate age is determined. These stacked tomographic profiles are compared to values of *V*_{s} estimated by applying the *V*_{s}-to-*T* parameterization to an oceanic plate cooling model^{29}. *H*_{1} is then calculated by computing the rms misfit between observed and calculated *V*_{s} profiles at depths of 75 and 100 km^{27}. Oceanic plate cooling models can only be used to minimize the difference between observed and calculated values of *V*_{s} at depths ≤125 km^{29}. Therefore, a second misfit function, *H*_{2}, is constructed at depths beneath the upper-thermal boundary layer (i.e., from >225 to 400 km). *V*_{s} values as a function of depth for the SL2013sv tomographic model are averaged across oceanic regions and compared with values of *V*_{s} estimated from a 1333 °C isentrope calculated for peridotite using material constants from ref. ^{42}. Both the oceanic plate cooling model and the temperature profile for convecting mantle are assumed to have a *T*_{p} of 1333 °C^{27,29,42}. Mantle with a *T*_{p} value of 1333 °C is considered to be at ambient temperatures within the final model. These empirically derived parameters are used to provide a prediction of shear-wave attenuation, *Q*^{−1} (ref. ^{73}). The calculated value of *Q*^{−1} can be compared to an estimate of seismic attenuation at depths >150 km beneath old oceanic floor (>100 Ma^{74}). The misfit between observed and predicted *Q*^{−1}, *H*_{3}, is calculated for these locations. Finally, the misfit between an assumed value of shear viscosity, \({\mathrm{log}\,}_{10}[{\nu }_{ref}]=3\times 1{0}^{20}\) Pa s, and the mean value of \({\mathrm{log}\,}_{10}[\nu ]\) determined beneath the thermal boundary layer, *H*_{4}, is calculated. The combined misfit, *H*, between observed and calculated values of *V*_{s}, is quantified by summing weighted versions of these four rms misfit functions, *H*_{1}–*H*_{4}, where

*w*_{1}–*w*_{4} are weighting coefficients with values of 10, 1, 2 and 2, respectively^{27}. See refs. ^{27} and^{47} for a detailed explanation of this *V*_{s}-to-*T* parameterization.

### Kolmogorov–Smirnov statistical test

A two-sample Kolmogorov–Smirnov test is used to calculate how likely it is that two cumulative distribution functions, CDFs, are drawn from the same reference distribution^{30}. The probability, *P*, is approximated by

where *D* is the Kolmogorov–Smirnov Test Statistic (i.e., the maximum magnitude difference between two CDFs). *D* can vary between 0 and 1. *m* and *n* are the number of samples in each CDF, which in this case are areal distribution of data (i.e., shear-wave velocity or lithospheric thickness) over the Earth’s surface and the subset of that distribution, which contains recent intraplate volcanism, respectively.

Database 1 is used to define locations of recent intraplate volcanism. Here, Database 1 is filtered to remove samples that are >10 Ma and samples that are >400 km from their predicted site of eruption. Note that we do not filter for MgO content. For ease of analysis, we have subdivided the Earth’s surface into 1° × 1° bins where *p* is the total number of bins. The subset of bins containing intraplate volcanic rocks is termed *q*. To avoid biasing this analysis by oversampling at high latitudes, the number of bins in each distribution, *p* and *q*, are weighted by the latitude of each bin (i.e., a bin closer to the equator covers a larger geographic area). *m* and *n* are calculated by weighting each bin, *p*_{i} and *q*_{i}, by \(\cos {\phi }_{i}\) where *ϕ*_{i} is the latitude of the bin such that

where \(\frac{\pi }{2}\) is the global average bin weighting. Thus, one bin at the equator in the distribution *p* will count for ≈1.6 bins while a bin at 80° north will count as ≈0.3 bins. The probability that 1° bins containing intraplate volcanic samples are randomly distributed so that the great majority lie above regions of negative Δ*V*_{s} is <1 in 10^{100}. The probability that bins containing intraplate volcanic samples are randomly distributed such that they lie above thin (<100 km thick) lithosphere is <1 in 10^{150}.

### Correlation coefficient

The quality of correlation between two datasets, for example between La/Sm and Δ*V*_{s}, can be determined using the Pearson product-moment correlation coefficient, *R*. In this study, we compare global databases subdivided into 1° × 1° bins. Since the geographic area of a given bin depends upon its latitude, the influence of each bin, *i*, on the value of *R* is once again assigned a weighting, \({w}_{i}=\cos {\phi }_{i}\), so that bins closer to the poles that cover smaller geographical areas have less effect on the value of *R*. The weighted Pearson product-moment correlation coefficient is defined as

where

The symbols *x*_{i} and *y*_{i} denote the *x* and *y* values of each bin. The intercept, *β*_{0}, and the slope, *β*_{1}, which define the linear best-fit to observed values are calculated using

To determine the statistical significance of these correlations, we can use either the population correlation coefficient or a look-up table of *t*-test critical values, which define the lowest value of *R* that is distinguishable from 0 for a given sample size. For these significance-limit calculations, we apply a weighting of 1 to all bins.

### Principal component analysis

Before principal component analysis is carried out, Database 1 is filtered so that only samples with 9 < MgO wt% <14.5, which are <10 Ma, which are <400 km from their predicted site of eruption, and which contain all eight incompatible trace elements, are included. These measurements are binned in the way described in the main text and any bins that have ≤5 samples are removed. Each bin is then normalized to average composition for each element within the filtered and binned version of Database 1. Principal component analysis is carried out using the statsmodels.multivariate.pca routine from the Python software package where the following settings are selected: number of calculated components = 3; standardize = “True”; and method = “NIPALS”.

To construct a synthetic dataset, elemental concentrations are calculated using the INVMEL model for a range of *T*_{p} values, lithospheric thicknesses, and mantle compositions. Random distributions of *T*_{p}, lithospheric thickness and *ε*Nd are generated using Gaussian distributions, each of which is defined by a mean and standard deviation of 1350 °C and 40 °C, 55 km and 8 km, and 5 and 2, respectively. Outer limits for *T*_{p}, lithospheric thickness and *ε*Nd are set at 1250–1450 °C, 35–75 km, and 0–10, respectively. Finally, resultant *T*_{p}, lithospheric thickness and *ε*Nd values are rounded to the nearest 5 °C, 1 km and 0.25, respectively. In each case, 150 values are generated and these values are combined to describe a series of input parameters for 150 INVMEL models. Concentrations of Ba, Nb, K, La, Nd, Zr, Sm and Yb calculated for these model runs are then modeled using principal component analysis. This procedure is repeated 99 times.

### Geochemical estimates of *T*
_{p}

*T*

The INVMEL-v12 geochemical model is used to calculate rare earth element (REE) concentrations generated by mantle melting for different values of *T*_{p} and lithospheric thickness^{41}. Concentration of each REE within an instantaneous melt, *c*_{l}, and concentration within the residue, *c*_{s}, are related to each other by two equations, which must be simultaneously solved. These equations are

where *X* is the total melt fraction that has been removed, \(\bar{D}\) is the bulk distribution coefficient for any given element within the solid assemblage, and \(\bar{P}\) is the bulk distribution coefficient for the melting assemblage^{75}. \(\bar{P}\) and \(\bar{D}\) vary as a function of depth since both the mineral assemblage present and individual mineral distribution coefficients are depth-dependent. To calculate *c*_{l} and *c*_{s} at each depth, the expressions in Eq. (8) are numerically integrated using a fourth-order Runga-Kutta scheme from the base of the melting interval, where *X* = 0, to the top of the melting interval, where \(X=X^{\prime}\). The average melt composition for all melt extracted from a rock at a single depth, *C*, is related to *c*_{l} by

The mean composition of all of the melt extracted from the melting interval 0–*h* is given by

\(X^{\prime} (z)\) is calculated from the equations of ref. ^{35} updated using revised parameter values from ref. ^{42}. \(X^{\prime}\) is a function of potential temperature, *T*_{p}, lithospheric thickness, *a*, the weight fraction of water present within the source region, \({X}_{{{\rm{H}}}_{2}{\rm{O}}}^{{\mathrm{bulk}}}\), and the modal clinopyroxene by mass of the peridotite source, *M*_{cpx}. \({X}_{{{\rm{H}}}_{2}{\rm{O}}}^{{\mathrm{bulk}}}\) is approximated from the concentration of Ce within the source region and \({X}_{{\rm{Ce}}}^{{\mathrm{bulk}}}\) is constrained by assuming that \({X}_{{{\rm{H}}}_{2}{\rm{O}}}^{{\mathrm{bulk}}}/{X}_{{\rm{Ce}}}^{{\mathrm{bulk}}}=200\)^{76}. *M*_{cpx} = 0.17 is fixed for all models.

Composition of the mantle source region used in this modeling procedure varies as a function of the average value of *ε*Nd calculated for each province. It is calculated by linearly mixing between a primitive mantle source, PM, and a depleted MORB mantle, DMM, in order to match the observed value of *ε*Nd (Supplementary Table 5^{77}). Note that varying REE element concentrations within the source region as a function of *ε*Nd also affects the value of \({X}_{{{\rm{H}}}_{2}{\rm{O}}}^{{\mathrm{bulk}}}\) for the source since we assume \({X}_{{{\rm{H}}}_{2}{\rm{O}}}^{{\mathrm{bulk}}}/{X}_{{\rm{Ce}}}^{{\mathrm{bulk}}}=200\). Exploiting the melting parameterization of ref. ^{42} and assuming an *ε*Nd value of 10, a value of *T*_{p} = 1312 ± 28 °C is required to generate 6.9 ± 2.2 km of oceanic crust at a mid-oceanic ridge, assuming a triangular melt geometry^{29}. We thus assume that 1312 °C is the ambient value of *T*_{p}. Within the mantle, the aluminous phase present varies as a function of depth between plagioclase, spinel and garnet. The mineral proportions used within each stability field are given in Supplementary Table 6. The plagioclase-spinel transition zone is set at 25–35 km and the spinel-garnet transition zone is set at 63–72 km^{2,37}.

The bulk distribution coefficient of the melting assemblage, \(\bar{P}\), is calculated to determine the concentration of REEs within the melt phase. The INVMEL model incorporates two melting regimes: one regime where all minerals listed in Supplementary Table 6 are present; and a second regime where the more fusible minerals are exhausted, leaving only olivine and orthopyroxene. The proportion of each mineral present within the source region by weight is given by *F*_{n}, where *n* = 1, 2, 3, 4, 5 and 6 refer to olivine, orthopyroxene, clinopyroxene, plagioclase, spinel and garnet, respectively. The point at which clinopyroxene and aluminous phases (i.e., plagioclase, spinel and garnet) become exhausted is referred to as *X*_{1}. *X*_{1} is assumed to be located where *X* = 0.18. The proportion of clinopyroxene and aluminous phase present is assumed to vary linearly between the initial proportion of each mineral present within the source, \({F}_{n}={F}_{n}^{0}\), and when the minerals are exhausted *F*_{n} = 0, between *X* = 0 and *X* = *X*_{1}. The proportion of olivine to orthopyroxene remains fixed and the concentration of these minerals within the source region is varied so that

where *N* is the total number of mineral phases.

To calculate the bulk distribution coefficient for a given element within the solid assemblage, \(\bar{D}\), the partition coefficients for each mineral, *D*_{i}, must be parameterized. Values of *D*_{i} for each REE in olivine, orthopyroxene and spinel are listed in Supplementary Table 8. In the mantle, partition coefficients necessarily vary as a function of pressure, temperature and mineral chemistry as sites within mineral lattices expand and contract. The equation for calculating the partitioning of an element with a charge *v*+ and a radius *r*_{i} entering into site *M* of a crystalline lattice is governed by

where *N*_{A} is Avogadro’s number, *R* is the gas constant, \({E}_{M}^{v+}\) is the elastic response of the site *M*, \({r}_{0(M)}^{v+}\) is the radius of the site and \({D}_{0(M)}^{v+}\) is the partition coefficient for an element with a charge *v*+ and a radius \({r}_{i(M)}^{v+}\)^{78}. The constants required to calculate *D*_{i} for REEs in clinopyroxene, plagioclase and garnet as a function of pressure, temperature and chemistry are listed in Supplementary Tables 8 and 9. Partitioning of REEs into garnet is governed by the fraction of pyrope present, PYR, and can be calculated from

where MgO_{mol}, Cr_{2}O_{4mol}, Al_{2}O_{3mol}, FeO_{mol}, MgO_{mol}, CaO_{mol} are the fraction by weight of each oxide within garnet divided by their respective molecular weights (Supplementary Table 7^{79}). Partitioning of REEs into clinopyroxene depends upon the proportions of Ca and Al within the clinopyroxene phase, \({x}_{Ca}^{M2}\) and \({x}_{Al}^{M1}\), respectively^{78}. These proportions are calculated using

where *x*_{i}, *i*_{mol}, *C*_{i} and *O*_{i} are the proportions of oxide *i* within clinopyroxene, fraction by weight of oxide *i* within clinopyroxene divided by its molecular weight, the number of cations in the oxide *i*, the number of oxygen atoms in the oxide *i* listed in Supplementary Table 7, respectively.

### Grid search procedure

Before comparing calculated REE concentrations, \({{\mathscr{C}}}_{m}\), with observed REE concentrations, \({{\mathscr{C}}}_{o}\), the observed concentrations are corrected for olivine fractional crystallization. For each sample, the proportion of olivine that must be added into the melt, *χ*_{ol}, to ensure that the melt is in equilibrium with olivine of 90% forsterite content, is calculated. The equations of ref. ^{4} are used, assuming Fe^{3+}/∑Fe = 0.15^{43}. The average *χ*_{ol} value for all samples within a given volcanic province is used and no REEs are assumed to be present within the cumulate olivine. Corrected observed REE concentrations, \({{\mathscr{C}}}_{c}\), are therefore given by

To identify that model that best fits the corrected REE concentrations, \({{\mathscr{C}}}_{c}\) is compared with \({{\mathscr{C}}}_{m}\) calculated for each pair of *T*_{p} and *a* values at 1 °C and 1 km intervals, respectively. *T*_{p} is varied between 1250–1550 °C and *a* is varied between 30 and 100 km. Certain samples do not have recorded concentrations for all 14 REEs. To ensure average REE concentrations are consistent, if <50% of samples record concentrations for a given element, it is omitted from the misfit calculation. Following this screening process, if the number of REEs that a province contains, *N*, is <7 then *T*_{p} and *a* are not recorded for that particular province. The root mean squared (rms) misfit, *H*, between \({{\mathscr{C}}}_{c}\) and \({{\mathscr{C}}}_{m}\) at each *T*_{p} and *a* pair is given by

where *σ*_{c} is the standard deviation of \({{\mathscr{C}}}_{c}\) for each province. The best-fitting pair of *T*_{p} and *a* values is taken to be the global minimum for *H*. If the value of *H* at the global minimum is >1 then *T*_{p} and *a* are not recorded for that province. Acceptable models are those where the value of *H* is <1.5 × the global minimum. The errors for each *T*_{p} and *a* estimate show this range of acceptable values.

## Data availability

All geochemical data analyzed during this study are included in this published article and are included and appropriately referenced within Supplementary Data File 1 and Supplementary Tables 1 and 2. Tomographic models, *V*_{s}-to-*T* models, and lithospheric thickness maps analyzed in this study can be accessed through the references provided (e.g., refs. ^{19,20,21,22,27,80}).

## Code availability

The INVMEL software package is available on request from D. McKenzie. All other codes used can be provided by the authors upon request.

## References

- 1.
McKenzie, D. & Bickle, M. J. The volume and composition of melt generated by extension of the lithosphere.

*J. Petrol.***29**, 625–679 (1988). - 2.
McKenzie, D. & O’Nions, R. K. Partial melt distributions from inversion of rare-earth element concentrations.

*J. Petrol.***32**, 1021–1091 (1991). - 3.
Herzberg, C. et al. Temperatures in ambient mantle and plumes: constraints from basalts, picrites and komatiites.

*Geochem. Geophys. Geosyst.***8**, 2006GC001390 (2007). - 4.
Lee, C. T. A., Luffi, P., Plank, T., Dalton, H. & Leeman, W. P. Constraints on the depths and temperatures of basaltic magma generation on Earth and other terrestrial planets using new thermobarometers for mafic magmas.

*Earth Planet. Sci. Lett.***279**, 20–33 (2009). - 5.
Hofmann, A. & Hart, S. An assessment of local and regional isotopic equilibrium in the mantle.

*Earth Planet Sci. Lett.***38**, 44–62 (1978). - 6.
Zindler, A. & Hart, S. Chemical geodynamics.

*Ann. Rev. Earth Planet. Sci.***14**, 493–571 (1986). - 7.
Hawkesworth, C. et al. Paraná magmatism and the opening of the South Atlantic.

*Geol. Soc. Spec. Publ.***68**, 221–240 (1992). - 8.
Kay, R. W. & Kay, S. M. Delamination and delamination magmatism.

*Tectonophysics***219**, 177–189 (1993). - 9.
King, S. D. & Anderson, D. L. Edge-driven convection.

*Earth Planet. Sci. Lett.***160**, 289–296 (1998). - 10.
Conrad, C. P., Wu, B., Smith, E. I., Bianco, T. A. & Tibbetts, A. Shear-driven upwelling induced by lateral viscosity variations and asthenospheric shear: a mechanism for intraplate volcanism.

*Phys. Earth Planet Interiors***178**, 162–175 (2010). - 11.
Niu, Y., Wilson, M., Humphreys, E. R. & O’Hara, M. J. The origin of intra-plate ocean island basalts (OIB): the lid effect and its geodynamic implications.

*J. Petrol.***52**, 1443–1468 (2011). - 12.
Priestley, K. & McKenzie, D. The thermal structure of the lithosphere from shear wave velocities.

*Earth Planet. Sci. Lett.***244**, 285–301 (2006). - 13.
Yamauchi, H. & Takei, Y. Polycrystal anelasticity at near-solidus temperatures.

*J. Geophys. Res. Solid Earth***121**, 7790–7820 (2016). - 14.
Schutt, D. & Lesher, C. Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite.

*J. Geophys. Res.***111**, 1–24 (2006). - 15.
Hirschmann, M. M. Partial melt in the oceanic low velocity zone.

*Phys. Earth Planet. Interiors***179**, 60–71 (2010). - 16.
Klein, E. M. & Langmuir, C. H. Global correlations of ocean ridge basalt chemistry with axial depth and crustal thickness.

*J. Geophys. Res.***92**, 8089–8115 (1987). - 17.
Dalton, C. A., Langmuir, C. H. & Gale, A. Geophysical and geochemical evidence for deep temperature variations beneath mid-ocean ridges.

*Science***344**, 80–83 (2014). - 18.
Argus, D. F., Gordon, R. G. & DeMets, C. Geologically current motion of 56 plates relative to the no-net-rotation reference frame.

*Geochem. Geophys. Geosyst.***12**, 1–13 (2011). - 19.
Ritsema, J., Deuss, A., Van Heijst, H. & Woodhouse, J. 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). - 20.
French, S., Lekic, V. & Romanowicz, B. Waveform tomography reveals channeled flow at the base of the oceanic asthenosphere.

*Science***342**, 227–30 (2013). - 21.
Schaeffer, A. J. & Lebedev, S. Global shear speed structure of the upper mantle and transition zone.

*Geophys. J. Int.***194**, 417–449 (2013). - 22.
Ho, T., Priestley, K. & Debayle, E. A global horizontal shear velocity model of the upper mantle from multimode Love wave measurements.

*Geophys. J. Int.***207**, 542–561 (2016). - 23.
Klöcking, M., White, N. J., Maclennan, J., McKenzie, D. & Fitton, J. G. Quantitative relationships between basalt geochemistry, shear wave velocity, and asthenospheric temperature beneath western North America.

*Geochem. Geophys. Geosyst.***19**, 3376–3404 (2018). - 24.
Gibson, S. A. & Geist, D. Geochemical and geophysical estimates of lithospheric thickness variation beneath Galápagos.

*Earth Planet. Sci. Lett.***300**, 275–286 (2010). - 25.
Plank, T. & Forsyth, D. Thermal structure and melting conditions in the mantle beneath the Basin and Range province from seismology and petrology.

*Geochem. Geophys. Geosyst.***17**, 1312–1338 (2016). - 26.
Priestley, K. & McKenzie, D. The relationship between shear wave velocity, temperature, attenuation and viscosity in the shallow part of the mantle.

*Earth Planet. Sci. Lett.***381**, 78–91 (2013). - 27.
Hoggard, M. J. et al. Global distribution of sediment-hosted metals controlled by craton edge stability.

*Nat. Geosci.***13**, 504–510 (2020). - 28.
Burgos, G. et al. Oceanic lithosphere-asthenosphere boundary from surface wave dispersion data.

*J. Geophys. Res.***119**, 1079–1093 (2014). - 29.
Richards, F. D., Hoggard, M. J., Cowton, L. & White, N. J. Reassessing the thermal structure of oceanic lithosphere with revised global inventories of basement depths and heat flow measurements.

*J. Geophys. Res.***123**, 9136–9161 (2018). - 30.
Kolmogorov, A. Sulla determinazione empirica di una lgge di distribuzione.

*Inst. Ital. Attuari, Giorn.***4**, 83–91 (1933). - 31.
Schilling, J.-G. In

*Mantles of the Earth and Terrestrial Planets*(ed. Runcorn, S. K.) 267–283 (Interscience Publication, 1967). - 32.
Kay, R. W. & Gast, P. W. The rare earth content and origin of alkali-rich basalts.

*J. Geol.***81**, 653–682 (1973). - 33.
Schilling, J.-G. Azores mantle blob: rare-earth evidence.

*Earth Planet. Sci. Lett.***25**, 103–115 (1975). - 34.
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). - 35.
Katz, R. F., Spiegelmann, M. & Langmuir, C. H. A new parameterization of hydrous mantle melting.

*Geochem. Geophys. Geosyst.***4**, 2002GC000433 (2003). - 36.
Hofmann, A. W. Mantle geochemistry: the message from oceanic volcanism.

*Nature***385**, 219–229 (1997). - 37.
Jennings, E. S. & Holland, T. J. B. A simple thermodynamic model for melting of peridotite in the system NCFMASOCr.

*J. Petrol.***56**, 869–892 (2015). - 38.
Niu, Y. & Green, D. H. The petrological control on the lithosphere-asthenosphere boundary (LAB) beneath ocean basins.

*Earth Sci. Rev.***185**, 301–307 (2018). - 39.
Huang, X.-L., Niu, Y., Xu, Y.-G., Chen, L.-L. & Yang, Q.-J. Mineralogical and geochemical constraints on the petrogenesis of post-collisional potassic and ultrapotassic rocks from western Yunnan, SW China.

*J. Petrol.***51**, 1617–1654 (2010). - 40.
Pilet, S., Baker, M. B., Müntener, O. & Stolper, E. M. Monte Carlo simulations of metasomatic enrichment in the lithosphere and implications for the source of alkaline basalts.

*J. Petrol.***52**, 1415–1442 (2011). - 41.
White, R. S., McKenzie, D. & O’Nions, R. K. Oceanic crustal thickness from seismic measurements and rare earth element inversions.

*J. Geophys. Res.***97**, 19683–19715 (1992). - 42.
Shorttle, O., Maclennan, J. & Lambart, S. Quantifying lithological variability in the mantle.

*Earth Planet. Sci. Lett.***395**, 24–40 (2014). - 43.
Ball, P. W. et al. Quantifying asthenospheric and lithospheric controls on mafic magmatism across North Africa.

*Geochem. Geophys. Geosyst.***20**, 3520–3555 (2019). - 44.
Lee, C.-T. A. Compositional variation of density and seismic velocities in natural peridotites at stp conditions: Implications for seismic imaging of compositional heterogeneities in the upper mantle.

*J. Geophys. Res.***108**, 1–20 (2003). - 45.
McKenzie, D. Constraints on melt generation and transport from U-series activity ratios.

*Chem. Geol.***162**, 81–94 (2000). - 46.
Takei, Y. Effects of partial melting on seismic velocity and attenuation: a new insight from experiments.

*Ann. Rev. Earth Planet. Sci.***45**, 447–470 (2017). - 47.
Richards, F. D., Hoggard, M. J., White, N. & Ghelichkhan, S. Quantifying the relationship between short-wavelength dynamic topography and thermomechanical structure of the upper mantle using calibrated parameterization of anelasticity.

*J. Geophys. Res.***125**, 2019JB019062 (2020). - 48.
Matthews, S., Wong, K., Shorttle, O., Edmonds, M. & Maclennan, J. Do olivine crystallization temperatures faithfully record mantle temperature variability?

*Geochem. Geophys. Geosyst.*https://doi.org/10.1029/2020GC009157 (2021). - 49.
Bassin, C., Laske, G. & Masters, G. The current limits of resolution for surface wave tomography in North America.

*EOS Trans AGU***81**, F897 (2000). - 50.
O’Neill, H. S. C. & Jenner, F. E. Causes of the compositional variability among ocean floor basalts.

*J. Petrol.***57**, 2163–2194 (2016). - 51.
Niu, Y. The meaning of global ocean ridge basalt major element compositions.

*J. Petrol.***57**, 2081–2103 (2016). - 52.
Pekeris, C. L. Thermal convection in the interior of the Earth.

*Geophys. J. Int.***3**, 343–367 (1935). - 53.
Hager, B. & Richards, M. Long-wavelength variations in Earth’s geoid: physical models and dynamical implications.

*Philos. Trans. R. Soc. Lond. Ser A, Math. Phys. Sci.***328**, 309–327 (1989). - 54.
Moucha, R. & Forte, A. M. Changes in African topography driven by mantle convection.

*Nat. Geosci.***4**, 707–712 (2011). - 55.
Hoggard, M. J., Austermann, J., Randel, C. & Stephenson, S. Observational estimates of dynamic topography through space and time, in

*Mantle convection and surface expressions*(eds Marquardt, H., Ballmer, M., Cottaar, S. & Konter J.)*AGU Geophysical Monograph Series*. (Washington, DC, 2021). https://doi.org/10.1002/9781119528609.ch15. - 56.
Parsons, B. & Sclater, J. G. An analysis of the variation of ocean floor bathymetry and heat flow with age.

*J. Geophys. Res.***82**, 803–827 (1977). - 57.
McKenzie, D., Jackson, J. & Priestley, K. Thermal structure of oceanic and continental lithosphere.

*Earth Planet. Sci. Lett.***233**, 337–349 (2005). - 58.
Madden, G. T.

*Paleocene Pycnodont Fishes from Jabal Umm Himar, Harrat Hadan area, Kingdom of Saudi Arabia*(US Geological Survey, 1983). - 59.
Çiner, A., Kosun, E. & Deynoux, M. Fluvial, evaporitic and shallow-marine facies architecture, depositional evolution and cyclicity in the Sivas Basin (Lower to Middle Miocene), Central Turkey.

*J. Asian Earth Sci.***21**, 147–165 (2002). - 60.
Wilson, J. W., Roberts, G. G., Hoggard, M. J. & White, N. J. Cenozoic epeirogeny of the Arabian Peninsula from drainage modeling.

*Geochem. Geophys. Geosyst.***15**, 3723–3761 (2014). - 61.
McNab, F., Ball, P. W., Hoggard, M. J. & White, N. J. Neogene uplift and magmatism of Anatolia: Insights from drainage analysis and basaltic geochemistry.

*Geochem. Geophys. Geosyst.***19**, 175–213 (2018). - 62.
Poore, H., White, N. & Maclennan, J. Ocean circulation and mantle melting controlled by radial flow of hot pulses in the Iceland plume.

*Nat. Geosci.***4**, 558–561 (2011). - 63.
Braun, J., Robert, X. & Simon-Labric, T. Eroding dynamic topography.

*Geophys. Res. Lett.***40**, 1494–1499 (2013). - 64.
Rovere, A., Stocchi, P. & Vacchi, M. Eustatic and relative sea level changes.

*Curr. Clim. Chang. Rep.***2**, 221–231 (2016). - 65.
Austermann, J., Mitrovica, J. X., Huybers, P. & Rovere, A. Detection of a dynamic topography signal in last interglacial sea-level records.

*Sci. Adv.***3**, e1700457 (2017). - 66.
Stephenson, S. N., White, N. J., Li, T. & Robinson, L. F. Disentangling interglacial sea level and global dynamic topography: Analysis of Madagascar.

*Earth Planet. Sci. Lett.***519**, 61–69 (2019). - 67.
Ding, X., Salles, T., Flament, N., Mallard, C. & Rey, P. F. Drainage and sedimentary responses to dynamic topography.

*Geophys. Res. Lett.***46**, 14385–14394 (2019). - 68.
Bird, P. Continental delamination and the Colorado Plateau.

*J. Geophys. Res.***84**, 7561–7571 (1979). - 69.
McKenzie, D. Some remarks on the development of sedimentary basins.

*Earth Planet. Sci. Lett.***40**, 25–32 (1978). - 70.
Davies, D. et al. Earth’s multi-scale topographic response to global mantle flow.

*Nat. Geosci.***12**, 845–856 (2019). - 71.
Lebedev, S. & Van Der Hilst, R. D. Global upper-mantle tomography with the automated multimode inversion of surface and S-wave forms.

*Geophys. J. Int.***173**, 505–518 (2008). - 72.
Jain, C., Korenaga, J. & Karato, S.-i. On the grain size sensitivity of olivine rheology.

*J. Geophys. Res.***123**, 674–688 (2018). - 73.
McCarthy, C., Takei, Y. & Hiraga, T. Experimental study of attenuation and dispersion over a broad frequency range: 2. The universal scaling of polycrystalline materials.

*J. Geophys. Res.***116**, 1–18 (2011). - 74.
Dalton, C. A., Ekström, G. & Dziewonski, A. M. Global seismological shear velocity and attenuation: a comparison with experimental observations.

*Earth Planet. Sci. Lett.***284**, 65–75 (2009). - 75.
Shaw, D. M. Trace element fractionation during anatexis.

*Geochimi. et Cosmochim. Acta***34**, 237–243 (1970). - 76.
Michael, P. Regionally distinctive sources of depleted MORB: Evidence from trace elements and H

_{2}O.*Earth Planet. Sci. Lett.***131**, 301–320 (1995). - 77.
McKenzie, D. & O’Nions, R. K. The source regions of ocean island basalts.

*J. Petrol.***36**, 133–159 (1995). - 78.
Wood, B. J. & Blundy, J. D. A predictive model for rare earth element partitioning between clinopyroxene and anhydrous silicate melt.

*Contrib. Mineral. Petrol.***129**, 166–181 (1997). - 79.
Van Westrenen, W., Blundy, J. D. & Wood, B. J. Effect of Fe

^{2+}on garnet–melt trace element partitioning: experiments in FCMAS and quantification of crystal-chemical controls in natural systems.*Lithos***53**, 189–201 (2000). - 80.
Priestley, K., McKenzie, D. & Ho, T. In

*Lithospheric Discontinuities*(eds Yuan, H. & Romanowicz, B.) chap. 6, 111–123 (American Geophysical Union, 2019).

## Acknowledgements

P.W.B. acknowledges support by Shell Research. We are grateful to K. Czarnota, I. Frame, L. Kennan, M.J. Hoggard, M. Klöcking, D. Lyness, F. M^{c}Nab, C.P.B. O’Malley, F. Richards, G. Roberts, O. Shorttle, and K. Warners for their help. D. McKenzie generously provided software for the INVMEL forward model. Y. Niu provided an insightful review. Figures were prepared using Generic Mapping Tools software. Database 1 was compiled with the aid of the GEOROC database (www.georoc.edu). Cambridge Earth Sciences contribution number esc.6010.

## Author information

### Affiliations

### Contributions

All authors jointly contributed to the conceptualization and design of this research project. P.W.B. collated the geochemical database, performed the geochemical modeling, carried out the statistical analysis, and drafted the figures. P.W.B. and N.J.W. wrote the paper with assistance and advice from J.M. and S.N.S.

### Corresponding authors

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks Yaoling Niu and the other, anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

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

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Ball, P.W., White, N.J., Maclennan, J. *et al.* Global influence of mantle temperature and plate thickness on intraplate volcanism.
*Nat Commun* **12, **2045 (2021). https://doi.org/10.1038/s41467-021-22323-9

Received:

Accepted:

Published:

## Comments

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