Subducted oceanic crust as the origin of seismically slow lower-mantle structures

Mantle tomography reveals the existence of two large low-shear-velocity provinces (LLSVPs) at the base of the mantle. We examine here the hypothesis that they are piles of oceanic crust that have steadily accumulated and warmed over billions of years. We use existing global geodynamic models in which dense oceanic crust forms at divergent plate boundaries and subducts at convergent ones. The model suite covers the predicted density range for oceanic crust over lower mantle conditions. To meaningfully compare our geodynamic models to tomographic structures, we convert them into models of seismic wavespeed and explicitly account for the limited resolving power of tomography. Our results demonstrate that long-term recycling of dense oceanic crust naturally leads to the formation of thermochemical piles with seismic characteristics similar to the LLSVPs. The extent to which oceanic crust contributes to the LLSVPs depends upon its density in the lower mantle for which accurate data is lacking. We find that the LLSVPs are not composed solely of oceanic crust. Rather, they are basalt rich at their base (bottom 100–200 km) and grade into peridotite toward their sides and top with the strength of their seismic signature arising from the dominant role of temperature. We conclude that recycling of oceanic crust, if sufficiently dense, has a strong influence on the thermal and chemical evolution of Earth’s mantle.


Introduction
Tomographic models reveal the presence of two large anomalous structures in the lowermost mantle. These large low-shear velocity provinces (LLSVPs) are surrounded by high shear-wave velocity regions that are associated with past subduction ( Fig. 1; Tanaka et al. 2009;Ritsema et al. 2011;and and Cottaar and Lekic 2016). The origin of the LLSVPs and their connection to mantle dynamics and the long-term chemical evolution of the planet remain unclear (see McNamara (2019) for a recent overview). Suggested origins include large-scale thermal upwellings (e.g., Schubert et al. 2004;Davies et al.2015;and Koelemeijer et al. 2017), dense piles of primordial and compositionally distinct mantle (e.g., Garnero and McNamara 2008;Burke et al. 2008), and accumulations *Correspondence: tjones@carnegiescience.edu 1 Department of Terrestrial Magnetism, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC, 20015 USA Full list of author information is available at the end of the article of recycled oceanic crust (e.g., Christensen and Hofmann 1994;Nakagawa et al. 2009Nakagawa et al. , 2010. Markedly different interpretations of LLSVPs persist because there is no unique cause of seismic velocity anomalies in the mantle, but there are many plausible ones. For example, the low-shear velocity anomaly (δ ln V s ≈ −2%) associated with LLSVPs may be equally well explained by an increase in temperature (Schuberth et al. 2009), a compositional change (McNamara and Zhong 2005), the presence of melt (Lee et al. 2010), or a combination thereof (Li and McNamara 2013;Ballmer et al. 2016). The elevated ratio between shear and compressional wavespeed anomalies within LLSVPs (Robertson and Woodhouse 1995;Masters et al. 2000;Romanowicz 2001;Moulik and Ekström 2016) is often regarded as evidence that they are thermochemical in origin, but it can also indicate the presence of post-perovskite without requiring a contribution from large-scale chemical variations (Koelemeijer et al. 2018). Measurements of solid-Earth tides (Lau et al. 2017) and (2020) (Ritsema, 2011) showing variations in shear-wave velocity, δ ln V s , with respect to the radial average. Both positive and negative velocity anomalies are highest toward the CMB. The color intensity is proportional to the velocity amplitude up to a maximum value of 2% normal mode splitting functions (e.g., Ishii and Tromp 1999;Trampert et al. 2004;Moulik and Ekström 2016;and Koelemeijer et al. 2016) suggest that the LLSVPs have a higher than average density across the bottom two thirds of their depth range. However, constraining the density from low-frequency data depends on the chosen set of normal modes, on how CMB topography is modeled, and on the theoretical simplifications made to model the data (e.g., Koelemeijer et al. 2017;Akbarashrafi et al. 2018).
In this study, we will test the hypothesis that LLSVPs are accumulations of dense recycled oceanic crust. We are motivated in particular by geodynamical studies that show that the long-term recycling of oceanic crust can significantly contribute to the formation of warm but dense pools of compositionally distinct material at the base of the mantle (Fig. 2; Christensen and Hofmann 1994;Brandenburg and van Keken 2007;Brandenburg et al. 2008;and Nakagawa and Tackley 2008). We will focus on how well thermochemical model simulations predict the seismological characteristics of LLSVPs by projecting the temperature and composition to seismic velocity. We will then test how well such velocity structures are resolved in the shear-wave velocity model, S40RTS (Ritsema et al. 2011), and the combined shear-and compressional-wave velocity model, SP12RTS (Koelemeijer et al. 2016).
As oceanic lithosphere subducts, three distinct components are recycled back into the mantle: a thin basaltic crust and a thicker base of harzburgite on top of the more primitive peridotite from which the crust was extracted. Due to convective stirring of these components, the mantle composition is an evolving mixture of the three endmember compositions. To simplify the description below, we will refer to end-member basaltic, harzburgitic, and peridotitic compositions as "basalt", "harzburgite", and "peridotite", respectively, and, where relevant, will use subscripts B, H, and P to identify them.
Mineral physics experiments indicate that the basaltic component is denser than ambient mantle through most of Earth's mantle with a relative density contrast δ ln ρ of up to 2-7%, except near the surface and in the uppermost lower mantle (Ringwood and Irifune 1988;(2020) 7:17 Page 3 of 16

Fig. 2
Plot of temperature and composition fields for Brandenburg et al. (2008) models for δ ln ρ of 0% (a, b) and 3.8% (c, d), where we use the non-linear radial projection of the models onto Earth coordinates. In both cases, recently subducted oceanic crust meets the CMB and the low viscosity upper mantle has been homogenized by convective stirring relative to the lower mantle. When δ ln ρ = 0%, oceanic crust is somewhat evenly distributed throughout the lower mantle. In contrast, when δ ln ρ = 3.8%, oceanic crust accumulates and the CMB forms sharp piles Ringwood 1990;Hirose et al. 1999;Aoki and Takahashi 2004;Hirose et al. 2005;Ricolleau et al. 2010;Tsuchiya 2011). Under these conditions, numerical experiments predict that the oceanic crust will accumulate at the core-mantle-boundary (CMB) (Christensen and Hofmann 1994;Davies 2002;Xie and Tackley 2004;Brandenburg and van Keken 2007). The presence of recycled oceanic crust as a major component of mantle plumes (Hofmann and White 1982) is empirical evidence that at least some portion of subducted material reaches the CMB.
The models presented in this study follow the approach of Christensen and Hofmann (1994) (herein CH94) who assumed that the basalt remains denser throughout the mantle and retains a constant density contrast with the ambient mantle. CH94 found that at sufficiently large density excess, the subducted oceanic crust accumulates at the CMB where it is swept into piles and mixes with the ambient mantle. They concluded that this process could not only account for important features of Earth's chemical evolution but may also explain seismic heterogeneity in the D" layer.
Due to computational limitations at the time, CH94 performed dynamical computations at a convective vigor that is significantly below that of present day Earth and used a scaling relationship to extrapolate their findings to models with present-day plate velocities. Brandenburg and van Keken (2007) explicitly showed that this extrapolation is an accurate approximation. Brandenburg et al. (2008) (herein BB08) further extended the study of CH94 by exploring the effects of crustal recycling in models with force-balanced plates and with Earth-like convective vigor. Their results confirmed the findings of CH94. They also showed that the formation and recycling of oceanic (2020)  crust could contribute to the HIMU-DMM-EMI diversity of Ocean Island Basalt (OIB) and Mid-Oceanic Ridge Basalt (MORB) in multiple isotope systems. The BB08 study confirmed that recycled dense oceanic crust forms thermochemical piles. These piles naturally form sharp edges, as observed in seismological studies (Ni et al. 2002;Ritsema et al. 1998), and they are broadly consistent with the size and shape of LLSVPs (Fig. 2 frames c and d). However, any assessment of similarity between thermochemical models and global tomographic models must take into account the limited spatial resolution of seismic tomography. As earlier studies have shown, shearwave velocity perturbations can be reduced, smeared, and spatially shifted due to the effects of uneven data coverage and damping (Ritsema et al. 2007;Schuberth et al. 2009). In this paper, we will map the temperature and compositional fields of the BB08 models into seismic velocity and will then account for the parameterization and regularization parameters of S40RTS (Ritsema et al. 2011) and SP12RTS (Koelemeijer et al. 2016). Our analysis of these seismologically filtered geodynamical models will include a comparison of their spectral content and radially averaged profiles to that of S40RTS and SP12RTS. We conclude that the recycling of oceanic crust, by itself or in addition to other processes, can explain the imaged seismic heterogeneity that is expressed by the LLSVPs as long as the oceanic crust has a moderate to high excess density with respect to the ambient mantle.

