Description: 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.

Description: The Large Low Velocity Provinces (LLVP) are two antipodal regions of reduced seismic velocity that extend about 800 km into the mantle from the core-mantle boundary. The LLVPs might affect the generation of plumes and organize large-scale plate motions. However – except for the reduced velocity – almost all properties of the LLVPs are the subject of vigorous debate. The LLVPs could simply be hot upwellings, or they could be chemically different from normal mantle. They could be a transient feature, exist since the Early Earth or be the result of continuous accumulation as a result of plate tectonics. To some extent, determining the density of the LLVPs could help to distinguish between these scenarios. However, most seismic methods are only weakly sensitive to density and so far both negative and positive density anomalies have been proposed based on seismology. A more direct means of assessing the density structure comes from inverting the gravity field. While density inversions are inherently non-unique, this can be somewhat alleviated by constraining the geometry of potential sources of the gravity anomalies. In this contribution, we use vote maps to constrain the geometry. A vote map is based on a collection of seismic tomographies and highlights areas of agreement between the seismic tomographies. We find that the LLVPs possess a slight positive density anomaly between 0.1 and 0.6 %. The variation results from how the lithosphere is treated, since we use both an isostatic model and seismically determined Moho depths, with the isostatic model resulting in smaller LLVP densities. The combination of increased density and reduced velocity can only be explained if the LLVPs are somewhat chemically different from ‘normal’ pyrolitic mantle. Using petrophysical data bases we estimate that an enrichment of 1-1.5% iron oxide content together with a temperature increase of 260 – 380 K with respect to an adiabatic temperature curve can explain the density increase and velocity decrease. Alternatively, the LLVPs would have to contain 40-60 % Mid-Oceanic Ridge Basalt and be 870 – 960 K hotter in order to explain our findings.

Description: Identifying the geodynamic processes that trigger the formation of new subduction zones is key to understand what keeps the plate tectonic cycle going, and how plate tectonics once started. Here we discuss the possibility of plume-induced subduction initiation. Previously, our numerical modeling revealed that mantle upwelling and radial push induced by plume rise may trigger plate motion change, and plate divergence as much as 15-20 My prior to LIP eruption. Here we show that, depending on the geometry of plates, the distribution of cratonic keels and where the plume rises, it may also cause a plate rotation around a pole that is located close to the same plate boundary where the plume head impinges: If that occurs near one end of the plate boundary, an Euler pole of the rotation may form along that plate boundary, with extension on one side, and convergence on the other. This concept is applied to the India-Africa plate boundary and the Morondova plume, which erupted around 90 Ma, but may have influenced plate motions as early as 105-110 Ma. If there is negligible friction, i.e. there is a pre-existing weak plate boundary, we estimate that the total amount of convergence generated in the northern part of the India-Africa plate boundary can exceed 100 km, which is widely thought to be sufficient to initiate forced, self-sustaining subduction. This may especially occur if the India continental craton acts like an “anchor” causing a comparatively southern location of the rotation pole of the India plate. Geology and paleomagnetism-based reconstructions of subduction initiation below ophiolites from Pakistan, through Oman, to the eastern Mediterranean reveal that E-W convergence around 105 Ma caused forced subduction initiation, and we tentatively postulate that this is triggered by Morondova plume head rise. Whether the timing of this convergence is appropriate to match observations on subduction initiation as early as 105 Ma depends on the timing of plume head arrival, which may predate eruption of the earliest volcanics. It also depends on whether a plume head already can exert substantial torque on the plate while it is still rising – for example, if the plate is coupled to the induced mantle flow by a thick craton.

Description: Glacial isostatic adjustment is dominated by Earth rheology resulting in a variability of relative sea-level (RSL) predictions of more than 100 meters during the last glacial cycle. Seismic tomography models reveal significant lateral variations in seismic wavespeed, most likely corresponding to variations in temperature and hence viscosity. Therefore, the replacement of 1D Earth structures by a 3D Earth structure is an essential part of recent research to reveal the impact of lateral viscosity contrasts and to achieve a more consistent view on solid-Earth dynamics. Here, we apply the VIscoelastic Lithosphere and MAntle model VILMA to predict RSL during the last deglaciation. We create an ensemble of geodynamically constrained 3D Earth structures which is based on seismic tomography models while considering a range of conversion factors to transfer seismic velocity variations into viscosity variations. For a number of globally distributed sites, we discuss the resulting variability in RSL predictions, compare this with regionally optimized 1D Earth structures, and validate the model results with relative sea-level data (sea-level indicators). This study is part of the German Climate Modeling initiative PalMod aiming the modeling of the last glacial cycle under consideration of a coupled Earth system model, i.e. including feedbacks between ice-sheets and the solid Earth.

