# A poorly mixed mantle transition zone and its thermal state inferred from seismic waves

Lauren Waszek, Benoit Tauzin, Nicholas Schmerr, Maxim Ballmer, Juan

Carlos Afonso

► To cite this version:

Lauren Waszek, Benoit Tauzin, Nicholas Schmerr, Maxim Ballmer, Juan Carlos Afonso. A poorly mixed mantle transition zone and its thermal state inferred from seismic waves. Nature Geoscience, 2021, 10.1038/s41561-021-00850-w. hal-03427247

HAL Id: hal-03427247

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

Submitted on 2 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
A poorly mixed mantle transition zone and its thermal state inferred from
seismic waves

L. Waszek¹,²*, B. Tauzin³,⁴, N. C. Schmerr⁵, M. D. Ballmer⁶,⁷, J. C. Afonso⁸

¹James Cook University, Physical Sciences, Douglas, QLD 4811, Australia

²New Mexico State University, Department of Physics, Las Cruces, NM 88003, USA

³Université de Lyon, UCBL, ENSL, CNRS, LGL-TPE, Villeurbanne, France

⁴Australian National University, Research School of Earth Sciences, Acton, ACT 2601,
Australia

⁵University of Maryland College Park, Department of Geology, College Park, MD 20742,
USA

⁶University College London, Department of Earth Sciences, London, WC1E 6BS, UK

⁷ETH Zurich, Institute of Geophysics, 8092 Zurich, Switzerland

⁸Macquarie University, Department of Earth and Planetary Sciences, North Ryde, NSW
2109, Australia

*Correspondence and requests for materials should be addressed to Lauren Waszek at:

lauren.waszek@cantab.net
The abrupt changes in mineralogical properties across the Earth's mantle transition zone impact significantly convection and thermochemical fluxes between the upper and lower mantle. While the 410-km discontinuity at the top of the mantle transition zone is detected with all types of seismic waves, the 660-km boundary is mostly invisible to underside P-wave reflections (P660P). The cause for this observation is debated. The dissociation of ringwoodite and garnet into lower-mantle minerals both contribute to the '660' visibility; only the garnet reaction favours material exchanges across the discontinuity. Here, we combine large datasets of SS and PP precursors, mineralogical modelling, and data-mining techniques to obtain a global thermal map of the mantle transition zone, and explain the lack of P660P visibility. We find that its prevalent absence requires a chemically unequilibrated mantle, and its visibility in few locations is associated with potential temperatures greater than 1800 K. Such temperatures occur in approximately 0.6% of Earth, indicating that the 660 is dominated by the ringwoodite decomposition, which tends to impede mantle flow. We find broad regions with elevated temperatures beneath the Pacific surrounded by major volcanic hotspots, indicating plume retention and ponding of hot materials in the mantle transition zone.
## Introduction

Many intra-plate volcanoes are proposed to result from hot narrow upwellings that originate from the lowermost mantle (1). Conversely, intermittent volcanic activity in the Pacific Ocean, combined with apparent motion of hotspots with respect to plate motion, indicates a shallower origin, perhaps in the mid-mantle or MTZ (2,3). Thus, whether large-scale upwelling encompasses the whole mantle, or occurs within semi-independent layers, remains debated (4,5,6). Solid mineralogical phase changes at 410 and 660 km depths in the MTZ play a key role in regulating mantle circulation. The ‘410’ and ‘660’ correspond primarily to the exothermic transition from olivine (Ol) to wadsleyite (Wd), and the endothermic dissociation of ringwoodite (Rw) into bridgmanite (Br) and ferropericlase (Fp) (7), with further mineralogical reactions adding complexity. These reactions typically result in a thin MTZ in hot regions and a thick MTZ in colder environments (8). The Rw-decomposition at the base of the MTZ acts to impede both cold downwellings and hot upwellings as it changes buoyancy across the phase change (6,8-15). As reactions produce further mineral assemblages around 410 and 660-km depth, some simulations predict lateral deflection and spreading of hot upwelling material (5,16). At high temperatures (greater than 2100K), another phase dissociation from majoritic garnet (Gt) to Al-bearing Mg-rich bridgmanite (Br) dominates at the base of the MTZ, and displays an opposite sign Clapeyron slope to the Rw transition (17). Unlike Rw-decomposition, the post-garnet (Gt) break-down promotes transfer of material from the lower mantle (18,19). Constraining the spatial extent of Gt-decomposition may therefore provide clues on upwelling streams of global mantle circulation.

Currently, the maximum temperature in the MTZ, and hence the global extent of Gt-decomposition dominated regions, is poorly constrained. Comprehensive information about the global temperature distribution in the MTZ can be obtained from the record of seismic
waves sampling the Earth's mantle between globally distributed sources and receivers (20,21). Translating seismic information into absolute mantle temperatures is difficult, as this requires comparing seismic observations with predictions from elastic properties at high pressure and temperature (p, T) and separating the effects of temperature from composition (e.g. 22,23,24). Furthermore, different seismic datasets have their own spatial sensitivity, coverage, and uncertainties that must be incorporated to avoid biases in interpretation (17,24-26).

Here, we combine global seismic datasets, mineral physics modelling, and data-mining techniques, to build a thermal model for the MTZ, constraining regions where the Gt-decomposition dominates over the Rw-decomposition. Progress over previous methodologies (21-23) includes an extended geographic coverage of seismic observations, and a joint analysis of data sensitive to both compressional and shear-wave velocity to alleviate trade-offs between temperature, composition, and chemical equilibration models. The mineral physics analysis is further integrated into an adaptive data-driven scheme for global multi-scale model restoration, incorporating full catalogues of physics-based synthetic data. The resulting thermal model reveals that hotspots in the southern Pacific are not directly located above hot regions of the MTZ, suggesting mediated mass transfer across the mantle

## Seismic observations

Seismic studies provide the first-order detailed constraints on the physical structures of the mantle. Nearly global seismic coverage of MTZ discontinuities is provided by precursors to the shear-wave phase SS and compressional equivalent PP, which travel through the mantle with one free surface reflection (Fig. 1a). Precursors (Fig. 1b) to SS and PP are generated by underside reflections from mantle discontinuities; their relative travel times are sensitive to discontinuity depths, and their amplitudes to the change in elastic properties across the discontinuity. The 410-km discontinuity is visible in both SS and PP data (17,20,21). In
contrast, the 660 is detected globally only in SS data ('S660S'), and typically absent from PP data ('P660P') (20,23).

In order to obtain global seismic constraints on the MTZ, we use manually-picked global datasets of SS (58,217 data) (27) and PP (136,512 data) (Supplementary Information SI1). The small-amplitude precursors are rarely detected in individual seismograms (Fig. 1b), thus we employ common mid-point stacking according to SS/PP bounce points (Fig. 1c,d). The individual data are corrected for 3D crustal and mantle velocity structure (SI2) using CRUST2.0 (28) and S40RTS (13), with PREM (29) as a reference model. Global stacks reveal the 410 signal in both datasets but a missing P660P (Fig. S1), in agreement with previous works (20,23,30). Data are then partitioned spatially. Fixed partitioning schemes introduce geographic biases that potentially mask small-scale structure. We therefore use a stacking approach based upon global tessellations of the Earth using Voronoi cells, which adjust their size to account for sharp changes in topography, noise, and data coverage (Methods A; Fig. S2-S8).

