Estimation of source, path, and site factors of S waves recorded at the S-net sites in the Japan Trench area using the spectral inversion technique

S-net is a large-scale ocean bottom (OB) network in the Japan Trench area, consisting of inline-type 150 observatories equipped with seismometers and pressure gauges. Among them, 41 observatories have been buried about one meter beneath the seafloors in the shallow water regions (water depth <1500 m). We analyzed the strong-motion data recorded at the S-net sites from earthquakes with magnitudes 3.5 < Mw ≤ 7 and focal depths < 70 km to understand the site amplification effect on the recorded motions. We used the spectral inversion technique and obtained some fundamental properties of the earthquake source spectra, path attenuation, and site factors from the horizontal-component S-wave portions of the recordings. We obtained that the source spectra followed the ω-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }^{-2}$$\end{document} source model generally well, and the estimated magnitudes were mostly within ± 0.3 magnitude units of the catalog magnitudes. The stress drops increased systematically with the focal depths, and the values for the shallow earthquakes in the Pacific Plate were higher than those for the interplate earthquakes with comparable focal depths. The path-averaged quality factors were generally frequency-dependent and were somewhat larger than those in the past studies. The peak site frequencies ranged between about 0.2 and 10 Hz, while the peak amplification factors ranged between 10 and 50. Even though the peak frequencies and amplification factors differed from site to site, the peak frequencies were mostly higher than 2 Hz at the outer trench stations, while they were lower than 2 Hz at many inner trench stations. The amplification factors at a few OB sites in the shallow water regions were comparable with the theoretical ones computed from the 1-D subsea model. The amplification factors at intermediate frequencies (~ 0.3 to 2 Hz) generally increased with P-wave travel time in the sediments estimated from the multi-channel seismic survey. At about 20% of the sites (mainly at the unburied stations), spurious site spectra were recognized at frequencies higher than about 4 Hz. If the site spectra between 4 and 10 Hz are required, using only the X-component records is recommended.


Introduction
S-net is a large-scale ocean-bottom network for earthquakes and tsunamis in the Japan Trench area established after the 2011 Tohoku-oki earthquake disaster. One of the main objectives of the establishment of the S-net was to enhance the Japan Meteorological Agency (JMA) earthquake early warning (EEW) and tsunami early warning (Aoi et al. 2020). The network consists of inlinetype 150 observatories divided into six segments, namely, S1, S2, S3, S4, S5, and S6 (Fig. 1). It transmits waveform data to the data center of the National Research Institute for Earth Science and Disaster Resilience (NIED) continuously. The observatories and cables in shallow water regions (water depth < 1500 m) have been buried one meter beneath the seafloors to protect them from fishing activities and increase ground coupling except for three stations, S1N08, S1N12, and S6N24. The total numbers of buried and unburied stations are 41 and 109, respectively.
Several previous studies reported that the average amplitudes of the long-period ground motions between about 1 and 20 s were larger at the ocean bottom seismograph sites in the Nankai Trough area (e.g., Nakamura et al. 2015;Kubo et al. 2019) and the Japan Trench area (Dhakal et al. 2021) compared to the average amplitudes at the strong-motion stations on land. Hayashimoto et al. (2019) constructed an equation to estimate the magnitude of the suboceanic earthquakes for EEW using the vertical component records at the seafloor sites to avoid the large amplification of the horizontalcomponent ground motions due to sediments and also to minimize the effect of the rotational noises during strong shakings. These studies have indicated that understanding the site amplification effects on the recorded motions is essential for the reliable estimation of magnitudes and ground motions for an EEW. Therefore, one of the main objectives of this paper is to evaluate the S-wave site amplification factors at the S-net sites in the Japan Trench area based on the recorded strong-motion data.
In the present study, we used the spectral inversion technique (e.g., Andrews 1986;Iwata and Irikura 1988), described at some length in a later section, to simultaneously separate the source, path, and site factors of S waves recorded at the S-net stations. Therefore, the analysis of the source and path factors was an integral part of the present study to validate them and infer the site amplification factors. It also must be emphasized that the study of source spectra of earthquakes is crucial to understand the source processes of earthquakes and predict strong-motions from future earthquakes (e.g., Aki 1967;Hanks 1979;Morikawa and Sasatani 2003;Allmann and Shearer 2009). In addition, as the seismic waves in the present study were observed in wide offshore area, it was of interest to examine the path-averaged quality factor of S waves in the offshore region as most past studies were based on recorded motions on land (e.g., Nakamura et al. 2006;Nakano et al. 2015).
In the next two sections, we describe the data and methods, respectively. Then, we present and discuss the source spectra of S waves. We present estimated magnitudes of the earthquakes and compare them with the catalog magnitudes. We also derive stress drops of the earthquakes and discuss them. Following the discussion of the source spectra, the path-averaged quality factors for S waves are presented and compared with several previous studies. Finally, we present and discuss the site amplification factors at the S-net sites in some details. It is expected that the results presented in this study provide a basis for more detailed studies in the future about developing ground-motion prediction models for EEW and other seismological applications.  (S-net). The six segments of the S-net (S1-S6) are denoted by different polygons, with station codes at the start and end of each segment. The short-dashed rectangle indicates the location of sites, namely, S3N26, S3N25, S3N24, S3N23, S3N22, S3N21, S3N20, S3N19, and S6N12 from left to right, respectively, discussed later. The NAM, PHS, and ST are shorthand for the North-American Plate, Philippine Sea Plate, and Sagami Trough, respectively. Right panel: location of the epicenters of earthquakes used in the present study. The depth to the Pacific Plate is marked with contours (thin dashed lines) at interval of 10 km. The 10 km contour coincides roughly with the Trench axes. The plate-depth data were taken from Hirose (2022)

