# An analysis of core-mantle boundary Stoneley mode sensitivity and sources of uncertainty

Alexander Robson, Harriet C.P. Lau, Paula Koelemeijer, Barbara Romanowicz

► To cite this version:

Alexander Robson, Harriet C.P. Lau, Paula Koelemeijer, Barbara Romanowicz. An analysis of core-mantle boundary Stoneley mode sensitivity and sources of uncertainty. Geophysical Journal International, 2022, 228, pp.1962-1974. 10.1093/gji/gab448. insu-03643035

HAL Id: insu-03643035

https://insu.hal.science/insu-03643035

Submitted on 20 Apr 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# An analysis of core–mantle boundary Stoneley mode sensitivity and sources of uncertainty

Alexander Robson,<fnref n="1" /> Harriet C.P. Lau,<fnref n="1" />,<fnref n="2" /> Paula Koelemeijer,<fnref n="1" />,<fnref n="3" /> and
Barbara Romanowicz<fnref n="1" />,<fnref n="3" />,<fnref n="4" />

Accepted 2021 October 27. Received 2021 August 4; in original form 2021 October 15

## SUMMARY

Stoneley modes are a special subset of normal modes whose energy is confined along the core–mantle boundary (CMB). As such, they offer a unique glimpse into Earth structure at the base of the mantle. They are often observed through coupling with mantle modes due to rotation, ellipticity and lateral heterogeneity, though they can be detected without such coupling. In this study, we explore the relative sensitivities of seismic spectra of two low-frequency Stoneley modes to several factors, taking as reference the fully coupled computation up to 3 mHz in model S20RTS. The factors considered are (i) theoretical, by exploring the extent to which various coupling approximations can accurately reproduce reference spectra and (ii) model-based, by exploring how various Earth parameters such as CMB topography, attenuation and S- and P-wave structures, and the seismic source solution may influence the spectra. We find that mode-pair coupling is insufficiently accurate, but coupling modes within a range of ±0.1 mHz produces acceptable spectra, compared to full coupling. This has important implications for splitting function measurements, which are computed under the assumption of isolated modes or at best, mode-pair or group coupling. We find that uncertainties in the P-wave velocity mantle model dominate compared to other model parameters. In addition, we also test several hypothetical models of mantle density structure against real data. These tests indicate that, with the low-frequency Stoneley mode spectral data considered here, it is difficult to make any firm statement on whether the large-low-shear-velocity-provinces are denser or lighter than their surroundings. We conclude that better constraints on long wavelength elastic mantle structure, particularly P-wave velocity, need to be obtained, before making further statements on deep mantle density heterogeneity. In particular, a dense anomaly confined to a thin layer at the base of the mantle (less than ~100–200 km) may not be resolvable using the two Stoneley modes tested here, while the ability of higher frequency Stoneley modes to resolve it requires further investigations.

**Key words:** Composition and structure of the core; Composition and structure of the mantle; Structure of the Earth; Surface waves and free oscillation; Theoretical seismology.

## 1 INTRODUCTION

Within the seismic spectrum, normal modes provide the strongest constraints on density structure within the Earth (Dahlen & Tromp 1998). Knowledge of density is important in determining the nature of mantle dynamics since buoyancy dictates the flow of mantle rock. Furthermore, the relationships between elastic parameters and density helps to identify the thermochemical nature of mantle heterogeneities (Karato 1993).

Seismic studies of density have focused on constraining the structure of the lowermost mantle, where two large, equatorial and antipodal low shear wave-speed structures beneath Africa and the central Pacific Ocean were first observed in seismic tomography over 30 yr ago (Dziewonski 1984). The existence of these two ‘large-low-shear-velocity-provinces’ (LLSVPs) has been demonstrated to be robust, consistently observed using a variety of data sets and methodologies (e.g. Lekic et al. 2012). With a shear wave-speed ($v_s$) reduction reaching 2–3 per cent, at lateral scales of thousands

<footnote n="1">¹Department of Earth and Planetary Science, University of California, Berkeley, CA 94720, USA. E-mail: hcplau@berkeley.edu</footnote>

<footnote n="2">²Department of Earth Sciences, Royal Holloway, University of London, Egham, Surrey TW20 0EX, UK</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris, 1 Rue Jussieu, 75005 Paris, France</footnote>

<footnote n="4">⁴College de France, 11 Pl. Marcelin Berthelot, 75231 Paris, France</footnote>
of kilometres, these LLSVPs cover roughly one quarter of the core–mantle boundary’s (CMB) surface area, and in some models, extend ~1000 km above the CMB. However, their origin, detailed morphology, stability over time, composition and role in mantle circulation, remain open questions (e.g. Forte & Mitrovica 2001; Mulyukova et al. 2015; Zhong & Rudolph 2015; Davaille & Romanowicz 2020). Their buoyancy structure is a topic of continued debate (e.g. Ishii & Tromp 1999; Kuo & Romanowicz 2002; Moulik & Ekström 2016; Koelemeijer et al. 2017; Lau et al. 2017). Some studies favour seismic wave-speed reductions as a result of increased temperatures alone, while others suggest the presence of compositionally distinct material within them (see the review by, e.g. McNamara 2019).

Resolving LLSVP density via seismological observations remains difficult. With the exception of reflection at impedance contrasts, body waves offer little sensitivity, making longer-period normal modes a more promising avenue to resolve the density structure of the deep mantle, at least at long wavelengths. Over 20 years ago, Ishii & Tromp (1999) first used this sensitivity in combination with free-air gravity constraints, and inverted for seismic wave speeds and density, proposing a density increase of a few percent within the LLSVPs. Shortly after, it was shown that the result of this inversion was highly non-unique (Masters et al. 2000; Romanowicz 2001) and heavily dependent on the wave-speed structure of the initial model (Kuo & Romanowicz 2002). Subsequent studies utilized probabilistic forward modelling approaches, again finding that normal mode observations favoured an increase in density within the LLSVPs, albeit the models were parametrized with thick depth layers due to computational constraints (Resovsky & Trampert 2003; Trampert et al. 2004). More recently, a variety of tomographic models generated to fit seismic observations and/or geodynamic constraints have found LLSVPs with increased densities, though the height, extent and magnitude of the denser regions vary across models (Simmons et al. 2010; Mosca et al. 2012; Moulik & Ekström 2016).

