S-wave attenuation structure in the Taiwan area and its correlation to seismicity

S-wave amplitude data from the Taiwan Telemetered Seismographic Network (TTSN) and the Central Weather Bureau Seismographic Network (CWBSN) are used in a tomographical inversion for a three-dimensional S­ wave attenuation structure beneath the Taiwan area. The data consist of three-dimensional, digital seismograms of 137 earthquakes with magnitudes between 3.0 and 5.0. The S-wave attenuation structure is obtained through the use of the damped least-squares method from attenuation time value (t*) and travel times. Attenuation time value is calculated by using the spec­ tral decay method under the assumption that t* is independent of frequency. S-wave spectra are corrected for a ro2 source spectrum model. The fre­ quencies of instrument response above 25 Hz or below 4 Hz usually have large uncertainty. Therefore, the least-squares line fitting is performed in the frequency band of 4 to 25 Hz. The results here show that the distribu­ tion pattern of the Q8 values in the Taiwan area correlates well with seis­ micity. In the upper crust (depths of between 0 to 15 km), the region of lower Q s area correlates with higher seismicity. However, in the lower crust and upper mantle (15 to 75 km), the area of higher seismicity correlates with the regions of higher Q s values. (


INTRODUCTION
The anelasticity of the Earth, as measured by the specific quality factor Q, has long been recognized as an important indicator of the composition, state and temperature of the deep interior. This seismic parameter is not well understood and is often difficult to measure. How ever, its measurement is essential for the prediction of ground motion produced by earth quakes or explosions, the determination of earthquake magnitude, the understanding of re gional geodynamics and the evaluation of seismic risk. Attenuation estimates are also very useful in characterizing certain physical properties of rock materials that cannot be inferred from seismic velocities.
Studies on the attenuation of seismic waves have occupied a significant portion of seismic literature over recent years (Der and McElfresh, 1976;Al-Shunkri et al., 1988;Kanamori, 1967a;1967b;Lee and Solomon, 1978;Solomon, 1972;1973;Teng, 1968). Attenuation measurements can carry much information regarding the physical state of the Earth's interior and are important in making corrections in source mechanism studies. Observations of lateral variations in Q may aid in determining regional differences in temperature, phase and compo sition.
The tectonic features of the Taiwan area are extremely complicated because of the nearby convergence of the Eurasian and Philippine Sea plates. The collision of the Philippine Sea and Eurasian plates as well as the spreading of the Okinawa Trough are the main reasons for the complex geological structures and high seismicity in this region (Figure 1). Roecker et al. (1987) applied a three-dimensional block inversion method to obtain the velocity structure of this area and reported the existence of velocity heterogeneity in the uppermost mantle. It is believed that the high-velocity zones are generally associated with the subducted Philippine Sea plate. Moreover, Chen (1993) investigated the corresponding attenuative property of P waves beneath this region. He pointed to the presence of a high-Q P zone associated with the subducted Philippine Sea plate in northeastern Taiwan and the subducted Eurasian plates in the south. Chen et al. (1996) determined the Q P structure in the Taiwan area and its correlation with seismicity. In the present study, seismic data as recorded by the TTSN and CWBSN are used to obtain a three-dimensional model of the Q8 structure beneath the Taiwan area. Figure 2 shows the distribution of the stations for these two networks. The CWBSN has 75 stations including (those of) the TTSN stations. These authors expect that the Q8 structure provide useful information for a better understanding of the material beneath this area. For this pur pose, the spectral decay of the seismic S phase is used as basic observed data. The three dimensional Q8 structure is parameterized by 9 layers bounded by horizontal interfaces, each of which is separated into a mesh of rectangular prisms by two orthogonal sets of vertical interfaces (see Figure 6). Totally, the study region is divided into 940 "blocks". The interfaces (including horizontal and vertical) are similar to those of the velocity model IV3 used by Roecker et al. (1987). The Q8 values in these blocks are estimated as the output of an inverse problem from the observed data.

SPECTRAL DECAY METHOD
Q is a dimensionless quantity that measures the amount of energy dissipated per radian as seismic waves propagate. In other words, a seismic wave loses 1/Q of its energy after traveling each radian owing to the anelasticity of the medium. Laboratory studies have previous shown Q to be independent of frequency (Knopoff, 1964). However, theoretically, Q varies signifi cantly with a function of frequency for a rock sample of uniform composition and grain size . Also, in a composite viscoelastic model with a range of relax ation times, Q varies smoothly with frequency and may be considered (to be) constant over a narrow frequency band (Kanamori and Anderson, 1977). The basic assumption that Q is inde pendent of frequency within the interested frequency band is often made in studies of seismic  Ho, 1986), showing the interaction of the Philippine Sea and Eur asian plates in Taiwan with the Ryukyu arc system developing to the east and northeast, and the Luzon arc system extending to the south.
The effect of Q on a seismic signal along a ray path describes the variation of the displace ment spectrumA(r,.f), namely: A(r,f) = Ao(r,f)e-n fi * , in which V is the wave velocity, Q is the spatial attenuation factor, and the integral is along the ray path s (Kanamori, 1967 a). The spectrum of displacement, observed at a site distance r from the epicenter, A(r,f), may be approximated by: where G(r,f) gives the effect of geometrical spreading that is usually assumed to be indepen dent of frequency, l(f) is the instrumental response, and S(f) is the displacement spectrum of radiation at the source. In a homogeneous medium, G( r,f) can be approximated by 1/r for body waves. At high frequencies, a widely accepted source spectrum model is the co 2 model, which assumes that S(f) is proportional to f-2 (Brune, 1970;Hanks, 1979). This model was used to study the attenuation of high-frequency waves from teleseisrnic sources (Der and McEl fresh, 1976;Lundquist and Cormier, 1980;Der et al., 1982;1985), and at regional distance (local earthquakes) (Modiano and Hatzfeld, 1982;. Since the primary interest here is not the absolute value oft* but rather its spatial variation within the crust, assuming a co 2 model should present no difficulties in drawing . conclusions. Under this assumption, AJ r, f) in Equation (1) becomes independent of frequency.
Since the instrumental response I(f) i s usually known for each component at all stations, the effect can be removed from the above spectrum A(r, f). In other words, Equation (1) can be rewritten, after taking the natural l ogarithm on both sides, as: lnA(r,f) = lnAo(r)-rrjt * .
The termA (r) in Equation (4 ) includes the geometric spreading effect, which is independent 0 of frequency. Since the spreading effect is frequency-independent, the value of A (r) only 0 affects the level of the spectrum but has no effect on the observed t*. If Q is frequency-independent, then t* is also frequency-independent. Thus, Equation ( 4) shows lnA(r,f) as a straight line in/ and can be differentiated with respect to/to obtain the t*: The corrected observations, in A(r,f) , is fitted to a straight line by using the least-squares method over a frequency range extending somewhat beyond t he comer frequencies. However, many studies have pointed out that Q is frequency-dependent in the crust under the continental regions, at least at frequencies near I Hz (for example : Mitchell, 1980), and the degree of the dependence varies regionally (Mitchell, 1981). But un fortunately, i t is almost impossible to investigate the frequency-dependency of Q with the data available in this study. Since the primary goal here is to find the spatial variations of attenuation in the upper crust of the Tai wan area:, the use of frequency-independent attenuation is assumed. For earthquakes used in this study, the comer frequencies range from about 1to4 Hz. A frequency band of 4-25 Hz is chosen in this study. The former is limited by the comer of the source frequency, while and the latter is limited by the uncertainty in instrument response above 40 Hz.

