Abstract
Determination of a response of the sea water column to teleseismic plane wave is important to suppress adverse effects of water reverberations in calculating receiver functions (RFs) using ocean-bottom seismometer (OBS) records. We present a novel non-linear waveform analysis method using the simulated annealing algorithm to determine such a water-layer response recorded by an OBS array. We then demonstrate its usefulness for the RF estimation through its application to synthetic and observed data. Synthetic experiments suggest that the water-layer response constrained in this way has a potential to improve RFs of OBS records drastically even in the high-frequency range (to 4 Hz). By applying it to data observed by the OBS array around the Kii Peninsula, southwestern Japan, we identified a low-velocity zone at the top of the subducting Philippine Sea plate. This zone may represent the incoming fluid-rich sediment layer that has been reported by active-source seismic survey.
1 INTRODUCTION
Scattered teleseismic phases, such as P-to-S conversion phases, have illuminated the Earth's interior because of their high sensitivity to seismic velocity discontinuities. Conventionally, such phases are extracted from teleseismic seismograms by deconvolving the radial or transverse component with the vertical-component record (Langston 1979). The resultant time-series are termed receiver functions (RFs). The RFs provide useful approximations of the Earth's impulse response.
Recent progress in passive seismic monitoring at offshore regions has offered more opportunities to conduct RF analysis using data of ocean-bottom seismometers (OBSs; e.g. Kawakatsu et al. 2009; Kumar et al. 2011; Akuhara & Mochizuki 2015; Audet 2016 and references therein). The application of RF analysis to OBS data, however, is not straightforward, especially in the high-frequency range where the ocean water layer yields strong reverberations dominating vertical-component records. Such reverberations violate the fundamental assumption of the RF analysis that the source wavelets can be approximated by vertical-component records.
Akuhara & Mochizuki (2015) developed a method to suppress the effect of the water reverberations, utilizing a linear filter that expresses the impulse response of the water layer, namely the water-layer filter (WLF). The inverse of WLF (IWLF) is used to eliminate the multiples. The WLF was formulated by assuming vertical incidence of seismic waveforms, requiring only two tuning parameters to fully describe the filter shape: a P-to-P reflection coefficient at the seafloor, R, and a two-way traveltime within the water layer, τ. Akuhara & Mochizuki (2015) determined these parameters by first estimating the source wavelet from on-land seismic station array data and then by deconvolving OBS vertical components with the source wavelet. Such an on-land array is, however, not always available at close distance from an OBS network.
This study aims to improve the WLF method to be more accessible and to demonstrate its usefulness. For this purpose, we first propose an alternative method to derive water-layer responses using only OBS data. We conduct non-linear waveform analysis to determine source wavelets and water-layer response simultaneously employing the simulated annealing (SA) method. We then show how the WLF method improves RFs using both synthetic and real observed data in a high-frequency range (up to 4 Hz).
2 NON-LINEAR WAVEFORM ANALYSIS FOR WATER-LAYER RESPONSE
We employ the SA algorithm (Kirkpatrick et al. 1983) to address this high-dimensional non-linear inverse problem. We perform a random walk with a constant step interval (1 per cent of the maximum amplitude of observed waveforms) to generate trial model parameters for source wavelets, while the other parameters are generated via random samplings within a certain continuous ranges, following previous studies (Chevrot 2002; Iritani et al. 2010, 2014; Tonegawa et al. 2013). We use the same tuning parameter as Iritani et al. (2010) through which the SA algorithm controls probability to accept trial parameters. The final results are obtained after 2000 random generations of each model parameter. Repeating this inversion eight times using different random seeds yield no significant change (Fig. S1, Supporting Information), suggesting that our inversion results most likely converge to the global minimum rather than local minimum. In the following result section, we show the averaged τ and R values over these eight estimates.
3 APPLICATION OF NON-LINEAR WAVEFORM ANALYSIS FOR OBS ARRAY DATA
We use data from the OBS array that was deployed around the Kii Peninsula from 2003 to 2007 (Fig. 1a). The array consists of OBSs with velocity sensors characterized by a flat frequency response above 1.0 Hz. The observation periods were different among observation sites, while 9–27 OBSs were deployed at the same time (Mochizuki et al. 2010). The orientations of the horizontal-component sensors were determined from the Rayleigh wave polarization (Akuhara & Mochizuki 2015). From continuous records of these OBSs, we extract P-wave records of teleseismic earthquakes (M > 6.0) at a distance range of 30°–90° and local deep (>300 km) events (M > 5.0) at a distance range of 3°–10°. We measure signal-to-noise ratio (SNR) on the vertical-component records while conducting automatic phase picking based on short-term average to long-term average ratio algorithm (Abt et al. 2010). The SNR is measured by taking 1 and 4 s long time windows before and after the phase pick, respectively. We adopt events with SNR ≥ 3.0 at eight or more OBSs for the inversion analysis to determine a water-layer response at each site. We obtained 118 events through this selection (Fig. 1b).
(a) Sites of ocean-bottom seismometers (OBSs) used in this study. Grey circles represent OBS sites. Thin lines show the water depth with the depth interval of 500 m. (b) Teleseismic and deep local events used in the inversion analysis.
After solving the inversion problem, we evaluate the quality of solutions by measuring cross-correlation coefficients (CCs) between synthetic and observed waveforms at each OBS (Fig. 2a). For the following discussion, we discard results of poorly fit seismograms with CC < 0.8.
Results of non-linear waveform analysis. (a) Observed (black) and synthetic (red) waveforms for the event with origin time specified on top of the panel. OBS names and cross-correlation coefficients are labeled on the left side of these traces. These traces are aligned in the order of OBS depths. (b) Source wavelet (red trace) and impulse response at each OBS site (black traces) from inversion analysis. Green thick bars represent the arrival times of the first water multiples expected from the OBS depths. (c) Estimations for τ and R values from different events at a single OBS site, LS11 (crosses). The star represents their averaged values with 2σ standard errors. (d) and (e) Correlation diagram between OBS depth and averaged (d) τ and (e) R values for each OBS. Error bars denote 2σ standard errors, but are scaled by 10 for τ.
The inversion results show good reproduction of observed waveforms even if the onsets of water reverberations are contaminated by coda of the direct arrivals (Fig. 2a). This insensitivity to coda contamination can be considered one of the advantages of this non-linear waveform fitting approach using the SA algorithm (e.g. Iritani et al. 2010). We also find consistency between our estimations for τ values and the expected timing of the reverberations from OBS depths measured by acoustic ranging (Fig. 2b).
Even at the same site, individual τ and R values from different events are scattered (Fig. 2c). If this scattering comes from the dipping seafloor, there should be some correlation between individual τ (or R) values and the locations of the incident point where P wave enters into the water column. However, we do not find any apparent correlation (Fig. S2, Supporting Information). We therefore assume that the scattered τ and R distributions mainly reflect random errors in the measurements.
Examining the results from all OBSs, we can see a strong correlation between the averaged τ values and the OBS depths (Fig. 2d). The slope of this linear trend is about 0.75 km s−1 that is consistent with the P-wave velocity of 1.5 km s−1 within the water layer. Our estimations for an average R value at each OBS site show large variation from 0.1 to 0.5 (Fig. 2e). The averaged value among all OBSs is 0.30, which is reasonable for the typical P-wave velocity and density of the marine sediment (Hamilton 1978).
4 RECEIVER FUNCTION ESTIMATION WITH INVERSE WATER-LAYER FILTER: SYNTHETIC TESTS
We calculate synthetic vertical- and radial-component records using the propagation matrix method (Haskell 1953) by assuming 1-D layered structure, composed of a water column, sediment and crust layers and a mantle half-space (Table S1, Supporting Information). In this calculation, a receiver is allowed to be placed beneath the free surface, so the effects of water reverberations are taken into consideration (e.g. Kumar et al. 2011). Although the formulation of the WLF (eq. 3) is based on the assumption of vertical incidence, here we employ oblique incidence with ray parameter of 0.06 s km−1 to check the validity of this assumption.
We compute two types of radial RFs by the deconvolution of the synthetic radial-component record with the synthetic vertical-component record pre-processed with IWLF and the original synthetic vertical-component record. The extended-time multitaper method (Shibutani et al. 2008) is incorporated through deconvolution operation. For the parameters of the IWLF, we use not only theoretical values derived from the layered model, but also slightly shifted values from the theoretical values. These shifted values aim to examine how poorly constrained parameters influence RF estimations. The resultant RFs are low-pass filtered to 4 Hz, and RFs estimated with IWLF are normalized by the factor of 1 + R for comparison. This amount of 1 + R corresponds to the energy of downward P-wave reflection at the seafloor which is to be removed from vertical-component records by the IWLF.
The synthetic tests show distinct advantages of the usage of the IWLF (Fig. 3a). The original RF shows many artificial peaks throughout the records which cannot be explained from the layered structure model (e.g. positive phase at 0.6 s and negative phase at 2.7 s). The RF processed with IWLF, on the other hand, shows no such artificial peaks but six isolated phases which can be interpreted as P-to-S conversion phase at the Moho (PsM) and the bottom of the sediment layer (Ps), and sediment-related reverberations (PpPs, PpSs, PsSs and PpPs + w, see Fig. 3b).
(a) Synthetic receiver functions (RFs). Three panels show results from different water-layer filter parameters. In each panel, red and grey traces are synthetic RFs calculated with and without the application of the inverse water-layer filter to vertical-component records, respectively. Blue solid and dashed bars represent theoretical arrival times of primary phases with positive and negative RF amplitudes, respectively. (b) Schematic illustration of primary phases which are recognizable on the synthetic tests.
We achieve the suppression of artificial peaks even when we use the shifted τ and R values. Determining the thresholds for acceptable shift amount is difficult because it depends on structure model and frequency range of low- or bandpass filter. From our experience, the acceptable shift amounts are about ±0.1 s and ±0.3, for τ and R values, respectively. For most stations, the standard errors for τ and R values are below these empirical values and thus these estimation errors would not cause a severe problem in RF estimation (Figs 2d and e).
5 RECEIVER FUNCTION ESTIMATION WITH INVERSE WATER-LAYER FILTER: APPLICATION TO OBS RECORDS
We computed radial RFs for two OBSs, LS11 and LS08, using the same method as in the previous section, regarding the deconvolution operation, filtering and normalization. Teleseismic or deep local events used here are selected with the same criteria described in Section 3, but we calculate SNR for 30 s long time windows before and after the theoretical P-wave arrival times that are expected from the IASP91 Earth model (Kennett & Engdahl 1991). About 50 event records with SNR ≥ 3.0 on vertical components are selected for each of these two OBSs (cf. about 20 event records are selected in Section 3).
We can see significant differences between the RFs resulting from the two methods (Fig. 4). We conclude that the WLF method successfully improved the RF estimations based on the following three observations: (1) reduction of artificial energy before the direct P-wave arrivals; (2) disappearance of negative amplitudes at zero lag time at LS08, which are also considered as artificial peaks and (3) improved coherence of RFs among different events (e.g. ∼1.6 s at LS08).
Radial receiver functions (RFs) calculated at two different sites, (a)–(e) LS08 and (f)–(j) LS11. (a) and (f) Backazimuth (orange circles) and ray parameter (green circles) distribution which are used to estimate RFs. (b) and (g) Radial RFs calculated without the application of IWLF to vertical-component records. (c) and (h) Radial RFs calculated with the application of IWLF. (d) and (i) Enlarged views of (c) and (h). Open circles represent the negative RF peaks. (e) and (j) Moveout pattern of the negative RF peaks along event backazimuths. Colour of each plot indicates the amplitudes of the peaks. Blue curves represent sinusoidal curve fitted to the plots.
For both OBSs, RFs after application of IWLF indicate strong coherent peaks around ∼1.5 s, although their polarities differ from each other. The small delay times suggest that velocity discontinuities exist at very shallow depth beneath the seafloor (≤1 km depth, if assuming typical velocity of the sediment). We interpret the positive amplitudes at LS11 as P-to-S conversion phase from the boundary between the unconsolidated sediment above and the old accretionary prism below (e.g. Tsuji et al. 2015). In contrast, the early negative phase at LS08 suggests the existence of inverted velocity contrast. A candidate for the cause of the negative amplitude phase is P-to-S conversion phase from a bottom simulating reflector, above which gas hydrated layer with high-velocity anomaly is located (e.g. Harris et al. 2013). It is important to note that the very shallow structure can alter the appearance of RFs drastically for the case of OBS records.
By applying the IWLF, we recover coherent signals with negative polarities at 3.1 and 3.3 s for LS08 and LS11, respectively. These negative signals indicate moveout patterns along event backazimuths characterized by a sinusoidal curve with 2π period (Figs 4e and j). These moveout patterns of the negative peaks suggest the existence of low-velocity zones (LVZs) with a northward dip. The depths of these LVZs are considered to be 14 km below sea level at LS08 and 16 km at LS11 if we assume the P- and S-wave velocities estimated by the previous 3-D tomographic analysis (Akuhara et al. 2013). Such depths are in good agreement with the depth of the subducting plate interface estimated by a previous study (Baba et al. 2002). We therefore interpret that these LVZs represent the top of the subducting Philippine Sea plate.
6 DISCUSSION
In southwestern Japan, LVZ along the plate interface has been identified by intense P-to-P reflection phases from active seismic sources and has been interpreted as fluid-rich sediment layer (Kodaira et al. 2002). Although further quantitative analysis is left for future studies, the LVZ revealed by this study may also represent the same fluid-rich layer as reported by Kodaira et al. (2002). Many RF studies have identified LVZs at the top of subducting plates worldwide and have interpreted them as fluid-rich (Bostock 2013 and reference therein). These findings have provided useful insight to better understand slow slip phenomena occurring at the downdip of seismogenic zones on megathrusts (Audet & Kim 2016 and reference therein). In contrast, for seismogenic zones, little information has been obtained by RF methods because of the limited analyses with OBS data. Our proposed method can extend the target region seaward and enables to access seismogenic zones.
We should note the non-uniqueness of our interpretation. As shown in the synthetic tests, a sediment layer beneath the seafloor causes various multiple phases on RFs, which may also be recognizable in the data (e.g. ∼5.0 s in Fig. 4c). Such multiples make interpretations of RFs difficult. Although several lines of evidence, including the moveout patterns and the consistency with the previous plate model, support our interpretation of the LVZ at the plate interface, more sophisticated analysis such as RF inversion analysis would be necessary for robust interpretation.
The frequency range we used in this study (<4 Hz) is higher than that used in ordinary RF studies (<1 Hz). Such high-frequency analysis is helpful to achieve high spatial resolution and thus the results become more comparable with those obtained by active-source seismic surveys. Most recently, detailed properties of a LVZ along the plate interface at shallow depths has been investigated by modeling P-to-P reflection waveforms of active seismic sources (e.g. Li et al. 2015). Our success in high-frequency RF estimation using OBS records provides another option to conduct such waveform modeling using passive seismic sources even for a wide depth range of the plate interface. This new analysis has potential to offer additional constraints on S-wave velocity or even anisotropic structure of the plate interface at seismogenic depth. These kinds of information are essential to discuss material and fluid content at the plate interface.
We thank Pascal Audet and an anonymous reviewer for their useful comments. This work was supported by JSPS KAKENHI Grant Numbers 14J10221 and 15K13558.
REFERENCES
SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this paper:
Figure S1. Histograms of 1σ standard deviations of (a) τ and (b) R over eight inversions using different random seeds.
Figure S2. (a) τ and (b) R values estimated from different events at a typical OBS, LS11. The size of crosses and circles show the amount of the deviations from the averaged values, and their locations represent the incident points of the first water reverberation into the sea water. Centre of each polar coordinate is identical to the OBS location, and dashed circle denotes intervals of 100 m.
Table S1. Layered structure and ray parameter used in synthetic tests.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.
