Geodynamic, geodetic, and seismic constraints favour 1 deflated and dense-cored LLVPs 2

21 Two continent-sized features in the deep mantle, the large low-velocity provinces (LLVPs), 22 influence Earth’s supercontinent cycles, mantle plume generation, and its geochemical bud- 23 get. Seismological advances have steadily improved LLVP imaging, but several fundamental 24 questions remain unanswered, including: What is the true vertical extent of the buoyancy 25 anomalies within these regions? And, are they purely thermal anomalies, or are they also 26 compositionally distinct? Here, we address these questions using a comprehensive range 27 of geophysical observations. The relationship between measured geoid anomalies and long- 28 wavelength dynamic surface topography places an important upper limit on the vertical 29 extent of long-wavelength LLVP-related density anomalies of ∼ 900 km above the core- 30 mantle boundary (CMB). Instantaneous mantle flow modelling suggests that anomalously 31 dense material must exist at their base to simultaneously reproduce geoid, dynamic to- 32 pography, and CMB ellipticity observations. We demonstrate that models incorporating 33 this dense basal layer are consistent with independent measurements of semi-diurnal Earth 34 tides and Stoneley modes. Our thermodynamic calculations indicate that the presence of 35 early-formed, chondrite-enriched basalt in the deepest 100–200 km of the LLVPs is most 36 compatible with these geodynamic, geodetic, and seismological constraints. By reconciling 37 these disparate datasets for the first time, our results demonstrate that, although LLVPs 38 are dominantly thermal structures, their basal sections likely represent a primitive chemical 39 reservoir that is periodically tapped by upwelling mantle plumes. 40 key questions around their vertical extent and composition remain un- resolved, with studies coming to divergent conclusions depending on the data they use. In 48 this study, we integrate a wide suite of recent geodynamic, geodetic, and seismological ob- servations with numerical modelling to gain a new, unified understanding of LLVPs. We find that they extend no more than 900 km above the core-mantle boundary, that the deep- est 100–200 km of these regions contain anomalously dense material, and that this material 52 most likely represents fragments of thick crust that formed early in Earth’s history ( 4 billion 53 years ago) and may have been capped by a meteorite-derived debris layer. These results suggest that, while LLVPs are mainly thermal features, their deepest portions represent a stable, early-formed chemical reservoir that is occasionally sampled by plumes of hot rock 56 ascending from the core-mantle boundary. and tide optimised by comparing an ensemble of thermodynamically self-consistent density models that cover a possible and to a suite of seismological, and geochemical constraints, we demonstrate that basal within the LLVPs likely comprises remnants of Hadean crust and chon- These results confirm that basal sections of LLVPs are potential Our work further suggests that lower if LLVPs were purely

. Geodynamic misfit as a function of input density and viscosity model. Observed and predicted dynamic topography power spectra of best-fit thermal model for TX2011 and S10 viscosity profile. Dark and light gray envelope = 99% and 50% confidence intervals for power spectrum of optimal spherical harmonic coefficients for oceanic residual depth measurements  mid-mantle (1000-2000 km), we also test R ρ = 0 in this region. In addition to using models 151 with vertically varying V S -to-density scaling factors, we construct a suite of thermochemi-152 cal models where chemical heterogeneity is represented as a density jump, ranging between 153 0.0-2.0%, between the LLVP interior and exterior. We generate density models using five 154 seismic tomographic models and perform instantaneous flow calculations using three mantle 155 viscosity profiles (Appendices A and B). Three key results emerge from this analysis. First, we find that acceptable fits to both layer (δρ c ≥ +0.8% for 13 of 15 tomographic and viscosity model combinations), but find 164 little to no excess density in the shallower 2000-2700 km depth range (δρ c ≤ +0.2% for 13 of 165 15 models; Figure 4d and 4e; Table S3). The thermochemical models also generally return The geodynamic inversions exhibit a preference for R ρ ∼ 0 throughout the mid-mantle, 174 which is too low for any plausible mantle composition and indicates that geodynamic ob-175 servables are incompatible with strong thermal buoyancy contributions from this depth.

176
Given that seismic tomographic models are dominated by l = 2 structure over the 1000-177 2000 km depth range, we explore this result further using associated sensitivity kernels for 178 instantaneous mantle flow. (e) Spectral amplitude of l = 2 VS anomalies from SEMUCB-WM1 tomographic model (French & Romanowicz, 2015). (f) Geoid kernel, K l N , coloured by geoid-to-VS anomaly correlation, rN , as a function of depth. (g) Dynamic topography kernel, K l A , coloured by dynamic topographyto-VS anomaly correlation, rA. (h) Geoid-to-topography ratio (GTR) kernel, coloured by rN .
Blue/red bands = values required to produce the observed GTR when thermal density anomalies are correlated/anti-correlated with the geoid.
The geoid-to-topography amplitude ratio (GTR) at l = 2 provides a crucial con-180 straint on the vertical extent of long-wavelength buoyancy anomalies associated with LLVPs.