This data mining approach confirms that the 410 is a global feature in both SS and PP datasets, as is the S660S. We detect significant deviations from global average depths on regional length-scales (Fig. 2, S3-4, Table S1, SI3). An apparent P660P signal is detected for only around 5% of the globe (Fig. 2a). However, synthetic tests reveal that many “detections” are likely artefacts (Methods B; Fig. S9, S10). Almost one third of the potential P660P patches are classified as artefacts of interfering phases, and half of the patches produce poor quality signals (Fig. 2a, S11, S12). A robust P660P signal is detected in only four patches (0.6% of the globe): three in the Pacific Ocean, and one at North Pole. Depth profiles along a cross-section through the Pacific Ocean (Fig. 2c-d) reveal narrow regions with single peak P660P signals. Elsewhere, the 660 in PP data is broad and/or extremely low amplitude. P660P signals may exist more widely, but would be below the resolution limit of our data, or
too weak to be detected within the noise. We estimate a detection threshold for P660P from
seismic uncertainties on amplitude and noise levels, depending on the presence or absence of
a P660P signal (Methods B). A conservative estimate for the amplitude detection threshold is
P660P/PP equal to 0.86%.

The S660S-S410S differential measurements and global seismic tomography models
show that the visibility of P660P is associated with seismic indicators for high mantle
temperatures (Fig. S13). While the global average MTZ thickness is 244.1 ± 6.2 km, regions
of robust P660P reflections display MTZ thicknesses of 238.9 ± 3.2 km (Fig. 2b). They also
display lower velocities in the S40RTS tomographic model, particularly at depths of 700 –
1000 km (Fig. S13c). Statistical t-tests confirmed that these differences are significant (p-
value less than 0.001).

**Mineral physics modelling**

For comparison to our data, we predict full seismic waveforms for SS and PP precursors
based on mantle mineral assemblages for a range of mantle geotherms (Fig. 3a-d, S14-19)
(Methods C, D). Phase assemblages are obtained via Gibbs free energy minimisation
calculations using the thermodynamic formalism from ref. 31. We use the mineral elastic-
property database from ref. 32. Properties are computed along pressure-temperature (P-T)
adiabatic gradients, with defining potential temperatures (T<sub>pot</sub>) associated to the extrapolation
of the adiabat to zero pressure (Fig. 3a-b). Compositions are based on the chemical analysis
from ref. 33.

We initially investigate the effects of temperature and composition on noise-free
synthetic waveforms (Fig. 3c-f, S18, S19). Subsequently, we account for the exact data
acquisition geometry by simulating each waveform in the database, incorporating event
depth, epicentral distance, Centroid Moment Tensor focal mechanism solutions (34), and
noise (Methods C). These realistic synthetic waveforms are processed the same way as the
observed data (spatial partitioning and stacking), to test predictions from compositional
models (Fig. S25, S26). The synthetic catalogues represent a significant improvement over
previous one-dimensional velocity profile calculations (20), revealing important effects of
data coverage and geometry, and highlighting the effects of temperature.

Compositional models are either an equilibrium assemblage (EA) of pyrolite (33), or
mechanical mixtures (MM) of various proportions of basalt (35) and harzburgite (36).
Modelling mantle composition in terms of these two components considers heterogeneity that
is introduced through partial melting at mid-oceanic ridges, and recycling of these
heterogeneities in the deep mantle through subduction (37,38). Mantle mixing occurs due to a
combination of two mechanisms: convective shearing, which deforms distinct blobs (stirring)
into thin streaks, and solid-state diffusion. These end-member processes result in two
compositional models, even for the same pyrolitic chemical bulk composition of the mantle
(39): EA is chemically homogeneous and similar to the rock type peridotite in the upper
mantle; MM is a solid-state mixture of basalt and harzburgite in chemical disequilibrium.
Despite being based on the same mantle bulk composition, the EA and MM models for a
basalt fraction *f* ≈ 0.2 exhibit differences in their mineral assemblages at a given pressure and
temperature (*p*-T), character of related mineral reactions (38), and hence elastic properties.

This modelling approach permits us to determine the thermochemical conditions of the
mantle for which (i) the SS and PP precursors are detectable, and (ii) the relationship between
MTZ thickness and temperatures.
156

A thermometer based on P660P visibility

157

We infer the “visibility” of P660P by comparing amplitudes of the synthetic waveforms with the detection capability of our observational data (i.e., P660P/PP greater than approximately 0.86%). This is initially performed with noise-free theoretical seismic responses, and second by using predictions for the same acquisition geometry as in the datasets, incorporating the effect of noise and interfering phases.

158

159

160

161

162

The velocity profiles for various geotherms are shown in Fig. 3a, and calculated $v_p$ gradients in the $p-T$ domain in Fig. 3b. Predictions from noise-free theoretical seismic responses for pyrolite EA (Fig. 3c,d) indicates that a visible P660P is associated with an increase in relative abundance of majoritic garnet (Gt) at high temperature (Fig. 3b), which acts to broaden and depress the 660 discontinuity to 720 km depth. This regime corresponds to potential temperatures above $T_{pot}=1600$ K (e.g. inset III in Fig. 3b). We find a similar temperature-dependent P660P visibility for MM models (Fig. 3f). A notable difference is the higher visibility threshold of $T_{pot}=1700-1800$ K for MM (Fig. 3f), slightly increasing with basalt fraction $f$. In this noise-free analysis, it is not possible to define an exact temperature threshold for P660P visibility, as this depends on composition and data noise, both of which are spatially variable.

163

164

165

166

167

168

169

170

171

172

173

A prior study (23) also identified a thermal dependence of P660P visibility, but found an opposite relationship, with weak reflection amplitudes for high-T geotherms (1900 K). We attribute this discrepancy to their use of Zoeppritz equations to describe the partitioning of seismic-wave energy at first-order interfaces, rather than considering full waveforms as we do here (Methods C1; Fig. S16).

174

175

176

177

178

The MTZ thickness estimated from SS data (e.g. Fig. 3c,e) decreases monotonically as $T_{pot}$ increases for $T_{pot}$ smaller than 1800-2100 K (depending on composition), but then increases again at very high temperatures (Fig. 3e), in contrast to a pure olivine model (22).

179

180
181

The main difference between EA and MM models occurs at *T*<sub>pot</sub> greater than 1800 K, with thinner MTZ predicted for MM (Fig. 3e).

182

183

We use our mineral-physics models to convert the observed MTZ thicknesses into a global field of vertically averaged MTZ temperature. Vertical thermal gradients would have a small influence on these measurements as the thickness is primarily controlled by the 410 topography, with influences of 0.1 km/K and 0.007 km/K for the 410 and 660 respectively (Fig. 3c; Methods D). Thus, temperature estimates based on MTZ thickness are blind to thermal gradients. We find that the global average *T*<sub>pot</sub> based on pyrolite EA is approximately 1730±50 K (Fig. 4a), in agreement with estimates from a smaller SS dataset (21). Peak-to-peak variations are ±300 K. However, the EA model for these temperatures predicts P660P reflections over around 45% of the globe (Fig. S25). Noting that this does not account for detection bias and interfering phases, the prediction is clearly at odds with the observed proportion of approximately 5.0% (Fig. 2a). Furthermore, the EA model cannot account for the thinnest MTZ from SS measurements (Fig. 3e).

184

185

186

187

188

189

190

191

192

193

194

195

