Refining Heterogeneities near the Core–Mantle Boundary Beneath East Pacific Regions: Enhanced Differential Travel-Time Analysis Using USArray

Abstract

Recent advancements in seismic data analysis have enhanced our grasp of the seismic heterogeneities near the core–mantle boundary (CMB). Through seismic tomography, persistent lower-mantle structures like the large low shear velocity provinces (LLSVPs) beneath the Pacific and South Africa have been identified. However, variations in the finer-scale features across different models raise questions about their origins. This study utilizes differential travel-time measurements from the USArray, operational across the contiguous United States from 2007 to 2014, to examine the impact of upper-mantle heterogeneities on tomographic models. By averaging the P-wave travel times and calibrating them with diffracted P-waves at the same stations, we mitigate the effects of shallow heterogeneities. The findings confirm that this method accurately maps the seismic anomalies beneath the USArray, consistent with other North American studies. Calibrated Pdiff travel-time data indicate significant anomalies in the mid-Pacific and Bering Sea and lesser anomalies in the northern Pacific, aligning with the global tomographic images. Moreover, the study highlights sharp travel-time variations over short distances, such as those across the northern boundary of the mid-Pacific anomaly, suggesting a chemically heterogeneous Pacific LLSVP.

1. Introduction

The advancements in seismic data acquisition and analytical technologies have significantly deepened our comprehension of the seismic heterogeneities near the core–mantle boundary (CMB). Seismic tomographic studies have illuminated persistent long-period patterns in lower-mantle structures, notably the large low shear velocity provinces (LLSVPs) beneath the Pacific Ocean and South Africa [1,2,3,4,5]. Despite the consistency in the long-wavelength features across different models [6], the resolution of the finer structural details by tomography has led to discrepancies among these models, raising questions about the origins of such differences. These inconsistencies are often attributed to variations in the datasets and inversion methodologies. Moreover, heterogeneities in the upper mantle near the seismic receivers can distort the interpretations of deeper mantle structures [7,8,9,10], particularly affecting our understanding of the fine-scale features within the lowermost mantle [9,11]. Earlier models [3,12,13] posited LLSVPs as colossal pile-like formations, a perspective gradually refined by later studies suggesting that these could be aggregates of upwellings [12,14]. This shift underscores the potential misinterpretation of small-scale plume clusters as artifacts from upper-mantle anomalies rather than features intrinsic to deep-mantle dynamics [7,15,16,17]. Prior research has explored various approaches for correcting receiver-side upper-mantle effects on travel-time anomalies. The strategies employed include differential travel-time measurements by comparing target phases to reference phases and isolating specific phase contributions [18,19,20,21,22,23] and adjustments based on pre-existing tomographic models as reference standards [9,24]. While these methods have achieved some success, their efficacy is often constrained by the choice of reference phases or models.

This study introduces a different approach to receiver-side travel-time corrections, leveraging differential datasets from teleseismic and diffracted events recorded by the USArray. By adjusting for travel-time delays across North America, we isolate the influences of receiver-side upper-mantle structures. The resulting travel-time anomalies from diffracted waveforms, after receiver-side upper-mantle corrections, reveal persistent variations attributed to lowermost mantle heterogeneities.

2. Data Measurements and Uncertainties

Our study collects global earthquake data recorded by the USArray network, encompassing TA, FA, and REF stations, capturing events with moment magnitudes ranging from 5.8 to 7.2 spanning from 2007 to 2014 (Figure 1 and Tables S1 and S2). We remove instrument responses from the waveform data and apply bandpass filtering within two distinct frequency ranges: 0.005–0.033 Hz and 0.03–0.2 Hz. For our travel-time analysis, we utilize P and Pdiff displacement waveforms recorded on the vertical component. Selection criteria for these waveforms include a high signal-to-noise ratio (SNR > 5), defined as the maximum amplitude of a P or Pdiff wave divided by the standard deviation of the background noise within a 30 s window prior to the P or Pdiff signal [9], and a simple source time function. The differential travel-time measurement involves cross-correlating (Figure 2) the P and Pdiff data with synthetic waveforms (Equation (1)). These synthetics are generated using the AXiSEM solver [25] based on the Preliminary Reference Earth Model (PREM) [26].

$$∆{{\mathrm{T}}_P}_j^i=∆{{\mathrm{T}}_{P\_obs}}_j^i-∆{{\mathrm{T}}_{P\_syn}}_j^i;∆{{\mathrm{T}}_{Pdiff}}_j^i=∆{{\mathrm{T}}_{Pdiff\_obs}}_j^i-∆{{\mathrm{T}}_{Pdiff\_syn}}_j^i$$

