Increased density of large low-velocity provinces recovered by seismologically constrained gravity inversion

The nature and origin of the two large low-velocity provinces (LLVPs) in the lowest part of the mantle remain controversial. These structures have been interpreted as a purely thermal feature, accumulation of subducted oceanic lithosphere or a primordial zone of iron enrichment. Information regarding the density of the LLVPs would help to constrain a possible explanation.

In this work, we perform a density inversion for the entire mantle, by constraining the geometry of potential density anomalies using tomographic vote maps. Vote maps describe the geometry of potential density anomalies according to their agreement with multiple seismic tomographies, hence not depending on a single representation. We use linear inversion and determine the regularization parameters using cross-validation. Two different input fields are used to study the sensitivity of the mantle density results to the treatment of the lithosphere. We find the best data fit is achieved if we assume that the lithosphere is in isostatic balance.

The estimated densities obtained for the LLVPs are systematically positive density anomalies for the LLVPs in the lower 800–1000 km of the mantle, which would indicate a chemical component for the origin of the LLVPs. Both iron-enrichment and a mid-oceanic ridge basalt (MORB) contribution are in accordance with our data, but the required superadiabatic temperature anomalies for MORB would be close to 1000 K.

1 Introduction

Seismology has systematically revealed more and more of the heterogeneity in the mantle. Interpretations of seismic images often aim at determining the density of features, because buoyancy is the driver of mantle dynamics. For instance, much research has focused on tracing sinking slabs in order to understand how subduction works (e.g. van der Meer et al., 2018). At greater depths such interpretation becomes more difficult, due to reduced resolution and lack of corroborating surface evidence. The most striking features of the lowermost mantle are probably the large low-velocity provinces (LLVPs) (Garnero et al., 2016).

The two LLVPs are antipodal regions of decreased seismic velocity that extend from the core–mantle boundary about 400–800 km into the mantle and cover about 25 % of the surface of the core–mantle boundary (Garnero et al., 2016). There is considerable debate regarding their nature and origin, with them being interpreted as a purely thermal feature (Koelemeijer et al., 2017), an accumulation of subducted oceanic lithosphere (Mulyukova et al., 2015) or iron-enriched zones (Ballmer et al., 2015; Deschamps et al., 2012). In addition, it has been suggested that plumes are formed preferentially at the edges of the LLVPs, based on the reconstructed locations of large igneous provinces and kimberlites (Burke et al., 2008; Conrad et al., 2013), which imply that the LLVPs are stable over long time periods.

The debate regarding the LLVPs is often centred on their density, because density allows one to distinguish between a purely thermal and a thermo-chemical origin. However, direct determinations of density using seismological and geodetic methods have led to contradictory results, indicating either positive (Ishii and Tromp, 1999; Lau et al., 2017; Moulik and Ekström, 2016; Trampert et al., 2004) or negative (Koelemeijer et al., 2017) density anomalies, with a strength between 0.5 % and 2.0 %.

Combining seismology and gravity to investigate the mantle is a well-established approach that typically includes dynamical modelling of mantle convection (Richards and Hager, 1984; Hager and O'Connell, 1979). Most commonly, seismic velocity anomalies are directly converted into density anomalies, and the vertical viscosity distribution is adjusted in order to fit the geoid (Richards and Hager, 1984). Using this technique, many models of viscosity and density inside the mantle have been derived (King, 1995, 2016; Steinberger, 2007, 2016; Richards and Hager, 1984; Hager et al., 1985).

While viscosity inversions can successfully explain the geoid with variance reductions of 80 % or more (King and Masters, 1992; Ricard et al., 1993), this approach relies on some assumptions that can be challenged. To derive a conversion factor from velocity to density anomalies based on mineral physics, often purely thermal anomalies are assumed (Karato, 1993). However, differences in terms of composition certainly play a role in the upper mantle/lithosphere (Griffin et al., 2009) and have been suggested for the lower mantle as well (Ballmer et al., 2017; Deschamps et al., 2012). For example, the velocity and density model Gypsum (Simmons et al., 2009) allows for a partial decorrelation of velocity and density.

Furthermore, the modelled geoid associated with a certain density distribution is partly based on the dynamic topography caused by this density distribution. Thus, viscosity inversion should also reproduce surface topography, since it is an important part of the total gravity effect. However, before the dynamic predictions can be compared to real topography, the isostatic part of topography must be removed using crustal models. While the agreement between predicted and non-isostatic topography has been increasing (Flament et al., 2013; Steinberger, 2016), there are still considerable differences (Molnar et al., 2015).

The justification for directly inverting for density instead of going through viscosity is that the gravity field is only affected indirectly by the viscosity distribution, through viscosity's influence on the deformation of the boundaries (surface and core–mantle boundary) (Hager and O'Connell, 1979). However, the gravity effect of any deformed boundary is independent of what caused the deformation. Thus, the gravity effect of the top boundary (i.e. topography) can simply be calculated and removed from the observed data, without needing any dynamic simulations or viscosity structures. In contrast, the core–mantle boundary (CMB) topography is not well known. While a variety of seismological techniques can be used to estimate CMB topography, there is no consensus regarding amplitude and patterns (see comparisons in Tanaka, 2010). Geodynamical simulations have shown that CMB deformation is sensitive to possible lateral viscosity variations in the lowermost mantle and the viscosity and density variations related to the presence post-perovskite (Deschamps et al., 2018; Deschamps and Li, 2019). Hence, CMB deformation cannot be accounted for a priori, but its possible effects must be considered a posteriori in the interpretation.

In this contribution, we explore an alternative approach for fitting the gravity field using mantle density structure. Instead of density–velocity conversion, we use tomographic vote maps (Lekic et al., 2012; Shephard et al., 2017) to provide geometrical constraints for a linear gravity inversion. The geometry is defined by the fast and slow zone derived from the vote maps, and each zone is assigned an unknown but constant density anomaly. These density anomalies are completely free parameters that are adjusted to fit the gravity field. As input data we use satellite gravity data, corrected for topography and lithosphere. For the latter correction we compare two different approaches. In the first approach it is assumed that the lithosphere is in perfect isostatic balance, and in the second we make use of a recently developed crustal model by Szwillus et al. (2019).

2 Methods and data

2.1 Gravity data

We use satellite gravity data from the global gravity field model GOCO05S (Pail et al., 2010) at the measurement height (225 km) of the GOCE satellite in its last mission phase (van der Meijde et al., 2015). Often, gravity data are mathematically downward continued to a lower height, but this introduces omission errors (Bouman et al., 2013), which we would like to avoid. In addition, the target density distributions are at such a great depth that downward continuation probably provides no additional sensitivity.