Performing calculations for MM with bulk pyrolitic composition (20% basalt, 80% harzburgite) resolves this discrepancy. It reduces the average *T*<sub>pot</sub> by around 100 K to 1630±50 K (Fig. 4a), in good agreement with petrological constraints from mid-ocean ridge volcanism (40), and shifts the visibility of the P660P to higher *T*<sub>pot</sub> (1700–1800 K). Thus, the thermal field for MM predicts P660P reflections over a smaller portion of the Earth (14.0%), better matching the data (Fig. S26). This computation for basalt fractions of 0.2 represents a lower temperature bound for P660P visibility, since it favours visibility at lower *T*<sub>pot</sub> than for higher basalt fraction. Fewer P660P observations would be predicted for higher basalt enrichments, in agreement with our global average basalt fraction estimate of approximately 0.35 (see next section “Effect of composition”; Figs. S22-S24).

196

197

198

199

200

201

202

203

204
A previous study (22) determined a preference for MM over EA thermal fields from SS data; however the discrimination between models relied on petrological thermo-barometry. Here, we show that complementary PP data provide a robust constraint. Our inferred temperatures generally agree with petrologic estimates for overlying hotspot temperatures, even though the former are systematically lower than the latter (Fig. S27), likely explained by lower seismic resolutions than typical plume radii.

Conversely, an investigation of thermochemistry using receiver functions (24) found that the mantle is better described by a patchwork of EA and MM composition, at odds with our preferred global MM model. This discrepancy is reconciled by considering the data resolution of both studies. Ps conversions sample continental areas beneath seismic stations with a maximum of approximately 500 km diameter in the MTZ, whereas the minimum lateral resolution of SS and PP-precursors is approximately 600 km. Conversion data are therefore sensitive to small regions in the mantle that have potentially been fully homogenized, while reflection data sample globally a mantle that corresponds to a mechanical mixture on the relevant longer wavelengths.

Along these lines, our modelling indicates that P660P visibility, in combination with MTZ thickness from SS precursors, can be used as a mantle thermometer. We propose that the marker for a very hot regime is thin-to-normal MTZ from SS data, combined with a visible P660P, and that normal temperature conditions ($T_{pot}$ of around 1600 K (40)) favour a widespread “absence” of P660P reflections.

225

## Effect of composition

Our mineral physics modelling further suggests that the global absence of P660P observations (Fig. 2a and 3f) is primarily explained by thermal rather than compositional

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228
effects. MM calculations predict that P660P amplitudes decrease with increasing basalt fractions (Fig. 3f, S18-S20). At high basalt fractions, it is essential to distinguish between signals of the Rw dissociation and of the deeper Gt breakdown (Fig. 3a-d, S17-S19). These signals are separable for basalt fractions smaller than a threshold of 0.7-0.8, at which point the shallower Rw signal becomes undetectable. In theory, the presence of P660P might therefore also constrain the amount of basalt enhancement near the 660-km discontinuity. This aspect is particularly appealing as the detection of basalt enrichment in the MTZ, i.e. a MM with $f$ greater than 0.2 (39), is consistent with previous evidence for accumulations of recycled oceanic crust at the base of the MTZ from observational and modelling studies (15,26,41).

Very high basalt fractions ($f$ greater than 0.6), however, are incompatible with our observed seismic amplitudes measurements (Figs. S20-S21). We incorporate SS and PP precursors amplitude information with MTZ thickness in a Bayesian inversion, to obtain radially averaged potential temperature and permissible basalt fractions in the MTZ (Methods E; Figs. S22-S24). Globally, the best-fitting solution occurs at $T_{pot} \approx 1625$ K and $f \approx 0.35$, with acceptable models ranging at 1450-1650 K and 0.18-0.45 (one standard deviation). For these compositions, the P660P visibility is fully controlled by the Gt dissociation (Fig. 3b) at high temperatures (Fig. 3f). In a similar inversion, locations with visible P660P (Fig. S24) favour higher temperatures (1750-1800 K), and further enhanced basalt fractions (0.4-0.45).

Overall, our data are therefore consistent with a moderately basalt-enriched MTZ, quantitatively agreeing with geodynamic models (42). We cannot rule out more extreme small-scale compositional variations with local reductions or enrichments in basalt fraction, or regions dominated by EA, especially below our resolution limit (approximately 600 km, Fig. S28). Similarly, compositional gradients across the MTZ may occur, but such variations do not significantly affect our temperature estimates and are virtually undetectable by
precursors due to their limited effects on discontinuity depths (0.5 km absolute depth perturbation for a 10% increase in basalt fraction; Fig. S24). The same holds true for compositional changes across the discontinuities. We estimate only a 6 K temperature difference for a 12% increase in *f*. Currently, the lack of thermodynamic data precludes considerations on the effects of minor components, such as volatiles or ferric iron (43,44).

259

260

## Implications for mantle convection

The observations of P660P are thus attributed to the broad Gt-decomposition dominating in the very hottest regions of the MTZ with $T_{pot}$ greater than 1800 K (Fig. 2-3). The dynamical effect of the Gt-decomposition with a positive Clapeyron slope promotes upward transfer of hot material across the 660 (6,18). Our findings for both P660P visibility and MTZ thickness suggest that such high-T and thus efficient transport are reached in only very few and small locations (Fig. 2a, 4a), detectable only in regions with the highest data coverage (Fig. S28), which the lower resolution S660S data cannot capture.

The P660P detections are primarily located beneath the Pacific Ocean, a region with a large concentration of hotspots (e.g. 45), but there is no direct correlation with individual hotspots (Fig. 4b). Our thermal field (Fig. 4b) reveals two broad provinces of elevated temperatures beneath the North (NPP) and South Pacific (SPP), separated by a central belt of “normal” temperature (CPB). The CPB includes the Pacific Superswell, a region of anomalously shallow seafloor associated with a geoid low and isotopic anomalies often ascribed to a low-density upwelling in the upper mantle (47). Interestingly, the edges of the NPP and SPP are correlated with the surface expression of 13 of the 24 Pacific hotspots (SI6, Fig. 4b, S29, S30). Many of these hotspots display intermittent activity (less than 30 Myrs), violating the classical Wilson-Morgan plume theory (3,46). Their locations at the edges of hot provinces are inconsistent with simple vertical material transfer through the MTZ.

269

270

271

272

273

274

275

276

277

278
Instead, it implies mediated mass transfer across the MTZ with (partial) failure and ponding of hot plume material for significant timescales (4,5) secondary thermal upwellings from the edges of ponding regions where thermal convective instability is advanced (6,47,48), and only the very hottest material (greater than 1800 K) rising straight through (5). Such a scenario implies an indirect link between plume generation zones in the lowermost mantle (1) and hotspot locations (14).

Since the Gt-decomposition is only prevalent at very high temperatures, the large hot regions with $T_{pot}$ between 1700-1800 K are dominated by the Rw dissociation, which impedes vertical flow (6,8,10). This outcome agrees with previous seismic studies, which demonstrate a mainly thermal control of the 660 (21,22). The Rw dissociation reaction is a good candidate to explain ponding of material near the 660 at the bottom of the MTZ, and other factors such as basalt content, subadiabaticity, and temperature-induced viscosity contrasts probably also play a key role in plume dynamics (4-6,15,37).

Thus, our results are consistent with a semi-layered style of mantle convection (49). This finding has important implications for heat and material cycling across the MTZ, and hence secular planetary cooling. Moreover, significant quantities of trapped hot materials imply that the MTZ and/or mid-mantle may be on average slightly warmer than predicted by geotherms that are rooted in the upper mantle. A dynamic barrier to flow further leads to imperfect mixing of materials (2,49). Geochemical constraints indicate that plumes carry compositionally-distinct material (3), hence plume stagnation implies that the bulk composition beneath the MTZ differs from that of the upper mantle (4,15,42,49). Integrating our results with additional geophysical constraints such as from geodynamical modelling and mineral physics experiments may shed further light on the dynamics of thermochemical plume segregation in the MTZ, in addition to models of mantle evolution.
Acknowledgments

