Global mantle structure from multifrequency tomography using P , PP and P -diffracted waves

core-diffracted (Pdiff) constrain lower study is the ﬁrst to invert a very large data set of Pdiff waves, up to the highest possible frequencies. This results in tomographic resolution matching and exceeding that of global S -wave tomographies, which have long been the models of choice for interpreting lowermost structure. We present three new global tomography models of 3-D isotropic P -wave velocity in the earth’s mantle. Multifrequency cross-correlation traveltimes are measured on all phases in passbands from 30 s dominant period to the highest frequencies that produce satisfactory ﬁts ( ≈ 3 s). Model DETOX-P1 ﬁts ≈ 2.5 M traveltimes from teleseismic P waves. DETOX-P2 ﬁts the same data, plus novel measurements of ≈ 1.4 M traveltimes of Pdiff waves. DETOX-P3 ﬁts the same data


I N T RO D U C T I O N
Core-diffracted P and S waves (Pdiff and Sdiff) are body waves that dive so deep that they sense the earth's core and are diffracted around it. Their abundance and relatively short wavelengths make them the most highly resolving wave type that extensively samples the lower third of the mantle. Since the lowermost mantle hosts the thermal boundary layer expected to generate mantle plumes, and also holds the oldest memory of subducting seafloor, the constraints added by core-diffracted body waves would be highly desirable for geodynamics and palaeo-reconstructions. Yet these waves have seen very little use in tomography. Large data sets of waveformbased Sdiff measurements have been inverted in only a few global shear-wave tomographies (Mégnin & Romanowicz 2000;Ritsema et al. 2004Panning & Romanowicz 2006, and the use of Pdiff waves has been much more limited. The only global P-wave tomographies to explicitly measure and model Pdiff traveltimes were by Wysession (1996), Kàrason & Van der Hilst (2001) and Li et al. (2008), who all used the same data set of 543 differential PKP-Pdiff measurements by Wysession (1996) and very approximative sensitivity kernels. Hence the lower third of the mantle has remained poorly illuminated especially by P waves.
The waveform-based measurement of Pdiff waves is challenging (e.g. Wysession 1996;Hosseini & Sigloch 2015;Euler & Wysession 2017) because high frequencies are selectively damped with increasing epicentral distance, resulting in poor signal-to-noise ratio. Even more challenging is the modelling of their wave propagation and non-ray-like sensitivities, requiring computationally expensive forward modelling techniques (Nissen-Meyer et al. 2007a,b;Liu & Tromp 2008;Igel 2017;Leng et al. 2019). As a result, the resolution of P-wave tomographies in the lower third of the mantle has been much poorer than in the middle third, the latter illuminated by teleseismic waves turning at mid-mantle depths. S-wave tomographies do not suffer from the same blind spot because they can illuminate the lowermost mantle with normal modes (e.g. Masters et al. 1996;Ritsema et al. 1999Ritsema et al. , 2011Durand et al. 2016;Koelemeijer et al. 2016), albeit at relatively low spatial resolution. In prior work most similar to ours, Sdiff waves have been measured and inverted in a waveform-based framework, although in smaller numbers and for longer wavelengths than used here [e.g. 1220 Sdiff and sSdiff measurements at periods >31.4 s in Mégnin & Romanowicz (2000); measurements from 253 earthquakes at periods >17 s in Bozdag et al. (2016)]. Our work aims at filling the P-wave resolution gap in the lower third of the mantle and achieving the highest resolution of any tomography, thanks to the relatively shorter wavelengths of Pdiff waves compared to Sdiff.
In Hosseini & Sigloch (2015), we presented a method to routinely measure finite-frequency traveltimes of Pdiff waves by crosscorrelating observed waveforms with synthetic seismograms across the broad-band frequency range. We measured a very large data set of 479 559 core-diffracted P-wave traveltimes, in frequency passbands ranging from 30.0 to 2.7 s dominant period. A 'prototomography' that projected Pdiff traveltime anomalies onto their core-grazing ray segments was able to recover major, known structural features along the core-mantle boundary (CMB) ( fig. 15 in Hosseini & Sigloch 2015). In the actual global tomographies presented here, Pdiff measurements have been extended to ≈1.4 M traveltimes (Table 1). Jointly with finite-frequency traveltimes of ≈2.5 M teleseismic P and ≈1.2 M PP phases, the Pdiff data are used in a broad-band, waveform-based inversion for P-velocity structure of the entire mantle.
Waveform-based tomography methods explicitly account for wavefront healing through physically realistic, finite-volume sensitivities (Fresnel zones) of seismic observables (such as traveltimes or amplitudes), away from the geometrical ray (Marquering et al. 1998;Dahlen et al. 2000;Hung et al. 2000;Montelli et al. 2004a;Tromp et al. 2005;Tape et al. 2007;Fichtner et al. 2008;Bozdag et al. 2016). The first global P-wave tomography that was not ray-based was presented by Montelli et al. (2004a), and used the more realistic, but computationally efficient 'banana-doughnut' sensitivity kernels derived by Dahlen et al. (2000). A method to measure finite-frequency traveltimes and amplitudes tailored to the banana-doughnut kernels was only developed later by Sigloch & Nolet (2006), so Montelli et al. (2004aMontelli et al. ( , 2006 used these kernels in conjunction with ≈88 000 pre-existing, correlation-based measurements of teleseismic P and PP waves at a dominant period of 20 s; together with 1.5 million analyst-picked, teleseismic P and PP traveltimes ) that were modeled by seismic rays. Obayashi et al. (2013) was the first, and so far only, global P-wave tomography obtained from self-consistent finite-frequency measurements and kernels (for >20 K teleseismic differential traveltime measurements at neighbouring stations in the Pacific hemisphere), which they complemented by 10 M analyst-picked teleseismic traveltimes ) modeled as rays. All other global P-wave tomographies have been based on ray theory (e.g. Dziewonski et al. 1977;Dziewonski 1984;Zhou 1996;Grand et al. 1997;Obayashi & Fukao 1997;Bijwaard et al. 1998;Kennett et al. 1998;Boschi & Dziewonski 2000;Kàrason & Van der Hilst 2001;Zhao 2004;Lei & Zhao 2006;Amaru 2007;Li et al. 2008;Simmons et al. 2012;Koelemeijer et al. 2016), benefiting from the simplicity and computational efficiency of ray-based modelling, and from many millions of analyst-picked, high-frequency traveltimes from stations deployed since 1960, as curated and distributed by the International Seismological Centre Weston et al. 2018). Since the sensitivities of Pdiff waves cannot be reasonably approximated by rays, analystpicked P arrivals at epicentral distances exceeding ≈97 • (the range beyond which Pdiff nominally emerges from teleseismic P) have not been included in tomographies.
In this paper, we present and discuss three global, waveformbased tomography models using very large data sets of traveltime measurements of P, PP and Pdiff waves by cross-correlation Table 1. Summary of the four different P-wave traveltime data sets used in the tomographic inversions for models DETOX-P1, P2 and P3. The first three are multifrequency cross-correlation traveltimes measured on P, PP and Pdiff seismograms in up to eight frequency passbands. Fourth row describes a supplementary set of analyst-picked P-wave arrival times from the ISC-EHB/EHB catalogues Weston et al. 2018 (Table 1). Section 2 describes the multiple-frequency crosscorrelation measurements (periods 30 to 3 s), and the analyst-picked arrival times carefully selected to close coverage gaps in the (relatively recent) waveform data. In order to account for diffraction of the deepest-diving P waves around the earth's core, Green's functions for all P phases are calculated by the fully numerical wave propagation tool AxiSEM  and Instaseis (van Driel et al. 2015), rather than by the approximate methods suitable for teleseismic (Chapman 1978;Sigloch & Nolet 2006) or triplicated waves (Fuchs & Müller 1971;Stähler et al. 2012). AxiSEM is efficient enough to afford modelling waveforms to shorter periods (1 s) than the data can routinely fit (3 s), that is, accepting a spherically symmetric background model enables full-waveform inversion that is not computationally limited. Section 3 describes the methodology for joint inversion of P, PP and Pdiff traveltimes. The finite-frequency sensitivity kernels in this study are based on the method by Dahlen et al. (2000), which for Pdiff waves is significantly more approximate than AxiSEM kernels (those are still in the debugging stage), but significantly less approximative than Pdiff sensitivities used previously (Wysession 1996; Kàrason & Van der Hilst 2001;Li et al. 2008). Section 4 appraises model resolution and a smearing artefact under the oceans caused by sparse PP paths. Section 5 describes and discusses low-velocity anomalies in the mantle (LLVPs and mantle plumes), and Section 6 does the same for fast-velocity anomalies (subducted oceanic lithosphere or 'slabs'). Sections 7 and 8 give further discussion and conclusions.

Multifrequency cross-correlation traveltimes
A finite-frequency traveltime anomaly is defined as the time delay that maximizes the correlation function of an observed seismogram relative to a reference waveform. For assembling a global data set of waveform measurements, using all available broad-band recordings around the globe, we cannot rely on correlating seismograms to neighbouring seismograms because this latter method depends on correlating clusters of highly similar waveforms (e.g. VanDecar & Crosson 1990), a condition only met for local and regional studies. On the scale of global seismic networks, the earthquake radiation pattern strongly modifies a seismogram's relative P, pP and sP amplitudes as a function of azimuth and distance. In sparsely instrumented regions, the data of which are particularly valuable, even seismograms from neighbouring stations may not be sufficiently similar to cross-correlate. Instead, we had to adopt a method that correlates each broad-band seismogram with a corresponding synthetic waveform, developed by Sigloch & Nolet (2006) for teleseismic P waveforms of shallow events, modified by Hosseini & Sigloch (2015) for Pdiff waves, and further applied to PP waves for this study. Compared to data-to-data correlation methods, datato-synthetic correlation requires significant additional computation and data processing, but it means that we leave no data behind: we try to use every event exceeding magnitude m b ≥ 5.8, regardless of depth; every available broad-band station that might have recorded the event; and to measure finite-frequency traveltimes across the entire broad-band frequency range (Hosseini-zad et al. 2012;Hosseini & Sigloch 2015).
If an event was recorded by N stations (at suitable P, Pdiff, and PP distance ranges), then each of the N synthetic broad-band seismograms consists of a forward-computed Green's function convolved with a broad-band source time function (STF). Green's functions are computed by the fully numerical wave propagation method AxiSEM ; van Driel & Nissen-Meyer 2014a;van Driel et al. 2015). It solves the wave equation up to frequency 1 Hz through the spherically symmetric velocity model IASP91 (Kennett & Engdahl 1991). IASP91 serves as the reference model relative to which the 3-D tomography will be computed. For measurement purposes, the reference model also includes the intrinsic attenuation values of the PREM model (Dziewonski & Anderson 1981), so that P-wave traveltime dispersion due to intrinsic attenuation (which is small but appreciable) is included in the synthetics, and thus properly accounted for upon cross-correlation with the observed waveforms.
In a pre-processing stage, the broad-band STF wavelet is estimated (jointly deconvolved) from a suitable subset N sub of the N seismograms and their N Green's functions; this high-quality subset consists of stations at teleseismic distances that are part of the permanent, international networks (network codes II, IU, G, GE). Some STF estimates are shown in Fig. 1(a), together with their focal mechanisms (beachballs) in Fig. 1(b). The STFs shown are typical in that they have slight imperfections compared to the theoretical expectation of a non-negative, compact pulse that is identical to zero prior to t = 0 and after the rupture has terminated (t > 3 − 15 s depending on event). Moderate oscillations before and after these times reflect the effect of noise on the deconvolution, consisting of noise in the waveforms (ambient noise before rupture, plus signal-generated 'noise' after rupture) as well as slight misalignments of the N sub seismograms on the (unknown) origin time prior to deconvolution. STF shapes are used as quality diagnostics, that is, events are not used for subsequent traveltime measurements when their STF estimates look badly non-physical. See Sigloch & Nolet (2006) for details.
The vast majority of earthquakes occur at shallow depths (<30 km deep). In such a case, the depth phases pP and sP cannot be ignored or windowed out of the seismogram because they arrive within a few seconds of the P pulse, and often exceed the amplitude of P itself. A pulse spacing of seconds is neither short nor long compared to the period bands in which we measure finite-frequency traveltimes (3-30 s), or compared to the source time functions (rupture durations) of events that are sufficiently large to be measured globally. We attempt to include all events that range between 5.8 and 7.5 m b in magnitude, which have typical rupture durations of Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 1. Measurement procedure for multifrequency traveltimes on the example of 13 typical earthquakes of magnitudes 5.4 to 6.5 M w and two earthquakes of magnitudes 7.5 and 7.9 M w in the southeast Asia. For each event, its depth, moment tensor, and broad-band STF are deconvolved from teleseismic recordings around the globe using the method of Sigloch & Nolet (2006). (a) Estimates of the 15 broad-band STFs. (b) Estimated focal mechanisms of the 15 events plotted as beach balls. Panels (d)-(g) show observed and predicted waveforms for the event marked red in panels (a) and (b), a magnitude 7.5 earthquake in Southern Sumatra (2009/09/30 10:16:09.2, at 0.72 • S, 99.87 • E, 82 km depth). (c) Transfer functions of our eight Gabor filters with center periods ranging from 30.0 to 2.7 s. The frequency axis is logarithmic. (d) Observed (dashed) and synthetic (solid lines) waveforms for the red event filtered to 30.0 s dominant period, using filter '1' in (c). The stations are located in the Middle East and Europe, bracketed by orange dashed lines in (b). Blue waveforms indicate teleseismic P waves, black waveform are core-diffracted P waves (Pdiff), all amplitude-normalized and aligned relative to the theoretical P/Pdiff arrival times. (e) The same as (d) but filtered to 15.0 s period (filter '3'). Panels (f) and (g) are like (d) and (e), respectively, except that stations are located in the Far East and North America, bracketed by purple lines in (b). The combination of teleseismic P and Pdiff waves covers a wide epicentral distance range, from 34 • to 140 • in (f) and (g). The processed time window includes the direct P phase and its depth phases pP and sP, which strongly change the waveform as a function of azimuth and filter period. Pdiff waveforms look increasingly 'low-passed' with increasing distance, a characteristic of the core diffraction process (Hosseini & Sigloch 2015).

