Detailed S-wave velocity structure of sediment and crust off Sanriku, Japan by a new analysis method for distributed acoustic sensing data using a seafloor cable and seismic interferometry
Express Letter
Shun Fukushima ORCID: orcid.org/0000-0002-9522-91441,2,
Masanao Shinohara2,
Kiwamu Nishida2,
Akiko Takeo2,
Tomoaki Yamada2 &
Kiyoshi Yomogida3
The S-wave velocity (Vs) structure of sediments and the uppermost crust in the landward slope of a subduction zone are important for determining the dynamics of the overriding plate. Although distributed acoustic sensing (DAS) measurements have improved the horizontal resolution of Vs structure in marine areas, the estimations have been limited to the uppermost sedimentary layers. In the present study, we applied seismic interferometry to DAS data of 13 h duration to image the sedimentary and crustal structure offshore of Sanriku, Japan with a spatial horizontal resolution of 2.5 km and > 3.0 km depth. We grouped the DAS data into 10 km long subarrays with 75% overlaps. We first applied a frequency-wavenumber filter to the DAS data to remove DAS instrumental noise and to allow effective extraction of surface waves from short-time records. We then applied a seismic interferometry method and estimated the phase velocities at each subarray. The estimated phase velocities of the fundamental-mode and first higher-mode Rayleigh waves were then used to determine one-dimensional Vs structures for each subarray. The resultant 2-D Vs structure was interpreted as representing sediments and crust. The upper sedimentary layers thicken seaward, while the entire sedimentary unit shows complex lateral variations in depth. The boundary between the sedimentary layers and the uppermost crust varies in depth from 1.8 to 6.8 km and is the deepest in the middle of the profile. Combining this result with the P-wave velocity (Vp) structure along the nearest survey line, determined in previous studies, allowed us to estimate Vp/Vs = 3.12, on average, for the lower sedimentary layers. Our method of applying seismic interferometry to marine DAS data broadens the techniques for estimating Vs and Vp/Vs structure of sedimentary layers and the upper crust across subduction zones. These results show that application of the frequency-wavenumber filtering and seismic interferometry to marine DAS data can estimate the Vs structure and the Vp/Vs structure, together with standard marine geophysical surveys of sedimentary layers and the upper crust across subduction zones.
Accurate S-wave velocity (Vs) structure of sedimentary layers and the uppermost crust in the landward slope of a subduction zone provide important information about the rheology of the overriding plate. Knowledge of the Vs structures in a shallow part should enhance the accuracy of a wide range of studies, such as hypocenter determination and consideration of rock properties. Since the P-wave velocity (Vp) structures can generally have high resolution, a reliable estimation of the S-wave velocities would provide an accurate estimation of the ratio of P- and S-wave velocities (Vp/Vs), which is useful information for the estimation of rock properties (Ayres and Theilen 1999) and fluid pressure properties (Kodaira et al. 2004).
In the subduction zone of the Japan trench, the recent acquisition of marine multichannel seismic surveys has provided high spatial resolution Vp structures of the sediment and upper crust (Miura et al. 2003; Takahashi et al. 2004). By contrast, the estimation of Vs structures has been limited. Consequently, seismic surveys in the Japan trench could provide estimations of the Vs structures based on the P to S conversion waves observed in data from Ocean Bottom Seismometers (OBSs) (Fujie et al. 2018). Although both seismic interferometry and receiver function methods are useful for estimating Vs with OBS records (Yao et al. 2011; Akuhara et al. 2020; Yamaya et al. 2021), obtaining spatially high-resolution Vs models has been difficult because the distance between OBS locations is typically greater than 6 km along active source seismic profiles and greater than 20 km apart for passive surveys.
In recent years, distributed acoustic sensing (DAS) measurements have been applied to seismic observations (e.g., Zhan 2020). DAS system enables us to measure strain or strain rate with very high spatial resolution over a long distance. The DAS measurement was performed on both the land and seafloor. For example, Dou et al. (2017) estimated shallow shear wave velocity structures by applying the ambient noise method to DAS records on land. In a marine area, Spica et al. (2020) estimated the Vs structure of shallow sediments (to 3 km depth) by applying a frequency–wavenumber (FK) analysis to the seafloor DAS measurements obtained off the Sanriku coast of Japan. However, their estimation of Vs structures was limited to a shallow sediments (< 3.0 km depth). Obtaining the Vs structure throughout the crust is important for deriving the Vp/Vs and the forearc rock properties across subduction zones. The phase velocity of the overtone or long period Rayleigh waves is useful for estimating Vs upper crustal layers, since the phase velocities of these waves are sensitive to Vs.
In the present study, we applied a seismic interferometry method to seafloor DAS records to extract the overtone or long-period (> 5 s) Rayleigh waves. We also proposed a practical method using FK filtering to extract surface waves in the low-frequency range from short-time records with FK filtering. We then estimated the phase velocities from the extracted Rayleigh waves. Finally, we obtained a 2-D Vs structural model of the sedimentary layers and the upper crust for the offshore area of the Sanriku coast of northeastern Japan.
Data and methods
In 1996, the Earthquake Research Institute of the University of Tokyo installed a seafloor seismic tsunami observation system that uses an optical fiber cable in the offshore area of Sanriku (Fig. 1). The length of the buried cable is approximately 120 km, and it extends between 0 and 47.7 km from the coast and covers an area with a sea depth from 0 to 2750 m (Shinohara et al. 2022). The system contains six spare fibers; therefore, DAS measurements have been conducted using a DAS interrogator unit from AP Sensing GmBH (Cedilnik et al. 2019) since 14 February 2019. In the present study, the DAS system recorded strain. The data were recorded at a temporal sampling frequency of 500 Hz, and the duration of the observations in this study was selected to be approximately 13 h. The sensing range, spatial resolution, and gauge length adopted in this study were 100 km, 5.1 m, and 40.79 m, respectively. Shinohara et al. (2019, 2022) and Spica et al. (2020) have described the cable setup and measurement quality of this DAS system.
Location of the seafloor cable containing optical fibers in the off-Sanriku region, Japan. Location of the optical seafloor cable in the off-Sanriku region, Japan. The black line indicates the route of the seafloor cable installed by the Earthquake Research Institute, the University of Tokyo. White triangles the indicate centers of the subarrays for calculation of the CCFs. Blue triangles indicate the centers at the distances from the coast of 39.94 and 67.48 km for the data processing demonstrated in Fig. 2 and 3. The red line shows a span of the estimated 2-D velocity profile shown in Fig. 4. Due to large incoherent noise near the coast and offshore areas, the spatial range of the present data processing was limited from 28 to 75 km from the coast. White lines indicate the survey lines of Takahashi et al. (2004)
Calculation of CCFs with an FK filter
Before calculating and stacking the cross-correlation functions (CCFs) of the DAS background noise records, we divided the entire dataset of the DAS array into 10 km-subarrays with a moving window of 75% overlap (Fig. 1). We decimated the original records by first reducing the sampling frequency from 500 to 2 Hz and then stacking 10 adjacent data to enhance the signal-to-noise ratios (SNRs) after filtering the anti-spatial aliasing. The resulting spatial interval of the data was then 51 m.
We computed the CCFs by first dividing continuous records into 10 min time window segments. We allowed for a 50% overlap of each segment to improve SNR. We computed the weighted average of the cross spectra between two channels, following Takagi et al. (2021) and Takeo et al. (2013, 2014). This processing was conducted in the frequency domain. We avoided any earthquake signals by calculating the mean power between 0.025 and 0.2 Hz in the frequency domain for each time window and discarding any data segments whose amplitude exceeded 10 times the amplitude of the previous time window. Takeo et al. (2013) provide more detailed explanations for the calculation of CCFs.
In general, the DAS instrumental noise has a coherent phase at the same time for all stations (Tribaldos and Ajo-Franklin 2021), and this source of noise particularly affects seismic interferometry studies. In the present study, we applied an FK filter to the DAS data to remove this zero-lag noise and to enhance the SNR before calculating and stacking the CCFs. The high spatial density of the DAS data enables the application of spatial Fourier transform; therefore, we performed a 2-D Fourier transform in the FK domain (Hudson et al. 2021; Atterholt et al. 2022). The zero-lag noise, which has an infinite phase velocity, is mapped to the zero wavenumber in the FK domain. We designed an FK filter \(w\left(f,k\right)\) as:
$$\begin{array}{c}w\left(f,k\right)=g\left(f,k,{c}_{\text{min}}\right)*\left(1-s\left(f,k,{c}_{\text{max}}\right)\right), \end{array}$$
where \(g\left(f,k,{c}_{\text{min}}\right)\) is a 10% Tukey window used to eliminate signals slower than the cut-off minimum phase velocity \({c}_{\text{min}}\). \(g\left(f,k,{c}_{\text{min}}\right)\) is defined as:
$$g\left( {k,f,c_{{\min }} } \right) = \left\{ \begin{array} {ll} 1,& {\text{for}}\, - 0.4k_{{\max }} < k < 0.4k_{{\max }} \\ 0, & {\text{for}}\,k < - 0.5k_{{\max }} \;{\text{or}}\,k > 0.5k_{{\max }} \\ \frac{1}{2}\left[ {1 + \cos \left\{ {\pi \left( {1 - \frac{k}{{0.1k_{{\max }} }}} \right)} \right\}} \right], &{\text{otherwise}} \\ \end{array} \right.,$$
where \({k}_{\text{max}}=4\pi f/{c}_{\text{min}}\). Furthermore, \(s(f,k,{c}_{\text{max}})\) is an 8-th order Kaiser filter that removes the signal faster than the maximum phase velocity \({c}_{\text{max}}\), including the zero-lag noise. In this study, we enhanced the surface wave signals by setting \({c}_{\text{min}}\) and \({c}_{\text{max}}\) to 0.35 km/s and 4.00 km/s, respectively.
Phase velocity estimation
In this section (Phase velocity estimation), we obtained the phase-velocity dispersion curves for the Rayleigh waves using the CCFs of the ambient noises recorded by the DAS observation. At the same time, we determined the 1-D Vs structures and errors of the phase velocities. In the next section (Inversion of the S-wave velocity structure), we provide more accurate estimates of the 1-D Vs structures using only the reliable frequency range of the dispersion curves (i.e., those with small errors) obtained in this section.
We adopted the spatial auto-correlation (SPAC) method to estimate the phase velocity (Aki 1957; Okada 2006; Nishida et al. 2008b), because the SPAC method represents the observed cross-spectra assuming a laterally homogeneous structure and a homogenous source distribution. According to Nakahara et al. (2021), the synthetic cross-spectrum for radial strain records (e.g., DAS) is defined as:
$$\begin{array}{c}{S}_{ij}^{\text{syn.}}\left(\upomega ,{c}_{\text{R}},{c}_{\text{L}}\right) = {A}_{\text{R}}\left(\upomega \right)\left[3{J}_{0}\left(\frac{\upomega }{{c}_{\text{R}}}{d}_{ij}\right)-4{J}_{2}\left(\frac{\upomega }{{c}_{\text{R}}}{d}_{ij}\right)+{J}_{4}\left(\frac{\upomega }{{c}_{\text{R}}}{d}_{ij}\right)\right]+{B}_{\text{L}}\left(\omega \right)\left[{J}_{0}\left(\frac{\upomega }{{c}_{\text{L}}}{d}_{ij}\right)-{J}_{4}\left(\frac{\upomega }{{c}_{\text{L}}}{d}_{ij}\right)\right], \end{array}$$
where \(\omega\) is the angular frequency, \({A}_{\text{R}}\) is the power spectrum for the Rayleigh wave, \({B}_{\text{L}}\) is the power spectrum for the Love wave, \({J}_{n}\) is the nth-order Bessel function of the first kind, \({d}_{ij}\) is the distance between ith and jth channels, \({c}_{\text{R}}\) is the phase velocity of the Rayleigh wave, and \({c}_{\text{L}}\) is the phase velocity of the Love wave. In Eq. 3, the first and second terms represent the Rayleigh and Love waves, respectively. In general, the Love waves have less energy than the Rayleigh waves at a period range below 0.1 Hz (Nishida et al. 2008a). Furthermore, as mentioned by Nakahara et al. (2021), Rayleigh waves predominate in the far field. The envelope of the Rayleigh wave term decays with an order of \({d}^{-1/2}\), while the Love wave term decays with an order of \({d}^{-3/2}\) in the far field. For these reasons, and given the small estimated contribution of the Love wave terms, we obtained our misfit function \(E\) as:
$$\begin{array}{c}L2\left(\omega ,{c}_{\text{R}}\right)=\sum_{ij}{\left(Re\left[{S}_{ij}^{\text{obs}}\left(\omega \right)\right]-{S}_{ij}^{\text{syn.}}\left(\omega ,{c}_{\text{R}}\right)\right)}^{2}, \end{array}$$
$$\begin{array}{c}E\left({\beta }_{l}, {h}_{l}\right)=\frac{1}{\omega }\int L2\left(\omega ,{c}_{R}\left(\omega ;{\beta }_{l};{h}_{l}\right)\right)d\omega , \end{array}$$
where \({S}_{ij}^{\text{obs}}(\omega )\) and \({S}_{ij}^{\text{syn.}}\left(\omega \right)\) are the observed and synthetic cross-spectra between ith and jth channels calculated from Eq. 3, respectively. \({\beta }_{l}\) and \({h}_{l}\) are the Vs and thickness at the lth layer, respectively. We followed the methods described by Takeo et al. (2022) and Yoshizawa & Kennett (2002) to estimate the phase velocity model. \({\beta }_{l}\) and \({h}_{l}\) are model parameters for minimizing \(E\left({\beta }_{l}, {h}_{l}\right)\) in Eq. 5. In this study, we adopted the simulated annealing algorithm of Goffe et al. (1994) as a global optimizer in our search for model parameters. We avoided numerical instability by constraining the Vs at each layer at greater than 80% of the value of the layer directly above it.
In practice, we estimated the phase velocity and evaluated the errors in the estimated phase velocity. First, we estimated the dispersion curves of the phase velocity from the cross-spectra for all pairs of channels for each subarray. We then obtained a 1-D model with six layers, where the deepest layer had an infinite thickness for each subarray. We adopted the scaling relationship between Vp, Vs, and density described by Brocher (2005). Next, searching by trial and error, we divided the subarrays into two groups according to the distances from the coast and frequency range. The dispersion curves of the phase velocity were estimated using the only fundamental mode of the Rayleigh wave at the frequency range between 0.08 and 0.50 Hz for one group that had distances ranging from 28 to 52 km from the coast. For the other group with distances ranging from 52 to 75 km, the fundamental mode (0.08–0.50 Hz) and the first higher mode (0.25–0.50 Hz or 0.08–0.50 Hz) of Rayleigh wave were used to obtain a dispersion curve. Due to large incoherent noise near the coast and offshore areas, the spatial range of the present data processing was limited from 28 to 75 km from the coast. Synthetic phase velocities were calculated using DISPER80 (Saito 1988).
We also used the bootstrap method (Efron 1992) to estimate errors. We aggregated station pairs randomly selected from all station pairs, allowing for overlap for a bootstrap sample. We then calculated the standard deviation of 100 dispersion curves estimated from each of the 100 bootstrap samples. The Vs structures estimated from the data of all pairs of channels for each subarray were used as the reference model.
Inversion of the S-wave velocity structure
The phase velocity of the Rayleigh wave estimated by seismic interferometry was used to obtain a 1 – D Vs isotropic model. The 1-D Vs structure was determined by minimizing the misfit function E, as defined by:
$$\begin{array}{c}E = \sqrt {\frac{1}{N}\sum {\left[ {\frac{{c_{{obs}} \left( \omega \right) - c_{{syn}} \left( {\omega ;\beta _{l} ;h_{l} } \right)}}{{c_{{err}} \left( \omega \right)}}} \right]^{2} ,} } \end{array}$$
where N is the number of estimated phase velocities, and \({c}_{obs}\), \({c}_{syn}\) and \({c}_{err}\) indicate the observed, synthesized, and uncertainty of the phase velocity, respectively. We obtained a stable result by selecting only a frequency band with phase velocities with an error lower than 0.1 km/s for the inversion analysis. We again set the Vs (\({\beta }_{l})\) and thickness (\({h}_{l})\) as the model parameters to minimize the value obtained from Eq. 6. We then used the simulated annealing algorithm method to obtain the optimal model parameters. The bootstrap average models estimated from the phase velocity estimation were used as the reference model. We estimated errors using the bootstrap method once again. Here, we aggregated a final dispersion curve randomly selected from the 100 dispersion curves estimated in the previous section (Takeo et al. 2013).
The 1-D Vs structures were estimated on the assumption that we could ignore the effect of lateral heterogeneity for profiles perpendicular to the seafloor cable. Previous seismic reflection and refraction surveys support this assumption since profiles perpendicular to the seafloor cable show little lateral heterogeneity to a depth of 10 km (e. g. Takahashi et al. 2004).
We computed the CCFs between all possible pairs of each subarray, and the resulting CCFs clearly showed a surface wave (Fig. 2). The group velocities of about 0.5 km/s were much slower than the surface waves from land observations in the same period range. The surface wave with low velocity corresponded to a special type of Rayleigh wave or Scholte wave (Stokoe et al. 1991; Bagheri et al. 2015) because the Love wave had less energy than the Rayleigh wave in the period range below 0.1 Hz. The SNR was much higher for the CCFs with the FK filter than without the FK filter. We used the method of Bensen et al. (2007) to estimate the SNR.
Calculated CCFs. Calculated CCFs. Each trace was normalized by its amplitude. A and B—CCFs at a distance of 39.94 km from the coast without and with an FK filter, respectively. The frequency range was between 0.17 and 0.50 Hz. C and D—CCFs without and with an FK filter at 67.48 km from the coast in the frequency range between 0.10 and 0.25 Hz. SNR indicates the signal-to-noise ratios calculated by Bensen's (2007) method. The dashed lines indicate an apparent velocity of 0.5 km/s. The FK filter improves the SNR values and enables effective extraction of surface waves
We selected the dispersion curve using the CCFs (Figs. 3-A and 3-E). Extracting the Rayleigh waves by interferometry enabled a stable estimation of the dispersion curves with errors less than 0.1 km/s by bootstrapping. Near the coast (Fig. 3A), the fundamental Rayleigh wave was dominant in the frequency range between 0.2 and 0.5 Hz. By contrast, in the case of an offshore subarray (Fig. 3E), both the fundamental (0.5–0.1 Hz) and the 1st higher mode (0.5–0.2 Hz) Rayleigh waves predominated. This result clearly indicates that the thickness of the low-Vs layer differs between the near-coast and offshore areas. We determined the Vs structure using only the phase velocities of the fundamental mode of Rayleigh waves in the near-coast area at distances between 28 and 52 km. In the area where the distance was greater than 52 km from the coast, both the fundamental and first higher modes of the Rayleigh waves were used for phase velocity estimation.
Phase velocity, Vs structures, and sensitivity kernel. A and E—Normalize residuals and estimation phase velocity at subarrays from the coast of 39.94 and 67.68 km, respectively. The error bars represent the uncertainties estimated by the bootstrap method. The white dots and cross makers indicate errors less than and more than 0.1 km/s, respectively. We used phase velocity only with an error less than 0.1 km/s for our inversion analysis. B and F—Vs structures at subarrays from the coast of 39.94 and 67.68 km. D and H White dots and green lines indicate the observed and estimated phase velocities, respectively. The errors were estimated by the bootstrap method. C Normalized Vs sensitivity kernel of the fundamental Rayleigh wave at 39.94 km from the coast. G and I Normalized Vs sensitivity kernels of fundamental and first higher mode Rayleigh waves at 67.68 km from the coast
We inverted the measured phase velocities into 1-D Vs structures (Fig. 3) below each subarray and calculated the 1-D Vs structures (Fig. 3B and F) and normalized Vs sensitivity kernels (Fig. 3C, G, and I) of the Rayleigh waves. Errors were calculated using the bootstrap method. A small error was estimated in the depth range where Vs values were slower than 2.0 km/s. The synthetic phase velocities calculated by DISPER80 using the 1-D Vs inversion results were consistent with the observations (Figs. 3D and H). Given that the sensitivity kernel had high sensitivity in the depth range of Vs slower than 2.0 km/s, our inversion provided a Vs structure with small errors at that depth range of Vs. By contrast, the layers of Vs greater than 2.0 km/s were not estimated well. However, the Vs sharply increases to those greater than 2.0 km/s at the boundary between the shallow slow layers and the layers below, and the mutual velocity ranges, with errors, are clearly separated. Therefore, the depth was reliably estimated for the discontinuity between the shallow slow layers and the layers below.
Finally, a 2-D Vs profile (Fig. 4) was obtained from the 1-D Vs structures at each subarray. Although each 1-D Vs structure at a given subarray was composed of six layers, we classified this 2-D profile into four units with a clear velocity discontinuity (Fig. 4). The velocities of some contiguous layers were coincident in terms of errors. In those cases, we unified the contiguous layers into one unit. The Vs of Unit-1 at depths of 0.7–2.6 km and Unit-2 at the depth of 0.9–4.0 km were slower than 0.6 km/s and approximately 0.7 km/s, respectively. The Vs of Unit-3 at depths of 1.2–6.8 km and Unit-4 layers below the Unit-3 were 1.1–1.8 km/s and greater than 2.0 km/s, respectively.
2-D Vs structure estimated from phase velocity information. 2-D Vs structure estimated from phase velocity information. The 2-D structure model consists of four units. Unit annotations are the average Vp/Vs value with errors. The errors were defined as double standard deviation during calculation of average. The white circles indicate the aftershocks off the Pacific coast of the Tohoku Earthquake determined by pop-up OBSs (Shinohara et al., 2012). The aftershocks within a box of 40 km in width on both sides along the fiber cable are projected onto this profile. Some aftershocks occurred in Unit4, whereas there was no seismicity observed in the sedimentary basin
We compared our 2-D Vs structure (Fig. 4) with the Vs and Vp structures obtained in other studies to interpret these units. The Vs values of Unit-1 and Unit-2 layers, at less than 0.6 km/s and about 0.7 km/s, respectively, were consistent with the results presented by Spica et al. (2020) with the same profile. Takahashi et al. (2004) provided the Vp structure estimated by refraction surveys along a profile parallel to our DAS cable, positioned northward, with a distance of about 20 km from the seafloor cable (Fig. 1). The Vp values of Unit-1 and Unit-2 in the present study were determined as 1.7 and 2.3 km/s, respectively, and these layers were interpreted as Neogene sediment by Spica et al. (2020). The thicknesses of both Unit-1 and Unit-2 increased seaward for a distance range from 28 to 47 km. This increase in thickness is also consistent with the results reported by previous studies. In distance range further from the coast between 47 and 75 km, however, the thicknesses of Unit-1 and Unit-2 showed very small variations.
The thickness of Unit-3 was 0.6 km at a distance of 28 km from the coast and increased with a horizontal distance. Unit-3 became the thickest (4.1 km) at a distance from 52 km from the coast. The thickness of Unit-3 decreased to 1.9 km at a distance of 65 km. At distances between 65 and 75 km, Unit-3 appeared to pinch out. These strong lateral heterogeneities in Unit-3 were also evident in the Vp model by Takahashi et al. (2004). The discontinuity observed between Unit-3 and Unit-4 can be interpreted as the boundary between sediments and the island arc uppermost crust consistent with the interpretation reported by Takahashi et al. (2004). Along the same Japan trench offshore Ibaraki and the south of the study area, Yamaya et al. (2021) obtained Vs values of the deepest sediment of approximately 2.1 km/s using records from pop-up ocean bottom seismometers. In contrast, Unit-3 in this study has smaller Vs compared to the results off Ibaraki. We attribute the difference in Vs values for the two regions to differences in porosities and lithologies in the lowermost sedimentary layers.
Although our Vs profile and the Vp profile of Takahashi et al. (2004) are 20 km apart (Fig. 1) and the spatial resolution of our results is much higher than that of Takahashi et al. (2004), we calculated the average value and errors of the Vp/Vs (Fig. 4) using the Vp estimated by Takahashi et al. (2004) to support our interpretation of the results. The errors were defined as twice the standard deviation during calculation of the average. The average Vp/Vs values of Unit-1 and Unit-2 were 5.00 and 3.34, respectively. The Vs in Unit-3 of 1.1–1.8 km/s was obtained from this study and the Vp in Unit-3 was estimated to be 4.2 km/s by Takahashi et al. (2004). As a result, the average Vp/Vs value in Unit-3 become 3.12. On the other hand, the Vp/Vs of Unit-4 was 1.98 from a Vs of > 2.0 km/s and a Vp of 5.0 km/s.
Because we found no previous reports on the Vp/Vs of marine sedimentary layers in our study area, we compared Vp/Vs with the values obtained in other marine areas for the similar geological ages. Vp/Vs values of sedimentary layers are known to vary depending on porosity of rocks, constitutive materials, and tectonic settings from place to place (e.g., Kodaira et al. 1996). We chose the ages of sediments as an indicator of the comparison rather than the mineral components, because porosity of medium primarily affects Vs of sediments and uppermost crust under marine environment. A comparison of the thickness of each Unit in our results to those of the layers obtained by Takahashi et al. (2004) indicates that Unit-1 and Unit-2 are Neogene sediments, and Unit-3 is Cretaceous sediments. The Vp/Vs of the layer deposited in Neogene was estimated at 5.10 in Norway (Kvarven et al. 2015) and 3.2–4.4 in the landward slope on the southern Japan trench (Yamaya et al. 2021). These values are comparable to our estimation of Vp/Vs in Unit-1 and Unit-2. Kodaira et al. (1996) reported that the Vp/Vs of the sedimentary layers in Norway ranged from 5 to 3, although their ages were unknown. By contrast, Kvarven et al. (2015) reported a Vp/Vs of 1.78 for the Cretaceous sediment layer in Norway. On the other hand, Yamaya et al. (2021) showed a Vp/Vs of 2.1 for the Cretaceous sediment layer beneath the landward slope of the southern Japan trench. Although our Vp/Vs values are believed to have some uncertainty due to lower spatial resolution of Vp results from the airgun-OBS surveys, the Vp/Vs for the Cretaceous sediment unit in our study region can be concluded to have large value obviously from a point of view of the errors estimated in this study. We therefore infer that the Cretaceous sediments in this study area is unconsolidated.
Now let us consider our result on a deep part or the crust. Using stationary seismic networks in land and sea by body-wave travel-time tomography, Matsubara et al. (2019) reported Vp/Vs of 1.6–1.9 for the upper crust below the northeastern Japan island arc. In our study, the Vp/Vs for Unit-4 seems to correspond to the uppermost part of the island arc crust from the point of view of their Vp/Vs Furthermore, aftershocks of the 2011 Tohoku-oki earthquake, as estimated by Shinohara et al. (2012), were located only in Unit-4 (Fig. 4), whereas no seismicity was experienced in the other three units. In summary, Unit-1, Unit-2, and Unit-3 consist of sedimentary layers, whereas Unit-4 represents the uppermost part of the igneous crust of an island arc.
Previous studies that estimated Vs structures using DAS records have been limited to elucidation of a shallow part of the sediment layers. For example, Spica et al. (2020) estimated the Vs structure to 3 km in depth. The application of the seismic interferometry with FK filtering to the DAS data now enabled the determination of phase velocities of surface waves whose periods are longer than those of previous studies. Using DAS data with high spatial resolution we were able to obtain an S-wave structure of lowermost sediment > 3.0 km depth and extending to 6.8 km depth where sediments were thickest.
The SPAC method is useful to estimate velocity structure using data from a seismic network; however, spatial aliasing due to a sparse distribution of stations causes uncertainty for the SPAC analysis. Because DAS measurement has a high-spatial density, spatial aliasing can be avoided. For example, we could prevent misidentification of a mode branch. Vs structures in marine area have conventionally been estimated by the observation of converted waves during seismic surveys or by seismic interferometry using spatially dense OBSs data. Our method offers a new approach for determining the Vs of shallow structures with high spatial resolution in marine areas using short-term DAS data and a seafloor cable. A shallow Vs structure in a marine area with high resolution by our method leads to a detailed distribution of Vp/Vs in a shallow region, which is useful for considering tectonics and rock properties.
The information of Vs is important for understanding the rock properties of the upper crust and sediments on the seafloor. We applied a seismic interferometry method to DAS records obtained during a recording period of 13 h by a seafloor cable installed off Sanriku, Japan. We found that applying the FK filter to the DAS data before calculating the CCFs effectively enhances the surface waves. The phase velocities of the Rayleigh waves were calculated by the SPAC method. Finally, we inverted the Vs structure from the phase velocities of the Rayleigh waves obtained from seismic interferometry. The Vs structures of the sediment layers and the upper crust were consistent with those determined by refraction surveys. The seismic interferometry method with DAS data collected through seafloor cables could estimate Vs structures with high spatial resolution.
The Vs structure of sediment layers has been widely explored using seismic interferometry and receiver function methods with OBS data. However, their spatial resolution is constrained by the spatial density of the installed OBSs. By contrast, a Vs structure with high spatial resolution can be estimated from short-term DAS records by applying the seismic interferometry method. Our results identified a new approach for estimating the heterogeneous Vs structure of sediment layers and the upper crust in subduction zones.
The DAS observations were performed as part of the Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research) by the Ministry of Education, Culture, Sports, Science and Technology of Japan. The cross-spectra data in this study are available at https://doi.org/10.5281/zenodo.5904941. The raw data supporting the conclusions of this study are available from the corresponding authors.
V S :
S-wave velocity
V P :
P-wave velocity
Vp/Vs:
Ratio of P- and S-wave velocities
Distributed acoustic sensing
FK:
Frequency–wavenumber
CCFs:
Cross-correlation functions
SNRs:
Signal-to-noise ratios
SPAC:
Spatial auto-correlation
OBSs:
Ocean bottom seismometers
Aki K (1957) Space and time spectra of stationary stochastic waves, with special reference to microtremors. Bull Earthq Res Inst 35:415–456
Akuhara T, Tsuji T, Tonegawa T (2020) Overpressured underthrust sediment in the Nankai Trough forearc inferred from transdimensional inversion of high-frequency teleseismic waveforms. Geophys Res Lett. https://doi.org/10.1029/2020GL088280
Atterholt J, Zhan Z, Shen Z, Li Z (2022) A unified wavefield-partitioning approach for distributed acoustic sensing. Geophys J Int 228:1410–1418. https://doi.org/10.1093/gji/ggab407
Ayres A, Theilen F (1999) Relationship between P-and S-wave velocities and geological properties of near-surface sediments of the continental slope of the Barents Sea. Geophys Prospect 47:431–441. https://doi.org/10.1046/j.1365-2478.1999.00129.x
Bagheri A, Greenhalgh S, Khojasteh A, Rahimian M (2015) Dispersion of Rayleigh, Scholte, Stoneley and Love waves in a model consisting of a liquid layer overlying a two-layer transversely isotropic solid medium. Geophys J Int 203:195–212. https://doi.org/10.1093/gji/ggv278
Bensen GD, Ritzwoller MH, Barmin MP, Levshin AL, Lin F, Moschetti MP, Shapiro NM, Yang Y (2007) Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys J Int 169:1239–1260. https://doi.org/10.1111/j.1365-246X.2007.03374.x
Brocher TM (2005) Empirical relations between elastic wavespeeds and density in the Earth's crust. Bulletin Seismol Soc Am 95:2081–2092. https://doi.org/10.1785/0120050077
Cedilnik G, Lees G, Schmidt, P, Herstrøm S, Geisler T (2019) Ultra-long reach fiber distributed acoustic sensing for power cable monitoring. In: Proceedings of the JICABLE'19 10th International Conference on Power Insulated Cables, Versailles, France, 23–27 June 2019. Available via DIALOG. https://www.apsensing.com/fileadmin/Publication%20Files/Cedilnik_et_al._-_2019-Jicable-Proceedings-E4-4-Ultralong-reach-DAS-.pdf. Accessed 17 May 2022.
Dou S, Lindsey N, Wagner AM, Daley TM, Freifeld B, Robertson M, Peterson J, Ulrich C, Martin ER, Ajo-Franklin JB (2017) Distributed acoustic sensing for seismic monitoring of the near surface: a traffic-noise interferometry case study. Sci Rep 7:11620. https://doi.org/10.1038/s41598-017-11986-4
Efron B (1992) Bootstrap methods: another look at the jackknife. Breakthroughs in statistics. Springer, Berlin, pp 569–593. https://doi.org/10.1007/978-1-4612-4380-9_41
Fujie G, Kodaira S, Kaiho Y, Yamamoto Y, Takahashi T, Miura S, Yamada T (2018) Controlling factor of incoming plate hydration at the north-western Pacific margin. Nat Commun 9:3844. https://doi.org/10.1038/s41467-018-06320-z
Goffe WL, Ferrier GD, Rogers J (1994) Global optimization of statistical functions with simulated annealing. J Econom 60:65–99. https://doi.org/10.1016/0304-4076(94)90038-8
Hudson TS, Baird AF, Kendall JM, Kufner SK, Brisbourne AM, Smith AM, Butcher A, Chalari A, Clarke A (2021) Distributed acoustic sensing (DAS) for natural microseismicity studies: a case study from Antarctica. J Geophys Res Solid Earth. https://doi.org/10.1029/2020jb021493
Hunter JD (2007) Matplotlib: a 2D graphics environment. Computing Sci Eng 9:90–95. https://doi.org/10.1109/MCSE.2007.55
Kodaira S, Bellenberg M, Iwasaki T, Kanazawa T, Hirschleber HB, Shimamura H (1996) Vp/Vs ratio structure of the Lofoten continental margin, northern Norway, and its geological implications. Geophys J Int 124:724–740. https://doi.org/10.1111/j.1365-246X.1996.tb05634.x
Kodaira S, Iidaka T, Kato A, Park J-O, Iwasaki T, Kaneda Y (2004) High pore fluid pressure may cause silent slip in the Nankai Trough. Science 304:1295–1298. https://doi.org/10.1126/science.1096535
Kvarven T, Mjelde R, Hjelstuen BO, Faleide JI, Thybo H, Flueh ER, Murai Y (2015) Crustal composition of the Møre margin and compilation of a conjugate Atlantic margin transect. Tectonophysics 666:144–157. https://doi.org/10.1016/j.tecto.2015.11.002
Matsubara M, Sato H, Uehira K, Mochizuki M, Kanazawa T, Takahashi N, Suzuki K, Kamiya S (2019) Seismic velocity structure in and around the Japanese Island arc derived from seismic tomography including NIED MOWLAS Hi-net and S-net data. In: Kanao M, Toyokuni G (eds) Seismic waves - probing earth system. London: IntechOpen. https://doi.org/10.5772/intechopen.86936
Miura S, Kodaira S, Nakanishi A, Tsuru T, Takahashi N, Hirata N, Kaneda Y (2003) Structural characteristics controlling the seismicity of the southern Japan Trench fore-arc region, revealed by ocean bottom seismographic data. Tectonophysics 363:79–102. https://doi.org/10.1016/S0040-1951(02)00655-8
Nakahara H, Emoto K, Nishimura T (2021) Extending the formulation of the spatial autocorrelation (SPAC) method to strain, rotation and tilt. Geophys J Int 227:287–302. https://doi.org/10.1093/gji/ggab217
Nishida K, Kawakatsu H, Fukao Y, Obara K (2008a) Background Love and Rayleigh waves simultaneously generated at the Pacific Ocean floors. Geophys Res Lett. https://doi.org/10.1029/2008gl034753
Nishida K, Kawakatsu H, Obara K (2008b) Three-dimensional crustal S wave velocity structure in Japan using microseismic data recorded by Hi-net tiltmeters. J Geophys Res Solid Earth. https://doi.org/10.1029/2007jb005395
Okada H (2006) Theory of efficient array observations of microtremors with special reference to the SPAC method. Explor Geophys 37:73–85. https://doi.org/10.1071/EG06073
Tribaldos VR, Ajo-Franklin JB (2021) Aquifer monitoring using ambient seismic noise recorded with distributed acoustic sensing (DAS) deployed on dark fiber. J Geophys Res Solid Earth. https://doi.org/10.1029/2020jb021004
Saito M. (1988). DISPER80: A subroutine package for the calculation of seismic normal-mode solutions, in Seismological Algorithms, 293-319, Academic, San Diego, Calif. Matsubara et al. (2019) is an online paper at https://www.intechopen.com/chapters/67965
Shinohara M, Machida Y, Yamada T, Nakahigashi K, Shinbo T, Mochizuki K, Murai Y, Hino R, Ito Y, Sato T, Shiobara H, Uehira K, Yakiwara H, Obana K, Takahashi N, Kodaira S, Hirata K, Tsushima H, Iwasaki T (2012) Precise aftershock distribution of the 2011 off the Pacific coast of Tohoku Earthquake revealed by an ocean-bottom seismometer network. Earth Plan Sp 64:8. https://doi.org/10.5047/eps.2012.09.003
Shinohara M, Yamada T, Akuhara T, Mochizuki K, Sakai S (2022) Performance of seismic observation by distributed acoustic sensing technology using a seafloor cable off Sanriku. Front Mar Sci. https://doi.org/10.3389/fmars.2022.844506
Shinohara M, Yamada T, Akuhara T, Mochizuki K, Sakai S, Hamakawa M, Kasajima T, Arioka T, Kubota S (2019) Distributed Acoustic Sensing measurement by using seafloor optical fiber cable system off Sanriku for seismic observation. OCEANS 2019 MTS/IEEE. Seattle, WA, pp 1–4. https://doi.org/10.23919/oceans40490.2019.8962757
Spica ZJ, Nishida K, Akuhara T, Pétrélis F, Shinohara M, Yamada T (2020) Marine sediment characterized by ocean-bottom fiber-optic seismology. Geophys Res Lett 47:e2020GL088360. https://doi.org/10.1029/2020GL088360
Stokoe KH, Gauer RC, Bay JA (1991) Experimental investigation of seismic surface waves in the seafloor. In: Hovem JM, Richardson MD, Stoll RD (eds) Shear Waves in Marine Sediments. Springer, Dordrecht, pp 51–58. https://doi.org/10.1007/978-94-011-3568-9_6
Takagi R, Toyokuni G, Chikasada N (2021) Ambient noise correlation analysis of S-net records: extracting surface wave signals below instrument noise levels. Geophys J Int 224:1640–1657. https://doi.org/10.1093/gji/ggaa548
Takahashi N, Kodaira S, Tsuru T, Park J-O, Kaneda Y, Suyehiro K, Kinoshita H, Abe S, Nishino M, Hino R (2004) Seismic structure and seismogenesis off Sanriku region, northeastern Japan. Geophys J Int 159:129–145. https://doi.org/10.1111/j.1365-246X.2004.02350.x
Takeo A, Forsyth DW, Weeraratne DS, Nishida K (2014) Estimation of azimuthal anisotropy in the NW Pacific from seismic ambient noise in seafloor records. Geophys J Int 199:11–22. https://doi.org/10.1093/gji/ggu240
Takeo A, Nishida K, Isse T, Kawakatsu H, Shiobara H, Sugioka H, Kanazawa T (2013) Radially anisotropic structure beneath the Shikoku Basin from broadband surface wave analysis of ocean bottom seismometer records. J Geophys Res: Solid Earth 118:2878–2892. https://doi.org/10.1002/jgrb.50219
Takeo A, Nishida K, Aoyama H, Ishise M, Kai T, Kurihara R, Maeda T, Mizutani Y, Nakashima Y, Nagahara S, Wang X, Ye L, Akuhara T, Aoki Y (2022) S-wave modelling of the Showa-Shinzan lava dome in Usu Volcano, Northern Japan, from seismic observations. Geophys J Int. https://doi.org/10.1093/gji/ggac111
Wessel P, Luis JF, Uieda L, Scharroo R, Wobbe F, Smith WHF, Tian D (2019) The generic mapping tools version 6. Geochem Geophys Geosyst 20:5556–5564. https://doi.org/10.1029/2019gc008515
Yamaya L, Mochizuki K, Akuhara T, Nishida K (2021) structure derived from multi-mode ambient noise tomography with dense OBS network at the Japan Trench. J Geophys Res Solid Earth 126:e2021JB021789. https://doi.org/10.1029/2021JB021789
Yao H, Gouédard P, Collins JA, McGuire JJ, van der Hilst RD (2011) Structure of young East Pacific Rise lithosphere from ambient noise correlation analysis of fundamental- and higher-mode Scholte-Rayleigh waves. C R Geosci 343:571–583. https://doi.org/10.1016/j.crte.2011.04.004
Yoshizawa K, Kennett BLN (2002) Non-linear waveform inversion for surface waves with a neighbourhood algorithm—application to multimode dispersion measurements. Geophys J Int 149:118–133. https://doi.org/10.1046/j.1365-246X.2002.01634.x
Zhan Z (2020) Distributed acoustic sensing turns fiber-optic cables into sensitive seismic antennas. Seismol Res Lett 91:1–15. https://doi.org/10.1785/0220190112
The authors thank Drs. M. Masuda, S. Tanaka and Messrs. T. Hashimoto, K. Miyakawa, T. Yagi of the Earthquake Research Institute, the University of Tokyo for technical support in the DAS observations. The collaboration with Fujitsu Laboratories Ltd. in the field of data acquisition is appreciated. We also acknowledge discussions we had with Drs. T.Akuhara and S.Takemura for this study. Some figures were made using the Matplotlib (Hunter 2007) and Generic Mapping Tools (Wessel et al. 2019).
This study was mainly supported by the Ministry of Education, Culture, Sports, Science, and Technology of Japan under the Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research). A part of this study was also funded by the Earthquake Research Institute, the University of Tokyo.
Department of Earth and Planetary Science, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo‑ku, Tokyo, 113-0033, Japan
Shun Fukushima
Earthquake Research Institute, The University of Tokyo, 1‑1‑1, Yayoi, Bunkyo‑ku, Tokyo, 113‑0032, Japan
Shun Fukushima, Masanao Shinohara, Kiwamu Nishida, Akiko Takeo & Tomoaki Yamada
Department of Natural History Science, Graduate School of Science, Hokkaido University, North 10 West 8 Kita-ku, Sapporo, Hokkaido, 060-0810, Japan
Kiyoshi Yomogida
Masanao Shinohara
Kiwamu Nishida
Akiko Takeo
Tomoaki Yamada
SF played a leading role in this study, including data processing, analysis, and completion of the manuscript. MS and TY led the temporal DAS observation and interpreted the results. KN, AT, and KY contributed to develop the analysis method and interpreted the results. All authors read and approved the final manuscript.
Correspondence to Shun Fukushima.
The authors declare that they have no competing interest.
Effectiveness of frequency-wave number (FK) filter for calculation of CCFs. The CCFs with FK filter ((B), (D) and (F)) are compared those without FK filter ((A), (C) and (D)). The CCFs at distances of 39.94 km, 54.95 km, and 64.97 km are shown in pairs of (A) and (B), (C) and (D), and (E) and (F), respectively. Each trace was filtered from 0.1 Hz to 0.2 Hz and was normalized by its amplitude. Seismic signals propagating with an apparent velocity range from 0.35 km/s to 4.0 km/s were selected by FK filtering. For the estimation of signal-to-noise ratio calculated by root mean square amplitude ratio between surface wave and noise data for the last 50s, the method of Bensen et al. (2007) was used. DAS instrument generally has a coherent noise has at the same timing in all stations (zero-lag noise) and the coherent noise interfere analyses for surfaces waves. We can remove the zero-lag noise by using FK filter due to difference of apparent velocity. Applying FK filter to DAS data before calculating CCFs effectively enhances surface waves in DAS records. Figure S2. Calculated cross correlation functions (CCFs)with a frequency-wave number filterat the frequency range between 0.08 and 0.50 Hz. Each trace was normalized by its amplitude. (A), (B), and (C) denoteCCFs at distancesof 29.93, 32.43, and 64.97km from the coast, respectively.The dashed lines indicate anapparent velocity of 0.5 km/s.Note that the seismic waves propagating with velocities of approximately 0.5 km/s are recognized in the CCFs. The group velocities of about 0.5 km/s were much slower than the surface waves in land for an identical period range. The surface wave propagating with a low velocity is estimated to correspond to a special type of Rayleigh wave or Scholte wave. Figure S3.:Normalizedresiduals and estimationofphase velocitiesin thesubarrays from the coast of29.93(A), 32.43(B), and 64.97(C).Bars indicate phase-velocity uncertainties, which wereestimated by a bootstrap method. The estimated phase velocitieswith errors less than0.1 km/s are indicated by white circles. Although the fundamental Rayleigh wave was dominant near the coast, both fundamental (0.5–0.1 Hz) and 1st higher mode (0.5–0.2 Hz) Rayleigh waves appeared far from the coast. The faster phase velocities were observed in areas closer to the coast. Figure S4. The obtained Vs structures at subarrays centered at distances at 29.93 km (A), 32.43 km (B), and 64.97 km (C) from the coast. The errors were estimated by the bootstrap method. Vs and thickness in each layer are set to be model parameters and minimized values of the misfit function between observed and synthetic dispersion curves in Equation 6
Fukushima, S., Shinohara, M., Nishida, K. et al. Detailed S-wave velocity structure of sediment and crust off Sanriku, Japan by a new analysis method for distributed acoustic sensing data using a seafloor cable and seismic interferometry. Earth Planets Space 74, 92 (2022). https://doi.org/10.1186/s40623-022-01652-z
S-wave velocity structure
Sediment and upper crust
Seismic interferometry
Seafloor cable
Frequency-wavenumber filter