LW acknowledges the Australian National University as host institution for the majority of
this research, as well as inspiration through discussion at the Cooperative Institute for
Dynamic Earth Research summer program. CIDER-II is funded as a "Synthesis Center" by
the Frontiers of Earth Systems Dynamics (FESD) program of NSF under grant number EAR-
1135452. We are grateful to Ian Campbell, Louis Moresi, and Eric Debayle for insightful
discussions and comments. We acknowledge support from a Discovery Early Career
Research Award (project number DE170100329), funded by the Australian Government
(LW), the National Science Foundation (NSF) under grant numbers EAR-1661985 and EAR-
1853662 (LW), the European Union's Horizon 2020 Research and Innovation programme
under the Marie Sklodowska-Curie grant agreement 793824 (BT), and Australian Research
Council grants LP170100233 and DP190102940 (JCA). Calculations were performed on the
ANU Terrawulf cluster, a computational facility developed with support from the AuScope
initiative. AuScope Ltd is funded under the National Collaborative Research Infrastructure
Strategy (NCRIS), an Australian Commonwealth Government Programme.

Author contributions

The study was conceived by LW and BT. The datasets were compiled by LW, partially
supervised by NCS. BT wrote the adaptive stacking and synthetic modelling codes, with
some contribution by LW. BT and LW ran the synthetic experiments. LW implemented the
analysis on observed data. All authors discussed the results and interpretations, and
contributed to the manuscript.
327 **Competing interests**

328 Authors declare no competing interests.
329 **Figure Legends**

330 **Figure 1.** SS and PP precursory reflections underneath the 410 and 660, and global data
331 coverage.

332 (a) Ray paths of SS, PP, and their precursors S410S, S660S, P410P, and P660P. (b) Examples of high-quality SS (event in Indian Ocean, 1 Jan 2016, station DRLN) and PP data (Fiji, 3 Jan 2017, DGAR). The precursory signals from reflections at the 410 and 660 are indicated. (c) SS bounce points of our data (Supplementary Information SI1). (d) PP bounce points.

337 **Figure 2.** Global observations of the 410 and 660 in SS and PP data.

338 (a) Classification of apparent detections of the 660 from PP data, based on common mid-
339 point stacking analysis. The contour of anomalies in black is associated with a detection
340 frequency in 50% or more of the 150 stacking parameterisations, and classified as the
341 threshold for observation. (b) Map of MTZ thickness from SS data. Robust P660P detections
342 are superimposed. (c) Depth profile taken along the Aa profile for SS data, calculated using
343 6° spherical caps. Seismicity at the North Honshu subduction zone is shown in black, and
344 seafloor elevation is included as an inset for reference. (d) The corresponding depth profile
345 for PP data reveals the effect of interfering phases near 160°E, and several actual detections
346 of the P660P (see Fig. S10 for modelling).

347 **Figure 3.** Effects of temperature on MTZ mineralogical phase transitions, and SS and
348 PP precursors.

(a) The velocity structure v_s extracted for an equilibrium assemblage of pyrolite (EA), along a set of adiabatic temperature gradients for different potential temperatures shows that the MTZ thickness varies non-linearly with temperature. Green lines and shading indicate first-
order or sharp gradients. Note that $v_p$ and $v_s$ have similarly-shaped profile. (b) Diagram of P-wave velocity gradients showing the effects of temperature on the 410 and 660. The grey background scales with the amplitude of the velocity gradient with the darkest regions associated with the sharpest and strongest gradients (insets). The shape of the gradients is a controlling parameter for reflected-wave amplitudes. (c) Set of synthetic SS seismograms calculated for a range of temperatures in a pyrolite EA mantle. (d) Similar set of synthetic PP seismograms. The peaks corresponding to the 410, 520, and 660 are labelled. Note that these waveforms are ideal and contain no noise. (e-f) Predictions of the MTZ thickness from SS data and P660P amplitude based on EA (black) and MM with variable $f$ (light-to-dark green), measured from synthetic data such as in (c) and (d). These predictions are compared with observations (histograms along the vertical axis). Note the sharp increase in MTZ thickness from SS data at 1800 K for a MM with $f=0.6$, as a consequence of the emergence of the garnet-bridgmanite reaction with larger signal amplitude (Fig. S18). Curves in panels (e-f) are sampled with a step of 100 K and obtained by linear interpolation. This adds a ±50 K uncertainty on the P660P visibility *Tpot* threshold.

**Figure 4. Thermal model of the MTZ and implications for global mantle circulation.**

(a) Histograms of global mantle temperatures in the MTZ based on seismic and mineral physics constraints. The positioning of the global average temperature and the absolute temperature threshold for the visibility of the post-garnet transition in PP data (grey and red bars) depend on the model of mantle composition (EA or MM). We use $f=0.2$ for MM, and colour shading indicates temperature. The P660P visibility threshold depends on the data noise level and composition (dashed red), but cannot extend far below around 1700 K otherwise we would over-predict the P660P visibility. (b) Composite model of the temperature in the MTZ based on S660S (colours) and P660P (black contours) observations with superimposed hotspot locations. The hotspot distribution is from ref. 45 as follows: 1.
377 Louisville; 2. Hawaii; 3. Easter; 4. Marquesas; 5. Pitcairn; 6. Austral-Cook; 7. Tahiti-Society.

378 NPP: North Pacific Province; SPP: South Pacific Province; CPB: Central Pacific Belt.
379

**Data availability**

380

Waveform data are publicly available from the IRIS Data Management Center (NSF grant EAR-1063471). Measurements of the 410 and 660 km discontinuities, the thermal map, and velocity models are available from the ISC Repository, DOI: 10.31905/7M3LMG8X.

381

382

383

**Code availability**

384

A downsampled database of synthetic seismograms for mineral physics models is available from https://zenodo.org/record/5512035. A simplified version of the adaptive stacking code is available from https://zenodo.org/record/5512805. The full database and software is available from Benoit Tauzin upon request.

385

386

387
1. Burke, K., Steinberger, B., Torsvik, T., Smethurst, M. Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the core–mantle boundary. *Earth. Planet. Sci. Lett.* **265**, 49-60 (2008).

2. Courtillot, V., Davaille, A., Besse, J., Stock, J. Three distinct types of hotspots in the Earth's mantle. *Earth. Planet. Sci. Lett.* **205**, 295-308 (2003).

3. Koppers, A., Staudigel, H., Pringle, M., Wijbrans, J. Short-lived and discontinuous intraplate volcanism in the South Pacific: Hot spots or extensional volcanism? *Geochem. Geophys. Geosys.*, **4**(10), doi:10.1029/2003GC000533 (2003).

4. Nakagawa, T., Buffett, B. Mass transport mechanism between the upper and lower mantle in numerical simulations of thermochemical mantle convection with multicomponent phase changes. *Earth. Planet. Sci. Lett.* **230**(1-2), 11-27 (2005).

5. Farnetani, C., Samuel, H. Beyond the thermal plume paradigm. *Geophys. Res. Let.*, **32**(7), doi:10.1029/2005GL022360 (2005).

6. Dannberg, J., Sobolev, S. Low-buoyancy thermochemical plumes resolve controversy of classical mantle plume concept. *Nat. Commun.* **6**, 6960 (2015).

7. Bina, C., Helffrich, G (1994). Phase transition Clapeyron slopes and transition zone seismic discontinuity topography. J. Geophys. Res. 99, 15853-15860.