Meanwhile, several geodynamic studies have disputed the view of LLSVPs as dense, compositionally distinct and stable structures. For example, Davies et al. (2012) generated a variety of thermal and thermochemical mantle convection models and concluded that purely thermal, buoyant, LLSVP models were sufficient to explain teleseismic S-wave traveltimes, while thermocompositional LLSVPs regularly overpredicted deep mantle shear wave-speed anomalies. In contrast, Jones et al. (2020) argued that material with chemical density contrasts of at least 3 per cent must exist to match the amplitudes of S-wave tomography model S40RTS. More recently, Davaille & Romanowicz (2020) explored the effects of compositional components within buoyant upwellings and found that it is possible to maintain overall buoyancy even with the inclusion of compositionally dense components.

Recently, Lau et al. (2017) used GPS measurements of Earth’s body tides for the first time, and performed a model-space search for a $v_s$-to-density scaling factor within and outside the LLSVPs, finding that best fitting models required an excess of density within the LLSVPs. In a similar fashion, Koelemeijer et al. (2017) used a newly collected data set of CMB Stoneley modes (Koelemeijer et al. 2013), a group of modes whose sensitivity is concentrated in the lowermost mantle and the uppermost core (Fig. 1), and thus well-suited for investigating deep mantle structure. Their model-space search methodology was applied to explore $v_s$-to-density scaling factors within and outside the LLSVPs, as well as the amplitude and sign of CMB topography. Koelemeijer et al. (2017) reported two groups of equally well-fitting models: one group, favouring dense LLSVPs, and the other, favouring buoyant ones. The buoyant

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Stoneley mode radial sensitivity: sensitivity kernels for isotropic P- and S-wavespeeds ($v_p, v_s$) and density, $\rho$, as a function of depth, computed in the radially anisotropic Preliminary Reference Earth Model (PREM, Dziewonski & Anderson 1981) for Stoneley modes $1S_{11}$ and $1S_{12}$ to structure of spherical harmonic degree, $s=2$. Note the low sensitivity in the upper and mid-mantle and lack of sensitivity below the uppermost outer core. The sensitivity kernels are normalized to plot on the same scale.</figcaption></figure>

LLSVP models were preferred due to dynamic arguments about the sign of the associated underlying CMB topography.

Here we revisit the question of Stoneley mode sensitivity to density heterogeneity in the deep mantle using full coupling calculations. Given their low energy near the Earth’s surface, the observation of CMB Stoneley modes often (though not exclusively) requires coupling with mantle modes in the presence of the Earth’s rotation, ellipticity and laterally heterogeneous (3-D) structure. As a result, several factors may play an important role, such as the coupling theory assumed and uncertainties on mantle model parameters.

Thus far, Stoneley mode-based studies have applied the so-called ‘splitting function’ approach (a derived data type involving intermediate processing steps of spectra; Giardini et al. 1988), that assumes coupling between Stoneley modes and neighbouring spheroidal modes is only of first order (Koelemeijer et al. 2013). Such truncated normal mode perturbation theory may introduce significant errors in the determination of Earth structure (Yang & Tromp 2015; Akbarashrafi et al. 2018; Jagt & Deuss 2021), and its effects on Stoneley modes has yet to be explicitly explored. While we do not directly calculate splitting functions, our results may hold implications for Stoneley mode splitting function measurements.

In an ideal world, we would couple all modes to as high a frequency as required, but several levels of approximation exist to ensure the problem remains tractable: self-coupling, mode-pair coupling, narrow-band (or group) coupling and full coupling (see Dahlen & Tromp 1998; Deuss & Woodhouse 2001).

Due to computational constraints, we focus on two low-frequency Stoneley modes: $1S_{11}$ and $1S_{12}$. Specifically, we compute synthetic normal mode spectra, based on full coupling of modes with frequencies up to 3 mHz within the framework of first-order perturbation theory, to investigate the effect of various levels of truncation in the coupling theory on the spectra of the Stoneley-fundamental mode
<table><caption>**Table 1.** Parameters for data processing: Stoneley mode observed centre frequencies (Koelemeijer *et al.* 2013), frequency window considered for spectra comparisons and time windows considered for spectra computations. Time windows, chosen to maximize signal, were also taken from this study.</caption><thead><tr><th rowspan="2">Stoneley mode</th><th rowspan="2">Observed centre frequency (μHz)</th><th colspan="2">Bandpass filter (mHz)</th><th colspan="2">Window (hr after event)</th></tr><tr><th>Low</th><th>High</th><th>Start</th><th>End</th></tr></thead><tbody><tr><td><sup>1</sup>S<sub>11</sub></td><td>2345.64 ± 0.41</td><td>2.330</td><td>2.365</td><td>20</td><td>60</td></tr><tr><td><sup>1</sup>S<sub>12</sub></td><td>2552.55 ± 0.09</td><td>2.540</td><td>2.580</td><td>10</td><td>60</td></tr></tbody></table>

pairs $^{1}$S$_{11}$ − $^{0}$S$_{15}$ and $^{1}$S$_{12}$ − $^{0}$S$_{17}$ (Table 1). We only consider these two Stoneley modes since they are comfortably located in the middle of this frequency band. We should note that this limits the extrapolation of our results to higher-frequency Stoneley modes.

We also assess the magnitude of errors due to uncertainties in mantle structure (elastic, anelastic and CMB topography). Finally, we compare the magnitude of all these effects against the magnitude of the signal generated by conceptual models of LLSVP density distribution. We conclude that low-frequency Stoneley modes can indeed offer constraints on the density of LLSVPs, however (1) the use of mode-pair coupling theory introduces significant errors, often larger than the signal being studied, with important implications for the use of low-frequency Stoneley mode splitting functions and (2) in order to better resolve the strength and spatial extent of density anomalies in the deep mantle, improved models of mantle velocity structure must first be constructed.

## 2 METHODOLOGY

### 2.1 Theory