100
K.  s; the durations of the example STFs in Fig. 1(a) ranges from 3 to 15 s. When one of these STFs convolves the tightly sequenced P-pP-sP spike train in the broad-band Green's function of a shallow event, the three spikes are blurred together into an inseparable composite in the synthetic seismogram.
Next the N broad-band observed and N broad-band synthetic waveforms are filtered to the eight frequency passbands shown in Fig. 1(c). We use Gabor filters in overlapping frequency bands, which are Gaussians in the log-frequency domain, with dominant periods ranging between 30.0 and 2.7 s. They have constant fractional bandwidth and good spectral concentration properties (Sigloch & Nolet 2006). In principle, we could measure up to 1 s dominant period (the corner period of the Green's functions) but the percentage of good correlations becomes very low beyond the <2.7 s-band, according to the success statistics for P and Pdiff of Hosseini & Sigloch (2015).
N cross-correlation functions are computed in each of the eight passbands. We obtain 8N finite-frequency traveltime anomalies, each of which is the time lag that maximizes a cross-correlation function; the corresponding maximum correlation value (crosscorrelation coefficient CC) serves as an important quality diagnostic. Panels (d)-(g) in Fig. 1 compare, for the red earthquake, its observed and synthetic P and Pdiff waveforms in different passbands after time-shifting each synthetic by the finite-frequency traveltime anomaly. Arriving 20-40 s after P, the pP and sP pulses are clearly separated in time only because this was a relatively deep event (82 km). Note the systematic changes of relative pulse amplitudes within each cluster; and the increasingly smoother appearance of the waveforms at large P-diffracted distances (loss of high-frequency content through diffraction). For all these reasons, this set of waveforms is too dissimilar for a data-to-data correlation strategy to be appropriate.
As misfit measure for tomography, cross-correlation traveltimes have very favorable properties. They are robust in that crosscorrelation or matched filtering is the optimal way of detecting a signal of known shape in white Gaussian noise (Sigloch & Nolet 2006, and signal processing references therein). This robustness is important in environments of high ambient noise, such as Pdiff waves at long distance ranges and high frequencies (Hosseini & Sigloch 2015). Cross-correlation traveltimes also retain a near-linear relationship to velocity model perturbations over a larger perturbation range than direct waveform differences (i.e. a samplewise L2 misfit) (Mercerat & Nolet 2013).
We estimated source time functions and applied the multiplefrequency (MF) measurement procedure to a global data set of earthquake sources. Source time function estimates were attempted for all earthquakes that occurred between 1999 and 2015 and exceeded magnitude m b ≥ 5.8. Many earlier events are also included (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999). Smaller events were attempted if they occurred in uncommon locations. Fig. 2(a) shows the resulting ≈3500 sources for which STFs and MF measurements were successfully obtained. Fig. 2(c) shows the 5389 broad-band stations that yielded successful MF measurements.
The large volumes of broad-band waveform data were retrieved from data centers that support FDSN web services (e.g. IRIS and ORFEUS data management centres) using the purpose-built, fully automated Python software obspyDMT (Hosseini & Sigloch 2017), which also executed instrument corrections to ground displacement, bandpass-filtering between 0.01 and 3.5 Hz, archiving on our local computer network, and automated updating of the data holdings when new waveforms became available on the servers of the data centres. Our obspyDMT software supports many other routine data retrieval, management and plotting tasks, on small and large sets of event-based and time-continuous seismograms. It is freely available at http://kasra-hosseini.github.io/obspyDMT [last accessed: August 2019] and described in Hosseini & Sigloch (2017).
About 13.8 million teleseismic P traveltimes were measured in eight frequency bands, as well as 32.5 million PP and 16.3 million Pdiff traveltimes. In order to reduce data redundancy and limit the computational cost of calculating sensitivity kernels, we used P measurements from only four frequency bands, PP measurements from only two bands, and Pdiff measurements from all eight bands. We only accepted measurements above the quality threshold of CC ≥ 0.8 for P and Pdiff, and CC ≥ 0.85 for PP, independent of frequency band. The standard deviation of the measurement error was assigned to each measurement according to its CC, ranging from ≈0.3 s for high-quality measurements to ≈1 s for low CC. The standard deviations were used to populate the data covariance matrix C d (eq. 2 in Section 3.3). Table 1 lists the number of accepted MF measurements on P, PP and Pdiff waveforms, which were subsequently used to generate tomography models DETOX-P1, -P2 and -P3. Fig. 3 shows histograms for the P, PP and Pdiff measurements of Table 1, after correcting the raw traveltime measurements for earth's ellipticity, topography and crustal structure (using crustal model CRUST2.0; Bassin et al. 2000). The corrections are described in more detail in Appendix A. Fig. 3(c) plots histograms of the corrected dT anomalies versus epicentral distance. Red dashed lines show the ray-theoretically predicted arrival times for an event at 0 km depth; the bulk of arrivals occur earlier because the sources have finite depths. Mantle heterogeneities (Oldham 1906) and uncertainties in epicentre location also contribute to deviations from the red line.
We exclude teleseismic P waves at distances >85 • where the core-reflected PcP phase arrives shortly after the direct P wave, thus interfering with the measurement of the teleseismic P-wave traveltime. At ranges moderately shorter than 85 • , this interference is too small to warrant exclusion (Sigloch & Nolet 2006). Although some global tomography models have included traveltimes of the corereflected P wave (PcP) (e.g. Amaru 2007;Simmons et al. 2012), or equivalently the ScS phase (e.g. Wysession et al. 1994), we explicitly decided to exclude PcP measurements. A secondary reason was the relatively low signal-to-noise ratio of PcP, due to small CMB reflection coefficients, especially at distances <70 • (Zhu & Wysession 1997). More importantly, Colombi et al. (2012) used 3-D wave-propagation simulations and sensitivity analysis to show that PcP waves are very sensitive to topography in the CMB (as opposed to volumetric heterogeneity near the CMB), whereas Pdiff measurements are less sensitive to undulations of the CMB, and most sensitive to volumetric structure. Given that CMB topography is unknown on the length scales relevant to us, we would not be able to correct for it. Instead parametrizing and inverting for CMB topography and volumetric structure jointly from PcP and Pdiff is a challenge left for the future. Here, we parametrize for volumetric structure only (like prior global body-wave models) and reject the PcP phase that is so sensitive to non-volumetric structure.
Core-refracted P waves (PKP) have been included in some global tomographies (e.g., Kàrason & Van der Hilst 2001;Amaru 2007;Li et al. 2008) and would have been a good addition here but have been left for the future. Both PKP and PcP waves tend to add high resolution in certain, regionally confined patches of lowermost mantle (Zhu & Wysession 1997), compared to the spatially more even and comprehensive coverage of Pdiff waves (Hosseini & Sigloch 2015).  Colour indicates average traveltime anomalies, averaged over all measurements regardless of event location, seismic phase and frequency, and before attempting to re-estimate origin times. (d) 6051 unique receiver locations selected from the ISC-EHB/EHB catalogues. This relatively small subset of all available ISC-EHB/EHB picks was selected to provide the most complementary coverage to MF data. (e) and (f) show the relative contributions of each data set. Earth's surface was subdivided into blocks of 4 • . The number of successful measurements was counted in each block for MF and ISC-EHB/EHB measurements; the total of the two numbers adds up to 100 per cent for each block (shades of red) that contains stations, else the blocks are white. The complementarity of the two data sets is clearly evident. MF measurements in (e) provide most of the information from the western Indian Ocean , East Africa (Nyblade et al. 2008(Nyblade et al. , 2011, Hawaii (Wolfe et al. 2009), French Polynesia (Barruol et al. 2002Suetsugu et al. 2005), offshore Cascadia (Toomey et al. 2014), Greenland and Antarctica. ISC-EHB/EHB stations in (f) greatly enhance the MF coverage and represent the main source of information for several regions, such as Brazil, India, Iran and Saudi Arabia. and Pdiff. We attempt to measure the teleseismic P phase in the distance range 32-85 • , PP in range 64-155 • and Pdiff from its emergence around 98 • to the maximum range with a clear waveform, typically to 160 • . (b) Histograms of multifrequency traveltime anomalies, after correction for earth's ellipticity, crustal structure, and topography. See Table 1 for the absolute numbers. (c) Distribution of successful MF traveltimes as a function of epicentral distance. Colour indicates the number (density) of measurements. Red dashed lines show the ray-theoretically predicted arrival times for an event at 0 km depth in the IASP91 reference model. We exclude teleseismic P measurements ranging between 85 • and the emergence of the Pdiff phase because the core-reflected PcP waveform interferes significantly with the direct P-wave at these distances.

Analyst-picked arrival times from the ISC
In order to achieve the best possible data coverage, we complemented our global multifrequency traveltimes by a subset of analystpicked P arrival-times from the EHB ) and ISC-EHB (Weston et al. 2018) catalogues. Multifrequency measurements depend on the availability of high-quality broad-band waveforms, which date from the 1990s onward. The International Seismological Centre (ISC) has been aggregating analyst-picked traveltimes that date back to the 1960s and older, from many geographical sites that have not subsequently been reoccupied by broadband instruments. Even today the ISC receives picks from many stations (broadband or other) that do not send their waveforms to international data centres.
EHB traveltimes are obtained from picked arrival times reported to the ISC and to the National Earthquake Information Center (NEIC), by reprocessing them using the Engdahl et al. (1998) algorithm, which includes iterative, non-linear relocation of earthquake sources and dynamic phase re-identification. The ISC recently released a new catalogue called ISC-EHB (Weston et al. 2018), which replaces the original EHB catalogue for events from 2000 to 2008 and in addition contains teleseismically well-constrained events from 2000 to 2014. We selected arrival times for 2000-2014 from the ISC-EHB catalogue, and for 1960-1999 from the EHB catalogue. They were corrected for the effects of crust, topography and ellipticity as detailed in Appendix A.
The ISC-EHB bulletin contains data for ≈140 000 events and many millions of measurements for various seismic phases. Our aim was to select the smallest useful subset of teleseismic P picks that would fill as many gaps as possible in our MF coverage (which are evident from the highly uneven source and station distributions of Figs 2a, c and e).
We considered P arrival picks at epicentral distances 30 • -90 • if the precision of their arrival time reading was better (smaller) than 0.1 s. A method was developed to homogenize the directional density coverage of measurements at source-receiver locations. The earth's surface was divided into blocks of 2 • . For each block that enclosed at least one source or one receiver, the number of rays and their directions were stored. When the number of rays from one of the four azimuthal bins (0-90 • , 90-180 • , 180-270 • and 270-360 • ) reached a threshold, no more rays would be accepted for that direction and block. Good azimuthal coverage is also desirable for each event because source (re)location is attempted during tomography. Hence, stations at a range of azimuths were selected for each event. The final selection included 5.7 million teleseismic P-phase picks, see Table 1.
Figs 2(b) and (d) show the distribution of sources and receivers for these ISC data. Comparison of Figs 2(e) and (f) reveals the benefit of adding EHB/ISC-EHB traveltime picks. They contribute the bulk of information for a number of regions, including Brazil, India, Iran and Saudi Arabia. Certain other regions, especially in the oceans, are constrained only by MF waveform data, often obtained from recent ocean-bottom seismometer deployments (Barruol et al. 2002;Suetsugu et al. 2005;Wolfe et al. 2009;Toomey et al. 2014;Suetsugu & Shiobara 2014;Stähler et al. 2016).

M E T H O D
We assume that the traveltime observations d i (i = 1, ..., N) are linearly linked to the discretized earth model m j (j = 1, ..., M) in the following explicit form (see Appendix B for the detailed description): (1) where m contains the P-velocity deviations dV P /V P from the reference earth model IASP91 on a suitably chosen inversion grid, and A links the traveltime observations in d to the discretized model parameters in m.