As a first step we calculate an ice-corrected Bouguer anomaly based on the topography of ETOPO1 (Amante and Eakins, 2009). Onshore we assume a density of 2670 kg m⁻³, and offshore we use a correction density of 1770 kg m⁻³ based on a reference crustal density of 2800 kg m⁻³ and a water density of 1030 kg m⁻³. In addition, the influence of ice is removed assuming a density of 917 kg m⁻³.

Since we are interested in the mantle density structure, the second step is to account for the gravity effect of the crust and/or lithosphere. Two inherently different approaches exist to achieve this. The first option is to assume some form of isostatic compensation and use this to account for the crust and/or lithosphere. Alternatively, the gravity effect of crustal models, such as Crust1.0 (Laske et al., 2013) or Litho1.0 (Pasyanos et al., 2014) can be calculated. The latter approach includes more geophysical data, but it also “propagates” errors of the crustal model into the mantle. In particular, all crustal models predict several kilometres of residual topography, even in the continents (Steinberger, 2016) and are thus not nearly in isostatic balance. As a consequence, the topography would need to be supported dynamically, which should lead to higher free-air anomaly values than observed (Molnar et al., 2015).

In light of this discrepancy, we explore two approaches. In the first approach, we assume that continental topography is compensated by crustal thickness variations, whereas ocean floor depth is compensated by variations in the mantle lithosphere density (the same isostatic model as in Szwillus et al., 2016, with a compensation depth of 120 km). The isostatic residual is then

$$\begin{array}{}\text{(1)}& \mathrm{isostatic}\mathrm{residual}=\mathrm{FA}-g_{\mathrm{topo}}-g_{\mathrm{iso}},\end{array}$$

where FA is the free-air anomaly (Fig. 1a), g_topo is the gravity effect of topography and g_iso is the gravity effect of the isostatic compensation masses.

Figure 1 Input data for the inversions. (a) Free-air anomaly at a measurement height of 225 km. (b) Isostatic residual gravity anomaly at 225 km height. This is obtained by correcting for the presence of topography and isostatic compensation masses. (c) Crustal residual gravity anomaly at 225 km height obtained by removing the gravity effect of the crustal model from Szwillus et al. (2019). (d) Residual topography after removing the isostatic contribution of the crustal model from Szwillus et al. (2019).

In the second approach, we use our recent crustal model based on kriging (Szwillus et al., 2019) and assume a density contrast between crust and mantle of 400 kg m⁻³ to estimate the gravity effect of varying crustal thickness. This gives the crustal residual:

$$\begin{array}{}\text{(2)}& \mathrm{crustal}\mathrm{residual}=\mathrm{FA}-g_{\mathrm{topo}}-g_{\mathrm{crust}},\end{array}$$

where g_crust is the gravity response of the crustal thickness variations.

The isostatic residual is very similar to the free-air anomaly (Fig. 1b), except in areas of high topography. In contrast, the crustal residual has a much higher amplitudes and patterns (Fig. 1c). For comparison we also computed the residual (non-isostatic) topography of the kriging crustal model (Szwillus et al., 2019), which anti-correlates with the crustal residual on large scales (Fig. 1d).

2.2 Geometry constraints

Additional constraints are required to overcome the inherent non-uniqueness of a density inversion. We use information from whole mantle seismic tomographies in the form of vote maps (Shephard et al., 2017) to define volumes which could potentially contain density anomalies. Vote maps were first introduced by Lekic et al. (2012) in the context of cluster analysis and then further developed by Shephard et al. (2017). A vote map is based on a collection of tomography models, and at each point it gives the number of tomographies that detect a significant anomaly at that location. The threshold that constitutes a significant anomaly is the standard deviation at that depth for each seismic tomography. Thus, even tomographies with vastly different amplitude ranges can be combined in a single vote map. We used the vote maps based on 17 S-wave tomography models as available on SubMachine (Hosseini et al., 2018) at a resolution of 1^∘ and 100 km depth spacing. There are two separate vote maps for positive and negative anomalies. Figure 2 shows an example of the vote maps at a depth of 2700 km, close to the core–mantle boundary.

Figure 2 Vote maps at 2700 km depth. (a) Vote map for fast (positive) anomalies. (b) Vote map for slow (negative) anomalies. Note the clear signature of the two LLVPs underneath Africa and the central Pacific.

The first step is to discretize the inversion problem by extracting regions of potential density anomalies from the vote map individually for each depth. Our underlying assumption is that density anomalies only occur in horizontally connected regions where at least N tomographies detect anomalous velocity. We choose N=10, because this is the cut-off that correlates with the strongest gradients of the vote maps. Each connected region is assigned a constant but unknown density anomaly. As a result of this process, we have a total of 1135 anomalous regions that are each described in terms of their depths, their pixels and unknown density anomaly values. These density anomalies are completely free parameters that will be adjusted during the inversion to fit the gravity field. Importantly, we impose no correlation between velocity and density.

2.3 Cross-validation inversion approach

Our inversion is a two-step approach. First, we forward calculate for each potential density anomaly region its characteristic gravity response. Second, the characteristic gravity responses of all of the potential density anomalies are placed into a matrix A and solved with a generalized Tikhonov regularization technique. The regularization parameters are chosen based on cross-validation.

The gravity effect of each potential density anomaly is estimated by approximating it as a collection of point masses. Each point mass represents a volume V, given by

$$\begin{array}{}\text{(3)}& V=r^{\mathrm{2}}\cos \left(\mathrm{lat}\right)(\mathrm{\Delta }\mathit{\varphi }{)}^{\mathrm{2}}\mathrm{\Delta }z,\end{array}$$

where r is the radius, lat is the latitude, Δϕ is the angular size of each pixel (1^∘) and Δz is the vertical resolution (100 km). The mass m is then found by multiplying the volume with the density anomaly ρ.

The gravity kernel K(P,Q) caused by a unit point mass located at $P=(r,\mathrm{lat},\mathrm{long})$ on a measurement point $Q=(r_i,{\mathrm{lat}}_i,{\mathrm{long}}_i)$ is given by

$$\begin{array}{}\text{(4)}& K(P,Q)=G{\displaystyle \frac{r_{\mathrm{s}}-r\cos \mathit{\gamma }}{{\left(r^{\mathrm{2}}+r_{\mathrm{s}}^{\mathrm{2}}-\mathrm{2}r{r}_{\mathrm{s}}\cos \mathit{\gamma }\right)}^{\mathrm{3}/\mathrm{2}}}},\end{array}$$

where γ is the great-circle distance between the point mass and the measurement point, and $G=\mathrm{6.67428}\times {\mathrm{10}}^{-\mathrm{11}}$ m³ (kg s²)⁻¹ is the gravitational constant. The total gravity effect of each potential density anomaly is then found by summing over the gravity effect of all its pixels.

