Near-Field Body-Wave Extraction From Ambient Seafloor Noise in the Nankai Subduction Zone

Ambient noise correlation is capable of retrieving waves propagating between two receivers. Although waves retrieved using this technique are primarily surface waves, the retrieval of body waves, including direct, refracted, and reflected waves, has also been reported from land-based observations. The difficulty of body wave extraction may be caused by large amplitudes and little attenuation of surface waves excited by microseisms, indicating that body wave extraction using seafloor records is very challenging because microseisms are generated in ocean areas and large amplitudes of surface waves are presumably observed at the seafloor. In this study, we used a unique dataset acquired by dense arrays deployed in the Nankai subduction zone, including a permanent cabled-network of 49 stations, a borehole sensor, and 150 temporary stations, to attempt to extract near-field body waves from ambient seafloor noise observed by multivariate sensors of broadband and short-period seismometers, differential pressure gauges (DPGs), and hydrophones. Our results show that P waves are extracted only in the DPG-record correlations at a frequency of 0.2–0.5 Hz, which can be seen up to a separation distance of two stations of 17 km with an apparent velocity of 3.2 km/s. At 1–3 Hz, P waves are observed only in the vertical-record correlations up to a separation distance of 11 km with an apparent velocity of 2.0 km/s. These velocity differences reflect the vertical velocity gradient of the accretionary prism, because the P waves at low frequencies propagate at relatively long distances and therefore the turning depth is greater. Moreover, the long-period and short-period P waves are observed at the slope and flat regions on the accretionary prism, respectively. To investigate the retrieved wavefield characteristics, we conducted a two-dimensional numerical simulation for wave propagations, where we located single sources at the sea surface above the flat and slope bathymetry regions. Based on our observations and simulations, we suggest that the retrieval of near-field body waves from ambient seafloor noises depends on the relative amplitudes of P and other surface waves in the ambient noise wavefield, and those are controlled by the subseafloor velocity structure, seafloor topography, and water depth.


INTRODUCTION
Ambient noise analysis applied to land-based seismic records has retrieved various wavefields propagating between two receivers (Wapenaar, 2004). Retrieved waves are mainly surface waves (e.g., Sabra, 2005;Shapiro et al., 2005), but body wave retrievals by correlating ambient noises from land-based observations have also been reported (Roux et al., 2005;Draganov et al., 2007;Zhan et al., 2010;Poli et al., 2011;Ryberg, 2011;Takagi et al., 2014). In particular, reflections from the 410 and 660 km discontinuities in the upper mantle could be detected from ambient noise records observed in Finland (Poli et al., 2012). Body waves propagating through the deep interior of the Earth, including core phases, have been extracted using globally distributed broadband seismometers (Nishida, 2013). Near the coastline, direct and refracted P waves could be detected by a dense array of seismometers deployed at Long Beach, California, and those waves can be used to estimate the three-dimensional (3D) velocity structure at shallow depth (Nakata et al., 2015).
The difficulty of extracting body waves is primarily caused by large amplitudes of surface waves observed at land stations, which correspond to microseisms excited by ocean swells in the ocean areas (Longuet-Higgins, 1950;Hasselmann, 1963). Retrievals of body waves may be owing to either large amplitudes of body waves excited by ocean swells near the coastline (Nakata et al., 2015) or low amplitude of surface waves in quiet regions that are away from the coastlines. On the other hand, for seafloor observations, if wave-wave interactions of ocean swell persistently excites microseisms including body and surface waves, seafloor sensors may capture such signals because they are close to excitation regions of the microseisms. Waves extracted from seafloor observations are primarily surface waves, including Rayleigh waves, Love waves, their higher modes (e.g., Takeo et al., 2013;Lin et al., 2016;Isse et al., 2019;Kawano et al., 2020), Scholte waves (Mordret et al., 2014), and ocean acoustically-coupled Rayleigh (ACR) waves (Ewing et al., 1957;Sugioka et al., 2001;Butler and Lomnitz, 2002;Butler, 2006). Here, ACR waves (or seismoacoustic modes) can be observed at frequencies of 0.5-5.0 Hz, which include higher modes of Rayleigh waves whose energies are distributed in the ocean and marine sediment, and have a propagation velocity slightly less than 1.5 km/s (Tonegawa et al., 2015). However, teleseismic body waves excited in the ocean areas can be observed at land stations (Gerstoft et al., 2006;Koper et al., 2010;Landès et al., 2010;Gualtieri et al., 2014;Farra et al., 2016;Nishida and Takagi, 2016). Although this means that most of the body wave energy is transmitted to the interior of the Earth and can be observed at distant stations, near-field body waves are also possibly observed under the conditions of dense seismic sensors deployed at the seafloor.
In the Nankai subduction zone, south of Japan, the Philippine Sea Plate (PHS) subducts northwestwards from the Nankai Trough, historically leading to megathrust earthquakes along the plate boundary. To investigate the seismic structure of the subduction zone, seismic exploration surveys have been conducted with dense survey lines, in which temporary ocean bottom seismometers (OBSs), each equipped with a hydrophone, have been deployed for refraction surveys (e.g., Nakanishi et al., 2018). Moreover, to monitor seismic activity in this region, a permanent cabled-network of seismometers and pressure gauges (Dense Oceanfloor Network System for Earthquakes and Tsunamis: DONET) Kawaguchi et al., 2015) has been installed on the accretionary prism where the bathymetry from the trough gradually becomes shallower landwards to a distance of 30-40 km, after which it is almost flat until 20-30 km before the coastline.
To explore the retrieval of body waves propagating between two seafloor receivers, we employed ambient noise records acquired by seismometers and pressure sensors of DONET and temporary OBSs. Such retrievals have potential for investigating the 3D seismic structure beneath the seafloor without natural and artificial seismic sources. In this study, we show near-field P-wave extractions at frequencies of 0.2-0.5 Hz from seafloor pressure gauges, mainly deployed at the slope bathymetry region, and those at frequencies of 1-3 Hz from seafloor seismometers, mainly deployed at the flat bathymetry region. Those extractions are further examined by 2D numerical simulations.