Geodynamic models
The BB08 models are finite element solutions of convective flow in an incompressible fluid at an infinite Prandtl number driven by thermal and chemical buoyancy. The models are based on a 2D cylindrical geometry with a reduced core radius to mimic the heat loss characteristics of Earth (van Keken 2001).
The models use a force-balanced plate approach (Gable et al. 1991) that yields average surface velocities and surface heatflow consistent with values observed today. This sets up energetically consistent convection with prescribed zones of divergence and convergence where oceanic crust forms and subduct, respectively. The resolution of the models varies from 50 km in the middlelower mantle to 13 km at both the surface and CMB, such that the thermal boundary layers, of characteristic thickness ∼100 km are adequately resolved. At divergent zones, a melting parameterization generates basaltic crust with a harzburgite residue from the peridotitic mantle to form oceanic lithosphere. When this lithosphere subducts at convergent zones, it sinks to the base of the mantle under negative thermal buoyancy before mixing with the surrounding mantle. We represent the peridotitic composition by tracers that are spaced, on average, at 5 km and that concentrate by a factor of 8 when the oceanic crust is formed such that a ∼10-km-thick basaltic crust is adequately modeled.
In non-dimensional form, the governing equations are the conservation of mass: ( 1 ) the conservation of momentum: and the conservation of heat: where u is the velocity vector, P the dynamic pressure, t time, T the temperature, andĝ the unit vector in the direction of gravity. C is the chemical composition, η the non-dimensional dynamic viscosity, and Q is the volumetric internal heating.˙ is the strain-rate tensor˙ = ∇u + ∇u T . Ra is the thermal Rayleigh number: where T is the assumed temperature contrast across the mantle and h is the thickness of the mantle. ρ 0 , κ 0 , α 0 , and η 0 are the reference values for density, thermal diffusivity, thermal expansivity, and dynamic viscosity, respectively. Rc is the compositional Rayleigh number: where ρ 0 is related to the density contrast between the oceanic crust and the depleted mantle (see below). For reference values, see Table 1. The models assume constant expansivity and diffusivity, and a temperature-and depth-dependent viscosity of the non-dimensional form where z is depth, b = ln(1000), and η z is 1000 in the lithosphere, 1 in the upper mantle, and 30 in the lower mantle. BB08 unfortunately introduced a number of notation errors that include typos in the equations (e.g., their equation (10) should read f = N/(N + 7M), not f = 1/(N + 7M)) and inverted the meaning of the term "extraction coefficient" where "retention coefficient" was intended. It also obscured the definition of the eclogite excess density (EED) as it was not fully defined. The EED used in their paper was defined as EED=Rc/Ra= ρ/(ρ 0 α 0 T) and ranged from 0 to 10%. The composition function C (Eq. 2) is 1 for peridotite, 0 for harzburgite (which has density ρ 0 ), and 8 for basalt. This means that the dimensional density contrast of the oceanic crust to depleted mantle ρ B−H = ρ B − ρ H can be found by multiplying EED by 8ρ 0 α 0 T. The more commonly used expression for relative excess density ∂ ln ρ = (ρ B − ρ P )/ρ P is with respect to ambient mantle found by multiplying EED by 7ρ 0 α 0 T. Using ρ 0 =4500 kg/m 3 , α 0 =3 × 10 −5 /K, and T=3000 K, we find ρ B−P =EED×2840 kg/m 3 . For EED=10% (the maximum value considered), we find a dimensional density of the oceanic crust to be ρ B =4824 kg/m 3 which translates to a relative density increase of δ ln ρ = 6.2%.
We select seven simulations that are similar to those in BB08 except that the models have been evaluated on a higher resolution finite element mesh (with 70 K vs. 46 K nodal points) and a larger number of tracers (3.7 M vs. 2.5 M). Between each model, we vary only the the density difference between basalt and harzburgite. This has some control over the dynamics, which leads to differences in the spatial distribution of chemical heterogeneity throughout the mantle.
Each simulation starts from snapshot of a purely thermal convection simulation at Ra=5 × 10 6 that was evolved for a sufficiently long time to let any of the transients associated with the initial conditions disappear. We then model oceanic crust formation and recycling for a period equivalent to 4.5 Byr and for various values for Rc (between 0 and 5 × 10 5 ). Piles of recycled crust build up over time for Rc > 0. Since we keep Rc/Ra fixed for each simulation, the accumulation of oceanic crust for the higher δ ln ρ models tends to make the convection slightly more sluggish towards the end of the model evolution (consistent with BB08).

