Title: Kilometer-scale structure on the core-mantle boundary at the source of the Hawaiian mantle plume

The lowermost mantle right above the core-mantle boundary shows a complex and heterogeneous landscape containing multiple poorly understood seismic features visible across a wide range of length scales. The smallest, but most extreme, heterogeneities yet observed are 15 'Ultra-Low Velocity Zones' (ULVZ), several of which have recently been linked to the base of mantle plumes. We exploit seismic shear waves that diffract along the core-mantle boundary to provide new insight into these enigmatic structures. We demonstrate that these waves have a strong frequency-dependent sensitivity to structure at different length scales above the core-mantle boundary, similar to the dispersive characteristics of surface waves. We measure a rare 20 core-diffracted signal refracted by a ULVZ at the base of the Hawaiian mantle plume at unprecedented high frequencies. This signal shows remarkably longer time delays compared to lower frequencies, indicating extreme internal variability within the Hawaiian ULVZ. Utilizing the latest computational advances in 3D synthetic waveform modeling, we are able to model this high frequency signal and constrain high-resolution structure on the scale of kilometers at the 25 core-mantle boundary, for the first time. Results reveal that the lowermost part of the Hawaiian ULVZ is extremely reduced in shear wave velocity, by up to -40%. 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.


Introduction
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 35 surface where they form hotspots and volcanoes. Over the past few decades, seismology has revealed that the landscape of the core-mantle boundary is highly heterogeneous, potentially analogous to the level of variability seen on the Earth's surface.
Ultra-Low Velocity Zones (ULVZs) represent the most extreme core-mantle boundary features 40 yet observed, generally showing shear-wave velocity reductions of 10% to 30% 1 . However, the small height of these structures -only tens of kilometers thick -means they are below the resolution limit of tomographic models, requiring higher frequency seismic waves to image them. Currently available seismic data has limited geometries that can be used to image these structures and therefore only a fraction of the core-mantle boundary is illuminated. Despite this, 45 a number of ULVZs varying in size, shape and velocity reduction have been reported 1 . Early studies interpreted ULVZs as sporadic and localized structure, but had little constraint on their lateral extent [2][3][4] . More recently four unusually large ULVZs linked to areas of hotspot volcanism have been mapped near Hawaii 5 , Iceland 6 , Samoa 7 and the Marquesas 8 . Unlike earlier studies, which suggested ULVZs represented small-scale structures, these new observations indicate the 50 presence of thin but wide ULVZs, on the order of 600-900 km across. These have been dubbed 'mega-ULVZs' in some literature. The large aspect ratio of mega-ULVZs suggest they are very dense compared to their surroundings 9 , which could be explained by iron enrichment 10,11 . 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 55 collected in hotspot lavas indicate a primordial material or a core-signature in the mantle plume 12 .
Three of the recently observed mega-ULVZs, beneath Hawaii, Iceland and the Marquesas, have been discovered using core diffracted shear-waves known as Sdiff 5,6 . Sdiff phases are observed at 60 epicentral distances over 100° from the earthquake epicenter, sometimes extending to up 140° after a long refraction path along the core-mantle boundary (Fig. 1B). Higher frequency Sdiff energy gets trapped in thin mega-ULVZs, becoming delayed and refracted. On a seismogram these guided waves appear tens of seconds after the main Sdiff phase, and for a cylindrical ULVZ show a hyperbolic delayed move-out as a function of azimuth 5 . We refer to these as Sdiff 65 postcursors. The travel time delays of Sdiff postcursors hold information on the size, shape, and average velocity reduction of the ULVZ they sample 5 . 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 Sdiff waves sample a reasonably large area of the coremantle boundary, the observation of Sdiff postcursor waves is still very rare. This is partly due to 70 the requirement of a continuous and dense array of recording seismic stations. The previous observations of mega-ULVZs can all be attributed to the deployment of the large-scale US Transportable Array 13 .

