Body Waves Retrieved From Noise Cross‐Correlation Reveal Lower Mantle Scatterers Beneath the Northwest Pacific Subduction Zone

Seismological studies have revealed heterogeneities at scales of 10–1,000 km in the lower mantle. For the first time, we demonstrate here that seismic interferometry of ambient noise can be used to detect faint scattered waves. We locate scattering structures at depths of ~900–1,000 km from P‐to‐P scattered wave signals in stacks of vertical cross‐correlograms from dense seismic arrays in northeast China. Independently, we show that this scattering structure produces S‐to‐P converted signals in the recordings of a deep earthquake. We constrain the thickness (H), the contrasts of the shear wave speed (δVs), and density (δρ) of the anomaly using 1‐D synthetics. Our modeling of both the reflected and converted signals indicates that H = 10–20 km, δVs = −7.2%, and δρ = 0.6%, which are similar to values reported previously and consistent with the interpretation that the scattering structure is a fragment of the basaltic oceanic crust associated with the subducted Izanagi plate.


Introduction
The lower mantle of the Earth is heterogeneous on various spatial scales. Global seismic tomography images have detected large-scale heterogeneities in the lower mantle that are believed to be related to subducted or delaminated slabs (Fukao & Obayashi, 2013;Tao et al., 2018). Small-scale heterogeneities in the mid-tolower mantle provide important insights into fundamental problems such as the circulation of mantle material, the mineral composition, and chemical evolution in the deep earth (Bentham & Rost, 2014;Brandenburg et al., 2008;Haugland et al., 2017;Kaneshima, 2016). Different seismic methods have been applied to detect scatterers/reflectors in the mid-lower mantle, such as studies of PP/SS/P′P′ precursors, receiver functions, and S-to-P conversions (Kaneshima, 2016). Reflectors have been detected globally in the depth range of 800-1,300 km, and they are located in seismic fast, slow, and neutral velocity domains (Waszek et al., 2018) with a lateral scale of 500-2,000 km. Converted phases from a sharp discontinuity at a depth of~1,000 km have been detected beneath Western Europe and Hawaii using receiver function analyses (Jenkins et al., 2017;Shen et al., 2003). In subduction regions, the S-to-P conversion phase of deep-sourced earthquakes has been widely applied and fine-scale heterogeneities with length scales on the order of~10 km have been mapped within a depth range of 800-1,800 km (He & Zheng, 2018;Kaneshima, 2019;Kawakatsu & Niu, 1994;Li & Yuen, 2014;Niu, 2014). The application of these methods and the inferred distribution of mid-mantle heterogeneities, however, are limited by the ray geometry confined by the earthquake distribution.
Over the past decade, it has become common to retrieve Green's functions by cross-correlating the continuous seismic noise or earthquake coda wave. The extraction of surface waves from ambient noise has been widely used to image the structure of the crust and the mantle, and theoretical studies have demonstrated the effectiveness of this method in diffuse fields with uniformly distributed noise sources (Lin et al., 2008;Luo et al., 2012;Shapiro & Campillo, 2004;Yang et al., 2007). Body-wave extraction is a greater challenge due to its weaker energy. With the deployment of dense seismic arrays and using the stacking strategy, however, body waves reverberating in the crust and reflecting off the 410-km and 660-km discontinuities have recently been successfully retrieved (Feng et al., 2017;Pedersen & Colombi, 2018;Phạm et al., 2018;Zhan et al., 2010). Even reflections from the core-mantle boundary (CMB) and the phases through the core have been reconstructed (Poli et al., 2015;Wang et al., 2015).
In this study, we apply the cross-correlation technique to seismic ambient noise recorded by dense seismic arrays in northeast China, and we retrieve the body-wave signals, which are confirmed to be reflected from scatterers at depths of~900-1,000 km beneath the Northwest Pacific subduction zone. To our knowledge, these are the first reported scatterers located in the mid-lower mantle retrieved using the seismic interferometry technique. Our results show great potential to improve the detection of small-scale heterogeneities in the lower mantle by taking advantage of the densely distributed broadband seismic arrays.