8. Christensen, U (1995). Effects of phase transitions on mantle convection. *Annu. Rev. Earth. Planet. Sci.*, **23**, 65-87.

9. Tackley, P., Stevenson, D., Glatzmaier, G., Schubert, G. Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the Earth's mantle. *Nature*, **361**, 699-704 (1993).

10. Bercovici, D., Schubert, G., Tackley, P. On the penetration of the 660 km phase change by mantle downflows. *Geophys. Res. Lett.*, **20**, 2599-2602 (1993).
413

414

11. Fukao, Y., M. Obayashi. Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. *J. Geophys. Res.*, **118**, 5920-5938 (2013).

415

416

417

12. Cao, Q., van der Hilst, R., de Hoop, M.V., Shim, S. Seismic imaging of transition zone discontinuities suggests hot mantle west of Hawaii. *Science*, **332**(6033), 1068-1071 (2011).

418

419

420

421

13. 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**(3), 1223-1236 (2011).

422

423

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

424

425

426

15. Ballmer, M.D., Schmerr, N.C., Nakagawa, T., Ritsema, J. Compositional mantle layering revealed by slab stagnation at ~1000-km depth. *Science Adv.*, **1**(11), doi:10.1126/sciadv.1500815 (2015).

427

428

16. Tosi, N., Yuen, D. Bent-shaped plumes and horizontal channel flow beneath the 660 km discontinuity. *Earth. Planet. Sci. Lett.*, **312**(3-4), 348-359 (2011).

429

430

431

17. Jenkins, J., Cottaar, S., White, R., Deuss, A. Depressed mantle discontinuities beneath Iceland: evidence of a garnet controlled 660 km discontinuity? *Earth. Planet. Sci. Lett.*, **433**, 159-168 (2016).

432

433

434

18. Hirose, K. Phase transitions in pyrolitic mantle around 670-km depth: implications for upwelling of plumes from the lower mantle. *J. Geophys. Res.*, **107**, B4, doi: 10.1029/2001JB000597 (2002).

435

436

19. Weidner, D.J., Wang. Y. Chemical- and Clapeyron-induced buoyancy at the 660 km discontinuity. *Journal of Geophysical Research: Solid Earth*, **103**(B4), 7431-7441 (1998).
437

20. Shearer, P., Flanagan, M. Seismic velocity and density jumps across the 410- and 660-kilometer discontinuities. *Science*, **285**, 1545–1548 (1999).

438

439

21. Houser, C., Williams, Q. Reconciling Pacific 410 and 660 km discontinuity topography, transition zone shear velocity patterns, and mantle phase transitions. *Earth. Planet. Sci. Lett.*, **296**(3-4), 255-266 (2010).

440

441

442

22. Ritsema, J., Xu, W., Stixrude, L., Lithgow-Bertelloni, C. Estimates of the transition zone temperature in a mechanically mixed upper mantle. *Earth. Planet. Sci. Lett.*, **277**(1-2), 244-252 (2009).

443

444

445

23. Deuss, A., Redfern, S., Chambers, K., Woodhouse, J. The nature of the 660-kilometer discontinuity in Earth's mantle from global seismic observations of PP precursors. *Science*, **311**(5758), 19-201 (2006).

446

447

448

24. Munch, F., Khan, A., Tauzin, B., van Driel, M., Giardini, D. Seismological evidence for thermo-chemical heterogeneity in Earth's continental mantle. *Earth. Planet. Sci. Lett.*, **539**, 1162-1169 (2020).

449

450

451

25. Schmerr, N., Garnero, E.J. Upper mantle discontinuity topography from thermal and chemical heterogeneity. *Science*, **318**(5850), 623-626 (2007).

452

453

26. Wu, W., Ni, S., Irving, J.C.E. Inferring Earth's discontinuous chemical layering from the 660-kilometer boundary topography. *Science*, **363**(6428), 763-740 (2019).

454

455

27. Waszek, L., Schmerr, N., Ballmer, M. Global observations of reflectors in the mid-mantle with implications for mantle structure and dynamics. *Nat. Commun.*, **9**(385) (2018).

456

457

28. Bassin, C., Laske, G., Masters, G. The current limits of resolution for surface wave tomography in North America. *EOS Trans AGU*, **81**, F897 (2000).

458

459

29. Dziewoński, A., Anderson, D. Preliminary reference Earth model. *Phys. Earth Planet. Int.*, **25**, 297-356 (1981).

460
461

462

30. Lessing, S., Thomas, C., Saki, M., Schmerr, N., Vanacore, E. On the difficulties of detecting PP precursors. *Geophys. J. Int.*, **201**(3), 1666-1681 (2015).

463

464

465

31. Connolly, J.A.D. Computation of phase equilibria by linear programming: a tool for geodynamic modeling and its application to subduction zone decarbonation. *Earth. Planet. Sci. Lett.*, **236**, 524-541 (2005).

466

467

32. Stixrude, L., Lithgow-Bertelloni, C. Thermodynamics of mantle minerals – II. Phase equilibria. *Geophys. J. Int.*, **184**(3), 1180-1213 (2011).

468

469

470

33. Baker, M.B., Stolper, E.M. Determining the composition of high-pressure mantle melts using diamond aggregates. *Geochimica et Cosmochimica Acta*, **58**(13), 2811-2827 (1994).

471

472

34. Ekström, G., Nettles, M., Dziewonski, A. The global CMT project 2004-2010: Centroid-moment tensors for 13,017 earthquakes. *Phys. Earth Planet. Inter.*, **200-201**, 1-9, (2012).

473

474

35. Gale, A., Dalton, C.A., Langmuir, C.H., Su, Y., Schilling, J. The mean composition of ocean ridge basalts. *Geochem. Geophys. Geosys.*, doi:10.1029/2012GC004334 (2013).

475

476

477

478

36. Afonso, J.C., Fernandez, M., Ranalli, G., Griffin, W.L., Connolly, J.A.D. Integrated geophysical-petrological modeling of the lithosphere and sublithospheric upper mantle: Methodology and applications. *Geochem. Geophys. Geosys.*, doi:10.1029/2007GC001834 (2008).

479

37. Ringwood, A. E. A model for the upper mantle. *J. Geophy. Res.*, **67**(2), 857-867 (1962).

480

481

38. Stixrude, L., Lithgow-Bertelloni, C. Geophysics of chemical heterogeneity in the mantle. *Ann. Rev. Earth. Planet. Sci.*, **40**, 569-595 (2012).

482

483

39. Xu, W., Lithgow-Bertelloni, C., Stixrude, L., Ritsema, J. The effect of bulk composition and temperature on mantle seismic structure. *Earth. Planet. Sci. Lett.* **275**, 70-79 (2008).
484

485

486

40. Herzberg, C. et al. Temperatures in ambient mantle and plumes: Constraints from basalts, picrites, and komatiites. *Geochem. Geophys. Geosys.*, **8**(2), doi:10.1029/2006GC001390 (2007).

487

488

41. Hirose, K., Fei, Y., Ma, Y., Mao, H.K. The fate of subducted basaltic crust in the Earth's lower mantle. *Nature*, **397**(6714), 53-56 (1999).

489

490

491

42. Yan, J., Ballmer, M.D., Tackley, P.J. The evolution and distribution of recycled oceanic crust in the Earth's mantle: Insight from geodynamic models. *Earth. Planet. Sci Lett.*, **537**, 1161-1171 (2020).

492

493

43. Thio, V., Cobden, L., Trampert, J. Seismic signature of a hydrous mantle transition zone. *Phys. Earth. Planet. Int.*, **250**, 46-63 (2016).

494

495