Location of the Hawaiian ULVZ 75
When Sdiff postcursors illuminate a ULVZ from a specific direction, determining its exact location has a degree of freedom along the diffracted path, and there are trade-offs between location, size, and velocity reduction in the modelling of these structures. The Hawaiian ULVZ was initially identified using Sdiff waves from deep earthquakes in Papua New Guinea recorded at the transportable array and other stations in the central United States 5 (Fig. S2). This left 80 uncertainty as to the exact location of the ULVZ in the NE-SW direction along the Sdiff 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 Sdiff waves that highlight the Hawaiian ULVZ in the N-S direction is now available. Five deep earthquakes from the Kermadec trench recorded in Alaska show postcursor energy caused 85 by the Hawaiian ULVZ ( Fig S3). Using this new data, we are now able to pinpoint the precise location of the Hawaiian ULVZ with unprecedented accuracy (Fig. 1A). The hyperbolic travel times of postcursors across this data set suggest the ULVZ is centered at 172.3°W and 15.4°Noffset further southwest from the Hawaiian Islands than previously thought 5 . Synthetic modelling of the waveforms from two directions indicate a highly axisymmetric cylindrical shape to the 90 ULVZ ( Fig. S2 and S3), similar to that proposed for the Icelandic ULVZ 6 .

High-frequency Sdiff Postcursors
Previously, Sdiff postcursors have been observed down to periods of 10s (frequencies ≤ 0.1 Hz), which limits the vertical resolution of imaged ULVZs to tens of kilometers (Fig. 1C). This is 95 insufficient to unravel the details of internal ULVZ velocity gradients, which requires the exploitation of higher-frequency observations. However, pushing Sdiff 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 exponentially with frequency 14,15 . In this study, we make the first high-frequency 100 observations of Sdiff postcursors down to 3s (≤ 0.3 Hz), allowing us to infer internal ULVZ structure on the order of kilometers.
We focus on the high-quality Sdiff data from a 2010 earthquake in the Papua New-Guinea area recorded in the contiguous United States (Fig. 1A in yellow). The Sdiff observations are filtered into two frequency bands for comparison, one from 10 to 20 s ('long-period'), and one from 3 to 105 6 s ('short-period'). The energy of the short-period phase is very weak; it is barely observable in the raw data (Fig. S1). To enhance the signal-to-noise ratio, we apply a sub-array phase-weighted beamforming technique for each station and its nearest 20 neighbors (Materials and Methods), which stacks the signals based on phase and directional coherency. Figure 2 shows the results of this phase-weighted stacking in separate time windows for the main Sdiff arrivals and postcursors, 110 each stacked along their respective incoming directions ( Fig. 2C and 2F). The long-period Sdiff stacked signals ( Fig. 2A) look similar to unstacked waveforms (Fig. S1), while the short-period postcursor stacked signals (Fig. 2D) emerge from what appears to be only background noise in the raw data (Fig. S1).
We observe that at long-periods the postcursor arrives at around 35 to 50s, varying with azimuth 115 (Fig. 2B). Strikingly, at short-periods the postcursors are significantly more delayed, arriving at 50 to 70s (Fig. 2E). Frequency dependent travel times with differences on the order of tens of seconds between long and short-period Sdiff postcursor waves has never before been documented, though a weak frequency-dependent dispersion of core-diffracted waves has been suggested based on global statistical analysis of ray parameters 16 . The main Sdiff phase and the 120 postcursor have different incoming directions (Fig. S5) due to 3D out-of-path effects. While the incoming backazimuth direction of Sdiff phases are slightly scattered ( Fig. 2C and 2F), they remain reasonably consistent across both long-and short-period measurements. Gradual deviation from 0° to 15° away from the epicenter backazimuth ( Fig. 2C and 2F) implies a strong bending effect that can only be explained by interaction with a structure of strong velocity 125 contrast (5). 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.