Next, the gravity response of all of the potential density anomalies are placed into a matrix A, where A_ij is the gravity effect of potential density anomaly j on measurement point i, assuming a density of 1. Using Eq. (4), A_ij becomes

$$\begin{array}{}\text{(5)}& A_{ij}=\sum _{k\in {\mathcal{R}}_j}K(P_k,Q_i)V_k,\end{array}$$

where ℛ_j is the set of indices belonging to potential density anomaly j; P contains all subsurface points and Q all measurement points.

The unknown density values of the potential density anomalies are placed in a vector ρ. In principle, $\mathit{\rho }={\mathbf{A}}^{-\mathrm{1}}\mathit{g}$, where g_i=g(Q_i) is the measured gravity value at point Q_i. However, in practice the solution will be unstable and non-unique, so that regularization is required. Here, we choose to minimize the following least-squares functional:

$$\begin{array}{}\text{(6)}& \mathit{\chi }(\mathit{\rho }{)}^{\mathrm{2}}=||\mathbf{A}\mathit{\rho }-\mathit{g}|{|}_{{\mathrm{\Sigma }}_{\mathrm{g}}^{-\mathrm{1}}}^{\mathrm{2}}+\mathit{\beta }\left|\right|\mathit{\rho }|{|}^{\mathrm{2}}+\mathit{\gamma }|\left|\mathbf{D}\mathit{\rho }\right|{|}^{\mathrm{2}}.\end{array}$$

The three terms relate to the data misfit, the magnitude of density variations and the vertical derivative of density anomaly. $\left|\right|\cdot |{|}_M$ is the quadratic form related to matrix M, and $\left|\right|\cdot \left|\right|=\left|\right|\cdot |{|}_I$ is the Euclidean norm. β and γ are the regularization parameters that enforce the minimization of absolute norm of the density anomalies and their vertical derivative respectively. In addition, the parameter β represents the ratio of data variance to model variance. D is a finite differencing matrix consisting of first-order forward vertical differences in density. Finally, Matrix Σ_g is an assumed correlation matrix of the gravity values. We use an isotropic correlation function with a correlation distance of 10^∘, based on a semivariogram analysis (Chilès and Delfiner, 2012). The main purpose of the matrix is to down-weigh observation points near the pole, which are otherwise over represented due to convergence near the pole.

The optimal ρ belonging to this formulation is

$$\begin{array}{}\text{(7)}& \widehat{\mathit{\rho }}={\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{\Sigma }}_{\mathrm{g}}^{-\mathrm{1}}\mathbf{A}+\mathit{\beta }\mathbf{I}+\mathit{\gamma }{\mathbf{D}}^{\mathrm{T}}\mathbf{D}\right)}^{-\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}{\mathbf{\Sigma }}_{\mathrm{g}}^{-\mathrm{1}}\mathit{g}.\end{array}$$

The regularization parameters β and γ exert critical control over the resulting density structure and the achievable data fit. The most popular methods to determine such parameters are L curves and cross-validation (Farquharson and Oldenburg, 2004). Here we use cross-validation to estimate the regularization parameter.

In k-fold cross-validation, the gravity data set is randomly split into k distinct index sets ℐ_j. The first set ℐ₁=𝒯 then becomes the training set and the remaining sets are combined to form the validation set 𝒱. The training set is then used to solve the inversion problem giving a solution ρ₁. Based on how well ρ₁ predicts the data in 𝒯 and 𝒱, the training and validation misfit is calculated. Each ℐ_𝒿 becomes the training set once, and this is repeated for many random partitions of the data.

By design, the training misfit is minimal when β and γ are zero (if there are no numerical instabilities). However, small values of β and γ also risk overfitting the data. By carrying out cross-validation for all combinations of β and γ in a reasonable range of values, the optimal β and γ corresponding to the minimum validation misfit can be found.

Apart from constraining the regularization parameters, this procedure also gives a bootstrap estimate (Efron and Gong, 1983) of the density structure. The collection of all the density anomaly estimates derived for all subsets of gravity data can be used to calculate a mean density anomaly value and its standard deviation for each potential anomaly, given certain regularization parameters. This is useful to determine how robust a single density anomaly is.

2.4 Boundary topography

The surface and CMB deformations implied by the density model recovered by inversion are important to judge the quality of the results. The surface topography predicted by the density model should roughly agree with the residual topography, because topography is strongly sensitive to density. The magnitude of the CMB deformation is required to assess how much it affects the gravity field fitting.

We determine the deformation of the upper and lower domain boundary caused by density structure as a postprocessing step after inversion. We use both an isostatic and a simple dynamic formulation to determine this deformation.

Let t_CMB(lat,long) be the CMB-topography as a function of longitude and latitude. Under the assumption of isostatic balance at the CMB, we find

$$t_{\mathrm{CMB}}(\mathrm{lat},\mathrm{long})=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{CMB}}}}\underset{\mathrm{CMB}}{\overset{\mathrm{CMB}+H}{\int }}\mathit{\delta }\mathit{\rho }(r,\mathrm{lat},\mathrm{long})\mathrm{d}r,$$

where Δρ_CMB is the density jump at the CMB, H is the height above the CMB where isostatic balance applies, and $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ is the mean density between r=CMB and $r=\mathrm{CMB}+H$.

We convert the density difference associated with undulating topography into a surface density σ=t_CMBΔρ_CMB and calculate its gravity response in spherical harmonics by

$$g=\sum _{l,m}{g}_l^m{Y}_l^m=\sum _{l,m}\mathrm{4}\mathit{\pi }G{\displaystyle \frac{l+\mathrm{1}}{\mathrm{2}l+\mathrm{1}}}{\left({\displaystyle \frac{r_{\mathrm{CMB}}}{r_{\mathrm{s}}}}\right)}^{l+\mathrm{2}}{\mathit{\sigma }}_l^m{Y}_l^m,$$

where ${\mathit{\sigma }}_l^m$ values are the spherical harmonic coefficients of the surface density, $Y_l^m$ is a spherical harmonic function and $r_{\mathrm{s}}=\mathrm{6371}+\mathrm{225}$ km is the radius of the observations.

Isostatic surface topography is calculated in the same way using the upper 300 km of the domain and a density contrast of 2670 kg m⁻³ at the surface.

To determine the dynamic surface and CMB topography we use the propagator matrix approach of Hager and O'Connell (1979), with free slip boundary conditions on the surface and CMB. The viscosity is layered and increases by a factor of 100 at 700 km depth. This is a simplification of the real viscosity structure, which is additionally affected by lateral variations due to temperature and – possibly – compositional differences.