Model selection and seismic wavespeed conversion
We first select a time slice from the past 2.5 Gyr of each simulation in which a broad thermochemical pile is present (for Rc > 0). Piles tend not to approach the size of LLSVPs earlier than about 2 Ga and are more robust for higher δ ln ρ. We also use a time-slice for Rc = 0 for comparison.
Two steps are required to create an Earth-like geometry from the 2D cylindrical models. We use a non-linear radial projection of the models onto Earth coordinates (i.e., the non-dimensional core radius 0.492 maps onto a CMB depth of 2885 km and the non-dimensional surface radius 1.492 corresponds to depth 0 km; van Keken (2001)). We then project a 180°section of the geodynamic model into an axisymmetric spherical shell representation that is representative of a fully 3D spherical model with the important constraint that the thermochemical pile is centered about the rotation axis. This effectively imposes the dome-like geometry characteristic of LLSVPs. Structures away from the rotation axis thus become ring-like. The 3D spherical model is then converted from temperature and composition to seismic wavespeed using the following approach. First, we convert potential temperature to absolute temperature by adding an adiabatic gradient for a pyrolite composition. This step is required because the geodynamic simulations assume the incompressible Boussinesq approximation.
Next, we map the pressure P, temperature T, and basalt fraction f to seismic shear-wave velocity V s and compressional-wave velocity V p using the method of Xu et al. (2008), where it is assumed that the mantle can be represented as a mechanical mixture of basalt and harzburgite. The bulk composition of the two end members are defined in terms of the proportions of the oxides SiO 2 , Na 2 O, CaO, FeO, MgO, and Al 2 O 3 .
We use the composition of basalt from Workman and Hart (2005) and the composition of harzburgite from Baker and Beckett (1999). The stable mineral assemblage and associated anharmonic seismic wavespeeds are computed for both basalt and harzburgite using Per-ple_X (Connolly 2005) and the mineral elastic parameter database of Stixrude and Lithgow-Bertelloni (2011).
The anharmonic seismic wavespeeds are corrected for frequency-dependent effects of attenuation using Maguire et al. (2016) attenuation model. The seismic wave speed V (V s or V p ) of intermediate compositions is taken as the linear combination of the end member compositions basalt and harzburgite: The velocity anomaly δ ln V is computed relative to the radially averaged seismic structure.

