The lowermost mantle right above the core-mantle boundary is highly heterogeneous containing multiple poorly understood seismic features. The smallest but most extreme heterogeneities yet observed are ‘Ultra-Low Velocity Zones’ (ULVZ). We exploit seismic shear waves that diffract along the core-mantle boundary to provide new insight into these enigmatic structures. We measure a rare core-diffracted signal refracted by a ULVZ at the base of the Hawaiian mantle plume at unprecedentedly high frequencies. This signal shows remarkably longer time delays at higher compared to lower frequencies, indicating a pronounced internal variability inside the ULVZ. Utilizing the latest computational advances in 3D waveform modeling, here we show that we are able to model this high-frequency signal and constrain high-resolution ULVZ structure on the scale of kilometers, for the first time. This new observation suggests a chemically distinct ULVZ with increasing iron content towards the core-mantle boundary, which has implications for Earth’s early evolutionary history and core-mantle interaction.

A new paper constrains a high-resolution ultra-low velocity zone (ULVZ) structure at the base of the Hawaiian mantle. The authors further propose this to be a chemically distinct ULVZ with increasing iron content towards the core mantle boundary.

The core–mantle boundary separates the Earth’s liquid iron-nickel outer core from the solid silicate mantle. Heat from the core powers convection in the mantle, driving hot buoyant upwellings known as mantle plumes, that rise from the core-mantle boundary to the Earth’s surface where they form hotspots and volcanoes. Over the past few decades, seismology has revealed that the mantle immediately above the core–mantle boundary is highly heterogeneous, potentially analogous to the level of variability seen at lithospheric and crustal levels near the Earth’s surface^{1}.

Ultra-Low Velocity Zones (ULVZs) represent the most extreme core–mantle boundary features yet observed, generally showing shear-wave velocity reductions of 10 to 30%^{2}. However, the small height of these structures - only tens of kilometers high - means they are below the resolution limit of tomographic models, requiring higher frequency seismic waves to image them. Currently available seismic data have limited earthquake-station geometries that can be used to image these structures and therefore only a fraction of the core–mantle boundary is illuminated. Despite this, a number of ULVZs varying in size, shape, and velocity reduction have been reported^{2,3}. Early studies interpreted ULVZs as sporadic and localized structures, but found little constraint on their lateral extent^{4–6}. More recently four unusually large ULVZs linked to areas of hotspot volcanism have been mapped near Hawaii^{7,8}, Iceland^{9}, Samoa^{10}, and the Marquesas^{11}. Unlike earlier studies, which suggested ULVZs represented small-scale structures, these new observations indicate the presence of thin but wide ULVZs, on the order of 600–900 km across. These have been dubbed ‘mega-ULVZs’ in some literature^{10,11}. The large aspect ratio of mega-ULVZs suggest they are very dense compared to their surroundings^{12}, which could be explained by iron enrichment^{13,14}. Increasing evidence suggests that the correlation between the geographic location of mega-ULVZs and surface hotspots may not be a coincidence, i.e. the isotopic tungsten anomalies collected in various ocean island basalts indicate a primordial material or signatures from core infiltration in the mantle plume^{15,16}.