The resulting kernels show that for the assumed viscosity structure, the influence of density anomalies below 500 km on topography is limited, even at very long wavelengths (Fig. 3a). The CMB kernels decrease much more slowly, even with a layered viscosity (Fig. 3b) and even at degree 10.

Figure 3 Topography and CMB kernels calculated using the propagator matrix approach (Hager and O'Connell, 1979). The black line is computed using a two-layer viscosity structure, the grey line uses a constant viscosity. The solid line is for spherical harmonic degree 2 and the dashed line for spherical harmonic degree 10.

3 Results

3.1 Inversion using isostatic residual

A grid search of β and γ values between 10⁻³ and 10³ using the cross-validation procedure leads to a minimum validation misfit of 5.32 mGal if β=0.03 and γ=1.0 (Fig. 4a). Although this minimum is not particularly well developed, we proceed with these values of β and γ, since neighbouring regularization settings give similar results. With these regularization parameters, the training misfit becomes 3.25 mGal, which is about 20 % of the rms of the input gravity field (96 % variance reduction).

Figure 4 Results of the inversion using the isostatic residual. (a) Validation misfit as a function of regularization parameter β. (b) Predicted field of the inversion results. (c) Spectral comparison of the input and predicted data. (d) Difference between input and predicted data.

The inversion is able to reasonably reproduce the main features of the isostatic residual gravity field (Fig. 4b). However, the amplitude of the residual field is systematically too small. A linear regression of the gravity field predicted from the model against the measured values gives a linear factor of 0.85, implying that the modelled anomalies are about 15 % too weak. This is probably a result of the regularization.

In the spectral domain, the predicted gravity field is able to reproduce the long-wavelength part up to spherical harmonic (SH) degree 10 (Fig. 4c). At shorter wavelengths, the predicted field is systematically too weak, going down to less than 10 % of the observed data at spherical harmonic degree 40. For comparison, the solution with no regularization ($\mathit{\beta }=\mathit{\gamma }=\mathrm{0}$) is able to reproduce the spectrum up to spherical harmonic degree 25 perfectly. Thus, wavelengths shorter than SH degree 25 simply cannot be produced, due to the insufficient spectral content of the vote maps.

The inverted density anomalies show no clear correlation with sign of the velocity anomaly (Sect. 4.2), which would suggest that many density anomalies are not caused only by temperature variations. The spatial distribution is highly complex, and not all density anomalies are easily related to expected features. However, some first-order trends can be identified (Fig. 5).

Figure 5 Density anomaly slices through selected depths of the inversion results using the isostatic residual.

The upper 100 km of the inversion results contain very little density variations (Fig. 5a). Some cratons (South America, South Africa, Eurasia) show a slight negative density anomaly (Fig. 5b). This is probably related to simplifications in our isostatic model and not representative of the actual density structure

At greater depths, subducted slabs would be expected in the density structure. This is overall reflected in our density inversion results, but not all slabs are resolved as positive density anomalies over their complete depth range. For example, the Andean slabs can be seen at a depth of 900 km as an anomaly of +6 kg m⁻³ (Fig. 5d) but not at 600 km depth.

The LLVPs have systematically a positive density anomaly of up to +10 kg m⁻³ (Fig. 5e and f). The top of the African LLVP appears at a depth of 1800 km and extends down to the CMB. The Pacific LLVP only appears at greater depths of 2200 km and is also continuous to the CMB. Both LLVPs are overlaid by slight negative density anomalies.

3.2 Inversion using crustal residual

The preferred regularization parameters derived by the cross-validation procedure are similar to those for the isostatic residual, β=0.1 and γ=1.0 (Fig. 6a). With these values, a validation misfit of 34 mGal and a training misfit of 23 mGal are achieved, which corresponds to about 40 % of the input signal. Thus, the fit for the crustal residual is considerably worse than for the isostatic residual, both in an absolute and a relative sense. The predicted and input gravity data only agree qualitatively in some areas (Fig. 6b), and large residuals remain (Fig. 6d).

Figure 6 Results of the inversion using the crustal residual. (a) Validation misfit as a function of regularization parameter β. (b) Predicted field of the inversion results. (c) Spectral comparison of the input and predicted data. (d) Difference between input and predicted data.

Overall the inverted density variations are much larger for the crustal residual than for the isostatic residual (Fig. 7). In the upper 300 km or so, this is expected, since the lithospheric density structure has not been accounted for by the crustal model. The ocean cooling trend is not imaged well by the vote maps, and as a result, only the Pacific and Indian oceans' spreading ridges are resolved as crude negative density anomalies (Fig. 7a). Some cratons (North America, Eurasia) show strong positive density anomalies (30 kg m⁻³), but this neither is consistent vertically nor applies to most geographical regions (Fig. 7a and b). Thus, the ability to resolve mantle lithospheric sources is limited and probably contributes to the poor fit of the data.

Figure 7 Density anomaly slices through selected depths of the inversion results using the crustal residual.

At depths of 600 and 900 km, some subduction-related density anomalies can be seen, but overall correlation is poor. Only in the lowest mantle a more consistent picture emerges, in particular the LLVPs are resolved as positive density anomalies (+40 kg m⁻³). As in the case with isostatic residual, the African LLVP appears to extend further upwards as a positive density anomaly than the Pacific LLVP (Fig. 7e and f).

4 Discussion

4.1 Resolved mantle structure

While the main features observed are the LLVPs, it is important to also scrutinize the results for the remaining mantle. A long-wavelength error in the upper-mantle density structures leads to an incorrect estimate of the density of the LLVPs. In the upper mantle above the transition zone (410–670 km depth) much of the seismic image can be linked to tectonic features: thick and cold lithospheric roots underlie the cratonic cores of the continents and can extend up to 300 km deep (Schaeffer and Lebedev, 2015); the oceans are characterized by a velocity trend related to cooling of the lithosphere (Priestley and McKenzie, 2006); and slabs formed during subduction sink through the mantle, sometimes stagnating at various depths and sometimes passing through the transition zone to the lowermost mantle (e.g. van der Meer et al., 2018).

The cratonic lithospheric mantle density is the result of a complex interaction between temperature, composition and mineralogy (Griffin et al., 2009; Fullea et al., 2009; Afonso et al., 2013). The density increase due to low temperature is probably (partially) offset by its depleted composition (the isopycnic hypothesis; Jordan, 1978). Indeed, a slightly positive density anomaly (+5 to 7.5 kg m⁻³) is recovered by the inversion with crustal residual in the cratonic cores of North America, central Europe, Siberia and East Antarctica. In contrast, the inversion with the isostatic residual does not show the cratons. The cratonic mantle keels affect surface topography and, hence, are removed by the isostatic correction.

The density of the oceanic mantle is expected to increase with age, following the cooling trend, which also explains seafloor subsidence. The inversion with crustal residual shows some decreased density associated with mid-oceanic ridges, but vote maps are unsuited to resolve the continuous gradient associated with the cooling of the plate. When using the isostatic residual, ocean floor subsidence and the cooling trend has already been accounted for during the isostatic correction and is thus not resolved by the inversion.

Subducted slabs are denser than surrounding mantle immediately after subduction. As the plate sinks, heat diffuses into the plate from the surrounding mantle, so that the plate slowly loses its negative buoyancy. How quickly a subducted plate thermally equilibrates depends on how effectively heat is transported to the slab by convection and conduction. Since some slabs seem to stagnate on mantle discontinuities, it is conceivable for a slab to loose all or most of its negative buoyancy. We compared our results with the slab depth contours from Slab 1.0 (Hayes et al., 2012).

Out of the 12 slabs contained in Slab1.0, only five are detected as positive velocity anomalies based on the vote maps (Fig. 8a) and only at depths greater than 200 km or more. These are Izu–Bonin, Kermadec–Tonga, Kamchatka–Kurils–Japan, South America and Sumatra, according to the terminology of Slab1.0.

Figure 8 Comparison with Slab1.0 (Hayes et al., 2012). The slab codes stand for the following: alu – Alaska–Aleutians; izu – Izu–Bonin; ker – Kermadec–Tonga; kur – Kamchatka–Kurils–Japan; mex – Central America; phi – Philippines; ryu – Ryukyu; sam – South America; sco – Scotia (South Sandwich); sol – Solomon Islands; sum – Sumatra; van – Santa Cruz Islands–Vanuatu–Loyalty Islands. (a) Average value of the velocity anomalies inside the Slab1.0 contour, where 1 = positive, −1 = negative anomaly. (b) Average density anomaly from inversion using the isostatic residual. (c) Average density anomaly from inversion of the crustal residual.

The reason might be that the tomography models are too heterogeneous in quality and resolution to resolve these relatively narrow features. Of course, vote maps are not the only way to extract features from seismic tomographies, and an alternative would be to use the “best” seismic tomography model available instead of a collection of models. But there are no clear criteria to decide which tomography result is the best, except that newer results are probably based on more and better data. Furthermore, there is some evidence that the differences between tomographic models partly reflect different regularization parameters, rather than different resolved features (Root et al., 2016; Root, 2020). Thus, vote maps only extract the features that are resolved by a majority of tomographic models

Still, out of the five detected slabs, four have positive density contrasts (+5 to +10 kg m⁻³) in the inversion with the isostatic residual (Fig. 8b). However, these density anomalies are not continuous over the entire depth range of the slabs and are limited to depths between 400 and 700 km. Similar results are obtained using the crustal residual as input, but the density anomalies are more sporadic and concentrated at a single depth of 400 km with stronger density variations of up to +30 kg m⁻³ (Fig. 8c).

Based on these results, we think that our density inversion gives acceptable results overall for most of the mantle. The binary classification scheme we use to extract regions from the vote maps leads to a very rough representation in the upper 300 km, which likely contributes to the poor data fit obtained using the crustal residual. In addition, we assume a constant density contrast over the vertical resolution of 100 km. Thus, the volume of relatively thin structures like slabs might be smeared out over a too large depth range.

4.2 Velocity–density relation

We compare the recovered density values with velocity values from the SMEAN2 tomography, an update of SMEAN (Becker and Boschi, 2002). SMEAN2 is itself an average of three tomography models and agrees reasonably well with the patterns of the vote maps. To correlate velocity and density, every potential density anomaly is assigned the average value of SMEAN2 over the area of the potential density anomaly.

Next, we calculated the distribution of the ρ∕v_s ratio at each depth. For this distribution, we used the variability in each density value that we calculated during the cross-validation procedure. In addition, we weighted the distribution by the area of each potential density anomaly. Two histograms were estimated for the fast and slow anomalies.

In the case of the isostatic residual, the median conversion factor changes significantly with depth (Fig. 9a). In the upper 300 km it is on average close to zero, starts to increase at depths greater than 300 km and reaches a maximum of ∼0.2 at 1000 km depth. Remarkably, the conversion factor then decreases rapidly with depth, becoming negative at around 1200 km depth. The distributions at the different depths are highly variable as well. At most depths, the scaling factor shows a multimodal distribution, so there is no single consistent velocity–density relation.

Figure 9 Distribution of density-anomaly-to-velocity-anomaly ratios as a function of depth. The blue line is the distribution for fast velocity anomalies, and the red line is for slow velocity anomalies. The yellow line gives the median of the ratio distribution. Note that the median sometimes appears to lie outside the distribution, which is due to heavy tails that are not visible on this scale. Furthermore, the histograms for slow and fast anomalies are scaled to have the same maximum height for easier visual comparison. (a) Scaling for results with isostatic residual. (b) Scaling for results with crustal residuals.

The velocity–density scaling relations we calculated could have the following causes. The small scaling values in the upper 300 km is due to the fact that we use the isostatic residual. The increased scaling factors in the transition zone depths could be related to undulations in the phase transitions depths, because the higher-pressure phases have both higher density and velocity. Thus, any undulation of the phase transition would strongly affect both density and velocity. The negative velocity-anomaly-to-density-anomaly ratios we calculated below 1200 km depth zone can only be explained by compositional variations. In Appendix A we calculate sensitivities of density and velocity to temperature, iron content and MORB (mid-oceanic ridge basalt) fraction changes. We find that the sensitivities to iron content and MORB fraction are fairly constant below the transition zone, whereas the temperature sensitivity slowly decreases with depth. Thus, the velocity and density values below the transition zone could be more indicative of compositional variations than above the transition zone.

The joint geodynamic–tomographic model Gypsum (Simmons et al., 2009) also recovered a depth-variable velocity anomaly to density anomaly scaling of similar magnitude as our results. In addition, they also obtained an increase towards a maximum at around 1000 km depth. However, while their results show a decrease below 1000 km, they did not obtain systematically negative scaling factors, because they regularized towards a thermal scaling factor, which is always positive.

However, the negative conversion factors might also be affected by unmodelled dynamic effects. Depending on the viscosity model, the geoid kernel can become negative at the longest wavelengths in the lower mantle (Steinberger, 2016). This is mainly due to the effect of dynamic topography at the surface and the CMB, which has the opposite sign of the density anomaly causing it. However, we do not think that this effect alone is sufficient to explain the sign change in the lower mantle. We use the gravity field, not the geoid, and the former is more sensitive to the dynamic surface topography, which is predicted to be small, based on our density model (see next section). Still, dynamic compensation of masses in the lower mantle on the CMB could bias our results in the direction of negative conversion factors (see also Sect. 4.4).

In the case of the crustal residual, the increased overall density variation leads to higher velocity anomaly to density anomaly scaling values (Fig. 9b). However, the distribution of the scaling values are so broad that there is not a clear change in the scaling with depth.

4.3 Surface topography

We calculated the contribution of the mantle to surface topography. This allows one to assess the quality of our results, because the predicted mantle topography contribution should be zero for the isostatic residual and equal to the residual topography in the crustal residual case.

In the isostatic case, the upper mantle contributes less than 250 m to surface topography (Fig. 10a). Contrary to expectation, there is a slight positive topographic contribution in some regions with thick lithosphere (North America, eastern Australia, East Antarctica) due to negative density anomalies in the upper mantle. The likely explanation is that we only use the Moho depth to compensate topography, but in reality there is an additional contribution from the lithospheric mantle. Thus, our isostatic Moho depths is too shallow, making the isostatic correction too positive, which is then corrected by the density inversion.

Figure 10 Isostatic contributions to boundary deformations from inverted density model. To facilitate comparison with the dynamic topography in Fig. 11 the results have been filtered to spherical harmonic degree 30. (a) Contribution of upper 300 km to surface topography for the inversion results based on the isostatic residual. (b) The same for the results based on the crustal residual. (c) Contribution of the lower 800 km of the mantle to topography at the core–mantle boundary based on the isostatic residual inversion. Positive values indicate a depression of the CMB. (d) The same for the crustal residual inversion.

These results are also confirmed by the dynamic topography calculation. With the two-layer viscosity model that we chose, the topography kernels are essentially zero except for the upper 300 km of the mantle. Hence, the predicted dynamic topography is very similar to the isostatic topography. Even at degree 2, where one would expect the strongest signal from dynamic topography, according to the kernels, there is hardly any topographic contribution, because the high-density LLVPs (the dominating degree 2 structure) are overlaid by negative density anomalies. The existence of a low-density layer on top of the LLVPs was also observed in a mantle convection simulation (Liu and Zhong, 2015).

The inversion results with the crustal residual reproduce the main feature of the residual topography of our crustal model to first order. The isostatic topography contribution of the upper 300 km is considerable (up to 1.5 km). In the North American and Siberian craton, the isostatic topography is mainly negative due to the strong influence of thick, dense lithosphere, but in other areas (e.g. South Africa) the topographic contribution is negative. In the oceans, the expected signal from ocean floor spreading can only be seen in the Pacific. This is because the vote maps only reflect the oceanic cooling trends in a very crude way, such that only the broad ridges in the Pacific can be captured using this technique.

Still, there is qualitative agreement between the residual topography of the crustal model we used and the topography contribution from the upper mantle. The remaining misfit is due to an imperfect crustal model and the poor representation of the lithospheric mantle in the vote maps. In a recent model (Afonso et al., 2019), it was found that topography and the gravity field can be explained by the crust and upper mantle (for wavelengths shorter than spherical harmonic degree 15), while still staying within the uncertainty in the global seismological model Litho1.0 (Pasyanos et al., 2014). This suggests that when using a crustal model to account for the upper Earth, one should jointly consider the uncertainties in the crustal model, possibly even co-inverting a crustal model.

4.4 Impact of CMB topography

We did not explicitly include the gravity effect of CMB deformation in our inversion. To estimate the impact of CMB topography, we use the results based on the isostatic residual gravity. If the lower 800 km of the mantle is isostatically balanced on the CMB and assuming a density contrast of 4500 kg m⁻³ at the CMB, the inverted density structure would lead to CMB undulations of ±1200 m. The CMB topography is dominated by the depressions underneath Africa and the Pacific that are caused by the increased density of the LLVPs. At satellite altitude, these undulations would cause a gravity response of ±20 mGal. Due to the great depth of the source, the gravity response is anti-correlated with the proper gravity effect of the density anomalies in the lower 800 km of the mantle. Thus, the gravity effect of isostatic CMB deformations would reduce the gravity effect associated with lower-mantle density anomalies from ±30 to ±15 mGal but does not affect the spatial pattern. Hence, the actual density anomalies in the lower mantle would have to be twice as large to produce the same gravity effect as without CMB deformation. Therefore, we expect that the inversion underestimates the true density variations in the lowermost mantle.

These results also hold if dynamic topography is calculated using a two-layer viscosity model at least at the wavelengths where the gravity signal of the CMB would be visible (compare Figs. 10d and 11d). However, several studies have found that the convection dynamics in the lowest mantle might be more complicated. Viscosity probably varies laterally, because of strong temperature or composition differences (Lassak et al., 2010; Deschamps et al., 2018). In addition, post-perovskite might have much lower viscosity than perovskite (Deschamps and Li, 2019). In some scenarios the CMB deformation decorrelates from the density distribution in the lowest mantle. These findings imply that there is an additional source of uncertainty for the LLVP densities that we estimated.

Figure 11 Dynamic contributions to boundary deformations from inverted density model and a two-layer viscosity model. The results are filtered to spherical harmonic degree 30. (a) Contribution of upper 300 km to surface topography for the inversion results based on the isostatic residual. (b) The same for the results based on the crustal residual. Note that some structures appear which were not seen in Fig. 10b. This is caused by strong density anomalies below 300 km depth, which were not considered in the isostatic calculation shown in Fig. 10b. (c) Dynamic core–mantle boundary topography based on the isostatic residual inversion. Positive values indicate a depression of the CMB. (d) The same for the crustal residual inversion.

Comparing our results with published CMB topography estimates based on seismological methods (Morelli and Dziewonski, 1987; Doornbos and Hilton, 1989; Sze and van der Hilst, 2003; Tanaka, 2010), we find that none of the models contain the first-order pattern of CMB depressions beneath Africa and the central Pacific, seen in our maps (Fig. 10c, d). In terms of the magnitude, the CMB topography based on the isostatic residual are much smaller than any seismic determination except Sze and van der Hilst (2003), whereas the crustal residual CMB topography has similar range as Doornbos and Hilton (1989) and Morelli and Dziewonski (1987).

4.5 Isostatic residual vs. crustal correction

The results obtained with isostatic and crustal residual disagree substantially. The correlation coefficient calculated for individual depth slices is typically less than 0.5 and even negative for some depths. The highest correlations are found in the depth ranges of 2500–2800 and 1500–2000 km and are mainly due to the influence of the LLVPs. In addition, the magnitude of the density anomalies is on average 4 times larger based on the crustal residual, which agrees with the relative magnitude of the input gravity fields.

Based on the fit to the data, the isostatic residual is preferred, because it achieves both higher absolute and relative data fit. However, the isostatic Moho depths disagree with the seismological determinations. Furthermore, the crustal residual lacks a correction for the oceanic cooling trend, which clearly contributes to the crustal residual. Thus, the better data fit preference for the isostatic residual is not as straightforward.

In any case, the isostatic and crustal residuals are extreme examples of how the crust and upper mantle can be accounted for in a gravity inversion. Our results clearly demonstrate that these different approaches lead to an enormous spread in terms of the recovered densities in the mantle. At the same time, signals from the deep Earth in terms of gravity or topography might affect modelling of the upper Earth. Commonly, high-pass filtering is used to remove the signal of the deep Earth (Bowin, 1991), but this is clearly insufficient, due to the spectral overlap (Root et al., 2015). Thus, a coupled approach that simultaneously considers the entire crust and mantle is required in order to properly model the gravity field and Earth's topography.

4.6 Implications for LLVP temperature and composition

Broadly speaking, the LLVPs could originate from any combination of increased temperature and compositional variation. After their first discovery, the LLVP were initially considered as a purely thermal feature: two “superplumes” that rise from the core–mantle boundary (Dziewonski, 1984). However, even if the LLVPs are purely due to temperature increase, a more likely explanation is that they are swarms of smaller plumes that are smeared due to the limited resolution of seismic tomographies (Schubert et al., 2004; Schuberth et al., 2009), with an additional influence from the stability field of post-perovskite (Koelemeijer et al., 2018).

In contrast to this isochemical view, some authors have proposed that the LLVPs are chemically distinct from the normal (pyrolitic) mantle. The chemical distinctiveness of the LLVP can either be accumulated over time or be a primitive reservoir that separated early in Earth's history (Deschamps et al., 2012). A likely process for accumulation is the separation of oceanic crust from subducted lithosphere (Mulyukova et al., 2015). However, these two hypotheses are not mutually exclusive, and a small amount of subducted oceanic lithosphere might be incorporated in pre-existing structure, as proposed in the basal melange (BAM) hypothesis (Tackley, 2012). If the LLVPs have a compositional component, they have to be intrinsically denser than surrounding pyrolitic mantle to have remained near the core–mantle boundary, despite very high temperatures. If the LLVPs are indeed old features, this would also explain the proposed spatial correlation between plume generation and the edges of the LLVPs over the last 200 million years (Burke et al., 2008).

Our inversion results indicate a positive density anomaly for the LLVPs both using isostatic and crustal residual. This would rule out a purely thermal origin of the LLVPs, since such a scenario would lead to negative density anomalies.

In order to test different scenarios using our inversion results, we make use of the petrological database of Stixrude and Lithgow-Bertelloni (2011). We construct a simple adiabatic model (see Appendix A) based on a pyrolitic composition for the mantle (Stixrude and Lithgow-Bertelloni, 2012) that serves as reference model. The adiabatic model is entirely self-consistent, and the only free parameter is the temperature at the top of the model, which we adjusted in order to fit the surface wave dispersion curves from PREM (Dziewonski and Anderson, 1981). Unlike previous thermochemical interpretations for the lower mantle (e.g. Deschamps et al., 2012), this corrects for the bias between the petrophysical database and PREM. Next, we applied first-order perturbations to the model and obtained sensitivities of shear wave velocity and density to changes in temperature, iron content and fraction of mid-ocean ridge basalt (MORB).

We assume a negative S-wave velocity deviation of 2 % for the LLVPs, since most mantle tomographies display roughly this amount of slowness. Our inversion results would place the density anomaly of the LLVP between 0.1 % (isostatic residual) and 0.3 % (crustal residual), but due to possible isostatic compensation at the CMB, (Sect. 4.4), the density anomalies could be twice as large (up to 0.6 %).

A temperature increase of 670 K leads to the required velocity reduction; however it would also entail a density change of −1 %, which would be incompatible with our findings. Likewise, a 2.6 % increase in iron content (without temperature change) fits the velocity reduction but leads to a density increase of 1.6 %. To fit both velocity and our lowest density estimate, a temperature increase of 380 K and an iron increase of 1.1 % is required, while our highest density estimate would require changes of 260 K and 1.6 % FeO respectively.

Adding a MORB fraction leads to higher required temperatures, because MORB is slightly faster than pyrolitic mantle at the depths of the LLVPs, according to our petrophysical calculations. Our lowest density estimate (+0.1 %) would require a temperature change of +870 K and a MORB fraction of 40 %, and the highest density estimate (+0.6 %) requires +960 K and 58 % MORB. Deschamps et al. (2012) also found high excess temperatures (+1500 K) were necessary to explain the shear wave velocity and bulk sound speed if the LLVPs are MORB-enriched. However, our interpretation is dependent on the petrological database of Stixrude and Lithgow-Bertelloni (2011), which is based on sparse measurements at extremely high pressures. Recent measurements of calcium–perovskite (Ca–Pv) (Gréaux et al., 2019; Thomson et al., 2019) imply that seismic velocities in MORB-enriched zones would be lower than thought previously. Thus, less extreme temperatures would be required to reconcile increased density and reduced velocity. Furthermore, the chemical composition of MORB is highly variable, especially its iron content (Deschamps et al., 2012), which could also strongly affect the impact of mixing a MORB fraction into a pyrolitic mantle.

Based on our results it is difficult to express preference for MORB or iron enrichment. Since MORB is introduced to the mantle by subduction, plate reconstructions place some constraints on the amount of MORB produced. Stixrude and Lithgow-Bertelloni (2012) estimate that the total amount of basalt input into the mantle corresponds to about 10 % of the volume of the mantle. As an absolute upper limit, assuming the LLVPs have a height of 800 km and cover 25 % of the core surface area, the volume would be around 4 % of the total mantle. Thus, to reach a MORB fraction of 60 % in the LLVPs, about a fourth of the total generated subducted basalt would need to accumulate.

5 Conclusions

In this paper we have presented how the gravity field and tomographic vote maps can be combined to estimate the density distribution inside the mantle. We have shown that the method is able to reasonably recover expected features in the mantle without requiring any information about the viscosity structure. Still, our recovered density structure leads to qualitative agreement between isostatic and residual topography. Furthermore, our results indicate that the LLVPs are slightly denser and hence chemically distinct.

In our first-order analysis, the nonlinear dependencies of velocity and density, the impact of melt, heterogeneity inside the LLVPs, and the possible presence of post-perovskite is neglected. However, our results can be reconciled with our present knowledge about rock properties at these extreme conditions.

One important difference compared to previous methods is that in our method density is free to vary independent of seismic velocity. While this gives the necessary freedom to the density inversion, it ignores the strong evidence for the important role of temperature. In the future, more precise petrophysical data could help to put constraints on the relative importance of temperature and composition. This would also help to reconcile our results with previous viscosity inversions.

Furthermore, our results show that the lithospheric mantle is critical to resolve disagreements between bottom-up and top-down methods. In fact, depending on how the lithosphere is treated, the inverted density anomalies can change by a factor of 4. Thus, a combined approach is necessary, where the uncertainties in seismic determinations of the lithosphere are considered in conjunction with signals resulting from deeper density anomalies.

Our inversion breaks down the mantle volume into discrete volumes with a constant density anomaly. This reduces the number of unknowns compared to a continuous inversion and means that less regularization is required. However, using vote maps to extract features from a collection of seismic tomographies is the most basic way to do this, and a more refined method (e.g. Fadel et al., 2015) could lead to better constrained volumes. In particular it might benefit an inversion if uncertainties could be placed on the size and position of the volumetric features to allow them to change during inversion.

Instead of relying on seismic tomographies or derived quantities like vote maps, it might be beneficial to stay closer to seismic data. The canonical choice would be seismic travel times and normal modes (as in Dziewonski and Anderson, 1981), but there are more uncommon seismic products that might provide additional constraints. For instance, core-diffracted waves are very sensitive to velocity near the CMB (Hosseini and Sigloch, 2015).

Appendix A: Adiabatic model

We derived a simple 1-D adiabatic model of Earth's mantle, based on a petrological database (Stixrude and Lithgow-Bertelloni, 2005, 2011). We assume a pyrolitic composition (45.1 % SiO₂, 38.1 % MgO, 8.0 % FeO, 3.3 % Al₂O₃, 3.1 % CaO) throughout the entire mantle.

The main relations – density ρ(P,T), specific heat capacity c_p(P,T) and heat expansivity α(P,T) – are functions of temperature and pressure and are calculated using the Stixrude database with PerpleX (Connolly, 2009). Pressure is purely hydrostatic as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial P}{\partial r}}=-\mathit{\rho }g.\end{array}$$