We generate synthetic spectra for Stoneley modes $^{1}$S$_{11}$ and $^{1}$S$_{12}$, applying the normal mode coupling code of Yang & Tromp (2015). This code allows full coupling of all modes within a specified frequency band, under the assumption of first order perturbation theory. We couple modes (with eigenfunction and eigenfrequency {$s'_{k}, \omega'_{k}$} associated with the eigenmode $k$) of the spherically symmetric, non-rotating, elastic and isotropic (SNREI) Earth model PREM (Dziewonski & Anderson 1981). Eigenmodes ($s'$) of a model perturbed from a SNREI earth model by rotation, anelasticity and lateral heterogeneity, are expressed as linear combinations of the SNREI eigenmodes as follows:

$$s' = \sum_{k} q_k s_k^0, \qquad (1)$$

where $q_k$ are weights that are dependent on the aspherical structure introduced by departing from the SNREI earth model. As described in Yang & Tromp (2015), these weights are determined by casting the perturbations away from the SNREI model as a (non-linear) eigenvalue problem, where the eigenvalues represent fully perturbed eigenfrequencies, $\omega'_{k}$ (where both $q_k$ and $\omega'_{k} \in \mathbb{C}$). Matrices in this problem are commonly known as splitting matrices, where asphericity breaks a $(2\ell + 1)$-degeneracy with $\ell$ representing the spherical harmonic degree for each mode $k$. For further details, see Yang & Tromp (2015).

The broad outline of this process is shown on the right-hand side of Fig. 2. The choice of asphericity (rotation, ellipticity and 3-D model) and coupling approximation (see Section 2.3.1) determines the entries in the splitting matrix. We solve the eigenvalue problem to determine the fully 3-D modes, {$s'_{k}, \omega'_{k}$}. These 3-D modes are used to form the Earth’s 3-D Green’s Function, which is convolved

with the earthquake source to produce the synthetic seismogram. This is subsequently processed following the same steps as applied to observed seismograms in order to obtain mode spectra in target frequency windows (see Section 2.2 and Fig. 2).

We will compare all our tests against a reference case that incorporates coupling through rotation, ellipticity and 3-D mantle structure using tomographic model S20RTS (Ritsema *et al.* 1999). This model only prescribes shear wave velocity variations and we assume the following scalings of $\partial \ln \rho / \partial \ln v_s = 0.3$ and $\partial \ln v_p / \partial \ln v_s = 0.5$ to describe the $v_p$ and $\rho$ structure (consistent with model S20RTS itself). The density scaling is compatible with a purely thermal hypothesis for the origin of mantle heterogeneity (e.g. Karato 1993). We assume no CMB topography and the 1-D attenuation model QL6 (Durek & Ekstrom 1996). For this reference case, we *fully couple* the modes up to 3.1 mHz (see Section 2.3.1). In each Section 2.3.1–2.3.6, we elaborate on specifics of the coupling level when different from this reference case.

All synthetic spectra were computed for the case of the 24 May 2013 Sea of Okhotsk ($M_w$ 8.3) event, a large and deep earthquake (~598 km depth) that, at the time of this study, was not incorporated in published splitting function data sets of Stoneley modes. For the reference case, we applied the Harvard CMT source solution (Ekström *et al.* 2012). In order to consider a realistic situation, we computed vertical component synthetic spectra for a collection of 55 stations for which high quality records are available for the 2013 Okhotsk earthquake (Fig. 3).

### 2.2 Data selection and processing

All observed seismograms available from Streckeisen STS-1 broadband seismometers for a period of 96 hr after the event were collected, individually inspected, and seismograms with glitches and/or excessive noise were removed. The records from the remaining 55 stations were kept (Fig. 3) and processed (Fig. 2, Table 1) for comparison to synthetics. These were obtained via the online Global Seismic Network databases (IU Network n.d.; IC Network n.d., DOIs noted in the reference list).

The left side of the flowchart (Fig. 2) describes how the raw seismic data were processed. Before the steps that are common to both observations and synthetics were performed, the instrument response was removed from the raw data. These steps included an initial detrending and application of a Hamming taper. Once the instrument response was deconvolved, the data were detided by removing a 2nd-order spline (see Beyreuther *et al.* 2010), resulting in a seismogram that can be directly compared to the corresponding synthetic seismogram. The open source processing toolkit *Obspy* (Beyreuther *et al.* 2010) was applied to perform these steps.

Both synthetics and observed seismograms were converted to spectra following the steps indicated by grey boxes in Fig. 2. The time windows considered for each Stoneley mode, before transforming to the frequency domain, consistent with the splitting function measurements by Koelemeijer *et al.* (2013), are listed in Table 1.

### 2.3 Details of synthetic tests

In the following Sections 2.3.1–2.3.6 we describe several synthetic tests performed to determine the sensitivity of the low-frequency Stoneley modes to various theoretical and model assumptions. We discuss their implications in Section 3.1.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Flow chart of method: the processing stages of seismic data (left-hand side) and the computation of synthetic spectra (right-hand side).</figcaption></figure>

### 2.3.1 Coupling approximations

Since the energy of Stoneley modes peaks around the CMB region (Fig. 1), coupling with mantle modes (whose sensitivity is concentrated within the mantle) is important to generate a strong surface observation of the Stoneley mode. To visualize the neighbourhood of modes that couple to $_1S_{11}$ and $_1S_{12}$ for each approximation of coupling, Fig. 4 shows spectra at two stations CASY and BILL, where the frequency brackets define the modes that are coupled for each approximation. We chose four levels of coupling: full (as in the reference case), both ±0.5 and ±0.1 mHz for narrow-band

coupling and mode-pair coupling. Full coupling up to 3.1 mHz includes coupling a total of 122 modes (including toroidal modes as well). As the width of the coupling band is narrowed, the number of modes coupled to $_1S_{11}$ and $_1S_{12}$ is reduced to 60 and 61 (respectively) for ±0.5 mHz, 14 and 12 modes for ±0.1 mHz and finally, two modes for mode-pair coupling. As an example, for station CASY (Fig. 4), we see that in the case of ±0.5 mHz coupling for $_1S_{11}$, the coupled spheroidal modes are modes whose SNREI eigenfrequencies lie in between and including those of $_0S_{11}$ and $_0S_{19}$. Once these selections of modes are determined, the splitting matrix can be constructed.
Sea of Okhotsk 2013

