Confirmation of a change in the global shear velocity pattern at around 1000km depth

S. Durand,1,2 E. Debayle,1 Y. Ricard,1 C. Zaroli3 and S. Lambotte3 1Laboratoire de Géologie de Lyon—Terre, Planètes, Environnement, CNRS, UMR 5276, École Normale Supérieure de Lyon, Université de Lyon, Université Claude Bernard Lyon 1, 2 rue Raphaël Dubois, Bâtiment Géode F-69622 Villeurbanne Cedex, France. E-mail: durand@uni-muenster.de 2Institut für Geophysik, Corrensstr. 24, D-48149 Münster, Deutschland 3Institut de Physique du Globe de Strasbourg, UMR 7516, Université de Strasbourg, EOST/CNRS, Strasbourg, France


I N T RO D U C T I O N
The question of the existence of a transition around 1000 km depth is not new, it was actually first proposed in the Earth model MODEL A (Bullen 1940(Bullen , 1942) ) in which, based on P and S changes of velocity gradient, the Earth was subdivided into seven shells including a transition zone (shell C) ranging from 410 to 1000 km depth and a lower mantle (shell D) ranging from 1000 to 2890 km depth (Bullen 1947).The limit between the lower and upper mantle was thus first settled at 1000 km depth.Then, certainly because of the success of the association of the 410, 520 and 670 km discontinuities with olivine phase changes (Ringwood & Major 1970;Ito & Takahashi 1989), the transition zone has been reduced to 410-670 km depth range (Dziewonski & Anderson 1981).However, there still remains a large number of seismic observations showing a change in seismic properties around 1000 km depth (Johnson 1969;Kawakatsu & Niu 1994;Le Stunff et al. 1995;Fukao et al. 2001;Deuss & Woodhouse 2002;Castle & van der Hilst 2003;Shen et al. 2003;Vanacore et al. 2006;Fukao & Obayashi 2013;French & Romanowicz 2015;Jenkins et al. 2016).The question is important to solve because it has strong implications on mantle dynamics and it controls the material exchanges between upper and lower mantle.
The origin of this seismic velocity change around 1000 km depth is debated.Various explanations have been proposed.It could be due to a change in composition (Anderson 2002;Ballmer et al. 2015;King et al. 2015;Marquardt & Miyagi 2015) generally associated with a viscosity change (Forte & Peltier 1991;King & Masters 1992;Kido et al. 1998;Rudolph et al. 2015).But all these assumptions suggest a global discontinuity, while seismic studies rather report observations at specific locations.It is thus now important to investigate whether this transition is related to some particular mantle structures, high-velocity slabs or slow velocity mantle plumes and Large Low Shear Velocity Provinces (LLSVPs), or if it is global.To do so, we develop in this study a new tomographic model of the Earth's mantle.
Seismic tomography is a powerful tool to investigate the dynamic of the Earth's interior.In the last 10 yr, the resolution of whole mantle seismic models has improved (Ritsema et al. 2011;French & Romanowicz 2014;Moulik & Ekström 2014) thanks to the continuously increasing number of seismic data and the improvement of theories, from ray approximation to more sophisticated finitefrequency theories (Dahlen et al. 2000).In this study, we present a new shear velocity tomographic model of the Earth's mantle, SEIS-GLOB2, which is largely based on new surface wave, normal mode and body wave measurements that have never been jointly inverted in a tomographic inversion.In a recent study, we published SEIS-GLOB1, a pure SV model of the entire mantle designed to study the large-scale structure across the D region (Durand et al. 2016).This study uses the same surface wave and normal mode data sets as Durand et al. (2016), but includes in addition a comprehensive data set of S, SS and ScS waves measured by Zaroli et al. (2010) which provides an optimal coverage between 750 and 2000 km depth.We also use normal mode cross-coupling data that allows us to constrain the large-scale odd spherical harmonic degrees in the whole mantle, in contrast to recent global tomographic models based on self-coupling of normal modes that are sensitive to the even degrees only (Ritsema et al. 2011;Moulik & Ekström 2014).Our body wave measurements combined with self-and cross-coupling normal mode measurements and with surface wave fundamental and higher modes measured by Durand et al. (2015) provide a data set ideally suited to study the mantle near 1000 km depth.By combining different seismic data, we then maximize the geographical coverage and the depth sensitivity.The lateral velocity variations are imaged using spherical harmonics developed up to the spherical harmonic degree 40, while for the variations in the radial direction we use a combination of 21 spline functions as in Ritsema et al. (2004).
In this paper, seismic anisotropy is neglected although both radial and azimuthal anisotropy should be present (Debayle & Ricard 2012).Both Rayleigh waves and spheroidal modes are only sensitive to vertically polarized seismic waves, SV.On the contrary, our body wave data could be affected by vertical anisotropy in the lithosphere and in D and therefore be partly SH.However, the body waves that we use have a quasi-vertical incidence so that they should not be very sensitive to the lithospheric radial anisotropy and we do not include waves traveling for a long distance in D such as S diff .We therefore consider that SEISGLOB2 is mostly an SV model.
In what follows, we first present our data set and in particular the data processing and selection that have been performed before the inversion.Then, we introduce the inversion scheme that is applied to build SEISGLOB2.Finally, we present our model SEISGLOB2 and highlight the change in velocity pattern around 1000 km depth.