The temperature is purely adiabatic, such that

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\partial T}{\partial r}}=-{\displaystyle \frac{T\mathit{\alpha }g}{c_p}}.\end{array}$$

The gravity acceleration at a specific depth is determined by the internal mass at that depth

$$\begin{array}{}\text{(A3)}& g\left(r\right)={\displaystyle \frac{G{M}_{\mathrm{int}}\left(r\right)}{r^{\mathrm{2}}}},\end{array}$$

and the internal mass is decreasing with depth according to

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\partial M_{\mathrm{int}}}{\partial r}}=+\mathrm{4}\mathit{\pi }{r}^{\mathrm{2}}\mathit{\rho }.\end{array}$$

This set of equations can be solved by finite differences if the pressure, temperature and gravity strength are specified at the top of the model domain. First, the ρ, c_p and α belonging to these P−T conditions are calculated using PerpleX (Connolly, 2009). Then, these values are used to update pressure, temperature and internal mass as follows:

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}P(r-\mathrm{\Delta }r)& =P\left(r\right)+\mathit{\rho }\left(r\right)g\left(r\right)\mathrm{\Delta }r\\ T(r-\mathrm{\Delta }r)& =T\left(r\right)+{\displaystyle \frac{T\left(r\right)\mathit{\alpha }\left(r\right)g\left(r\right)}{c_p\left(r\right)}}\mathrm{\Delta }r\\ M_{\mathrm{int}}(r-\mathrm{\Delta }r)& =M_{\mathrm{int}}\left(r\right)-\mathrm{4}\mathit{\pi }{r}^{\mathrm{2}}\mathit{\rho }\left(r\right)\mathrm{\Delta }r.\end{array}\end{array}$$

