We investigate broadband SPdKS waveforms from earthquakes occurring beneath Myanmar. These paths sample the core–mantle boundary beneath northwestern China. Waveform modeling shows that two ∼250 × 250 km wide ultra‐low velocity zones (ULVZs) with a thickness of roughly 10 km exist in the region. The ULVZ models fitting these data have large S‐wave velocity drops of 55% but relatively small 14% P‐wave velocity reductions. This is almost a 4:1 S‐ to P‐wave velocity ratio and is suggestive of a partial melt origin. These ULVZs exist in a region of the Circum‐Pacific with a long history of subduction and far from large low‐velocity province (LLVP) boundaries where ULVZs are more commonly observed. It is possible that these ULVZs are generated by partial melting of mid‐ocean ridge basalt.

Ultra‐low velocity zones (ULVZs) are thin patches of material with reduced seismic wavespeeds sitting on top of the core–mantle boundary (CMB). Most efforts at ULVZ characterization infer either a partially molten (e.g., Williams and Garnero, 1996; Revenaugh and Meyer, 1997) or Fe‐enriched compositional origin to the ULVZs (e.g., Muir and Brodholt, 2015; Wicks et al., 2017; Krier et al., 2021). Previous efforts have suggested that ULVZs may be preferentially found beneath hot spot volcanoes (e.g., Jenkins et al., 2021; Cottaar et al., 2022) or along large low‐velocity province (LLVP) boundaries (e.g., Li et al., 2017; Fan et al., 2022; Lai et al., 2022). But ULVZs are also found in the Circum‐Pacific region (Havens and Revenaugh, 2001; Thorne et al., 2019, 2020, 2021)—a region inferred to be associated with the resting place of subducted slabs (e.g., Grand, 2002). Recent studies suggest that different mineral phases have the potential to transport water to the CMB (see Walter, 2021 for a review). In this case, water released could introduce melting in the mantle or even react with outer core Fe to produce a dense low‐velocity phase (Mao et al., 2017). Mid‐ocean ridge basalt (MORB) in subducted slab crust could potentially melt when heated up at the CMB, serving as another possible source of melt (e.g., Andrault et al., 2014). Thus, a natural place to search for ULVZ genesis could be at the margins of subducted materials.

We use the SPdKS seismic phase to further investigate a ULVZ located beneath northwestern China identified previously (Thorne et al., 2020, 2021). The SPdKS phase is an S wave that intersects the CMB at the critical angle for P‐wave diffraction (Fig. 1). These legs of P diffraction are sensitive to structures along the CMB. Here, we analyze SPdKS data crossing this region using a waveform modeling approach and demonstrate that the seismic waveforms for these data are consistent with a partially molten ULVZ origin. In addition, geodynamic modeling constrained by the past 200 million years of subduction (Seton et al., 2012) suggests that MORB material is likely being transported to the CMB in this region. If this MORB material is undergoing melting, ULVZs could be generated here.

Seismic data were collected as a part of two previous studies (Thorne et al., 2020, 2021). This collection consisted of over 270,000 high‐quality seismograms for events occurring between 1990 and 2017 with event depths ≥ 75 km in the epicentral distance range from 90° to 130°. Details on data collection and processing steps are given in Thorne et al. (2020) but include standard processing steps such as (1) removing instrument response, (2) rotating to the radial component, and (3) integrating to displacement.

We selected events that crossed a probable ULVZ target zone in East Asia (Thorne et al., 2020, 2021) and identified seven events from 1995 to 2017 with Mw5.8 and depths > 75 km with simple source time functions (see Table S1, available in the supplemental material to this article). We searched for additional recent events, but no high‐quality events occurred in this region since 2017. For these events, we obtained 2307 total records. The location of events and receivers used are shown in Figure 2. Three of these events had many records spanning both a wide azimuthal and epicentral distance range, which we focus on here: (1) 3 September 2009 (19:51), (2) 13 April 2016 (13:55), and (3) 24 August 2016 (10:34). Distance profiles for these events are provided in the supplemental material (Figs. S1–S10).