Mw 8.3, z = 598 km

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Source and receivers: location and mechanism of the 2013 Sea of Okhotsk deep earthquake (red and white beach ball) and distribution of very broad-band (VBB) stations considered in this study (yellow and red inverted triangles). Only the 55 yellow stations have been used in the comparison with real data.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Coupling approximations: synthetic full-coupling spectra showing different considered coupling levels marked for ¹S₁₁ at station CASY, Casey, Antarctica and ¹S₁₂ at station, BILL, Bilibino, Russia, computed for the 2013 Sea of Okhotsk (Mw 8.3) earthquake. Colour bands mark the extent of the different levels of coupling considered in this study, with modes included in mode-pair coupling in pink, modes included in ±0.1 mHz coupling in blue and modes included in ±0.5 mHz coupling in green. These colours are also reflected in the brackets above the spectra. Note that the list of labelled modes is non-exhaustive, only indicating some of the more prominent modes present.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Test models: (a) depth dependent scaling factors for $\partial \ln v_p/\partial \ln v_s$ (solid lines) and $\partial \ln \rho/\partial \ln v_p$ (dashed lines, which overlap for the two models) for mantle elastic models. (b) The $Q_\mu$ models – QL6 (Durek & Ekstrom 1996), Qm1 (Widmer et al. 1991), $Q_{\text{PREM}}$ (Dziewonski & Anderson 1981). In (a) and (b), the reference models are S20RTS and QL6. (c) the reference excess density field, $\partial \ln \rho$, based on S20RTS (Ritsema et al. 1999) is shown at the base of the mantle, with the maximum variations corresponding to the colourbar of $\epsilon = \pm 0.8$ per cent. (d) The test LLSVP density models, shown at the base of the mantle, with $\epsilon = \pm 2$ per cent, $\pm 10$ per cent for the ‘Dense Tall’ and ‘Dense Basal’ LLSVP models, respectively.</figcaption></figure>

<table><caption>**Table 2.** Source parameters for each of the three source solutions based on the Harvard CMT reference solution (Ekström et al. 2012), the United States Geological Survey solution (USGS, 2013) and solution E35 from Abubakirov et al. (2015).</caption><thead><tr><th></th><th>Half-duration<br>(s)</th><th>Latitude<br>(°)</th><th>Longitude<br>(°)</th><th>Depth<br>(km)</th><th><i>M</i><sub>rr</sub></th><th><i>M</i><sub>θθ</sub></th><th><i>M</i><sub>φφ</sub><br>(× 10<sup>27</sup> Nm)</th><th><i>M</i><sub>rθ</sub></th><th><i>M</i><sub>rφ</sub></th><th><i>M</i><sub>θφ</sub></th></tr></thead><tbody><tr><td>CMT</td><td>35.7</td><td>54.89</td><td>153.22</td><td>598.1</td><td>-16.7</td><td>3.82</td><td>12.8</td><td>-7.84</td><td>-35.7</td><td>1.55</td></tr><tr><td>USGS</td><td>35.7</td><td>54.89</td><td>153.22</td><td>598.1</td><td>-14.0</td><td>3.10</td><td>10.8</td><td>-7.90</td><td>-35.4</td><td>0.80</td></tr><tr><td>E35</td><td>16.0</td><td>54.89</td><td>153.22</td><td>640.0</td><td>-12.9</td><td>4.10</td><td>8.70</td><td>-11.3</td><td>-37.7</td><td>0.60</td></tr></tbody></table>

### 2.3.2 Elastic mantle structure

To test the effect of uncertainties in mantle elastic structure on Stoneley mode spectra, we also calculated synthetic spectra adopting SP12RTS (Koelemeijer et al. 2016) as an alternative velocity model. Differences between SP12RTS and the reference model (S20RTS) are small in $\partial \ln v_s$ (Koelemeijer et al. 2016) due to similarity in the data and methods used in model construction, though the reference model S20RTS does not include the augmented data set of SP12RTS. However, in $\partial \ln v_p$ these differences are more pronounced as, unlike S20RTS, $\partial \ln v_p$ was inverted for independently in the construction of SP12RTS. As a result, in addition to exploring the effects of uncertain $\partial \ln v_s$ structure, this comparison simultaneously addresses errors due to assuming fixed (and usually constant) $\partial \ln v_p/\partial \ln v_s$ scaling, which is common in tomographic models. For comparison, Fig. 5(a) shows the scaling factors $\partial \ln v_p/\partial \ln v_s$ and $\partial \ln \rho/\partial \ln v_s$ as a function of depth for both reference model S20RTS and for SP12RTS.

### 2.3.3 Anelastic mantle structure

Anelasticity introduces two effects: dispersion results in the frequency dependence of the elastic moduli [where $v_s \rightarrow v_s(\omega)$, $v_p \rightarrow v_p(\omega)$], and attenuation adds an imaginary part to a mode's eigenfrequency (whose real part is also perturbed due to dispersion). For our reference anelastic structure, we assume the attenuation, $Q^{-1}$, model ‘QL6’ (Durek & Ekstrom 1996). To explore Stoneley mode spectra sensitivity to $Q^{-1}$ structure, we also produce spectra

assuming QM1 (Widmer et al. 1991) and the Q structure of PREM (Dziewonski & Anderson 1981). Together, these represent the most commonly adopted radial attenuation models, reflecting the range of variability in long period-derived Q models, all of which are broadly consistent (e.g. Cammarano & Romanowicz 2008; Romanowicz & Mitchell 2015). Fig. 5(b) shows each model for shear attenuation ($Q_\mu$).

### 2.3.4 LLSVP density structure

