Extraction of Mantle Discontinuities From Teleseismic Body‐Wave Microseisms

Ocean swell activities excite body‐wave microseisms that contain information on the Earth's internal structure. Although seismic interferometry is feasible for exploring structures, it faces the problem of spurious phases stemming from an inhomogeneous source distribution. This paper proposes a new method for inferring seismic discontinuity structures beneath receivers using body‐wave microseisms. This method considers the excitation sources of body‐wave microseisms to be spatially localized and persistent over time. To detect the P‐s conversion beneath the receivers, we generalize the receiver function analysis for earthquakes to body‐wave microseisms. The resultant receiver functions are migrated to the depth section. The detected 410‐ and 660‐km mantle discontinuities are consistent with the results obtained using earthquakes, thereby demonstrating the feasibility of our method for exploring deep‐earth interiors. This study is a significant step toward body‐wave exploration considering the sources of P‐wave microseisms to be isolated events.


Introduction
Microseisms are random seismic wavefields with a frequency range of 0.05-0.50Hz excited by ocean swell activities (e.g., Nishida, 2017) and can be categorized into primary microseisms (PMs; <0.1 Hz) and secondary microseisms (SMs; >0.1 Hz).PMs are excited by topographic coupling between surface ocean gravity waves and seismic waves (Hasselmann, 1963), whereas SMs are excited by the nonlinear effects of surface ocean gravity waves (Longuet-Higgins, 1950).Although Rayleigh and Love waves dominate the SMs, teleseismic P-waves have also been observed (e.g., Gerstoft et al., 2008).It has long been understood that these random wavefields represent noise in earthquake seismology.
In the late 2000s, seismic interferometry (SI) turned the noise into a signal elucidating the Earth's internal structures (e.g., Snieder & Larose, 2013).SI is a technique for extracting wave propagations between seismic stations by calculating the cross-correlation of random seismic wavefields of station pairs.Because the random wavefield is excited by distributed noise sources, only noise sources within the stationary phase regions contribute constructively to the cross-correlation function.Because this assumption is more valid for surface waves, surface-wave exploration became widely used in the late 2000s (e.g., Shapiro & Campillo, 2004).

10.1029/2023GL105017
2 of 8 prohibiting the accurate reconstruction of wave propagations (e.g., Li et al., 2020;Pedersen & Colombi, 2018).Such studies showed that for estimating seismic structures using P-wave microseisms, considering their source locations is important.Back projection before cross-correlation is one of the possible solutions for improved extraction (e.g., Liu & Shearer, 2022).
This study employs a strategy that is different from SI. Instead of assuming distributed sources, P-wave microseisms are assumed to be spatially isolated events (within several hundred kilometers) with persistent excitation (∼6 hr).We developed a new method for estimating seismic structures using P waves from the centroid locations of P-wave microseisms, focusing on the P-s conversion beneath the receivers.To detect this conversion, we generalized a receiver function method for earthquake data (e.g., Langston, 1979) to P-wave microseisms.Receiver function analysis using earthquakes sometimes involves narrow azimuthal coverage from inhomogeneous hypocenters.Utilizing P-wave microseisms in the generalized receiver function (gRF) analysis can improve the azimuthal coverage and broaden the illuminated area.
This study targets mantle discontinuities at 410/660 km depth.The structure of mantle discontinuities has been inferred by several methods (Kind & Li, 2015), such as the precursor of PP and SS phase (e.g., Shearer & Masters, 1992), ScS reverberation (e.g., Revenaugh & Jordan, 1991), and receiver functions (e.g., Tonegawa et al., 2005).We will discuss the feasibility of our method by comparing them with previous studies.