We reviewed each trace in the epicentral distance range from 108° to 115° and categorized them as either anomalous or not anomalous (Thorne et al., 2019, 2020). In typical preliminary reference Earth model‐like (PREM) records in this distance range, either SKS exists alone, or at epicentral distances roughly greater than 111°, SPdKS also starts to emerge from the shoulder of SKS with a lower amplitude than SKS. Thus, in nonanomalous waveforms, we observe either 1 or 2 arrivals; and when we observe 2 arrivals, the second arrival has a lower amplitude than the first. Here, records were identified as anomalous if they had an additional arrival (e.g., 2 or 3 arrivals when only 1 or 2 are expected) or if the second arrival emerged with a larger amplitude than the first. The locations of the anomalous and nonanomalous ray paths on the CMB are shown in Figure 2a,b. The anomalous records are concentrated in two azimuthal bands with relatively normal waveforms located between them. All anomalous waveforms overlie areas of high seismic velocity in tomographic models on both source and receiver sides of the path (Fig. 2a). The concentrations of anomalous waveforms on the source side of the path are consistent with ULVZ‐likelihood calculations of Thorne et al. (2021; see Fig. 2b). For the concentration of anomalous waveforms to the east, a second high‐likelihood zone exists near the Pd path exit, perhaps indicating complex ULVZ structure on the source side of the path.

On the receiver side of the path, Pd segments are widely distributed (Fig. 2a). Following the interpretation of 3D synthetic seismograms by Thorne et al. (2021), it is likely that anomalous structure giving rise to these anomalous waveforms is concentrated near the source. Otherwise, anomalous structures would be spread out over several thousand kilometers beneath North America (see figs. 3–5 in Thorne et al., 2021).

We sorted records into 6° wide azimuthal bins, which divided areas of anomalous and nonanomalous waveforms yet also provided enough waveforms in each bin to obtain quality waveform stacks. Records were sorted into 1° epicentral distance bins, band‐pass filtered with corners between 0.025 and 1.0 Hz, and stacked using the adaptive stacking algorithm of Rawlinson and Kennett (2004).

Example waveforms for an event occurring on 13 April 2016 are shown in Figure 3. This figure shows waveforms for three different azimuthal bands that are representative of the waveform behavior encountered in this region. In Figure 3a, seismic traces for the western cluster of anomalous waveforms are represented in the azimuth band from 354° to 360°. The waveform stacks (blue) from epicentral distances 111°–115° show clear waveform differences from the PREM synthetics (green). In Figure 3b, data for the central region are shown (azimuths 0°–6°). These waveforms are similar to the PREM predictions. The eastern region, azimuths 24°–30°, shows a high degree of complexity starting at an epicentral distance of 111°. Three distinct arrivals are evident between 112° and 116°, and the waveforms are consistently distinct from PREM. In all the three azimuthal bins waveforms are matched well by PREM predictions for distances less than ∼110°, indicating a simple source time function. For distances greater than 120°, the waveforms are also similar to PREM. At large epicentral distances, SPdKS waveforms commonly look similar to PREM, which has been previously noticed in global studies (Thorne and Garnero, 2004). This similarity may be related to wavefront healing of the long Pd path around ULVZs but is not evident in all locations (e.g., beneath the Coral Sea Jensen et al., 2013), which may provide important evidence toward ULVZ size. Waveform stacks compared to PREM are shown for all azimuthal bins in the supplemental material (Figs. S11–S18) as well as individual waveforms in the most highly anomalous regions (Figs. S19, S20).

Data stacks were compared to synthetic seismograms computed with the PSVaxi technique (e.g., Vanacore et al., 2016; Krier et al., 2021). The primary model space is computed for boxcar‐shaped models (Fig. 3d) on a regular grid.

  • ULVZ thicknesses (h) ranging from 5 to 45 km in 5 km increments.

  • ULVZ length (l) in the great circle arc direction of 0.75°, 1.5°, 3.0°, 6.0°, 9.0°, and 12.0°.

  • ULVZ edge position (l1) from 5.5° to 29.5° in 1.5° increments.

  • S‐wave velocity reductions (δVS) of −10%, −20%, −30%, −40%, and −50% with respect to the PREM.

  • P‐wave velocity reductions (δVP) of −10%, −20%, −30%, −40%, and −50% with respect to the PREM, noting that we never allow δVP reductions to be larger than δVS.

Additional models exist in our model space for source‐side ULVZ geometries that have been computed for various studies (Jensen et al., 2013; Thorne et al., 2013,2019; Krier et al., 2021), giving a total of over 21,000 ULVZ models. SPdKS has minimal sensitivity to density (ρ), and in this model space, we fix it to ρ=+10% with respect to PREM (e.g., Rost et al., 2005). All models are computed at a 500 km source depth and are shifted to a common source depth (see Thorne et al., 2019).