DATA AND COMPUTATION
Digital seismograms from 137 (from 1987 to 1994) local events in the Taiwan region are used to obtain the body-wave spectra. A total of 1523 t* values are obtained . The magnitudes of these earthquakes range from 3.0 to 5.0. Figure 3 shows the location of these earthquakes and the seismic stations used in this study. The first 2.0 second starting from t he onset time of the S-waves of the seismograms are analyzed in term of S phase. The Fast Fourier Transform ( FFT ) method is used to study the S phase spectra. A frequency band from 4 to 25 Hz is selected in this study. The spectrum for frequencies lower than the comer-frequency and eigen frequency are rather unreliable, while that for those higher than 25 Hz has a very small ampli tude due to a high cut-off filter i n the recording system. The slope for the Fourier spectrum versus frequencies is determined by using the method of least squares. The procedure for determining the spectrum and calculating the attenuation of body waves is basically the same as that used by Hough et al. (1988) and Al-Shukri et al. (1988). The only difference between the present study and that of Hough et al. (1988) is that a Parzen window is used here, whereas a cosine window before spectral analysis was used in the earlier study. It is found that the Parzen window produces a smoother spectrum than the cosine window does. Consequently, it enables the elimination of the need for smoothing the spectrum before fitting the general trend with the least-squares method. The smoothing of the power spectrum also eliminates the effect of variance in the instrument response. Figure 4a shows a sample of one seismogram used for calculating the t* value and S wave window. The truncated data by the time window is shown in Figure 4b. Figure 4c shows the spectra calculated from the data in Figure 4b