Station Data
The continuous records used in this study consist of threecomponent seafloor motions and pressure fluctuation and were acquired by broadband seismometers and differential pressure gauges (DPGs) at each station of DONET and 4.5 Hz short-period sensors and hydrophones at each temporary station. The broadband seismometers (Guralp CMG-3T) of DONET were buried 1 m below the seafloor, and have a flat velocity response from 100 Hz to 360 s (e.g., Nakano et al., 2013). The response of the hydrophone decreases from 2 Hz to lower frequencies (e.g., Tonegawa et al., 2015). Also used were the three components of a broadband seismometer deployed in a borehole at a depth of 900 m from the seafloor (Kopf et al., 2011;, at which lower noise levels than those at the seafloor are observed due to the amplitude decay of persistently propagating surface waves. The sampling rates of all sensors are decimated to 40 Hz. The stations in DONET are installed in the eastern (DONET1) and western (DONET2) part of the Nankai subduction zone with a station spacing of 15-20 km, whereas the temporary stations are distributed along five lines that cover the central to eastern part of the subduction zone with a station spacing of 5 km (Figure 1). The observation periods of DONET and temporary stations are 2011∼present and September-December of 2011, respectively. We prepared two datasets of the continuous records to retrieve wavefields of lower (0.2-0.5 Hz) and higher frequency (1-3 Hz) components using an ambient noise analysis.

Cross-Correlation Function Calculation
Cross-correlation functions (CCFs) were calculated with ambient noise records in which energetic signals including earthquakes were suppressed by the following log-normal shaped function.
where σ 2 and the time length, T 400 s. When σ > 1 in Equation 1, the peak of the function is located forward and its tail amplitude is still preserved at the end of the time length (Supplementary Figure S1). The time length, T, was determined from the durations of relatively large earthquake signals at 0.2-0.5 Hz observed by DONET. Supplementary Figure S1A shows a waveform example for a deep earthquake with a magnitude of 5.7, a depth of 410 km (earthquake catalog from United States Geological Survey), and an epicentral distance of 3.67°. The coda amplitudes can be observed for a duration of ∼400 s. If the coda portions of earthquakes are longer than 400 s and still have large amplitudes (see Equations 2, 3 for amplitude criteria), Equation 1 is repetitively applied to the rest of the coda. A cosine taper with a time window of 20 s was also applied to both the edges of the function. The root-mean-squared (RMS 1hour ) amplitudes were calculated with 1-h continuous records at a frequency of 0.2-0.5 Hz. When amplitudes in the 1-h record at 0.2-0.5 Hz, A max , exceed five times the RMS 1hour , the raw records were divided by the following function with a time shift of 80 s from the time of A max : The reason of the time shift is that the maximum amplitude of the log-normal shaped function (Equation 1) is approximately 80 s from the starting time. In addition, when moderate-sized earthquakes occur in Japan, the large and small amplitudes of the S and P waves within an S-P time of <80 s are observed, and the S wave amplitude often corresponds to A max . Equation 2 with a time shift of 80 s and Equation 3 can effectively suppress the large amplitudes of such earthquake signals. An example of the suppression of the deep earthquake signals is shown in Supplementary Figure S1D. The CCFs were calculated using a time window of 200 s with spectral whitening (Brenguier et al., 2007). The CCFs were stacked over 2 months. When the RMS of one segment (200 s) in the processed records was less than 0.3 times the RMS 1hour , the CCFs of those segments were discarded. Previous studies calculated CCFs using continuous records with one-bit normalization where amplitudes are equalized while preserving their polarities. Using the continuous records acquired by a borehole sensor and a nearby seafloor sensor (KMD16, Figure 1) with a horizontal separation distance of 3.7 km, we compared the resulting CCFs ( Figure 2). Without one-bit normalization they show P waves and ACR waves at lag times of ±2 s and +4 s, respectively, while with one-bit normalization the CCFs do not show any clear signals. The reasons for the absence of clear signals remains elusive, but we suppose that larger and smaller amplitudes in the ambient seafloor noise are dominated by different wavefields, as will be discussed in Discussion. In this study, we preserved amplitude information in the continuous records with suppressing energetic signals when calculating CCFs. Because CCFs calculated with spectral whitening measure the phase difference between two time series (e.g., Prieto et al., 2009;Tonegawa et al., 2009), the phase difference obtained in the CCFs reflects portions of the time series with relatively large amplitudes. Figure S2 shows CCFs for 4 × 4 components at 0.2-0.5 Hz. Most of the component combinations only show ACR wave propagations with a propagation velocity of 1.5 km s −1 or surface waves with slower velocities, whereas pressure-pressure (P-P) CCFs show body wave propagations. In Figure 3, the reference station is located at a lower latitude, so that signals in the positive lag time represent waves propagating northwards. P wave propagations are observed up to a separation distance of 17 km in the positive lag time along the travel time curve of the P waves estimated from the velocity structure (Tonegawa et al., Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 610993 2017). Here, we constructed a one-dimensional (1D) velocity model averaged over all of the 1D velocity profiles obtained from the DONET stations below the seafloor (Tonegawa et al., 2017). The P wave can be observed in the positive lag time of the CCF for the station pair of KMB06 and KMB08 ( Figure 3D). The turning depth of P waves reaches up to 6 km from the seafloor (Figure 3), which samples the plate boundary in the shallow subduction zone. Figure 4 shows P-P and vertical-vertical (Z-Z) CCFs at 1.0-3.0 Hz. Although P waves were observed in the P-P CCFs at 0.2-0.5 Hz, in the P-P CCFs at 1.0-3.0 Hz only ACR waves were observed. Instead, P waves were retrieved in the positive and negative lag times of the Z-Z CCFs at 1.0-3.0 Hz, and they reach up to 11 km in separation distance of two stations. At the middle frequency range of 0.5-1.0 Hz, we did not extract body waves (Supplementary Figure S3), and hence focus on body wave retrieval at the frequency bands of 0.2-0.5 Hz and 1.0-3.0 Hz in the subsequent sections.

Supplementary
To confirm the stability of the obtained CCFs, we compared them for stacking periods of 3 months, 1 month, and 2 weeks in the low frequency band (Supplementary Figure S4). Because the observation period of Dataset 2 was almost 2 months, we did not confirm the stability of Dataset 2. The CCFs were almost stable over all of the stacking periods; hence, we discuss the characteristics of the retrieved waves based on a 2 months stacking period.

DISCUSSION P Wave Retrieval
We measured the apparent velocities of P waves in high and low frequency bands using a slant stack technique. Given apparent velocities, the theoretical travel times can be calculated using the distances between two stations. The apparent velocity was varied between 1 km/s and 5 km/s, with an increment of 0.1 km/s. We averaged the absolute amplitudes of the CCFs within a 0.5 s time window from the theoretical travel time. The averaged values were stacked for the separation distances between two stations: less than 10 and 17 km for the high and low frequency bands, respectively, because large amplitudes were obtained in the distance ranges (Figures 3, 4). This process was performed for positive and negative lag times. We obtained velocities of 2.0 km/s in the high frequency band ( Figure 4) and 3.2 km/s for positive lag times in the low frequency band (Figure 3). This velocity difference and the fact that P waves in the low frequency band propagate over relatively long distances indicate that the turning depths of these P waves were relatively shallow and deep, respectively, and the apparent velocities reflect the P-wave velocity structure (Vp) at these depths. Indeed, the travel time curve gradient of the P waves gradually increased as a function of distance ( Figure 3B), which reflects the vertical velocity gradient of the accretionary prism. Seismic exploration surveys of the entire Nankai accretionary prism have also reported a Vp of 2.0-4.0 km/s in the marine sediments and a gradual increases with depth of the Vp structure at shallower depths in the accretionary prism (e.g., Kodaira et al., 2002;Takahashi et al., 2002;Nakanishi et al., 2008). In the higher frequency band, because the turning depth is relatively shallow, the turning P waves are simply observed at stations with separation distances less than 11 km. For the low frequency band, because the turning depth is close to the plate boundary, if seismic velocity discontinuities are present near the bottom of the prism, refracted and head waves can be generated. Indeed, the presence of a low velocity layer has been reported at the bottom of the accretionary prism in the southern DONET1 region (Park et al., 2010;Kamei et al., 2012;Akuhara et al., 2020). In this study, for the low frequency observations from DONET1, the retrieved P waves may contain such multiple P phases.
In order to investigate the region where P waves were retrieved at 0.2-0.5 Hz, we selected P-P CCFs that contain P waves by cross-correlating between the reference CCF and individual P-P CCFs, which is a similar approach to that of a previous study (Nakata et al., 2015). The reference CCF was calculated by stacking the P-P CCFs of all station pairs within separation distances of 10-17 km along the travel time curve of the P waves shown in Figure 3B (Tonegawa et al., 2017). We calculated the cross-correlation coefficients (CC) between the reference CCF and individual P-P CCFs with a time window of ±2 s from the travel time curve. If CC > 0.6, we consider that the CCF possibly contains P waves, and plot the pair in the map (red Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 610993 lines in Figure 5A). As a result, pairs showing P waves are primarily concentrated at the southern part of DONET, where the seafloor slope to the trough is formed. To explore the regions where P waves are observed at 1.0-3.0 Hz, we aligned Z-Z CCFs with the separation distance of two stations of 5 km along Lines A-D (Figure 1). Because the separation distances in the temporary OBS array were equal, aligning the CCFs with 5 km separation distances along lines allowed us to compare the retrieved waves with the geological setting. The result shows that P waves were primarily retrieved at the flat bathymetry region, while signals were weak at the slope bathymetry region (Figure 6). This contrasts with the retrieval of P wave at 0.2-0.5 Hz. Note that the retrieved P waves in the Z-Z CCFs are different from the ACR waves detected in the P-P CCFs using hydrophone records in a previous study (Supplementary Figure S5) (Tonegawa et al., 2015), in which the P waves in the Z-Z CCFs were faster.
Previous studies that retrieved body waves at land stations speculated that body waves are converted from the Rayleigh waves due to the heterogeneous structure of the Earth's upper crust (Roux et al., 2005), and that there are specific structures where body wave energy is trapped and scattered, such as low velocity layers, basins, topography, and heterogeneity (Zhan et al., 2010). Although the retrieved P waves retrieved in our observations may have included contributions from the correlation of P waves trapped and scattered by heterogeneities within the accretionary prism toe, in which the original P waves were generated by wave-wave interactions at the sea surface, the wavelength of the P wave at low frequency was 15 km, for a frequency of 0.2 Hz and a propagation velocity of 3.0 km/s. This wavelength appears to be long for sufficient seismic wave scattering. Therefore, we conducted numerical simulations to evaluate the contribution of the original P waves associated with wave-wave interactions to the P waves extracted from our observations.

Wave Propagation From Numerical Simulation
In this section, we confirm whether the characteristics of the waves retrieved from ambient seafloor noises can be reproduced by a simple 2D numerical simulation for the ocean-solid Earth system. Since it appears that P retrieval is related to the slope and flat regions in the bathymetry, we compare wave propagations depending on frequency (0.2-0.5 Hz and 1-3 Hz), source locations (slope and flat), and component (vertical velocity and stress τ zz ). We used a 2D finite difference method with a rotated-staggered grid for second order approximations in time and space (Saenger et al., 2000). The calculation has been performed in the displacement-stress scheme with an absorbing boundary condition (Clayton and Engquist, 1977), and converted the vertical displacement to vertical velocity seismograms to compare with the τ zz component. We applied vertical single forces with Ricker wavelets for central (maximum) Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 610993 frequencies of 2.0 Hz (2.9 Hz) and 0.35 Hz (0.5 Hz), respectively, to produce wavefields at high and low frequencies. The grid size is 10 × 10 m and 20 × 20 m for high and low frequencies, respectively. Stations are set at the seafloor within a distance ranging between 0 and 200 km with an interval of 1 km. The seismic velocity structure at the subseafloor is based on a Vp structure resulting from a seismic exploration survey (Figures 1, 7A) (Nakanishi et al., 2008). The 1D Vp profile at the northern edge of the model is extended landward (Figure 1). The Vs and density were derived from Vp using empirical relations (Brocher, 2005). A 1D profile (depth dependent) of the Vp is applied to the sea water (Munk, 1974), while Vs and density are 0 and 1.02 g/cm 3 , respectively. The calculation was unstable in the cases where Vp/Vs of the accretionary prism at shallow depths was large and the small horizontal-scale bathymetry was complex. To avoid these problems, we set Vp/Vs 2.5 when Vp/Vs > 2.5 and applied a horizontal distance moving average of 10 km to the seafloor. The minimum Vs is 0.56 km/s because the minimum Vp in the sea water is approximately 1.4 km/s (Munk, 1974). The source locations were set to horizontal distances of 80 and 110 km at the sea surface, which correspond to the slope and flat bathymetry regions, respectively ( Figure 7A).
The numerical simulation results at high frequency shows that direct and subsequent P waves are propagating within the subseafloor structure, and were produced by multiple reflections of P waves in the sea water (e.g., red arrow in Figure 7B). ACR waves can also be observed after the apparent velocity of 1.5 km/s (e.g., blue arrow in Figure 7C). At low frequency, in addition to P waves, subsequent P waves, and ACR waves, the propagation of Rayleigh waves with an apparent velocity of ∼0.5 km/s was also extracted (e.g., orange arrow in Figure 7D). After the direct and subsequent P waves propagate to a horizontal distance of 30-40 km ( Figure 7B), their amplitudes at both frequency bands are decayed at greater distances. Although our observations only show P wave propagations of 10-15 km in horizontal distance, such an amplitude decay as was obtained in the numerical simulation is one of the reasons for the absence of P wave propagations at farther distances.
We consider that the retrieved waves in our observation can be linked to the relative amplitudes of the waves that emerged in the numerical simulation. In our observation, P waves could be retrieved at a frequency of 1-3 Hz at the stations on the flat bathymetry. In the numerical simulation for the flat bathymetry, the amplitudes of direct and subsequent P waves in the vertical velocity component are larger than those in τ zz , and the amplitudes of the ACR waves are large in τ zz ( Figure 7B).
Because we preserved the amplitude information when calculating CCFs, large amplitudes of those P and ACR waves may be correlated and emphasized in the Z-Z CCFs and P-P CCFs, respectively. Moreover, the amplitudes at the coda part of the ACR waves are large in the slope region ( Figures 7B,C). This may hinder the retrieval of P wave at the slope region in our observation.
For the low frequency wavefields in our observation, P waves could be extracted at the slope region and mainly propagate northwards. In the numerical simulation for the bathymetric slope, direct and subsequent P waves propagating to shallow water depths have larger amplitudes than those propagating to deep water depths (arrow in Figure 7E), while the amplitudes to both directions are comparable in the flat bathymetry region ( Figure 7D). This asymmetric radiation pattern of the P amplitudes may cause the azimuthally-asymmetric extraction of the retrieved P waves in our observation. Another important issue of the low-frequency P retrieval is that P waves could not be extracted in the observed Z-Z CCFs, although they have large amplitudes in the vertical velocity component in the numerical simulation ( Figure 7E). This is because Rayleigh wave amplitudes significantly emerged in the vertical velocity component, compared with τ zz (Figure 7E), and ACR and Rayleigh waves propagate long distance with relatively FIGURE 6 | P wave retrieval at 1-3 Hz. (A) Z-Z CCFs for pairs whose separation distance is 5 km are aligned along line A (Figure 1). Red triangles represent P wave retrievals, and red arrow indicates the location where the P wave is retrieved, which corresponds to the flat bathymetry region. Right panel shows the bathymetry along the line A. Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 610993 9 little attenuation, compared with P waves. Thus, it is considered that cross-correlating the vertical velocity components at the seafloor tends to result in the retrieval of Rayleigh and ACR waves. For the flat region in the simulation, because the Rayleigh waves have large amplitudes ( Figure 7D) and propagate from the flat to slope regions, this effect also supports the idea that ACR and Rayleigh waves dominated the observed Z-Z CCFs at this frequency.
The ambient noise wavefield at the seafloor contains body waves and multiple modes of surface waves. Since their amplitudes vary with the geological setting, including subseafloor velocity structure, seafloor topography, and water depth, the summed amplitudes of the waves contribute to seafloor motions below the seafloor and pressure variations above the seafloor. This may result in wavefield differences across the seafloor in terms of amplitude and propagation velocity. We therefore considered that the wavefield differences in the solid Earth and ocean controls the extracted wavefield in the CCFs obtained from seismometers and pressure sensors in this study.

Cross-Correlation Functions With and Without One-Bit Normalization
The one-bit normalization is capable of equalizing the energies of incoming waves from heterogeneously distributed sources, and useful for cross-correlating ambient noise records (Brenguier et al., 2007). We here compared the obtained CCFs with and without one-bit normalization at higher and lower frequencies in more details. Supplementary Figure S6A shows the P-P CCFs at 0.2-0.5 Hz with one-bit normalization. The propagations of the P and ACR waves are almost comparable to those in the P-P CCFs without one-bit normalization (Figure 3). On the other hand, the signal-to-noize ratios (S/Ns) of the P and ACR waves in the Z-Z CCFs at 1-3 Hz with one-bit normalization (Supplementary Figure S6B) are significantly lower than those without one-bit normalization ( Figure 4A). This feature is consistent with the example shown in Figure 2, in which P and ACR waves can be constructed in 1-h CCFs without one-bit normalization.
Although the reason for such low S/Ns in the retrieved waves remains elusive, one possible reason may be related to the relative amplitudes between the wavefield right after excitation by wavewave interaction and the wavefield where the excitation occurred some time ago. The wavefield associated with wave-wave interaction contains primary waves due to the wave-wave interaction with relatively large amplitudes and secondary waves that contain multiply scattered waves of the primary waves. For example, the primary waves include direct and subsequent P waves and ACR waves, and the secondary waves are dominated by ACR waves, particularly at higher frequencies, scattered by both small-scale heterogeneous structures below the seafloor and complicated seafloor topography. If we calculate CCFs without one-bit normalization, the CCFs possibly contain the effects of the primary waves. However, if we calculate CCFs with one-bit normalization, because the amplitudes of primary and secondary waves are equalized, correlating those waves presumably requires a larger amount of ambient noise records. This may also be related to the rapid convergence to a robust CCF when non one-bit normalization was applied (Seats et al., 2012), although the time window used for calculating the CCFs did not overlap in the present study, as in the case of Seats et al. (2012).
To evaluate this, CCF constructions from ambient noise in various marine geological settings are required. In particular, because the relative amplitudes of the ambient noise wavefields at the seafloor may be controlled by the subseafloor velocity structure, seafloor topography, and water depth, it is necessary that the retrieval of body waves and the convergence from multivariate components is analyzed with and without one-bit normalization in various marine settings. Such experiments may also be applicable to body-wave extractions from land-based three-component seismometer records.

Acoustically-Coupled Rayleigh or S Wave?
At 0.2-0.5 Hz, clear signals can be traced in the positive and negative lag times along the travel time curve of the S wave up to a separation distance of two stations of 23 km (Figure 3). At separation distances greater than 23 km, the apparent velocity of the relatively weak signals is close to 1.5 km/s and slower than the travel time curve of the S wave, which corresponds to ACR waves. If the waves observed at distances up to 23 km corresponds to S waves, the pressure gauges observe the pressure response associated with the seafloor displacement corresponding to S waves from the subseafloor structure, including incident S, reflected P, and reflected S waves.
However, it appears that these waves also correspond to ACR waves. Our numerical simulation indicates that ACR waves have higher phase velocities than the group velocity of ∼1.5 km/s (arrow in Supplementary Figure S7), which reflect the shear wave velocity at the deeper part of the accretionary prism. It is therefore considered that the observed waves corresponding to the S-wave travel time curve are the results of cross-correlating the ACR waves with high phase velocities. The reason for the apparent velocity change at a distance of 23 km may be that the wave propagations at distances greater than 23 km primarily reflect the seismic structure beneath the flat region because relatively longer separation distances of two stations can be used in this region.

CONCLUSION
We present the retrieval of near-field P waves from ambient noise records observed at the seafloor, using data from a permanent cabled network (DONET) and temporary stations. The following are the major findings on the characteristics of P retrievals.
(1) P waves could be extracted in the P-P CCFs at 0.2-0.5 Hz in the slope bathymetry region, and they propagate to shallow water depths.
(2) P waves could be extracted in the Z-Z CCFs at 1-3 Hz in the flat bathymetry region. (3) Our numerical simulations indicate that the relative amplitudes among the P, ACR, and Rayleigh waves and their attenuations are important for the extraction of the waves. (4) The relative amplitudes of these waves are controlled by the marine geological setting, including the subseafloor seismic velocity structure, seafloor topography, and water depth.
The propagation distance of P wave is only up to 15-20 km, but their bottoming depth reaches up to 5-6 km below the seafloor. This indicates that the P wave can reach the plate boundary in the shallow Nankai subduction zone. If such retrievals can be realized in various shallow subduction zones, the retrieved P waves can be used to investigate seismic structures by, e.g., tomographic approaches. It is expected that more details on near-field body wave extractions from ambient seafloor records will be investigated under various conditions in seismic seafloor experiments.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/ restrictions: Seafloor seismometer data are available from JAMSTEC upon request. The data that support the findings of this study are available from the corresponding author upon reasonable request. Requests to access these datasets should be directed to Takashi TONEGAWA, tonegawa@jamstec.go.jp.

AUTHOR CONTRIBUTIONS
TT processed the data, and drafted the manuscript. TK and EA acquired the data underlying this study. All authors contributed to interpretation and the final manuscript.

FUNDING
This work was supported by JSPS KAKENHI Grant Number 19H04632 in Scientific Research on Innovative Areas "Science of Slow Earthquakes."