We first compare data stacks in each azimuthal bin to synthetic seismograms from our precomputed model space. We create an empirical source wavelet for each event by stacking SKS waveforms in the epicentral distance range from 90° to 105°. Then, we do a grid search through triangle and truncated triangle functions such that when it is convolved with the PSVaxi synthetics we match the SKS waveforms for that event (see Thorne and Garnero, 2004 for a detailed description of generating the empirical source). All data stacks are aligned such that the SKS arrival starts at zero seconds. Then, we find the best shift for the synthetic seismogram to match the data stack using cross‐correlation. For each synthetic‐data pair, we compute the cross‐correlation coefficient (CC) for the best‐shifted seismogram over the interval from −5 to +30 s. The time window was chosen to include possible ULVZ‐related precursors to SKS (i.e., SPKS) and additional postcursors that exist in many of the ULVZ models. We also integrate the synthetic and data trace over the same interval and compute an integral goodness of fit (IGF) measurement defined as
(1)IGF=1.0|530data(t)dt530synthetic(t)dt|530synthetic(t)dt.
We define our goodness of fit (GF) as the CC multiplied by the integral (IGF). This is averaged over all distances for which we have data stacks, to get a single number for each synthetic model:
(2)GF(syntheticmodel)=1ni=1nIGF(i)×cc(i),
in which we have n data stacks in the azimuth window.

We use the results from the comparisons with the boxcar models as a starting point to compare data with more complicated trapezoidal ULVZ models (Fig. 3e). We introduce an additional length parameter (l3) that can lie anywhere in between the l1 and l2 edges. We allow the height to be adjustable at each angular location and have three height parameters: h1 (at the l1 edge), h3 (at l3 in between l1 and l2), and h2 (at the l2 edge). We allow any of the height parameters to be zero, allowing for a wide variety of ULVZ shapes (a sample is shown in Fig. 3f).

We computed a GF for each synthetic ULVZ model compared with each event, for each 6° azimuthal band. The highest quality event is the 13 April 2016 event, so we focus the discussion on this event, with waveform comparisons for the other two primary events shown in the supplemental material (Figs. S33–S36). We found that the models with the highest GF are consistent for the three regions, which for this event we define as: (1) the western cluster (azimuths from 348° to 360° and 0° to 6°), (2) the central cluster (azimuths from 6° to 18°), and (3) the eastern cluster (azimuths from 18° to 36°).

For both the eastern and western clusters, the best‐fit model has ULVZ model parameters of δVS=50%, δVP=50%, h = 45 km, l = 6°, and l1=5.5°. This model has elastic parameters at the extremes, but with an edge position of l1=5.5°, the point where Pd initiates on the CMB is near the far edge (the l2 edge) of the ULVZ, and there is little interaction of the SPdKS wavefield with the ULVZ. This model predicts interesting wavefield complexity similar to observed data; but the extreme properties of this model predict unrealistically large‐amplitude postcursors that we do not observe in any of our data if the ULVZ is moved such that the SPdKS wavefield has more interaction with the ULVZ (see Figs. S21, S22).

We also investigated histograms of ULVZ properties for the top‐fitting models. The histograms reveal that both eastern and western clusters are consistently explained by models with mean values of δVS=50%, δVP from −20% to 0%, h = 10 km, l = 1.5°–3°, and l1 from 8.5° to 10°. The central cluster is dominated by small‐length models with δVS=10% and δVP=0%, which are the smallest P‐ and S‐wave velocities in the model space and are identical to PREM predictions for lengths of l = 0.75° and edge positions l1=5.5°. This indicates ULVZ‐like structure in the eastern and western clusters and a PREM‐like structure in the central cluster (see histograms in Figs. S40–S46).

We refined the model further by conducting a grid search in the vicinity of the starting model with smaller step sizes between parameters. Our refined models showed ULVZ parameters in the eastern cluster with δVS=55%, δVP=14%, h = 8 km, l = 2.5°, and l1=8.5°. For the western cluster we obtained: δVS=56%, δVP=14%, h = 8 km, l = 1.75°, and l1=9.0°. Figures S23 and S24 show waveform comparisons for these ULVZ models compared to data. These models provide a good fit in the early portions of the waveform that includes SKS and SPdKS but do not fit the complexity of the SPdKS postcursor wavefield.