(1)

where i and j are the event and station index. To ensure the reliability of our findings, we retained only differential measurements exhibiting high cross-correlation coefficients (CC) greater than 0.7 for further analysis. High CC values indicate a strong similarity in waveform shape between the observed data and synthetic waveforms. This stringent selection criterion mitigates potential ambiguities arising from complex source rupture behaviors, thereby enhancing the clarity and accuracy of our travel-time analysis.

2.1. Source-Side Ambiguity

The $∆{{\mathrm{T}}_{P\left(Pdiff\right)}}_j^i$ we calculate can be viewed as reflecting the cumulative effects along their respective raypaths (Figure 1a). To specifically address receiver-side travel-time anomalies, we adjust each event’s differential time measurements by subtracting the mean differential time of individual event (Equation (2))

$$\begin{gathered} ∆{{\mathrm{T}}_{P_c}}_j^i=∆{{\mathrm{T}}_{P_{obs}}}_j^i-∆{{\mathrm{T}}_P}_{mean}^i;∆{{\mathrm{T}}_{Pdiff\mathrm{c}}}_j^i=∆{{\mathrm{T}}_{Pdiff}}_j^i-∆{{\mathrm{T}}_{Pdiff}}_{mean}^i, \\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}∆{{\mathrm{T}}_P}_{mean}^i={\displaystyle \frac{\left({\sum }_{j=1}^n∆{T_P}_j^i\right)}{\mathrm{n}}};∆{{\mathrm{T}}_{Pdiff}}_{mean}^i={\displaystyle \frac{\left({\sum }_{j=1}^n∆{T_{Pdiff}}_j^i\right)}{\mathrm{n}}} \end{gathered}$$

(2)

This approach aims to isolate anomalies associated with the receiver, assuming that any effects from source-side structures or earthquake mislocation errors are nullified due to the similarity in raypaths near the source region. Our dataset spans a wide range of earthquake depths, from shallow to deep. Considering the depth phase’s potential impact on differential time measurements due to upward raypaths near the source region, we conducted a comparative analysis of travel-time patterns from both a shallow and a deep event. For the deep event, separating the direct wave from the depth phase is straightforward for both frequency bands (Figure 3). However, for the shallow event, the proximity between the direct phase and the depth phase makes separation challenging (Figure 3). By computing differential time measurements over varying window lengths—one that exclusively encompasses the direct phase wavelet and another that includes both the direct and depth phase wavelets—we ascertain that the inclusion of the depth phase does not significantly distort the travel-time patterns (Figure 3). This outcome, depicted in our results, suggests that depth phase inclusion is permissible for events where separation from the direct wave energy is unfeasible.

2.2. Azimuthal Variations

After adjusting for the mean of the differential travel-time measurements to mitigate source-side uncertainties, residual variations in travel-time anomalies ($∆{{\mathrm{T}}_{{P\left(Pdiff\right)}_c}}_j^i$) could still be contaminated by mid- or lower-mantle structures, particularly for raypaths traversing longer distances. Figure 4a displays the $∆{{\mathrm{T}}_{P_c}}_j^i$ across various azimuths from different seismic events during 2008 and 2009. The amplitude of each travel-time plot varies because we subtracted only the mean measurement for each event. Despite differences in station distribution between these two years, overlapping stations in the networks reveal generally consistent travel-time patterns across most regions, suggesting similarity in raypaths near the receivers. However, some regions display divergent travel-time patterns between events, likely reflecting variations in mid- or lower-mantle raypaths from events arriving from different azimuths to the USArray.

To address potential discrepancies arising from these variations, we refined our methodology by integrating travel-time measurements from approximately homogeneously distributed events across diverse azimuths (Figure 5). Specifically, 78 events were used in the low-frequency analysis, while 76 events contributed to the high-frequency investigation, revealing only slight discrepancies in the event datasets across frequencies (Supplementary Tables S1 and S2). We limited the inclusion of events with epicentral distances beyond 90° in both frequency analyses, thereby minimizing the influence of deep-mantle heterogeneities on our measurements. This integration is crucial for ensuring that data from different events are robustly combined, with all travel-time measurements for each event referenced to the same station, TA.R11A, as shown in Equation (3). Subsequently, we average the travel-time measurements for individual stations across all events (Equation (4)), ensuring that each station recorded at least seven measurements from different events, to construct maps of receiver-side P-wave travel-time anomalies (Figure 4b and Figure 6). This comprehensive coverage of receiver-side travel-time anomalies beneath the USArray enables us to effectively remove the influence of receiver-side upper-mantle effects from our diffracted-wave travel-time measurements [27].