Modeling of the ULVZ Basal Layer
Detailed waveform modelling of the Hawaiian ULVZ based on long-period postcursor data has been shown in (5). We refine their preferred model of a 20km tall cylindrical ULVZ of radius 455km with velocity reduction of 20% Vs, to include details of finer internal structure based on our short-period observations and updated location. Initially we explore the model space of a 135 simplified cylindrical ULVZ using computational-cheap ray-based modelling (Fig. S9, S10), followed by 2.5D modelling 15 down to 3s (Materials and Methods). This reveals that the shortperiod postcursor observations can be explained by an ~2km thick layer with extreme velocity reduction (-40% Vs) at the base of the ULVZ, or by the presence of a less anomalous, but wider spread, basal-layer (Fig. S11). 140 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 1s achievable for the first time 17 . We compute 3D synthetics using a 20km thick ULVZ model with an extreme 2km basal layer of -40% dVs 145 based on our 2.5D modeling, as well as three additional models for comparison (Fig. 3C): a uniform ULVZ of -20% dVS, a gradient ULVZ varying from -10% dVS at the top to -30% dVS at the bottom, and a two-layered ULVZ with an upper 10 km at -10% dVS and a bottom 10 km at -30% dVS (Fig. 3C). The latter three models have an equivalent velocity reduction when integrated vertically across the ULVZ. Broadband (Fig. 3B) and long-period data ( Fig. S12A) are 150 unable to discriminate between these four models. Short-period data however, show strong differences in the travel times of the modeled Sdiff postcursors (Fig. 3A), with the uniform ULVZ model showing the weakest delay times and the two-layered model showing the strongest delay times. Both the gradient ULVZ and the ULVZ with a 2km extreme basal layer show good fits to the observed dispersion across both frequency bands. While not presenting a unique solution, 155 these 3D high-frequency synthetics do demonstrate that we have the ability to resolve a kilometer-scale basal-layer structure at the core-mantle boundary, and that strong velocity layering or gradients explains the observed dispersion.

Implications for Earth's early evolution history and core-mantle interaction 160
This unprecedented record of extreme seismic velocity reduction in the basal-layer of the Hawaiian mega-ULVZ sheds new light on the nature of these features and the complex processes happening at the core-mantle boundary. The steep thermal boundary just above the core is likely to explain some of the change in velocity with depth, but thermal effects can, at most, only explain several percent velocity reduction. The axisymmetric shape of the ULVZ suggests a 165 natural dynamical link between a cylindrical plume upwelling and the partial melting that might occur at its base 6 . However, the extreme shear velocity reductions of up to 40% that we observe would require a melt fraction likely above the percolation threshold 18 . Additionally, these melts are likely to be iron-rich and denser than the solid 19 , causing melts to drain out 18 .
Our velocity constraints instead suggest a compositionally distinct mega-ULVZ containing increasing iron content with depth. Iron-rich post-perovskite 20 or iron-rich magnesiowüstite 11 have been proposed as candidate minerals that could form significant proportions of solid-state ULVZs. Mixing iron-rich magnesiowüstite with bridgmenite would suggest a model with ~15% magnesiowüstite at the top to ~80% magnesiowüstite at the base 11 . Strongly iron-enriched 175 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 21 . Alternatively, the material represents a remnant of an ancient global basal magma ocean 22 . In this case, the vertical varying ULVZ structure might holds clues to changing levels of iron fractionation at different stages of magma ocean crystallisation. 180 An alternative source of iron enrichment in the lowermost mantle is through interactions with the Earth's core. The more diffuse velocity transition beneath mega-ULVZs might serve as a unique source area of mass transfer from the core. This is supported by observations of anomalous, potentially core-related, isotopic signals in hotspot lavas 12,24,25 , which imply the Hawaiian ULVZ 185 is not a closed reservoir, but is likely to be entrained in small amounts into the plume 26 . Our observation also seems to fit the proposed manifestations of hydrous minerals. In this scenario the vertically varying structure of the ULVZ might represent a chemical gradient zone, where a reaction series between water and iron is happening 27 .