The added complexity of the postcursor wavefield could be fit by an arrival that emerges from the l2 ULVZ edge. To explore this possibility, we generated a series of sensitivity tests (Figs. S25–S30). An additional postcursor is generated at the correct travel time to match our observations if the l2 edge is located around 20°, but an additional negative polarity arrival gets generated from the l1 edge, which destructively interferes with the SPdKS wavefield. Thorne et al. (2021) shows that SPdKS waveforms in the eastern cluster could be interacting with two distinct ULVZs. We tested models with two boxcar ULVZs located on the source side. These types of models fit the character of the waveforms well (see Fig. S31). But similar to our previous experiments, a negative polarity arrival is generated at the l1 edge that degrades the overall GF.

To test if we can fit the complexity of the waveforms without introducing a negative polarity arrival in the middle of the SKS and SPdKS wavefield, we tested a series of trapezoidal models (see Fig. 3e,f). We computed an additional 1500 models with fixed values of δVS=55% and δVP=14% based on our boxcar model results. We allowed all three thickness parameters (h1, h3, and h2) to be 0, 5, 10, or 15 km. We allowed l1 to be 8° or 9°; l2 to be 10°, 12°, 14°, or 16°; and l3 to be 12°, 14°, 16°, 18°, 20°, or 22°.

For the western cluster, we find that the model with the highest GF has length parameters of l1=8°, l3=10°, and l2=12°. Corresponding height parameters are h1=15  km, h3=10  km, and h2=5  km. The waveform fits are shown in Figure 4a,b. For the eastern cluster, the model with the highest GF has length parameters of l1=8°, l3=10°, and l2=12°. Corresponding height parameters are h1=15  km, h3=10  km, and h2=5  km. The waveform fit between this synthetic prediction and waveforms from the eastern cluster is shown in Figure 4c,d. These models fit the anomalous SKS and SPdKS portions of the wavefield at distances around 112° but underestimate the amplitude of the third arrival for the eastern cluster. A model with l1=8°, l3=10°, l2=14° and h1=0  km, h3=10  km, and h2=10  km visually appears to fit the eastern cluster data better with respect to the third arrival (Fig. S32) but still gives amplitudes that are too low in the central portion of the wavefield. The locations of the inferred ULVZs based on the highest GF trapezoidal models are plotted in Figure 2b.

We analyzed SPdKS seismic waveforms interacting with the ULVZ structure beneath northwestern China. These data reveal two thin (∼10 km) ULVZ patches on the order of 250 × 250 km, separated by relatively PREM‐like waveforms in between. We best model these waveforms with a trapezoidal ULVZ model (in 2.5D) to capture some of the waveform complexity characteristic of these data, but some of the waveforms contain greater complexity than displayed in the models we computed (see highly anomalous waveforms in Figs. S19, S20). It is possible that some of this complexity is caused by boundary interactions with 3D ULVZ structure because these ULVZs have finite extent in the lateral direction. The 3D nature of these ULVZs is supported by analyzing the absolute SKS arrival times for the study region (Figs. S37, S38). Here, we observe SKS arrival times that are delayed in the eastern and western azimuths concerning the PREM model but are fast in the central azimuthal band. The largest SKS delays occur for the eastern cluster and are as large as 5 s. Some of the SKS delay for this cluster can be attributable to low‐seismic velocities on the receiver side of the path as is demonstrated in models of seismic tomography (Fig. S39); however, these models only predict roughly 2 s of delay, so a ULVZ as we predict in this location could account for the additional delay. The largest delays are observed for azimuths commensurate with the eastern and western clusters at the smallest latitudes, which also supports the ULVZ model occurring further to the south as inferred in this study (see Fig. 2). The 3D modeling efforts of small‐scale ULVZs suggest that three arrivals could arise out of 3D ULVZ structure with an azimuthally dependent moveout (Thorne et al., 2021). However, it is not clear from these data that such moveout exists in the third arrival. As fully 3D waveform modeling becomes more approachable future efforts should model this in full 3D (see e.g., Krier et al., 2021).

The S‐ and P‐wave velocity reductions we recovered for these ULVZs are −55% and −14%, respectively. This gives an S‐ to P‐wave velocity ratio of nearly 4:1, which is suggestive of a partially molten origin to these ULVZs (Williams and Garnero, 1996). A recent study has demonstrated that iron‐rich magnesiowüstite could also produce a 3:1 velocity ratio as an alternative to the partial melt hypothesis in certain cases (Dobrosavljevic et al., 2019). However, this is only indicated for relatively small S‐wave velocity decreases. For S‐wave velocity decreases with magnitudes larger than 30%, as we infer in this study, a 1:1 to 2:1 velocity ratio is predicted for the iron‐rich compositional case. Modeling in 3D geometries may give different elastic parameters than those we recovered here because, in the 2.5D simulations, the energy is focused differently than in 3D (Krier et al., 2021). But in considering over 20,000 2.5D ULVZ models the best‐fitting models predominantly had 3:1 or 4:1 S‐ to P‐wave velocity ratios, and thus these ULVZs are excellent candidates for a partially molten origin.

