The teleseismic data provided in this study obtained from the IRIS global network, these data were recorded at a total of 28 broadband OBSs deployed on Tamu Massif of Shatsky Rise by the Woods Hole Oceanography Institution (WHOI). Our data consists of the teleseismic events recorded by the hydrophone component of each station between July and September 2010, all the earthquake events used in the study are based on the earthquake catalog published by the United States Geological Survey (USGS), and these event were selected with the following criteria: (a) earthquakes with M_w ≥ 4.8; (b) epicenter distance between 30° and 90°. We excluded earthquakes with epicenter distances less than 30° in order to keep the near-vertical incidence of P-wave. Earthquakes with epicenter distances between 90° and 120° were excluded due to the ray paths of these events bypass at the core-mantle boundary and being sensitive to the core. After selected, we have obtained a data set with raw hydrophone component records of a total of 155 seismic events (Figure 2; Supplementary Table S1).

The hydrophone component of the OBS was chosen for this study for two reasons. First, if the OBSs are poorly coupled to the substrate when deployed on the soft seafloor sediments, resulting in low signal-to-noise ratios of the recorded three-component signals, particularly the horizontal component recordings, makes it difficult to apply the receiver function which uses passive source data, although it has been widely used in the detection of continental lithospheric discontinuities. While the hydrophone component of OBS is suspended in water and is not affected by the coupling condition, the signal-to-noise ratio of the recorded signal thus is higher than the three-component signals and provides more useful information. Second, compared with the traditional component, uses the hydrophone component alone for autocorrelation calculation, reducing the step of rotating the component in the conventional seismic method, it alleviates the difficulty of usage and simplifies the processing steps. In this study, a total of 28 broadband OBSs were deployed along two survey lines across Tamu Massif (Figure 1). Line A extends through the entire Tamu Massif from west to east, while line B extends from south to north, covering the central part of the volcano. The distribution of these two lines allows us to obtain the structural changes information of the oceanic crust in the different directions of Tamu Massif. The thin sediment beneath most of the stations has less influence on the arrival of the reflection signal, ensuring the accuracy of the results.

In this study, we apply a novel geophysical method proposed by Dong et al. (2022) to calculate the autocorrelation of the teleseismic P-wave coda recorded in the hydrophone component of OBSs. P-wave coda is the oscillations of the ground continuing long after the slowest theoretical primary wave arrival time, this part corresponds to the recorder energy after the passage of all primary waves. Since the previous studies, it has been clearly proven that coda waves are an effective resource in extracting information about the characteristics of wave paths, it contains a great deal of information from underground (Herraiz and Espinosa, 1987). Based on the characteristics of the P-wave coda, we can calculate its autocorrelation to detect the lithospheric discontinuities.

After obtaining the original waveform records, we convert the original format to SAC format, and rewrite the channel header information. Then, we remove some non-zero values and linear trends of the data, which are caused by the nature of the instrument, in order to preserve the calculation and analysis of the data. In the following step, we select a window around the theoretical primary wave arrival time in seismogram, the window length is a generic parameter adapted to a particular type of application (Phạm and Tkalčić, 2017). We used the AK135 global velocity model (Kennett et al., 1995) to calculate the theoretical arrival times of P-waves for different seismic events, and intercepted waveforms 5 s before and 45 s after the theoretical arrival time of the selected seismic events, this window length containing the P-wave coda data required for our study, the multiple reflection signals from the lithospheric discontinuity surface and reflections from the deeper lithosphere. Meanwhile, we use Hanning windows to taper off both ends of the waveform to avoid spectral domain artifacts during the time domain to frequency domain conversion and calculation.

Due to Earth filtering, high-frequency signals are quickly absorbed and attenuated by geodesic filtering and low-frequency signals become dominant. To address the problem of low-frequency bias in the records, we adopted spectral whitening to amplify the high-frequency signal content while reducing the low-frequency signal content. The variation of the width of the spectral whitening window we set affects the distribution and energy of the illusion only and does not interfere with the signal from the reflecting interface, thus not affecting our identification of the real signal. Finally, we use a zero-phase band-pass filter of 0.5–3 Hz to improve the sharpness of the reflection signal (by removing the very long period signals) and to avoid spurious effects caused by the unexpected amplification of very high frequency noise due to the spectral whitening.