Data
We prepared the initial earthquake data set based on the moment tensor catalog of F-net, NIED. Then, we obtained 10 min of continuously recorded accelerograms at the S-net stations, beginning from 1 min before the earthquake origin time, from more than 1500 earthquakes, which occurred between February 2016 and October 2021. The moment magnitudes (M w ) of the earthquakes ranged between 3.8 and 7.0, and focal depths were shallower than 70 km. S-net sensors are housed inside cylindrical pressure vessels and record waveforms in three mutually perpendicular directions. However, the sensor axes are not aligned in the horizontal and vertical directions. In this study, conversions of the waveform data from the sensor axes to the horizontal and vertical directions were carried out following the procedures in Takagi et al. (2019). The present study used waveform data from two horizontal directions, namely, X and Y. The X direction coincides with the direction of the long axis of the cylindrical pressure vessel, and the Y direction is perpendicular to the X direction.
In the present study, all records were processed uniformly. The mean of a 1-min pre-event noise window was subtracted at each time step for the S-net records. In this paper, the mean or average implies the arithmetic mean of the data under consideration. For the landstation records, the mean of 10 s pre-event noise window was subtracted if the noise window was available; otherwise, the mean of the whole record was subtracted. Then, fourth-order low-cut filtering was applied to suppress low-frequency noises at 0.07 Hz. The S-wave arrival times were picked manually. After several trial analyses, the following magnitude-dependent S-wave time windows were selected: 8 s for M w < 4.5, 12 s for 4.5 ≤ M w < 5, 16 s for 5 ≤ M w < 5.5, and 20 s for M w ≥ 5.5. The one second time windows just before and after the defined time windows for the S waves were cosine-tapered. Then, zeroes were padded in the rear end, such that the window length remained 40.96 s for all recordings. The equal length allows one to evaluate the Fourier spectra at identical frequency points between all recordings. To evaluate the signal-to-noise ratios (SNRs), noise time histories of equal lengths with the corresponding S-wave time histories were selected from the pre-event parts of the records, and cosine tapering and zero paddings were applied to the noise time histories similar to that for the S-wave parts as explained above. Then, Fourier spectra were computed and smoothed using a Parzen window of 0.2 Hz for both S-wave and noise parts of equal durations of 40.96 s. Finally, the SNRs were obtained as the ratios of mean values of the smoothed spectra of two horizontal components between the S-wave and noise parts.
Final records were selected after imposing the following criteria. The vector peak ground accelerations (PGAs) of the original records were between 5 and 50 cm/s 2 . The maximum PGA value of 50 cm/s 2 was selected to avoid nonlinear site response and minimize the effect of the rotational noises in the records (e.g., Dhakal et al. 2017;Nakamura and Hayashimoto 2019;Takagi et al. 2019). The focal depths were shallower than 70 km to avoid the complicated attenuation effect for the deeper earthquakes. The epicentral distances were between 20 and 200 km. The choice of the distance limits was arbitrary, but it was expected that the minimum distance of about 20 km reduced the effect of the location errors for the records at small epicentral distances. The records were selected if the SNRs were larger than three at all frequencies between 0.1 and 20 Hz. At least three records exceeding the above threshold SNRs were available at each station and for each earthquake.
A preliminary analysis of the focal and centroid moment tensor depths estimated by the JMA showed that the two depth parameters differed largely for many events, even though the magnitudes of the events were smaller than about 5.5. The F-net, NIED, moment tensor solution obtains the depths of the events keeping the epicentral locations the same as that of the JMA epicentral locations. We found that the JMA centroid depths and NIED moment tensor depths were similar. In this study, we used the epicentral location from the JMA and depth from the F-net moment tensor solution to calculate the source-to-site distance used in the spectral inversion described in the next section. Hereafter, unless stated otherwise, the term 'distance' is used to mean the sourceto-site distance, and the 'depth' in the sense of depth obtained in the F-net NIED moment tensor solution. The location of the S-net stations and the epicenters of the selected events used in this study are depicted in the left and right panels of Fig. 1, respectively. The final data set contained 6326 recordings from 605 earthquakes. The general distributions of the data set in terms of the magnitude and distance, magnitude and depth, S-wave arrival time and distance, and the peak ground acceleration and distance are shown in Fig. 2a-d, respectively.

