# Constraining Jumps in Density and Elastic Properties at the 660 km Discontinuity Using Normal Mode Data via the Backus Gilbert Method

Harriet C. P. Lau, Barbara Romanowicz

► To cite this version:

Harriet C. P. Lau, Barbara Romanowicz. Constraining Jumps in Density and Elastic Properties at the 660 km Discontinuity Using Normal Mode Data via the Backus Gilbert Method. Geophysical Research Letters, 2021, 48, pp.97-110. 10.1029/2020GL092217. insu-03590053

HAL Id: insu-03590053

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

Submitted on 26 Feb 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
# Geophysical Research Letters

## RESEARCH LETTER

10.1029/2020GL092217

**Key Points:**

* We use recent normal mode center frequency data to constrain the elastic/density properties of the mantle 660 km discontinuity

* We find that acceptable range of jumps in P-wave speed and density fall outside that of standard seismic reference models

* Our data preclude the global discontinuity being as shallow as 650 km depth

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

H. C. P. Lau,
hcplau@berkeley.edu

**Citation:**

Lau, H. C. P., & Romanowicz, B. (2021). Constraining jumps in density and elastic properties at the 660 km discontinuity using normal mode data via the Backus-Gilbert method. *Geophysical Research Letters*, *48*, e2020GL092217. https://doi.org/10.1029/2020GL092217

Received 21 DEC 2020
Accepted 18 APR 2021

© 2021. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution-NonCommercial-NoDerivs
License, which permits use and
distribution in any medium, provided
the original work is properly cited,
the use is non-commercial and no
modifications or adaptations are made.

# Constraining Jumps in Density and Elastic Properties at the 660 km Discontinuity Using Normal Mode Data via the Backus-Gilbert Method

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

Harriet C. P. Lau¹ ID and Barbara Romanowicz¹,²,³ ID

¹Department of Earth and Planetary Science, University of California, Berkeley, CA, USA, ²College de France, Paris, France, ³Institut de Physique du Globe de Paris, Paris, France

**Abstract** We apply the Backus-Gilbert approach to normal mode center frequency data, to constrain jumps in P, S, bulk-sound speed and density at the “660” discontinuity in the earth’s mantle (~650–670 km depth). Different 1-D models are considered to compute sensitivity kernels. When using model PREM (Dziewonski & Anderson, 1981, Physics of the Earth and Planetary Interiors, 25, 297–356. doi:10.1016/0031-9201(81)90046-7) as reference, with a “660” at 670 km depth, the best-fitting jumps in density, P- and S-wave speeds range from (5.1–8.2)%, (5.3–8.0)%, (5.0–7.0)%, respectively, so the PREM values lie outside the ranges of acceptable density and P wave speed jumps. When shifting the depth of “660” to 660 km, the density and S wave speed jumps increase, while the P-wave speed jump decreases. Normal mode data do not support a global transition at 650 km depth. The density jumps are closer to those of pyrolite than PREM, while our bulk-sound wave speed jumps suggest a larger garnet proportion at “660.”

**Plain Language Summary** Within the mantle, at around 650–670 km depth, rocks undergo a phase change due to an increase in temperature and pressure, from a relatively less dense mineral phase of ringwoodite to the denser phases of bridgmanite and periclase. This transition often marks the boundary between the upper and lower mantle. Switching between these two phases holds many dynamical consequences, when it comes to the convection of the mantle. We use low frequency seismic data known as normal modes, vibrations of the entire planet, to revisit the elastic and density properties of the mantle at this transition. We find that changes in the density and P-wave speed across the boundary are different to what most standard seismic models have. We also find that the transition cannot be as shallow as 650 km. These conclusions bear consequences for those who research the convection and composition of the mantle and global seismologists.

## 1. Introduction

Phase transitions that occur throughout the mantle greatly affect mantle dynamics and their precise location can provide information about the thermal and compositional variations within the earth. In this study, we focus on the so-called “660” discontinuity (hereafter, 660) which has been observed to occur at depths between ~650 and 670 km, attributed to the transition between the mineral phases ringwoodite/ spinel at lower pressures to bridgmanite and oxides at higher pressures (e.g., Birch, 1952; Frost, 2008; Ring- wood, 1991; Shearer, 2000; Shim et al., 2001; Stixrude & Lithgow-Bertelloni, 2005). Both seismology and mineral physics have provided important insights into the nature of 660.

Such a phase transition will produce sharp jumps in seismic wave-speed, manifested by seismological observations of reflected phases such as precursors to short period P’P’ phases (e.g., Xu et al., 2003), precursors to SS (e.g., Shearer, 2000) and PP phases (e.g., Deuss et al., 2006), and converted phases as detected in receiver function studies (e.g., Andrews & Deuss, 2008). The depth of the sharp jumps in wave speed listed in seismic reference, spherically symmetric (1-D) models vary from 670 km (for PREM, Dziewonski & Anderson, 1981) to 650 km (for STW105, Kustowski et al., 2008), with a currently preferred value of 660 km. At the global scale, the topography of this discontinuity reaches upto ±30 km (Andrews & Deuss, 2008), with somewhat larger excursions locally, for example, in subduction zones (e.g., Niu & Kawakatsu, 1995). Observed jumps in S-wave speed, $\Delta v_s$, P-wave speed, $\Delta v_p$, and in density, $\Delta \rho$, range from 4.5%–10.1%, 2.5%–5.6%, and 4.2%–10.2%, respectively (Estabrook & Kind, 1996; Castle & Creager, 2000; Kennett & Engdahl, 1991; Montagner
& Anderson, 1989; Morelli & Dziewonski, 1993; Shearer & Flanagan, 1999). Along with inherent trade-offs between the different physical parameters, the complicated nature of seismic signals observed across the boundary itself must contribute to the wide range of seismically observed jumps (Andrews & Deuss, 2008).