Three of the recently observed mega-ULVZs, beneath Hawaii, Iceland, and the Marquesas, have been discovered using core-diffracted shear waves known as S_{diff}^{7,9,11}. S_{diff} phases are observed at epicentral distances over 100°, and sometimes extending to up 140° after a long refraction path along the core–mantle boundary (Fig. _{diff} energy at higher frequencies and shorter wavelengths (on the order of the ULVZs height) which propagate closer to the core–mantle boundary, can get trapped in thin mega-ULVZs, becoming delayed and refracted. On a seismogram these guided waves appear tens of seconds after the main S_{diff} phase, and for a cylindrical ULVZ show a hyperbolic delayed move-out as a function of azimuth^{7}. We refer to these as S_{diff} postcursors illustrated in Fig. _{diff} postcursors hold information on the size, shape, and average velocity reduction of the ULVZ they sample^{7}. Here we demonstrate that the frequency content and dispersive properties of these phases also contain information on the vertical internal velocity structure of ULVZs. Although S_{diff} waves sample a reasonably large area of the core–mantle boundary, the observation of S_{diff} postcursor waves with a hyperbolic move-out is still quite rare (though evidence of S_{diff} postcursors without this characteristic, which may be caused by other lower mantle structures, have been observed across the Pacific^{3,11}). This is partly due to the requirement of a continuous and dense array of seismic stations. The previous observations of mega-ULVZs can all be attributed to the deployment of the large-scale US Transportable Array^{17}.

_{diff} waves at 96°, 100°, 110°, and 120° for 1D Earth model PREM^{49}. The dashed lines from top to bottom mark the 410 km, 660 km discontinuity, and 2791 km (100 km above the core–mantle boundary). _{diff} ray paths on the background tomography model SEMUCB_WM1 at 2791 km depth^{50}. Beachballs of events plotted in different colors including 20100320 (yellow), 20111214 (green), 20120417 (red), 20180910 (purple), 20180518 (brown), 20181030 (pink), 20161122 (gray), stations (triangles), and ray paths of S_{diff} waves at pierce depth 2791 km in the lowermost mantle used in this study. The event used in short-period analysis is highlighted in yellow. Proposed ULVZ location is shown in black circle. Dashed line shows cross-section plotted in

Sketch shows the wavefront propagating through a 900 km diameter cylindrical ultra-low velocity zone (ULVZ) (red) and caused delayed long-period (blue) and short-period (cyan) S_{diff} postcursors.

In this work, we exploit a new compilation of postcursors for seismic shear waves that diffract along the core–mantle boundary that sample the Hawaiian ULVZ. We find strongly delayed postcursors at higher frequencies than previously observed. The postcursors have strong frequency-dependent sensitivity to structure at different length scales above the core–mantle boundary. Utilizing the latest computational advances in 3D waveform modeling, we model this high-frequency signal and illuminate the internal ULVZ structure on the scale of kilometers at the core–mantle boundary.

When S_{diff} postcursors illuminate a ULVZ from a specific direction, the exact location can be anywhere along the diffracted path, and there are trade-offs between location, size, and velocity reduction in the modeling of these structures. The Hawaiian ULVZ was initially identified using S_{diff} waves from deep earthquakes in Papua New Guinea recorded at the transportable array and other stations in the central United States^{7} (Fig. _{diff} path, and calls into question the reliability of the simplified cylindrical model proposed to fit the data. With the redeployment of the transportable array to Alaska starting in 2014, increased data coverage of S_{diff} waves that highlight the Hawaiian ULVZ in the N-S direction is now available. Four deep earthquakes from the Kermadec trench recorded in Alaska show postcursor energy caused by the Hawaiian ULVZ (Fig ^{7}. Although the exact shape of the ULVZ is not well-constrained given the current data coverage, synthetic modelings of the waveforms from two directions indicate a cylindrical shape is a good first approximation of the Hawaiian ULVZ (Figs. ^{9}.

Previously, S_{diff} postcursors have been observed down to periods of 10 s (frequencies ≤ 0.1 Hz), which corresponds to a 500 km wide horizontal resolution (based on Fresnel zone half width) and limits the vertical resolution of imaged ULVZs to tens of kilometers. This is insufficient to unravel the details of internal ULVZ velocity gradients, which requires the exploitation of higher frequency observations. However, pushing S_{diff} observations to higher frequencies is challenging: arrival amplitudes are weaker, background noise from ocean waves is louder (peaking at 0.14 Hz), and the computational costs for modeling full 3D synthetic data increases as power of frequency^{18,19}. In this study, we make the first high-frequency observations of S_{diff} postcursors down to 3 s (≤0.3 Hz), allowing us to infer internal ULVZ structure on the order of kilometers (Fig.

We focus on the high-quality S_{diff} data from a 2010 earthquake in the Papua New-Guinea area recorded in the contiguous United States (Fig. _{diff} observations are filtered into two frequency bands for comparison, one from 10 to 20 s (‘long-period’), and one from 3 to 6 s (‘short-period’). The energy of the short-period phase is very weak; it is barely observable in the raw data (Fig. _{diff} arrivals and postcursors, each of which is stacked along their respective incoming directions (Fig. _{diff} stacked signals (Fig.

The time axis is aligned by the S_{diff} travel time predicted by PREM^{49}. The beamforming is performed by the predicted S_{diff} slowness at 8.32 s/°. _{diff} slowness is varied by ±5%.

We observe that at long-periods the postcursor arrives at around 35 to 50 s, varying with azimuth (Fig. _{diff} postcursor waves have never been documented before, though weak frequency-dependent dispersion of core-diffracted waves has been suggested based on global statistical analysis of ray parameters^{20}. The main S_{diff} phase and the postcursor have different incoming directions due to 3D out-of-path effects (Figs. _{diff} phases is slightly scattered (Fig. ^{7}. Combining observations of postcursor travel-time delays and backazimuth deviations, suggests that while waves at the two periods sample similar geographical regions, the seismic velocities at the different length scales above the core–mantle boundary they are sensitive to, differ significantly.

Detailed waveform modeling of the Hawaiian ULVZ based on long-period postcursor data has been shown in the previous study^{7}. We refine their preferred model of a 20 km tall cylindrical ULVZ of radius 455 km with a shear velocity reduction of 20%, to include details of finer internal structure based on our short-period observations and updated location. Initially we explore the model space of a simplified cylindrical ULVZ using computationally cheap ray-based modeling (Figs. ^{19} down to 3 s (“Methods”). This reveals that the short-period postcursor observations can be explained by an ~2 km thick layer with extreme velocity reduction (40%) at the base of the ULVZ, or by the presence of a less anomalous, but wider spread, basal layer (Fig.

While it is still very challenging to simulate full 3D ULVZ synthetics at the high frequencies we explore here, a recent method development combining wavefield injection with AxiSEM3D makes full 3D global ULVZ synthetics down to 1 s achievable for the first time^{21}. We compute 3D synthetics using a 20 km thick ULVZ model with an extreme 2 km basal layer of 40% shear velocity reduction, based on our 2.5D modeling (Mesh details in Fig. _{diff} postcursors (Fig.

Synthetics are made for a uniform ULVZ model of 20 km with −20% _{diff} travel time predicted by PREM^{49}. The waveforms of the short-period postcursor are filled to emphasize them. Note that the gradient model (blue) and 2 km extreme basal model (purple) nearly overlap. The faded yellow shading indicates the time delay ranges for the long- and short-period postcursors as observed in the real data. (C) Cartoons of four different ULVZ models tested in the synthetics.

This unprecedented record of extreme seismic velocity reduction in the basal layer of the Hawaiian mega-ULVZ (Fig. ^{9}. However, the extreme shear velocity reductions of up to 40% that we observe would require a melt fraction likely to be above the percolation threshold^{22}. Additionally, these melts would likely be iron-rich and denser than the solid^{23}, causing them to drain out^{22,24}.

Our velocity constraints instead suggest a compositionally distinct mega-ULVZ containing increasing iron content with depth. Iron-rich post-perovskite^{25} or iron-rich magnesiowüstite^{14} have been proposed as candidate minerals that could form significant proportions of solid-state ULVZs. Mixing iron-rich magnesiowüstite with bridgmanite would suggest a model with ~20% magnesiowüstite at the top to ~70% magnesiowüstite in the extreme basal layer^{14}.

ULVZs have been hypothesized to represent remnants of an ancient global basal magma ocean^{26}. If this is the case, the vertical varying ULVZ structure observed holds clues to changing levels of iron fractionation at different stages of magma ocean crystallization. Alternatively, strongly iron-enriched compositions in the lowermost mantle may trace Earth’s early impact history, when iron-rich impactor cores mixed with the silicate mantle and accumulated at the core–mantle boundary^{27}. An alternative source of iron enrichment is through interactions with the Earth’s core. The extreme anomalous composition on the CMB might allow compositional mass transfer from the core^{28}. It is possible that the extreme basal layer in our model represents the first observation of crystallization products formed by core exsolution—a process that is suggested to have driven the early dynamo^{29–31}. Core products infiltrating the mantle is supported by observations of anomalous, potentially core-related, isotopic signals in hotspot lavas^{15,16,32}, which imply the Hawaiian ULVZ is not a closed reservoir, but is likely to be entrained in small amounts into the plume^{33}.

The observation of this extreme anomaly in the deepest mantle, not only pushes the boundary on the degree of chemical heterogeneity present, but also implies a strong variation in heat flux across the CMB, suggesting greater potential for mechanical and electromagnetic coupling^{28}. These are important boundary conditions needed to understand the convective flows generating the geodynamo and mantle plumes^{34–37}.

All the data for this study are obtained from the IRIS Data Management Center. We select seven events with depths from 12 to 413 km and seismic moment magnitude from 5.9 to 7.8 that sample the Hawaiian ULVZ. The detailed source information is shown in the Table ^{38}. Horizontal components are rotated to the radial and tangential orientations. In this study we analyze the tangential component of the seismogram, because the SH component of the S_{diff} wave propagates further and attenuates less than the SV component and thus has stronger energy. Low-quality noisy traces are manually identified and removed. We use a zero-phase fourth-order Butterworth band-pass filter for filtering. Examples of the tangential component S_{diff} data for the 2010 event filtered in 10–20 s and 3–6 s period bands are presented as a function of azimuth are shown in Fig.

Events used in this study. Source information is obtained from the global CMT catalog.

Event | Region | Date and time (UTC) | Latitude (°) | Longitude (°) | Depth (km) | Magnitude (M |
---|---|---|---|---|---|---|

20100320 | New Ireland Region, P.N.G. | 2010-03-20 14:00:50 | −3.32 | 152.33 | 413 | 6.6 |

20111214 | Eastern New Guinea Reg., P.N.G. | 2011-12-14 05:04:59 | −7.49 | 146.83 | 133 | 7.1 |

20120417 | Eastern Neaw Guinea Reg., P.N.G. | 2012-04-17 07:13:49 | −5.66 | 147.16 | 209 | 6.8 |

20180518 | South Of Kermadec Islands | 2018-05-18 01:45:31 | −34.67 | −178.21 | 14.3 | 6.1 |

20181030 | North Island, New Zealand | 2018-10-30 02:13:39 | −39.07 | 174.94 | 226 | 6.1 |

20161122 | North Island, New Zealand | 2016-11-22 00:19:43 | −40.79 | 177.58 | 12.0 | 5.9 |

20180910 | Kermadec Islands Region | 2018-09-10 04:19:02 | −31.91 | −179.13 | 119 | 6.9 |

(

The previous modeled location of the Hawaiian ULVZ sits just southwest of the surface hotspot, centered roughly between 172.5 W and 162.5W^{7}. With the diffracted data from events in Papua New Guinea recorded at the transportable array in the central USA, which propagate mainly from west to east, the ULVZ location was well-constrained latitudinally (Fig. _{diff} postcursors to the previous events recorded in the central USA (Fig. _{diff} paths in two almost orthogonal directions (using 20100320 and 20180910 events), the ULVZ is located further to the southwest than previously thought, centered around 172.3°W and 15.4°N (Fig. ^{39} and compared to data from all events, as shown in Figs. ^{40} and the S_{diff} radiation patterns. The reduction of the main phase energy is largely underestimated in the synthetic data.

While some of these events have a hint of a high-frequency further delayed postcursors, none were of convincing quality to allow interpretation. However, the overall noisy data and presence of depth phases, do not allow us to definitively confirm the absence of high-frequency postcursors from this dataset.

The sensitivity kernel of a wave illustrates which part of the Earth affects the observed waveform. Some numerical software packages (e.g. SPECFEM3D_GLOBE^{41,42}) can now calculate the finite-frequency sensitivity kernels for specific phases. However, calculations for kernels at higher frequencies (i.e. up to 0.33 Hz) are still very challenging given current computational resources. Here we approach the high-frequency sensitivity kernels with a more heuristic analytical method. We note that the tangential component of the guided postcursor shear waves we observe near the CMB are analogous to surface Love waves, as they both have free stress boundaries in the SH system. We apply the theory previously developed for Love waves in a vertically heterogeneous medium^{43} to velocity profiles at the CMB in order to provide an estimate of the S_{diff} sensitivity kernel at different frequencies. We assume the wave coming from the mantle side is a quasi-plane wave of a specific frequency and wavenumber. We then extend the wave propagation from the lower 200 km to the CMB using the propagator matrix method. Making use of the boundary condition at the CMB, we obtain the eigen wavenumber for each specific frequency and then transform them into a sensitivity kernel. Fig.

We also estimate the first Fresnel zone widths by computing the width at the core–mantle boundary for which arrivals at the stations will constructively interfere (i.e. arrive within wave period/4). For a 10–20 s period we estimate the half width to be 200–300 km, and for a 3–6 s wave 100–200 km. Note that these widths are less than the radius of the ULVZ.

Beamforming is an array method used to measure the incoming direction and slowness of a signal as it passes through an array^{44}. We use beamforming not only to determine the direction of the incident wave, but also to enhance the energy of the original signal by stacking data using the measured incident direction to align specific phases. The beamforming stack (B) is given by:_{i}(_{i} is the distance vector to the reference station.

First, we apply this beamforming stacking on the original array data in order to determine the most likely incoming angle. Unlike other body waves, the S_{diff} phase has a fairly constant slowness value for the grazing distances as long as the velocity variations are negligible at the diffraction exit points at the core–mantle boundary. We fix the absolute value of S_{diff} slowness ^{45}, and search the incident direction in a range of −50 to 50 degrees with respect to the incoming angle of the great circle path from the event epicenter provided by the TauP software^{46}. We apply this procedure for each station by forming a subarray consisting of the nearest 20 stations, or all stations within 4 degrees epicentral distance.

Figure

We note that 3D out-of-path effects and velocity variations at the turning point of the seismic ray have the potential to cause slowness variations leading to significant mislocation in plotted slowness vectors. To test whether this has an impact on our beamforming results we perform a sliding-window F-K analysis of the S_{diff} main arrival and postcursor allowing slowness to vary from predicted values. An example of this test using a subarray of multiple stations surrounding seismic station TA.T33A is shown in Fig.

We also retested our beamforming result allowing for greater variation in slowness of ±10% (slowness differences 0.8 s/deg) and compared this with our initial ±5% slowness variance calculations, as shown in Fig. _{diff} signal returns highly similar results, compared to our original ±5% variance analysis. We find that while there is some small tradeoff between the backazimuth, arrival time, and the absolute slowness value, the influence of absolute slowness has a minimal influence on results.

We produce final waveform stacks by multiplying the linear stack with a phase-weighted stack using the observed incidence angle measurements of the main phase and postcursor, such that_{i} represents the time series of _{i} is the distance vector for _{k} denotes the instantaneous phase obtained from Hilbert transform of the original data series _{i}(^{47}. We again stack over sub-arrays where

Because of the differences between our main arrival and postcursor phases, we have to split our stacked seismogram into two parts: the first window uses the measured main arrival incoming direction and the second window uses the measured postcursor incoming direction. The final stacked waveforms for all stations are shown in Fig.

The dispersive nature of the observed postcursors is difficult to observe in seismograms in the time domain. Here we use the wavelet spectrum, which is well-suited for non-stationary seismic data, to show the energy distribution of signals in both the frequency domain and the time domain^{48}. These spectra illustrate the stronger phase dispersion in the postcursor compared to the main S_{diff} arrival. We perform a continuous wavelet transform using the complex Morlet wavelet on the stacked time series from 20 s before to 120 s after the predicted S_{diff} arrival time. Spectra of linear stacks and phase-weighted stacks based on main S_{diff} and postcusor incoming angles are compared in Fig.

Although calculating wave propagation through a 3D ULVZ model is a complicated process, we first estimate the basic properties of the base of the ULVZ using simple calculations of postcursor time delay. Along profile that samples the ULVZ center, multi-pathing effects are negligible. This provides us with a simple geometrical relationship that links the properties of the ULVZ with the arrival time delays:

Next, we use an approximate ray-based method and the travel-time and backazimuth constraints from the refracted short-period postcursors to further estimate the properties of the lowermost part of the ULVZ structure. We only predict the horizontal wave propagation using the horizontal slowness and assume a full decoupling of the vertical propagation. The transmission of the wave is calculated using Snell’s law only where the ray enters and exits the ULVZ boundary. Figure _{diff} energy in the real data from 50° to 65° azimuth which is absent in this simplified simulation that only provides a first order estimate of S_{diff} postcursor behavior.

With this computationally efficient tool, we can explore a large parameter space of radius and velocity reductions for the ULVZ, before computing more expensive full waveform synthetics. We implement a grid search varying the ULVZ velocity reduction from 15 to 50% in steps of 5%, and its radius from 355 km to 855 km in 50 km steps. Figure

In the travel-time misfit plot (Fig.

We build on initial ray-tracing based results to further explore the causative effect of the ULVZ lateral structure using AxiSEM 2.5D synthetics down to 3s^{18}. These synthetics allow us to capture the full finite-frequency sensitivity of the diffracted phases but are more computationally feasible in terms of exploring the model space than full 3D simulations. 2.5D models assume azimuthal symmetry and are only accurate in the event-station plane, such that out-of-plane effects cannot be captured. Thus, we compare our synthetic results to the waves that have propagated through the center of the ULVZ without refracting, i.e. the postcursors waves with minimal delay times. These waves show a travel-time deviation of 12 s between the long- and short-period postcursors. We construct two end-member models for the mega-ULVZ structure: (1) a 455 km radius mega-ULVZ model with an extremely reduced basal-layer of −40% _{diff} waves by varying its thickness from 0 to 5 km (Fig.

Using our estimates of ULVZ properties from the above analyses, we create models for full 3D waveform modeling. A new hybrid method has recently been developed using a combination of wavefield injection and AxiSEM3D, which allows computation of the global wavefield with small-scale 3D internal heterogeneity down to periods of 1s^{21}. This method stores the wavefield information at the injection boundary and then uses AxiSEM3D to compute the wavefield inside and outside the boundary. As demonstrated in the 2.5D synthetics, our S_{diff} observations are more sensitive to vertical rather than lateral structure. Thus, we explore several variations of the R455 model with internal vertical structure. With a dense vertical mesh to capture our input model on the kilometer scale, each of such synthetics takes 128 nodes, with 64 cores per node, 96 h to finish on the Cambridge HPC cluster. Testing the wider R655 ULVZ scenario would require a larger injection boundary, and thus even more computing resources.

In Fig.

We construct four different ULVZ models for comparison: a uniform ULVZ of −20% _{diff} phase. The postcusor for the two-layer model arrives much later at roughly 60–70 s. The postcursor of the gradient model and the 2 km extreme model overlap at arrival times around 50–60 s—similar to arrivals seen in our high-frequency observations. All postcursors in these four models show a dispersive nature, to varying extents, between 3 s and 10 s (Fig. _{diff} phase, unlike its associated postcursor, is not dispersive.

The authors acknowledge the IRIS (Incorporated Research Institutions for Seismology) data center for providing all the data. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service. We thank D. Al-Attar for input on the kernel calculation, Carl Martin for the wavefront simulations that Fig.

Z.L.: Conceptualization, Methodology, Software, Formal Analysis, Writing—Original Draft. K.L.: Methodology, Software. J.J.: Visualization, Writing—Review & Editing, Supervision. S.C.: Conceptualization, Methodology, Writing—Review & Editing, Supervision, Funding acquisition.

The raw seismic waveform data are archived and available from the IRIS (Incorporated Research Institutions for Seismology) data center (

Code of simulating the 3D synthetics waveform modeling in this study is available at

The authors declare no competing interests.

SUPPLEMENTARY INFO

Peer Review File

The online version contains supplementary material available at