Description: We provide the model results of the manuscript "Glacial-isostatic adjustment models using geodynamically constrained 3D Earth structures" (Bagge et al. 2020, Paper) including the (1) predicted relative sea-level and (2) applied sea-level data. The predicted relative-sea level is calculated with the VIscoelastic Lithosphere and MAntle model VILMA (Klemann et al. 2008, 2015, Martinec et al. 2018, Hagedoorn et al. 2007, Martinec & Hagedoorn 2005, Kendall et al. 2005). The glacial-isostatic adjustment models uses different Earth structures (3D, 1D global mean and 1D regionally adapted; Bagge et al. 2020, Paper; Bagge et al. 2020, https://doi.org/10.5880/GFZ.1.3.2020.004) and ice histories (ICE-5G, Peltier 2004; ICE-6G, Peltier et al. 2015, Argus et al. 2014; NAICE, Gowan et al. 2016) resulting in 44 3D models, 54 1D global mean models and 162 1D regionally adapted models. For more information on model description and input data see Bagge et al. (2020, Paper) and Bagge at al. (2020, https://doi.org/10.5880/GFZ.1.3.2020.004). The provided output data include (1a) the global distribution of predicted relative-sea level at 14 kilo years before present as ensemble range of the 3D GIA models for three ice histories as netCDF files, (1b) the predicted relative-sea level at eight locations at 14 kilo years before present for all models as ASCII file and (1c) the predicted relative sea-level for the deglaciation period for all models as ASCII files. Eight locations include Churchill, Angermanland, Ross Sea (Antarctica), San Jorge Gulf (Patagonia), Central Oregon Coast, Rao-Gandon Area (Senegal), Singapore and Pioneer Bay (Queensland, Australia). (2) The about 520 applied sea-level data provide information on time, relative sea-level and type of sea-level data. They are extracted for the eight locations from the GFZ database using SLIVisu (Unger et al. 2012, 2018) and provided as ACSII files.

Description: We provide 18 3D Earth structures on a global grid. This supplementary material of the manuscript includes (1) 18 netcdf files of the 3D Earth structures and (2) 72 figures that visualize the lithospheric thickness, lateral average viscosity of the asthenosphere, transition zone and upper mantle for all 18 3D Earth structures. The Earth structures were derived from seismic tomography models (Schaeffer & Lebedev 2013, 2010 update of Grand 2002) and, under consideration of geodynamic constrains, transferred to viscosity (Steinberger 2016, Steinberger & Calderwood 2006). The 18 Earth structures vary in conversion from seismic velocity to viscosity. Detailed description of the procedure can be found in the corresponding manuscript (Bagge et al. 2020a), where the Earth structure data were applied to the glacial-isostatic adjustment model VILMA (Klemann et al. 2008, 2015, Martinec et al. 2018) to predict the relative sea-level during the last deglaciation. The netCDF files are provided on a Gaussian grid of 256x512 grid points. Each Earth structure consist of 167 layers, while lateral variations in Earth structure are considered for 114 layers between surface and 870 km depth and radially symmetric layers are considered for 50 layers from 870 km to the Earth’s core. The Earth structure is given as logarithmic viscosity in log10[Viscosity(Pa s)]. To visualize the global 3D structures, we calculated the lithospheric thickness and average viscosity of the asthenosphere, upper mantle and transition zone. The lithospheric thickness is defined as minimum depth with a viscosity 〈 10^23.5 Pa s, the asthenosphere is defined between the base of the lithosphere and 225 km depth, the upper mantle between 225 km and 410 km and the transition zone between 410 km and 670 km depth.

Description: Glacial-isostatic adjustment (GIA) is the key process controlling relative sea-level (RSL) and paleo-topography. The viscoelastic response of the solid Earth is controlled by its viscosity structure. Therefore, the appropriate choice of Earth structure for GIA models is still an important area of research in geodynamics. We construct 18 3D Earth structures that are derived from seismic tomography models and are geodynamically constrained. We consider uncertainties in 3D viscosity structures that arise from variations in the conversion from seismic velocity to temperature variations (factor r) and radial viscosity profiles (RVP). We apply these Earth models to a 3D GIA model, VILMA, to investigate the influence of such structure on RSL predictions. The variabilities in 3D Earth structures and RSL predictions are investigated for globally distributed sites and applied for comparisons with regional 1D models for ice center (North America, Antarctica) and peripheral regions (Central Oregon Coast, San Jorge Gulf). The results from 1D and 3D models reveal substantial influence of lateral viscosity variations on RSL. Depending on time and location, the influence of factor r and/or RVP can be reverse, for example, the same RVP causes lowest RSL in Churchill and largest RSL in Oregon. Regional 1D models representing the structure beneath the ice and 3D models show similar influence of factor r and RVP on RSL prediction. This is not the case for regional 1D models representing the structure beneath peripheral regions indicating the dependence on the 3D Earth structure. The 3D Earth structures of this study are made available.

Description: The existence of mantle plumes was first proposed in the 1970s to explain intra-plate, hotspot volcanism, yet owing to difficulties in resolving mantle upwellings with geophysical images and discrepancies in interpretations of geochemical and geochronological data, the origin, dynamics and composition of plumes and their links to plate tectonics are still contested. In this Review, we discuss progress in seismic imaging, mantle flow modelling, plate tectonic reconstructions and geochemical analyses that have led to a more detailed understanding of mantle plumes. Observations suggest plumes could be both thermal and chemical in nature, can attain complex and broad shapes, and that more than 18 plumes might be rooted in regions of the lowermost mantle. The case for a deep mantle origin is strengthened by the geochemistry of hotspot volcanoes that provide evidence for entrainment of deeply recycled subducted components, primordial mantle domains and, potentially, materials from Earth’s core. Deep mantle plumes often appear deflected by large-scale mantle flow, resulting in hotspot motions required to resolve past tectonic plate motions. Future research requires improvements in resolution of seismic tomography to better visualize deep mantle plume structures at smaller than 100-km scales. Concerted multi-proxy geochemical and dating efforts are also needed to better resolve spatiotemporal and chemical evolutions of long-lived mantle plumes.

Description: Most hotspots, kimberlites, and large igneous provinces (LIPs) are sourced by plumes that rise from the margins of two large low shear-wave velocity provinces in the lowermost mantle. These thermochemical provinces have likely been quasi-stable for hundreds of millions, perhaps billions of years, and plume heads rise through the mantle in about 30 Myr or less. LIPs provide a direct link between the deep Earth and the atmosphere but environmental consequences depend on both their volumes and the composition of the crustal rocks they are emplaced through. LIP activity can alter the plate tectonic setting by creating and modifying plate boundaries and hence changing the paleogeography and its long-term forcing on climate. Extensive blankets of LIP-lava on the Earth’s surface can also enhance silicate weathering and potentially lead to CO2 drawdown (cooling), but we find no clear relationship between LIPs and post-emplacement variation in atmospheric CO2 proxies on very long (〉10 Myrs) time-scales. Subduction is a key driving force behind plate tectonics but also a key driver for the long-term climate evolution through arc volcanism and degassing of CO2. Subduction fluxes derived from full-plate models provide a powerful way of estimating plate tectonic CO2 degassing (sourcing) and correlate well with zircon age frequency distributions through time. This suggest that continental arc activity may have played an important role in regulating long-term climate change (greenhouse vs. icehouse conditions) but only the Permo-Carboniferous icehouse (~330-275 Ma) show a clear correlation with the zircon record.

Description: The spatial distribution of the geochemical domains hosting recycled crust and primordial (high‐3He/4He) reservoirs, and how they are linked to mantle convection, are poorly understood. Two continent‐sized seismic anomalies located near the core‐mantle boundary—called the Large Low Shear Wave Velocity Provinces (LLSVPs)—are potential geochemical reservoir hosts. It has been suggested that high‐3He/4He hotspots are spatially confined to the LLSVPs, hotspots sampling recycled continental crust are associated with only one of the LLSVPs, and recycled continental crust shows no relationship with latitude. We reevaluate the links between LLSVPs and isotopic signatures of hotspot lavas using improved mantle flow models including plume conduit advection. While most hotspots with the highest‐3He/4He can indeed be traced to the LLSVP interiors, at least one high‐3He/4He hotspot, Yellowstone, is located outside of the LLSVPs. This suggests high‐3He/4He is not geographically confined to the LLSVPs. Instead, a positive correlation between hotspot buoyancy flux and maximum hotspot 3He/4He suggests that it is plume dynamics (i.e., buoyancy), not geography, which determines whether a dense, deep, and possibly widespread high‐3He/4He reservoir is entrained. We also show that plume‐fed EM hotspots (enriched mantle, with low‐143Nd/144Nd), signaling recycled continental crust, are spatially linked to both LLSVPs, and located primarily in the southern hemisphere. Lastly, we confirm that hotspots sampling HIMU (“high‐μ,” or high 238U/204Pb) domains are not spatially limited to the LLSVPs. These findings clarify and advance our understanding of deep mantle reservoir distributions, and we discuss how continental and oceanic crust subduction is consistent with the spatial decoupling of EM and HIMU.