The reference model (S20RTS) assumes a purely thermal origin of heterogeneity, as reflected in the positive scaling of $\partial \ln \rho/\partial \ln v_s$. Under this assumption, the LLSVPs represent regions of thermally driven upwellings of lower density than their surroundings. The density field at the base of the mantle in this model is shown in Fig. 5(c). To test the sensitivity of low-frequency Stoneley mode spectra to the density distribution within the LLSVP regions, we explored two additional conceptual models in which the LLSVPs are denser than average: (1) ‘Dense Tall’ and (2) ‘Dense Basal’. The former model is that found in Lau et al. (2017). In that study, a whole distribution of models were found to fit body tide data and here we choose to use the mean of the best fitting models. This model features mild excess density within the LLSVPs (defined by the -0.65 per cent shear wave-speed anomaly contour as in Torsvik et al. 2006) extending up to 600 km above the CMB. Outside this region, we applied a scaling of $\partial \ln \rho/\partial \ln v_s$ such that the overall mass of the Earth is conserved. For the ‘Dense Tall’ and ‘Dense Basal’ models the peak surrounding buoyancies at the base of the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Box and whisker diagrams of the distribution of synthetic residuals associated with each synthetic test case, arranged in the same order as in Fig. 7.</figcaption></figure>

<table><caption>**Table 3.** Median and mean residuals for each model of the quantities defined in eq. (2). The top 11 entries are relative to the synthetic reference model (i.e. the synthetic residual), arranged in the same order as in Fig. 7, whereas the bottom four entries are relative to real data (i.e. the data residual).</caption><thead><tr><th></th><th>Varied<br>Property</th><th>Median<br>Residual<br>(per cent)</th><th>Mean<br>Residual<br>(per cent)</th></tr></thead><tbody><tr><td rowspan="11">Synthetic</td><td>±3 km CMB topography</td><td>1.99</td><td>3.27</td></tr><tr><td>Q-PREM attenuation Model</td><td>2.52</td><td>4.45</td></tr><tr><td>QM1 attenuation Model</td><td>1.95</td><td>3.72</td></tr><tr><td>E35 source solution</td><td>21.31</td><td>43.29</td></tr><tr><td>USGS source solution</td><td>7.17</td><td>10.53</td></tr><tr><td>'Dense Basal' LLSVP</td><td>61.57</td><td>233.11</td></tr><tr><td>'Dense Tall' LLSVP</td><td>59.46</td><td>200.00</td></tr><tr><td>SP12RTS velocity model</td><td>21.59</td><td>74.09</td></tr><tr><td>Pair coupling</td><td>59.43</td><td>337.64</td></tr><tr><td>±0.1 mHz band coupling</td><td>1.99</td><td>3.27</td></tr><tr><td>±0.5 mHz band coupling</td><td>1.50</td><td>3.00</td></tr><tr><td rowspan="4">Data</td><td>'Dense Tall' LLSVP</td><td>59.46</td><td>200.00</td></tr><tr><td>'Dense Basal' LLSVP</td><td>61.57</td><td>233.11</td></tr><tr><td>S20RTS</td><td>61.24</td><td>227.12</td></tr><tr><td>SP12RTS</td><td>61.50</td><td>207.89</td></tr></tbody></table>

mantle reach ~-1 and ~-2 per cent, respectively. The latter model is based on the total excess mass within the LLSVPs of the 'Dense Tall' model, with this excess mass confined to a 100-km-thick basal layer within the LLSVPs (e.g. Romanowicz 2017). Surrounding the LLSVPs in this depth range, we again choose $\partial\ln \rho/\partial\ln v_s$ such that total mass is conserved. For both models, above the depths where material is present, the models are identical to the reference case. The density structure just above the CMB of these two models is displayed in Fig. 5(d).

### 2.3.5 CMB topography

With little consensus on CMB topography (Koelemeijer 2021), we designed a hypothetical model to explore the sensitivity of Stoneley mode spectra to CMB topography. We required that the spatial pattern of the topography matches that of density at the base of the mantle in the reference model (S20RTS), using a linear scaling factor. This assumes that a buoyant region of lowermost mantle material will correspond to an upwelling, resulting in an upward deflection of the CMB and that the reverse is true for dense regions. The scaling is defined such that the largest absolute amplitude of
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Example Stoneley mode spectra and residuals: a comparison of test cases to the synthetic full-coupling reference (black line) for 1S11 (at station CASY, Casey, Antarctica; left-hand panel) and 1S12 (at station BILL Bilibino, Russia; right-hand panel). Absolute values of the residual between each test and the reference are shown filled by different colours, with a label and a percentage describing the integral under the absolute residual relative to that of the reference. Note that the residual for the ‘Dense Basal’ model is small and we return to the significance of this in Section 3.2.</figcaption></figure>

topography is ±3 km. This is approximately consistent with observations of CMB topography at long wavelengths (Koelemeijer 2021).

### 2.3.6 Source parameters

In addition to the reference Harvard CMT source solution (Ekström et al. 2012), we considered source parameters as given in the US Geological Survey catalogue (hereafter, ‘USGS’; USGS, 2013) as well as those of Abubakirov et al. (2015) (hereafter ‘E35’). The reference (Harvard CMT) and USGS solutions are broadly similar in their construction, both utilizing teleseismic body waves alongside intermediate to long-period surface waves, to compute a centroid moment tensor. In contrast, the E35 solution was obtained using only regional waveforms (epicentral distance 8–25°) between 125 and 250 s. The parameters for each solution are listed in Table 2. Differences between the solutions are small in terms of the radiation pattern, but larger in terms of depth and source duration.

# 3 RESULTS AND DISCUSSION

## 3.1 Insights from synthetic tests

To quantify how the results of each synthetic test compare to our reference case (full coupled synthetics up to 3.1 mHz), we calculated the residuals between each test case (syn) and the reference (ref) amplitude spectra across each station and mode pairs, where $Y^i$ is

the associated amplitude value. The $i$th residual is then defined as

$$ \text{residual}_i = \left| \frac{Y_{ref}^i - Y_{syn}^i}{Y_{ref}^i} \right|, \qquad (2) $$

where $i$ encapsulates a given mode-station pair at the sampled frequency for each mode frequency band (Table 1) and $||X||$ denotes taking the absolute value of quantity $X$. This results in a distribution of $N$ residuals for each test case, where $N = (55 \times N_1) + (55 \times N_2)$, accounting for 55 stations and the number of sampled frequency points within each Stoneley mode frequency band, $N_1$ and $N_2$.