Methods
Spectral inversion method is a simple yet powerful technique for estimation of either site amplifications, path attenuation, and source parameters of earthquakes or their combinations and has been used in several previous studies (e.g., Andrews 1986;Iwata and Irikura 1988;Castro et al. 1990;Kato et al. 1992;Yamanaka et al. 1998;Satoh and Tatsumi 2002;Bindi et al. 2004;Tsuda et al. Dhakal et al. Earth, Planets and Space (2023) 75:1 2006;Oth et al. 2011;Ren et al. 2013;Nakano et al. 2015;Klimasewski et al. 2019;Fletcher and Boatwright 2020). In the present study, the observed S-wave acceleration Fourier spectra were expressed in the form shown in the following equation: where O ij (f ) is the observed S-wave Fourier amplitude spectrum at the jth station for the ith earthquake at frequency f , S i f is the source spectral amplitude of the ith earthquake at frequency f , G j f is the site amplification factor at the jth station at frequency f , Q s f is the average quality factor for S wave along the propagation path at frequency f , R ij is the distance between the jth site and ith event, and V s is the average S-wave velocity along the propagation path.
In the present study, the Fourier spectral amplitudes were computed as the mean values of the two horizontal components, as described in the previous section. The spectral amplitudes between 3 and 20 Hz were resampled with interval of approximately 0.1 Hz to save the computational time, and the results were obtained at 294 frequencies between ~ 0.0732 and 20 Hz. Equation (1) was linearized by taking base-10 logarithms and was solved by a least-square method after rearranging the terms and using the constrained conditions explained in the next paragraph (see Additional file 1 for details). The least-square solution obtained in this study may be expressed, as shown in the following equation (e.g., Searle 1971): where x is a solution set of unknown parameters (the source, path, and site factors), A is the design matrix of known values for the source, path, and site terms, and b is the set of observations for a given frequency. The T and −1 in Eq. (2) indicate the usual transpose and inverse matrix operations, respectively.
Reasonably, as the method was used by Andrews in 1986, a constrained condition is required in the inversion due to a trade-off between independent parameters. For example, Iwata and Irikura (1988) used the constrained condition that the site amplification factors were greater than two at all sites, while Yamanaka et al. (1998) and Nakano et al. (2015) used the theoretical amplification factors evaluated at a reference site. In the present study, we used the approach employed by Yamanaka et al. (1998) and Nakano et al. (2015). The theoretical amplification factors at two reference sites on land were used as constrained conditions. The use of one reference site could be sufficient for the constrained condition. However, using two or more reference sites provides more robust results if the reference sites are the rock sites. The theoretical amplification factors at the two reference sites, namely, IBRH14 and MYGH11, were calculated based on the structure estimated by the inversion of the surface-to-borehole spectral ratios (SBSRs) (see Fig. 1 for the location of the sites). Simulated annealing technique (e.g., Yamanaka 2005) was used to invert the spectral ratios for velocity models. The borehole sensors at the two sites were set up at the depths of 100 and 207 m in the layers having S-wave velocities of 3.2 and 2.78 km/s, respectively, according to the PS-logging data (see Additional file 2). The observed SBSR at the two sites (IBRH14 and MYGH11) computed from several earthquake records and theoretical SBSR using the velocity models obtained from the inversion of the observed SBSR are shown in Fig. 3a, d, respectively. The velocity models are shown in Fig. 3b, e, respectively. The theoretical amplification factors, used as constraints in the spectral inversion, are plotted in Fig. 3c, f, respectively. The layer parameters used to compute the theoretical values of the SBSR and site amplification factors are presented in the Additional file 2. The plots show that the observed spectral ratios are mostly one at frequencies lower than 2 Hz, suggesting that the shallow sediments produce small or no amplification at those frequencies. The theoretical amplification factors of about two at those frequencies show the free surface effects. In contrast to the small amplification at the lower frequencies at the reference sites, amplification factors of about 10 was seen at frequencies around 10 Hz at both sites. The estimated V s values at the depths of borehole sensors were about 2.8 km/s, close to the PS-logging values mentioned previously. Hence, the amplification factors obtained at the other sites from the solution of Eq. (1) might be considered close to the absolute site amplification factors above the seismic basement of ~ 3 km/s (e.g., Nakano et al. 2015). The V s value was 3.45 km/s in Nakano et al. (2015) referred to as the seismic basement. We used V s equal to 3.5 km/s in Eq. (1), and the source, path, and site factors obtained in the present study are described in the next sections.