44. Schulze, K. et al. Seismically invisible water in Earth's transition zone? *Earth. Planet. Sci. Lett.*, **498**, 9-16 (2018).

496

497

498

45. Müller, R.D., Royer, J.Y., Lawver, L.A. Revised plate motions relative to the hotspots from combined Atlantic and Indian Ocean hotspot tracks. *Geology*, **21**(3), 275-278 (1993).

499

500

46. Kim, S.S., Wessel, P. New global seamount census from altimetry-derived gravity data. *Geophys. J. Int.*, **186**(2), 615-631 (2011).

501

502

47. Davaille, A. Simultaneous generation of hotspots and superswells by convection in a heterogeneous planetary mantle. *Nature*, **402**, 756-760 (1999).

503

504

48. King, S., Ritsema, J. African hot spot volcanism: small-scale convection in the upper mantle beneath cratons. *Science*, **290**(5494), 1137-1140 (2000).

505

506

49. Kellogg, L.H., Turcotte, D.L. Mixing and the distribution of heterogeneities in a chaotically convecting mantle. *J. Geophys. Res. Solid Earth* **95**, B1: 421-432 (1990).
507

50. Hofmann, A. Mantle geochemistry: the message from oceanic volcanism. *Nature*, **385**,

508

219-229 (1997).

509

510
511

# Methods

512

## A. Adaptive Stacking (AS)

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

The PP and SS waveforms require dedicated processing to extract weak reflections from the 410 and 660 (amplitude of 0.4-6% the reference phase) from the noise. The typical approach is grouping the data in bins at the Earth's surface and stacking the waveforms. Here we refine this approach, called adaptive stacking (AS), which extracts information about the interface depth and amplitude of reflections with uncertainties. The Earth is partitioned into cells with irregular shapes called Voronoi cells. This parameterisation adapts its shape to account for sharp variations in topography. The initial step randomly distributes Voronoi nuclei throughout the globe. The distribution is uniform on the sphere (51), and the initial number of cells (2000) chosen so the surface area covered by an individual cell is smaller than the area of the first Fresnel zone of a single datum. The data are then grouped into the Voronoi cells, according to the location of the reflection points of the reference phase. The grouped time series are converted to depth using PREM (29) after correction for 3D structure, weighted by signal-to-noise ratio (SNR), and stacked. The maximum stacked-amplitude within a specific depth range surrounding a discontinuity (410 or 660) is defined as the discontinuity depth. The depth windows are 350-480 and 380-440 km for the 410 from SS and PP data respectively, and 600-730 km (SS) and 610-750 km depth (PP) for the 660. These ranges were chosen based on synthetic experiments (Fig. 3), and an analysis of the detection performance in the observed datasets. In particular, they account for interference effects from reflections at the 520-km discontinuity, plus the deeper signals from the post-garnet phase disasocciation (Fig. 3). Bootstrap resampling provides errors on depth and amplitude estimates (52), which are estimated from one standard deviation of a distribution of 100 samples. In the AS, Voronoi cells with post-stacking standard deviations larger than 10-km are removed from the global partition. This procedure iterates until convergence, i.e. when all
cells have a standard error below 10-km, producing an optimised tessellation for each run. A
visual output of the iterative process is presented in Figure S2. A single run produces patterns
with sharp edges, and the spatial resolution of the topography variations is influenced by the
cell geometry and therefore the initial distribution of the Voronoi cells. Thus, we repeat the
process 150 times to obtain measurements across a variety of cell distributions and
geometries. We use the outputs of these 150 runs to generate smooth average spatial
distributions over a grid of 0.5x0.5-degree latitude and longitude, producing detailed global
maps of discontinuity depths and MTZ thickness (Fig. 2, S3-5).

544

## B. Observations of P660P phases and consistency of their timing and slowness

AS cannot be implemented for the P660P because it is weak or absent globally. An absent
detection might arise either from the structure of velocity gradients near 660-km depth (Fig.
3a,b,d), or noise in the data. We confirmed with synthetic experiments that topography on the
P660P has limited effect at the frequencies of our data (SI5.3). To avoid a null observation
due to high noise, we minimise its influence by stacking P660P over the parameterisation
constructed for the P410P phase. Because this parameterisation is optimal for stacking
constructively for the 410-km discontinuity topography in terms of data coverage and quality
(by suppressing the noise), we expect that it will be similarly successful for the 660-km
discontinuity topography. We validated this approach via the SS data with a result shown in
Figure S5. Topographies obtained from stacking the S660S phase based on the S410S and the
S660S Voronoi-parameterisations are in good agreement (correlation coefficient R=+0.88).

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

For P660P detection, we implement a criterion to identify locations with a clear single pulse:
Voronoi cells with a standard error for depth of greater than 10 km are defined as “null-
observations”. We retain locations in which P660P is detected with at least a 50% probability
across the 150 iterative procedures (Fig. 2a, S9, S4d, S6d).
After global detection of the candidate P660P reflections (Fig. 2a), we performed slant-stacking on subsets of the PP data in these regions (Fig. S9) for quality evaluation. We show vespagrams for positive and spurious-detections in Figure S11. When observed (patches 1-4), the phases have a relative timing of -125 s and a differential slowness of -0.3 s/deg with respect to the PP phase, with no obvious interfering phases. These values are consistent with P-wave reflections underneath the 660, as demonstrated by synthetic modelling (Methods C).

However, in many cases (patches 5-14), the detections are ambiguous because the slowness is poorly resolved and/or the presence of large-amplitude interfering phases.

Systematic analysis of the depth-epicentral distance sections (Figure S12) reveals that the epicentral distance coverage for these patches is quite inhomogeneous, favouring the contribution of interfering phases that mimic P660P (typically PKP; Fig. S12-3, S12-4).

Synthetic modelling from mineral physics (see Methods C) helps distinguish true P660P detections because the synthetic data incorporate the same acquisition geometry as observed data. We compare observed vespagrams (Fig. S11) and depth-distance sections (Fig. S12) to predictions using a pyrolitic compositional model, and two reference potential temperatures 2000 K and 1400 K. These temperatures were selected as extrema of P660P visibility. A temperature of 1400 K does not generate a detectable P660P signal, although interfering phases may simulate such an arrival if epicentral distance distribution is poor.

Conversely, a temperature of 2000 K consistently generates a large P660P signal. These experiments show that with good epicentral distance coverage, the data should be able to discriminate between these two temperature regimes. Our criteria for a robust detection (Fig. 2a) are therefore: an even, broad distribution of epicentral distance in the time-distance profiles without strong interfering phases, a signal at the relative predicted time and slowness in the observed and 2000 K synthetic vespagrams, and no signal in the 1400 K vespagrams.
Bootstrap resampling is used to obtain an estimate of the 95% confidence level, to determine
whether a signal is significant.

We estimate the average noise level for the 660 depth window via the mean standard
deviation on the average trace from bootstrap resampling. The average P660P/PP amplitude
ratio is 1.42±0.23%, with an average noise of 0.26±0.11%. At locations with no P660P we
find noise of 0.65±0.21%. A conservative estimate for the minimum detectable P660P
amplitude is therefore 0.86%.

### C. Seismic Modelling

C.1. *Frequency response to complex velocity gradients.* Zoeppritz equations are commonly
used to compute reflection/conversion coefficients for first-order seismic interfaces
(20,23,53,54). This approach shows limitations when the discontinuity extends over a
significant depth interval. For the MTZ phase transformations, the use of
reflection/transmission coefficients for first-order interfaces is valid only at long period
(55,56). To demonstrate the effect of the complex shape of velocity gradients on PP
underside reflections (i.e. PP precursors), we reconstructed their frequency response by sub-
dividing the elastic gradient in thin homogeneous layers, then computing the effective elastic
response using a plane-wave matrix algorithm (57-59). We modified an algorithm from ref.
58 to calculate the generalized reflection-transmission coefficients (56,60).