$$∆{{\mathrm{T}}_{P_{rt}}}_j^i=∆{{\mathrm{T}}_{P_c}}_j^i-∆{{\mathrm{T}}_{P_c}}_{R11A}^i;∆{{\mathrm{T}}_{{Pdiff}_{rt}}}_j^i=∆{{\mathrm{T}}_{{Pdiff}_c}}_j^i-∆{{\mathrm{T}}_{P_{diffc}}}_{R11A}^i$$

(3)

$${Avg\left(P\right)}_j={\displaystyle \frac{\left({\sum }_{i=1}^{nevt}∆{{\mathrm{T}}_{P_{rt}}}_j^i\right)}{\mathrm{n}\mathrm{e}\mathrm{v}\mathrm{t}}};{Var\left(P\right)}_j={\displaystyle \frac{{\sum }_{i=1}^{nevt}({∆{{\mathrm{T}}_{P_{rt}}}_j^i-{Avg\left(P\right)}_j)}^2}{nevt-1}}$$

(4)

3. Results

3.1. Receiver-Side Upper-Mantle Structures

By averaging the differential travel times from a geographically and temporally diverse dataset of earthquakes recorded by the USArray, we have gained profound insights into the receiver-side upper-mantle travel-time anomalies Avg(P) across North America. In Figure 4b, we first compare our high-frequency results with the tomographic model US-SL-2014 [7], focusing on the western United States. The consistency between our travel-time patterns and the tomographic model indicates that our averaged travel-time map captures the first-order features of receiver-side structural heterogeneities. The results also align well with local tectonic features and igneous activities, which we believe are closely related to the subduction processes of the Juan de Fuca Plate off the west coast and the Farallon Plate beneath the inland areas. For instance, the low-velocity features in the southern Rocky Mountains (RMs) and the Basin and Range (B&R) area are attributed to the delamination of the Farallon Plate, which causes the upwelling of the hot asthenosphere, leading to crustal thinning [28,29]. Figure 6 further presents the Avg(P) across the entire US, yielding insights that align closely with the known geophysical and geological features of North America [30]. The delineated transition area, spanning −104° to −100° longitude, highlights a clear demarcation between the active tectonics of the Pacific plate margin with positive travel-time anomalies and the stable continental interior of the North American craton with significant negative travel-time anomalies (Figure 6). Furthermore, the subtle travel-time anomalies observed in the northeast–southwest direction in the eastern United States (Figure 6) could offer avenues for exploring the mechanisms behind lithospheric base delamination and small-scale mantle upwelling [31]. These processes, implicated in the formation of topographic and volcanic features, suggest a dynamic interaction between the lithosphere and underlying mantle that warrants further investigation. The distinct Avg(P) near the Mississippi River estuary presents an intriguing case of a frequency-dependent seismic response. The presence of positive anomalies in the high-frequency data, absent in the low-frequency analysis (Figure 6), may reflect localized heterogeneities or scattering phenomena that are more pronounced at higher frequencies. This discrepancy underscores the importance of incorporating a multi-frequency approach to seismic analysis to fully capture the spectrum of subsurface features. The variance maps (Figure 6) reveal a geographic pattern of data reliability, with lower variance in the west and higher variance in the east, culminating along the eastern coastline. The 2014 station relocation to Alaska introduces a caveat to our analysis, potentially affecting the data quality and filtering efficacy.

3.2. Lower-Mantle Structures After the Corrections