LINEARIZED INVERSION FOR S-WA VE ATTENUATION
A least-squares inversion scheme is developed here to determine a three-dimensional image of attenuation within a specific volume of the Earth beneath a local or regional seismic net work. Data used for this technique are the observed t* values calculated from S-waves. The inversion procedure starts by separating the volume of the crust and upper mantle beneath Taiwan into a number of layers, each layer being subdivided into a number of rectangular blocks; in each block both the velocity and a.!\enuation are assumed to be constant. The attenu ation within each block is denoted as Q i , where j is an index of a certain block. The discretization permits the calculation of the initial t* for each ray _ \hat traverses the starting model. For a crustal model with discrete values of velocity and Q i in each block, Equation (2) becomes: where Q1l is the unknown value of attenuation in each block For N observations, a set of equations can then be formed as : where oQ j 1 is the unknown attenuation perturbation ( Q;;6 -Qj -l ) for the jth block through which rays have traveled for a particular period of time. Note that tii in Equation (8) is the theoretical travel time in the jth block from the ith observation. Equation (8) can be rewritten in matrix form as: where Y is the Ot * vector, Tis the travel-time matrix, and q is the unknown vector of attenu ation. A stochastic least-squares solution to Equation (8) is given by where q is the approximate solution of q 'e 2 is the damping parameter with a positive value, and I is the identity matrix. Equation (9) is an overdetermined system due to the expected errors in the observations and t* ab useing the spectral decay technique. The stochastic least squares solution is preferred because it is more stable than the undamped solution. The resolu tion matrix for the solution is given by: (11) and the covariance matrix by: where O' � is the variance of error in the data. The standard errors of the solution are equal to the square roots of the diagonal elements of C. Here, the initial 9 2 is set to 10. Whenever an inverted Q8 are negative, adding a larger 0 2 to the corresponding Q block is needed until that inverted Q8 become positive. But the larger the 0 2 is set, the more overestimated are the Q8 values. The average of e 2 in this study is approximately 1000. The above scheme allows for the final solution to be obtained interactively and for the examination of the convergence of the solution at each step. The essential step of the forward solution for setting up Equation (8) is the calculation of the travel-time matrix. The travel-time traced by a ray in each individual block is calculated by using a ray-tracing algorithm in which blocks are assumed to be homogeneous and rays within these blocks are characterized by approximately straight-line paths. Thurber and Ellsworth (1980), using numerical integration through a heterogeneous medium, obtained an excellent approximation to the actual travel time. In the present study, a numerical integration procedure similar to the one used by Thurber and Ellsworth (1980) is employed to calculate the travel time. First, the ray is classified into a number of segments having equal lengths. Then, the travel time for each segment is calculated, and the times of all segments belonging to the same block are added together to form the total travel time in that block. If a ray hits a block or a layer boundary, the length of the segment is modified such that it is the distance between the end of the last segment and the boundary. The travel time for this modified segment is calcu lated accordingly. In this study, the S-wave velocity model inverted by Roecker et al. (1987) is used to calculate the travel times. The earthquakes used in this study are relocated using the velocity model mentioned above. The event with an error of relocation (ERH and ERZ) be yond 5 Km is discarded.
The inversion procedure is tested extensively. Both homogeneous and heterogeneous models are first used to generate synthetic data for t*. In order to be certain that the inversion program is correct, some different percentages of error are added to t* and the inversion is performed to obtain Q values. The velocity and Q model for testing is shown Figure Sa. The numbers denote the block number. Its corresponding velocity and Q values are shown within the parentheses. The "true" travel times and "true" t* can be calculated from this model. Subsequently, some variances are added to the "true" t* and "true" travel times to examine the amounts that attribute to the inverted Q values. Figure 5b shows the results. It is found that the difference between the inverted and the "true" Q is proportional to the percentage of input errors. The results indicate that the inversion program here is proper. These data are then inverted by using a number of different starting models. In the foregoing calculations, the initial value of Q8 is set at 0. Then the average Q8 values are estimated and used to check the dependence of the results for initial values. An average Q8 of 0.0005 is obtained. The results indicate that the initial values do not influence the evaluated Q8 values to any great extent. In all cases, with a sufficient number of observations, solutions using error-free synthetic data  adequately match with the expected solutions. The errors of the final solution are usually less than the level of the data errors.