Data Processing
All the data for this study is obtained from the IRIS Data Management Center. We select eight events with depths from 12 to 413km and seismic moment magnitude from 5.9 to 7.8 that sample the Hawaiian ULVZ. The detailed source information is shown in the Table S1. We remove the instrument response by convolving the original data with the inverted instrument spectrum using ObsPy 1 . 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 Sdiff 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 Sdiff data for the 2010 event filtered in 10-20s and 3-6s period bands and presented as a function of azimuth are shown in Fig. S1.

Update on Location of Hawaiian ULVZ
The previous modeled location of show Hawaiian ULVZ sits just southwest of the surface hotspot, centered roughly between 172.5W and 162.5W 2 . 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. S2), but a degree of uncertainty remained as to the exact longitudinal location due to lack of crossing data at the time of publication. The recent deployment of the Alaska transportable array provides a new direction to illuminate the Hawaiian ULVZ using diffracted phases propagating from south to north. We identify five new events from the Kermadec Islands that are recorded in Alaska which show similar hyperbolic Sdiff postcursors to the previous events recorded in the central USA (Fig. S3). The symmetry of the postcursors suggests the Hawaiian ULVZ has an axisymmetric structure likely quasi-cylindrical in shape. In this case, the least delayed postcursors represent waves that have propagated through the center of the cylinder have not been refracted out of plane and thus show no additional travel-time delay due to extra path length. Based on the intersection of the minimally-delayed Sdiff paths in two almost orthogonal directions, the ULVZ is located further to the southwest than previously thought, centered around 172.3°W and 15.4°N. Long-period synthetics for this new location are computed using coupled-CSEM 3 and compared to data from all events, shown in Fig. S2 and S3. Results show a good fit to the delay time and move-out of postcursors, although there is some variation between observed and modeled amplitudes, likely due to inaccuracies in the assumed sources and therefore Sdiff radiation patterns