This process is repeated iteratively from the top of the model domain down to the core–mantle boundary, with a step size Δr of 1 km.

We begin this integration at a depth of 80 km, with a pressure of 2.5 GPa and g=9.81 m s⁻². In order to determine the temperature at the top, we relied on PREM (Dziewonski and Anderson, 1981). At first we attempted to directly minimize the difference between the velocities in PREM and the modelled velocity from PerpleX. However, this approach was unsuccessful, because the depths of the main discontinuities in PREM (410, 660 km) are incompatible with the Stixrude database.

Instead, we relied on Rayleigh wave dispersion curves to choose the temperature at the top of the model. Using Mineos (Masters et al., 2019) we calculated the phase velocity of the fundamental Rayleigh wave between frequencies of 0 and 50 mHz for PREM and for the adiabatic model, for temperatures at the top between 1300 and 1800 K in steps of 50 K. We found that a top temperature of 1550 K gives the best fit to the PREM dispersion curves (Fig. A1a), if equal weight is given to phase velocities at all periods. In particular, we find that the T=1550 K model fits the dispersion curve data with an average relative accuracy of 0.3 %, which is similar to the fit with which PREM fits the data used in its construction. Thus, our model is equivalent to PREM with respect to the dispersion curve data. During these calculations we also found that the presence of a post-perovskite phase completely prevents fitting the long-period dispersion curves, so we excluded post-perovskite from the PerpleX calculations.

