Global influence of mantle temperature and plate thickness on intraplate volcanism

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. Here, the authors compile a global geochemical database of Neogene-Quaternary intraplate volcanism. By comparing the distribution and composition of these rocks with tomographic models they show that intraplate volcanism can be used to constrain upper-mantle structure at the time of eruption.

M antle convection drives plate tectonics, redistributes heat, cycles chemical species, and generates dynamic topography at the Earth's surface. Constraining spatiotemporal 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 shearwave 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 presentday 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 subaerial 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 colocation 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  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 lithosphereasthenosphere 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 midoceanic 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 rethickens 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 shearwave 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: ∂P , ν r , E a , V a and ∂T s ∂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, log 10  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°b ins 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 ϕ i where ϕ i is the latitude of the bin such that where π 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 productmoment 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 ϕ i , so that bins closer to the poles that cover smaller geographical areas have less effect on the value of R. The weighted Pearson productmoment 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 . 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 and X 0 the total melt fraction at the current and previous time step, respectively; D is the bulk distribution coefficient for any given element within the solid assemblage, and P is the bulk distribution coefficient for the melting assemblage 75 .
P and 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 0 . 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 0 ðzÞ is calculated from the equations of ref. 35 updated using revised parameter values from ref. 42 . X 0 is a function of potential temperature, T p , lithospheric thickness, a, the weight fraction of water present within the source region, X bulk H 2 O , and the modal clinopyroxene by mass of the peridotite source, M cpx . X bulk H 2 O is approximated from the concentration of Ce within the source region and X bulk Ce is constrained by assuming that X bulk H 2 O =X bulk Ce ¼ 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 bulk H 2 O for the source since we assume X bulk H 2 O =X bulk Ce ¼ 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, 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 0 n , 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 if n < 3; p n ¼ F n ð1 À ∑ N n¼3 p n Þ if X < X 1 ; D ¼ ∑ N n¼1 F n D n and P ¼ ∑ if X ≥ X 1 ; where N is the total number of mineral phases.
To calculate the bulk distribution coefficient for a given element within the solid assemblage, 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 À4πN A E vþ M 1 2 r vþ 0ðMÞ ðr i À r vþ 0ðMÞ Þ 2 þ 1 3 ðr i À r vþ 0ðMÞ Þ where N A is Avogadro's number, R is the gas constant, E vþ M is the elastic response of the site M, r vþ 0ðMÞ is the radius of the site and D vþ 0ðMÞ is the partition coefficient for an element with a charge v+ and a radius r vþ iðMÞ 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 PYR ¼ MgO mol ð1 À  Table 7 79 ). Partitioning of REEs into clinopyroxene depends upon the proportions of Ca and Al within the clinopyroxene phase, x M2 Ca and x M1 Al , respectively 78 . These proportions are calculated using x M2 Ca ¼ x Ca ; ð19Þ 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, C m , with observed REE concentrations, 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, C c , are therefore given by To identify that model that best fits the corrected REE concentrations, C c is compared with 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 C c and C m at each T p and a pair is given by where σ c is the standard deviation of 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.