Density structure of Earth's lowermost mantle from Stoneley mode splitting observations

Advances in our understanding of Earth's thermal evolution and the style of mantle convection rely on robust seismological constraints on lateral variations of density. The large-low-shear-wave velocity provinces (LLSVPs) atop the core–mantle boundary beneath Africa and the Pacific are the largest structures in the lower mantle, and hence severely affect the convective flow. Here, we show that anomalous splitting of Stoneley modes, a unique class of free oscillations that are perturbed primarily by velocity and density variations at the core–mantle boundary, is explained best when the overall density of the LLSVPs is lower than the surrounding mantle. The resolved density variations can be explained by the presence of post-perovskite, chemical heterogeneity or a combination of the two. Although we cannot rule out the presence of a ∼100-km-thick denser-than-average basal structure, our results support the hypothesis that LLSVPs signify large-scale mantle upwelling in two antipodal regions of the mantle.

a-c, 2-parameter search without CMB topography variations for a, All modes of the KDR13 and DRH13 data sets. b, All modes of the old HT96 and RR98 data sets, which were used in previous studies. c, All modes of the HT96 and RR98 data sets, using updated measurements from DRH13. d-f, Similar as a-c but for the 3-parameter search including CMB topography variations. As in Fig. 3 in the main text, the probability of density models is shown for different values of H, the scaling factor between lower mantle density variations and CMB topography, which is indicated on the right of each row. Negative values of H correspond to dynamically feasible models. Figure 8: Probability and L2-norm changes for lowermost mantle modes. a-b, Probability values and c-d L2-norm changes compared to SP12RTS for modes sensitive to shear-wave, compressional-wave and density structure in the lowermost mantle. Probability values and L2-norm changes are shown for the two best fitting models of Table 1, which feature light LLSVPs (solid blue and red lines in a, c). We also present the fit of the same models when the sign of RLL is reversed (dashed blue and red lines in b, d).

Supplementary Note 1: Splitting function amplitudes
Splitting functions are measured using an iterated damped least squares inversion of free oscillation spectra. To robustly constrain the density scaling factor R in the lowermost mantle, it is crucial to have accurate estimates of seismic anomaly amplitudes, in addition to their shape 7 . Therefore, much care has been taken in the development of the splitting function data set, to ensure that the splitting function amplitudes are independent of the damping value of the inversion 8 . The norm damping value has been varied by several orders of magnitude and we have ensured to pick damping values for which the results do not change significantly when lowering the damping further.

Supplementary Note 2: Splitting function predictions for other Stoneley modes
Supplementary Fig. 4 shows observed and predicted splitting function maps for additional Stoneley modes, similar to Fig. 2 in the main text. Again, predictions for mantle model SP12RTS ( Supplementary Fig. 4c, 4h and 4m) underestimate the splitting function amplitudes, which are reduced further by including dense LLSVPs ( Supplementary Fig. 4d, 4i and 4n). A better fit to the observed splitting functions ( Supplementary Fig. 4b, 4g and 4l) is only obtained for a positive value of RLL, i.e. light LLSVPs ( Supplementary Fig. 4e, 4j and 4o).  Fig. 5f, bottom) improve the fit compared to the default density model ( Supplementary Fig. 5e, centre). 2S16 and other Stoneley modes, show a strong trade-off between the CMB topography and LLSVP density scaling factor, making it difficult to constrain both separately.  Supplementary Fig. 6. This density model shows ~0.9 % lower densities for the LLSVPs and ~1.1 % higher densities for the surrounding regions in the "Ring around the Pacific". It should be noted that this model is not necessarily unique as a range of models fit the splitting function data within their uncertainties, but all feature low density LLSVPs.