DATA
SEISGLOB2 includes the same surface wave and normal mode data sets as our previous model SEISGLOB1 to which we add about 100 000 S-wave traveltimes measured by Zaroli et al. (2010).In this section, we describe the three types of data incorporated in SEISGLOB2.

Phase velocity maps
In order to constrain the upper-mantle structure, we use the unique data set established by Durand et al. (2015) including some 22 000 000 Rayleigh wave phase velocity measurements with their associated errors.These measurements have been obtained applying an automated waveform inversion approach (Debayle & Ricard 2012) on some 534 359 seismograms, in the period range 40-360 s for the fundamental mode and up to the fifth overtone.
From these measurements, performed on individual seismograms, a smooth phase velocity map is produced applying a continuous regionalization algorithm (Tarantola & Valette 1982;Montagner 1986;Debayle & Sambridge 2004).This step consists in combining the measurements using an a priori Gaussian covari-ance function where an horizontal correlation length controls the smoothness of the obtained map.For more details about the method used to build the phase velocity maps, the reader should refer to Durand et al. (2015).Fig. 1 shows examples of these phase velocity maps for the fundamental mode at 40 s and the fifth overtone at 50 s.
The phase velocity maps are then decomposed into spherical harmonics Y t s (θ, φ) up to the degree 40.The obtained spherical harmonic coefficients, σ k st , where (s, t) are the spherical harmonic degrees and orders, and where k refers to a given mode, constitute the data for the inversion with depth.We also compute the associated uncertainties on the spherical harmonic coefficients that are used as data error in the inversion.Table 1 summarizes the number of measurements used for Rayleigh surface waves.This data set provides by itself a very good lateral resolution in the upper mantle, ranging from about 600 km in the shallow upper mantle to about 1200 km in the transition zone (Debayle et al. 2016), and can constrain the velocity structure down to around 1200 km depth (Durand et al. 2015).

Normal mode self-and cross-coupling coefficients
Surface wave data are not sufficient to constrain the whole mantle, so this data set is completed with spheroidal modes that are sensitive to the whole mantle.They consist of self-and cross-coupling structure coefficients.These coefficients are denoted c kk st , where k and k refer to the two modes that are coupled and s and t to the spherical harmonic degrees and orders.In the case of self-coupling k = k , while in the case of cross-coupling k = k .These self-coupling, respectively cross-coupling, structure coefficients represent the spherical harmonic coefficients of the splitting function, f kk (θ, φ), such that This function f kk (θ, φ) can be understood as a map of the frequency perturbations of a normal mode frequency due to the presence of some aspherical elastic heterogeneities inside the Earth.Self-coupling data are only sensitive to the even degrees of the mantle structure, while cross-coupling data can constrain either the even or the odd degrees.It is therefore important to combine both self-and cross-coupling structure coefficients.Fig. 2 shows examples of splitting functions built from measurements done by Deuss et al. (2013).The top row displays the observed splitting function for self-coupling of the two spheroidal modes, 9 S 10 and 6 S 15 , while the bottom row displays the splitting function of the cross-coupling of the same two spheroidal modes.The self-coupling of modes 9 S 10 and 6 S 15 only constrains the even degrees of the mantle structure, while the cross-coupling of these two modes is sensitive to the odd degrees.
For this study, we compile published structure coefficients taken from various studies (Smith & Masters 1989;Resovsky & Ritzwoller 1998a;Masters et al. 2000;Deuss et al. 2013;Koelemeijer et al. 2013).The data from Koelemeijer et al. (2013) include Stoneley modes characterized by a sensitivity confined near the CMB and therefore are important to constrain the structure of the deepest layers.Moreover SEISGLOB2, as well as our previous model SEISGLOB1, includes the recent cross-coupling coefficients catalogue of Deuss et al. (2013).This represents a significant increase compared to the nine pairs of modes included in the work of Resovsky & Ritzwoller (1998a), 19 yr ago, which before SEIS-GLOB1 was the only study that incorporated cross-coupling coefficients.We exclude all the modes that are sensitive to the inner  core, since it is known that they are characterized by an anomalous splitting only explained when inner core radial anisotropy is introduced (Morelli et al. 1986;Woodhouse et al. 1986;Tromp 1993).
We only include the structure coefficients up to degree 8 (s = 8, t ≤ s in eq. 1) because we consider that the measurements are not accurate enough above degree 8.It yields a total of 184 measured spheroidal modes, 158 for self-coupling structure coefficients and 26 for cross-coupling structure coefficients.A summary of the number of coefficients used in the inversion is presented in Table 2. Finally, each measurement is provided with an uncertainty that will be used in the inversion.We also provide in Tables S1 and S2 in the Supporting Information exhaustive lists of the mode self-and cross-coupling data, respectively, that have been used.