These values are summarized in Fig. 6, with Table 3 listing the median and mean residuals for each model over the entire data set (first 11 entries). We emphasize that these residuals do not yet assess any fit to real data, which we reserve for Section 3.2. These residuals are also shown in Fig. 7 as a function of frequency for two stations. From these summary values and figures, it is clear that mode-pair coupling, by a significant margin, produces the largest residuals from our reference case of full coupling for S20RTS.

Widening the coupling frequency range reduces the residual, as is expected. Restricting the coupling between modes to within ±0.5 mHz of the target modes generally introduces negligible error in the context of other uncertainties (Figs 6 and 7). Even the narrow-band coupling truncated to within ±0.1 mHz of the target mode (which includes only 14 and 12 modes for 1S11 and 1S12, respectively) results in relatively accurate spectra (Fig. 7).

Aside from mode-pair coupling, the largest source of discrepancy among those explored here is the choice of mantle velocity model (orange residuals in Figs 6 and 7), with the use of SP12RTS
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Examples of comparisons of observed and predicted spectra showing the effect of density at stations ABPO (Madagascar), OBN and ARU (Russia), AAK (Kyrgyzstan), LVC (Chile), ADK (Alaska): (a) $|S_{11}|$, windowed between 20 and 60 hr and (b) $|S_{12}|$, windowed between 10 and 60 hr. The y-axes ($A$ and $\phi$) show normalized amplitude and phase. The normal mode frequencies (as predicted in PREM; Dziewonski & Anderson 1981) within each frequency window are also marked (dotted vertical lines). For each column of panels, horizontal frequency axes are shared.</figcaption></figure>

representing a median difference of ~22 per cent (Table 3). This reflects strong sensitivity to the $v_p$ structure, in agreement with early forward modelling of these modes by Koelemeijer et al. (2013, see their fig. 5). The $v_p$ structure is generally less well constrained at long wavelengths than $v_s$ as a result of uneven data sampling of the mantle, although recent models are showing more consistent results (Koelemeijer 2021). Meanwhile, differences due to the choice of 1-D attenuation model and the addition of CMB topography are very small, and negligible compared to other factors.

The uncertainty due to the source varies by several percent from station to station (not shown). Overall, as expected, the USGS source produces a significantly smaller difference with the reference CMT

source than the E35 solution of Abubakirov et al. (2015, Fig. 6). An accurate source depth is clearly important, although it should be said that here the use of splitting functions has the advantage of averaging over many sources.

The mean residual values between hypothetical density structures (the buoyant) reference, the ‘Dense Tall’ and ‘Dense Basal’ LLSVP models are consistent with the sensitivity kernels of $|S_{11}|$ and $|S_{12}|$. A relatively large mean residual between synthetics for the ‘Dense Tall’ LLSVP model and the reference is evident (Table 3 and Figs 6, 7). As expected, the mean residual between the reference and the ‘Dense Basal’ LLSVP model is extremely small, in line with the suggestion that a dense anomaly might be invisible to Stoneley
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Examples of comparisons of observed and predicted normalized spectra showing the effect of the background mantle model at stations ABPO (Madagascar), OBN and ARU (Russia), AAK (Kyrgyzstan), LVC (Chile), ADK (Alaska): (a) $|S_{11}|$, windowed between 20 and 60 hr and (b) $|S_{12}|$, windowed between 10 and 60 hr, similar to Fig. 8.</figcaption></figure>

mode data if confined to a sufficiently thin layer at the base of the mantle (Romanowicz 2017), which is the Stoneley mode blind spot for basal density variations.

## 3.2 Comparison with real data

Having considered various contributing factors to the synthetic residual discussed thus far, we now briefly compare our synthetic tests with real data.

Fig. 8 shows examples of individual spectra observed at several stations compared with synthetic spectra associated with each density model (calculated with full coupling of modes up to 3.1 mHz).

Different models fit the data better at different stations. For example, at stations ARU and ABPO (Fig. 8a) all models except the ‘Dense Tall’ model indistinguishably overestimate the spectrum amplitude for the mode pair $0S_{15}-1S_{11}$, while this density model underestimates the amplitude by a similar amount. Similar variability is seen for the pair $0S_{17}-1S_{12}$ (Fig. 8b). Finally, inspecting the other spectra in Fig. 8, the differences between all density models are rather small for mode pair $0S_{17}-1S_{12}$.

We now explore whether a better constrained $v_p$ model might improve the fit to these data. Fig. 9 is analogous to Fig. 8, where we compare our reference case for S20RTS with synthetics for the SP12RTS model and real data. Again, an inconsistent picture
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Box and whisker diagrams of the distribution of data residuals associated with each test case.</figcaption></figure>

emerges, where in some cases, for example at ADK, SP12RTS performs significantly better than our reference, whereas at OBN the opposite is true. We acknowledge that these only represent a small subset of spectra and so to explore this more holistically, we once more computed the residuals of the spectra, except that the reference ('ref' terms in eq. 2) is now associated with real data, and our reference model will be part of the collection of synthetic ('syn' term in eq. 2) tests under consideration. Table 3 and Fig. 10 display these results in an analogous fashion as before, comparing data against the four models plotted in Figs 8 and 9.

These misfits are substantially larger than most of the synthetic contributions explored in Section 3.1, as expected when working with real data. Based on these results, it seems impossible to distinguish between the two mantle models S20RTS and SP12RTS, given the similar distribution of data residuals (Figure 10). It should be kept in mind, however, that these tomographic models have been developed using a range of data with complementary sensitivity to seismic structures throughout the mantle (Ritsema et al. 1999; Koelemeijer et al. 2016).

## 3.3 Implications for constraining deep mantle density using Stoneley modes