Another recent study found a ULVZ adjacent to our westernmost ULVZ using the Sdiff seismic phase (location is indicated in Fig. 2) and also inferred a large S‐wave velocity decrease of between 32% and 40% with a fixed 20 km thickness (Wolf et al., 2024). Trade‐offs exist between ULVZ thickness, S‐wave velocity reduction, and size, and hence the ULVZ imaged by Wolf et al. (2024) may have a similarly large velocity reduction if they cut their ULVZ thickness by half commensurate with our modeling. It is unclear if the ULVZ we observe in the western cluster is a distinct ULVZ or an extension of the ULVZ imaged by Wolf et al. (2024) as we lack sufficient overlapping SPdKS coverage.

Typically, partially molten ULVZs are expected to exist in the hottest part of the lower mantle, such as the interiors of the LLVPs (Li et al., 2017). However, the ULVZs in this study occur in a region far from the LLVPs and beneath long‐term subduction. Hansen et al. (2023) used tracers to track the subducted materials in global mantle convection models and showed that regions beneath subduction outside the LLVPs in the lowermost mantle could contain widespread subducted materials. Here, we further analyzed the convection model of Hansen et al. (2023). We introduce passive tracers at 300 km depth beneath subduction regions during four different time periods of 200–150, 150–100, 100–50, and 50–0 Ma. Once introduced to the model domain, the tracers are advected with mantle flow to the deeper interiors (Movie S1). We find that subducted materials sink vertically at most depths; but they spread out laterally in the lowermost mantle. The majority of subducted materials that reach the base of mantle in our study region at the present day are those that are introduced to the model domain between 150 and 100 Ma at the western Pacific subduction zones (Fig. 5a–c, panel b in Movie S1). The model shows that our study region is generally colder than regions beneath the central Pacific and Africa where the LLVPs are located (Fig. 5d) and contains widespread subducted materials in the lowermost 100 km of the mantle (Fig. 5e,f, Movie S2).

We note that the mantle convection model discussed earlier does not distinguish between subducted MORB and subducted lithospheric mantle. Mineral physics experiments have suggested that subducted MORB has higher intrinsic density than the surrounding mantle (e.g., Hirose et al., 2005) and could segregate from subducted slabs when it reaches the lowermost mantle (e.g., Tackley, 2011; Li, 2023). The subducted MORB is at the top of subducted slabs and could even directly contact the CMB before crustal segregation occurs (Li, 2023). During its migration along the CMB, the MORB could be shaped into accumulations with variable shapes (Li, 2023).

We show additional evidence that ULVZs can be found beneath areas of subduction in the Circum‐Pacific. The most anomalous SPdKS waveforms observed in this region show as much complexity as those observed for a ULVZ beneath northern Mexico, also in a region of long‐standing subduction (Thorne et al., 2019). In both locations, waveform modeling suggests 4:1 S‐ to P‐wave velocity reductions, which could indicate a partially molten ULVZ. The ULVZs in our study region could be explained by partial melting of subducted MORB. The MORB has a relatively low melting temperature and may be partially molten at the CMB (e.g., Hirose et al., 1999; Andrault et al., 2014) and cause a significant reduction of seismic velocities (e.g., Williams and Garnero, 1996).

All seismic recordings used in this study are available in the repository (doi: 10.7278/S50d-7n1m-4fdp). The supplemental material include plots of data used in this study, synthetic ULVZ predictions, and comparisons between data and synthetics.

The authors acknowledge that there are no conflicts of interest recorded.

The authors acknowledge the University of Utah Center for High Performance Computing (CHPC) for computing resources and support. M. T. and M. F. were partially supported by the National Science Foundation (NSF) Grant Numbers EAR‐1723081 and EAR‐2139966. M. F. was also supported through the University of Utah Undergraduate Research Opportunities Program (UROP). M. L. was supported by NSF Grant Number EAR‐2216564. Tomography models were obtained from the SubMachine utilities (Hosseini et al., 2018). The authors thank Deputy Editor‐in‐Chief Vera Schulte‐Pelkum, Stuart Russell, and an anonymous reviewer for constructive comments that improved this article.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.