This study aims to refine our understanding of the lowermost mantle structures by mitigating the effects of receiver-side upper-mantle contributions. Unlike those conventional approaches that rely on corrections derived from pre-existing tomographic models, our corrected travel-time measurements are calculated from the averaging of all the $∆{{\mathrm{T}}_{P_{rt}}}_j^i$ values obtained from earthquakes homogenously distributed across the USArray. This method ensures minimal dependence on model-specific biases, offering a more objective assessment of the deep-mantle dynamics. To validate whether our corrected travel-time residuals ($∆{\mathrm{T}}_{cor}=∆{{\mathrm{T}}_{{Pdiff}_{rt}}}_j^i-Avg(P)$) are associated with the response from the diffracted raypaths at the CMB, we examined three events occurring in proximate regions but varying in depth. Figure S1 presents the high-frequency $∆{\mathrm{T}}_{cor}$ for these events, displaying a remarkable consistency across both the station projections (top row) and CMB path mappings (bottom row). This uniformity suggests that our method successfully isolates and suppresses those heterogeneities that are external to the diffractive path. Further analysis involved generating pairwise plots to examine the relationships between the variables in the travel-time correction process and the resultant residual values, as shown in Figure S2. Both the low- and high-frequency data exhibit a general linear relationship between the diffracted travel-time measurements and $∆{\mathrm{T}}_{cor}$. Conversely, a weak correlation between the average P-wave travel time used for the corrections and the $∆{\mathrm{T}}_{cor}$ suggests that our approach effectively nullifies the influence of shallow structural heterogeneities in the station area. Figure 7 maps out the $∆{\mathrm{T}}_{cor}$ along their diffraction paths at the CMB for both low- and high-frequency bands. While both frequency ranges reveal similar overarching characteristics, a detailed examination of the high-frequency data uncovers more nuanced structural features. The depiction of the low-velocity body in the central Pacific appears to be more fragmented and is interspersed with high-speed anomalies in the high-frequency analysis.

4. Discussion

To enhance the interpretability of our travel-time analysis, we partitioned our measurements into 2° × 2° grid boxes, calculating the mean value and standard deviations within each grid to delineate the primary travel-time patterns across our study region. A broad concordance emerges, particularly within the Pacific region (Figure 8), while comparing our findings with a range of published global P-wave and S-wave models [5,19,32,33]. The maps display key seismic heterogeneities, with the mid-Pacific low-speed body (MP), the north Pacific high-velocity body (NP), and the Bering Sea low-velocity anomaly (BS) prominently featuring in both our results and the global models [3,4,5]. High-frequency grid-averaged results reveal that the MP is subdivided into smaller, more distinct structural units. This suggests that what are often perceived as monolithic large low shear velocity provinces (LLSVPs) may actually be composites of more discrete heterogeneities, as previously observed beneath the western Pacific regions [34] or as suggested by recent geodynamic simulations [15,17]. However, a notable divergence is the absence of the low-velocity feature extending southward from the southern Bering Sea in our grid-averaged map, a structure commonly noted in many global models. This difference could be attributed to the elongated diffractive path of seismic waves traversing the north Pacific and the nearly perpendicular orientation of the north–south low-velocity anomalies relative to the direction of wave propagation. Such an orientation and path length might limit these anomalies’ impact on differential travel times, leading to their omission in our averaged results.

In our analysis of high-frequency data, the grid-averaged map delineates the northern boundary of the central Pacific heterogeneity [23,35,36,37,38,39,40], consistent with the contours highlighted by both the S-wave and P-wave models [32] and aligning closely with the P-wave boundary identified by Frost and Rost [36] through forward modeling (indicated by the arrow in Figure 9). To assess the boundary’s sharpness, we employed a double-differential travel-time calculation, ddt_AB = |dt_A − dt_B|, where A and B are neighboring stations from earthquakes with great-circle paths approximately parallel to the boundary. This analysis was particularly concentrated on ddt_AB measurements from station pairs orthogonal to the wave propagation direction. Given the unequal distribution of the stations, we only consider the station pairs with inter-station distances below 3 degrees. Figure S3 presents the $∆{\mathrm{T}}_{cor}$ and ddt_AB anomalies for the seismic event 201112140504A, displaying the data in diffractive path representations. A notable segregation into high- and low-velocity regimes is observed around a 48° azimuth, pinpointing the boundary of the Pacific LLSVP within this angular range. Consistent findings across different seismic events further reinforce this delineation. Interestingly, we identified some larger ddt_AB values within the high-velocity regions indicated by tomographic models, located significantly away from the LLSVP boundaries (Figure S3a). If the tomographic images in this region are well-resolved, then this anomaly suggests a nuanced complexity within the high-velocity zones, diverging from the expected uniformity. Alternatively, if the resolution of the tomographic models is limited in this region, our results may indicate a shift in the boundary of the Pacific LLSVP compared to previous studies [32,36]. Additionally, when plotting the $∆{\mathrm{T}}_{cor}$ against the azimuth angle for a specific seismic event (Figure 10), a distinct shift in the travel-time variability is evident. Pairs of stations at azimuths less than approximately 52 degrees exhibit greater travel-time variation than those at larger angles, indicating a pronounced difference in the seismic wave behavior across different subsurface structures.