We compute the autocorrelation of the whitened waveform in the time domain, and retain a half of the symmetric autocorrelogram having positive time lags. After calculating the P-wave coda autocorrelation, we normalize autocorrelation results according to the amplitude and then stack the normalized results by linear stacking to improve the signal-to-noise ratio. Then the stacked P-wave coda autocorrelation results of each station are arranged according to projection points to obtain the P-wave coda autocorrelation amplitude curve profile of Tamu Massif. The oceanic crust structure and Moho reflectors beneath the volcano are then identified through the reflected signals on the profile.

For the results with only two signal peaks, we first exclude the reflected signals from the water surface, then identify the reflected signals from the oceanic crust and characterize these signals to depict an overview of the structure of the oceanic crust beneath the Tamu massif. Taking station A16 as an example (Figure 4), this station is located on the flanks slope of Tamu Massif, where the water depth is quite deep and the wave propagation time in the water is relatively long. We arrange the autocorrelation results for different seismic events according to the ray parameter, which indicates the unique horizontal slowness corresponding to each earthquake event, then identify structure information by recognizing reflected signals on the ray parameters. In the autocorrelation results for station A16, the signal reflected from the water surface is clearly separated from the signal reflected, which from the oceanic crust by stacking the autocorrelation results for each ray parameter, thus two signals can be easily distinguished from each other. The water depth of station A01 is 3,755 m and the wave propagation speed in water is 1.5 km/s, thus the theoretical two-way travel time of the reflected signal from the water surface should be at about 5.00 s. The polarity should be negative after the reflection from the water surface (Dong et al., 2022). In Figure 4, we can see an obvious negative phase with a two-way travel time of about 5.0 s, which is consistent with the theoretical characteristics of the water surface reflection signal, so we identify it as the reflection signal from the water surface. In analogy, the two-way travel time of the water surface reflections can be calculated based on the deployment depth of each OBS in the cruise record of Korenaga and sager (2010). Combining with the polarity of phases, we can distinguish two signals that are close in time, thus other reflection signals could be analyzed more accurately. In Figure 4, it can be seen that after linear stacking of the results, another reflection signal is present at ∼2.0 s, which indicates the presence of a reflection interface at this depth, and this interface is observed in several ray parameters.

Different with Figures 5, 6, the two-way travel times of the reflection signals from the internal structure of oceanic crust is not consistent in different ray parameters generally distribute in 1.8∼2.2 s (Figure 4). After stacking, it can be identified that the depth of the structure reflector is located at the two-way travel time of ∼1.9 s beneath the Tamu massif. Autocorrelation results for most of the OBS stations at the Tamu Massif show a similar pattern, that is the reflected signals recorded on different ray parameters are scattered. On the contrary, station A01 (Figure 5) and A07 (Figure 6), demonstrate that the times of reflection signals from internal structures in the oceanic crust are relatively consistent at most ray parameter. The reflection signals from this structure are seen at the two-way travel time of ∼2.4 s and can be identified without stacking. In the profile of line A, reflected signals from internal oceanic crust structure can be seen in the stacked autocorrelation results of most stations. These signals vary in depth, they are all located at 1.5–2.5 s (Figure 7) and are generally parallel to the terrain.

We consider this signal as reflections from a layer composed of alternating eruptions of denser, higher velocity massive lava and sparser, lower velocity pillow lava flows, which have less density and lower velocity compared to the lower oceanic crust, with a strong acoustic impedance contrasts between them and thus able to generate a reflection signal.

On the profile, except for the reflection signals of the layer composed of lava flows and the water surface, some stations have also recorded a signal from other reflectors. Taking stacked autocorrelation results of station A01 as an example (Figure 5), there is a reflection signal at 3.6 s that has not yet been recognized, in addition to the reflection signals from the water surface at 7.2 s and the lava flow at 2.3 s. This reflection signal has the same characteristics as the two identified reflections in the autocorrelation results, with a clear amplitude at the two-way travel time. In addition, the signal keeps clear when adapting different widths of the Hanning window or bandpass range, indicating an interface at the corresponding depth. After distinguishing the reflection signals received by these stations from this interface and projecting them on the line A profile of Tamu Massif (Figure 7), it can be seen that the interface is situated a few seconds beneath the Tamu Massif on both sides and deepen overall as it approaches the middle. Signals from these reflectors are seen at about 3.6 s below station A01, at about 4.5 s below station A05, at about 3.7 s for A21 (Figure 7), while these reflectors are not observed at the remaining stations. This suggests that this interface is not quite continuous or the interface is too weak to produce a consistent reflective signal.