Efforts to combine the mineral physics and seismological approaches aim to tie physical causes to observed seismic properties. For example, by applying equations of state derived from mineral physics, assuming a mantle of adiabatic pyrolite composition, Cammarano et al. (2005) showed that wave-speed jumps that satisfy seismic reference models lay toward the higher end of permissible values from mineral physics constraints. The depth of the discontinuity also provides insights into the nonpyrolitic components (e.g., ilmenite, garnet) within the transition zone (e.g., Ishii et al., 2018; Vacher et al., 1998; Wang et al., 2006).

In this study, we revisit the estimation of globally averaged $\Delta v_s$, $\Delta v_p$, $\Delta v_b$, and $\Delta \rho$ across 660 by applying Backus-Gilbert based methods (Backus & Gilbert, 1970; Masters & Gubbins, 2003; Pijpers & Thompson, 1992) to an extensive recent normal mode catalog (Deuss et al., 2013; Roult et al., 2010).

## 2. Data and Methodology

### 2.1. Data

We use the normal mode center frequencies (and uncertainties) compiled by Robson and Romanowicz (2019) which is based on a combination of the Reference Earth Model catalog (Laske, Dziewonski, & Masters, 2001), observations from Deuss et al. (2013) and radial modes from Roult et al. (2010). The data are provided in Supplementary Table 1.

### 2.2. Methodology

In a spherical elastic and isotropic earth model, the eigenfrequency $\omega_k$ of any isolated normal mode multiplet (denoted by the index $k$) has distinct sensitivity kernels to $v_s$, $v_p$, and $\rho$ structure across the mantle and to topography of any discontinuity, $d_i$. In the framework of first order perturbation theory, the fractional change in eigenfrequency may be expressed as (e.g., Woodhouse & Dahlen, 1978):

$$ \frac{\delta \omega_k}{\omega_k} = \int_0^a \left[ M_{v_p}^k(r) \frac{\delta v_p}{v_p}(r) + M_{v_s}^k(r) \frac{\delta v_s}{v_s}(r) + M_{\rho}^k(r) \frac{\delta v_{\rho}}{v_{\rho}}(r) \right] dr + \sum_i M_{d_i}^k \delta d_i, \quad (1) $$

where $r$ is the radius and $r = a$ is the surface, $M_{v_p}^k$, $M_{v_s}^k$, $M_{\rho}^k$ are the sensitivity kernels of mode $k$ to perturbations in $v_p$, $v_s$, and $\rho$, respectively. $M_{d_i}^k$ is the sensitivity to topography, $d$, on the $i$th discontinuity.

By considering a linear combination of Equation 1 over a set of modes $k$, we obtain:

$$ \sum_k c_k \frac{\delta \omega_k}{\omega_k} = \int_0^a \left[ K_{v_p}(r) \frac{\delta v_p}{v_p}(r) + K_{v_s}(r) \frac{\delta v_s}{v_s}(r) + K_{\rho}(r) \frac{\delta v_{\rho}}{v_{\rho}}(r) \right] dr + \sum_i K_{d_i} \delta d_i \quad (2) $$

where $K_X = \sum_k c_k M_X^k$ for parameter $X$ which, in this study, $X = v_p, v_s$ or $\rho$. The coefficients $c_k$ may be determined such that $K_X$ is designed to enhance the sensitivity of the weighted observations (left-hand side of Equation 2) to a specific region within the mantle and a specific parameter $X$, while simultaneously reducing the sensitivity to other parameters, $Y$, $Z$, and $d_i$. If the weights $c$ are successfully determined, in the most ideal case $K_X$ will be only nonzero across the region of interest, and $K_Y$, $K_Z$ and $K_{d_i}$ will be zero everywhere. We will refer to both the weighted kernels and data as composite kernels and data. An additional condition required of the composite kernel is that it should be unimodular:

$$ \int_0^a K_X(r) dr = 1. \quad (3) $$

This is the essence of the Backus-Gilbert methodology. Finding the best combination of data, that is, finding $c$, requires solving an inverse problem and thus carries with it the same regularization issues as in typical geophysical inverse problems.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Composite kernels and inversion results. (a) Target kernels for study (blue, *T*), a half-Gaussian beneath 660; orange, *T*<sub>s</sub>, a half-Gaussian above, and black, *T*<sub>full</sub>, a full-Gaussian straddling both beneath and above 660). (b–d) The resulting kernels to enhance sensitivity to ρ, v<sub>p</sub>, and v<sub>s</sub>, respectively. (In the inset panels, the resulting kernels for the other parameters chosen to be muted). The top panels show relative sensitivity to topography on each discontinuity. (e–g) Resulting perturbations in parameters *X*, where *X* is ρ, v<sub>p</sub>, and v<sub>s</sub> in panels (a, b, c), respectively. The circles denote how each kernel samples the background PREM model, while the crosses are the resulting perturbations when applying the composite data. The color intensity of the bar represents the kernel density (see panels b–d). The width of each bar is the uncertainty in the result (Equations 14 and 15). The gray lines display the PREM profile. (h) A schematic depiction of the test models we consider. The perturbations above and below 660 are guided by solutions shown in panels (e–g), and we consider these to increase or decrease linearly from the background PREM model in two cases: 100 km above and below the discontinuity (dashed lines) and 200 km above and below the discontinuity (not depicted). The thick black line is the scaled kernel that is associated with the black bars in panels (e–g). This is used as a criterion to exclude test models that do not fit the constraint from the black crosses in (e–g). (i) Schematic diagram of how we adjust PREM for a shallower discontinuity at both 660 and 650 km depth, preserving the percent jump value to that of the standard PREM model. Panels (e–h) share the same vertical axes.</figcaption></figure>