Data and Processing
We chose the Northwest Pacific subduction zone as our study region. A stagnant slab has been imaged lying subhorizontally in the mantle transition zone (MTZ) (Fukao & Obayashi, 2013). The topography of the 660-km discontinuity has been investigated in detail by numerous studies of converted/reflected waves (Ai et al., 2003;Li et al., 2008;Liu et al., 2015). In addition, source-sided S-to-P conversion waves have been investigated and have confirmed the existence of heterogeneities in the lower mantle within a depth range of 900-1,200 km (Kaneshima, 2009;Li & Yuen, 2014;Niu, 2014). The deployment of numerous temporary and permanent seismic stations in northeast China provides an unprecedented opportunity to retrieve body waves related to deep mantle structures or small-scale heterogeneities in the lower mantle from ambient noise.
The continuous seismic records used in this study were primarily provided by the NECsaids Array (Wang et al., 2016), deployed intermittently in northeast China from August 2010 to September 2017. The array comprises 60 stations, primarily along with two nearly orthogonal dense profiles, while the other 25 stations were installed in the surrounding area. In addition, we collected contemporaneous continuous records from stations in the China Digital Seismic Network (CEAarray) (Zheng et al., 2010) and NECESSArray (Tao et al., 2014). A total of 125 stations were used in this study with the average station spacing estimated to be~30 km, with station spacing as small as 15-20 km just adjacent to the E-W profile shown in Figure 1a. The densely distributed stations increase the probability of extracting body-wave reflections using the noise correlation technique.
Our basic seismic data processing followed the procedures described by Bensen et al. (2007) and . For each station, the vertical components of the continuous seismic waveforms were cut into 1-day-long segments after being decimated to five samples per second. Then, all the daily segments were band-pass filtered between 0.1 and 0.5 Hz after pre-processing, for example, removing the data mean and trend, the deconvolution of the instrumental response, the time domain normalization, and the spectral whitening. The positive and negative lags were folded for each cross-correlation to obtain the so-called symmetric component.
In individual station-pair cross-correlation curves, the surface wave is always the dominant phase, while the body-wave signals can barely be identified (Figure 2a). We only retained station-pairs with interstation distances less than 200 km to reduce the possible influence of surface waves on the P410P phase, and we normalized each cross-correlation curve with its maximum amplitude over the following 80-240 s (Figure 2b) after muting the time window of 0-80 s. We conducted normal move-out corrections with IASP91 model (Kennett & Engdahl, 1991) as a reference, same to the scheme employed by Feng et al. (2017), to enhance the coherence of the known reflection phases of P410P and P660P (Figure 2c). To observe variations in the deeply reflected phases, the corrected zero-offset traces were stacked via common mid-point stacking within different circular areas around the latitudes 43-44°N. We defined the radius of the stacked bins as 0.5°as applied by Feng et al. (2017). We also defined the mid-point of the station-pair as the reflection point. The corrected zero-offset traces with reflection points in the same bin were stacked using the nonlinear phase weighted stacking (PWS) method (Schimmel & Paulssen, 1997).