Body wave traveltimes
In order to refine the model in the mid-and lower mantle, we finally complete the data set with body wave traveltimes measured by Zaroli et al. (2010) at 34 s of period.The body wave data set contains 47 007 S, 42 174 SS and 11 480 ScS delay times that are given with respect to the radial model IASP91 (Kennett & Engdahl 1991) and that have been corrected from the effect of the crust using the model CRUST2.0(Bassin et al. 2000), the Earth's ellipticity and the anelastic effect from PREM model (Dziewonski & Anderson 1981), as detailed in Zaroli et al. (2010). Fig. 3(left) shows the body wave ray density from the surface down to the CMB.It clearly shows that this data set offers an ideal coverage down to 2250 km depth.

Crustal considerations
The strongest elastic heterogeneities are located in the crust (e.g.Ricard et al. 1996).There are two possibilities to take the crust into account.We can either choose an a priori crustal model and invert for the mantle structure only (Ishii & Tromp 2001;Ritsema et al. 2011) or we can invert for both crust and mantle structures (Resovsky & Ritzwoller 1998b;Durand et al. 2016).Since the body wave delay time measurements have already been corrected from the crust by Zaroli et al. (2010) using the CRUST2.0model (Bassin et al. 2000), we choose to correct all the data from crustal contributions using the same model.To do so, we did compute the frequency and phase velocity perturbations associated with CRUST2.0 for each mode and subtract this effect from our initial normal mode and surface wave data.The importance of such corrections is illustrated in Fig. 1 where we show phase velocity maps before and after crustal corrections.For surface waves, these corrections can be significant as shown for the fifth overtone at 50 s.

Forward problem
The chosen parametrization for the inversion could be either a grid that can be adapted to the data coverage (see for instance Zaroli et al. 2015) or a basis of smooth functions such as spherical harmonics and splines.Since both normal mode and surface wave data can be easily expressed in terms of spherical harmonics (Dahlen & Tromp 1998), we choose to use the same formalism for body wave data.It is also less computationally costly than using a grid of the same resolution.The seismic properties are also radially expanded using n = 21 spline functions, z n (r).In the following, we also choose to neglect the discontinuity topography effect (Ishii & Tromp 2001;Ritsema et al. 2011), except for the Moho which has been taken into account in the crustal corrections (see Section 2.4).Therefore, we can finally write where m n st = [δα/α, δβ/β, δρ/ρ] n st represents the perturbations of the seismic properties relative to a 1-D reference model, M(r ) the P velocity, S velocity and density sensitivity kernels (Woodhouse 1980) computed in the 1-D reference model and finally, c 0 , U 0 , ω 0 and β 0 are the phase velocity, group velocity, angular frequency and shear velocity in the reference 1-D model, respectively.The 1-D reference model is chosen to be PREM (Dziewonski & Anderson 1981) which implies to change the reference of the body wave delay times which were given with respect to IASP91 (see Section 2.3).PREM is a suitable reference model for long-period S waves as it simultaneously provides a good fit to normal mode central frequencies, surface and body waves.For this reason, PREM has been used as reference in a number of previous global tomographic studies (e.g.Ritsema et al. 2011;Chang et al. 2015).A drawback of PREM is the presence of a 220 km discontinuity which may not be a global feature.Kustowski et al. (2008) built a new 1-D reference Earth model, STW105, which is continuous at 220 km.However, since they do not include normal modes in the inversion, it is not clear that their model agrees with normal mode central frequencies.
On Fig. 3(right), we display examples of effective sensitivity for the various S waves used in the inversion.It is noteworthy that even though we use the ray theory to compute the body wave delay times, because our parametrization is finite (up to spherical harmonic degree 40 and 21 splines), we cannot resolve the rays as 'rays' so they appear larger as displayed on these panels.Our parametrization thus leads to significantly wider sensitivity kernels than the ray itself, which is similar to the width of finite-frequency 'banana doughnuts' (Dahlen et al. 2000), even though the origins of the sensitivity enlargements are different (physics of propagation for the latter, parametrization of the inversion for the former).This can partly explain why the results obtained with finite-frequency inversions until today give very similar results to those obtained with ray theory.examples of effective sensitivity kernels for each kind of body wave.We use the ray theory but because our parametrization is not infinitely refined once projected, the ray sensitivity is larger than the theoretical infinitely thin ray path.