To expand upon this, we introduce the concept of a target kernel, *T*, as introduced by Pijpers and Thompson (1992), whose methodology we closely follow (though they considered only one free parameter). *T* will be designed such that it follows the shape of the desired sensitivity. Here, we will explore three kernels: (1) a narrow Gaussian centered at 660 (solid black line, Figure 1a) which is defined as:
$$ T_{\text{full}} = \frac{1}{\Lambda} \exp\left(-\left(\frac{r-r_0}{\Delta}\right)^2\right) \qquad (4) $$

where $\Lambda$ is chosen so that the area under $T_{\text{full}}$ is 1, $\Delta$ is the characteristic width of the Gaussian centered at $r = r_0$. The remaining two target kernels are truncated versions of this Gaussian, one where the kernel is identical to the full Gaussian above 660, but is zero below 660, $T_+$ (orange kernel, Figure 1a) and the other where the kernel is identical to the full Gaussian below 660, but zero above, $T_-$ (blue kernel, Figure 1a). The truncated kernels, $T_-$ and $T_+$, will provide estimates on either side of the 660, which we will use to determine new jump constraints for $\rho$, $v_p$, and $v_s$. The full Gaussian, $T_{\text{full}}$ will provide an overall constraint across the 660 boundary when testing synthetic models in Section 3.2.

In order to determine $c$ such that the resulting $K_X$ is as similar to $T$ as possible, we minimize the following expression:

$$ \Phi = \int_{0}^{a} \left[ (K_X - T)^2 + K_Y^2 + K_Z^2 \right] dr + \sum_{i} K_{d,i}^2 + \mu \sum_{ij} E_{ij} c_i c_j, \qquad (5) $$

where $E$ is the covariance matrix of data errors and $\mu$ is its corresponding trade-off parameter. In Pijpers and Thompson (1992), they used an additional constraint to ensure that the area under $K$ is unity. Here, since $T$ is designed to be this way, we do not include this additional constraint. Minimizing $\Phi$ with respect to $N$ unknowns (where $N$ is the number of normal mode center frequencies considered) yields $N$ linear equations which have the form:

$$ \sum_j \left[ \int_0^a M_X^i M_X^j + M_Y^i M_Y^j + M_Z^i M_Z^j dr + \sum_d M_{d,i} M_{d,j} + \mu E_{ij} \right] c_j - \int_0^a M_X^i T dr = 0. \qquad (6) $$

and may be written in matrix form as:

$$ A_c = v \qquad (7) $$

where $c = [c_1, c_2, c_3, ..., c_N]$ and vector $v$ is

$$ v_i = \int_0^a M_X^i T dr \qquad (8) $$

for $i = 1, 2, ..., N$. The elements of the $N \times N$ symmetric matrix $A$ are:

$$ A_{ij} = \int_{0}^{a} M_{X}^{i} M_{X}^{j} + M_{Y}^{i} M_{Y}^{j} + M_{Z}^{i} M_{Z}^{j} dr + \sum_{d} M_{d,i} M_{d,j} + \mu E_{ij}, \qquad (9) $$

where $i, j = 1, 2, ..., N$. Since $\int_0^a T dr = 1$, an estimate of the quantity of interest, $\tilde{X}$, may be obtained as follows. (Note that in the following expressions we make explicit any dependence on $r$ and that the “~” symbol denotes any value integrated over $r$.)

$$ \sum_k c_k \frac{\delta \omega_k}{\omega_k} = \int_0^a \frac{\delta X}{X_0}(r) K_X(r) dr = \left( \frac{\delta X}{X_0} \right), \qquad (10) $$

where $X_0$ is the unperturbed depth profile of parameter $X$ from which the kernels $K_X$ were determined. Taking $\delta X = X - X_0$, we may write
$$1 + \sum_k c_k \frac{\delta\omega_k}{\omega_k} = \int_0^a \frac{X}{X_0}(r) K_X(r) dr \qquad (11)$$

(the value of 1 arising from Equation 3). The radially dependent solution is thus:

$$X(r) = X_0(r) \left[ 1 + \frac{\delta X}{X_0}(r) \right]. \qquad (12)$$

However, as Equation 10 does not yield radially dependent $\delta X/X_0$, our solution, $\tilde{X}$, is an approximation of this (i.e., $\tilde{X} \approx X$) over the radial range for which the kernel is nonzero (or nonnegligible):

$$\tilde{X}(r) = X_0(r) \left[ 1 + \left( \frac{\delta X}{X_0} \right) \right], \qquad (13)$$

where $r$ is the radius of interest.

The uncertainty in the estimate of $\tilde{X}$, $\epsilon$, is due to two sources: data error, $\epsilon_{obs}$, and contamination from imperfections in the composite kernels, $\epsilon_{con}$, since in practice, they will not be fully zero where desired. The contribution from errors in observation is given by

$$\epsilon_{obs}^2 = c \cdot E \cdot c, \qquad (14)$$

and $\epsilon_{con}$ is due to nonzero contributions from $K_Y$, $K_Z$ and $K_{d,i}$ (Masters & Gubbins, 2003) but also differences between $T$ and $K_X$. This may be estimated by the following expression