For those signals with two-way travel times >2.5 s, we consider them to be reflections from the Moho. Similar to the receiver function, probably because the energy of the secondary Moho is faint and weaker than the noise energy, the reflection signal of 2p_mp is invisible in results, thus we focus on identifying the signals from the primary reflection phase of the Moho. The results are agreement with previous results obtained from various geophysical methods (Sager et al., 2010; Korenaga and Sager, 2012; Zhang et al., 2015; Zhang et al., 2016; Chen et al., 2022), that is, the Moho reflectors move from shallow to deep from the sides toward the middle, the thickness of the oceanic crust becomes thicker (Figure 7). Because of the thick crust in the middle of the volcano, it is difficult to receive the reflected signals from the Moho reflectors in the middle part.

Before recognizing the received signal’s the source, we took steps to confirm the signal’s reliability, ensuring that the signal we analyze in the discussion following is a true reflected signal. To ensure that our target signal is preserved as much as possible and to minimize interference from other frequency signals, filtering is a crucial step in this work. Our expected detection range is the lithosphere, and at the beginning of the study, because the bandpass filter was set too low, the reflection signals from deeper depths were mixed with the expected signals, making it impossible to distinguish and identify the signals useful for this study. When the bandpass filter was set too high, there were too many high frequency signals that contain information that is not useful to us, and we could not identify which was the expected signal. After extensive testing, we set the bandpass filtering in the range of 0.5∼3.0 hz. Within that bandwidth, the reflected signal from the oceanic crustal structure beneath Tamu massif is most prominent. We also adjusted the window of spectral whitening and the length of the cut time window of the received signal, and can see that the waveform of the real reflected signal does not change after doing these operations, so we can confirm which are the valid signals in the autocorrelation results.

We propose that the reflection signal 1.5–2.5 s below the sea floor of Tamu Massif is from a layer composed of alternating eruptions of denser, higher velocity massive lava and sparser, lower velocity pillow lava flows, which have less density and lower velocity compared to the lower oceanic crust, judged mainly in combination with previous Integrated Ocean Drilling Program (IODP) core and borehole logging data (Sager et al., 2011), MCS profiles interpreted by synthetic seismogram for IODP (Sager et al., 2013; Zhang et al., 2015; Zhang et al., 2016), OBS refraction signals (Korenaga and Sager, 2012). Evidence from synthetic seismograms of core and borehole logging data from ocean drilling suggests that the basement of Tamu Massif is dominated by dense massive lava flows with pillow lava flows interspersed (Sager et al., 2013), with similar features beneath the basement of some other oceanic plateaus such as the Ontong Java Plateau, the Manihiki Plateau and Kergulen Plateau (Rotsein, 1992; Rotstein et al., 1992; Inoue et al., 2008; Spitzer et al., 2008; Pietsch and Uenzelmann, 2015). Results from IODP and synthetic seismograms show that massive lava flows interstratify with pillow lava flows above the oceanic crust of the Tamu Massif, supporting the possibility that the upper oceanic crust portion of the Tamu Massif can constitute a mixed layer of massive and pillow lava flows. The velocity structure provided by the OBS refraction also demonstrates that the seismic wave velocity above the oceanic crust is low and increases with depth, and eventually has a significant difference with the seismic wave velocity in the middle of the oceanic crust, and we find the interface of the seismic wave velocity discontinuity that may be caused by this layer through the teleseismic P-wave coda autocorrelation. The depth of this layer can be seen at about 1.5–2.5 s of two-way travel time beneath the sea floor. They are generally parallel to the seafloor terrain, may suggesting that they may formed during rapid eruptions of viscous low magma.