Reparameterization and the seismic resolution operator
To meaningfully compare geodynamic models with tomographic models, we must account for the parameterization and limited resolution of the latter. This process is illustrated in Fig. 3. We begin by reparameterizing our geodynamically predicted wavespeed model onto the basis of S40RTS. Lateral variations of δ ln V s are described using spherical harmonics up to degree and order 40, and depth variations by 21 vertical splines. Following the approach of Ritsema et al. (2007), we define the resolution operator to be R = G † G, where G is the operator of the seismic forward problem and G † is its generalized inverse. We then modify our reparameterized model of seismic wavespeed (the input model m IN ) by multiplying it with R to obtain a filtered output model, m OUT , as if imaged by tomographic inversion: where R is a complete description of the variable spatial resolution of S40RTS. We compute R for the same damping parameter as that used in S40RTS.
To explore how both δ ln V s and δ ln V p are recovered after tomographic filtering, we repeat this processes using R from the joint V s and V p tomography model SP12RTS. SP12RTS uses a similar model parameterization to S40RTS but only up to spherical harmonic degree 12. Since SP12RTS incorporates information from both V s and V p , the model vector m consists of two parts: where S and P are model vectors describing dlnV s and dlnV p , respectively. Therefore, Eq. 8 becomes: Jones et al. Progress in Earth and Planetary Science (2020) where the diagonal blocks, R SS and R PP , describe how V s and V p map into V s and V p , respectively, while the offdiagonal blocks, R SP and R PS , describe how V s and V p map into each other.