The result for P-wave reflections underneath various types of structures for the 660
km discontinuity is shown in Figure S16. We show examples of two elastic models for a
pyrolitic composition and 1200 and 2100 K adiabatic temperature gradients, and of the
*iasp91* Earth's reference model (61). The *iasp91* structure is a simple first-order discontinuity
with approximately 6% *v*<sub>p</sub> and *v*<sub>s</sub> jumps and 7.5% density contrast (Fig. S16a). For this
structure, the reflection coefficient reaches 4% in the typical 8-40 s period range of analysis,
but shows little variation with frequency. Conversely, the responses in frequency of mineral physics-based models display significant deviation from those of a first-order interface. At long periods, the elastic gradients appear as a first-order discontinuity to seismic waves. The reflection coefficient is maximum for the 2100 K elastic profile because the overall contrast in elastic properties is larger for this model. The 1200 K model is the least efficient for reflecting long-period P-waves.

For the two mineralogical models, the shape of the response – a low-pass filter – is typical of that from an extended elastic gradient: the larger the gradient, the lower the cut-off filter frequency (Fig. S16b). For SS precursors and other phases, the response of these structures is instead flat within the band-pass of analysis of the data (not shown here). Particularly noteworthy is that the roll-off of the filter for these structures for PP occurs at the lower limit of the typical frequency range for analysis of PP precursors. Thus, the reflective characteristics of the 660 for PP-like waves have complex fluctuations depending on the frequency content of the data and the wave slowness. In the frequency range of analysis, and for the same jump in elastic properties, the reflections have lower-amplitudes for gradual boundaries than for a single first-order interface. The 2100 K model is more reflective than the 1200 K model (1% for the 1200 K model and 2% for the 2100 K model, versus 4% for iasp91).

### C.2 A Database of Full-Waveform Synthetic Seismograms.
Accurate seismic modelling is crucial to test the effects of the processing on time-domain signals, and the influence of stacking seismic data with different source mechanisms and source-receiver geometry. The acquisition geometry controls in particular interference patterns from other seismic phases. In addition, seismic modelling allows us to quantitatively investigate the effect of temperature on the waveforms. We use a reflectivity algorithm (62,63) to model the
complete waveform on synthetic vertical and transverse components for each seismogram in
our databases, performed in one-dimensional (1D) Earth's seismic models.

Instead of simulating wave-propagation for the 2036 PP and 1666 SS sources, we
generate a library of seismograms for five elementary moment tensors (Fig. S14) taken as
sources located at depths from 2.5 km to 75 km by step of 2.5 km. Seismograms are
reconstructed from a linear combination of these elementary seismograms (64). An example
of a reconstruction and comparison to observed data is presented in Figure S14. We
confirmed for arbitrary moment tensors that the reconstructions provide the same waveforms
as when computing the seismograms using these moment tensors (Fig. S14). In each
simulation, 200 receivers are located at distances from 80° to 179.5°, every 0.5°. This reduces
the number of simulations to 150 per seismic models. In the library of elementary
seismograms, the sampling precision in the event depths and epicentral distances introduce an
error in the reconstruction. However, the effect is small as the measurements are ultimately
performed for stacked seismograms. We add random noise to the synthetics (Fig. S15), from
a uniform distribution within 1-30% the amplitude of the reference phase. This procedure
reproduces the noise effect in the observed data. The noise is filtered for the same frequencies
before application. The synthetic databases are then processed in a similar way to the
observed data (AS, with the same procedure for retaining locations with P660P signals) for
comparison.

D. Elastic Models from Mineral Physics

The reference compositional model in our study is pyrolite, which is a model composition for
the Earth's upper mantle (37). Its phase transformations are in qualitative agreement with
observed seismic discontinuities (Fig. 3a). The composition is based on the analysis from ref.
33 with weight percentages of oxides: 45.5% SiO₂, 3.98% Al₂O₃, 7.18% FeO, 38.3% MgO,
and 3.57% CaO. Qualitative comparison with other studies using slightly different pyrolitic compositions does not show obvious difference in the phase diagram (22,32,39). Phase assemblages are obtained via Gibbs free energy minimisation with the code Perple_X (31), using the thermodynamic formalism and mineral elastic properties database ($v_p$, $v_s$, and density) from ref. 32. This computation supposes that the chemistry of the mantle is fully equilibrated (Equilibrium Assemblage; EA). We test other models based on mechanical mixtures (MM) of harzburgite and basalt (the properties of these two lithologies are also computed in the same manner). A mixture of these two harzburgite-basalt compositions is supposed to represent at first order mantle material that has undergone complete differentiation at mid-oceanic ridges and is cycled back into the mantle at subduction zones.
We test various mechanical mixtures parameterized in terms of fraction of basalt. The overall physical properties of the assemblages are shown in Figure S17, computed from a weighted average of the constituents (Voigt average). We find differences in elastic properties computed with Reuss, Voigt, and Voigt-Reuss-Hill averages to be insignificant.

Properties are computed along pressure-temperature adiabatic gradients, with corresponding potential temperatures associated with the extrapolation of the adiabat to the Earth's surface. The pressure step size is sampled so that the elastic profiles are discretized with approximately 4 km depth interval. This allows us to accurately model frequency-dependent waveforms within the frequency range of analysis. We compute 1D elastic profiles for 19 adiabats with potential temperatures ranging from 600 K to 2400 K in intervals of 100 K (Fig. 3b and S17). The shear and compressional seismic attenuations are fixed to the attenuation profiles in PREM (20). A difference in the viscoelastic behaviour should be expected from accounting for grain boundary sliding through extended Burgers' models (65), but synthetic experiments (e.g. 66) showed that this effect is small relative to PREM attenuation. The overall effect of increased attenuation is mainly seen for shear waves (i.e.
SS) at high temperatures, with less impact on PP waveforms. We simulate teleseismic wave-
propagation in each of these 19 models, building the libraries of seismograms as described
above to generate synthetic data corresponding to every seismogram in the observed
databases.

Theoretical PP and SS waveforms for the reference pyrolitic composition (EA) and
various mixtures of basalt and harzburgite (MM) are shown in Figure 3c-d and S18-S19,
respectively. The trends in the thermal response of the P660P reflection are consistent across
the compositional models, with low reflection amplitudes at intermediate-high temperatures
(below 1700-1800 K), and higher amplitudes above this temperature threshold (Fig. S19). In
comparison, SS data do not show obvious temperature-dependent variations in reflection
amplitudes (Fig. S18, S22). Increasing the basalt fraction decreases the overall impedance
contrast across the phase changes between olivine polymorphs, and increases the impedance
contrast across transitions in the garnet system (Fig. S17). This behaviour is mapped into the
amplitudes of pulses of PP and SS waves reflected at the velocity gradients (Fig. S18-S19).

Uncertainties in predicted MTZ thickness and depth gradients for seismic velocities derived from thermodynamic calculations have been discussed in a number of works (e.g.41,67-70, and references therein). These studies concluded that the main sources of intrinsic error (i.e. those not related to uncertainties in T and P) come from uncertainties in bulk composition, stability of mineral phases and elastic properties of some end-member minerals, especially the garnet and pyroxene end-members; the latter two being generally more important. According to previous results (69,70), we estimate uncertainties (one standard deviation) of the order of ±5 km for the MTZ thickness and approximately 10% for the depth gradients (dv_p/dP).