Sensitivity kernels
Each row in A of eq. (1) contains the 3-D sensitivity kernel K T of one measurement numerically integrated onto the inversion grid (see Fig. 4, right-hand column). The use of sensitivity kernels allows us to properly distribute the sensitivity of a given seismic phase with a specific frequency content into the mantle. As we have two types of traveltime observations (Section 2), two groups of sensitivity kernels are also calculated: frequency-dependent and ray-based kernels. Dahlen et al. (2000) developed a method to efficiently calculate sensitivity kernels K T using paraxial approximation in conjunction with ray theory and the Born approximation (Born 1926). Using this method as implemented by Tian et al. (2007b), we computed a 3-D frequency-dependent sensitivity kernel for each cross-correlation traveltime observation in our data set (i.e. ≈ 5M sensitivity kernels). Fig. 4 shows P, PP and Pdiff kernels.
Due to the inherent limitations of ray-theory, Pdiff kernels cannot be accurately calculated using the method of Dahlen et al. (2000). Here, approximated Pdiff kernels were computed instead, as shown in Fig. 4. Computation of approximated Pdiff kernels is described in more detail in Appendix C. We benchmarked these Pdiff sensitivity kernels against wavefield-based kernels calculated by AxiSEM ) and MC Kernel  for various epicentral distances, event depths and relevant frequencies. The approximated kernels give a good representation of the 'true' core-diffracted sensitivities. This approximation decreased the computational cost of one sensitivity kernel by at least three to four orders-of-magnitude compared to wavefield-based kernel calculations, that is, 3 s kernel -1 versus 18 000 s kernel -1 (a wavefield with 0.1 Hz highest frequency was used for the latter). We will incorporate wavefield-based kernels once work on the necessary software package 'MC Kernel' is complete ) and optimized for large data sets.
We use ray theory (hypothetically infinite frequency) to efficiently compute the sensitivity kernels of onset times selected from ISC-EHB/EHB bulletins. Rays are numerically integrated on the inversion grid as shown in Fig. 4 (bottom row) for P.

Model parametrization
We sample the unknown velocity structure of the mantle using an irregular distribution of points connected into space-filling tetrahedrons, with linear interpolation of dV P /V P between the four vertices of each tetrahedron. The spatial distribution of model parameters limits the maximum resolution retrievable in different regions. Therefore, the interpolation supports should be sufficiently dense to not limit the expression of the data's full information content, anywhere in the mantle. On the other hand, very dense parametrization . 3-D sensitivity kernels for cross-correlation traveltimes of a teleseismic P wave at 60 • distance, of a PP wave at 120 • , and a Pdiff wave at 120 • . The bottom row shows a ray-theoretical sensitivity kernel used to model an analyst-picked, teleseismic P arrival time from the ISC-EHB/EHB catalogue and at 60 • distance. The ray-theoretical kernels produce very sparse matrix rows, with only 0.05 per cent non-zero entries on average. In case of cross-correlation sensitivity kernels (the first three rows), the dominant period of the Gabor filter is 21.2 s. Left-hand column shows strictly 2-D sections through the mantle in the plane of wave propagation; right-hand column shows the same kernels numerically integrated on the tetrahedral inversion grid. P and PP kernels are computed by the method of Dahlen et al. (2000). PP kernels were limited to distances <155 • because the paraxial approximation is not valid near antipodal epicentral distances (Tian et al. 2007b). Since diffracted wave kernels cannot be calculated with the paraxial method, Pdiff kernels are approximated (see text). Each kernel fills one sparse row of the inversion matrix A. The fraction of non-zero entries (i.e. tetrahedral grid nodes of non-zero sensitivity) is on average 1.9 per cent for P kernels, 4.1 per cent for PP kernels and 2.9 per cent for Pdiff kernels.
increases the computational cost for both kernel calculations and inversion. Irregular tetrahedral meshes have the advantage of being able to accommodate different sampling densities in a straightforward manner (Sambridge et al. 1995;Sambridge & Gudmundsson 1998;Nolet & Montelli 2005). The mesh can thus be rendered denser selectively in places where the data offers increased resolving power, for example, densely instrumented regions.
We use the DistMesh tool (Persson & Strang 2004) and the Quickhull method of Barber et al. (1996) to generate non-overlapping tetrahedrons based on Delaunay criteria. DistMesh generates tetrahedral meshes of high quality as measured by the regularity of the tetrahedrons (large ratios of volume to average facet length cubed). To ensure that the hull nodes always enclose the earth's surface, the radius of the convex hull is set to r + 50 km where r = 6371 km is the earth's radius.
A method was developed to generate spherical tetrahedral meshes with adaptive cell sizing based on the kernel coverage with the following steps (Hosseini 2016): first, a (quasi-)homogeneous inversion grid with 200 km edge-length tetrahedrons was generated, and the sensitivity kernels of a relatively large, representative subset of the data were calculated. The absolute values of the kernels were then projected and summed at each node. This is similar to a column sum plot of A, and we assume that it is a good measure of resolving power of the underlying data. In the second step, a new inversion grid was generated in which the number of interpolation supports was increased for the regions with large sensitivity kernel values (≈80 km edge-length in the best-sampled regions). This was achieved by assigning the desired edge-length to each node of the grid according to its sensitivity kernel value. After some iterations, the final inversion grid with 396 501 nodes was generated as shown in Fig. 5. DETOX-P1, P2 and P3 models use the same parametrization. Fig. 6 compares two versions of the same P-wave kernel (60 • epicentral distance), once discretized onto a finely meshed part of the inversion grid under the northern Pacific rim (left-hand column), and once projected onto a coarsely meshed region beneath the southernmost Pacific (right-hand column). The primary features of the kernel are rendered even on the coarse grid.

Model space (m vector)
Due to rank deficiency in A of eq. (1), an infinite set of equally satisfactory solutions exits for m (Nolet 2008;Menke 2018). We use a combination of smoothing and damping operators to regularize the problem: n is a tunable parameter to control the strength of norm damping. c −1/2 m acts as smoothing operator and s is a smoothing factor that can be used to control the intensity of the smoothing operator. The data covariance matrix C d is assumed to be diagonal (i.e. uncorrelated errors) with assigned a priori uncertainties σ 2 i as entries (Section 2).
Entries of m in eq. (2) are the 3-D solutions for isotropic compressional P-wave velocity structure dV p /V p with additional unknowns to account for hypocentral parameters (origin time, longitude, latitude and depth). A priori uncertainties were assigned to each component as follow: 1 per cent for dV p /V p , 10 km for latitude, longitude and depth of hypocentre corrections and 2.0 s for event origin time corrections. Fig. 7 schematically shows the system of equations to solve: d r = Gm in which d r and G are d and A in eq. (2) with appended rows for regularization.
In DETOX-P3, A contains ≈10.7M data rows and ≈606K columns with ≈ 54 721M non-zero elements. Since every kernel is close to zero in most areas of the mantle, the system of equations is sparse (0.8 per cent in total), that is, ≈ 54 721M non-zero elements in A out of ≈ 6 510 212M elements. On average, the fractions of non-zeros for multifrequency measurement rows were ≈1.9 per cent for P, ≈4.1 per cent for PP, and ≈ 2.9 per cent for Pdiff. This for picked arrival times was ≈0.05 per cent.
Many of the rows in A of eq. (2) are dependent (N M), and the measurements in d have errors. As a result, the system of equations cannot be solved explicitly, and instead, we minimize the L 2 norm ||Gm − d r || 2 2 with LSQR algorithm, a conjugate-gradient type method for solving large, sparse, linear equations (Paige & Saunders 1982b;Nolet 1985). Similar to conjugate gradient least squares (CGLS) method, row-based LSQR uses a recursive scheme that monotonically reduces the residual vector (with no actual matrix inversion). In theory and in the absence of machine accuracy problems, LSQR is algebraically equivalent to applying CG to the normal equations G T Gm = G T d r but with better numerical properties in practice (Paige & Saunders 1982a,b). As the model is parametrized on an irregular tetrahedral mesh, the approximate volume associated with each vertex is computed and used to weight model parameters in the LSQR inversions (Zaroli et al. 2015;Aster et al. 2018). For this, we assume that each node is enclosed in a hypothetical tetrahedron whose facet length equals to the average length of all facets adjoining the node.
We employ a two-step approach to simultaneously invert for structure and for effects of source mislocation and origin time errors (Nolet 2008). First, an inversion with very low regularization parameters was performed to identify outliers. The measurements that deviated more than 3σ from the general distribution were labeled as outliers and removed from the subsequent inversions. Secondly, inversions with variable regularization parameters were done. To ensure convergence, the results were obtained after 10 000 LSQR iterations; however, with LSQR most of the convergence is achieved within a small number of iterations (500-1000 iterations, except for the inversions with very low regularization parameters).

Selection of preferred models
By changing the intensity of regularization parameters ( n and s ) in eq. (2), the trade-off between the residual data misfit in terms of reduced chi-squared χ 2 red (absolute chi-squared divided by the number of data N): and the L 2 model norm ||m|| 2 2 as a measure of model complexity can be analyzed. Fig. 8 shows an example of this trade-off (Lcurve) for DETOX-P2; two cross sections at 1300 and 2800 km depths are plotted through the tomography models obtained by constant ratio damping smoothing = 0.25 but with different amplitudes of damping and smoothing parameters. By increasing the intensity of regularization, models show lower norms ||m|| 2 2 and higher χ 2 red , as expected.
Selection of smoothing and damping parameters is normally done by assuming that the models close to the bending of the L-curve are the best compromise between minimizing both χ 2 red and ||m|| 2 2 . Models far from the bending are either dominated by noise ('vertical' part of the curve with high ||m|| 2 2 values, e.g. χ 2 red = 0.93 in Fig. 8) or the data is poorly exploited in them ('horizontal' part Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 5. Global tetrahedral grid used for inversion. Adaptive gridding improves the poor numerical conditioning of the inverse problem, which is caused by extremely uneven source and receiver distributions. Irregular tetrahedra are easily adapted to the expected spatial resolution of the tomographic model, as estimated a priori by comparing the cumulative kernel sensitivities of grid nodes (a proxy for 'goodness of sampling'). Right panel shows average distance between one vertex and all its neighbours as a function of vertex depth. Node spacing or tetrahedral facet length varies between ≈80 km in the upper mantle under highly sampled regions (e.g. North America) and ≈400 km in the core. The volume of a typical tetrahedron is 5-7 times smaller than the volume of a cube of equal facet length. 396 501 nodes form 2 338 710 tetrahedra, with grid connectivity statistics summarized by the histogram.
with high χ 2 red and low ||m|| 2 2 values, e.g. χ 2 red = 2.0 in Fig. 8). See Zaroli et al. (2013) for detailed discussion. However, there is a subjectivity inherent to the choice of models close to the bending of the L-curve due to imperfect quantification of data uncertainties, especially their correlations.
We experimented extensively with the damping smoothing ratios using the whole or only a subset of data vector, and a set of models close to the bending of L-curves was chosen as preferred solutions to the inverse problem and will be discussed in Sections 4, 5 and 6. Changes in the model amplitudes due to regularization can affect the interpretation of tomography models. As an example, conversion of velocity structure to temperature field, one of the inputs for geodynamic modellings, is directly influenced by the amplitudes of anomalies. In Appendix E, the effects of regularization on the amplitudes of our tomography models are discussed.