Effects of reparameterization and tomographic filtering
We first look at the effects of the filtering processes itself on δ ln Vs. Figure 3 illustrates how reparameterization and R modify the input model to an output model with the resolution of S40RTS. Only the broadest scale of seismic heterogeneity of the input model is preserved following the filtering process. The highest δ ln V s anomalies are in the uppermost mantle, and the relative position of high and low seismic velocity features are well preserved in the lower mantle.
The most pronounced effects are amplitude reduction and smoothing. These effects are accentuated where data coverage is sparse. This is for example of the case beneath the African plate (Fig. 4, left column), where δ ln V s anomalies in the same filtered geodynamic input model are smaller and smoother than beneath the Pacific plate (Fig. 4, right column). Consequently, cold subducted lithosphere and the thermochemical piles are more dominant features of the lower mantle in our Pacific crosssections. Where data coverage is particularly poor, low δ ln V s anomalies disappear entirely after filtering. For example, subducted lithosphere is imaged along the CMB in the Pacific cross-section with δ ln V s ∼ 0.5% (Fig. 4, δ ln ρ=3.1%), but the same structure is absent in the African cross-section.

Effects of δ ln ρ
The spatial patterns of δ ln Vs for all filtered dynamic models and along radial segments of S40RTS are illustrated in Fig. 4. In the isochemical convection case, plume clusters (Fig. 4, bottom row) are imaged as regions of low δ ln V s in the lower mantle, but these are significantly smaller and of weaker amplitude than LLSVPs in S40RTS. Furthermore, regions of high δ ln V s are absent from the lowermost mantle.
By comparison, regions of strongly negative δ ln V s for the thermochemical convection cases are broader and of stronger amplitude. Their size does not reflect that of the basaltic crust piles. Instead, the margins of the low δ ln V s regions map onto the thermal anomaly associated with the basaltic crust piles. This implies that LLSVPs should not be interpreted as pure oceanic crust. Rather, under the parameters explored here, they are basalt rich at their base (bottom 100-200 km) and grade into peridotite toward their sides and top, which can reach mid-mantle depths.
The range of δ ln V s in all thermochemical simulations is comparable to that of S40RTS after tomographic filtering is applied. However, the low velocity structures resemble LLSVPs in their size and strength only when δ ln ρ ≥ 3.1% so that enough basaltic crust accumulates at the CMB.
The thermochemical piles are generally flanked by moderate-to-high δ ln V s anomalies, corresponding to cold lithosphere that has recently subducted. In these models, the subducted lithosphere sweeps the hot, dense material into piles. This is consistent with S40RTS, and other tomographic models, that show LLSVPs surrounded by a ring of high seismic velocity (e.g., Dziewonski et al. 1977;Grand et al. 1997;and Kennett et al. 1998).
In contrast to S40RTS, subducted lithosphere in a number of the filtered models is coherently imaged through the entire mantle. This indicates that subducted lithosphere in those cross sections has a higher temperature anomaly or are broader than subducted slabs in Earth's mantle. Our models have important limitations that may explain this difference. We force an axisymmetric spherical symmetry on the models causing the downwellings in the geodynamical models to be projected as ring-like features which likely are recovered with higher amplitude and geographical extent than subduction zones in the Earth. Additionally, on Earth, subduction is one-sided and occurs at a range of angles while in the geodynamic models, subduction is two-sided and the slabs sink almost exclusively vertical.
Given that our models are two-dimensional and have dynamic plates, we can neither draw meaning from geographical correlations nor impose past plate motions. We therefore compare the similarity of cross sections of S40RTS and filtered geodynamic models in two ways. First, we compare the heterogeneity spectra of δ ln V s by computing the power spectral density (PSD) at each depth (Figs. 5 and 6). This gives the power as a function of wavenumber k that represents the number of wavelengths per 360 • and thus are analogous to spherical harmonic degree. The heterogeneity spectrum is normalized by the maximum of the PSD at each depth. Second, we compare the root-mean-square (RMS) of δ ln V s as a function of radius (Figs. 8 and 9).
In the lower mantle, the spectral characteristics of our filtered models are similar to the two S40RTS slices (Figs. 5a and 6a) when δ ln ρ ≥ 3.1. The dominance of power at lower k in such cases is a consequence of the long-wavelength structure of the thermochemical piles.
For lower values of δ ln ρ, power is somewhat randomly spread over wavenumbers k < 20 beneath the African plate and k < 40 beneath the Pacific plate, since basaltic piles do not form in as coherent a fashion as at larger density contrasts. Spectral characteristics diverge between the filtered models and S40RTS in the upper mantle, except when δ ln ρ = 6.2%, which is the highest value used in Power spectral density plots for S40RTS and the tomographically filtered models as imaged beneath the Pacific plate. Spectral power is normalized to the maximum power at each depth and plotted as a function of wavenumber and radius. High power at low wavenumbers characterize the lower mantle in S40RTS and the filtered geodynamic models when δρ ≥3.1%. This is due to the formation of large thermochemical structures at this density contrast. Below this value, power in the lower mantle is more randomly distributed across high and low wavenumbers due to a lack of broad-scale structure. The upper mantle of all geodynamic models contains more power at higher wavenumbers compared to S40RTS this study. At the scale of our interpretation, the spectral characteristics of LLSVPs in S40RTS do not vary longitudinally. This is illustrated in plots of the spectral content through the Pacific and African LLSVPs spanning 30 • (Fig. 7). Similar correlations can be observed in the RMS δ ln V s profiles for our filtered models (Figs. 8 and 9). The S40RTS profile is given for comparison. As before, models with δ ln ρ ≥ 3.1% show good correlation to S40RTS in the lower mantle. The larger amplitudes close to the thermal boundary layers and small peaks around the mid-to-upper mantle in S40RTS are poorly matched.
Lastly, we compute the expected δ ln V s /δ ln V p ratio (S/P) of our models after filtering with the resolution operator of SP12RTS. S/P is calculated by dividing the RMS of δ ln V s by the RMS of δ ln V p at each depth. Figure 10 illustrates S/P as a function of radius for the purely thermal case, for the thermochemical cases δ ln ρ=3.1% and δ ln ρ=6.2% and for SP12RTS along the same radial segment to which filtering is applied (Fig. 10). Model slices are taken along the same longitude as Pacific mantle sections above but are centered at 15 • latitude where P-wave coverage is higher. In each case, we include the effect of post-perovskite (pPv), whose transition from bridgmanite (brg) increases the amplitude of Vs and thus S/P.
In all filtered models, the increase in S/P throughout the lower mantle is comparable to SP12RTS. In the profile shown, the magnitude of S/P of SP12RTS is higher than average due to the presence of the Pacific LLSVP. S/P is weak in the purely thermal case and strengthens in the thermochemical cases as δ ln ρ is increased. When δ ln ρ=6.2%, the magnitude of S/P far exceeds the maximum values observed in SP12RTS. Peak S/P consistently occurs at ∼ 2500-km depth in our filtered models, which is ∼ 100 km deeper than in SP12RTS. The same result was (2020) 7:17 Page 10 of 16

