Crustal attenuation characteristics of S-waves beneath the Eastern Tohoku region, Japan

An inversion method was applied to crustal earthquakes dataset to find S-wave attenuation characteristics beneath the Eastern Tohoku region of Japan. Accelerograms from 85 shallow crustal earthquakes up to 25 km depth and magnitude range between 3.5 and 5.5 were analyzed to estimate the seismic quality factor Qs. A homogeneous attenuation model Qs for the wave propagation path was evaluated from spectral amplitudes, at 24 different frequencies between 0.5 and 20 Hz by using generalized inversion technique. To do this, non-parametric attenuation functions were calculated to observe spectral amplitude decay with hypocentral distance. Then, these functions were parameterized to estimate Qs. It was found that in Eastern Tohoku region, the Qs frequency dependence can be approximated with the function 33 f1.22 within a frequency range between 0.5 and 20 Hz. However, the frequency dependence of Qs in the frequency range between 0.5 and 6 Hz is best approximated by Qs (f) = 36 f0.94 showing relatively weaker frequency dependence as compared to the relation Qs (f) = 6 f2.09 for the frequency range between 6 and 15 Hz. These results could be used to estimate source and site parameters for seismic hazard assessment in the region.


Introduction
The study region is confined within 38°N to 39.5°N and 140.5°E to 142.5°E, and contains eastern part of Tohoku region in Japan's Honshu Island as shown in Fig. 1. Eastern Tohoku mainly comprises Iwate and Miyagi prefectures. The region is situated over one of the most active tectonic regimes of the world where Pacific plate is subducting under the Eurasian plate, generating numerous earthquakes in this area.
There is a close connection between the topography of the Tohoku region and the distribution of active faults. The topography of the Tohoku region is characterized by alternate mountain ranges and lowlands running in a northsouth direction. The primary active faults of the Tohoku region lie at the boundaries of the mountains and the lowlands.
The crustal depth beneath Tohoku is estimated to be 25-30 km (Kaminuma and Aki 1963;Zhao et al. 1992) and many of the large destructive inland earthquakes are known to have occurred at these shallower depths. The shallower crustal earthquakes tend to result in a greater intensity of surface shaking and often cause the greatest loss of life and damage to property. In engineering applications like seismic hazard analysis, S-wave phases of ground motions are of most interest, as most of the seismic energy from the earthquake source is radiated in the form of S-waves. Therefore, quantification of the S-wave attenuation characteristics is necessary.
The records of strong earthquakes can be inverted to obtain the characteristics of attenuation, radiation and propagation pattern of seismic waves and to construct models of ground behaviour at various stations (Pavlenko and Wen 2008). In this study, linear inversion analysis following a non-parametric approach has been performed to investigate the frequency dependence of attenuation of S-waves (Q s ).

Dataset
The dataset comprises 315 three-component earthquake records, from 85 earthquakes recorded at 6 K-NET (Kyoshin network) stations installed in the Eastern Tohoku region of Japan in the framework of National Research Institute for Earth Science and Disaster Prevention (NIED). The earthquake records were obtained over a period of nine years from Dec-2003 to May-2012 with magnitudes M w ranging between 3.5 and 5.5 and hypocentral distances between 35 and 125 km. The data have been sampled at 100 Hz. The data represent no distance dependence in magnitude as shown in Fig. 2.
The site details of the K-NET stations used in this study are listed in Table 1. These stations were initially installed with K-NET95-type accelerometers which were later replaced with K-NET02-type accelerometers in 2004 (Okada et al. 2004). Since, the crustal thickness of Eastern Tohoku region ranges in between 25 to 30 km, the hypocentral depth of the selected events was restricted not to exceed 25 km. Moreover, only those earthquakes were selected which were recorded by at least three stations. The details of epicentral coordinates of the events used and the hypocentral distances to the recording stations are provided in Appendix Table 2. The ray paths between different source-station pairs in the region are shown in Fig. 3. There are many crossing ray paths between different source station pairs ensuring reliable estimates of homogeneous attenuation model.