Scaling seismic parameters
In order to reduce the computing time and because of a lack of resolution, a priori scaling factors between (β, α) and (β, ρ) are generally imposed, such that d ln (α) = ν α d ln (β) and d ln (ρ) = ν ρ d ln (β).This simplification was first applied by Li et al. (1991) who took ν α = 0.5 that they deduced from forward modeling of the splitting function of some normal modes and ν ρ = 0.25 according to Anderson et al. (1968).Later, theoretical calculations from Karato (1993) led to similar values for ν ρ .Resovsky & Ritzwoller (1998b) also impose scaling factors but used different values for the upper and lower mantle.In SEISGLOB2, we choose to impose depth independent scaling factors, taking ν α = 0.55 and ν ρ = 0.2 (Ishii & Tromp 2001;Durand et al. 2016).We have shown in Durand et al. (2016) that changing the scaling laws does not affect significantly the shear velocity model.

Inverse problem
Eqs (2)-(4) show that the forward problem is linear and can be written as d = Gm where d = [c kk st , σ k st , δt i ] is the data vector, and m = m n st the inverted parameter vector.The inverse problem is solved using a damped least-squares scheme (Tarantola & Valette 1982) and the final 'best-fitting' model m, hereafter called SEISGLOB2, can be expressed as follows: where m 0 is the a priori model, C p the a priori parameter covariance matrix, C d the data error matrix and 'T' denotes the transpose matrix.The matrix C p is chosen diagonal.A constant value, such as is the Kronecker symbol), would allow the high degrees to be as strong as the low degrees which is not what we observe from the amplitude spectrum of tomographic models.Therefore, we choose to impose a covariance that decreases with the harmonic degree such that C p (nst, n s t ) = σ 2 p s x δ(n − n )δ(s − s )δ(t − t ).We choose x = −0.1 and σ p = 3 × 10 −4 .This very slight a priori decrease of the amplitude with s is much smaller than what is found after the inversion which proves that our choice does not arbitrarily damp the short wavelengths or affect the final power spectrum of the heterogeneities.
The data covariance matrix C d is diagonal and filled with the data uncertainties.However, as the data (normal modes, surface waves and body waves) are of very different origins and in different numbers we rescale the uncertainties such that the final misfit for each subset is somewhat similar.There is certainly some arbitrariness in our choice (see Table 3).