Fig. 6
Power spectral density plots for S40RTS and the tomographically filtered models as imaged beneath the African plate. Spectral power is normalized to the maximum power at each depth and plotted as a function of wavenumber and radius. Same basic features as for mantle beneath the Pacific except that S40RTS has more power at low wavenumbers in the upper mantle as well as in the lower mantle found by Koelemeijer et al. (2018), who proposed two scenarios to resolve the offset: either the Clapeyron slope of the brg-pPv transition shallows at high temperatures or the entire stability field of pPv occurs at shallower depths.

Discussion
Our results support the hypothesis that long-term oceanic crust formation and recycling leads naturally to the formation of relatively dense thermochemical piles in the lowermost mantle with characteristics similar to the LLSVPs in S40RTS. We find the best match if the density contrast between basaltic crust and ambient mantle is around 3% or larger. Note that the precise threshold will depend upon parameters that were not varied in this study but are nonetheless uncertain, such as the temperature dependence of viscosity and the mantle's radial viscosity structure (e.g., Lau et al. 2016). This interpretation implies a distinct thermal and chemical evolution for the seismic structure of the lower mantle. Piles of basaltic crust have two main attributes that influence V s : their long residence time at the CMB and their distinct composition. The latter modifies the recovered V s directly through its effect on the elastic parameters. The former is more important. Initially cold, subducted lithosphere has a high V s . As it warms, basaltic crust segregates from harzburgite and accumulates at the CMB. Stabilized by an excess density, these accumulations become hotter than the ambient mantle and develop a broad thermal halo. The resulting thermochemical structure has seismic characteristics similar to LLSVPs after accounting for tomographic filtering.
Although both former basaltic crust and depleted harzburgite are eventually returned to the upper mantle, the latter has a much shorter residence time in the lower mantle. This has an interesting implication for the evolution of mantle structure. While regions of the lower mantle become increasingly enriched with the accumulation of basalt, the upper mantle becomes increasingly depleted with the addition of harzburgite. Such processes will contribute to the development of the modern depleted upper mantle alongside the extraction of continental crust; a hypothesis also proposed by Campbell (2002). Our model predicts that LLSVPs are not fixed features of the mantle in several respects. For one, they are not merely accumulations of basaltic crust but mixtures of basalt, harzburgite, and peridotite, with the concentration of basalt increasing with depth and changing in time. Their shape and size are robust although they become less so with decreasing δ ln ρ. With lower δ ln ρ, it becomes easier for subduction to displace the accumulations of oceanic crust. There is no clear indication of the robustness of LLSVPs one way or the other; therefore, a range of dynamic behavior remains plausible.
The influence of long-wavelength thermochemical piles remains pronounced after filtering to the even lower spatial resolution of SP12RTS. However, in the presence of pPv, purely thermal models and thermochemical models produce an increase in lower mantle S/P comparable to SP12RTS. The increase is stronger in thermochemical models and proportional to the density of oceanic crust. For the highest δ ln ρ examined here, the amplitude of S/P is unrealistically large. Thus S/P, while not indicative of large scale compositional variation, as demonstrated here and first by Koelemeijer et al. (2018), restricts the maximum density contrast of thermochemical piles. Future work should examine the effect of basaltic crust on S/P in the absence of pPv and over the full range of experimentally predicted δ ln ρ.
Oceanic crust recycling also provides an explanation for seismic properties of the lower mantle on a much smaller scale. Haugland et al. (2018) mapped recycled oceanic crust as P velocity anomalies in the mantle. Their analysis predicts wave-front deflections by fragments of recycled crust. The magnitude of such distortions is sufficient to explain PKIKP precursors between epicentral distances of 130-142°in both arrival time and magnitude when δ ln ρ ≈ 4% (see Haugland et al. (2018); their Figure 6).
Our study highlights the need for large-scale chemical variation to explain the seismic properties of LLSVPs. This is in disagreement with previous work that has shown thermochemical piles to overpredict the δ ln V s amplitudes of LLSVPs (Bull et al. 2009;Davies et al. 2012Davies et al. , 2015. We suggest that a difference in how thermochemical piles form in each model could explain the discrepancy. In the models of (Bull et al. 2009;Davies et al. 2012Davies et al. , 2015, thermochemical piles form from an initial dense layer above the CMB. The volume of chemically anomalous material is chosen to be ∼ 3% of the mantle so that it is consistent with volumetric estimates for regions with a pronounced decrease in V s relative to V p (Hernlund and Houser 2008). In our models, thermochemical piles form by the steady accumulation and subsequent warming of dense oceanic crust. After filtering, the resulting V s anomaly is consistent with the size of LLSVPs. However, less than half the size of the anomalous region is chemically distinct from the ambient mantle, and only the bottom ∼ 100 km approaches pure oceanic crust. Thus, the influence of composition on the seismic properties of thermochemical piles is weakened for two reasons: a smaller portion of their size is chemically distinct and the portion (2020) 7:17 Page 12 of 16 Fig. 8 RMS δ ln V s as a function of radius for S40RTS slices (orange) and dynamic models (blue) after tomographic filtering beneath the Pacific plate. RMS δ ln V s is poorly matched for all models in the upper mantle. When δρ ≥3.1%, the amplitude of RMS δ ln V s in the lower mantle can be matched except for the small inflection at the CMB. Overall, RMS δ ln V s is smoother for the filtered models than for S40RTS that becomes diluted due to mixing with the ambient mantle. Indeed, without dilution, accumulations of oceanic crust are predicted to yield positive V s anomalies (Deschamps et al. 2012). Deschamps et al. (2012) found that the sensitivity of V s to oceanic crust in the lowermost mantle is mostly positive and that in order to induce a δ ln V s of -2.0%, typical of LLSVPs, oceanic crust would have to be 1600 K above ambient mantle temperature. They concluded that this is unrealistically high and that LLSVPs must therefore not be piles of pure oceanic crust. In contrast, we find that piles consisting of a mixture of oceanic crust, peridotite, and harzburgite require an excess temperature of only 400-700 K to explain the highest (negative) δ ln V s amplitudes of LLSVPs.
Although we find that thermal heterogeneity plays a dominant role in LLSVP formation, our study is also in contrast with previous work suggesting that the LLSVPs could be purely thermal (Schubert et al. 2004;Schuberth et al. 2009;Bull et al. 2009;Davies et al. 2012Davies et al. , 2015Koelemeijer et al. 2018) as we find that some dense compositional component is required. This discrepancy may be explained by the organizing effect of plate-motion history. In previous studies where past plate motions are imposed at the model surface, the formation of longwavelength and strong amplitude velocity structures is attributed to the organizing force of Earth's particular tectonic history and form whether compositional variations are present or not. The influence of past plate motions cannot be explored in our 2D geometry.
We cannot exclude the importance of primordial heterogeneity that may still be sequestered near the base of the mantle (e.g., Tackley 1998;and Labrosse et al. 2007). In fact, geochemical arguments based on mantle noble gas chemistry (e.g., Allègre et al. 1996;and Albarède 1998) and studies of short-lived isotope systems (e.g., Boyet and Carlson 2005;and Rizo et al. 2016) strongly suggest that primordial heterogeneity has been preserved in the Earth's mantle. From these studies, it is not clear however what the required volume for preservation is. We share the view, suggested in previous studies (e.g., Tackley 2012;Li et al. 2014;and Ballmer et al. 2016), that there likely is a combination of both primordial chemical heterogeneity and recycled oceanic crust that contributes (2020) 7:17 Page 13 of 16 Fig. 9 RMS δ ln V s as a function of radius for S40RTS slices (orange) and dynamic models (blue) after tomographic filtering beneath the African plate to the observed structure of the LLSVPs. Teasing out their relative contributions will require further interdisciplinary work combining geodynamics and geochemistry (as in CH94; van Keken and Ballentine 1998, 1999); Xie and Tackley 2004;Jackson et al. 2018) with further seismological analyses as, for example, in this paper or Haugland et al. (2018) . Clear improvements to the models that can be made include (i) compressible convection, (ii) better representation of the spherical geometry, and (iii) more flexible implementation of plate tectonics. As for (i), the use of full compressible mantle convection (e.g., Nakagawa et al. 2010;and Bossmann and van Keken 2013) will allow for the self-consistent implementation of phase changes and realistic equation of state. This will allow for an improved and more consistent comparison with the seismological observations and models. As for (ii), while we expect eventually to be able to perform full 3D simulations such as in Nakagawa et al. (2010) and Barry et al. (2017), the required high resolution and time evolution over the age of the Earth makes the use of 2D geometries still highly attractive. In BB08, a cylindrical "shrunken core" approach was used to better represent the heat loss properties of the Earth. The "spherical annulus" geometry of Hernlund and Tackley (2008) provides a more natural approach that, while 2D, represents the additional curvature by mapping of the equations onto a spherical cross section that retains the appropriate scaling between CMB radius and Earth's radius. As for (iii), we use a formulation for introducing tectonic plate into the models that is based on defining the plate geometry and computing the velocity of each plates by integrating the stress that the underlying convection exerts on them. This approach, first introduced by Gable et al. (1991), has the advantage over kinematically prescribing plates (e.g., Christensen and Hofmann 1994;and Davies et al. 2012) that the formulation is energy conservative and that one can reach Earth-like convective vigor without significant computational issues. An interesting alternative is the use of the "yield-stress" rheology (e.g., Tackley 2000) which has been demonstrated to lead to plate-like behavior without the requirement of prescribing the geometry of the plates. While it is still a question whether plate-like velocities can be modeled with this approach, recent advances using continental rafts suggest Earth-like convective vigor may be in reach (e.g., Arnould et al. 2018;Coltice et al. 2013;and Rolf et al. 2018).

Conclusions
In summary, we have shown that oceanic crust formation and recycling as modeled by BB08 can satisfactorily explain the LLSVPs in both the size and amplitude of the shear-wave anomalies as long as oceanic crust is denser than ambient mantle by 3% or more. The LLSVPs formed in our models are not composed solely of oceanic crust. Rather, they are basalt rich at their base (bottom 100-200 km) and grade into peridotite toward their sides and top. The strength of their seismic signature arises from warming over billions of years. These models appear to be well suited as base models for further exploration of the interplay between geodynamical and geochemical processes that lead to the long-term chemical and thermal evolution of the Earth. These models have been shown to match the present-day surface velocities and heat flow characteristics of the Earth, match the required outgassing efficiency of the crust and mantle as evidenced from 40 Ar, and match to a reasonable extent the HIMU-DMM-EM1 characteristics observed in OIBs and MORBs (BB08). This latter point has been further demonstrated by including the Lu-Hf isotope system (Jones et al. 2019).