Data processing and analysis
The records were corrected for instrumental response and baseline. Instrument correction was only applied on the records from K-NET95-type accelerometers for frequencies larger than *15 Hz, as the response of these accelerometers slightly decrease for frequencies larger than this limit. The baseline was corrected by subtracting the average of all the points of the record. S-wave portion of both the horizontal components of the recorded accelerograms were analyzed. For each horizontal component of the records, the S-wave window was selected such that it starts 1 s before the S-wave onset and ends when 80 % of the total energy of the record is reached as shown in Fig. 4. Typical window lengths range between 5 and 12 s. If in any case the window's length increases over 20 s, it was fixed to a maximum duration of 20 s to avoid having too much coda energy in the analyzed time window (Oth et al. 2011). The beginning and the end of S-wave window were tapered with a 5 % cosine taper.
A plot showing the build-up of Arias intensity with time is known as a Husid plot (Fig. 4) and it serves to identify the interval with the arrival of majority of the energy.
The Fourier amplitude spectra (FAS) were computed for each window and smoothed around 24 frequency points equidistant on logarithmic scale between 0.5 and 20 Hz using the running mean filter. This smoothing technique is optimal for reducing random noise while keeping the sharpest step response (Smith 1999). The smoothing bandwidth is determined by trial-and-error. Having tested various smoothing bandwidths, the following criteria were adopted: for each frequency point, a bandwidth equal to 0.3 of an octave was considered; if this bandwidth was smaller than 0.5 Hz, then 0.5 Hz was used as the smoothing  bandwidth. The signals to be analyzed have to satisfy stringent signal-to-noise ratio (SNR) requirements. Only data points with SNR greater than three were included in the final dataset. Figure 5 shows the raw and smoothed FAS of the chosen S-wave window from an accelerogram recorded at station IWT007. At lower frequencies, the database is slightly sparser due to SNR constraints as shown in Fig. 6.

Method
The generalized inversion technique (GIT) has been widely applied to crustal earthquake datasets (Andrews 1986;Castro et al. 1990). As a result of inversion, the frequency dependence of the attenuation of seismic waves as well as source characteristics and site response can easily be found (Parolai et al. 2000;Bindi et al. 2006).  These dataset from Eastern Tohoku region allows us to make stable and reliable estimates of attenuation properties as it covers a reasonable range of magnitudes, distances and focal depths. Moreover, each event in the dataset was recorded by at least three stations, creating multiple crossing ray paths from sources to stations. Under these conditions, the effect of in-homogeneities in whole-path attenuation is thought to be effectively averaged out and GIT could successfully be applied to our dataset.
In this study, a non-parametric approach was employed (Castro et al. 1990), which considers attenuation phenomenon to be a smooth function, decreasing with distance. In non-parametric approach, the inversion is performed in two steps. Since, the attenuation character-isticsÂ ij ðf ; r ij Þ are determined from the first step, we restrict the procedure to the first step of inversion. However, a subsequent inversion can be performed to separate source and site effects (Andrews 1986;Oth et al. 2009).