Data
We used seismograms from 690 seismic stations (Figure 1b) of the high-sensitivity seismograph network (Hi-net: https://doi.org/10.17598/NIED.0003;Okada et al., 2004) deployed by the National Research Institute for Earth Science and Disaster Resilience (NIED).The vertical and horizontal components of the velocity meters with a natural frequency of 1 Hz were deployed at the bottom of the boreholes of the stations.We analyzed the data after eliminating the coherent periodic noise originating from the logger (Takagi et al., 2015) and correcting for the instrumental response in the time domain (Maeda et al., 2011).
To detect the P-s conversion, we analyzed teleseismic P-wave microseisms from 0.10 to 0.25 Hz.This frequency range was determined by the narrow spectral content of P-wave microseisms (e.g., Figure 3 of Nishida & Takagi, 2016) due to the acoustic resonance of the ocean (Gualtieri et al., 2014).Assuming that the excitation  (Nishida & Takagi, 2022).The contour lines show the equidistant circles at every 10° from the N.HNSH station (34.7283°N 134.2744°E,Okayama Prefecture, Japan) of Hi-net.(b) Blue triangles denote the Hi-net stations used.(c) The ray paths of P410s (blue line) and P660s (green line) at the epicentral distance of 60° with AK135 (Kennett et al., 1995).The star denotes the P-wave source location, whereas the red triangle denotes the station.

10.1029/2023GL105017
3 of 8 sources were persistent in time but localized in space, they were approximated using a persistent vertical single force at the centroid location.Nishida and Takagi (2022) constructed a centroid single-force catalog of P-wave microseisms from 2004 to 2020 using Hi-net data.To construct the catalog, they developed an autofocusing method that utilizes information on both the slowness and wavefront curvature.Although the catalog also includes PP and PKIKP events, the number of events with high signal-to-noise ratios is smaller than that of P events.Accordingly, we mainly focused on P events.Each event in the catalog includes the event date, centroid location, centroid single force, P-wave beam power, and median absolute deviation of the beamforming result.We chose 5,780 P events with high signal-to-noise ratios.The signal-to-noise ratio was defined as the ratio of the P-wave beam power to the median absolute deviation of the beamforming result.The centroid locations of the P-wave microseisms were classified as the Northern Atlantic Ocean, Northern Pacific Ocean, and Southern Pacific Ocean (Figure 1).The histogram of epicentral distances revealed peaks at 30-40° and 90-100° (Figure S1a in Supporting Information S1).

Calculation of gRFs of P-Wave Microseisms
The gRF analysis for an earthquake was generalized for persistent P-wave microseisms.The source time function of P-wave microseisms is typically several hours to several days (e.g., Zhang et al., 2010), whereas those of most earthquakes are shorter than several minutes.This section describes the procedure for generating gRFs using P-wave microseisms.The centroid location of the P-wave microseisms was fixed at 6 hr.Six-hour-long waveforms were split into 1,024-s time windows.We selected time windows with low mean squared amplitudes of 0.05-0.10Hz (i.e., less than 7.6 × 10 4 nm 2 /s 2 ) and 0.10-0.20 Hz (i.e., less than 2.5 × 10 4 nm 2 /s 2 ) to avoid contamination from local earthquakes and local Rayleigh waves.In this frequency range, locally excited Rayleigh waves dominate the records in Japan because teleseismic Rayleigh waves attenuate with propagation owing to scattering and intrinsic attenuation.Indeed, the relative amplitude of P waves to Rayleigh waves is anti-correlated with the mean squared amplitude in the 0.1-0.2Hz range in Japan (Takagi et al., 2018).Although data selection according to the H/V ratio at inland stations is feasible for enhancing teleseismic P-wave microseisms (e.g., Pedersen et al., 2023), we chose data with a low mean squared amplitude.
For each time step, the Fourier spectrum of the source time function of the incident P-wave P(ω) was estimated by stacking vertical components all over the stations, along with the theoretical P-wave travel time.Stacking the vertical components of all stations enhanced the relative amplitude of body waves to Rayleigh waves.The theoretical P-wave travel time T ip was calculated using the tauP toolkit (Crotwell et al., 1999) for the 1-D velocity model AK135 (Kennett et al., 1995), with correction for the 3-D velocity model (0-60 km depth) beneath Japan (Nishida et al., 2008), assuming a vertical incident.The Fourier spectrum of the source time function of the jth time step is given by where N is the number of seismic stations and, Z ij is the vertical component of the jth time window at the ith station.The gRF of P-wave microseisms RF i at frequency ω is evaluated by minimizing the squared difference S i (ω) between the radial component R ij (the ith station and jth time window) and RF i ⋅ P j as, where Nt is the number of time steps.At each frequency, from the condition that the partial derivatives of S i with respect to the receiver function are zero, the radial gRF of P-wave microseisms RF i can be calculated as where w is the 5% water level to avoid instability.The numerator is the power spectrum of the incident P-wave   * , whereas the denominator is the cross-spectrum of the incident P-wave P j and the radial component R ij .In the case of only a one-time window, this equation is equivalent to the typical receiver function, including the water level.The vertical component of the gRF is calculated in the same manner, replacing R ij (ω) by Z ij (ω) in Equation 3.
All the radial gRFs for all events were linearly binned-stacked to enhance the signal-to-noise ratio.The width of the distance bin was 25 km.All radial gRFs RF i were normalized by the peak amplitude of the vertical gRFs ZF i , and aligned at the peak time of ZF i (e.g., Zhang et al., 2010).Quality control of the gRFs was performed before stacking, rejecting those with large amplitudes of 50-500 s before the P peak.The radial gRFs using PP events were calculated in the same manner, while those using PKIKP events were calculated in a slightly different manner (Supporting Information S1).Because the number of events was less than the number of P events, results are shown in Supporting Information S1 (Figures S3-S6).

Binned Stack of gRFs
The binned stack of radial gRFs shows the P-s converted waves at the mantle discontinuities with an epicentral distance of 3,000-11,000 km (Figure 2).The relative arrival times of these waves were consistent with the theoretical travel times of AK135.
The ringing shape of P waves varies with the epicentral distance because the narrow-band spectral content of SMs is affected by the ocean depth in the source area because of the resonance of the water layer (e.g., Gualtieri et al., 2014).The amplitude variation of vertical gRFs before the P peaks depends on the stacking number (Figure S2 in Supporting Information S1).The binned stack of the vertical gRFs shows no clear PP-and P-reflected waves at the mantle discontinuities (PP410P/PP660P).
Possible reasons for the emergence of P-s converted phases in the radial components and the lack of these reflection phases in the vertical components are as follows: First, the reflection phases being included in the source time function because of a small difference between the PP410P and PP660P slowness and that of P (∼0.0005 s/ km with an epicentral distance of 80°); and second, the reverberation at the bounce point causing destructive interference for the reflected phases.However, the P-s converted phases were free of these reasons.

1-D Depth Migration
All the radial gRFs were depth-migrated with the 1-D structure of AK135.Assuming that the radial gRFs consist of P-s converted waves, all gRFs with epicentral distances greater than 30° were stacked along the theoretical travel time of the P-s converted waves at depths of 200-1,000 km.Travel time was calculated using the tauP toolkit (Crotwell et al., 1999).Four groups (all events, Northern Atlantic events, Northern Pacific events, and Southern Pacific events) of gRFs were depth-migrated (Figure 3).At a shallow depth (200 km), results were 10.1029/2023GL105017 5 of 8 characterized by large amplitudes, which can be explained by contamination by the ringing of the P waves.All results showed common positive peaks around the depths of 410 and 660 km.The depth of peak amplitude was somewhat deep because of the difference in the crustal structure between Japan and AK135.
The amplitudes of the P-s waves were compared to previous seismic observations.The P-s transmission coefficients estimated using Japanese broadband stations are 2.7 ± 0.8% and 6.0 ± 0.7% (Kato & Kawakatsu, 2001), while using the PREM parameters (Dziewonski & Anderson, 1981), they are 1.5% and 4.6% (assuming the near vertical incident with slowness 0.06 s/km, based on Equation 4. 3.34 of Saito, 2016).In this study, we obtained the P-s transmission coefficients of the 410-and 660-km discontinuities as 0.8%-1.8%and 0.7%-1.5%,respectively.These values are smaller than those reported in previous studies, with the P-s transmission coefficient at 660 km depth being smaller than that at 410 km depth.Waveform distortion due to the high water level during deconvolution may decrease the amplitude of P-s converted waves.However, waveform distortion alone cannot explain the smaller amplitude of P660s compared with that of P410s.A plausible reason is the sharpness of the 410-and 660-km discontinuities, that is, the frequency dependence of the P-s transmission coefficients because this study used a relatively higher frequency band (0.10-0.25 Hz) than those in previous studies (0.005-0.200Hz).The widths of the mantle discontinuities may be related to the water content of the transition zone (e.g., Helffrich & Wood, 1996).Our results are consistent with the observed frequency dependence of the P-s conversion at the 660 km discontinuity reported previously (Tonegawa et al., 2005).Another possible reason is that the depression of the 660-km discontinuity (a ∼40 km depression in southwest Japan; Tono et al., 2005;Tonegawa et al., 2006) may make P660s incoherent for stacking.

3-D Depth Migration
The migration of all the radial gRFs with the 3-D structure was calculated.Assuming that the radial gRFs consisted of P-s converted waves, all the gRFs with epicentral distances greater than 30° were projected to the 100-1,000 km depth domain.Figure 4 shows a cross-section using these depth-projected gRFs.At each point on the cross-section, the gRFs were stacked with a horizontal distance smaller than 1°.The image of the cross-section shows the horizontally continuous 410-and 660-km discontinuities.The 410-km discontinuity has a gap where the Pacific plate (red line in the figure) crosses the 410 km depth.The 660-km discontinuity is depressed where the Pacific slab stagnates above the 660-km discontinuity.However, the Pacific plate is not evident in the cross-section because the spatial stacking length is too large to image the subducting plate surface (Kawakatsu & Yoshioka, 2011).10.1029/2023GL105017 6 of 8 The cross-section from the resultant gRFs was compared with the stacked image from conventional receiver functions using an earthquake (Tonegawa et al., 2006).Because the locations of both the cross-section and frequency band are similar, the cross-section of Tonegawa et al. ( 2006) is adequate for reference.The signatures of the mantle discontinuities were common in both images, demonstrating the reliability of this method.

Applicability of This Method to Other Seismic Arrays
The relative amplitude of body waves to surface waves of microseisms is one of the factors determining the signal-to-noise ratio in this method.The waveforms of the Hi-net stations are dominated by surface waves excited by the sea around Japan, whereas those of stations in the continental region are less dominated by surface waves owing to scattering decay (e.g., Vinnik, 1973).An inland array is preferable; however, our method still requires more than 500 stations to accurately extract the incident P-wave microseism.It also requires long-term observations, with time periods of over 10 years, and a sufficient number of events.Because dense observations have become popular in the last decade, our method could be feasible for extracting P-s converted waves using seismic networks on the continent.

Conclusions
We performed a generalized receiver function analysis of P-wave microseisms.Contrary to the SI assumption, this method assumes that the sources of body-wave microseisms are persistent in time and spatially localized.This strategy constitutes a solution to the spurious phase problem of body-wave extraction using SI.The resultant migration image showed the extraction of the P410s/P660s from the radial gRFs.The stacked image from the extracted P-s converted waves was consistent with the images from the conventional receiver function analysis using earthquakes.The results demonstrate the feasibility of our method using body-wave microseisms for exploring the Earth's deep interior, while this method can also be applied to other modern arrays.This study is a significant step toward seismic exploration while considering the sources of P-wave microseisms to be isolated events; it has the potential to extract information not only from the receiver side but also along the path, including the source side.

Figure 1 .
Figure 1.Centroid distribution of the mantle P-wave microseisms used in this study.(a) Red dots denote 5,780 P-wave source locations (Nishida & Takagi, 2022).The contour lines show the equidistant circles at every 10° from the N.HNSH station (34.7283°N 134.2744°E,Okayama Prefecture, Japan) of Hi-net.(b) Blue triangles denote the Hi-net stations used.(c)The ray paths of P410s (blue line) and P660s (green line) at the epicentral distance of 60° with AK135(Kennett et al., 1995).The star denotes the P-wave source location, whereas the red triangle denotes the station.

Figure 2 .
Figure 2. Binned stack of all the radial gRFs of all seismic stations.Each stacked trace is normalized by the maximum value.The peak at 0 s is the P-wave.Two dashed lines show the theoretical travel times of P410s and P660s from AK135.The figure shows these phases along each line.The right panel shows the number of radial gRFs used in the binned stack.

Figure 3 .
Figure 3. (a) Depth migration result of the radial gRFs converted from 200 to 1,000 km depth.Red, blue, green, and black lines indicate the results using all the events, events in the Northern Atlantic, events in the Northern Pacific, and events in the Southern Pacific, respectively.Peaks exist around the depths of 410 and 660 km.(b) Enlarged view in the 370-450 km depth range in panel (a).(c) Enlarged view in the 600-720 km depth range in panel (a).

Figure 4 .
Figure 4. Cross-section (A-A′) of the depth-projected radial gRFs.(a) Locations of the cross-section and stations (blue triangles) used in this study.(b) Stacked image of the A-A′ cross-section from 100 to 1,000 km depth.The red line denotes the depth of the top of Pacific slabs (Nakajima & Hasegawa, 2006; Nakajima et al., 2009).The shallow part near A′ is muted due to the poor ray path coverage.(c) Locations of the cross-section and stations (black triangles) of Tonegawa et al. (2006).(d) Stacked images from receiver function analysis of earthquakes (Tonegawa et al., 2006).This result used a 0.16 Hz low-pass Gaussian filter.OMH refers to the oceanic Moho, and LBP shows the lower boundary of the Pacific slab.Panels (c) and (d) were modified from Figures 1a and 6b of Tonegawa et al. (2006).