181
In Figures 5a and b, we show the l = 2 components of observed non-hydrostatic geoid 182 height anomalies and water-loaded dynamic topography, which yield an estimated GTR 183 of ∼ 0.21 ± 0.07. These deflections must be caused by l = 2 density anomalies, with the 184 strongest corresponding shear-wave velocity (V S ) anomalies found within the LLVP regions, 185 the mantle transition zone, and the asthenosphere (Figure 5e). These V S anomalies are 186 anti-correlated with the observed geoid and dynamic topography, with the exception of the 187 transition zone, where V S anomalies correlate with the geoid but remain anti-correlated or 188 become decorrelated with dynamic topography (Figure 5f-g).

189
Individual l = 2 sensitivity kernels for the geoid, dynamic topography, and GTR   Goodness-of-fit to semi-diurnal body tide constraints is calculated following the method- level. By contrast, the best-fitting thermochemical density model based on the same to-232 mographic input, but with chemical heterogeneity in the base of LLVPs, yields statistically significant outcomes (95.8% significance level). 234 We predict Stoneley mode splitting functions by adapting the methodology of Koele-  Secondly, by calculating the instantaneous mantle flow associated with each model, CMB 239 deflections are dynamically consistent with each LLVP density structure rather than scaled. 240 We find that the misfit between observed and predicted Stoneley mode splitting functions 241 is ∼20% lower for the optimal TX2011-based thermochemical density model compared with 242 its equivalent thermal model (Table S4; Figure S9). This conclusion appears to contradict compared to the SP12RTS model adopted in that study (Text S1.4).    Table A1; for more discussion see   The relatively modest excess density below 2700 km recovered in our initial geody-  Figures S12, S13, and S15). Nevertheless, while each endmember composition can generate 296 density models that satisfy the 2-4% excess density threshold for long-term chemical hetero-  Table A1).

328
For all the reasons listed above, we conclude that the most likely candidate for the shifting the ULM-LLM boundary to 2800 km (see Text S1.2).

454
The second class of density models are created to investigate likely chemical compo- [0., 0.1, ..

481
In the upper 300 km of the mantle, density structure is identical to the first class 482 of models. Below 400 km, densities are taken directly from the thermodynamically self-483 consistent parameterisation described above, whilst between 300 km and 400 km depth, 484 densities derived from the two parameterisations are smoothly merged by taking their 485 weighted average, as described for the first class of models. Since optimal thermal and 486 thermochemical density models recovered from geodynamic inversions consistently find that 487 R ρ (LMM) ∼ 0, density anomalies in the 1000-2000 km depth interval are set to zero for all 488 models, although by using a high-pass filter to remove degree-two structure, we also test the 489 effect of including only small-scale density anomalies in this depth region, described further 490 in the Supporting Information (Text S1.3; Figure S14). ship between geodynamic observables (dynamic topography, geoid and CMB topography) 505 and laterally varying density anomalies across the mantle. We impose free-slip surface and 506 CMB boundary conditions. For each assumed viscosity profile, the resulting sensitivity ker-507 nels vary as a function of depth and the spherical harmonic degree under consideration.

508
Dynamic topography δA lm can then be determined using where K l A is the dynamic topography kernel, r is radius, ∆ρ 0 is the density difference where K l N is the geoid kernel, g R⊕ is surface gravity and γ is the gravitational constant.

516
CMB topography, δC lm , is determined according to where K l C is the CMB topography kernel and ∆ρ C is the density difference between the We assess model performance using a combined misfit function to assess compatibility with geoid, dynamic topography and excess CMB ellipticity constraints. Following previous studies (Steinberger & Calderwood, 2006; N. A. Simmons et al., 2009), we define the misfit to geoid and dynamic topography based on variance reduction (VR), a proxy for the proportion of observed signal explained by a given model prediction. Geoid misfit, χ N , is defined to be equivalent to 1 − VR N , where VR N represents geoid variance reduction, and is calculated globally using where N lm terms represent spherical harmonic coefficients of observed (subscript o) and predicted (subscript c) geoid, and l max = 30 is the maximum spherical harmonic degree.
Dynamic topography misfit, χ A , is defined analogously to χ N (i.e., χ A = 1−VR A ). However, since accurate residual depth measurements only exist at specific oceanic locations, rather than compare spherical harmonic coefficients, we instead determine this value in the spatial domain according to where A i terms are predicted and observed dynamic topography at N A = 2278 geographic locations (Hoggard et al., 2017), and values are weighted by the surface area of the 1 • bin associated with each data point in order to correct for latitudinal variation in sampling density. Since excess CMB ellipticity is defined using a single spherical harmonic coefficient, rather than using a variance reduction-based misfit definition, we use the expression In Figure 2, we present optimal results for the S10 viscosity profile (Steinberger et where N S = 9 is the number of individual Stoneley modes investigated, the second sum-573 mation term includes only even degree terms, where s max is the maximum order. In most 574 calculations s max = 2; however, we also test the impact of setting s max to the maximum 575 degree at which splitting function measurements are available for a particular mode, as well 576 as the consequences of adopting different misfit criteria (Text S1.4; Table S4) and ignoring 577 CMB topography (Text S2.3; Figure S12). 578 We combine χ S and χ G to yield a joint total misfit function, χ T , using