SEISGLOB2 and data misfit
The model is presented Fig. 4 at various depths.For comparison, we also display SEISGLOB1 (Durand et al. 2016) in order to emphasize the gain in resolution provided by the use of a comprehensive body wave data set in SEISGLOB2.As expected, we discern very well in the upper mantle the imprint of surface structures such as cratons, oceanic ridges and some hotspots, like Hawaii where a circular low-velocity anomaly is clearly visible down to 670 km depth.As we go deeper into the transition zone, the subduction slabs become discernable and can be followed down to 1000 km depth for some of them, such as the Kermadec, Tonga, Java, South American and Mediterranean subductions.This is in complete agreement with the recent results of Zaroli (2016) who inverted the same S and SS body wave traveltimes, without the surface waves and normal mode measurements.At 1000 km depth, in agreement with previous models (e.g.Grand 1994;Van der Voo et al. 1999;Li et al. 2008) but even more strikingly, fast anomalies underline the closure of the Thethys (from under the Mediterranean sea to New Zealand) and the Farallon subduction (west of the Pacific-America border).Note that the maps between 670 and 1500 km contain shorter scale heterogeneities compared to the shallower and deeper maps.Finally, at the base of the mantle we find the LLSVPs below Pacific and Africa surrounded by a belt reminiscent of the positions of Mesozoic subductions (Ricard et al. 1993).
In Table 3, we present the initial and final misfits for each kind of data.Our scaling of the different data sets favours slightly the  surface and body waves which are better fitted by SEISGLOB2 than the normal modes.However, we show in Fig. 5 that SEISGLOB2 predicts accurately both surface waves and normal mode observations.This comforts us with the idea that our data weighting is a good compromise that captures the body wave information on in-termediate and short wavelengths, while remaining consistent with the long-wavelength structure constrained by the normal modes.Finally, even if we cannot entirely reject a dominance of even low degrees as we include more self-coupling than cross-coupling coefficients, we argue that the use of cross-coupling coefficients should , where s is the harmonic degree, x and A constants, which has been used to fit the spectrum of SEISGLOB2 (red), S40RTS (cyan) and SEMUCB-WM1 (blue).The red shaded area displays the zone ranging from 670 to 1500 km depth where the SEISGLOB2 spectrum is the flattest.make SEISGLOB2 more reliable than previous models for the low and odd degrees.
We also provide in Fig. S1 in the Supporting Information the correlations of SEISGLOB2 with other recent models [S40RTS (Ritsema et al. 2011), S362WMANI+M (Moulik & Ekström 2014), SEMUCB-WM1 (French & Romanowicz 2014)].SEISGLOB2 shows a good agreement with other seismic models, as good as among these models.Correlations are above the 95 per cent confidence level for degrees up to 20 in the shallow upper mantle.In the transition zone and lower mantle, the agreement between global models generally ranges between the 66 per cent and 95 per cent confidence level up to at least degree 10.The correlation increases for degrees lower than 10 at the base of the mantle.S40RTS appears to be the model closest to SEISGLOB2, which is expected since they share the same theoretical framework and used data of the same kind.However, there are many different steps in a global tomographic inversion and it is difficult to pinpoint exactly the origin of the discrepancies between two models using similar approaches.Moreover, even though they use similar data we include a larger self-coupling data set and cross-coupling measurements that are absent in S40RTS which can already lead to significant differences with S40RTS (Durand et al. 2016).In Fig. S2 in the Supporting Information, we show our final 1-D model, those obtained by various tomographic studies and PREM.The main differences between these models are located in the upper mantle.Our final 1-D model is very close to PREM and includes the 220 km discontinuity, while S362WMANI+M and SEMUCB-WM1, which do not use PREM as a reference model, are continuous at 220 km.