After identifying the reflection signals from the water surface and lava flows in the autocorrelation results, reflection signals from other interfaces deeper below the sea floor can be identified consequently. Reflection signals indicate a discontinuous interface distributed on both sides of Tamu Massif and the interface signal gradually disappears as moves toward the middle of the volcano. The energy and the quality of the interface reflected signal varies in stations, with some stations receiving clear signals and some slightly blurred, for which we infer that the signal comes from the Moho reflectors at the crust-mantle boundary. Although there is no independent verification that this reflector is Moho, but we can comparing our results with the oceanic crust models of previous multi-channel seismic surveys and gravity investigation (Korenaga and Sager, 2012; Zhang et al., 2015; Zhang et al., 2016; Chen et al., 2022). We found some stations that the two-way travel time of the signal is consistent with Moho reflectors detected by others, and Moho distribution pattern consistent with the previous observations that oceanic crust grows thicker as it approaches the middle of the volcano (Figure 8). In Figure 8, we project some OBS stations which have observed the Moho to the corresponding location of Zhang et al. (2016) MCS reflection profile and compare with their results, the red dot represents the reflection signal we received from Moho, the projected signal two-way travel time is consistent with the Moho as interpreted by Zhang et al. (2016) based on the multi-channel seismic reflection profile.

For A07 station, except for the reflection signals from water and lava flows, a signal is shown at 4.2 s (Figure 6). We do not interpret it as Moho because its arrival is shallower than that of A05, which is inconsistent with the gradual deepening trend of Moho as interpreted by previous researchers (Korenaga and Sager, 2012; Zhang et al., 2015; Zhang et al., 2016; Chen et al., 2022). Hence, for the signal that is shallower than Moho, we infer that it may be the transition from the upper crust to the lower crust, or a reflection caused by the intrusion within the lower crust.

Oceanic plateaus are less studied than continental since they are more difficult to access and submerged under the water, making it difficult for researchers to collect data. Even though it might be challenging to collect data, it is still essential to investigate the oceanic plateau. Scholars have attempted to research the oceanic plateaus from a variety of approaches in order to understand and learn more about the evolutionary formation mechanisms of these plateaus, in Shatsky Rise, there have been some geophysical studies here (Sager et al., 2010; Korenaga and Sager, 2012; Zhang et al., 2015; Zhang et al., 2016; Sager et al., 2019; Huang et al., 2020; Chen et al., 2022). These researches have collected spatial information and stratigraphic structure of the Shatsky Rise region by different geophysical methods. In this paper, teleseismic P-wave coda autocorrelation was applied in this region, and these results provide new constraints on the structure of the oceanic crust beneath the Tamu Massif.

In our autocorrelation results, we have retrieved the structure of the oceanic crust of the Tamu Massif. Within the volcano, a layer composed of denser, higher velocity massive lava and sparser, lower velocity pillow lava flows is found, although this signal has not been identified in active source seismic data (Sager et al., 2010; Zhang et al., 2015; Zhang et al., 2016). We infer that the cause of this situation is the difference in resolution among geophysical methods. Active source seismic using to shorter wave length, the obtained seismic profiles usually have a higher resolution and can obtain the delicate structure of the underground. In contrast, the P-wave of passive sources has a longer wavelength and can reflect a larger range of velocity variations, although at lower resolution. Therefore, we infer that the strata of Massive lava and pillow lava flows in the Tamu Massif are deeper beyond the resolution of MCS profiles, but can be detected by our approach.

Although the two-way travel time of the reflected signals from the water surface shows strong consistency, the two-way travel time of the signals from the internal structure of the oceanic crust shows some variations in some of the autocorrelation results with different ray parameters (Figure 4), and these signals can recover clear reflections after stacking. But in some other stations, there is a high consistency among the autocorrelation results of the signals from the internal structure of the oceanic crust in different ray parameters, such as A07 (Figure 7), those strong and clear lava flow reflection signals can be identified in most of the autocorrelation results of seismic events without linear stacking. We deduce the following two explanations as to why some stations failed to record the interface of alternating lava flows: 1) The eruption process is not stable and homogeneous, the characteristics of the interface thus are not consistent, resulting in some parts of the interface being clearer and some parts being blurred. 2) Due to the short period of time of instrument deployment, the data used in this study might not with sufficient quality to show the reflections.