Supplementary Note 5: Comparison with the data of previous studies
Our analysis, suggesting light LLSVPs, contradicts previous normal mode density models that obtained higher densities for the LLSVPs 9-12 . However, these models were constructed with fewer data and without the new Stoneley mode measurements. It has been suggested that significant trade-offs between upper and lower mantle density structure existed in these studies and that a range of models with both positive and negative scaling factors fit the data equally well [13][14][15] . At the same time, it was suggested that reliable density models may be obtained in future when more normal mode data are available. We repeat our twoparameter model space search for the subset of modes available in the HT96+RR98 data set (Supplementary Table 1) using both the original and updated splitting measurements, reported in Supplementary Fig. 7. For the original HT96+RR98 measurements 5-6 , we indeed find best fitting models with negative LLSVP density scaling factors, i.e. RLL = −0.7 to −1.3 ( Supplementary Fig. 7b). Using updated measurements for these modes from the DRH13 data set, we also find negative scaling factors for the LLSVPs: RLL = −2.4 to −3.1 ( Supplementary Fig. 7c). On the other hand, the Stoneley modes (Fig. 3c in the main text) find positive scaling factors for the LLSVPs with values between 1.7 and 1.9. In contrast, a model space search using all available modes of the DRH13+KDR13 data set [3][4] ( Supplementary Fig. 7a) produces a very broad region of best fitting models, with both positive and negative values of RLL. In the three-parameter model space search, where we include CMB topography variations ( Supplementary Fig. 7d-f), we observe similar patterns as in Supplementary Fig. 7a-c. For every value of H, both positive and negative values of RLL are found for all modes of the DRH13+KDR13 data sets (Supplementary Fig. 7d. The original measurements of the HT96+RR98 data sets (Supplementary Fig. 7e) also have difficulty resolving the sign of RLL for different values of H. However, the updated measurements for the HT96+RR98 data sets (Supplementary Fig. 7f) show a similar pattern to the Stoneley modes ( Fig. 3f in the main text), except that the shift to the different model class happens for negative values of H already, resulting in dynamically feasible models with dense LLSVPs.
We believe that by focusing on the Stoneley modes, we are able to extract the signal originating from the deep mantle. This philosophy is similar to the approach in body wave studies of using core-diffracted and reflected waves and not upper mantle phases to study the CMB region. Supplementary Fig. 7 also illustrates that using the HT96+RR98 selection of normal modes masks the signal due to the lowermost mantle, as the variations in probability values are very small (Supplementary Fig. 7b and c). In addition, none of the models in the entire model space search pass our observability criterion 16 due to the large uncertainties in the original splitting function measurements. Thus, we confirm the previous findings of dense LLSVPs using the HT96+RR98 data set. However, we also agree with subsequent studies that these data were insufficient to constrain the density structure of the lowermost mantle robustly. In contrary, we state that our Stoneley mode splitting data with unique sensitivity to the lowermost mantle are best fitted by light LLSVPs.

Supplementary Note 6: Fit to other lower mantle sensitive modes
To investigate whether light LLSVPs (as preferred by the Stoneley modes) are incompatible with the splitting functions of other modes sensitive to the lowermost mantle, we examine the fit of these modes to the best fitting density models listed in Supplementary Table 2. We define a mode to be significantly sensitive to the lowermost mantle if the sensitivity to structure below 2500 km depth is at least half the maximum sensitivity in the mantle.
Supplementary Fig. 8 shows the probability values and L2-norm changes with respect to the fit for SP12RTS for the resulting 26 modes. Most non-Stoneley modes sensitive to the lower mantle show small changes in the L2-norm ( Supplementary Fig. 8c-d), with only modes 0S2, 0S7 and 14S9 strongly preferring denser LLSVPs (dashed lines). However, the Stoneley modes have much higher L2-norm values for dense LLSVPs (Supplementary Fig. 8d), whereas they decrease very pronouncedly for light LLSVPs (Supplementary Fig. 8c). This discrepancy is likely to be due to unmodelled structure in the mid mantle. Other lower mantle sensitive modes have broader sensitivity kernels in the lower mantle and are therefore affected by structure in the mid mantle that is unaccounted for. On the other hand, the Stoneley modes are limited in sensitivity to depths near the CMB and will not be affected by unmodelled structures at mid mantle depths. Additionally, our inferences of light LLSVPs are based on the range of possible models that fit the splitting functions within their uncertainties, not on an L2-norm. Inspection of the probability values of these 26 modes indicates that modes such as 0S2 (which has been used to infer dense LLSVPs in recent studies 17 ) cannot be used to distinguish between light and dense LLSVPs, as the probability values are the same. 0S2 is the longest period normal mode, and consequently very difficult to measure, with only 78 usable spectra (8 earthquakes) available for the splitting function measurement 18 . The associated uncertainties in the splitting function coefficients are therefore very large, and almost any density model can be fitted within these uncertainties. As comparison, Stoneley modes are typically measured using 2000-3000 spectra from 90 earthquakes, giving rise to much smaller uncertainties 3 . Furthermore, the average probability of these 26 modes is higher for models with light LLSVPs (0.38 and 0.36 in Supplementary Fig. 8a) than for models with dense LLSVPs (0.25 and 0.25 in Supplementary Fig. 8b) or for SP12RTS (0.26).

Supplementary Note 7: Effect of default density scaling factor
The inversion for shear-and compressional-wave velocity variations in SP12RTS assumes R  Fig. 9b).

Supplementary Note 8: Testing trade-offs with shear-wave velocity structure
To test possible trade-offs between velocity and density structure for the Stoneley modes, we have repeated our model space search for two different models of velocity structure.
Model SP12RTS_low_D is a lower damped version of SP12RTS and hence features larger velocity amplitudes, about 1.4 times as large as in SP12RTS (Supplementary Fig. 10). This model allows us to investigate whether dense LLSVPs can be obtained when the shear-wave velocity amplitudes are larger, as the sensitivity kernels of the Stoneley modes for shearwave velocity and density are similar in the lowermost mantle ( Supplementary Fig. 3).
Supplementary Fig. 11 to 13 are equivalent to Fig. 2 to 4 in the main text, but we now use SP12RTS_low_D to describe the velocity structure in the entire mantle. The predicted splitting function maps (compare Supplementary Fig. 11 and Fig. 2  In Supplementary Fig. 15 to 18 we show the same results while using SP12RTS_P_scaled instead of SP12RTS. The splitting function maps ( Supplementary Fig. 14) now show slightly lower amplitudes, but the differences with SP12RTS are smaller than when the shear-wave velocity structure was changed ( Supplementary Fig. 11). Supplementary Fig. 15 (compare to

Supplementary Note 10: Velocity-density trade-offs
As Supplementary Fig. 11 to 16 demonstrate, we consistently find positive values of RLL with the amplitude depending primarily on the CMB topography scaling factor. Using model SP12RTS_low_D has a larger effect than using model SP12RTS_P_scaled, confirming the notion that the Stoneley modes respond more to variations in shear-wave velocity than Pwave velocity. Naturally, a full model space search of both velocity and density with all uncertainties taken into account is preferable over these experiments. However, they indicate that low densities are strongly preferred by the Stoneley mode data for reasonable models of lowermost mantle velocity structure.

Supplementary Note 11: Presence of a lower dense layer
We have fixed the depth, below which we vary the density scaling factors, to 2500 km depth, corresponding approximately to the anomalous region found in recent lower mantle tomographic clustering analyses 9 . Similar density models are obtained for mode 2S16 when we change this depth (RLL between 3 and 4 and RSR between 1 and 2), with the best fit to the splitting function measurements found for depths between 2400 and 2600 km (Supplementary Table 3