The effect of varying lower mantle density structure is small, compared to the effect of mode-pair coupling, as illustrated for Stoneley modes $_1S_{11}$ and $_1S_{12}$ in Figs 6 and 7. This is consistent with recent work by Al-Attar et al. (2012) and Akbarashrafi et al. (2018) for non-Stoneley modes. Moreover, the distributions of individual data residuals are similar between the two density models and S20RTS and SP12RTS (Fig. 10). Other studies have stressed the importance of wide-band coupling of modes such as Deuss & Woodhouse (2001), Irving et al. (2008) and Yang & Tromp (2015). The results herein are in support of these finds, but since we perform an in-depth analysis only on two modes, coupling across ~0.1 mHz seems to be sufficient when considering the uncertainties in other factors. The wider implication is that considerable improvements must be made in our Earth models and theoretical considerations in order to improve our ability to discern density structure in the deep mantle.

It is important to note that the work presented here is based on two low-frequency Stoneley modes for one earthquake. Generalizing these results to other, similarly deep and large earthquakes is likely warranted, but extrapolating these results to higher frequency Stoneley modes such as $_2S_{16}$ and $_3S_{26}$ must be done with caution.

The sensitivity of these Stoneley modes is more focused at the CMB and their ability to resolve mantle density structure substantially better, as demonstrated by Koelemeijer et al. (2013, see their fig. 5). While these higher frequencies modes are also affected by P-wave velocity variations in the mantle, the effect on the inferred density structure is small (Koelemeijer et al. 2017). Nonetheless, the effect of theoretical considerations on higher frequency Stoneley modes requires further investigations. Particularly, it will be important to study whether narrow band coupling (coupling 10–20 modes) will be computationally feasible and sufficiently accurate, as it is shown here to be the case for low-frequency Stoneley modes.

## 4 CONCLUSIONS

In the context of trying to resolve the long wavelength density structure at the base of the mantle, we have explored two factors that affect the prediction of Stoneley mode spectra, focusing on low-frequency Stoneley modes $_1S_{11}$ and $_1S_{12}$. The first is associated with approximations regarding the extent of mode coupling, where we demonstrate that mode-pair coupling produces considerably inaccurate spectra. Narrow band coupling across ~0.1 mHz produces spectra with satisfactory accuracy, when considering other sources of uncertainty. The second factor is associated with effects of uncertainty on 3-D mantle structure. Here, we conclude that improved constraints on 3-D mantle $v_p$ structure are an important step for enabling the use of low-frequency Stoneley modes to determine other structures of interest, such as the density of the LLSVPs. In particular, we show that if excess density within the LLSVPs is confined to a thin basal layer, these low-frequency Stoneley mode spectra cannot detect the presence or absence of such a layer. Extending our analysis to higher frequency Stoneley modes is key to ascertain their use in studies of lower mantle density.

## ACKNOWLEDGEMENTS

This work was supported by NSF grant EAR-1923865. PK acknowledges support from a Royal Society University Research Fellowship (URF\R1\180377). Raw VHZ data were downloaded from the IRIS data centre (www.iris.edu). We thank network operators of the IRIS GSN, Geoscope, GEOFON for maintaining high quality very broad-band stations, from which data (doi.org/10.7914/SN/IC, doi.org/10.7914/SN/IU) were used in this study. We thank Dr Rudolf
Widmer-Schnidrig and an anonymous reviewer for their construc-
tive comments that improved this manuscript.

DATA AVAILABILITY

The software used for all calculations within in this paper will be
shared on reasonable request to the corresponding author.

REFERENCES

Abubakirov, I.R., Pavlov, V.M. & Titkov, N.N., 2015. The mechanism of the deep-focus, Sea of Okhotsk earthquake of May 24, 2013 as inferred from static displacements and broadband seismograms, *J. Volcanol. Seismol.*, **9**(4), 242–257.

Akbarashrafi, F., Al-Attar, D., Deuss, A., Trampert, J. & Valentine, A.P., 2018. Exact free oscillation spectra, splitting functions and the resolv-ability of Earth's density structure, *Geophys. J. Int.*, **213**(1), 58–76.

Al-Attar, D., Woodhouse, J.H. & Deuss, A., 2012. Calculation of normal mode spectra in laterally heterogeneous Earth models using an iterative direct solution method, *Geophys. J. Int.*, **189**(2), 1038–1046.

Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y. & Wassermann, J., 2010. ObsPy: a Python toolbox for seismology, *Seismol. Res. Lett.*, **81**(3), 530–533.

Cammarano, F. & Romanowicz, B., 2008. Radial profiles of seismic atten-
uation in the upper mantle based on physical models, *Geophys. J. Int.*
**175**(1), 116–134.

Dahlen, F.A. & Tromp, J., 1998. *Theoretical Global Seismology*, Princeton Univ. Press.

Davaille, A. & Romanowicz, B., 2020. Deflating the LLSVPs: bundles of mantle thermochemical plumes rather than thick stagnant “piles”, *Tectonics*, **39**(10), e06265, doi:10.1029/2020TC006265.

Davies, D.R., Goes, S., Davies, J.H., Schuberth, B.S.A., Bunge, H.P. & Ritsema, J., 2012. Reconciling dynamic and seismic models of Earth's lower mantle: the dominant role of thermal heterogeneity, *Earth planet. Sci. Lett.*, **353**, 253–269.

Deuss, A. & Woodhouse, J.H., 2001. Theoretical free-oscillation spectra: the importance of wide band coupling, *Geophys. J. Int.*, **146**, 833–842.

Durek, J.J. & Ekstrom, G., 1996. A radial model of anelasticity consistent with long-period surface-wave attenuation, *Bull. seism. Soc. Am.*, **86**(1A), 144–158.

Dziewonski, A.M., 1984. Mapping the lower mantle: determination of lateral heterogeneity in P velocity up to degree and order 6, *J. geophys. Res.*, **89**, 5929–5952.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**(4), 297–356.

Ekström, G., Nettles, M. & Dziewoński, A.M., 2012. The global CMT project 2004-2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200**, 1–9.

Forte, A.M. & Mitrovica, J.X., 2001. Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data, *Nature*, **410**(6832), 1049–1056.

Giardini, D., Li, X.-D. & Woodhouse, J.H., 1988. Splitting functions of long-
period normal modes of the Earth, *J. geophys. Res.*, **93**(B11), 13 716–13
742.