Predictions for amplitudes of P410P, S410S, and S660S are shown in Figure S20, and
P660P amplitudes in Figure 3f, and compared with global observations. The observed S410S,
S660S, and P660P amplitude distributions all fall within the range predicted by basalt fractions of 0-0.6. Observed P410P amplitudes are slightly larger than predicted. This discrepancy between observations and theory may arise from variations in chemistry, uncertainties in the thermoelastic mineral properties, inaccuracies in tomographic corrections, or in the seismic forward modelling.

### E. Joint inversion for temperature and composition

In order to determine the temperature and composition that best match our combination of SS and PP measurements, we perform statistical inversions of the MTZ thickness, and amplitudes of S410S, S660S, P410P, and P660P, exploring the entire model space of ($T_{pot}, f$). These data have very different sensitivity to $T_{pot}$ and $f$. Figure S22 shows that the S410S, S660S, and P410P amplitudes depend almost entirely on mantle composition, whereas for 0-0.6 basalt fractions, the MTZ thickness is controlled by temperature. P660P amplitudes have a more complicated dependence on both $T_{pot}$ and $f$. We also show in Figure S22 the range of properties that satisfies each data individually, for regions away from clear P660P observations (i.e. 95% of the Earth's surface).

These seismic observables are then used to perform a probabilistic inversion for ($T_{pot}, f$) (Fig. S23, S24). For this, we consider (i) locations away from a clear P660P signal (amplitudes below the noise level of 0.65-0.21=0.44%), and (ii) regions with clear P660P observations. We use a Bayesian formulation, the posterior distribution on model parameters given the data d being such that:

$$p(m|d) \propto p(d|m) p(m),$$

where $p(m)$ is the a priori distribution on model parameters ($T_{pot}, f$), and $p(d|m)$ is the likelihood of observing the data d (SS, PP precursor amplitudes or the MTZ thickness) from ($T_{pot}, f$). We use uniform priors within [0;1] and [1000;2400] K for $f$ and $T_{pot}$. The likelihood
represents the misfit between observed $d^{obs}$ and predicted data d:

$$p(d|m) = \exp(-1/2 (d^{obs}-d)^2 / \sigma^2),$$

$d^{obs}$ being either the average MTZ thickness or SS/PP precursor amplitudes over the considered regions, and $\sigma$ representing one standard deviation on spatial variations over the region. Figures S23 and S24 show individual likelihoods for each data taken separately, with the highest probability centred on the range of properties that satisfies the data (Fig. S22).

The joint posterior probability distribution (Fig. S23b, S24b) is then computed by $p(m|d) = \prod p_i(m|d_i)$, where i is the index of the dataset (amplitude or thickness). These posterior distributions are approximately 2D Gaussian functions with a well-defined ($T_{pot}, f$) maximum and a small positive trade-off. Best-fitting solutions and standard deviations on $T_{pot}$ and $f$ are measured from the marginal posterior distributions (Fig. S23b, S24b).
746 Methods References

747 51. Marsaglia, G. Choosing a point from the surface of a sphere. *The Annals of Mathematical Statistics*, **43**(2), 645-646 (1972).

748

749 52. Efron, B., Tibshirani, R. An introduction to the bootstrap. *Chapman and Hall, New York*. doi:10.1007/978-1-4899-4541-9 (1991).

750

751 53. Zoeppritz, K. On the reflection and penetration of seismic waves through unstable layers. *Göttinger Nachrichten* **1**, 66-84 (1919).

752

753 54. Aki, K, Richards., P.G. *Quantitative seismology: Theory and methods*. New York, p 801 (1980).

754

755 55. Bostock, M.G. Seismic waves converted from velocity gradient anomalies in the Earth's upper mantle. *Geophys. J. Int.*, **138**(3), 747-756 (1999).

756

757 56. Tauzin, B., Pham, T.S., Tkalčić, H. Receiver functions from seismic interferometry: A practical guide. *Geophys. J. Int.*, **217**(1), 1-24 (2019).

758

759 57. Kennett, B.N.L. *Seismic wave propagation in stratified media*. ANU Press, p 288 (1983).

760 58. Chapman, C.H. Yet another elastic plane-wave, layer-matrix algorithm. *Geophys. J. Int.*, **154**(1), 212-223 (2003).

761

762 59. Chapman, C.H. *Fundamentals of seismic wave propagation*. Cambridge university press (2004).

763

764 60. Ma, Y., Wang, R., Zhou, H. A note on the equivalence of three major propagator algorithms for computational stability and efficiency. *Earthquake Sci.* **25**(1), 55-64. (2012)

765

766

767 61. Kennett, B.N.L., Engdahl, E. Traveltimes for global earthquake location and phase identification. *Geophys. J. Int.* **105**(2), 429-465 (1991).

768

769 62. Fuchs, K., Müller, G. Computation of synthetic seismograms with the reflectivity method and comparison with observations. *Geophys. J. Int.*, **23**(4), 417-433 (1971).

770
771

63. Müller, G. The reflectivity method: a tutorial: *J. Geophysics*, **58**, 153-174 (1985).

772

64. Kikuchi, M., Kanamori, H. Inversion of complex body waves—III. *Bull. Seismo. Soc. Am.*, **81**(6), 2335-2350, (1991).

773

774

65. Faul, U.H., Jackson, I. The seismological signature of temperature and grain size variations in the upper mantle. *Earth. Planet. Sci. Lett.*, **234**(1-2), 119-134 (2005).

775

776

66. Munch, F.D., Khan, A., Tauzin, B., Zunino, A., Giardini, D. Stochastic Inversion of P-to-S Converted Waves for Mantle Composition and Thermal Structure: Methodology and Application. *J. of Geophys. Res.: Solid Earth*, **123**(12), 10-706 (2018).

777

778

779

67. Connolly, J.A.D., Kerrick, D.M. Metamorphic controls on seismic velocity of subducted oceanic crust at 100–250 km depth. *Earth Planet. Sci. Lett.*, **204**, 61-74 (2002).

780

781

68. Afonso, J. C. et al. 3-D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle I: A priori information and geophysical observables. *J. Geophys. Res. Solid Earth*, **118**, 2586-2617 (2013).

782

783

784

69. Stixrude, L., Jeanloz, R. Constraints on Seismic Models from Other Disciplines - Constraints from Mineral Physics on Seismological Models. 10.1016/B978-0-444-53802-4.00022-1 (2015).

785

786

787

70. Connolly, J. A. D., Khan, A. Uncertainty of mantle geophysical properties computed from phase equilibrium models. *Geophys. Res. Lett.*, **43**, 5026-5034 (2016).

788
<figure><img src="image_2.png" /><figcaption>a</figcaption></figure>

SS bouncepoints

<figure><img src="image_3.png" /><figcaption>c</figcaption></figure>

PP bouncepoints

<figure><img src="image_4.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>c</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_10.png" /></figure>
<caption>a</caption>

Vs elastic profiles

<figure><img src="image_11.png" /><figcaption>b</figcaption></figure>

Vp gradients

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

SS data

Mantle *Tpot* (K)

<figure><img src="image_13.png" /><figcaption>d</figcaption></figure>

PP data

Mantle *Tpot* (K)

<figure><img src="image_14.png" /><figcaption>e</figcaption></figure>

MTZ thickness from SS data

Mantle *Tpot* (K)

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

P660P visibility

Mantle *Tpot* (K)

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

Compositional
models

Fraction of basalt in Mechanical Mixture

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>b</figcaption></figure>

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