The resulting temperature curve for the preferred model is nearly linear but shows a distinct kink below the 660 km discontinuity, due to the different properties of perovskite and a slower temperature increase at greater depths due to the decrease in thermal expansivity α (Fig. A1b). The temperature at the bottom of the mantle is roughly 2350 K. The pressure curve is also nearly linear but is slightly bent due to the increase in density with depth.

We then applied first-order perturbations in terms of temperature, iron content and mid-oceanic ridge basalt (MORB) fraction to the adiabatic model. The vertical resolution of our density models is 100 km, so we applied the perturbation over the same depth range. To determine the sensitivity to temperature variations, we simply used our existing lookup table from PerpleX, while for compositional variations we determined new phase equilibria and corresponding rock properties with FeO content increased by 1 %. For the MORB we proceeded somewhat differently, because the MORB is likely not in phase equilibrium with the surrounding mantle rocks, due to the long timescale of chemical diffusion (Stixrude and Lithgow-Bertelloni, 2012). Thus, we determined the phase equilibrium of a pure MORB and then calculated velocities and densities as volume averages of the MORB fraction and the surrounding mantle. These results can be used together to determine how sensitive v_s and density are to changes in temperature, iron content or MORB fraction (Fig. A2).

Figure A1 Results of the adiabatic model. (a) Rayleigh temperatures at the top of the adiabatic model (80 km). (b) Temperature and pressure profiles for the preferred model with T_top=1550 K.

Figure A2 First-order perturbation curves for deviations from a adiabatic, hydrostatic model with pyrolitic composition. The red and blue background shading indicates negative and positive changes of the respective quantity. (a) Relative change for a temperature change of 100 K. (b) Relative change for including 1 % more FeO. (c) Relative change for a mixture of 90 % pyrolite and 10 % MORB.