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

1.   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.

Figure  1.  Map of Japan showing Iwate and Miyagi prefectures which constitute Eastern Tohoku

DownLoad: Full-Size Img PowerPoint

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 north-south 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).

2.   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.

Figure  2.  Distribution of magnitudes with hypocentral distances

DownLoad: Full-Size Img PowerPoint

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.

Table  1.  Location of the seismic stations used in this study

 | Show Table

DownLoad: CSV

Table  2.  Earthquake coordinates and hypocentral distances

 | Show Table

DownLoad: CSV

Figure  3.  Map of the Eastern Tohoku region that showing the distribution of the epicentres (stars) and the stations (triangles)

DownLoad: Full-Size Img PowerPoint

3.   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.

Figure  4.  Accelerogram of one of the earthquakes recorded at the station IWT007 (lower panel), highlighting the S-wave duration and the corresponding Husid plot (upper panel)

DownLoad: Full-Size Img PowerPoint

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.

Figure  5.  Fourier amplitude spectrum (FAS) of the selected S-wave window from an accelerogram recorded at station IWT007. Raw and smoothed acceleration spectrums are shown by solid grey and solid black lines, respectively. Noise spectrum (dashed black line) is also shown for the same event to illustrate signal to noise levels

DownLoad: Full-Size Img PowerPoint

Figure  6.  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

DownLoad: Full-Size Img PowerPoint

4.   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 characteristics 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).

4.1   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) describes seismic attenuation with distance, and 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. is not supposed to have any specific shape and implicitly 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, 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 at some reference distance r₀.

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 and . This is the reason that '^' is used to distinguish them from the 'uncontaminated' path and source terms and S_i (f), respectively. contains, in fact, an average value of the site amplifications of all stations that recorded the ith earthquake. Thus, 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 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₁ is used to impose at the reference distance r₀. Since the dataset contains most of the events with hypocentral distance greater than 40 km, reference distance r₀ is set to 40 km. By setting at r₀ = 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₂ 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 one for each bin and the values of , one for each earthquake i. In this way, the unknown values which reduce the deviation between observations and model predictions can be evaluated in the form of sum of the squared residuals. The residual between the observed spectral amplitudes and those predicted by , are interpreted as site effects. The attenuation curves evaluated at different frequencies are shown in Fig. 7.

Figure  7.  Attenuation profiles observed at 9 of the 24 frequencies studied (solid lines) by fitting one function to the entire dataset (dashed lines). The functions are normalized to zero (in logarithm) at the reference distance r₀

DownLoad: Full-Size Img PowerPoint

4.2   Seismic quality factor (Q_s)

The attenuation function 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 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 (β) 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

where and m = {[πfr/ Q_s(f)β]·log₁₀(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.

Figure  8.  Attenuation curves corrected for geometrical spreading [log₁₀ A(f, r) - log₁₀ G(r) versus r] and fitted straight line (least squares fit) at 9 of the 24 frequencies studied

DownLoad: Full-Size Img PowerPoint

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.

Figure  9.  Estimated Q_s-model for Eastern Tohoku obtained from the least square fit of Q_s between 0.5 and 20 Hz

DownLoad: Full-Size Img PowerPoint

5.   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.

Figure  10.  Comparison of Q_s for crustal earthquakes in Japan

DownLoad: Full-Size Img PowerPoint

6.   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://creativecommons.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.