Source spectra
The obtained acceleration source spectra, S i f , from the inversion were multiplied by ω −2 , where ω is the circular frequency equal to 2π f , to obtain the displacement source spectra. The plots of displacement source spectra are shown in Fig. 4 for nine example events (see Table 1 for details about the events). The ω −2 source spectral model shown in Eq. (3) for far-field S-wave spectra (e.g., Aki 1967;Brune 1970Brune , 1971 were fitted with the displacement spectra, searching for the flat spectral level, , and the corner frequency, f c , which minimized the misfit function defined by Eq. (4) (e.g., Satoh and Tatsumi 2002) for each earthquake separately: where D(f ) is the displacement spectral amplitude at frequency, f : where S f is the source spectral amplitude obtained from inversion at frequency f , D f is the model spectral amplitude at frequency f , and df is the frequency width, f i+1 − f i . Considering the smaller SNRs at lower frequencies and simple forms for the source spectral model, the summation in Eq. (4) was taken over three frequency ranges based on the magnitudes of the earthquakes: 0.2-10 Hz for M w ≤ 5, 0.1-10 Hz for 5 < M w ≤ 6.0, and 0.07-10 Hz for M w > 6. (3) The fitted source spectral models for the nine example events mentioned above are also shown in Fig. 4. It was found that the ω −2 model fits the source spectra generally well within the range of frequencies under consideration. Even though the spectra were fitted for frequencies up to 10 Hz, the fits were reasonably good up to 20 Hz for the crustal and interplate earthquakes (Fig. 4b, c), while the models underestimated the spectra to some extent at frequencies over 10 Hz for the intraslab earthquakes (Fig. 4d). The earthquake types are discussed later.
We estimated magnitudes of the earthquakes from the source spectra obtained in this study and compared them with the catalog magnitudes. First, we calculated the seismic moment, M o , using Eq. (5) (e.g., Brune 1970;Iwata and Irikura 1988;Nakano et al. 2015), and then obtained the magnitude (M w ) using Eq. (6) (Hanks and Kanamori 1979): where ρ is density, V s is S-wave velocity in the source region, R is distance (unit reference distance of 1 km in this study), is flat level of displacement source spectrum, R θφ is the average radiation coefficient, and P is the energy partition ratio. Following Nakano et al. (2015), we used the V s values of 3.6 km/s for the crustal earthquakes and 4.0 km/s for the other earthquakes. Similarly, the densities of 2700 and 3000 kg/m 3 were used for Location map of example events and their source spectra obtained from the spectral inversion. a Focal-mechanism plots coincide with the locations of corresponding epicenters, and the colors of the compression parts indicate the corresponding depths to the foci. b Black lines denote the displacement source spectra obtained from the inversion and red lines show the fitted lines using the ω −2 source model. The labels, C1, C2, C3 means three crustal events, whose catalog magnitudes are indicated with letter M, respectively. The estimated magnitudes based on the flat level of the ω −2 source models are given inside the parentheses. c, d Similar to b, but for the interplate and intraslab earthquakes, respectively. See Table 1 for details about the events Table 1 Source parameters of the earthquakes, whose epicentral locations and source spectra are shown in Fig. 4 Earthquake type (Event code as used in Fig. 4 the source regions of the crustal and other earthquakes, respectively. The values of 0.63 and 1/ √ 2 were used for R θφ and P, respectively, as we used the mean value of the two horizontal components. These values were identical to those used by Nakano et al. (2015). A comparison of the F-net catalog and estimated magnitudes is shown in Fig. 5. We found that the two groups of magnitudes generally agree well, despite the fact that the catalog magnitudes were estimated using much lower frequency ground motions recorded by broadband seismic stations. The differences were mostly smaller than 0.3 magnitude units except for a few events.
Stress drop is an important source parameter used in the simulation of high-frequency ground motions and bears special interest in seismology (e.g., Boore and Atkinson 1987). Here, we discuss the stress drop values assuming simple circular fault models following the previous studies (e.g., Oth et al. 2010;Nakano et al. 2015). First, the source radius for each event was calculated using Eq. (7), and then the stress drop was calculated using Eq. (8) (Brune 1970(Brune , 1971): where r is radius of the source, V s is S-wave velocity, and f c is corner frequency: where �σ is stress drop, and M 0 and r are defined in Eqs.
(5) and (7), respectively. Expressing the seismic moment and radius in Nm and m, respectively, Eq. (8)  The Unclassified (inverted triangles, green) means the earthquakes, whose tectonic types were not ascertained. b Plots of stress-drop residuals against depth. The different symbols correspond to the different earthquake types as indicated in a above. The filled symbols denote the mean residuals for the corresponding earthquake types computed in the interval of 10 km and the vertical bars denote the range of one standard deviation. The mean residuals are shifted slightly in horizontal directions from the centers of the corresponding intervals for clarity stress drop values in units of bar equal to 10 5 Nm −2 . The stress-drop values for the earthquakes used in this study are plotted as a function of depths in Fig. 6. It was found that the stress drops generally increased with the depths (Fig. 6a). The logarithms of stress drops were fitted with depths assuming a linear relationship, which is given in the following equation: where �σ is stress drop in unit of bar as abovementioned, and D is depth in unit of kilometer. The stressdrop residuals obtained as the logarithmic ratios between the results obtained from Eqs. (8) and (9) are shown as a function of depths in Fig. 6b.
To get an idea about the variation of the stress drops between the earthquakes, the earthquakes were divided into four types based on their tectonic origins: interplate (IP) earthquakes, intraslab (S) earthquakes, crustal earthquakes in the Pacific Plate (PAC-C), and crustal earthquakes in the North American Plate (NAM-C). The earthquake types were decided based on their hypocentral locations (JMA), focal mechanisms obtained from the F-net NIED moment tensor solutions, and plateboundary data compiled by Hirose (2022). Some earthquakes could not be classified easily and were grouped as unclassified. It was found that the mean stress-drop residuals for the shallow earthquakes in the Pacific Plate (depths < 30 km) were larger than the mean residuals for the interplate earthquakes within the corresponding depth ranges. The difference was quite large (almost a factor of four, ~ 0.6 on the base 10 log scale) for depths shallower than 10 km. Allmann and Shearer (2009) reported that the stress drops for the intraplate earthquakes in the oceanic lithosphere were twice as high as those for the interplate earthquakes based on a global data set. For depths deeper than 30 km, the differences between the interplate and intraslab earthquakes were minor. The mean residuals for the crustal earthquakes (NAM-C) were small or near zero except for the depth range between 30 and 40 km, where the mean residual was noticeably larger than that of other groups of earthquakes. The mean residuals for the unclassified earthquakes were close to zero except for the depths between 10 and 20 km, in which the mean residual was larger than those from the other groups of earthquakes. However, the number of earthquakes was just two in the depth range, and it is likely that the events belonged to the intraslab type.
The stress drops obtained in Nakano et al. (2015) were slightly larger for the intraslab earthquakes than the interplate ones for the deeper events. The stress-drop values in Nakano et al. (2015) and current study were (9) log 10 (�σ ) = 1.2768 + 0.0124D quite comparable. For example, the average values for the interplate and intraslab earthquakes at a depth of about 40 km were about ~ 50-60 bars in Nakano et al. (2015, Fig . 17), which were close to the value of about 60 bars obtained in the present study. However, we found that the average values of the stress drops for the crustal earthquakes were slightly larger in the present study than those in Nakano et al. (2015). The spatial distributions of the stress drop for the interplate, intraslab, and crustal earthquakes are depicted in Additional file 3. The plots do not show a noticeable lateral variation of the stress drops but depict that the stress drops, on average, increase with the depth of the events, as discussed above.

Path factors
In the present study, we simply assumed the R −1 factor for geometrical attenuation. In this section, we discuss the path attenuation in terms of the quality factor, Q s , for The values obtained in the present study are denoted by circles. The red dashed and solid lines denote the relationships for the shallower (0-30 km) and deeper (30-60 km) parts of the fore-arc region in the north-east Japan (Nakamura 2009). The green dashed line, green solid line, and cyan-colored line denote the relationships for fore-arc region corresponding to the Japan Trench subduction zone for crustal (C), plate boundary (B), and intraslab (I) earthquakes, Region 2, in Nakano et al. (2015). The black dashed line denotes the relationship reported in Oth et al. (2011) for subcrustal earthquakes, and the blue line denotes the relationship reported in Satoh and Tatsumi (2002) for north-east Japan S wave. The Q s values obtained in the present study are plotted as a function of frequency in Fig. 7. The Q s values estimated in several previous studies using similar techniques and earthquake events from similar regions (e.g., Satoh and Tatsumi 2002;Nakamura 2009;Oth et al. 2011;Nakano et al. 2015) are also depicted in the plots for comparison. One major difference between the previous and present studies is that the results in the previous studies were derived using the records primarily from the stations on land as the seafloor observatories were limited or unavailable. As a result, the average paths covered by the seismic waves between the previous and present studies may have varied significantly. The Q s values obtained in this study showed the following features. The values generally increased linearly with frequencies up to 3 Hz, became almost constant between 3 and 10 Hz, and rose gently above 10 Hz (see Fig. 7). The relation between the frequencies and Q s values obtained in the present study is given in Eq. (10) for frequencies between 0.3 and 3.3 Hz. The values between 3.3 and 10 Hz may be considered constant, equal to the value at 3.3 Hz. The values at other frequencies may be regarded as similar to those to be explained shortly.

Nakamura (2009) obtained 3-D attenuation structures
beneath the Japan Islands. The region numbers 4 and 13 in Nakamura (2009) correspond to the 0-30 km and 30-60 km depth ranges in the fore-arc region of northeast Japan. The 30-60 km region overlaps with the present study area to some extent. The Q s values for the two regions in Nakamura (2009) are shown by the dashed and solid red lines in Fig. 7. The applicable frequency ranges in Nakamura (2009) were between 1 and 10 Hz. The plots show that Q s values obtained in this study were generally larger than those for the 30-60 km depth ranges in Nakamura (2009) in the applicable frequency ranges, though the differences were smaller near 1 and 10 Hz. Nakano et al. (2015) obtained three sets of Q s values corresponding to the three types of earthquakes: crustal, plate boundary (or interplate), and intraplate (or intraslab), denoted by the letters C, B, and I, respectively, following the notations in their paper. The Q s values for the region number two in their study, close to the region in the present study, are shown for the different earthquake groups in Fig. 7. The applicable frequency ranges in Nakano et al. (2015) were between 0.5 and 20 Hz for the interplate and intraslab earthquakes, while they were between 0.5 and 5 Hz for the crustal earthquakes. The Q s values at lower and higher frequencies in the present study were close to the Q s values for the intraslab earthquakes in Nakano et al. (2015), while the Q s values (10) Q s = 310f 1.12 at frequencies around 4-5 Hz were close to the crustal earthquakes in Nakano et al. (2015).
The Q s values plotted in Fig. 7, marked with Oth et al. (2011), were for subcrustal earthquakes from wide onshore areas, the larger part of which faces the Pacific Ocean and Philippine Seas. The applicable frequency ranges were between about 0.5 and 20 Hz in Oth et al. (2011). The Q s values plotted in Fig. 7 for Satoh and Tatsumi (2002) were for subduction earthquakes in northeast Japan and were valid between 0.2 and 20 Hz. The Q s values in Oth et al. (2011), Nakanao et al. (2015 for interplate earthquakes (Nakano et al. (2015) B in Fig. 7), and Satoh and Tatsumi (2002) were noticeably smaller than the values in the present study, except for Satoh and Tatsumi (2002) at frequencies higher than about 10 Hz. In summary, the average Q s values obtained in this study are either similar to or larger than those in the previous studies mentioned above. The larger values in the present study are generally as expected, because the seismic rays travel a considerable portion of the high Q Pacific Plate (e.g., Umino and Hasegawa 1984) as the S-net sites are closer to the Pacific Plate compared to the networks in the other studies. Even though the path-averaged Q s values differed between the present and other studies to some extent, the slopes of the fitted lines between the frequencies and Q s values were around unity at most frequencies. Future studies are required to image the detailed Q s structure in the seafloor regions.

Site spectra
We discuss the site amplification factors obtained from the spectral inversion in this section. First, we present the site amplification factors at three KiK-net sites (stations on land), namely, FKSH19, IWTH23, and IWTH27 (see Fig. 1 for the location of the sites). The site amplification factors at the three sites are shown in Fig. 8. The mean SBSRs computed from several earthquake records at the sites are also shown in the figure. The borehole sensors were set up in the layers having V s values of 3.06, 2.2, and 2.79 km/s at the FKSH19, IWTH23, and IWTH27, respectively, according to the PS-logging data. Therefore, the SBSRs may be considered equivalent to the site amplification factors at the sites with respect to a reference rock site. The peak frequencies of the spectral ratios were approximately 3.27, 16.06, and 7.34 Hz at the three sites, FKSH19, IWTH23, and IWTH27, respectively. The peak frequencies and their amplitudes matched well between the spectral ratios and the amplification factors at all the sites. The overall shapes of the curves for the spectral ratios and amplification factors were similar. The values of amplification factors at frequencies lower than about 2 Hz were smaller than two, suggesting no amplification effect by the sediments at the frequencies. The above comparison suggests that the obtained site amplification factors at the three sites were reasonable.
Hereafter, we describe the site spectra at the S-net sites. The site spectra at nine example sites are plotted in Fig. 9 (red lines) along with the individual-earthquake site spectra (grey lines) and theoretical amplification factors (blue lines) at each site computed from the J-SHIS (Japan Seismic Hazard Information Station) velocity model (Fujiwara et al. 2009(Fujiwara et al. , 2012. The sites were from the S3 segment except for one site and were located almost in the east-west direction (see Fig. 1 for the location of the sites). The sites were S3N26 (water depth 128 m), S3N25 (230 m), S3N24 (849 m), S3N23 (1220 m), S3N22 (1645 m), S3N21 (2779 m), S3N20 (5225 m), S3N19 (5591 m), and S6N12 (6111 m), from west to east directions, respectively. The panels a to i in Fig. 9 are arranged in the order of increasing water depths. The curves for the amplification factors from the inversion (red lines) showed multiple peaks and were generally different from site to site. The shallow water sites (a, b, c) showed relatively smaller amplifications (about ten or less) at lower frequencies compared with the values at the deeper water sites (d, e, f, g, h) (about ten or over) except for the deepest site, S6N12, panel i. The S6N12 site, among the nine sites, showed a peak amplification value of about 20 at approximately 6.68 Hz. The spectra at most sites showed a decreasing trend with frequencies higher than about 10 Hz. We explain later that the peak amplification factors at some sites were influenced to some extent by factors irrelevant to the site geology.
The theoretical amplification factors computed for the vertical incidence SH wave using the J-SHIS velocity model above the seismic basement (V s 3.2 km/s) are comparable with the results from the spectral inversion over wider frequencies at the two shallow-water sites, S3N26 (water depth 128 m) and S3N25 (230 m), as shown in Fig. 9a, b. Meanwhile, the theoretical values underestimated the site amplification factors from the inversion at the other sites. The differences generally increased with the increase of water depths. A similar tendency was seen at a few other shallow and deep-water sites (see the Additional file 4). The velocity profiles used to compute the amplification factors at the example sites are shown in the Additional file 4. While the total thicknesses of sedimentary layers above the V s 3.2 km/s layer were between about 1.5 (S6N12, S3N24) and 4 km (S3N19), the thicknesses of layers with V s < 1 km/s were between about 0.1 (S6N12) and 0.9 km (S3N20). The sites located toward the land from the Japan Trench axis showed greater thickness of the sediments than the sites located farther offshore from the trench axis. The obtained site amplification factors at the sites were generally consistent with the distribution of thickness of sedimentary layers in the J-SHIS model. Large amplification factors were obtained at the thicker sediment sites, such as the S3N19, S3N20, and S3N21, where the thicknesses of the sedimentary layers were about 3-4 km (see the Additional file 4). The J-SHIS model in the oceanic region was constructed using the limited geophysical survey data (Fujiwara et al. 2012) and lacked information for shallow sediments with velocities smaller than 600 m/s. Thus, the differences in site amplification factors between the inversion and J-SHIS model may be attributed to the preliminary nature of the velocity model. Moreover, the S-wave time windows selected in this study may have included the 3-D effects of the sedimentary layers around the sites (e.g., Dhakal and Yamanaka 2012). Despite such limitations, the similar amplification factors between the spectral inversion results and the theoretical values at some of the shallow water sites suggest that the site amplification Fig. 9 Comparison of the site amplification factors (red lines) obtained from the spectral inversion with the theoretical amplification factors (blue lines) using the J-SHIS subsurface velocity model. The grey lines denote the site factors corresponding to each earthquake deduced by dividing the original spectra by the source spectrum of corresponding earthquake and path factors. The plots a, b, c, d, e, f, g, h, and i correspond to sites, namely, S3N26, S3N25, S3N24, S3N23, S3N22, S3N21, S3N20, S3N19, and S6N12, indicated in each panel (see Fig. 1 for the location of the sites). The number of earthquake records (N) used in the spectral inversion and the depth to the seafloor (D) are indicated for each site. See Additional file 4 for the J-SHIS model at the sites mentioned above factors obtained in this study may serve as a basis for the validation and improvement of the velocity models in the offshore region. Sawazaki and Nakamura (2020) analyzed the spectral ratios of coda portions of the ground motion records between the horizontal Y and X components and reported that the spectral ratios show a characteristic shape 'N' with peaks and troughs at about 7 and 13 Hz, respectively. They showed that the two frequencies coincided roughly with the natural vibration periods of the cylindrical pressure vessels in the Y and X directions, respectively. They suggested that the natural vibrations were induced due to the poor coupling between the seabed and cylindrical pressure vessels that housed the S-net sensors. The 'N-shaped' spectral ratios were conspicuous at several unburied seismometer sites. To avoid or minimize the effect of the natural vibrations of the pressure vessels on the estimated site spectra, we performed the spectral inversion of only the horizontal X-component spectra, and comparisons of the source, path, and site spectra were made with the values described previously (obtained using the mean spectra of horizontal X and Y components). The moment magnitudes, corner frequencies, stress drops, and Q s values derived from the two inversions are plotted in Fig. 10. We found that the source spectra and path factors from the Fig. 10 a Comparison of the estimated magnitudes (Mw) between the inversions of the X-component spectra and the mean spectra of X and Y components. b, c Similar to the plot a, but for the corner frequencies (Hz) and stress drops, respectively. d Comparison of the quality factors between the two inversions plotted against the frequencies two inversions were essentially the same in a statistical sense.
Example comparisons of the site spectra from the two inversions at nine selected sites, locations being depicted in Fig. 11, are shown in Fig. 12. In addition, shown in Fig. 12 are the mean spectral ratios between the horizontal Y and X component records (Y/X ratios). The site spectra and the Y/X ratios at three sites, namely, S2N14, S2N15, and S2N16, are shown in the panels a, b, and c, respectively. The water depths at the sites were 162, 264, and 874 m, respectively. The three sites were buried sites. It can be seen that the Y/X spectral ratios are near one, and the site spectra from the two inversions (the red and blue curves denoting the results for the mean spectra of two components and X-component spectra only) are very similar. Similarly, the panels d, e, f, g, h, and i in Fig. 12 show the site spectra and Y/X spectral ratios for the other six sites, where the sensor houses were not buried. The water depths ranged between 2247 and 6230 m. The site spectra from the two inversions showed similitude in shape and peak values. The Y/X spectral ratios had obvious troughs at frequencies near 10 Hz or higher at all the sites, but the peaks were not so prominent except for the two sites, namely, S2N21 and S2N22 (see the panels h and i in Fig. 12). The peak frequencies of the site spectra and the Y/X spectral ratios were generally different. However, we found that the site spectra at 30 sites (20% of the total stations) obtained from the inversion of the mean spectra of X and Y components were noticeably contaminated at frequencies that corresponded with the peaks and troughs of the Y/X spectral ratios. The list of 30 sites, along with the peak frequencies and other properties, is provided in the Additional file 5 (see Fig. 11 for the location of the sites). Out of the 30 sites, the example plots of the site spectra from the two inversions and the Y/X spectral ratios are shown in Fig. 13 for the selected nine sites. The peak amplitudes of the mean Y/X spectral ratios ranged between about 3 and 9, and the peak frequencies were between about 4 and 9 Hz at the 30 sites. The upper three panels (a, b, c) in Fig. 13 depict the spectra for the three sites, for which the maximum values of the Y/X spectral ratios were about 3; the sites in the middle three panels (d, e, and f ) had the maximum Y/X spectral ratios of about 4. The last three panels (g, h, and i) had the maximum Y/X spectral ratios of about, 7, 7.5, and 9, respectively. The plots depict that the site spectra from the two inversions were generally similar except for the frequencies at and around the peaks and troughs of the Y/X spectral ratios. We found that the difference between the two inversion results becomes conspicuous when the peak value of the Y/X spectral ratios is about 3 or larger.
It is interesting to see how the amplification factors from the spectral inversion were spatially distributed over different sites for a given frequency to get a general picture of the variation in the site spectra. The spatial distributions of the site amplification factors for three frequencies, approximately 0.33 Hz, 1 Hz, and 5 Hz, are shown in Fig. 14b-d, respectively. The amplification factors at 5 Hz are plotted from the analysis of X-component records only to avoid the spurious spectra discussed above. Figure 14a shows the two-way travel time in the sedimentary layers for the P wave obtained from the multi-channel seismic (MCS) reflection survey (Nishizawa et al. 2022). The outer trench sites, part of the S6 segment, mostly showed smaller values (~ 0.5 s) for the travel times. Similarly, the twoway travel times were smaller than 1 s at several sites close to the coast and in the S1 segment. At the other sites, the two-way travel times were mostly larger than 1 s. The longest two-way travel time was about 2.5 s beneath the S2N21 site (see Fig. 14a). The longer travel times indicate thicker low velocity sediments compared with the smaller ones. The site amplification Fig. 11 Red circles denote the locations of the sites for which the site spectra are plotted in Fig. 12. The site codes are written near the site symbols. Red triangles denote the locations of the sites, where the site spectra at higher frequencies were spurious (see Fig. 13) factors for 0.33 and 1 Hz, depicted in Fig. 14b, c show the patterns roughly similar to the distribution of the two-way travel times shown in Fig. 14a. In contrast, the amplification factors for 5 Hz were larger than the ones for the lower frequencies at many sites in the S6 segment, where the two-way travel times were smaller. The median values of the amplification factors over the S-net sites were about 10 for the three frequencies, Fig. 12 Example comparison of the site spectra obtained from the spectral inversions of mean spectra of X-and Y-component records (red lines) and spectra of X-component records (blue lines). The plots a, b, c, d, e, f, g, h, and i correspond to sites, namely, S2N14, S2N15, S2N16, S2N17, S2N18, S2N19, S2N20, S2N21, and S2N22, as indicated in each panel. See Fig. 11 for the location of the sites. The black lines denote the mean spectral ratios between the horizontal Y-and X-component records. The number of earthquake records (N) used for each site and the depth to the seafloor (D) are indicated in each plot Dhakal et al. Earth, Planets and Space (2023) 75:1 Fig. 13 Example comparison of the site spectra obtained from the spectral inversions of mean spectra of X-and Y-component records (red lines) and spectra of X-component records (blue lines) at sites, where the spurious site spectra were recognized. The plots a, b, c, d, e, f, g, h, and i correspond to sites, namely, S4N26, S6N12, S4N10, S5N17, S5N15, S4N09, S6N16, S5N14, and S1N06, as indicated in each panel. See Fig. 11 for the location of the sites. The black lines denote the mean spectral ratios between the horizontal Y-and X-component records. The number of earthquake records (N) used for each site and the depth to the seafloor (D) are indicated in each plot Fig. 14 a Two-way travel time for P wave from the base of sediment to the seafloor estimated by multi-channel seismic (MCS) reflection survey (Nishizawa et al. 2022). The color filled squares and circles denote the sites located, within and beyond 5 km from the survey line, respectively. The open squares denote the sites, where the travel times were not estimated. b, c, and d Amplification factors at frequencies of approximately 0.33, 1, and 5 Hz, respectively discussed above. The smaller travel times and larger amplification for the higher frequencies might suggest the thickness of the sediment to be relatively small in the outer trench sites. Azuma et al. (2019) obtained sediment thicknesses based on travel-time differences between P and PS converted waves at the sedimentary basement below the S-net sites and suggested that the sediment thicknesses in the landward of the trench tend to be thicker than those estimated by two-way travel times of the deepest reflectors on the MCS profiles (Nishizawa et al. 2022). We plan to integrate the travel-time information, site spectra, and the available subsea models to reconstruct the improved velocity models in our future studies. Finally, the relationships between the two-way travel times and amplification factors at the frequencies of 0.33, 0.5, 1, 2, 5, and 10 Hz are presented in Fig. 15a-f, respectively. The amplification factors from the two inversions, using the mean spectra of X and Y components and X-component spectra only, showed similar trends with the two-way travel times as expected. The plots for the 0.33, 0.5, 1, and 2 Hz indicated that the amplification factors generally increased with the two-way travel times, albeit the small values of correlation coefficients (~ 0.25). We found that the correlation coefficients decreased with the increase in frequency over 2 Hz and was approximately zero at 10 Hz. The mean amplification factors were approximately 10 and 2 at frequencies of 5 and 10 Hz, respectively. These latter results may suggest that the thicker sediment layers dampen the higher frequency components more than the lower ones due to the lower Qs values in the sediments (e.g., Boore and Smith 1999).

Fig. 15
Relationships between the two-way travel times for P wave and site amplification factors at frequencies of approximately 0.33, 0.5, 1, 2, 5, and 10 Hz (panels a-f), respectively. The red and black lines indicate the linear fits between the travel times and logarithms of amplification factors. The numerals in parentheses indicate the values of correlation coefficients between the corresponding data sets. Note the difference in the scale of the vertical axis for 10 Hz f compared to the scales for other frequencies. The amplification factors are lower than five for 10 Hz at most sites

Conclusions
We analyzed the strong-motion records at the S-net sites following the earthquakes with magnitudes 3.5 < M w ≤ 7 and focal depths < 70 km and evaluated the source, path, and site spectra of S waves for understanding the S-wave site amplification factors at the sites. We used the spectral inversion technique to obtain the source, path, and site spectra simultaneously from the horizontal-component S-wave portions of the recordings.
The source spectra were fitted with the conventional ω −2 source model searching for the flat levels of the spectra and corner frequencies, and it was found that the source spectra followed the ω −2 source model generally well for most earthquakes up to 10 Hz. The estimated magnitudes from the flat levels of the source spectra were mostly within ± 0.3 magnitude units of the F-net NIED catalog magnitudes. The stress drops increased with focal depths, which is consistent with the previous studies in the region. The stress drops for the shallow earthquakes in the Pacific Plate were higher than those for the interplate earthquakes with comparable focal depths, while the values were similar for events with depths greater than about 30 km. The pathaveraged quality factors increased with frequency up to about 3 Hz, remained almost constant between about 3 and 10 Hz, and then increased gently up to 20 Hz. The values were either similar to or larger than those reported in the previous studies.
The peak site frequencies were between about 0.2 and 10 Hz, and the peak amplification factors ranged between about 10 and 50. The peak frequencies and amplification factors generally differed from site to site. However, a moderate regional pattern was observed. The peak frequencies were mostly higher than about 2 Hz at the outer trench sites, while they were lower than 2 Hz at many inner trench sites. The amplification factors at a few shallow-water sites close to the coast were comparable with the theoretical ones computed from the J-SHIS subsea model. The amplification factors from the inversion were larger in the regions, where the sedimentary layers were thicker in the J-SHIS model. However, the amplification factors were considerably underestimated by the J-SHIS model at many sites over wide frequencies. The amplification factors at intermediate frequencies (~ 0.3 to 2 Hz) generally increased with the P-wave travel time in the sediments estimated from the multi-channel seismic reflection survey. All the above results, by and large, suggested that the site spectra were reasonable. However, it was found that the site spectra included dominant peaks at frequencies between about 4 and 10 Hz at about thirty unburied stations, resulting from the coupling problem between the instrument vessels and seabed sediments.
We performed the spectral inversion using only the spectra from the horizontal X-component records and found that the spurious peaks at those sites were eliminated, while the source spectra and path factors were very similar to those from the joint use of both X-and Y-component records. We recommend avoiding the Y-component records if the site spectra between about 4 and 10 Hz are required at the unburied stations for any application. Moreover, the site spectra over about 10 Hz may not be reliable due to the higher modes in the Y directions and vibrations of the sensor vessels in the X directions.
The results discussed in this paper may be used in the prediction of ground motions for EEW, engineering design of offshore structures, and so on. It is expected that the abovementioned results may also serve as a basis for more detailed future studies regarding the source properties of the subduction zone earthquakes, quality factor of the crust and mantle, and improvement of the subsurface velocity models in the region.