Sensitivity Kernel for Sdiff Wave at CMB
The sensitivity kernel of a wave illustrates which part of the Earth affects the observed waveform. Some numerical software packages (i.e. SPECFEM_Global 4,5 ) can now calculate the finitefrequency sensitivity kernels for specific phases. However, calculations for kernels at higher frequencies (i.e. up to .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 both have free stress boundaries in the SH system. We apply the theory previously developed for Love waves in a vertically heterogeneous medium 6 to velocity profiles at the CMB to provide an estimate of the Sdiff 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 200km 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. S4 shows the eigensolutions of displacement and traction for 3s-, 10s-, and 20s-period waves. From the plot, we see that shorter periods have sensitivity closer to the CMB, with the 3s-period showing significant sensitivity to the lowermost 5 km. These kernels have guided our proposed velocity models used to fit our multi-frequency observations.

Beamforming Analysis
Beamforming is an array method used to measure the incoming direction and slowness of a signal as it passes through an array 7 . 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: where H represents the Hilbert transform on the stacked original data time series , u is the slowness vector as function of the incoming angle !, and 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 Sdiff 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 Sdiff slowness at 8.323 s/deg, as given by the IASP91 model 8 , 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 9 . We apply this procedure for each station by forming a sub-array consisting of the nearest 20 stations, or all stations within 4 degrees epicentral distance. Fig. S5 shows examples of the resulting beamforming phase stacks. We find two separate energy peaks in the time range of main arrival and postcursor respectively. These two peaks locate at different beamforming directions, suggesting these two signals pass through the array with different incident angles. The record is kept only if we observe one dominant local maximum in the estimated arrival time window. The backazimuth directions and travel times for the maxima in our beamforming stacks are shown in Fig. 2BC and EF in the main text. To estimate the uncertainty of our fixed slowness, we repeat the procedure with slowness values varied by ±5%. The differences in measured backazimuth and travel time are given as error bars on our measurements.

Waveform Stacking
We produce final waveform stacks by multiplying the linear stack with a weighted phase stack using the observed incidence angle measurements of the main phase and postcursor, such that: where represents the time series of -th station of an array, is the distance vector for i-th station with respect to the reference station, u is the measured directional slowness vector for the reference station obtained from beamforming, and denotes the instantaneous phase obtained from Hilbert transform of the original data time series . Weighting factor governs the sharpness of the noise reduction. To balance data distortion and signal coherency, we use in our processing, as recommended in 10 . We again stack over sub-arrays where N=20 nearest stations, or all the stations within four degrees epicentral distance.
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. 2A and C in the main text, with the two parts split by a dashed line. This procedure is applied for the two different frequency ranges explored, and significantly aids in bringing out a coherent postcursor signal in the 3-6 period range. For comparison to the stacked waveforms shown in Fig. 2, the raw data of event 20100320 is shown in Fig. S1.

Wavelet Spectrum of the Stacked Seismograms
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. These spectra illustrate the stronger phase dispersion in the postcursor compared to the main Sdiff arrival. We perform a continuous wavelet transform using the complex Morlet wavelet on the stacked time series from 20s before to 120s after the predicted Sdiff arrival time. Spectra of linear stacks and phase-weighted stacks based on main Sdiff and postcusor incoming angles are compared in Figure S6. The phase-weighted stacks have a better signal-noise-ratio. The spectra show that most of the signal energy is distributed above 3s. We also observe a gap in energy near the period of 6s in the data from event 20100320. It is based on these observations, that we define the two frequency ranges used in this study to demonstrate the dispersive nature of these phases: short-period at 3-6s, and long-period at 10-20s.

Clustering of Beamforming Results
The results from our beamforming backazimuth measurements show some scatter, thus we make use of a clustering algorithm (DBSCAN -density-based spatial clustering of applications with noise) to help threshold our short-period observations 11 . With no priori information on the data, the DBSCAN algorithm identifies the core sample and the outliers by measuring the neighborhood distance and density. If a point is closely connected to its topological neighbors, it will be treated as trustworthy. Conversely if a point is not reachable by many of the core points in the set within a given covering distance, it will be discarded as noise. This provides us a robust way of filtering measurements without any subjective manual data picking. The result of applying the DBSCAN algorithm to the observed deviations in backazimuth for the short-period postcursor is shown in Fig. S7. Our following misfit calculations only use the measurements deemed robust based on this method.

First Estimate of Basal ULVZ Properties
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: where r is the radius of the structure, dv is the fractional velocity reduction, and v is the background velocity at the CMB. This relationship is plotted in Figure S8 for the least delayed long-period and short-period postcursors. Results reveal that either extreme slow material as the base of the ULVZ or a ULVZ that is much wider at its base could explain the arrival times of delayed postcursors. If the bottom part of the ULVZ remains the same radius at 455 km, the velocity reduction would be up to 30%. If we otherwise assume the velocity reduction to remain constant at 20%, then the material should be spread much more widely with a radius of 700 km.

Ray-based Modeling of the ULVZ
Next, we use a ray-approximated 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 considered only where the ray enters and exits the ULVZ boundary. Figure S9 shows the refraction pattern caused by the cylindrical ULVZ and the predicted travel-times and backazimuths as a function of azimuth from the event.
With this computationally efficient tool, we are able to 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 355km to 855km in 50km steps. Figure S10 shows how well different combinations of ULVZ properties fit the short-period postcursor observations. Misfits are calculated based on the arrival times, the backazimuths, and a combination of the two. We choose to use the L1 norm in the misfit calculation, as it is more robust and resistant to outliers for a small number of samples.
In the travel-time misfit plot (Fig. S10A), we see the expected tradeoff between the velocity reduction and the ULVZ radius. The travel-time misfit suggests the cylindrical ULVZ to be either small, with an extreme velocity reduction (radius 405km, -50% dVs), or to be larger, with a more modest velocity reduction (radius 855km, -20% dVs). The backazimuth misfit also shows a tradeoff, but with the opposite orientation (Fig. S10B): requiring the ULVZ to be either small with slight velocity reduction (radius 455km, -15% Vs) or to be widespread and extremely reduced (radius 855km, -45% Vs). The combination of these two misfits normalized by each minimum value, creates a joint misfit which employs the two opposing trade-offs. This mistfit identifies the best-fitting intermediate model with a radius between 605km to 755km and a velocity reduction between 25% to 30% (Fig. S10C). The best-fit result peaks at radius 655km δ t~2 r v dv 1− dv and velocity reduction at 30%. Misfits for a ULVZ of constant radius of 455 km are less good, but would predict a velocity reduction of 40-45%.

2.5D Synthetic Exploration of Basal Layer Thickness
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 12 . 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 455km radius mega-ULVZ model with an extremely reduced basal-layer of -40% dVS, referred to as R455; 2) a similar sized ULVZ including a 655 km wide and -30% dVS basal-layer reflecting the lateral spreading of the mega-ULVZ material, referred to as R655. We analyze how the height of the basal-layer in both models interacts with Sdiff waves by varying its thickness from 0 to 5km (Fig. S11). We find that for model R455 the best fitting solution suggests a 2 km thick extreme basal-layer, while for model R655 a basal-layer thickness slightly above 2 km is suggested.