On the existence of a velocity spectrum pattern change around 1000 km depth
Fig. 6 displays the spectra computed for SEISGLOB1 (Durand et al. 2016), SEISGLOB2, S40RTS (Ritsema et al. 2011), S362WMANI+M (Moulik & Ekström 2014), and SEMUCB-WM1 (French & Romanowicz 2014).SEISGLOB1 and SEISGLOB2 share a parametrization in spherical harmonics and in splines, except that SEISGLOB1, only constrained by long-period surface waves and normal modes, is limited to degrees lower than 20.Fig. 6 shows all the information added mostly in the lower mantle thanks to the use of a large body wave data set in SEISGLOB2.Then, it shows that SEISGLOB2 has the highest amplitudes for larger harmonic degrees around 670-1200 km depth which suggests an accumulation of heterogeneities with horizontal wavelengths between about 4000 km (degree 10) and about 1300 km (degree 30).This feature is centred on 800-1200 km depth but extends from 670 to 1500 km depth.In order to quantify this observation, we then fit the spectrum of SEISGLOB2 and of the other models with a function of the form A s x , where s is the harmonic degree, x and A two constants.Then, we display in Fig. 6, the exponent x with depth.The spectrum of S362WMANI+M is too far from the function A s x to be displayed in Fig. 6.Moreover, it only contains degrees up to 20 and is thus less comparable to other models which contains degrees up to 40.A closer to zero exponent produces a flatter spectrum, meaning that high degrees have similar energy as low degrees.It is noteworthy that actually all spectra appear to be the richer in short-wavelength components around 1000 km depth.For SEISGLOB2, the flattest part is highlighted by a red shaded area which extends from 670 down to 1500 km depth, while for SEMUCB-WM1 it is from 670 down to 1100 km depth and for S40RTS it is around 800 km depth.It is noteworthy that we find an a posteriori exponent lower (x ranges from −0.2 to −0.6) than what was imposed in the a priori covariance matrix (x = −0.1)meaning that the flattening of SEISGLOB2 spectrum does not come from our a priori choices and is rather imposed by the data.
It can also be noted that SEISGLOB2 spectrum appears different from other models in the 250-660 km depth range.This difference is also seen in the correlations (Fig. S2, Supporting Information).This originates from our surface wave data set and this point is discussed in Debayle & Ricard (2012).These authors argued that it is due to the short radial correlation length required to fit the surface wave waveforms used to produce the phase velocity maps.
We also provide in Fig. S3 in the Supporting Information the radial correlation functions of SEISGLOB2 and of the same three tomographic models.It shows that for all the models there is a high degree of radial correlation through the lower mantle below 1000 km depth.Then, we observe that for SEISGLOB2 that there is a weak negative correlation (in blue) between S-wave velocities in the depth ranges 500-1000 km and 1200-2870 km.A weak or negative correlation layer is also observed in the radial correlation functions of S40RTS and SEMUCB-WM1.This suggest that in these models, the structures above and below 1000 km are decoupled.In S362WMANI+M, a weakly negative correlation is observed between the upper mantle and the structure below 1000 km.It is however difficult to estimate if this observation is required by the data, as the authors decouple the structure above and below the 670 km discontinuity a priori.All these observations suggest that there is a change in the global velocity pattern occurring around 1000 km depth in SEISGLOB2 but also in the other recent tomographic models.

The structure of the LLSVP and the behaviour of the slabs in the lower mantle
We now propose to check whether this observation of a change in the velocity spectrum pattern around 1000 km depth affects mantle structures such as slabs, LLSVPs or mantle plumes.We thus present in Figs 7 and 8 some cross-sections obtained in SEISGLOB2 and compare them with the other recent tomographic models (Ritsema et al. 2011;French & Romanowicz 2014;Moulik & Ekström 2014).
Figs 7 and 8 reveal that SEISGLOB2 and S40RTS agree very well which certainly comes from the fact that they are built in a similar way and they are based on the same kind of data (see also Fig. S1, Supporting Information).Even though SEMUCB-WM1 is based on a waveform inversion it also shares similar longwavelength patterns with SEISGLOB2.However, it also displays these vertical conduits across the whole mantle below hotspots reported in French & Romanowicz (2015) and that are absent in other models.S362WMANI+M is the most different and even the longwavelength structures are significantly different from the three other models.
In Fig. 7, we focus on the mantle below the Pacific.On these cross-sections, the Pacific LLSVP is very visible.The strongest low-velocity signature is located within the D layer but it is noteworthy that the LLSVP actually extends farther up in the lower mantle to almost 1000 km depth.In order to highlight this feature, we have overprinted on the profiles the contour line −0.1 per cent.On the profiles in the right-hand column, and conspicuously in SEISGLOB2 and S40RTS, we can also identify the Pacific slab subducting below the Tonga corresponding to a fast anomaly which extends laterally near 1000 km depth in SEISGLOB2.
On Fig. 8, we now focus on the mantle below Atlantic and Africa.We observe on these cross-sections the African LLSVP.It also extends to shallow depths, up to 1000 km depth.The cross-section in the middle column crosses the Afar and Tanzania hotspots.It has been suggested that the African hotspots are underlaid by a low-velocity anomaly in the upper mantle connected to the African LLSVP via a tilted low-velocity anomaly (Ritsema et al. 1999).However, SEISGLOB2, in agreement with S40RTS, suggests that the Afar and Tanzania hotspots are underlaid by two distinct lowvelocity anomalies in the upper mantle.The tilted low-velocity anomaly seems to only connect the Tanzania upper mantle to the south Atlantic CMB, while the low velocity beneath Afar extends down to about 900 km depth.However, although a low-velocity region is observed in the lower mantle beneath the Afar hotspot in all models, there is no clear evidence for a continuous vertical low-velocity anomaly linking the CMB to the surface.
Finally, SEISGLOB2 and S40RTS display a low-velocity anomaly beneath the Icelandic hotspot (right-hand column) crossing the upper mantle down to 1000 km depth.This anomaly is tilted toward the north in both models and seems to stagnate near 1000 km depth.In SEMUMCB-WM1, this anomaly is more deeply rooted in the lower mantle, still with a lower amplitude below 1000 km depth.However, the three models suggest that the Icelandic hotspot conduit is continuous at least down to 1000 km depth.
We also present in Fig. 9 tomographic cross-sections located on some subduction zones.In order to facilitate the identification of subducting slabs, we also represent some earthquakes by grey dots.For the majority of the slabs displayed in Fig. 9, they cross the transition zone and flatten around or below 1000 km depth.It is only for the Philippines and Mariana subductions that the slabs appear to be confined into the transition zone.Finally, it is noteworthy that on the cross-section across Java (Fig. 9, bottom), where we clearly see the slab crossing the transition zone and then partly extending laterally around 1000 km depth, Vanacore et al. (2006) have found P to S conversions reflecting the existence of a discontinuity around 1000 km depth.This is also in agreement with the slab imaging proposed by Fukao et al. (2001) showing the slab stagnating below the transition zone.
A common feature that we observe on all the cross-sections in Figs 7-9 is the existence of some structural changes around 1000 km.It is either where some subducting slabs spread laterally or where some mantle plumes seem to change behaviour and where the LLSVPs stop.These observations of a modification of the broad scale structure near 1000 km depth and of a flatter spectrum with shorter wavelengths heterogeneities in the depth range 670-1500 km, suggest that a significant change in mantle dynamics occurs in the upper half of the lower mantle rather than across the upper-lower mantle transition at 660 km.