Non-parametric approach
The dependence of the spectral amplitudes U (f, r) at frequency f and on distance r may be written as where U ij (f, r ij ) is the spectral amplitude (acceleration) from the ith earthquake at the jth station resulting, r ij is the source-site distance, the non-parametric attenuation function (NAF)Â ij ðf ; r ij Þ describes seismic attenuation with distance, andŜ i ðf Þ is a scaling factor depending on the size of the ith event. With this approach, the attenuation function remains independent of the size of an event.Â ij ðf ; r ij Þ is not supposed to have any specific shape and implicitly Example of distribution of hypocentral distances at the two selected frequencies, 0.5 and 5.54 Hz, respectively. At 0.5 Hz, the database is slightly sparser than at 5.54 Hz due to signal-to-noise ratio constraints contains all effects leading to attenuation along the travel path (geometrical spreading, anelasticity, scattering, etc.). Based on the idea that these properties vary slowly with distance,Â ij ðf ; r ij Þ is constrained to be a smooth function of distance i.e., with a small second derivative (Castro et al. 1990) and to take the valueÂ ij ðf ; r 0 Þ ¼ 1 at some reference distance r 0 . The representation given by Eq.
(1) does not include a factor related to the site, and hence the site effects are necessarily absorbed both inÂ ij ðf ; r ij Þ andŜ i ðf Þ. This is the reason that '^' is used to distinguish them from the 'uncontaminated' path and source termsÂ ij ðf ; r ij Þ and S i (f), respectively.Ŝ i ðf Þ contains, in fact, an average value of the site amplifications of all stations that recorded the ith earthquake. Thus,Ŝ i ðf Þ is not the true source spectrum of the ith event.
Equation (1) can be easily linearized by taking the logarithm: Equation (2) represents an over-determined system of the form Ax = b, where b is the data vector containing the logarithmic spectral amplitudes, x is the vector containing the model parameters, and A is the system matrix relating the two of them. The distance range was subdivided into nine distance bins (N D ), each 10 km wide, and the value of A ij ðf ; r ij Þ was computed in each bin. In matrix formulation, Eq. (2) takes the following form.
The left part of the system matrix in Eq. (3) contains the factors related to the attenuation parameters, whereas the right-hand side reflects those related to the source terms. The system matrix also includes rows relevant to constraints. The weighting factor w 1 is used to imposeÂ ij ðf ; r 0 Þ ¼ 1 at the reference distance r 0 . Since the dataset contains most of the events with hypocentral distance greater than 40 km, reference distance r 0 is set to 40 km. By settingÂ ij ðf ; r 0 Þ ¼ 1 at r 0 = 40 km is, in fact, equal to assuming that there is no attenuation over that distance from the source. Hence, there is a cumulative attenuation effect over these 40 km that is impossible to resolve, and therefore, the Q s -model derived from the slopes of the attenuation functions, only reflects the attenuation characteristic over the remaining part of the travel path (i.e., more than 40 km away from the source). The weighing factor w 2 is implemented to achieve monotonically decaying attenuation curves with reasonable degree of smoothness to suppress the site-related effects and yet preserve variations of the attenuation characteristics with distance (Oth et al. 2009).
At each of the 24 selected frequencies, an inversion was performed in a least-square sense and a solution x = (A T A) -1 A T b for a numerically stable system was computed (Menke 1989). As a result of successful inversion, the modal matrix gives the NAFsÂ ij ðf ; r ij Þ one for each bin and the values ofŜ i ðf Þ, one for each earthquake i.  evaluated in the form of sum of the squared residuals. The residual between the observed spectral amplitudes and those predicted by U ij ðf ; r ij Þ ¼Â ij ðf ; r ij Þ ÁŜ i ðf Þ, are interpreted as site effects. The attenuation curves evaluated at different frequencies are shown in Fig. 7.