$$\epsilon_{con}^2 = \int_0^a |K_X - T| \parallel \epsilon_X | + |K_Y - T| \parallel \epsilon_Y | + |K_Z - T| \parallel \epsilon_Z | dr + \sum_i |K_{d,i} - T| \parallel \epsilon_d |, \qquad (15)$$

where $\epsilon_X$, $\epsilon_Y$, $\epsilon_Z$, and $\epsilon_d$ are uncertainties in the parameters $X$, $Y$, $Z$, and $d$, respectively, and $|\cdot|$ denotes taking the absolute value.

# 3. Results

## 3.1. Composite Kernels

We calculated the sensitivity kernels for each mode to each parameter $M$ according to Dahlen and Tromp (1998), using the software package MINEOS (https://geodynamics.org/cig/software/mineos/), adopting the widely used seismic reference model PREM (Dziewonski & Anderson, 1981) as reference. In Section 4.1, we explore the effect of the choice of 1-D reference model.

We solve for three separate sets of three kernels $T_{full}$, $T_-$, $T_+$ (Figure 1a). In the first set, we enhance sensitivity to $\rho$ and suppress sensitivities to $v_p$, $v_s$, and $d_i$ (Figure 1b). In Figures 1c and 1d we show the analogous resulting composite kernels for enhancing sensitivity to $v_p$ and $v_s$, respectively.

While not perfect, the overall shape of the composite kernels $K_X$ capture $T$ very well, with small amounts of noise in all composite kernels (i.e., $K_Y$ and $K_Z$), which will add to the uncertainty in the estimate.

## 3.2. New Estimates on Jumps

In Figures 1e–1g, we show the results obtained by solving Equation 7 successively for each set of composite kernels. In each figure, it is important to focus on the shift between the circle and cross, rather than the shift from the PREM profile. Focusing on the truncated-Gaussian solutions (orange and blue) and the width of the kernel density bar, it can be seen that for all parameters, the available normal mode data better constrain the parameter beneath 660 (which is at 670 km in PREM) than above. For $v_s$, the PREM values satisfy those of the composite data well (crosses overlap the circles in Figure 1g). However, the composite data call for
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Resulting jumps across 660. (a) Black bars (Backus-Gilbert solutions) depict the jumps based on the results Figures 1e–1g. The length of the vertical gray lines here corresponds to the width of the colored bars in the latter figure. The blue boxes denote synthetically produced models tested against the data, perturbing PREM above and below the discontinuity over a length scale of 100 km (Figure 1h). These ranges are less than the gray bars as models that did not satisfy the constraint imposed by $\mathcal{T}_{\text{full}}$ (condition 17) were culled. Symbols represent values from different seismic reference models and from molecular dynamics (MD) calculations. (b) Table of jumps across 660 all listed as percentages. The references for each are as follows: (1): Dziewonski and Anderson (1981); (2): Kennett et al. (1995); (3): Kustowski et al. (2008); (4): Matsui (2000).</figcaption></figure>

a smaller density jump (Figure 1e), and for a larger jump in $v_p$ (Figure 1f) than in PREM, with shifts both above and below 660.

When one considers averaged values of PREM both above and below 660, PREM satisfies the composite data for all parameters. This is indicated by the black circles falling on the black crosses in all panels (Figures 1e–1g).

To explore this further, we tested how these new estimates perform in reproducing the composite data. We perturbed PREM both above and below its 670 km discontinuity with our new estimates. The values of $v_s$, $v_p$ and $\rho$ were linearly interpolated to the PREM background model across two length scales: 100 and 200 km. A schematic depiction of this perturbation is shown in Figure 1h where we show test models for the 100 km length scale. These perturbations to the background model may change the nature of the eigenfrequencies (Dahlen & Tromp, 1998) but recalculating the eigenfrequencies with these updated models and testing them against the data is one way to ensure that these models have not changed too much.

We produced many perturbed PREM models ($X + \delta X$), choosing values of parameters above and below the discontinuity within the uncertainty shown in Figures 1e–1g, applying all possible combinations of $\rho$, $v_p$, and $v_s$. This resulted in 46,656 models for each length scale tested. To further scrutinize these models, and before confronting them with the data, we tested whether these models met the constraint provided by $\mathcal{T}_{\text{full}}$ (shown in a scaled version in Figure 1h). This condition simply requires that the value $\bar{X}$ must lie within the horizontal span of the associated black colored bar where

$$ \bar{X} = \int_{0}^{a} (X + \delta X) \mathcal{T}_{\text{full}} dr \qquad (16) $$

and $(X + \delta X)$ is the perturbed model tested. That is,

$$ (\tilde{X}_{\text{full}} - \varepsilon) \le \bar{X} \le (\tilde{X}_{\text{full}} + \varepsilon). \qquad (17) $$

The result of this additional condition is that none of the models perturbed across a length scale of 200 km were able to meet the constraint, whereas for the 100 km length scale, 10,000 models satisfied the constraint. Suggesting that any perturbation from PREM cannot be too wide.

For all these results, listed in Figure 2b, the associated jumps are displayed as percentages in Figure 2a, where the final models tested are shown by the blue boxes. We note that after the initial culling of models,
$\Delta\rho$ and $\Delta v_p$ are significantly different from their respective PREM values. We subject this culled subset of models to an additional test, as follows. For all models, we predict the set of composite data and define the chi-squared misfit, $\chi^2$, as

$$ \chi^2 = \sum_i \frac{(\Omega_i^{\text{mod}} - \Omega_i^{\text{obs}})^2}{\sigma_i^2}. \qquad (18) $$

We note that $\Omega_i$ is the composite datum from each inversion performed where $i = [\rho^{-}, \rho^{+}, v_{p}^{-}, v_{p}^{+}, v_{s}^{-}, v_{s}^{+}]$ and $\Omega_i = \sum_{k=1}^{N} c_k^i \omega_k$. Indeed, for each parameter $i$ enhanced, a different set of coefficients $c_k^i$ is obtained, where $k$ is the index of a mode. Each inversion is accompanied by a composite uncertainty, $\sigma_i$, weighted in the same manner.

We present the misfit reduction, $\gamma$, in Figures 3a–3d, where $\gamma$ is the ratio of $\chi^2$ calculated from the test models to $\chi^2$ calculated by PREM. These models correspond only to the perturbed models over a length scale of 100 km, since the 200 km models did not satisfy the $T_{\text{full}}$ constraint (Equation 17). They span the blue boxes shown in Figure 2a. Trade-offs between wave-speeds $\Delta v_s$ and $\Delta v_p$ are shown in Figures 3a and 3c, while trade-offs between wave-speeds $\Delta v_s$ and $\Delta v_p$ are shown in Figures 3b and 3d. Each row corresponds to models at two fixed density jumps, where $\Delta\rho$ is 5.1% and 8.2% for panels (a, b) and (c, d), respectively. These values correspond to the minimum and maximum $\Delta\rho$ values in the blue boxes of Figure 3.

In the models tested, the misfit was reduced (i.e., $\gamma < 1$) for a significant portion of the models. We see a general preference for lower $\Delta v_s$ across the range tested and preference for higher $\Delta v_p$ across this range, though these two parameters display some covariance. A trade-off also exists between $v_s$ and $v_p$, whereby larger values of $\Delta v_s$ are paired with smaller values of $\Delta v_p$.

We performed F-tests for all cases at the 99% level of significance (solid red line). For $\Delta v_s$, we see that we reach levels of 99% significance in the region of ~ (4.5–7.5)% for a $\Delta\rho$ of 5.1% and ~ (5.0–7.5)% for $\Delta\rho$ of 8.2%. As such, we report that these composite data provide revised estimates of these jumps at a 99% significance value of: $\Delta\rho = (5.1–8.2)\%$, $\Delta v_s = (5.0–7.0)\%$, $\Delta v_p = (5.3–8.0)\%$ and $\Delta v_p = (4.0–9.5)\%$, though these ranges are all correlated and should not be taken at face value. The original PREM values are 9.7%, 6.7%, 4.7%, and 3.4%, respectively (see also Figure 2b).

As a further point of comparison, impedance contrasts across 660 are well constrained by body wave studies (e.g., Shearer & Flanagan, 1999). When comparing $v_s$ and $v_p$ impedance contrasts with SS and PP precursors (as summarized in Deuss, 2009), PREM, pyrolite, and piclogite show significantly higher contrasts in both $v_s$ and $v_p$ (~0.13–0.16 for both) than body wave-derived estimates (~0.08–0.11 in $v_s$ and ~0.05–0.08 in $v_p$). In Figure 11 of Deuss (2009), the trade-off between these two quantities is highlighted clearly. Our $v_s$ and $v_p$ contrasts are 0.08–0.14 and 0.09–0.16, respectively, aligning much closer to PREM and mineral physics values, though do span values close to the upper ends of body wave inferences. It is not clear why the lower frequency normal modes might see a stronger impedance. Perhaps the lower resolution of the normal mode data may be sensitive to velocity gradients either side of the discontinuity which may result in the impression of a larger jump. However, further investigation is required.

# 4. Discussion

## 4.1. Effect of Background Model

In both Figures 2a and 3, we have overlain our results with the corresponding jumps in PREM (gray stars) and the seismological reference models “ak135” (orange circles, Kennett et al., 1995) and “STW105” (yellow circles, Kustowski et al., 2008). While for PREM, the largest differences are seen in $\Delta\rho$ and $\Delta v_p$, our results are in much closer agreement with the other two seismic models in $\Delta\rho$ and with ak135 for $\Delta v_p$. However, our estimate for $\Delta v_s$, while in agreement with PREM and STW105, is substantially different from that of ak135 (which is more than ~2% different). These differences may reflect the fact that that STW105 and PREM were derived from much more similar datasets than ak135, which consisted mainly of short period travel time data sampling the earth’s upper mantle beneath continental areas. We note that no single model shows an obvious consistency with our new estimates for all three parameters.
<figure><img src="image_5.png" /></figure>

**Figure 3.** Misfits from synthetic tests and their trade-offs. Panels (a–d) display contour plots of the misfit reduction, $\gamma$, for models that span the blue boxes in Figure 2a (i.e., perturbed models over a 100 km length scale). Each row displays models at fixed $\Delta\rho$ values of 5.1% (panels a, b) and 8.2% (panels c, d). These $\Delta\rho$ values span the full range of models that satisfied condition (17). Left panels display trade-offs between $\Delta v_s$ and $\Delta v_p$ and right panels display trade-offs between $\Delta v_s$ and $\Delta v_b$. Panels (e–h) show $\gamma$ for test models that meet the criterion described by condition (17), assuming different background models: the standard version of PREM (blue circles, also shown in panels a–d, where dark blue circles distinguish a subset of these models for which $\Delta\rho$ is the minimum possible density jump of 5.1% highlighting the directionality of trade-offs), PREM in which the discontinuity has been artificially shifted upward to 660 km (orange circles, see Figure 1e), and assuming the background model ak135 (orange shaded range).

PREM and ak135 provide a good representation of the range of background models available given the differing nature of the datasets used to produce them. However, another notable difference between these reference models is that PREM places the 660 at a depth of 670 km, whereas ak135 places this discontinuity at 660 km. Here we test how robust our results are to the choice of background model, and specifically, the depth of the 660 in the reference model.
Starting from PREM, we artificially adjusted the discontinuity to two depths: 650 and 660 km. On the upper side of the discontinuity, we kept the PREM parameter values down to the new discontinuity depth. On the underside, we imposed the density and wave speeds that conserve the original percent jump of PREM for each quantity, in order to isolate the effect of the discontinuity depth. We then linearly interpolated to the original PREM values at 670 km depth (Figure 1e). For each case, we produced new composite kernels, $K_+$, $K_-$, and $K_{\text{full}}$ and repeated the analysis of Sections 2 and 3.

No model with the discontinuity fixed at 650 km satisfied the condition (17), whereas where the depth of discontinuity was 660 km, 9,460 models satisfied this constraint. We show the misfit reduction for the latter distribution of models in Figures 3e–3g (dark orange circles). For comparison, we display the results from the standard analysis (blue circles). By elevating the discontinuity to a shallower depth of 660 km, the trade-offs between $\Delta v_p$ and $\Delta v_s$ switch direction, requiring an increase of $\Delta\rho$ and $\Delta v_s$ and a reduction of $\Delta v_p$. Furthermore, it seems that a global average depth of 650 km for the 660 is too shallow to satisfy normal mode data. In line with this, mineral physics experiments of, for example, Shim et al. (2001) and Ishii et al. (2018), also indicate a deeper depth is favored.

These trade-offs between a shallower discontinuity and increase in both $\Delta\rho$ and $\Delta v_s$ and a reduction of $\Delta v_p$ seem to be consistent if we repeat the entire exercise with ak135 (Kennett et al., 1995). We find values corresponding to those of the gray bars in Figure 2a of (5.1–12.8)%, (1.4–4.3)%, and (6.5–9.5)% for $\Delta\rho$, $\Delta v_p$, and $\Delta v_s$, respectively. The same culling exercise that reduced the gray bars to the dark blue boxes in Figure 2a (for a perturbation length scale of 100 km) did not result in significant changes to these ranges.

For $\Delta\rho$, the ak135 estimate is similarly poorly constrained relative to the gray bars, but shifted to higher density jumps. For $\Delta v_p$ and $\Delta v_s$, the span of the ak135 estimates are roughly two-thirds of the gray bars (compare Figures 3e–3h with the black bars in Figure 2a). In Figures 3e–3h (orange circles), the misfit reductions of the modified PREM model (with a depth of discontinuity of 660 km, as in ak135) reduce toward the ranges spanned by the ak135 result for all parameters except for $\Delta v_b$.

These differences between the Backus-Gilbert solutions based on PREM and ak135 are not trivial, and illustrate the strong nonlinearity of the problem, combining the effects of the depth of the discontinuity and of the jumps in the three parameters considered. Since ak135 was constrained by a very different type of data, these results likely represent an unrealistic “worst case scenario” when applied to normal mode center frequency data.

## 4.2. Physical Implications

The characteristics of the 660 phase boundary have had much influence on the conceptual picture of mantle convection. Given that the negative Clapeyron slope implies that the transition shifts to higher pressures at colder temperatures, the idea that this transition is a barrier to general mantle circulation has been proposed extensively to satisfy geochemical constraints on mantle heterogeneity (e.g., Allgre, 1997) and explored dynamically with consideration to how parameters such as mantle viscosity might be affected (e.g., van Keken & Ballentine, 1998). While the picture of mantle convection continues to evolve with better seismic imaging techniques and more sophisticated modeling capabilities, many of these approaches require a background seismic model, and, in many cases, PREM is the model of choice.

As such, revising the globally averaged characteristics above and beneath the 660 in the light of recent data is important for both the geodynamical and seismological communities. Paired with increasingly accurate measurements from the mineral physics literature, a more complete picture of the physical characteristics of this region will be reached.

The mineral physics estimates of Matsui (2000) for $\Delta\rho$ and $\Delta v_b$ of the model mantle compositions of pyrolite and piclogite (Bass & Anderson, 1984; Ringwood, 1962) are included in Figures 2a and 3a–3d (green crosses and lines, respectively). Our estimates of $\Delta\rho$ are more in line with these mineral physics estimates relative to PREM, being closer to a pyrolitic composition. This is further visualized in Figure 3e, where misfit reductions point toward the estimate for pyrolite. However, we do not improve the $\Delta v_b$ fit to either of these possible compositions (Figure 3h). It seems from our results in Section 4.1 that these same conclusions stand whether we consider a discontinuity at 670 km as in PREM, or 660 km as in ak135.
Nonolivine components of the upper mantle, and in particular, the presence of ilmenite, may affect such jumps across 660 (e.g., Vacher et al., 1998). More recently, Ishii et al. (2018) explored the transition of ring-woodite to garnet and magnesiowüstite, as did Wang et al. (2006). The latter found that the large velocity jumps (in NE Asia) may involve a larger fraction of garnet transforming to perovskite. They also explored how these compositions could affect the velocity gradients surrounding the discontinuity. Our normal mode study cannot resolve these gradients. Indeed, while our new estimates bring PREM closer to mineral physics estimates for model mantle minerals, these still lie on the higher end of estimates for adiabatic pyrolite (Cammarano et al., 2005).

## 5. Conclusion

We have used the Backus-Gilbert method to find a combination of normal mode center frequency data that enhances sensitivity to just above and below the 660 discontinuity, for density, P-wave speeds, S-wave speeds. We have determined the best-fitting ranges of jumps in these parameters when assuming PREM as a background model (Figure 2b). There is significant covariance between these parameters. The corresponding PREM value for our $\Delta\rho$ lies above this range, the $\Delta v_p$ lies below this range, and $\Delta v_s$ lies within this range. When shifting the depth of 660 to 660 km brings out additional trade-offs, resulting in a range of acceptable models that span larger values of $\Delta\rho$ and $\Delta v_s$, and smaller values of $\Delta v_p$. In these calculations, we also found that the normal mode data do not support a globally averaged phase transition depth as shallow as 650 km depth.

Our results produce a range of values for $\Delta\rho$ and $\Delta v_b$ that are generally higher than those estimated by mineral physicists for the pyrolite model, and in particular even higher than PREM for $\Delta v_b$, supporting the possibility of a larger proportion of garnet in the transformation to perovskite. Still, the density jump of PREM is at the high end of the acceptable models resulting from our study, which may be important for geodynamicists modeling global convection. Finally, the inability to obtain a consistent result when using ak135 as a reference model may reflect frequency dependence of structure and/or the presence of significant lateral variations around the 660.

## Data Availability Statement

Compiled normal mode center frequencies data are available through Deuss et al. (2013), Roult et al. (2010), and Laske, Dziewonski, & Masters (2001).

## Acknowledgments

The authors thank two anonymous reviewers for providing constructive comments on the original manuscript. This work was supported by NSF grant EAR-1923865. They thank Dr Alexander Robson for sharing the compiled normal mode center frequencies (where the data are available through Deuss et al., 2013, Roult et al., 2010, and Laske, Dziewonski, & Masters, 2001). H. Lau thanks David Al-Attar (University of Cambridge) for some fruitful discussions related to Backus-Gilbert methods.

## References

* Allgre, C. J. (1997). Limitation on the mass exchange between the upper and lower mantle: The evolving convection regime of the earth. *Earth and Planetary Science Letters*, 150, 1–6. https://doi.org/10.1016/S0012-821X(97)00072-1

* Andrews, J., & Deuss, A. (2008). Detailed nature of the 660 km region of the mantle from global receiver function data. *Journal of Geophysical Research*, 113, B06304. https://doi.org/10.1029/2007JB005111

* Backus, G., & Gilbert, F. (1970). Uniqueness in the inversion of inaccurate gross Earth data. *Philosophical Transactions of the Royal Society of London Series A*, 266, 123–192. https://doi.org/10.1098/rsta.1970.0005

* Bass, J. D., & Anderson, D. L. (1984). Composition of the upper mantle: Geophysical tests of two petrological models. *Geophysical Research Letters*, 11, 229–232. https://doi.org/10.1029/GL011i003p00229

* Birch, F. (1952). Elasticity and constitution of the Earth's interior. *Journal of Geophysical Research*, 57, 227–286. https://doi.org/10.1029/JZ057i002p00227

* Cammarano, F., Deuss, A., Goes, S., & Giardini, D. (2005). One-dimensional physical reference models for the upper mantle and transition zone: Combining seismic and mineral physics constraints. *Journal of Geophysical Research*, 110, B01306. https://doi.org/10.1029/2004JB003272

* Castle, J. C., & Creager, K. C. (2000). Local sharpness and shear wave speed jump across the 660-km discontinuity. *Journal of Geophysical Research*, 105, 6191–6200. https://doi.org/10.1029/1999JB900424

* Dahlen, F. A., & Tromp, J. (1998). *Theoretical global seismology* (p. 1025). Princeton University Press.

* Deuss, A. (2009). Global observations of mantle discontinuities using SS and PP precursors. *Surveys in Geophysics*, 30, 301–326. https://doi.org/10.1007/s10712-009-9078-y

* Deuss, A., Redfern, S. A. T., Chambers, K., & Woodhouse, J. H. (2006). The nature of the 660-kilometer discontinuity in Earth's mantle from global seismic observations of PP precursors. *Science*, 311, 198–201. https://doi.org/10.1126/science.1120020

* Deuss, A., Ritsema, J., & van Heijst, H. (2013). A new catalogue of normal-mode splitting function measurements up to 10 mHz. *Geophysical Journal International*, 193, 920–937. https://doi.org/10.1093/gji/ggt010

* Dziewonski, A. M., & Anderson, D. L. (1981). Preliminary reference Earth model. *Physics of the Earth and Planetary Interiors*, 25, 297–356. https://doi.org/10.1016/0031-9201(81)90046-7
Estabrook, C. H., & Kind, R. (1996). The nature of the 660-kilometer upper-mantle seismic discontinuity from precursors to the PP phase. *Science*, **274**, 1179–1182. https://doi.org/10.1126/science.274.5290.1179

Frost, D. J. (2008). The upper mantle and transition zone. *Elements*, **4**, 171–176. https://doi.org/10.2113/GSELEMENTS.4.3.171

Ishii, T., Kojitani, H., & Akaogi, M. (2018). Phase relations and mineral chemistry in pyrolitic mantle at 1600-2200°C under pressures up to the uppermost lower mantle: Phase transitions around the 660-km discontinuity and dynamics of upwelling hot plumes. *Physics of the Earth and Planetary Interiors*, **274**, 127–137. https://doi.org/10.1016/j.pepi.2017.10.005

Kennett, B. L. N., & Engdahl, E. R. (1991). Traveltimes for global earthquake location and phase identification. *Geophysical Journal International*, **105**, 429–465. https://doi.org/10.1111/j.1365-246X.1991.tb06724.x

Kennett, B. L. N., Engdahl, E. R., & Buland, R. (1995). Constraints on seismic velocities in the Earth from traveltimes. *Geophysical Journal International*, **122**, 108–124. https://doi.org/10.1111/j.1365-246X.1995.tb03540.x

Kustowski, B., Ekström, G., & Dziewoński, A. M. (2008). Anisotropic shear-wave velocity structure of the Earth's mantle: A global model. *Journal of Geophysical Research*, **113**, B06306. https://doi.org/10.1029/2007JB005169

Laske, G., Dziewonski, A., & Masters, G. (2001). Reference earth model [data set]. Retrieved from https://igppweb.ucsd.edu/~gabi/rem.dir/surface/smodes.list

Masters, G., & Gubbins, D. (2003). On the resolution of density within the Earth. *Physics of the Earth and Planetary Interiors*, **140**(1-3), 159-167. https://doi.org/10.1016/j.pepi.2003.07.008

Matsui, M. (2000). Molecular dynamics simulation of MgSiO₃ perovskite and the 660-km seismic discontinuity. *Physics of the Earth and Planetary Interiors*, **121**, 77-84. https://doi.org/10.1016/S0031-9201(00)00161-8

Montagner, J.-P., & Anderson, D. L. (1989). Constrained reference mantle model. *Physics of the Earth and Planetary Interiors*, **58**, 205-227. https://doi.org/10.1016/0031-9201(89)90055-1

Morelli, A., & Dziewonski, A. M. (1993). Body wave traveltimes and a spherically symmetric P- and S-wave velocity model. *Geophysical Journal International*, **112**, 178-194. https://doi.org/10.1111/j.1365-246X.1993.tb01448.x

Niu, F., & Kawakatsu, H. (1995). Direct evidence for the undulation of the 660-km discontinuity beneath tonga: Comparison of Japan and California array data. *Geophysical Research Letters*, **22**, 531-534. https://doi.org/10.1029/94g03332

Pijpers, F. P., & Thompson, M. J. (1992). Faster formulations of the optimally localized averages method for helioseismic inversions. *Astronomy and Astrophysics*, **262**, L33-L36.

Ringwood, A. E. (1962). A model for the upper mantle. *Journal of Geophysical Research*, **67**, 857-867. https://doi.org/10.1029/JZ067i002p00857

Ringwood, A. E. (1991). Phase transformations and their bearing on the constitution and dynamics of the mantle. *Geochimica et Cosmochimica Acta*, **55**, 2083-2110. https://doi.org/10.1016/0016-7037(91)90090-R

Robson, A. J. S., & Romanowicz, B. (2019). New normal mode constraints on bulk inner core velocities and density. *Physics of the Earth and Planetary Interiors*, **295**, 106310. https://doi.org/10.1016/j.pepi.2019.106310

Roult, G., Roch, J., & Clévédé, E. (2010). Observation of split modes from the 26th December 2004 Sumatra-Andaman mega-event. *Physics of the Earth and Planetary Interiors*, **179**, 45-59. https://doi.org/10.1016/j.pepi.2010.01.001

Shearer, P. M. (2000). Upper mantle seismic discontinuities. Geophysical monograph series (117, pp. 115-131). Washington DC: American Geophysical Union. https://doi.org/10.1029/GM117p0115

Shearer, P. M., & Flanagan, M. P. (1999). Seismic velocity and density jumps across the 410- and 660-kilometer discontinuities. *Science*, **285**, 1545-1548. https://doi.org/10.1126/science.285.5433.1545

Shim, S.-H., Duffy, T. S., & Shen, G. (2001). The post-spinel transformation in Mg2SiO4 and its relation to the 660-km seismic discontinuity. *Nature*, **411**, 571-574. https://doi.org/10.1038/35079053

Stixrude, L., & Lithgow-Bertelloni, C. (2005). Thermodynamics of mantle minerals - I. Physical properties. *Geophysical Journal International*, **162**, 610-632. https://doi.org/10.1111/j.1365-246X.2005.02642.x

Vacher, P., Mocquet, A., & Sotin, C. (1998). Computation of seismic profiles from mineral physics: The importance of the non-olivine components for explaining the 660 km depth discontinuity. *Physics of the Earth and Planetary Interiors*, **106**, 275-298. https://doi.org/10.1016/S0031-9201(98)00076-4

van Keken, P., & Ballentine, C. (1998). Whole-mantle versus layered mantle convection and the role of a high-viscosity lower mantle in terrestrial volatile evolution. *Earth and Planetary Science Letters*, **156**(1), 19-32. https://doi.org/10.1016/s0012-821x(98)00023-5

Wang, Y., Wen, L., Weidner, D., & He, Y. (2006). SH velocity and compositional models near the 660-km discontinuity beneath South America and northeast Asia. *Journal of Geophysical Research*, **111**, B07305. https://doi.org/10.1029/2005JB003849

Woodhouse, J. H., & Dahlen, F. A. (1978). The effect of a general aspherical perturbation on the free oscillations of the Earth. *Geophysical Journal International*, **53**, 335-354. https://doi.org/10.1111/j.1365-246X.1978.tb03746.x

Xu, F., Vidale, J. E., & Earle, P. S. (2003). Survey of precursors to P'P': Fine structure of mantle discontinuities. *Journal of Geophysical Research*, **108**(B1), 7-1. https://doi.org/10.1029/2001JB000817