Discussion
The existence of a lower mantle transition has already been proposed based on seismic observations.Its depth varies with the studies from 800 to 1200 km (Jenkins et al. 2016).The observation of such discontinuity in the lower mantle first comes from studies led in sub-duction zones (Kawakatsu & Niu 1994;Fukao et al. 2001;Deuss & Woodhouse 2002;Fukao & Obayashi 2013).It was first proposed to be the signature of the end of dissolution of garnets in the subduction slabs.However, this continuous transition is not expected to imply very significant velocity perturbations and should occur before 800 km deep (Ricard et al. 2005).Recently, Marquardt & Miyagi (2015) also show that the strength of ferropericlase significantly increases at pressures around 20-65 GPa, corresponding to 670-1000 km depth, which could cause a increase of viscosity and explain slab stagnation in the upper lower mantle.However, some studies also reported observations of this transition below continents and hotspots (Le Stunff et al. 1995;Deuss & Woodhouse 2002;Shen et al. 2003).Thus, the idea that this transition could be global has emerged.This study confirms that we can actually observe this global change in the spectrum of the recent global tomographic models, not only in SEISGLOB2.Now, if there is a global velocity pattern change around 1000-1200 km depth, we still have to understand its origin.It could be related to a continuous increase of viscosity with depth.Rudolph et al. (2015) indeed find that they best explain the observed geoid for viscosity profiles displaying an increase around 1000 km depth.The viscosity increase also enables them to explain the decrease in correlation for layers above and below this depth in SEMUCB-WM1.Since we saw that the radial correlation function of SEISGLOB2 is also characterized by a stronger decorrelation around 410-1000 km depth (see Fig. S3, Supporting Information), then it would also be a good explanation for our model.Such a viscosity increase would disrupt the slab geometry by inducing a thickening or a folding of slabs (Ricard et al. 1993;Ribe et al. 2007).This transition could also be the signature of a global and modest compositional change with depth.Indeed, Ballmer et al. (2015) show through numerical simulations that a moderate increase in MORB content in the lower mantle compared to the upper mantle (8 per cent) may be sufficient to produce this kind of stagnant slabs.Ballmer et al. (2015) then show that it is possible to produce and maintain such a chemical difference through the thermochemical convection of the mantle including the effect of the Clapeyron slopes of the olivine phase transitions.In their study, they also do not exclude the fact that this chemical change could be combined with a viscosity increase but they argue that a viscosity increase alone is not sufficient to explain slab stagnation.

C O N C L U D I N G R E M A R K S
In this paper, we present our new shear velocity tomographic model SEISGLOB2 based on surface wave, body wave and normal mode data.We observe in the spectrum of SEISGLOB2, a global change in the velocity pattern around 1000 km depth which extends from 800 down to 1500 km depth.In this depth range, the spectrum is flatter and contains a greater amount of heterogeneities with horizontal wavelengths between 1300 and 4000 km.We also confirm some changes in the structures of LLSVP, mantle plumes and subduction slabs around 1000 km depth.These observations strongly confirm previous ones on the existence of a global transition at around 1000 km depth.This observation suggests either the presence of a global chemical change or/ and a viscosity jump around 1000 km depth, but its precise origin is not yet known.This has a strong impact on the mantle geodynamics and implies that future numerical modeling should try to understand and reproduce this observation.Finally, SEISGLOB2 will be made freely available to the community.