To better characterize the spatial location of the boundary between the high- and low-velocity anomalies near the CMB, we have refined our analytical approach by implementing a weighting system for the observed travel-time variations, as detailed in Equation (5):

$$\mathrm{W}{\mathrm{d}\mathrm{d}\mathrm{t}}_{\mathrm{A}\mathrm{B}}={\mathrm{f}}^{\prime }\left(\mathrm{A}\mathrm{z}\right({\mathrm{M}}_{\mathrm{A}\mathrm{B}}\left)\right)\times {\mathrm{d}\mathrm{d}\mathrm{t}}_{\mathrm{A}\mathrm{B}}$$

(5)

where $\mathrm{A}\mathrm{z}\left({\mathrm{M}}_{\mathrm{A}\mathrm{B}}\right)$ represents the azimuth angle at the midpoint ${\mathrm{M}}_{\mathrm{A}\mathrm{B}}$ for each station pair and ${\mathrm{f}}^{\prime }$ signifies the first derivative of the best-fitting curves (depicted as red lines in Figure 10a) that model the $∆{\mathrm{T}}_{cor}$ as a function of azimuth. This derivative serves as a weighting factor that adjusts the $∆{\mathrm{T}}_{cor}$ based on the azimuthal dependence, enhancing the sensitivity of our measurements to changes along the hypothesized boundary zones. Figure 10b displays the weighted lateral travel-time variations for the seismic event 201112140504A against its diffractive path. The result highlights that the most pronounced travel-time variations cluster around an azimuth of 48°. We further explored how variations in the inter-station distances influence the detection of boundary features across seismic events. Figure S4 depicts the lateral travel-time variations mapped against their diffractive paths for event 201112140504A, incorporating a range of inter-station distances. The results show that consistent boundary features were observed at azimuths of 48° across different lateral spacings. This consistency was preserved even at the closest inter-station distance tested, approximately 1.1°, translating to a horizontal separation of about 66 km at the CMB. The spatial characteristics of the travel-time anomalies correlate closely with the dynamic simulation results of Heyn et al. [41]. Their study investigated the viscosity differences between denser thermochemical piles and the surrounding high-speed materials in the deep mantle. Their simulations predicted that, under a thermochemical evolution model—characterized by a higher internal viscosity within the piles—a short-period, asymmetric topographic depression would form around the peripheries of large-scale low-velocity bodies. Such a depression typically extends approximately 80 to 120 km in width and reaches about 4 km in depth. Further supporting evidence comes from the PcP observations along the northern boundary of the central Pacific anomaly by Frost and Rost [36]. They highlighted the significant impact of subducting slab materials on adjacent large-scale low-velocity zones, pointing to a structural orientation defined by an approximate 26° angle between contrasting velocity zones over a horizontal extent of 120 km. Drawing parallels between this evidence and our observational data, we hypothesize that the composition of the Pacific LLSVP is chemically distinct from the surrounding mantle.

5. Conclusions

Leveraging the unprecedented coverage of the USArray, operational across the United States from 2007 to 2014, this study introduces a differential travel-time observation method that significantly enhances our understanding of the deep mantle. This approach, exploiting the array’s dense and widespread deployment, is specifically engineered to minimize the impacts of upper-mantle heterogeneity, thereby highlighting the contributions that originate from the lowermost mantle. Via the analysis of average P-wave travel-time anomalies at both high and low frequencies, we have exposed the fundamental structural characteristics of the station areas, uncovering a complex tapestry of heterogeneous features beneath the United States. The differential travel-time residuals, observed at 30s and 5s periods, align with the foundational layers of the LLSVPs, revealing an intricate array of small-scale structures that have remained elusive in the traditional global seismic imaging tomography models. Our method further incorporates weighting schemes into our observations, allowing us to precisely delineate fine-scale structures such as the boundaries of LLSVPs that often evade detection by conventional seismic tomography. An interesting finding from our detailed examination includes the pronounced and stable boundary features at the northern boundary of the Central Pacific, spanning approximately 66 km. This observation lends substantial support to the hypothesis that the Pacific anomaly consists of structures with distinct chemical properties divergent from the surrounding mantle. Moreover, our regional analyses reveal the simultaneous presence of anomalous high- and low-speed travel times within the deep mantle, highlighting the crucial role of subducted slab materials in influencing mantle dynamics. However, the processes by which these materials contribute to upward thermal convection remain insufficiently understood, posing significant questions that require further interdisciplinary research to fully elucidate the complex interactions at play within Earth’s mantle.