Full 3D Sdiff Synthetics
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 a small-scale 3D internal heterogeneity down to periods of 1s 13 . 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 Sdiff observations are more sensitive to the vertical rather than lateral structure. Thus we explore several variations of the R455 model with internal vertical structure. With a careful treatment of the meshing on the kilometer scale, each of such synthetics takes 128 nodes, with 64 CPU per node, 24 hours 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.
We construct four different ULVZ models for comparison: a uniform ULVZ of -20% dVS a gradient ULVZ varying from -10% dVS at the top to -30% dVS at the bottom, and a twolayered ULVZ with an upper 10 km at -10% dVS and a bottom 10 km at -30% dVS, and a 20km thick ULVZ model with an extreme 2km basal layer of -40% d based on our 2.5D modeling. In the main paper we show it is difficult to distinguish these models in broadband data, although dispersive postcursor energy can be observed (Fig. 3B). Fig. S12 shows the waveform filtered into the two period bands of interest. The 10s-20s waveforms are highly similar for all models, demonstrating that waves at these periods are unable to differentiate between models containing different kilometer-scale internal structure. Data filtered for the 3s-6s band show strong differences between the modeled waveforms. The uniform model has a postcusor around 40s after the main Sdiff phase. The postcusor for the two-layer model arrivals much later at roughly 60-70s. The postcursor of the gradient model and the 2km extreme model overlap at arrival times around 50s-60s -similar to arrivals seen in our high-frequency observations. All postcursors in these four models show a dispersive nature, to varying extents, between 3s to 10s (Fig. S13). The uniform and two-layer models show the weakest and strongest dispersion respectively, while the gradient and 2km extreme model show moderate dispersion, comparable to our observations. It should be noted that the main Sdiff phase, unlike its associated postcursor, is not dispersive.  On the left we show event radiation patterns (blue) with azimuth range for which data is shown highlighted (yellow). Next is the ray-path geometry shown between the entry (purple) and exit (green) points to the lowermost 100 km of the mantle, the event location (yellow star) and station locations (white triangles). Background shows tomography model SEMUCB_WM1 at 2800 km 14 . The right panels show the observed waveform data (black, left) and CSEM synthetics (red, right) filtered from 10s to 20s. Yellow colored bands highlight the Sdiff postcursor.