RESULTS AND DISCUSSION
The Qs values for each block are inverted. The distribution of the Q8 values obtained in each block is shown in Figure 6. The boundaries of each block in every layer are also shown to the right of the same figure. Generally speaking, the standard errors in Q8 are less than 10 (for an example, see Figure 6). To investigate the characteristics of the Q8 structures, the contour lines for different layers are also constructed and presented in Figure 7. The earthquakes lo cated in each layer are also plotted in the same panel with the Q8 structure. For the top layer (depths of between 0 and 5 km), the Q8 values range from 32 to 1471 (Fig. 6a). About 66% of the Q8 values are less than 200. Two Q8 highs are found: one at the central part of the Central Range, and the other near the Taichung Basin, The lateral variation of attenuation in the top layer indicates that the alluvial layer, which covers the southwest part of the region, seems to have little effect on the Q5 values. It is known that the thickness of alluvium gradually increases to the southwest and reaches a maximum thickness of about 1 km in the extreme southwestern region of Taiwan. However, if the alluvium has a measurable effect on Q8 a trend of decreasing Q8 values should be parallel to the increase in thickness.
In fact, however, Figure 6a shows no such trend. Moreover, the model shows no differ ence in attenuation between the northern part and the southwestern part of Taiwan where the alluvium is about 1 km thick. This indicates that the alluvial cover has only a very small effect on the Q8 values in the blocks. Earthquakes with magnitude greater than 4.0 located in this layer are plotted in Figure 7a. About 85% of the events occurred in zone :::: 200 of the Qs values ( Figure 7a). In high Q8 zones (greater than 600), only 5% of the earthquakes occurred.
In layer 2 (depths of between 5 and 10 km), the Q8 values range between 44 and 805 ( Figure 6b). About 50% of the Q8 values are less than 180. The average for this layer is 150.
The lowest value occurs in a region of intense earthquake activity, while the highest values -occur in the outer blocks far away from the active region. Four local Q8 highs are in alignment with the north-south direction parallel to the Coastal Range. More than 96% of events located in this layer are in Q8 zones lower than :::; 400 (Figure 7b).
In the depth range of 10 to 15 km (layer 3), the Q8 values range from 88 to 649 (Figure 6c) with an average of 230. About 55% of the Q5 values are less than 200. The standard error of each block is slightly less than 5. About 60% of the earthquakes are located in zone S400 of Q5• While 30% of the events occur located in Q8 zones lower than 200 (Figure 7c).
In the depth range of 15 to 25 km (layer 4), the Q8 values are distributed from 58 to 1200 (Figure 6d). With approximately 25% of the blocks having Q8 values of less than 200. About 55% of the blocks have Q8 values between 400 and 500. The average of this layer is 300. The average standard error for this layer is 4. Two significant Q8 highs align parallel to the Central Layer 1 (0 -5 km) Layer 2 (5 -10 km) Range. There is a Q8 low off the Ilan Plain. The other side of the Ilan Plain is Q8 high. Almost 50 % of earthquakes occurred in 600 of Q8 zones (Figure 7d) . In the depth range of 25 to 35 km (layer 5), the Q8 values range between 55 and 925 ( Figure 6e) and reveal an average of about 400. About 50% of the blocks have Q8 values less than 400. The standard error is slightly less than 6. The Q8 values increase gradually from northwest to southeast. Three significant highs parallel to the eastern coast of Taiwan. The local Q8 low off the Ilan Plain still exist. More than 93% of earthquakes occurred in Q8 zones lower than 600. About 80% of earthquakes were located in Q8 zone between 400 and 600 In the depth range of 35 to 50 km (layer 6), most of the Q8 values are greater than those in the upper layers, with a range of 140 to 263 1 (Figure 6f), and an average of about 500. About 43 % of the blocks have Qs values less than 600. Only 10% of the blocks have values less than 300. The standard error is slightly less than 5. Th ree local Qs highs can be found in this layer.
They are located near Hsinchu, Tai tung and the off the Ilan Plain. The local lows off the Ilan Plain in layer 4 and 5 are repla\ed by Q5 local high. About 7 5 % of the earthquakes occurred in Qs zones lower than 800 (Figure 7f). In the depth range of 50 to 75 km (layer 7), the Q8 values range between 135 and 2066 ( Figure 6g). The Qs values on the land are slightly greater than 400. The average of this layer is about 550. About 60% of the blocks have Q5 values less than 400. The standard error is slightly less than 2. One local high almost totally covers the Puli Basin. The local Q6 high off the Ilan Plain still exist. 60% of earthquakes occurred in Q8 zones lower than 800 (Figure 7g).
In depths greater than 75 (layers 8 and 9), the lateral changes in the Q5 values are not considered because only a few regions of the half space were passed by rays.
In general, low Q8 may be caused by weakness and fracture. Obviously, in the Taiwan area higher seismicity is mostly associated with low Q8 in the upper crust. However, in the upper mantle, higher seismicity has high Q8 • This author suggests that in the low Q8 areas, the material may be in a partial melting state, which always causes flow and the elimination of the possibility of earthquakes.
The standard errors shown in Figure 6 are estimated by the regression deviation of t*, the variation of the velocity model and the variance matrix in the inversion procedure (see Equa tion (12). As presented in Equation (2), the Q values are inverted by determining the travel times and velocities. The velocity model used in this study is assumed to be "true". Actually , the accuracy of the velocity model directly affects the ray paths, the travel times or both of them. The variations in travel times are not proportional to the variations in the velocity model. Therefore, to verify and estimate the variance contributed by the velocity model is not easy. It (Continued) means that the uncertainty of Q5 cannot be judged directly from the variance of the velocity model.

CONCLUSIONS
The results of the present study indicate that Q5 structures in the Taiwan area have a good correlation with seismicity. In the upper crust (at depths between 0 to 15 km), the regions of lower Q5 are consistent with higher seismicity. However, in the lower crust and upper mantle (between 15 to 75 km), the regions of higher Q5 values are consistent with higher seismicity.
Liang for some drawings. This work was supported by the National Science Council, ROC under Grants NSC84-21l1-M-003-009.