A C K N O W L E D G E M E N T S
This work was supported by the French ANR SEISGLOB ANR-11-BLANC-SIMI5-6-016-01, the European Research Council within the framework of the SP2-Ideas Program ERC-2013-CoG, under ERC grant agreement number 617588 and the German Deutsch Forschungsgemeinschaft under the DFG project HAADES DU1634/1-1.SEISGLOB2 can be downloaded on the website of S. Durand: http://earth.uni-muenster.de/∼durand.

Figure 1 .
Figure 1.Examples of surface wave data for the fundamental mode at 40 s (left) and the fifth overtone at 50 s (right).Top row: phase velocity perturbation maps.The reference velocity is indicated on the top of the map.Middle row: phase velocity perturbation maps corrected from crustal effects.Bottom row: associated uncertainties normalized to 0.05 km s −1 .The green lines show the plate boundaries.The two curves represent the normalized shear wave velocity sensitivity kernels for each mode.

Figure 2 .
Figure 2. Examples of splitting functions built from the catalogue of Deuss et al. (2013).Top row: splitting functions for the self-coupling of the modes 9 S 10 and 6 S 15 .Bottom row: splitting function for the cross-coupling of these two modes 9 S 10 -6 S 15 .The colour scale represents the frequency perturbations in µHz with respect to those predicted by the reference model PREM.The normalized shear wave velocity kernels are also represented besides the maps.The spherical harmonic degrees constrained by these modes are indicated below each map.

Figure 3 .
Figure 3. Left: ray density (i.e.number of rays per area of 4• × 4 • ) at various depths intervals obtained from the body wave data(Zaroli et al. 2010).Right: examples of effective sensitivity kernels for each kind of body wave.We use the ray theory but because our parametrization is not infinitely refined once projected, the ray sensitivity is larger than the theoretical infinitely thin ray path.

Figure 4 .
Figure 4. Shear velocity anomalies of SEISGLOB1 (left) and of SEISGLOB2 (right) at various depth with respect to PREM (Dziewonski & Anderson 1981).The green lines show the plate boundaries and black circles the principal hotspots.

Figure 5 .
Figure 5.Comparison of some normal mode self-and cross-coupling splitting functions with those predicted by SEISGLOB2 and of some observed phase velocity maps with those predicted by SEISGLOB2.

Figure 6 .
Figure6.Five left-hand panels: logarithm of the amplitude spectra of SEISGLOB1, SEISGLOB2 and of three recent tomographic models.Right-hand panel: exponent of the function As x , where s is the harmonic degree, x and A constants, which has been used to fit the spectrum of SEISGLOB2 (red), S40RTS (cyan) and SEMUCB-WM1 (blue).The red shaded area displays the zone ranging from 670 to 1500 km depth where the SEISGLOB2 spectrum is the flattest.

Figure 7 .
Figure 7.Comparison of cross-sections below Pacific between SEISGLOB2 (top row) and three recent tomographic models: S40RTS (second row), SEMUCB-WM1 (third row) and S362WMANI+M (fourth row).We have overprinted the contour line −0.1 per cent and the depths lines at 410, 670 and 1000 km depth (dashed black lines).The maps displays where the cross-sections have been done and the blue and red circles enable to locate along the cross-section.The pink circles represent the hotspots and the green curves represent the plate boundaries.

Figure 8 .
Figure 8. Same as Fig. 7, but below Atlantic and Africa.

Figure 9 .
Figure 9. Example of slab imaging with SEISGLOB2.In order to highlight the slabs, we have overprinted some earthquakes in grey dots and the contour lines −0.1 per cent and +0.5 per cent.Most of them stop around 1000-1200 km depth except for the Mariana subduction where the slab is more confined into the transition zone.

Table 1 .
Table of the number of phase velocity measurements and finally the number of inversed coefficients.

Table 2 .
Table of the spheroidal mode structure coefficients used in this paper, and the corresponding studies.

Table 3 .
Table of initial and final misfits and of the variance reduction percentage as function of the data type.