X-Phase Retrieved From Cross-Correlations
Most stacked cross-correlation curves show clear phases around the times of 100 and 150 s, and some have another signal around 200-220 s (Figures 2e and 2f). We show two typical examples with the stacked bin located adjacent to detected clusters of mid-mantle scatterers near the Western Japan Sea (Regions II and III in Figure 1a). The detected mid-mantle scatterers are distributed sporadically in this region at depths of 900-1,200 km (Li & Yuen, 2014;Niu, 2014), and this information can be used to facilitate the detection and confirmation of new signals. For comparison, we chose another area (Region I in Figure 1a) at a similar latitude but far from the hinge of the slab (Li et al., 2008) to avoid the potential complicated effects, such as the significant depression of the 660 or the effects caused by the slab interface of dipping portion. The two peaks appearing at~100 and~150 s in the stacked waveforms ( Figure 2f) correspond to the body waves reflected at the 410-and 660-km discontinuities, referred to as P410P and P660P, respectively, which have been detected and explored in recent studies (Feng et al., 2017;. The slowness is around zero for both phases in the vespagrams calculated by Radon transformation (Schultz & Gu, 2013) (Figures 2d, S4b, S5b, and S6b), confirming the phases are reflected at the deep mantle discontinuities. However, one intriguing and unknown signal, which we label the X-phase occurs at~218 and~200 s for the II and III curves, respectively, with amplitudes 22% and 12% of the P410P, one of the reference phases used in our study.
Because we use the vertical components to calculate the cross-correlations, it is likely that the X-phase is a P-polarized body wave. The most reasonable candidate is a PP phase, which is reflected at a discontinuity or scatterer in the deeper mantle. We estimate the 95% confidence interval by bootstrap resampling 70% of the traces per stack (Figures 2e, S4c, S5c, and S6c), which indicate that the peaks of X-phase are significant. Meanwhile, the slowness being close to 0 (Figures S5b and S6b) also confirms that the X-phase is a reflected wave. Based on the IASP91 model, we calculated the reflection depths of the two phases to be~1,008 and 930 km, which is consistent with scatterers previously reported at depths of 987 and 937 km (Li & Yuen, 2014;Niu, 2014) located both at the edges of the red and orange circles in Figure 1a. In the following section, we indicate and give new evidence for the origin of the X-phase using teleseismic S-to-P converted waves (the SxP phase).

Scatterers Confirmed by SxP Converted Phase
Source-sided S-to-P conversion waves are widely used to detect mid-lower mantle heterogeneity beneath subduction regions. The large-scale seismic array stacking technique enhances the coherent later arrivals, Regions with P-to-P phases reflected by the scatterers in the deep mantle retrieved from ambient seismic noise are marked with red and orange circles (regions II and III, respectively). Region I, which is far from the hinge of the subducting slab, was selected as a reference. The black crosses and gray ovals mark the locations of the scatterers detected by Li and Yuen (2014) and Niu (2014) which facilitates the identification of the mid-mantle scatterers (e.g., Kawakatsu & Niu, 1994). We searched for deep earthquakes and teleseismic records that could be used to map landward scatterers in the lower mantle using the S-to-P converted phase. We found an earthquake that occurred on 5 April 2013 recorded by 98 stations in Europe (Code: CH, FR, G, GR, GU, IV, MN, NI, OE, and SI) with epicentral distances between 75°and 80°that was appropriate for this purpose (Figure 1b). We aligned the direct P phase to 0 s and filtered the data between 0.2 Hz and 1.0 Hz (Figure 3a). Clear phases can be observed in the window between 26 and 29 s. We then performed a Radon transformation (Schultz & Gu, 2013) to obtain a vespagram (Figure 3b). In the calculated slowness-time plot, energy consistently arrives at~27.1 s following the direct P wave. The slowness of this phase is slightly larger, with a value of 0.02 s/deg, than that of the direct P wave, indicating a subhorizontal structure beneath the source. An array beam-forming analysis of the data shows that the back azimuth of the S-to-P phase deviates −1°from the great-circle path defined by the source and the receiver (Figure 3c), almost the same direction as the direct P wave. The 1-D structure modeling explains the observed slowness and travel time adequately (Figure S1), indicating a relatively simple 1-D structure beneath. The conversion depth is estimated to be~900 km, with the conversion points located 36 km southeast of the center of the orange circle ( Figure 1a).
Due to the consistency in the locations of the conversion points of the S-to-P waves and the reflection points of the X-phase retrieved from the ambient noise, it is reasonable to speculate that the X-phase is caused by a structure of one or two spatially very close scatterers in the mid-mantle at a depth of~900 km. Based on this assumption, we took advantage of the combination of reflected and converted phases obtained from the ambient noise and teleseismic records to constrain the P and S velocities simultaneously.

Scatterer Structure Constrained by Waveform Modeling
Due to the sensitivity of the SdP phase to Vs, we first used forward waveform modeling to put constraints on the shear velocity and thickness of the anomalous structure. We stacked 26 traces with signal-to-noise ratios greater than 15 for direct P waves for which the SxP phase could be clearly identified from the records ( Figure S2). We then calculated theoretical seismograms of the vertical displacement using the 1-D reflectivity method (Wang, 1999). We used S660P converted at the 660-km discontinuity as the reference phase for normalization. The density and P wave velocity do not significantly affect the SxP waveform, as has already been addressed by Haugland et al. (2017). Based on the stacked waveform, it is reasonable to assume a slow Vs anomaly layer at a depth of approximately 900 km (Haugland et al., 2017;Niu, 2014). We varied Vs, allowing the depth to vary between ±10 km, and generated synthetic seismograms. We then compared the observed stacked signal to synthetic waveforms for different δVs and thickness values. Our calculations indicate that an anomaly with a thickness >20 km broadens the pulse width significantly and that the best match has a thickness of~10 km. We then fixed the thickness and varied the shear velocity. The preferred model had a δVs value of approximately −7.2% ( Figure S2d) with the center located at a depth of 895 km ( Figure S1).
We then modeled the X-phase located beneath the orange circle, which was also partly sampled by the teleseismic conversion phases; this allowed both Vp and the density to be constrained simultaneously. This

Geophysical Research Letters
region has the largest number of stacking traces of 361 and thus is the most appropriate candidate for waveform analysis. Because of the temporal normalization and spectral whitening in data processing, the absolute amplitude of the X-phase cannot be directly used (Lin et al., 2011). We used both P410P and P660P as reference phases with the velocity structure around the MTZ taken from the published literature (Li et al., 2013;Wang & Niu, 2010). We fixed δVs to be −7.2% and then changed the Vp, density, and thickness values of the layer to match the observed reflection phases. Figure 2f shows that the synthetic seismogram, the black trace, matches the observations well. Our preferred model has δVp = 0.2% and δρ = 0.6% with a thickness of~20 km. Meanwhile, a relatively sharp 410-km discontinuity with a thickness of~13 km and a broad 660-km discontinuity with a thickness of~76 km are required to match the arrival time and amplitude of the P660P and P410P phases. The structure around the MTZ is consistent with triplication investigations of a thick 660-km discontinuity (Li et al., 2013;Wang & Niu, 2010).
Note that, due to the tradeoff between the density, velocity, and thickness of the velocity anomaly layer, a suite of models is available to fit the X-phase. Based on the Zoeppritz equations, the impedance change is approximately equal to the sum of δVp and δρ when those values are small. A linear change in δVp and δρ will generate a reasonably good fit to the observations ( Figure S3), making it impossible to resolve the exact values of δVp and δρ simultaneously using just one trace. Nevertheless, the P wave tomography image shows that the Vp anomaly at depths of around 1,000 km is approximately 0.2% (Fukao & Obayashi, 2013;Kaneshima, 2019). Based on this inference, we estimate that the density of the anomalous layer is 0.6% larger than the ambient mantle. Because the frequency content of the cross-correlations is predominantly 5-10 s, we cannot resolve layers thinner than 15 km. Therefore, the final preferred model that best fits the X-phase is δVp = 0.2%, δρ = 0.6%, and δVs = −7.2% with a thickness of approximately 20 km centered at a depth of 918 km.

Discussion and Conclusions
Body waves, including the waves reflected at the upper mantle discontinuities and the lower mantle scatterers, are retrieved simultaneously from seismic ambient noise. We notice that the P410P are remarkably similar in all the stacked traces. Because of the use of tapering in the data processing, which decreases the signal smoothly to zero or near zero at ends, we have compared the stacked cross-correlations for cases with and without tapering in Figure S7. The similarity of the stacked waveforms obtained from both processing excludes the potential effect of tapering on the consistency of the observed P410P. We also compared the stacked cross-correlations with those with the interstation distances less than 100 km for region III. In that distance range, the surface coda wave does not affect the P410P phase that arrives later. The results ( Figure S7c) show a coherent feature with that obtained from all the station-pairs within 200 km. The above two systematic comparisons and the bootstrapping uncertainty estimate reconfirm that the P410P is a robust feature in the stacked curves, indicating a relatively simple and consistent 410-km discontinuity in these regions.
The P660P, however, appears at different times and its magnitude varies considerably. We suspect a complicated interaction between the slab and the mantle material around here (Wang et al., 2020), which contributes a lot to the spatial difference of the P660P observed in our study. Previous studies using triplication waveform modeling have found that the 660 discontinuity beneath this study region is depressed and broad with a transitional thickness of 40-50 km (Li et al., 2013;Wang & Niu, 2010). Receiver function analysis shows that the deconvolved Ps converted signals vary substantially and that some show multiple phases (Tian et al., 2016). A radically changed P660P and consistent P410P waves are also observed in published results conducted for the North China Craton,~800 km southwest to our study region (Feng et al., 2017). The amplitude ratio of the P660P may be as high as 2, even in adjacent areas with~40% overlap. In our study, the red region shown in Figure 1a, which does not exhibit the P660P at~160 s, is just located within 50 km to the east of the hinge between the subducting slab and the 660-km discontinuity (Li et al., 2008). The "P660P" sampling in this region appears at an earlier time~134 s, corresponding to phases reflected at depth~570 km, which may be related to the fast velocity anomaly or the upper boundary of the cold slab (Wang et al., 2020) ( Figure S8). The variation of the travel time and amplitude of the P660P and P410P phases will improve the understanding of the geodynamic processes in this region; however, this is beyond the scope of our study and will be left for future investigation.
We identified two anomalies located closely in the lower mantle via simultaneous waveform modeling of two different phases. Due to the limitations of ray sampling, the SxP conversion points were located 30-40 km southeast of the center of the P-to-P reflection points (Figure 4a) with an average depth~23 km shallower than that of the X-phase. We noticed that the P velocity in the crust and upper lithosphere was slow with a δVp value ranging from approximately −1.5% to −2.0% (Fukao & Obayashi, 2013;Liu et al., 2017).
Assuming an average 2% slower velocity layer with a thickness of 100 km, the estimated mid-mantle reflection depth of the P phase should be 4 km shallower than the estimate based on the IASP91 model. Considering the longer period used for the X-phase (~5-10 s) and the~1-5-s period used for the SxP phase, we argue that the two anomalies may share a similar source of heterogeneity or manifest two very spatially close heterogeneities; however, we cannot discriminate the detailed geometry using the available data.
We identified another region with an obvious X-phase in the stacked waveforms (Figure 2f) of the region II marked by the red circle (Figures 1a and 2f). The reflection depth was calculated to be 1,008 km according to the IASP91 model. Niu (2014) investigated the same region in detail via an array analysis of the SdP phases recorded by USArray stations, and the scatterers were estimated to be located at 987 km depth after 3-D migration. This depth discrepancy arises partly from the model we used for the transformation from the time domain to the depth domain. Due to partial area-overlapping of the scatterers detected by Niu (2014) with the reflection region of the ambient noise samples, we propose that both may originate from a similar scattering source at a mid-mantle depth range of 980-1,000 km. The detailed geometric relationship between these scatterers needs to be further investigated.
Fine-scale scatterers in the mid-lower mantle have been identified beneath subduction zones, with most located in the depth range of 800-1,800 km (e.g., Haugland et al., 2017;Kaneshima, 2016Kaneshima, , 2019Li & Yuen, 2014;Niu, 2014). In such studies, the S wave velocity varies significantly between −2% and −12% (Haugland et al., 2017;Kaneshima, 2016), and the density anomaly is generally proposed to be higher than that found in this study, with the largest values of 2-9% detected in the Mariana subduction zone (Niu et al., 2003). Due to the insensitivity of the SdP phase to Vp, studies on the P velocity are limited with variations between −1% and 1% (Kaneshima, 2019). Our results indicate that small-scale heterogeneities in the mid-mantle can be manifested in both converted and reflected waves. Simultaneous forward waveform  Fukao and Obayashi (2013). The symbols are the same as those used in Figure 1a. Profile AB was selected to pass through the center of the circular regions II and III. (b) Cross section along the profile AB. The solid and dashed white lines denote branches of the S and P waves, respectively. The orange and red rectangles denote the locations of the scatterers revealed by the ambient noise; the blue rectangle represents the scatterer detected by the S-to-P phase analyzed in our study, while the gray rectangle represents that reported by Niu (2014). modeling indicates that the detected heterogeneities are located in a depth range of 900-1,000 km beneath the area west of the Japan Sea, with H = 10-20 km, δVs = −7.2%, and δρ = 0.6%.
Considering an anomaly with a large Vs but smaller Vp variation, we argue that the small-scale heterogeneities that produce the observed high-amplitude SP and PP phases are associated with the mid-ocean ridge basalt (MORB) subducted into the lower mantle. At depths below 800 km, the transition in silica from ferroelastic post-stishovite to the CaCl 2 -type phase (Tsuchiya, 2011) will give rise to substantial differences in elasticity, which may be observed in the converted and reflected phases in our study. The scatterers are distributed at different depths and show variable strengths, which may indicate the complex composition of the mantle material or the different degrees of the phase transitions, which vary with the composition, temperature, and compression conditions. We propose that the detected scatterers may originate from a compositionally differentiated reservoir related to the ancient subducted basaltic crust of the Izanagi plate (Li & Yuen, 2014), which subducted completely beneath Southern Japan 50-60 Ma ago (Seton et al., 2012).
In summary, we retrieved body waves reflected from scatterers in the mid-lower mantle for the first time using seismic noise correlations. The phase reflected from scatterers at depths of~900-1,000 km that appeared in the noise correlation curves may be partly due to our stacking scheme and to the relatively large number of correlation curves. We explored and confirmed the PP-waves reflected from scatterers in the mid-lower mantle using the teleseismic SdP phases. We obtained a velocity anomaly generated by the reflected and converted signals simultaneously. Our study illustrates the exciting potential to explore the fine-scale structure in the mid-lower mantle, especially in regions far from subduction zones, which can be used to constrain the mineralogy and dynamic processes of the deep Earth.

Data Availability Statement
Seismic data from the CEA stations are provided by Data Management Center of China National Seismic at Institute of Geophysics, China Earthquake Administration (https://doi.org/10.11998/SeisDmc/SN, http:// www.seisdmc.ac.cn) and the cross-correlation functions can be downloaded via https://doi.org/10.6084/ m9.figshare.12708362. Earthquake waveforms of the stations in Europe are downloaded from the GFZ Metadata Editor (http://eida.gfz-potsdam.de/webdc3/). The Northeast China Extended SeiSmic Array (NECESSArray) data were downloaded through the Incorporated Research Institutions for Seismology (https://doi.org/10.7914/SN/YP_2009). Waveforms of the NorthEast China Seismic Array to Investigate Deep Subduction (NECsaids) data are deposited in the Seismic Array Laboratory, Institute of Geology and Geophysics, Chinese Academy of Sciences (https://doi.org/10.12129/IGGSL.Data.Observation, http:// www.seislab.cn) and the cross-correlation functions of seismic ambient noise obtained in this study can be downloaded via https://doi.org/10.6084/m9.figshare.12708362. The Seismic Array Laboratory will make the NECsaids Array data publicly available from October 2021 (3 years after the completion of the NECsaids project). Sac2000 and GMT (Wessel et al., 2013) are used in basic data processing and plotting of some figures.