Seismic quality factor (Q s )
The attenuation functionÂ ij ðf ; r ij Þ combines the effect of the geometrical spreading and the anelastic attenuation in the same function. By using the non-parametric form of attenuation function, it becomes very easy to test it against any assumed geometrical spreading without repeating earlier computations. The attenuation termÂ ij ðf ; r ij Þ can be parameterized in terms of frequency-dependent quality factor Q s (f) and geometric spreading G (r). Thus, we can express the attenuation function asÂ The average S-wave velocity (b) estimated in the region is 3.2 km/s, measured from the S-P arrival times (Zhao et al. 1992;Kurahashi and Irikura 2011). The Q s estimates are also sensitive to the choice of geometrical spreading function. Considering body wave propagation in an infinite homogeneous medium, the geometric spreading function G (r) was chosen to be 1/r. The amplitudes are normalized to 40 km as most of the observed spectral amplitudes start at 40 km. Thus, for each frequency f, Eq. (4) is linearized correcting the empirical attenuation functions by the effect of geometrical spreading G (r) = 40/r ij and taking the logarithm. Thus, Eq. (4) is written as The functions are normalized to zero (in logarithm) at the reference distance r 0 where aðrÞ ¼ log 10Âij ðf ; r ij Þ À log 10 Gðr ij Þ and m = {[pfr/ Q s (f)b]Álog 10 (e)}. It is important to note that m is the resulting slope of the linear least-squares fit obtained between 40 and 120 km for each frequency analyzed. Therefore, Eq. (5) takes the following form: It becomes evident that by correcting the attenuation functions for geometrical spreading and plotting versus distance, Q s (f) can be evaluated from the slope of a linear least-squares fit. Figure 8 shows the attenuation functions corrected for the geometrical spreading and the regression used to estimate Q s for different frequencies.
A homogeneous attenuation model for the studied region from a linear fit of the determined values of Q s over the selected frequency band of 0.5 to 20 Hz takes the form Q s (f) = 33 f 1.22 as shown in Fig. 9. The regression error of the fit is shown by vertical error bars and can be significantly reduced by choosing multiple Q s models over different frequency bands.

Results and discussion
The obtained NAFs show that in general the attenuation curves are well constrained. At higher frequencies (f [ 12 Hz), initially the change in rate of amplitude decay is fast but gradually slows down beyond 70 km distance. A strong frequency dependence of Q s is found between the frequency range of (0.5 and 20 Hz) and is given by the relation Q s (f) = 33 f 1.22 . Below 6 Hz, Q s can be best approximated by Q s (f) = 36 f 0.94 showing relatively weaker frequency dependence as compared to relation Q s (f) = 6 f 2.09 for the frequency range between 6 and 15 Hz. This change in frequency behaviour of Q s can be attributed to the localized discontinuities within the crust which are present in the form of faults or densely fractured zones in the region (Umino et al. 1990). A strong attenuation in Tohoku region is quite obvious, as evident from several past studies. The Q s estimates found in this study are comparable to those obtained for the Noto region in Japan (Itoi et al. 2008). Both functions are obtained from the crustal earthquake datasets. The slopes of the functions are quite similar, leading to the determination of almost identical Q s (f)-model as shown in Fig. 10. The derived Q s (f)-model is also comparable to those of crustal earthquakes at other regions in Japan.

Conclusions
Frequency-dependent seismic quality factor Q s (f) has been evaluated by applying non-parametric inversion method to crustal earthquakes dataset. A strong attenuation in Eastern Tohoku region has been observed between (0.5 and 20 Hz) and is given by the relation Q s (f) = 33 f 1.22 which is close to the background value in the region. This strong attenuation can be attributed to structural heterogeneities created by large earthquakes in the region. The estimates of Q s could be further refined by adopting a bilinear Q s (f)-model. The frequency dependence of Q s in the frequency range between 0.5 and 6 Hz can be best approximated by Q s (f) = 36 f 0.94 showing relatively weaker frequency dependence as compared to the relation Q s (f) = 6 f 2.09 for the frequency range between 6 and 15 Hz. The observed estimates of Q s could provide basic input for determining source and site parameters (Parolai et al. 2000). Therefore, more realistic estimates of strong ground motion parameters could be obtained for assessment of seismic hazards in the region.

Acknowledgments
The study was carried out as a part of author's M.Sc Research under the project: ''Strengthening of Earthquake Engineering Center'', funded by Higher Education Commission, Government of Pakistan. The data used in this study were downloaded online from K-Net (Kyoshin Network) website: http:// www.kik.bosai.go.jp/, managed under the framework of National Research Institute for Earth Science and Disaster Prevention (NIED). The author wishes to thank Professor Apostolos Papageorgiou for his valuable guidance and support to carry out this research.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix
See Table 2.