IC Network, n.d. doi.org/10.7914/SN/IC, GSN Dataset for IC network.

Irving, J. C.E., Deuss, A. & Andrews, J., 2008. Wide-band coupling of earth's normal modes due to anisotropic inner core structure, *Geophys. J. Int.*, **174**(3), 919–929.

Ishii, M. & Tromp, J., 1999. Normal-mode and free-air gravity constraints on lateral variations in velocity and density of earth's mantle, *Science*, **285**(5431), 1231–1236.

IU Network, n.d. doi.org/10.7914/SN/IU, GSN Dataset for IU network.

Jagt, L. & Deuss, A., 2021. Comparing one-step full-spectrum inversion with two-step splitting function inversion in normal mode tomography, *Geophys. J. Int.*, **227**(1), 559–575.

Jones, T.D., Maguire, R., van Keken, P.E., Ritsema, J. & Koelemeijer, P., 2020. Subducted oceanic crust as the origin of seismically slow lower-mantle structures, *Progress in Earth and Planetary Science*, **7**(17), 255–267.

Karato, S.-I., 1993. Importance of anelasticity in the interpretation of seismic tomography, *Geophys. Res. Lett.*, **20**, 1623-1626.

Koelemeijer, P., 2021. Towards Consistent Seismological Models of the Core-Mantle Boundary Landscape, Marquardt, Ballmer, Cottaar & Konter *Mantle Convection and Surface Expressions*, Wiley-Blackwell. 229–255 doi:10.1002/9781119528609.ch9

Koelemeijer, P., Deuss, A. & Ritsema, J., 2013. Observations of core-mantle boundary Stoneley modes, *Geophys. Res. Lett.*, **40**(11), 2557–2561.

Koelemeijer, P., Ritsema, J., Deuss, A. & van Heijst, H., 2016. SP12RTS: a degree-12 model of shear- and compressional-wave velocity for Earth's mantle, *Geophys. J. Int.*, **204**(2), 1024–1039.

Koelemeijer, P., Deuss, A. & Ritsema, J., 2017. Density structure of Earth's lowermost mantle from Stoneley mode splitting observations, *Nat. Commun.*, **8**, 15241.

Kuo, C. & Romanowicz, B., 2002. On the resolution of density anomalies in the Earth's mantle using spectral fitting of normal-mode data, *Geophys. J. Int.*, **150**(1), 162–179.

Lau, H.C.P., Mitrovica, J.X., Davis, J.L., Yang, H.-Y. & Al-Attar, D., 2017. Tidal tomography constrains Earth's deep-mantle buoyancy, *Nature*, **551**, 321–326, doi.org/10.1038/nature24452.

Lekic, V., Cottaar, S., Dziewonski, A. & Romanowicz, B., 2012. Cluster analysis of global lower mantle tomography: a new class of structure and implications for chemical heterogeneity, *Earth planet. Sci. Lett.*, **357**, 68–77.

Masters, G., Laske, G. & Gilbert, F., 2000. Matrix autoregressive analysis of free-oscillation coupling and splitting, *Geophys. J. Int.*, **143**(2), 478–489.

McNamara, A.K., 2019. A review of large low shear velocity provinces and ultra low velocity zones, *Tectonophysics*, **760**, 199–220.

Mosca, I., Cobden, L., Deuss, A., Ritsema, J. & Trampert, J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, *J. geophys. Res.*, **117**(B6), B06304, doi:10.1029/2011JB008851.

Moulik, P. & Ekström, G., 2016. The relationships between large-scale variations in shear velocity, density, and compressional velocity in the Earth's mantle, *J. geophys. Res.*, **121**, 2737–2771.

Mulyukova, E., Steinberger, B., Dabrowski, M. & Sobolev, S.V., 2015. Survival of LLSVPs for billions of years in a vigorously convecting mantle: Replenishment and destruction of chemical anomaly, *J. geophys. Res.*, **120**(5), 3824–3847.

Resovsky, J. & Trampert, J., 2003. Using probabilistic seismic tomography to test mantle velocity-density relationships, *Earth planet. Sci. Lett.*, **215**(1-2), 121-134.

Ritsema, J., van Heijst, H.J. & Woodhouse, J.H., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland, *Science*, **286**(5446), 1925–1928.

Romanowicz, B., 2001. Can we resolve 3D density heterogeneity in the lower mantle? *Geophys. Res. Lett.*, **28**(6), 1107–1110.

Romanowicz, B., 2017. Geophysics: the buoyancy of Earth's deep mantle, *Nature*, **551**(7680), 308–309.

Romanowicz, B. & Mitchell, B., 2015. 1.25 - Deep earth structure: Q of the Earth from crust to core, in *Treatise on Geophysics*, 2nd edn, pp. 789–827, ed Schubert, G., Elsevier.

Simmons, N.A., Forte, A.M., Boschi, L. & Grand, S.P., 2010. GyPSuM: a joint tomographic model of mantle density and seismic wave speeds, *J. geophys. Res.*, **115**(B12),

Torsvik, T.H., Smethurst, M.A., Burke, K. & Steinberger, B., 2006. Large igneous provinces generated from the margins of the large low-velocity provinces in the deep mantle, *Geophys. J. Int.*, **167**(3), 1447–1460.

Trampert, J., Deschamps, F., Resovsky, J. & Yuen, D., 2004. Probabilistic tomography maps chemical heterogeneities throughout the lower mantle, *Science*, **306**(5697), 853–856.
Widmer, R., Masters, G. & Gilbert, F., 1991. Spherically symmetric attenuation within the Earth from normal mode data., *Geophys. J. Int.*, **104**(3), 541–553.

Yang, H.-Y. & Tromp, J., 2015. Synthetic free-oscillation spectra: an appraisal of various mode-coupling methods, *Geophys. J. Int.*, **203**(2), 1179–1192.

Zhong, S. & Rudolph, M.L., 2015. On the temporal evolution of long-wavelength mantle structure of the Earth since the early Paleozoic, *Geochem. Geophys. Geosyst.*, **16**(5), 1599–1615.

## SUPPORTING INFORMATION

Supplementary data are available at *GJI* online.

S1.pdf

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.