A P P R A I S A L O F G L O B A L P -T O M O G R A P H Y M O D E L S D E T O X -P 1 , P 2 A N D P 3
This section contains technical appraisals of our global tomography models generated from teleseismic P traveltimes (DETOX-P1), P and Pdiff (DETOX-P2), or P, Pdiff and PP (DETOX-P3). Subsequent sections will discuss results in terms of earth structure.
All three DETOX models can be accessed and interactively visualized via the SubMachine web portal (http://submachine.earth.ox. ac.uk/, last accessed: Aug. 2019). We invite the reader to visualize Figure 6. The effect of numerically integrating the quasi-continuous sensitivity kernels onto the relatively coarse inversion mesh, which is also spatially variable. We compare two versions of the same P-wave kernel (60 • epicentral distance), once discretized onto a finely meshed part of the inversion grid (northern Pacific rim, left-hand column), and once projected onto a coarsely meshed region (southernmost Pacific, right-hand column). Top panels show the density (average spacing) of mesh vertices at 100 km depth, with the two wavepaths superimposed as black lines with coloured dots. Middle two panels show mesh vertex densities in the whole-mantle cross-sections defined by the two wavepaths. Bottom panels show the two discretized kernels along the same cross-sections. The primary features of the kernel, that is, its 'doughnut hole' (zero sensitivity along the ray path) and its second Fresnel zone of opposite polarity (blue), are rendered even on the coarse grid. The continuous kernels are symmetric along the central ray path but are rendered slightly asymmetrically on the irregular mesh. Kernel amplitudes appear higher on the coarse grid (bottom right) than on the fine grid (bottom left) because the continuous kernel's total sensitivity, which is identical in both cases, is mapped onto fewer vertices in the case of a coarse grid region. See supplementary Fig. D1 in Appendix D for more mesh density maps. our models themselves as they digest the remainder of this study. SubMachine is a collection of web-based tools for visualization, analysis and quantitative comparison of global-scale, volumetric (3-D) data sets of the subsurface . It also provides functionality to overlay reconstructed plate boundaries, coastlines, Ultra-Low Velocity Zones, hotspot locations, and/or shear wave splitting.

Sampling improvements from adding PP & Pdiff measurements
We define the kernel coverage for a given node of the inversion grid as the sum of absolute values of all sensitivity kernels that sensed that node, that is i |A ij | for the jth node. Thus cumulative sensitivity kernels in an area (as compared to other areas) can be thought of as fairly representative of the ability to resolve structural features, even though this measure ignores the extent to which kernels crisscross, and at which angles (both are important for good resolution). Fig. 9 shows kernel coverage at several depths and for three data sets: only teleseismic P measurements as in DETOX-P1 (left-hand column); P and Pdiff as in DETOX-P2 (middle); and P, Pdiff and PP measurements as used in DETOX-P3 (right-hand column). Before computing the column sums, rows of A were scaled according to a priori uncertainties assigned to their corresponding measurements [the data covariance matrix C d in eq. (2), Section 3.3].
Body waves sample the upper mantle only in the vicinity of sources and receivers, which accounts for their very heterogeneous cumulative sensitivities at 100 and 600 km depth. For DETOX-P3, the bouncing points of PP waves significantly enhance kernel coverage at these depths because they also occur in regions without earthquakes and stations, which includes the oceans and high latitudes. Coverage becomes more homogeneous in the mid Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 7. Schematic of the linear system of equations d r = Gm to be solved by regularized least-squares minimization min||Gm − d r || 2 2 . Matrix G consists of kernel matrix A (one kernel per row) and 2 × M appended rows for norm damping and smoothing (in yellow). Vector d r is populated by traveltime anomalies d (N = 10 740 922 elements in DETOX-P3; 9 570 176 elements in DETOX-P2; 8 190 884 elements in DETOX-P1) and 2 × M appended zeros corresponding to regularization rows. Number of successful records of each data type/phase used in the actual inversion (after removing the outliers) differs only marginally between DETOX-P1, P2 and P3 models (see Table 1). Model vector m consists of P-velocity anomalies dV P /V P (M = 396 501 entries) and appended elements for source corrections c (in purple, four elements per event: hypocentre and source origin time). mantle (1200 km) for all three data sets. These depths are constrained mainly by teleseismic P waves. The southern hemisphere remains undersampled.
The cumulative sensitivity of DETOX-P1 starts to significantly fall behind that of DETOX-P2 and -P3 in the lower third of the mantle (2000-2800 km), where Pdiff waves contribute massively. Inclusion of this thus neglected wave type has increased cumulative sensitivity by factors of 100-1000, compared to only teleseismic P waves. (Section 7 will return to the question how the lowermost mantle has been illuminated (at all) by existing global P models, in which core-sensitive waves have not been systematically used.)

Resolution tests
For resolution tests, synthetic data vectors d s yn were constructed based on synthetic input models m inp : Solutions were then obtained through the same inversion procedure and parameter settings as for the observed data. Comparisons between input and output models permit an appraisal of resolution for our global models at different depths. Fig. 10 compares the resolution of DETOX-P1 and DETOX-P2 using as test functions 3-D Gaussians dV P /V P = (dV P /V P ) center exp ( − r 2 /w 2 ) with radii w = 800 km and (dV P /V P ) center = ±3 per cent. r is distance from the center of each 3-D Gaussian velocity anomaly.
In the upper mantle, the very variable spatial resolution is controlled by the uneven distribution of sources and receivers. Wellresolved areas mainly coincide with instrumented continents and with regions of high seismicity, as well as PP bouncing points in the case of DETOX-P3 (see Fig. F1 in Appendix F for a resolution test comparing DETOX-P2 and DETOX-P3). The amplitudes of anomalies in the upper mantle are decreased significantly by regularization, especially in sparsely instrumented regions. In the mid and lower mantle, input patterns are well resolved in both inversions. In the lowermost mantle, input patterns are retrieved very well in DETOX-P2 due to the inclusion of Pdiff measurements. Fig. 11 compares the DETOX models side-by-side at five different depths. At shallow depths (100 and 600 km), DETOX models show similar structures especially in well instrument regions, e.g., North America, Europe, Southeast/East Asia, Australia and East Africa. The most pronounced slow anomalies are imaged in the west of North and South America, East African Rift, western and central Europe, mid-Atlantic ridge and Tibetan plateau. More localized low-velocity regions beneath several hotspot locations (e.g. Iceland, Hawaii and the Samoa-Tahiti-Society Islands group) are apparent in all models as well. At 100 km depth, DETOX models reveal highvelocity features beneath various subduction zones and cratonic regions (e.g. western Australia, eastern Europe, Congo and Siberia). As expected, DETOX-P1 and P2 show patchy anomalies in the oceans due to the lack of earthquakes and seismic instruments. In fact, the only data sources in oceans are sparse island stations and ocean bottom seismometers (see Fig. 2 for the source-receiver distribution, and Fig. 9 for the global kernel coverage). DETOX-P3, on the other hand, images more structures in the oceans, for example, low-velocities beneath mid-Atlantic ridge and North-tocentral Pacific, and high-velocities in North Asia. Unfortunately, part of the low-velocity structures in oceans are artefacts that arise due to downward smearing of lithospheric structure; see Section 4.4 for discussion.

Side-by-side comparison of the DETOX models
From the upper mantle transition zone to mid-mantle depths, long and narrow high-velocity structures appear beneath North and South America, southern Eurasia, New Zealand, Alaska and Japan in all models. At 1200 km depth, there is a strong agreement between the DETOX models on these sharp high-velocity features. The main difference between the models is on low-velocities imaged in North-East Pacific Ocean in DETOX-P3 model.
At 2200 and 2800 km depth profiles of Fig. 11, DETOX-P1 either shows no or very low-amplitude features. This model contains only teleseismic P waves, and although deep-diving P waves in the far teleseismic range start to sense the CMB, the majority of them sample the mantle at depths ≤1800 km. The sharp decline in the global data coverage at these depths (see Fig. 9) results in very low resolution models in the lower-third of the mantle. Nevertheless, high velocities beneath Central America and Eastern Asia are detectable. On the other hand, thanks to core-sensitive Pdiff measurements, various high and low velocity features are imaged in DETOX-P2 and -P3 at depths ≥2000 km which will be discussed in Sections 5 and 6. Fig. 12 shows decompositions of the three DETOX models into spherical harmonic coefficients using SHTools (Wieczorek & Meschede 2018). The comparison shows how the addition of Pdiff data (going from DETOX-P1 to DETOX-P2) greatly intensifies dV P /V P structure imaged on the lower half of the mantle. Degree 2 structure is dominant only in the lowermost 400 km, and directly above the CMB, degree 3 is equally strong as degree 2. Addition of PP data in DETOX-P3 intensifies dV P /V P structure in the upper half of the mantle, including at very long wave lengths, which points Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 8. L-curves quantify the trade-off between model norm ||m|| 2 2 (a proxy for model complexity) and data misfit χ 2 red as the severity of regularization is varied. The effects on the solution of increasing regularization (from left to right along the L-curve) are compared visually for DETOX-P2, at 1300 km (top) and 2800 km depth (bottom). The ratio damping smoothing = 0.25 is the same for all solutions. Regions of nominally best data fit (very low χ 2 red values) are undesirable because those solutions also attempt to fit the measurement noise and become very detailed and 'grainy' to do so. The preferred DETOX-P2 model (framed in black) is more conservative and "simple" in appearance, its moderately lower data fit of χ 2 red = 1.23 indicating that noise is being rejected. Very strongly regularized solutions produce overly smooth/damped models by overriding valid information contained in the data (high misfits χ 2 red ). (If measurement uncertainties σ i were perfectly known, we would simply accept the theoretically optimal mode χ 2 red = 1 (dashed red line) as the preferred solution. Our subjective preference for χ 2 red = 1.23 implies that the σ i were underestimated slightly.).
Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 9. Comparison of how models DETOX-P1, P2 and P3 sample the mantle at different depths. Plotted in gray shades are cumulative sensitivity kernel values (sum over columns of kernel matrix A: i |A ij | for the jth node). Dark shades indicate that many kernels traverse the region, which is a moderately good proxy for tomographic resolution. (A better proxy would take into account the incidence angles of kernels.) Note that the colour scale is logarithmic. Compared to DETOX-P1, which contains only teleseismic P waves, DETOX-P2 and P3 show much better sampling below 2000 km, thanks to Pdiff waves. DETOX-P3 also includes PP waves, which mainly improve sampling under the oceans in the upper half of the mantle.
to large areas of upper mantle beneath the oceans featuring slower structure than indicated by the reference model, as discussed in Section 4.4.

The effects of including PP waveform measurements
The difference between DETOX-P2 and DETOX-P3 is that the latter includes PP traveltimes (cf. Fig. 3  Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 11. Comparison of tomography models DETOX-P1 (using P waves), DETOX-P2 (P and Pdiff) and DETOX-P3 (P, Pdiff, PP) at different depths. The inversion grid and regularization procedure are the same for all models, and the resulting data fits are almost the same (χ 2 red = 1.15, 1.23, 1.20 for DETOX-P1, -P2 and -P3, respectively). DETOX-P2 and P3, which include Pdiff waves, produce much more structure in the lowermost mantle than DETOX-P1, consistent with the resolution tests of Figs 10 and F1. DETOX-P3, the only model to contain PP-waves, produces significant additional structure in the upper 1500 km of the Pacific and Indian Oceans (see Section 4.4). (Fig. 12). This is expected because the addition of data constraints in regions sampled only by PP (large parts of the oceans, especially the Pacific, see Fig. 9 at 100, 600 and 1200 km depth) pushes back against the tendency of norm damping to render these regions almost structureless in DETOX-P2. However, the most pronounced intensification of structure appears at the longest spatial scales (harmonic degrees 0-4, see Fig. G2), which was unexpected. At degree zero, which indicates static offsets from reference model IASP91, DETOX-P3 differs from DETOX-P2 mainly between 500-900 km depth. In the map renderings of Fig. 11, the difference in the longest wavelengths is expressed in a pattern that, at 600 and 1200 km depth, vast tracts of Pacific Ocean, and Indian Ocean to some extent, appear 'redder' (seismically slower) in DETOX-P3 than in DETOX-P2. These slow anomalies are robust under regularization.
Do these additional, large-scale velocity anomalies in DETOX-P3, which are always slower-than-average and located beneath the oceans (concentrated between 500 and 900 km depth), correspond to actual velocity anomalies in the mantle, or might they represent Figure 12. Spectral decomposition of the three DETOX models into spherical harmonic functions, as a function of depth. Colour shades show the magnitudes of the spherical harmonic coefficients as a function of harmonic degree on the (logarithmic) x-axis, and of depth in the mantle on the y-axis. The logarithmic colour scale is chosen such that every colour increment signifies and increases in spectral power of a factor of √ 2. The physical meaning of intensifying colour shades for DETOX-P2 and -P3 is the same as in the maps of Figs 9, 10 and 11: the magnitude of dV P /V P increases as more data are adding constraints on mantle structure that was previously pushed towards zero (the a priori model) by the action of norm damping. In supplementary Figs G1 and G2, we compare these DETOX power spectra to those of other global tomography models. an artefact, software bug, or waveform measurement problem? True structure is plausible, given that these mantle volumes are sampled only by PP waves. The reference model IASP91 (Kennett & Engdahl 1991), relative to which DETOX-P3 appears slow under the oceans, was constructed from P-wave picks of major phases reported to the ISC between 1964 and 1987. These were overwhelmingly teleseismic P picks, with PP quite noisy and not very numerous. The limited range of teleseismic P means that they sample mantle structure mostly beneath continents, a spatial bias shared by IASP91. If the upper 1500 km of mantle beneath oceans and continents featured systematically differing P-wave velocities, this might be transpiring only now that PP measurements are becoming abundant, especially high-quality waveform measurements. We think that truly slower Pvelocity structure beneath oceans is the most likely explanation, and that therefore DETOX-P3 is a more complete and 'better' mantle model than DETOX-P2, as it should be.
Independent checks of this assertion are not straightforward. Other global tomography studies do not seem to have reported the difference between excluding versus including PP waves. To get indirect hints, we attempted to check the harmonic degree-zero structure of different global models that included PP or SS waves, and whether more of them might indicate slower-than-reference velocities at 500-900 km depth. There are only two other global Pmodels using finite-frequency data: Montelli et al. (2004b) agrees that the mantle between 500 and 900 km is slower than indicated by IASP91 (degree 0 structure in Fig. G2), which is encouraging. Obayashi et al. (2013) reported their 3-D model relative to a custombuilt reference model, which we could not easily compare. The same is true for one of only two finite-frequency S-wave models (Zaroli et al. 2015). The other S-wave model, by Montelli et al. (2006), is slower than IASP91 at 900-1500 km, but the P-and S-velocities of IASP91 may well have different biases. Other global tomographies are even less comparable to ours in terms of data used. At this stage, no clear conclusions can be drawn from these comparisons.
When we first noticed the large-scale, high-amplitude, slowvelocity differences between DETOX-P2 and -P3, we suspected issues with crustal and or topography corrections at the PP bounce points, which would smear downward into the mantle along sparse PP paths. In the oceans, the differences between reality (thin oceanic crust plus water layer) and the reference model (35 km of continental crust, no water layer) are very large (Fig. A3), and they are sampled twice by the up-down reflecting PP wave. Fig. A2 shows the combined crust-plus-topography correction times of all our PP measurements, plotted at their bounce points. For oceanic bounce points, the reference model typically predicts PP arrival times that are delayed by 2-3 s relative to the actually measured arrival times. A mistake in implementing the crustal and topography corrections at the bounce points could therefore translate into strong artefacts that would resemble the slow, oceanic anomalies present in DETOX-P3, because (wrong) crustal structure would be smeared downward along the steeply incident PP bouncing legs.
We have checked and re-checked our implementation of the correction computations and found no mistakes. It also seems implausible that the values of the crustal and bathymetry corrections indicated by model CRUST2.0 could be wrong enough to produce artefacts of this magnitude, given that oceanic structure is relatively simple and predictable. To get a sense for the magnitude of spurious mantle artefacts produced by a worst-case mistake, we designed a test, reported in detail in Appendix A2 and Fig. A4. Similar to a DETOX-P3 resolution test, it is an inversion where the synthetic 'data' consist of only the (sign-reversed) crust-plustopography correction times δT = −δT corr = −(δT crust + δT topography ) while retaining the full P+Pdiff+PP kernel matrix. The resulting, inverted mantle "model" should consist purely of those artefacts of downward smearing that would be introduced into an inversion of the real DETOX-P3 data if it had accidentally omitted all crust and topography corrections. These artefacts are shown in Fig. A4, and they are indeed highly correlated with the slow structures that appear in DETOX-P3 but not DETOX-P2 (highlighted by ellipses), and have roughly the same magnitude of dV P /V P . The anomalies smear down near-vertically to >1200 km because almost everywhere under the oceans, PP waves sample too sparsely to produce abundant crossings of their steeply incident bouncing legs.
From this test result, we cannot conclude that there is a problem with crustal/topography corrections for PP bounce points. Due to the steeply incident geometries of PP, the same inversion result would be produced by slow-velocity anomalies that were truly present in the (upper) mantle beneath PP bounce points, that is, at depths not too much deeper than the crust. This sensitivity ambiguity of PP makes it nontrivial to assess the validity of the differences between DETOX-P2 and -P3. We reiterate our conclusion that most likely the difference has a real structural cause-at some depth range between the asthenosphere and 1500 km, most likely between 500 and 900 km, dV P /V P in the oceanic mantle is significantly slower Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 than indicated by "continental" mantle reference model IASP91because this finding could plausibly have remained undetected until now and because we have thoroughly vetted our traveltime correction procedure for implementation mistakes.
In the following, we interpret whole-mantle sections in areas of suspected mantle plumes, with a focus on demonstrating the benefits of Pdiff data. In order to remove any lingering doubts about PP from the discussion, we will interpret DETOX-P2 instead of DETOX-P3, given that any vertical smearing of slow structure along sparse PP paths could suggest the presence of plume conduits where there aren't any.

Large low velocity provinces (LLVPs)
In the lower third of the mantle, the southern hemisphere is dominated by large-scale slow velocity anomalies (Figs 13, 14 and 16). At long spatial wavelengths, these structures have long been prominent in global S-wave tomographies (e.g. Woodhouse & Dziewonski 1989;Ritsema et al. 1999Ritsema et al. , 2011Montelli et al. 2006;Houser et al. 2008;Simmons et al. 2010;Lekic et al. 2012;Auer et al. 2014;Cottaar & Lekic 2016;Koelemeijer et al. 2016;Garnero et al. 2016), which have stressed their 'degree 2' partitioning into two 'Large Low Shear Velocity Provinces (LLSVPs)' centred under southern Africa and the South Pacific, and separated by belts of fast-velocity anomalies. Since our Pdiff inversions now also show these structures very clearly, we instead adopt the term 'Large Low Velocity Province (LLVP)', thus dropping the qualifier 'Shear'.
However, a two-partite structure is relatively less dominant in our results because much more spatial detail and subpartitioning are apparent. The resolution tests of Fig. 13 show that the structural details discussed below should be readily resolved.
In Fig. 14, the African LLVP at 2800 km depth is imaged as an assemblage of three intensely slow patches centred under southwest Africa (lat -13.5 • , lon 12.2 • ); under northwest Africa (lat 29 • , lon -2 • ) with an extension under Europe; and under the Equatorial Atlantic (Ascension Island region, lat -5 • , lon -10 • ). Their approximate diameters are 1000, 1000 and 600 km, respectively. Also considered part of the African LLVP is a long, NW-SE striking tongue of moderately slow material that stretches from Madagascar to Antarctica (lon ≈110 • ); as well as a cluster of moderately slow anomalies beneath Arabia, Afar and the northwestern Indian Ocean. Our tomography shows these slow patches clearly delineated from each other, that is, separated by (relatively narrower) regions of seismically fast or almost neutral mantle. An additional slow patch of this character is the so-called 'Perm anomaly' under Russia, which appears as a distinct feature in the vast majority of global models (Lekic et al. 2012).
Under South America, DETOX-P2 and -P3 observe an additional group of slow anomalies near the CMB that have not appeared as clearly in prior work. In the southwestern Atlantic offshore Brazil, Uruguay and Argentina, a very slow patch is centred on 33 • S, 50 • W, (diameter ≈1000 km) and connects to moderately slow anomalies beneath much of South America. This area is well resolved according to the resolution tests of Figs 13 and 10. While inspection of prior global tomographies shows only vague evidence of these patches ), a number of regional studies of the lowermost mantle seem to have picked up on parts of these features (Wysession et al. 2001;Fisher et al. 2003;Hung et al. 2005). The most interesting aspect about these newly imaged anomalies is that they indicate a relatively smooth longitudinal transition between the African LLVP patches described earlier and slow anomalies beneath the southeastern Pacific (discussed below), which have been counted among the Pacific LLVP by prior work. The discovery of South American slow anomalies thus blurs the distinction between the two African and Pacific 'superplume' domains, a division that had been considered solid, and a crucial observation for geodynamic convection modelling to match.
A map view centred on the Pacific at 2800 km depth (Fig. 16) reveals two large areas of intensely slow anomalies. In an almost north-south striking belt from the Mexican west coast to ≈40 • S, dV P /V P appears almost uniformly slow by more than 1 per cent. This slow CMB belt roughly underlies the East Pacific spreading ridge although there is no suggestion of an upward connection to this near-surface feature. Further west, a vast area beneath the equatorial and southern Pacific is occupied by half a dozen intensely slow patches set against a moderately slow background. This description fits the longitudes between ≈150 • W (Hawaii) and ≈135 • E (eastern Australia and Papua-New Guinea). The very slow patch just west of Hawaii coincides with an Ultra-Low Velocity Zone first described by  from Sdiff waveform analysis. The overall region has been considered to form the core of the Pacific LLVP, especially since it was more evident in many tomographies than the previously discussed, N-S striking belt beneath the eastern Pacific. The latter is however very prominent and distinct in the DETOX models, and is located at a similar distance from the Pacific LLVP 'core' as from the newly imaged anomalies under South America. This underscores the longitudinal continuity of slow CMB patches around the southern hemisphere.
In summary, seismically slow anomalies in the lowermost mantle are concentrated at southern and equatorial latitudes and can be described as moderately slower-than-average background, onto which are superimposed a dozen or more intensely slow patches of typically 600-1400 km diameter. By contrast, seismically fast structure (subducted lithosphere) is concentrated under the northern hemisphere, and under Antarctica and the Southern Ocean (see Section 6).
The DETOX-P2 and -P3 models revise current views of the lowermost mantle in that we observe an almost globe-spanning belt of intensely slow patches at southern and equatorial latitudes, rather than two clearly separated LLVPs beneath Africa and the Pacific. The only major perturbation of this belt occurs at the longitudes of eastern Asia, where a massive accumulation of fast anomalies under the northern hemisphere spills southward into the Indian Ocean, nearly severing the spatial continuity between the slow, Madagascarto-Kerguelen tongue (African LLVP) and the westernmost part of the Pacific LLVP under Australia.
High resolution in the lowermost mantle directly results from the Pdiff measurements (see Sections 4.1 and 4.3); no structure is imaged in DETOX-P1 with teleseismic P measurements only. The resolution tests of Fig. 13 show that the input models can be retrieved very well in the lower-third of the mantle (also see Fig. 10). The outlines of inputs have been preserved, and many small-scale features have been resolved. However, amplitude recovery is spatially variable. The amplitudes of the input models are underestimated in the outputs, particularly ≈ -60 • latitude in South Atlantic and south of Africa due to the lack of data and regularization that we use to suppress the effects of noise in the data. Figure 13. DETOX-P2 at 2800 km depth: preferred model (top right) and three resolution tests. Input 1 consists of Gaussian spheres spaced by 30 • , with peak anomalies of 3 per cent. Input 2 consists of cylindrical Gaussian anomalies dV P /V P = (dV P /V P ) center exp ( − r 2 /w 2 ) with radii w = 400 km, i.e., synthetic mantle plumes inserted at the locations of intraplate hotspots (magenta triangles in top right map). Input 3 is the same except that w = 300 km. The output results of tests 2 and 3 suggest that plume clusters (e.g. in the southern Pacific or offshore western Africa) do not blur together completely at 2800 km depth, and hence that internal structure of the Pacific and African LLVPs should be resolved. The actual tomography indeed shows significant subdivisions within the two LLVPs, which could be the footings of individual mantle plumes.

Mantle Plumes
Figs 14 and 16 show 10 vertical whole-mantle cross sections (panels A to P), which are chosen to intersect the mantle beneath volcanic hotspots, the presumed locations of deep mantle plumes. In each case, we compare DETOX-P2, which contains Pdiff data, to DETOX-P1, which is estimated from teleseismic P traveltimes only. The limits of the African and Pacific LLVPs are shown in panels A-P as black dashed lines labeled 'AL' and 'PL', respectively. They were traced subjectively in DETOX-P2 rather than by any rigorous algorithm, since their purpose is merely to give a visual reference highlighting the differences between DETOX-P1 and DETOX-P2.
Each section A-P is accompanied by a resolution test for a vertical, cylindrical mantle plume (slow columnar anomaly) in Figs 15 and 17.
The DETOX inversions include data from several temporary seismic arrays that were dedicated to improving the sampling of the (upper) mantle beneath certain hotspots. This includes stations around La Reunion island (Barruol & Sigloch 2013;Stähler et al. 2016 Fig. 14. Cross-sections through Ascension, Iceland, Afar, Kerguelen, Azores and Canary, from surface to CMB. The input patterns are 3-D cylinders with a Gaussian lateral velocity profile dV P /V P = (dV P /V P ) center exp ( − r 2 /w 2 ) with radius w = 300 km and (dV P /V P ) center = −3 per cent; centered on hotspot locations marked by pink triangles on the map. r is radius from the center of each 3-D symmetric Gaussian velocity anomaly. Colours indicate dV P /V P with respect to IASP91. For each hotspot, the left panel shows the input pattern, and right panel shows the output pattern retrieved by the DETOX-P2 inversion, using realistic regularization. Beneath Iceland and Afar, the columns are recovered very well at almost all depths. Input patterns below Azores and Canary are well resolved but the amplitudes are slightly underestimated. Even though the two columns are very close to each other in the input model, DETOX-P2 data set has sufficient resolution to resolve two separate cylinders in this region. Below Ascension, the outlines of the input pattern can be resolved but amplitudes in the mid mantle are underestimated and some lateral smearing is present. Note that section "D" was chosen to track the SE-striking 'tongue' of the LLVP (see the main text and Fig. 14) and does not cut through the central axis of the vertical, test cylinder, where dV P /V P reaches its maximum. See Fig. F5 for a vertical cross-section that runs through the surface location of the Kerguelen hotspot.  Fig. 17 for resolution tests beneath these hotspots. Dotted ellipses in H and K indicate likely artificial downward smearing of shallow heterogeneities beneath Hawaii and Tahiti. In panel G, high-velocity anomalies beneath Central America (subducted lithosphere, labelled 'Cocos' and 'Trans-Americas') are visible at all mantle depths. This well-sampled region is one of the few where even DETOX-P1 shows structure in the lowermost mantle (see Fig. 11), but the superior detail of DETOX-P2 at these depths is particularly clear in this comparison.  Fig. 15 for explanations of plotting style. The input pattern beneath Galapagos is well resolved. For all other regions, the amplitudes of dV P /V P are underestimated, to variable degrees as a function of depth. For Hawaii, amplitudes are least well retrieved in the mid-mantle, which mirrors the appearance of structure actually imaged by DETOX-P2 under Hawaii (in Fig. 16).
Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 In Fig. 14(A) beneath the Ascension Island area, DETOX-P2 shows a continuous connection of slow anomalies from the surface down to the African LLVP. The structure is not vertical; it tilts to the right (NE-ward) in the upper 1500 km, and back to the left (SW-ward) in the lower half of the mantle, where dV P /V P is also weaker. Interestingly, the columnar resolution test of Fig. 15 shows a very similar output, of a plume 'kinked' to the right (NE-ward) and back, with weak recovery of dV P /V P in the lower half of the mantle. DETOX-P1 closely resembles DETOX-P2 in the upper 1200 km. As expected, it fades at deeper depths due to poor sampling combined with the action of norm regularization. Only a hint of the LLVP is present, in stark contrast to DETOX-P2, which shows a particularly intense low-velocity patch beneath the Ascension area, as discussed earlier.
Beneath Iceland (Fig. 14B), DETOX-P2 images a continuous low-velocity conduit extending from the surface to the CMB, but again not exactly vertical. Excellent seismic illumination of the area means that a vertical conduit should be very well resolved at all mantle depths, according to the resolution test, Fig. 15(B). Hence the upwelling under Iceland indeed seems to originate from the lowermost mantle but is positively not vertical, which is consistent with prior work Montelli et al. 2006). DETOX-P1 is very similar to DETOX-P2 down to depths exceeding 2000 km which is another indication of the excellent seismic illumination of the area. Again the LLVP is apparent only in DETOX-P2 and is not present vertically beneath Iceland. Rather, the upwelling is rooted in a slow anomaly beneath southern Greenland, which is expressed relatively weak in dV P /V P and may or may not form part of the LLVP. For the upper mantle, the DETOX models and almost all other pertinent tomographies agree on the existence of an extended and seismically very slow anomaly beneath Iceland (e.g. Ritsema et al. 1999;Hung et al. 2004;Montelli et al. 2006;Rickers et al. 2013). The role of additional low-velocity anomalies in the upper 1000 km, east of Iceland and abutting the prominent East European craton (dark blue), remains up for speculation, for example that the plume flow may have encountered resistance to vertical ascent into the upper mantle (French & Romanowicz 2015).
Beneath East Africa, the Afar area and Arabia in Fig. 14(C), a vast tract of the upper 1000 km of mantle are filled with very slow anomalies. By comparison, the LLVP is shifted SW-ward and is located beneath South Africa. A possible connection is made by a slow anomaly tilting up and NE-ward between 2000 and 1000 km depth, a puzzling situation that has been resolved previously (e.g. Ritsema et al. 1999Ritsema et al. , 2011Hansen et al. 2012). DETOX-P2 images a second very slow patch on the CMB beneath Arabia, but no clear upward connection to the Afar area, whereas the resolution test of Fig. 15(C) indicates that a vertical plume beneath Afar would be quite resolvable.
The long LLVP tongue that strikes south-eastward from southern Africa through the Southern Ocean to Kerguelen is sectioned by Fig. 14(D). Under southern Africa, the LLVP rises from the CMB to depths as shallow as ≈1000 km-very well resolved by DETOX-P2 and partly also by DETOX-P1. Its elongate appendage into the Southern Ocean remains mostly confined to the lowermost 500 km except beneath the Kerguelen area, where the LLVP fills the lowermost 1000 km. Kerguelen is cut by two sections in Fig.  F5, which show slow anomalies extending upward to the Kerguelen hotspot. The resolution tests indicate that a vertical plume would be confidently resolved by DETOX-P2.
The Canary and Azores areas are sectioned by Fig. 14(E), which cuts tangentially along the western African margin. DETOX-P2 shows slow anomalies in the upper mantle beneath the two hotspots, which extend downward across all depths, albeit not in a simple, vertical geometry. The LLVP is sharply delineated in the lowermost 500 km but is shifted SE-ward relative to shallower slow anomalies. The columnar resolution tests of Fig. 15(E) indicates that DETOX-P2 has good resolution beneath these islands although d V p /V p amplitudes are underestimated.
Shifting focus to hotspots of the Pacific basin, a NE-striking cut through Easter Island and the Galapagos archipelago reveals a complex picture in Fig. 16(G). According to the resolution test of Fig. 17(G), a vertical plume beneath Galapagos would be very well resolved, except for moderate underestimation of its amplitude between 700 and 1500 km depth. The actual tomography shows slow anomalies only in the upper 1000 km beneath the hotspot. The lowermost mantle is filled by fast anomalies ('Trans-Americas slab'), which are so robust that they are picked up even by DETOX-P1, and are agreed upon by the vast majority of global tomography models . Hence the Galapagos upwelling cannot originate from the CMB vertically beneath the hotspot. Nolet et al. (2019) recently reported a 200-300 km thin, almost vertical mantle plume beneath Galapagos, from the surface to ≈1900 km depth, but no deeper. Their study is consistent with our conclusion in that they also rule out an origin from the CMB (though without explicitly tying to a slab), and that our tomography would likely miss a thin plume between 1000 and 1900 km depth because it does not include the MERMAID float seismograms (Sukhovich et al. 2015) that constituted the major novelty of the Nolet et al. (2019)

study.
A vertical plume beneath Easter Island would be detected confidently only in the uppermost and lowermost mantle, whereas intermediate depths would be afflicted by artificial NE-ward kinking of the plume conduit (like for Ascension Island) and by severe attenuation of dV p /V p amplitudes-see Fig. 17(G). The actual tomography in Fig. 16(G) resembles the resolution test output in that the lowermost and uppermost 700 km beneath Easter are intensely slow, and that the upper-mantle anomaly tilts down towards the northeast, with relatively weak amplitude. All this would be the expected imaging outcome if the actual plume conduit were vertical, but our result is obviously not sufficient to confirm the vertical plume scenario. Fig. 16(M) shows an almost perpendicular cut through Easter in which the plume conduit looks more nearly vertical but also interrupted, its 'kink' protruding out of the viewing plane. Again this result resembles the output of the vertical plume resolution test, Fig. 17

(M).
Hawaii is cut by two almost perpendicular sections in Figs 16(H) and (K). The resolution tests, Figs 17(H) and (K), show that a vertical plume would be confidently resolved by DETOX-P2 in the upper 1000 km and the lower 1000 km, but would appear very strongly attenuated throughout the middle third of the mantle. The actual tomography result in Fig. 16(K) closely resembles the test output, and Fig. 16(H) bears a reasonable resemblance, both of them showing strong slow anomalies in the upper and lower 800-1000 km, and weak slow anomalies at intermediate depths. (Slow anomalies enclosed by ellipses in (H) should be ignored. They are ray-shaped and very likely represent downward smearing artefacts, caused by lack of crossing ray paths beneath Hawaii.) Hawaii's upper mantle anomaly appears to be centred slightly southwest relative of the hotspot, and the imaging results are consistent with a vertical whole-mantle plume beneath this location, which in the lowermost mantle hosts a large ULVZ .
The hotspot archipelago of French Polynesia/Tahiti is also sectioned by Fig. 16(H). It shows very slow anomalies in the upper 1000 km and the lowermost 500 km. Again this resembles the resolution test in Fig. 17(H), which shows that a vertical whole-mantle conduit would be difficult to resolve at depths between 1000 and ≈2200 km (although note that the section misses the central axis of the vertical cylinder, giving the test output an overly pessimistic appearance).
The Marquesas archipelago is sectioned by Fig. 16(P) and shows surprisingly strong slow anomalies, given that the columnar resolution test of Fig. 17(P) indicates severe attenuation of dV P /V P at all depths except the lowermost 500 km. The actual tomography does show an intensely slow LLVP, but also relatively strong anomalies in the uppermost mantle, and moderately slow dV P /V P to below 1000 km. All of these are consistent with a whole-mantle plume in this poorly sampled area. In the regional tomography model of Obayashi et al. (2016), large-scale slow anomalies from the CMB to 1000 km depth are imaged underneath the South Pacific superswell consistent with DETOX-P2.
In this section, we have stressed tests for vertical, whole-mantle plume conduits because this represents the most influential and long-standing hypothesis for mantle structure beneath intraplate hotspots. In a few cases, the addition of P-diffracted waves to conventional teleseismic data sets was shown to confidently rule out the null hypothesis of a vertical, substantial, whole-mantle plume because seismically fast or neutral anomalies can now be robustly imaged near the CMB beneath the hotspot (e.g. under Galapagos, Iceland, Afar, Azores-Canary). In other cases, both the upper and lowermost mantle beneath a hotspot are observed to be seismically very slow, and resolution tests indicate that lack of slow structure imaged at intermediate depths would be expected from a lack of sampling at these depths (e.g. Ascension, Easter, Hawaii, Tahiti, Marquesas). This leaves open the possibility of a vertical, substantial, whole-mantle plume, and in fact strengthens this option compared to prior work, by confidently confirming a slow lowermost mantle. The investigation of more complex plume models (nonvertical, intermittent, thin, etc.) will be the subject of future work (e.g. Tsekhmistrenko et al. 2018).

S E I S M I C A L LY FA S T A N O M A L I E S , T H E R E M N A N T S O F PA L A E O -S U B D U C T I O N
The DETOX models feature a large number of seismically fasterthan-average anomalies that are well resolved and delineated, and which are broadly consistent with prior work. In Fig. 11, fast (blue) anomalies in the lithosphere (depth 100 km) coincide with cratons, shields and other old continental lithosphere (the oceans are largely unsampled) and also with active subduction zones (e.g. Andes, Sunda). By 600 km, the continental lithospheric signature is fading in favour of slabs beneath active subduction zones. Below these depths, fast anomalies are widely accepted to represent the thermal anomalies associated (only) with subducted lithosphere. At 1200 km, subducted lithosphere is concentrated in two vast belts, one beneath the Americas (corresponding to the accretions of the American Cordilleras) and one beneath the Tethyan mountain chains, from the Alps to the Southwest Pacific. Fig. 18 presents a 3-D rendering of seismically fast anomalies between 600-1800 km depth in DETOX-P1, a representation that highlights the continuity of these two huge slab complexes in the mid-mantle, laterally and in depth. Their spatial correlation with the two major orogenies of Mesozoic times has played a decisive role in the acceptance of whole-mantle convection (Grand et al. 1997;Bijwaard et al. 1998;Van der Voo et al. 1999). Detailed geological interpretations of these mid-mantle slab complexes have inferred approximate sinking rates of 1 cm/a ( Van der Voo et al. 1999;van der Meer et al. 2010;Sigloch & Mihalynuk 2013van der Meer et al. 2018), from which one can extrapolate that the mantle holds a memory of at least 300 Myr of subduction.
By 2200 km depth in Fig. 11, DETOX-P1 has started to fade and to diverge from DETOX-P2 and -P3, and fast anomalies are less obviously structured in linear belts. The lowermost mantle (2800 km) presents 'cleaner' patterns once again, with one enormous accumulation of slab under the eastern half of Asia, which must represent the Paleozoic assembly of Eastern Asia (e.g. Şengör et al. 1993;Domeier 2018). Secondary belt-like anomalies, still thousands of kilometres long, are imaged beneath the Americas, the Atlantic, the North Pacific and the Southern Ocean/Antarctica. Visually the linear fast anomalies in the lowermost mantle appear significantly widened compared to, for example, at 1200 km, but this is deceptive because the area of the CMB is only about a quarter as large as the earth's surface area.
In DETOX-P2/P3, the resolution gains for slabs in the lowermost mantle (compared to DETOX-P1 or existing P-wave tomographies) are similar to the gains discussed for seismically slow patches near the CMB. This can be appreciated on the example of Fig. 19, which presents a zoom into the American hemisphere, where voluminous slabs are present from the surface to the CMB.
The DETOX models include finite-frequency traveltimes for most of the USArray deployment (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), resulting in similarly high resolution in the upper half of the mantle as that achieved by dedicated USArray tomographies (Tian et al. 2009;Schmandt & Humphreys 2010;Sigloch 2011;Pavlis et al. 2012;Burdick et al. 2014Burdick et al. , 2017. Down to 1500 km depth, DETOX-P1 and DETOX-P2 differ very little in Fig. 19, but below that the differences accumulate, with slabs at 2200 km in DETOX-P1 appearing like faint, smeared shadows of those in DETOX-2. These resolution differences are all the more significant because a major reorganization of slab locations is evident from 1500 to 2200 km. Indeed, slabs above ≈2000 km have been used extensively to infer palaeogeography and even to reinterpret the geologic record (Grand et al. 1997;Grand 2002;Ren et al. 2007;Sigloch et al. 2008;van der Meer et al. 2010;Sigloch 2011;Pavlis et al. 2012;Sigloch & Mihalynuk 2013van der Meer et al. 2018;Chen et al. 2019), whereas only rough proposals have been made for deeper slabs (van der Meer et al. 2018).
The lowermost mantle beneath Central America is one of the best sampled regions in this depth range, thanks to dense seismic networks in North America recording abundant seismicity from the Andes. Even DETOX-P1 robustly picks up on high-velocity structure beneath Central America at 2800 km depth, from deep-diving teleseismic waves alone (Fig. 19 bottom right). This slab is equally evident in other global P-wave tomographies (e.g., Grand et al. 1997;Montelli et al. 2006;Amaru 2007;Houser et al. 2008;Li et al. 2008;Simmons et al. 2012), the regional finite-frequency tomography of Hung et al. (2005), and array seismology studies (e.g., Fisher et al. 2003;Thomas et al. 2004;Hutko et al. 2006;Kito et al. 2007). However, DETOX-P2 shows this same slab in much crisper definition, while also picking up on a whole network of elongated, neighbouring slabs, which offers the prospect of reconstructing a network of continuously closing plate boundaries, of the same character as observed on the surface today.
A similar point can be made on the mantle beneath Europe and surrounding regions (Fig. 20). The upper 900 km show no appre-Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure 18. Seismically fast anomalies in the mid mantle (600-1800 km), as imaged by DETOX-P1. 3-D isosurfaces enclose structures that are dV P /V P = +0.25 per cent faster than ambient mantle at any given depth; slow and seismically neutral features are not rendered. Colour signifies depth and changes every 200 km. The shallowest structures (at 600-1000 km) appear in the foreground in turquoise and greenish colours; behind them lies structure at 1000-1200 km in yellow, 1200-1400 km in red, etc. The map is clipped at absolute latitudes ≥75 • . Fast anomalies at sublithospheric depths are generally accepted to represent subducted (cold) oceanic lithosphere. The bulk is located beneath the Americas, corresponding to the Mesozoic-Cenozoic Cordilleran orogenies, and beneath the southern margins of Eurasia, corresponding to the Tethyan orogenies of the same ages.
ciable difference between DETOX-P1 and -P2, whereas a wholemantle section through the Aegean subduction zone convincingly shows (in DETOX-P2) that the Aegean slab does not extend to the CMB, which is instead occupied by slow anomalies. This could not have been inferred from DETOX-P1, which fades below ≈1800 km depth. At 100 km beneath Europe, sharply delineated boundaries between fast and slow anomalies are prominent beneath known, lithosphere-cutting boundaries, most prominently the Tornquist-Teisseyre Zone, which separated the Precambrian Eastern European Platform from younger terranes that accreted in Phanerozoic times (e.g. Zielhuis & Nolet 1994;Zhu et al. 2012). At 600 km depth, Mediterranean slabs are similarly sharply delineated with slab-free regions to their north. This shows how teleseismic waves can yield quite impressive lateral resolution under regions that are densely instrumented at the surface. An analogous case could be made for North America in Fig. 19, where the upper mantle is mapped out even better by USArray.

D I S C U S S I O N
DETOX-P1 features very low dV P /V P amplitudes at 2800 km depth because illumination is based only on teleseismic wavepaths of less than 85 • and 90 • epicentral distance in finite-frequency measurements and analyst-picked arrival times, respectively. Such waves sample the lower few hundred kilometers only marginally (Fig. 9, bottom left), and norm regularization acts to suppress dV/V structure in undersampled areas (Fig. 10, bottom right). By comparison, DETOX-P2 shows strong anomalies because it includes Pdiff traveltimes, and these anomalies are robustly imaged (Fig. 10, bottom center). However, the other three global P-models in Fig. 21 show comparably strong dV P /V P anomalies as DETOX-P2, even though none of them explicitly included a large Pdiff data set (nor normal modes). There are at least three-and-a-half possible explanations.
Most likely the models have included P picks up to ≈97 • , the distance at which the Pdiff range nominally starts in ray theory (for shallow events). However, numerical experiments show that a P-waveform containing periods between 30 s and 3 s displays diffractional character already at distances exceeding ≈85 • (Hosseini & Sigloch 2015), which means that nominally 'teleseismic' P picks between ≈85 • and ≈97 • are actually measured on seismograms that already sense the core, with very non-ray-like sensitivities. Tomographies including these waves have truly sampled the lowermost mantle, but this sampling is modeled by a very crude approximation (rays), translating into a commensurate image blurriness. DETOX-P2 and -P3 are expected to contain a milder version of this modelling blurriness because its Pdiff sensitivity kernels, while not rays, are still based on Dahlen's approximate method  for teleseismic waves Hung et al. 2000;Tian et al. 2007b), plus a heuristic 'flattening' interpolation across the CMB. Hence additional gains in resolution are expected once sensitivity kernels will be calculated by the AxiSEM methods , as synthetic seismograms for the measurements already are.
Moreover, the apparent illumination of the lowermost mantle can be increased by choosing coarse inversion grid cells that span large depth ranges (e.g. a tetrahedral node spacing of 600-800 km in the lower mantle in Montelli et al. 2006). A truly teleseismic wave sampling the shallow part of such a cell will determine dV/V for the entire cell in absence of deeper diving waves. Thus shallower structure gets projected into the deepest mantle, which may or may not be readily apparent in the end result. In Fig. 21 for example, PRI-P05 at 2800 km shows the imprint of individual tetrahedra (straight division lines between red and blue domains), which hint at this 124 K. Hosseini et al. Figure 21. Comparison of global P-and S-wave tomographies at 2800 km depth. DETOX-P1 and DETOX-P2 are compared to three other P-models: PRI-P05, the first global finite-frequency model, which also contained a large set of EHB picks (Montelli et al. 2004b); MITP08, computed from EHB P-picks (Li et al. 2008); and GAP-P4, the only other global finite-frequency model (Obayashi et al. 2013). Also shown are global S-wave models SEISGLOB2 (Durand et al. 2017), S40RTS (Ritsema et al. 2011), and SEMUCB-WM1 (French & Romanowicz 2014). Colours indicate velocity perturbations with respect to the original reference models (rather than one common model). Colour saturation values are 1 per cent for P models and 2 per cent for S models. Bottom: Vote maps (Shephard et al. 2017;Hosseini et al. 2018) for 2800 km depth, generated from 11 global P-wave and 14 global S-wave models (listed below). Left map shows model agreement on seismically fast structure (slabs), right map on seismically slow structure (LLVPs). Red and black areas indicate that majority of models agree that fast or slow structure is present; greenish-yellow means that about half of the models agree. Models seem to agree easily on the outlines of the African LLVP, and somewhat less on the Pacific LLVP. Visual comparison to the eight individual models also shows that the apparent consensus in the vote map (which ignores amplitudes) glosses over rather variable appearances in individual models. The fast-velocity vote map shows much less agreement than the slow vote map, and individual models are obviously more variable in their blue patches. Seismically fast structure, has always been more challenging to image in the lowermost mantle, a major motivation for adding more solid constraints from Pdiff waves. The 11 global P models included in the votemaps: DETOX-P2, PRI-P05, MITP08, GAP-P4, UU-P07 (Amaru 2007), Hosseini2016 (Hosseini 2016), GyPSuM-P (Simmons et al. 2010), HMSL-P06 (Houser et al. 2008), SP12RTS-P (Koelemeijer et al. 2016), SPani-P (Tesoniero et al. 2015), LLNL G3Dv3 (Simmons et al. 2012); and 14 global S models: S40RTS, SEMUCB-WM1, GyPSuM-S (Simmons et al. 2010), HMSL-S06 (Houser et al. 2008), PRI-S05 (Montelli et al. 2006), SP12RTS-S (Koelemeijer et al. 2016), SPani-S (Tesoniero et al. 2015), S20RTS (Ritsema et al. 1999), S362ANI+M (Kustowski et al. 2008), SAVANI , SAW642ANb (Panning et al. 2010), SEMum (Lekić & Romanowicz 2011), SGLOBE-rani (Chang et al. 2015), TX2015 (Lu & Grand 2016). phenomenon. Large grid cells are a pragmatic approach to improving the conditioning of the inverse problem, but of course create an attendant blurriness about the true spatial origins of traveltime anomalies, which may not include the lowermost mantle despite the model's appearance. (Finally, sufficient relaxation of regularization parameters would always result in strong dV/V anomalies in the undersampled lowermost mantle, but we have no reason to expect this instability in any of the models shown.) In summary, most earlier global P-models did actually sample the lowermost mantle, but each of the above effects results in decreased image resolution, for which the cumulative effect is difficult to assess. By adding a large Pdiff data set that samples the lowermost mantle extensively rather than marginally, using proper sensitivities, small grid cells and conservative regularization, these above effects should be minimal in the DETOX models, and the appraisal of imaging results by resolution tests is as straightforward as for other mantle regions.
In shear-wave tomography, constraints on the lowermost mantle derive mainly from normal mode data (Ritsema et al. 1999(Ritsema et al. , 2011Moulik & Ekström 2014;Koelemeijer et al. 2016;Durand et al. 2016Durand et al. , 2017. The availability of normal mode splitting measurements to constrain dV S /V S explains why S-velocity models have been the preferred avenue for investigating lowermost mantle structure and dynamics, perhaps most prominently reflected by the term 'Large Low Shear Velocity Province (LLSVP)'. On the downside, normal modes constrain only spatial scales of many hundreds to thousands of kilometers, and resolve only structure of even degrees (2, 4, 6, etc.), unless mode coupling is taken into account. Hence almost all normal mode inversions have been insensitive to that half of mantle structure which cannot be rendered by even-degree spherical harmonic functions-a large and non-intuitive blind spot. Specifically, one has to entertain the possibility that two pronounced, separate LLSVPs have featured so prominently in S-wave models because degree-2 structure is what the method is sensitive to. Where mode coupling has recently been taken into account (Durand et al. 2016, the dominant degree 2 structure of the LLSVPs fades in favour of a joint dominance of degrees 2 and 3. This result points in the same direction as the emergence of additional slow patches under South America in our Pdiff inversions, outside the previously described LLSVP boundaries ('the plume generation zone'; Burke et al. 2008).
Comparison of the models in Fig. 21 shows that the South American low-velocity anomalies at 2800 km that are seen so clearly in DETOX-P2 (and vaguely also in P-models PRI-P05 and MITP08), are not apparent in the S-wave models (except possibly SEMUCB-WM1). Sensitivity and resolution issues are a likely explanation, but on the other hand, the mantle's true V p anomalies need not perfectly covary with its V s anomalies, given our uncertainties about rock compositions and phase transitions in the lowermost mantle (e.g. Masters et al. 2000;. Hence demonstrating robust differences between P-and S-wave models is of high scientific interest, but challenging due to the completely different sensitivity patterns of diffracted body waves versus normal modes. The most straightforward and transparent approach would be side-by-side and joint Pdiff-Sdiff inversions, which have not been attempted. A few global S-wave models (e.g. Mégnin & Romanowicz 2000;Ritsema et al. 2011;Koelemeijer et al. 2016) have already complemented their normal mode data by core-diffracted S waves. These Sdiff data sets have been much smaller than the 1.4 M Pdiff traveltimes used in the DETOX models, and only a few global S models have been published that uses Sdiff measurements but not normal modes to constrain the lowermost mantle.

C O N C L U S I O N S
In providing very detailed and robust resolution of the lower third of the mantle, our new models DETOX-P2 and DETOX-P3 remove a major blind spot of existing global P-wave tomographies. This is achieved by measuring and modelling ≈1.4 M multifrequency traveltime measurements of core-diffracted P waves (Pdiff), where previous work included no more than a few hundred data. Inclusion of Pdiff data also stabilizes structures retrieved in the midmantle, which are constrained primarily by teleseismic P and PP waves. The latter add valuable constraints to DETOX-P3, but the shortcoming (relative to DETOX-P2) is that PP waves have also introduced significant low-velocity artefacts in the upper ≈1800 km beneath oceans, which a user of DETOX-P3 needs to be aware of (or use DETOX-P2 instead). Forward modelling for the >5 M multifrequency measurements is done by full wave propagation up to 1 Hz frequency and is not computationally limited, thanks to the efficiency of the AxiSEM software ; van Driel & Nissen-Meyer 2014a,b;van Driel et al. 2015), which is based on the adoption of a spherically symmetric background model. Future gains in resolution are expected once sensitivity kernels will also be calculated by AxiSEM rather than the more approximate method of Dahlen et al. (2000), Hung et al. (2000) and Tian et al. (2007b).
The most significant new findings arguably concern the arrangement of slow and very slow velocity anomalies in the lowermost mantle. Rather than just two Large Low Velocity Provinces under Africa and the Pacific, we resolve their subdivisions into at least a dozen very slow patches of 600-1400 km diameter, embedded in a background of moderately slow anomalies that occupies much of the southern and equatorial latitudes. The first clear resolution of additional slow patches beneath South America closes one of the two longitudinal gaps between the African and Pacific LLVPs. This leaves only one gap beneath the Indian Ocean, where a vast pile of subducted slabs under Eastern Asia protrudes southward. Hence we would characterize low-velocity anomalies in the lowermost mantle as an almost globe-circling chain of intensely slow patches beneath the southern hemisphere, set against a moderately slow background occupying much of that hemisphere.
The DETOX-P2 and P3 models support or are consistent with substantive whole-mantle plumes beneath Iceland, Hawaii, Afar, Ascension, Kerguelen, Azores, Canary, Easter, Galapagos, Tahiti and Marquesas. The plumes are rooted in the well-resolved LLVPs, but their upward geometries are not always near-vertical.
The DETOX models also sharpen up subducted slab piles in the lower two thirds of the mantle. These anomalies tend to look more narrowly elongated than in prior work, especially in the lowermost mantle. This is reassuring because subducted slabs are expected to reflect an interconnected network of palaeo-trenches, analogous to the present-day character of such plate boundaries.
The DETOX models are freely available and can be explored interactively via the SubMachine website, together with two dozen other global tomography models (http://submachine.earth.ox.ac.u k/, last accessed: August 2019).

A C K N O W L E D G E M E N T S
We thank editor Gabi Laske, as well as Michael Wysession and an anonymous reviewer for their careful and constructive reviews. K.H., K.S. and A.M. were supported by funding to K.S. from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 639003 DEEP TIME) and from a Philip Leverhulme Prize awarded by The Leverhulme Trust. M.T. was funded by the NERC DTP in Environmental Research award NE/L002612/1, and the Clarendon Award from St. Edmund Hall College, University of Oxford. K.H. acknowledges the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (http://www.lrz.de). We thank Mitch Mihalynuk for discussions which included the suggestion of the name "DETOX". K.H. thanks Simon Stähler and Ceri Nunn for providing feedback on an early version of this paper. Our data processing made extensive use of the ObsPy toolkit (Beyreuther et al. 2010;Megies et al. 2011;Krischer et al. 2015). The IRIS Data Management Center provided access to waveforms, metadata, and derived products. All authors acknowledge discussions within the COST Action ES1401 "TIDES".

A P P E N D I X A : T R AV E LT I M E C O R R E C T I O N S A1 Crust, topography and ellipticity corrections
Synthetic seismograms are calculated on the spherically symmetric earth model of IASP91 in which neither effects of the crust nor ellipticity and topography are incorporated; only a layer with constant thickness and velocity is considered as a crust. To account for the known deviations from spherically symmetric earth models and to bring the predicted time closer to that of a 3-D Earth, corrections should be applied to the computed traveltimes (Nolet 2008): where (T obs − T pred ) is the measured traveltime anomaly (Section 2) and δT corr is the correction term.
We use CRUST2.0 (Bassin et al. 2000) to correct for the 3-D crustal structure. Fig. A1 shows this for Pdiff as an example. To Figure A1. Map of crustal traveltime corrections for Pdiff waves based on crustal model CRUST2.0. These values depend only on receiver location (not epicentral distance) because all Pdiff waves arrive at the same (critical) incidence angle. This map can also give a rough impression of the regional variation of crustal corrections for (non-critical) P and PP waves (Bassin et al. 2000).  consolidate the ellipticity effects, the earth is approximated by an ellipsoid. In case of Pdiff, ellipticity correction is calculated based on the parts of the ray that goes downward and upward to/from CMB (Kennett & Gudmundsson 1996;Nolet 2008).
Crust, topography and ellipticity corrections were applied to all traveltime observations using the method of Tian et al. (2007a) and the corrected measurements were added to the d vector of eq. (2).

A2 PP traveltime corrections
Crustal correction is crucial in global tomography as it weakens the crustal smearing or leakage of shallow heterogeneities into the deeper structure. In DETOX tomographic inversions, it is formulated as (Tian et al. 2007a): where t BG is calculated according to the incident angle and the length of the ray in the crustal layer of IASP91 (Kennett & Engdahl 1991), and t 3D is computed on 2 • × 2 • crustal model CRUST2.0 (Bassin et al. 2000). Fig. A2 shows average crustal corrections δT crust of the bouncing points of all PP measurements in our data set on a grid with 360 × 360 blocks. In Fig. A2(b), average crustalcorrection values lower than 2 s (i.e. |δT crust | <2 s) are masked out. Except for Tibet and Andes regions, almost all corrections with high absolute values are beneath oceans. This is due to the large differences between the background model IASP91 and the crustal model CRUST2.0 in oceanic regions. The IASP91 background model is a spherically symmetric velocity model based on the arrival time picks of major seismic phases reported to the ISC between 1964 and 1987 (Kennett & Engdahl 1991). Constrained mainly by teleseismic P picks, IASP91 is geographically biased towards mantle beneath well instrumented regions, that is, beneath the continents. By construction, it features a 35-km-thick layer of continental crust everywhere (Kennett & Engdahl 1991). Fig. A3 shows the discrepancies between CRUST2.0 and IASP91 in terms of crustal thickness, which are largest under the oceans. Synthetic seismograms for the finite-frequency correlations are computed through IASP91. If we had somehow made a mistake in implementing crustal corrections prior to tomography, that is, in not properly 'substituting' the traveltimes predicted through IASP91 crust by those predicted though CRUST2.0 crust, then wrong traveltime anomalies would enter the inversion procedure.
We designed a worst-case test to explore the outcome for the case that we had not corrected at all for the differences between IASP91 and CRUST2.0. The expectation is that the erroneous, crustal traveltime signal would be artificially smeared downwards along ray paths in sparsely sampled regions. This effect might be extreme beneath the oceans, where steeply incident PP bouncing legs tend not to cross and where the discrepancy between IASP91 and CRUST2.0 is extreme. Fig. A4 shows the outcome of this 'resolution test' and compares it to the structures observed in the DETOX-P2 and DETOX-P3 models.
The outputs of 'crustal-correction resolution analysis' (Fig. A4, right-hand column) are generated by setting the raw measurements (T obs − T pred ) in eq. (A1) to zero. In other words, the data vector in the inversion contains only δT = −δT corr = −(δT crust + δT topography ). In an ideal setup with abundantly crossing wave paths at shallow depths, the output of this test should show structure only in the upper 100 km depth. Indeed there is very little downward smearing under North America or Europe, for example. By contrast, downward smearing under the oceans is extreme, producing strong low-velocity artefacts down to > 1200 km depth. Two of the largest of these low-velocity artefacts are highlighted in Fig. A4 because they also occur in DETOX-P3 (which includes PP paths) but not in DETOX-P2 (which does not include PP paths). Their dV P /V P amplitude is comparable in the test output and in DETOX-P3. From this test, we do not conclude that we have implemented the crustal and topography corrections wrong (because we have checked them Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure A4. DETOX-P2 (left-hand column) and DETOX-P3 (middle column) tomography models compared with a resolution test (right-hand column) at different depths. The synthetic 'data' in the input model consist of only the correction terms for topography and the crustal model, that is, −δT corr = −(δT crust + δT topography ). PP crustal corrections are very large in oceanic regions because the spherically symmetric reference model IASP91 features a continental lithosphere, and bouncing segments of a PP kernel sample this 'wrong' lithosphere twice. Indeed there is no indication of similar artefacts beneath continents: note for example the non-correlation between the synthetic test and DETOX-P3 beneath the Tibetan Plateau, India, North-East Africa, Europe, Japan and western Pacific, western North America and Canadian shield at 100 km depth. In these regions, there is also agreement between DETOX-P2 and DETOX-P3, that is, structure is constrained not only by PP phases. See Fig. A2 for average crustal corrections at PP bouncing points and Fig. A3 for difference in crustal thicknesses between CRUST2.0 and IASP91. Refer to Section 4.4 for discussion. very carefully). Rather we think that the earth's mantle at relatively shallow levels beneath the oceans (probably between 500 and 900 km depth, see Section 4.4) is truly slower than indicated by IASP91, but that this signal is artificially spread over a much larger depth range by the steep geometries of PP legs bouncing under the oceans. An anomalously slow layer in the upper mantle would manifest itself in a similarly smeared fashion in DETOX-P3 as wrong crustal corrections manifest themselves in the resolution test. Refer to Section 4.4 for further discussion.

A P P E N D I X B : L I N E A R I N V E R S I O N
We assume that the traveltime observations d i (i = 1, ..., N) are linearly linked to the earth model m(r) by the traveltime kernels K T i (r ) in the following explicit form with discrete data and continuous model function: where K T i (r ) is the sensitivity kernel of ith data (Section 3.1), and the integral is over the whole earth's volume V. In practice, the volume is limited to the region where the amplitude of the sensitivity kernel is significant. This equation is in continuous form with ∞-unknowns, and first it should be discretized into a finite number of unknowns through model parametrization or generation of the inversion grid (Section 3.2). Inversion grid (linearly) interpolates the continuous kernels onto the discrete nodes of the mesh; therefore, the unknown function m(r) can be written as a linear combination of M known 'basis functions' (h j ): where m j (j = 1, ..., M) are the finite-number of unknowns in the inversion. By substituting eq. (B2) into B1: in which: A ij (i, j = 1, ..., M) embodies the physics of the problem, and it contains partial derivatives of the data with respect to the model parameters (see eq. B3). A ij links one measurement d i to all the model parameters m j according to the projected sensitivity kernel K i onto the inversion grid. It is calculated and assembled for all the measurements before the inversion (Section 3.1). We also add source-correction terms to m to simultaneously invert for hypocentral parameters (origin time, longitude, latitude and depth) and perturbations in velocity (refer to Section 3.3).
The above equation can be written in matrix form as: The method to solve this equation is described in Sections 3.3 and 3.4.

A P P E N D I X C : A P P RO X I M AT E D P D I F F S E N S I T I V I T Y K E R N E L S
Computation of the approximated sensitivity kernels for Pdiff traveltime measurements comprises two steps: (i) Dynamic ray tracing following the method of Tian et al. (2007b) to compute the geometrical spreading factors and second derivatives of the traveltime along the wavefront. As Pdiff is not a ray geometrical phase, three ray segments are instead computed for each Pdiff wave path: one downward ray to the CMB, one upward ray from the CMB, and one ray that moves along the CMB.
(ii) Compute sensitivity kernels along the three ray segments using the method of Dahlen et al. (2000). Fig. C1 shows two example Pdiff kernels at 120 • and 140 • epicentral distances. These approximated sensitivity kernels will be replaced by 'real' wavefield-based kernels once work on the necessary software package 'MC Kernel' ) is complete and optimized for large data sets as used in this study. Fig. D1 shows the spatial resolution of the tetrahedral inversion grid used in all three DETOX models, at various depths. Grid refinement is data-driven, adapting to the expected spatial resolution of the tomographic model locally (Section 3.2). The spacing of tetrahedra vertices is determined a priori from the cumulative kernel sensitivities at each vertex, which serves as a proxy for goodnessof-sampling and is shown in Fig. 9. The resulting grid is sufficiently dense to not limit the expression of the data's full information content anywhere in the mantle-compare to Fig. 6, which renders a P-wave kernel on a coarsely versus finely gridded region.

A P P E N D I X E : R E G U L A R I Z AT I O N A N D A M P L I T U D E S O F T O M O G R A P H Y M O D E L S
To better characterize the effects of regularization on amplitudes, Fig. E1 compares the mean of absolute amplitudes as a function of depth in four tomography models of Fig. 8. We excluded the model with χ 2 red = 0.93 as it is clearly contaminated by noise. The amplitudes of the solution with χ 2 red = 1.11 are higher than the other three models in all depths. Fig. E1(right-hand panel) plots the mean-amplitude ratios |m depth |/|m re f | with respect to the preferred DETOX-P2 model with χ 2 red = 1.23. Amplitudes of the more-regularized model with χ 2 red = 2.0 are 30-50 per cent of the reference model with χ 2 red = 1.23 (see Fig. E1, right-hand panel). Note that the changes are not linear in depth.
Although these values cannot be generalized to all tomography models, they show the magnitudes to which the models can be affected by the regularization in a realistic example.

A P P E N D I X F : O T H E R R E S O L U T I O N T E S T S
Figs F1-F5 show various resolution tests for DETOX models; refer to the caption of each figure for explanation.

A P P E N D I X G : C O M PA R I S O N O F G L O B A L T O M O G R A P H Y M O D E L S I N T H E S P H E R I C A L H A R M O N I C D O M A I N
Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure C1. Two examples of Pdiff sensitivity kernels at epicentral distances of 120 • and 140 • . The finite-frequency sensitivity kernels in this study are computed based on the method of Dahlen et al. (2000), which for Pdiff waves is significantly more approximate than AxiSEM kernels, but significantly less approximate than the Pdiff sensitivities used by prior work. Figure D1. Maps showing the average distance between adjacent vertices in the adaptive tetrahedral inversion grid, at 100, 600, 1200, 1800 2200 and 2800 km depth. Figure E1. Typical model amplitude (magnitude of dV P /V P ) in DETOX-P2, as a function of depth and severity of regularization. Left-hand panel: mean of absolute velocity anomalies |m depth | in four tomography models along the L-curve of Fig. 8 (χ 2 red = 1.11, 1.23, 1.51, 2.0). |m depth | at a specific depth is defined as is the number of model parameters at that depth profile. Values in the upper mantle are so small because norm damping (in body-wave tomography) acts to allocate non-zero dV P /V P only to the vicinity of sources and receivers, which are sparse in most regions. Right-hand panel: ratios |m depth |/|m re f |, where the denominator expresses the anomalies of the preferred model (hence the orange line is constant at value 1). |m re f | at a specific depth is defined as V P ) i is the velocity anomaly in the preferred model with χ 2 red = 1.23 and M d is the number of model parameters at that depth profile. Note that regularization acts to change the amplitudes in a non-constant manner as a function of depth, especially in the upper 500 km.
Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure F1. Resolution test comparing DETOX-P2 with DETOX-P3 at different depths. Left-hand column shows the input test pattern, consisting of Gaussian spheres (3-D Gaussian functions) spaced by 30 • , with peak anomalies of 3 per cent. Noiseless synthetic traveltimes were created, and the solutions for DETOX-P2 (middle column) and DETOX-P3 (right-hand column) were obtained by the same inversion procedure and parameter settings as for the real data. Figure F2. DETOX-P2 resolution test beneath the circum-Arctic region (left-hand column) and the circum-Antarctic region (right-hand column) at different depths. Technical setup as in caption to Fig. F1. The test output shows excellent near-CMB resolution beneath the circum-Arctic region, but significant resolution caveats beneath the circum-Antarctic region. Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure F3. Resolution test comparing DETOX-P2 with DETOX-P1 at different depths under Eurasia. Left-hand column shows the input test pattern, consisting of Gaussian spheres (3-D Gaussian functions) spaced by 15 • , with peak anomalies of 3 per cent. For each resolution test, the input anomalies are located at several depths to fill the whole mantle as opposed to put anomalies at one depth at a time (e.g. Li et al. 2008). Noise-free synthetic traveltimes were created, and the solutions for DETOX-P2 (middle column) and DETOX-P1 (right-hand column) were obtained by the same inversion procedure and parameter settings as for the real data. Figure F5. DETOX-P2 at 2800 km centred on the Kerguelen hotspot, and two intersecting cross-sections. Colours indicate P-velocity anomalies with respect to IASP91. Note that only Section 'd' runs through the surface location of the hotspot (pink triangle), on which the cylindrical input anomaly is centred. Section 'D' was chosen to track the SE-striking 'tongue' of the LLVP (as in the main text) and does not cut through the test cylinder axis, where dV P /V P reaches its maximum. Hence both input and output look more damped in 'D' than in 'd', but the actual test results are very similar in both sections: a whole-mantle plume could be resolved, with some caveats on amplitudes. Refer to Section 5.2 for discussion. Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure G1. Comparison of global P-and S-wave tomographies in the spherical harmonic domain as a function of depth. DETOX models are compared to three other P-models: PRI-P05 (Montelli et al. 2004b); MITP08 (Li et al. 2008); and GAP-P4 (Obayashi et al. 2013). Also shown are global S-wave models SEMUCB-WM1 (French & Romanowicz 2014), S40RTS (Ritsema et al. 2011), and SEISGLOB2 (Durand et al. 2017). Colour shades show the magnitudes of the spherical harmonic coefficients as a function of harmonic degree on the (logarithmic) x-axis, and of depth in the mantle on the y-axis. The logarithmic colour scale is chosen such that every colour increment signifies and increases in spectral power of a factor of √ 2. Refer to Fig. 12 for discussion and to Fig. 21 for a visual comparison of these models at 2800 km depth.
Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020 Figure G2. Comparison of global P-and S-wave tomographies in the spherical harmonic domain as a function of depth. Like Fig. G1, except that the x-axis is linear and that the degree zero, which indicates static offsets from reference model, is also plotted for each model. Downloaded from https://academic.oup.com/gji/article-abstract/220/1/96/5571093 by University of Oxford user on 19 March 2020