Prediction of Peak Ground Acceleration in Southwestern Taiwan as Revealed by Analysis of CHY Array Data

1 Center for General Education, Ching Yun University of Technology, Chung-Li, Taiwan, ROC * Corresponding author address: Prof. Jen-Kuang Chung, Center for General Education, Ching Yun University of Technology, Chung-Li, Taiwan, ROC; E-mail: jkchung@cyu.edu.tw Empirical regional attenuation relationships for the amplitude of S and Lg waves were achieved by performing regressions over a large number of accelerograms, collected from TSMIP stations in southwestern Taiwan over a ~7-year period. From estimations of effective ground-motion duration, utilized to characterize the source term and path effects, a piecewise continuous duration function relating to hypocentral distance was derived. For simplicity, Brune’s ω -square model was used in this article. Assuming a high-frequency approximated decay parameter of κ0= 0.03 sec and static stress drop ∆σ = 100 bars for the averaging characteristics of the source spectra, the attenuation curves of peak horizontal accelerations that relate to the hypocentral distance for the given moment magnitudes were effectively obtained using random vibration theory. However, the current model cannot be arbitrarily generalized for the prediction by individual specific earthquake and does not include extended faults based on the source assumption applied. Due to sensitive variation of peak ground motion predictions (strongly associated with the stress drop and attenuation factors), it is suggested that the source parameters and the sites chosen should be localized into several subsets (as far as this is possible) through analyses of historic earthquakes and strong-motion records. In practice, this work should be valid in predicting peak ground motions well when simulating a characterized event. (


INTRODUCTION
The sources of damage caused to man-made structures due to strong ground motions by earthquakes are complicated.Besides design and quality of construction, a clear understanding of ground response at soil sites during an earthquake is helpful in the prediction of damage distribution.For instance, serious damage to buildings relating to ground motion followed a particular pattern after the 1999 Chi-Chi earthquake (M L 7.3) in Taiwan (Shin et al. 2000).As shown in the accelerograms recorded on the free ground surface, stronger ground shaking appears clearly at regions well beyond the vicinity of the Chelungpu fault, for example, the Taipei metropolitan area and the southwestern region.This phenomenon indicates that local geological conditions are possibly the most critical factors in the analysis of ground motion.However, nonlinear amplification at sediment sites appears somewhat unpredictable and is more pervasive than seismologist once though (Aki 1993;Beresnev and Wen 1996).This difficult problem implies that in practice, a forward simulation of ground-motion response ought be an efficient modeling method without foreknowledge of every individual effect on a motion.
In waveform modeling of ground motion, a deterministic approach is generally applied by fitting shape and amplitude of every primary pulse in the seismograms.Such methods have been successfully used for long-period waveforms recorded at regional distances.Another approach, however, must be utilized for describing high-frequency ground motion over shorter distances, since heterogeneous earth and non-uniform fault rupturing make seismograms more complex.To establish such relationships with the degree of damage caused by earthquakes, the peak ground acceleration (PGA) induced by seismic waves with a frequency range of 1~10 Hz is generally simulated by statistical method instead of forward theoretical or numerical modeling.When the available recordings are sufficient in number, in point of view of statistics, regression techniques are used to predict peak ground motion in a target region.For example, Tsai and Bolt (1983) derived the attenuation forms of peak ground acceleration and velocity using SMART-1 array data recorded from six moderate earthquakes.Tsai et al. (1987) deduced the attenuation characteristics of peak ground acceleration from accelerograms recorded by the SMA-1 strong-motion network spread over the entire Taiwan area.Based on the accelerograms recorded on free-field hard-rock sites from 49 earthquakes, the attenuation relations in Kanai's form (Kanai 1961) and Campbell's form (Campbell 1981) of peak ground acceleration were also proposed by Ni and Chiu (1991).In their results, the attenuation behavior of peak amplitude with distance may have poor correlation with site effects of seismic waves because data used are limited in specific condition (almost classified to the hard-rock sites of the SMA-1 network, otherwise, to a very localized spread of sites as in SMART-1 array).However, with the advantages in quality and number of strong-motion data collected by the Taiwan Strong Motion Instrumentation Program (TSMIP) network operated by the Central Weather Bureau since 1992 (Shin 1993), Huang (1995) regressed the attenuation of peak acceleration by adopting a form of Campbell's attenuation using a total of 526 digital accelerograms.Using more available strong-motion recordings, other studies (e.g., Tan and Loh 1996;Shin 1998;Liu 1999) on the regression of peak ground motion for the whole Taiwan area, with particular reference to local-site response in metropolitan regions, have been proposed over the past decade.
An alternative approach for performing the prediction of ground motion needs a source model and a predictive relationship, which allows the estimation of the specific ground-motion parameter as a function of magnitude, distance from the source, and frequency.An empirical attenuation function can be obtained by regressing a relatively simple functional form containing an average crustal attenuation parameter and a geometrical spreading function.The dynamic analysis of structures requires an estimate of signal duration, and this is the reason why the measures of ground motion based on some average of the motion over time are sometimes used in seismic hazard.Adopting the ways (outlined in the following context) described in the papers of Raoof et al. (1999) and Malagnini et al. (2000), an estimate of regional attenuation was carried out in the present study using 158 sets of strong-motion data collected by the CHY array of TSMIP network during the period 1993 to 2000.The advantage of their method over the regression analyses is to take into account the duration parameter as a function of frequency and distance, which is more effective in structural engineering applications.Then a statistical tool called the random vibration theory (RVT) (Cartwright and Longuet-Higgins 1956) was used to estimate the peak motion of a random time history, given the empirical attenuation parameter, source spectrum, and its duration in time.Such a stochastic method is generally based on the assumption that the amplitude of ground motion at a site can be specified with a random phase spectrum such that the motion is distributed over a duration related to the earthquake magnitude and to the distance from source (Boore 1983(Boore , 1986)).The purpose of ground-motion prediction is to establish an adequate building code for earthquake engineering and hazard assessment based on possible peak ground motion at a given distance from a scenario large earthquake.Comparing the predictive peak motion to the observations, local site effect can also be examined for further modification in different regions.

REGIONAL ATTENUATION RELATIONSHIP
Ground motion at a particular site is influenced by source, travel path, and local site conditions.Following the statements in Raoof et al. (1999), Malagnini et al. (2000), and Malagnini and Herrmann (2000), the logarithm of the peak values of ground motions is represented in a general form for a predictive relationship as the following: log ( , ) log ( , ) log ( , , ) log ( ) where a r f ( , ) is the peak amplitude of the narrow bandpass-filtered ground acceleration recorded at the hypocentral distance r.The term E r f ref ( , ) represents the seismic source scaled to the reference distance r ref , which can be arbitrarily chosen according to the distribution of data used P r r f ref ( , , ). describes the crustal attenuation in the region, and S f ( ) is a site term.The instrument factor is eliminated when the motions are recorded by seismometers with the same response in the range of the frequency of interest.
By separating the spectrum of ground motion into source, path, and site components, the model can be easily modified to account for specific situations or to account for improved information about particular aspects of the model.To obtain the empirical estimate of the path term in equation ( 1), a coda normalization technique (Aki 1980;Frankel et al. 1990) is used to eliminate the source and site terms by taking the ratio of the peak amplitude carried by the S or the Lg waves in the filtered time histories to the rms amplitude of the stable seismic coda.The essence of this method is that the site and excitation effects are identical for both of the S waves and coda waves.Therefore, the coda-normalized amplitude can be written in the form: log ( , ) log ( , , ) log ( , , ) where A reference time window within the coda waves is chosen for calculating the rms amplitude in the definition of Parseval's theorem.However, the arbitrariness in the choice of the length of time window is removed by forcing log ( , , ) P r r r f ref ref

=
=0.The preliminary esti- mate of the empirical path term accounts for geometrical spreading, attenuation (combining intrinsic and scattering attenuation), and the general increase of duration with distance due to wave propagation and scattering at each hypocentral distance for each sampling frequency.Based on the results of the coda normalization method, the Q model of the S/Lg waves propagating through the crust will be inverted by giving an empirical geometrical spreading model.

STRONG MOTION DATA SET AND DATA PROCESSING
The data used in this study were collected by the TSMIP network of the Central Weather Bureau.The network consists of over one thousand digital accelerographs deployed over nine metropolitan areas and several tens of buildings for structural safety evaluation (Shin 1993).Most of the free-field strong-motion stations are located on the sedimentary ground of elementary schools.Each station equips with a three-component force-balanced accelerometer, either an A-900 or IDS-3602 in the early stage of this program, and with a 16-bit resolution and a 96db dynamic range, or other types of accelerometers more recently.Through the testing of instrument parameters, the basic specifications of both types of instruments are almost the same, except the sampling rate.
To ensure the sampling quality of propagation characteristics in the target region, the averaged distribution of sites is a basic requirement.A total of 117 TSMIP stations (the first three characters of station code are labeled with "CHY") located in Yun-lin, Chia-yi, and Tainan Counties were chosen for searching available accelerograms (Fig. 1).Their site conditions have been categorized into four classes, namely, B, C, D, and E based on the geologic and geomorphologic data (Lee et al. 2001), and are listed in Table 1.It is obvious that most of the sites at the western flood plain are classified to class E (soft soil).Meanwhile, only four sites located at the Central Range are recognized as class B (rock sites).Subsurface soil properties are critical factors in site effect studies of earthquake ground motion.According to the available P-S logging profiles under the project established by the CWB and National Center for Research on Earthquake Engineering (NCREE), the shear wave velocities of strata in the upper 30 m at some of the TSMIP/CHY stations are ranging from 190 m s 1 − to 230 m s 1 − without strong dependence on the site class.However, it can be expected that site conditions at some stations close to mountains may differ from that of soft sedimentary sites, mainly due to variations in soil properties and thickness.With a criteria of local magnitude greater than 4.0 and focal depth shallower than 50 km, there were 146 earthquakes prior to the 1999 Chi-Chi earthquake to be sorted, which triggered at least one accelerometer of the CHY array within the source-receiver distance of 150 km.Another set of 12 events was selected from the aftershocks of the Chi-Chi earthquake and the October 22, 1999 Chia-Yi earthquake ( M L 6.4) sequences with magnitude greater than 6.0.The epicenters of these earthquakes are shown in Fig. 1.Even though the number of earthquakes used is enough, it is noticed that the azimuthal coverage of propagation direction of seismic waves is mostly concentrated in the eastern side of the target region due to the distribution of seismicity.Figure 2 demonstrates the distribution of hypocentral distances for the selected 2860 three-component accelerograms.All these waveforms have no abnormal spikes caused by unstable instrument output or unknown reasons.The records containing overlapped wavetrains due to multiple events with very close origin time were also abandoned during data selection.Thus, a total of 5720 horizontal acceleration time histories were gathered as the dataset.
Table 1.The site classification of TSMIP/CHY array based on Lee et al. (2001) Each accelerogram was narrow-bandpass filtered around each of the following sampling frequencies of 0.2, 0.25, 0.33, 0.4, 0.5, 1.0, 2.0, 3.0, 4.0, and 5.0 Hz.The peak amplitude of each filtered time history, carried by direct S or Lg waves windowed at a length of 5 sec (as indicated in some records shown in Fig. 3), was picked out.Meanwhile, the time window for calculating the rms amplitude of coda waves was 10 sec.Using the coda normalization technique described in the previous section, the attenuation functions for peak amplitude at each sampling frequency were estimated by choosing a reference hypocentral distance of 40 km to remove the effect due to the arbitrariness in the choice of coda duration.The results, demonstrated in a log-to-log space, are shown in Fig. 4. To avoid biases from large outliers, an L1norm minimization (Bartels and Conn 1980) was performed.The medians (denoted as black diamonds at some specified distances in the figure) can be considered as the nodes of the piecewise empirical attenuation function.It is noted that the logarithm of peak amplitude of ground motion is normalized to zero at the reference distance for each frequency.By visual inspection of the curves, the peak motions recorded at larger distance present relatively slower decay with the hypocentral distance.This phenomenon can be caused by the discrepancy of attenuation characteristics between the S (or Lg) waves and the coda waves, and perhaps by site amplification, although the assumption for applying coda normalization method has been taken into account.A multi-parameter inversion can be used to extract the site term by modeling equation (1).However, unless detailed information about the geology underlying the sites is well investigated, the ground-motion scaling relationships obtained from further regression are similar to their estimates computed by using the coda normalization technique (Malagnini and Herrmann 2000).Therefore, no more inversions were carried out to obtain empirical attenuation functions in this study.The effects produced by site condition and radiation of seismic sources are not explicitly separated and will be merged into the estimations of geometrical spreading and the anelastic models.

PREDICTION OF PEAK GROUND ACCELERATION
A useful and simple method for simulating ground motion is to combine functional descriptions of ground motion's amplitude spectrum with a random phase spectrum, as represented in equation (1).Basically, the essence of the stochastic method is that the motion is distributed over a duration related to the earthquake magnitude and to the distance from the source (Hanks and McGuire 1981;Boore and Atkinson 1987).A rapid way of obtaining the peak motions (response spectra, peak acceleration, peak velocity, peak displacement, etc.) is to use random vibration theory, as proposed by Cartwright and Longuet-Higgins (1956), when given knowledge of the Fourier amplitude spectra of time history and signal duration.It is particularly useful for simulating the high-frequency ground motions of most interest to engineers.In order to use this theory for the extremes of transient earthquake ground motions, the definition of duration used in determining rms amplitude could be calibrated (Boore and Joyner 1984).A computer code (Boore 2000) was used in this study, one which demands several input parameters including seismic source spectrum, a crustal attenuation function (anelastic Q and geometrical spreading models), and a distance-dependent duration function.The empirical anelastic attenuation function is modeled using the following functional form [equation ( 13) in Malagnini et al. 2000 where r is the hypocentral distance; the reference distance r ref is 40 km used here and β = 3.5 km s 1 − (Yeh and Tsai 1981) is the average shear-wave velocity in the source regions.g r ( ) is a piecewise linear geometrical spreading function.For larger distances, the peak amplitude of motion is mainly carried by higher mode surface waves (Lg) and the geometrical attenuation is generally expected to be proportional to 1/ r (Shin and Herrmann 1987).Taking the essence of the expression as equation ( 4), however, it may be possible that the geometrical spreading effect and the coefficient of anelastic attenuation could exchange with each other (Atkinson 1989).Similar to the multi-segment models of geometrical spreading supposed in Atkinson and Mereu (1992) and Atkinson and Boore (1995), when considering dominant wave types at various distances, a bilinear geometrical spreading function shown in Fig. 5   for the Taiwan region (Sokolov et al. 2000) and was also adopted by Roumelioti and Beresnev (2003), based on data averaging of the variety of site conditions from rock sites to soft soils of different thickness.Based on equation ( 4) and the supposed geometrical spreading function g r ( ) , the average Q values can be estimated by giving the distances and the coda-normalized amplitudes at sampling frequencies.Neglecting the abnormal bias at 0.2 Hz, these discrete anelastic attenuation factors were matched with a bilinear function shown as Fig. 6, which has different behavior crossing the frequency of 0.4 Hz.The frequency-dependent Q model, as the input parameter of the simulation using random vibration theory, is: A Qs value of about 200, by analyzing the 1~5 Hz seismic waves propagating through the central Taiwan region, was proposed by Wang (1987).A more detailed Qs model was constructed by Chen (1998), which represents increasing value from 130 to 400 with increasing depth of stratum from the surface to 35 km deep in the of interest frequency band of 4~25 Hz.Basically, there is no obvious discrepancy between the function in equation ( 5) and the above two models.
In addition, distance-dependent ground motion duration is an important function for use of stochastic method simulations.Although the Fourier amplitude spectrum of ground motion is not dependent on duration, duration is a function of the path, as well as the source that is related to the inverse of the corner frequency.The definition of duration of ground motion in this study is given as the width of the time window that limits the 5%~75% portion of the seismic energy following S-wave arrival.From the results performed for each accelerogram at each sampling frequency, the distance-dependent duration function can be represented by a connected series of straight-line segments derived by the L1-norm optimization.Figure 7 shows the durations varied with hypocentral distance for the regional data at different frequencies.Larger scatterings of the durations are observed at lower frequencies, also at larger hypocentral distances, which are induced by shallow geologic conditions, as well as by seismic noise characteristics at specific sites.Again, an L1-norm minimization was performed to obtain the averages at sampling hypocentral distances from the large outliers.The reduction in duration at lower frequency at a distance of 150 km is caused by truncation at the end of seismograms, due to a typically programmed 15 seconds for post-event memory as the ground motion amplitude is lower than the threshold level (in general, 0.2% of the full scale of 2 g) (Liu et al. 1999).The empirical distance-dependent duration function shown in Fig. 8, obtained by averaging the estimates for the frequency range of 1~5 Hz (considered as the major band to perform the recorded peak acceleration), was given for predictive modeling of ground motion decay through random vibration theory.The function, as with that of Atkinson and Boore's (1995), can be matched with a bilinear one, which is represented by the heavy line in the figure .The important parameter required in calculations using random vibration theory is the shape of the source spectrum.The most common model is the ω -square model (Aki 1967;Brune 1970).The scaling of the spectra is determined by specifying the dependence of the corner frequency on seismic moment.Analyses from Huang and Yeh's (1999) study indicate that the source spectra for some events in the Chiayi-Tainan region cannot be explained by the ω -square model.However, the earthquakes used here occurred not only in the southwestern area, but also in central and eastern Taiwan.A typical ω -square model with a single corner frequency was still adopted in this study.
For modeling the shape of source spectra, the stress drop ∆σ of the seismic source is a critical parameter in the prediction of peak ground motion.Previous studies' published results vary widely.From a nonlinear, least-square regression based on the functional form of Brune (1970), Ou and Tsai (1993) calculated an average stress drop of about 150 bars for the events of Chia-yi and Tainan regions.By analyzing the source spectra of S waves, which were determined using records of eighteen earthquakes occurring in the same area with local magnitude of 2.8 ≤ M L ≤ 5.8, an average stress drop of 125 bars was obtained (Huang and Yeh 1999).In their study, a stress drop of about 60 bars was estimated from the spectral level in a lower Fig. 7.The 5%~75% duration distribution for the regional data at different frequencies.By an L1-norm minimization process, the averaged effective durations were computed at different distances indicated by the diamonds.
frequency range, while the value of 600 bars was required to interpret high-frequency amplitudes.
According to the assumption of a point, circular fault model, the pulse width of the sourcetime function is a key factor in the determination of stress drop.Much higher estimates, in the kilobars, have also been obtained using different source-time durations for the 1993 Tapu earthquake ( M L = 5.8) (Shin 1995;Huang et al. 1996).As discussed by Huang and Yeh (1998), the essential differences among these results may be attributed to the pulse width of the source-time function used.Moreover, readings of the low-frequency flat levels and corner frequencies from displacement spectra are usually done by eye, making results subjective and, consequently, strongly affecting the estimation for stress drop (Moya et al. 2000).Recently, Huang et al. (2002) also suggested high dynamic stress drop proportional to the magnitude of the aftershocks of 1999 Chi-Chi earthquake.Their estimates are ∆σ = 991 bars for M L = 6.4 aftershock and ∆σ = 831 bars for M L = 6.0 one, respectively, by calculating the integrals of squared velocities and squared displacements.In fact, high stress drop only implies a rapid slip velocity of dislocation (Atkinson and Beresnev 1997;Beresnev 2001).In the other word,  the stress drop only serves as a proxy for the source radius in the relationship between the radius and the corner frequency (Beresnev 2001).Moreover, most of the previously proposed high stress drops, which were obtained by fitting high-frequency levels of observed and modeled spectra, could probably be exaggerated due to local site effects which reveal in the high frequencies.Consequently, those very high estimations of stress obtained from some special stations and/or earthquakes, even in the same target region, are not suitable requirements for developing a generic attenuation model.Through the Lg spectral analyses of the short-period and strong motion records from 865 earthquakes (1.28 ≤ M L < 5.88) in Taiwan, for the cases of seismic moment M 0 greater than 10 22 dyne-cm, a determination of stress drop in the range from 10 to 100 bars was derived by Chiang (1994), which is comparable with the value for California's earthquakes derived by Brune (1970).In Chiang's result, a bilinear relationship between local magnitude and seismic moment was supposed by adding seismic moments of 54 earthquakes ( M L ≥ 5) published in the PDE catalog.Among these values with extreme diversity, a typical, constant stress drop of 100 bars was quoted in this article because the source distribution in Chiang (1994) is similar to the one selected here.
For most applications, it is advisable to represent the effects of the path by a simple function, as in equation ( 4) that accounts just for geometrical spreading and the attenuation.However, in the strictest sense, the modification of seismic waves by local site conditions is part of the path effect.Unless in the situation that the specific-site effect including nonlinear behavior on the ground motion can explicitly be determined, many of the simulations from the stochastic method are intended to be used for the prediction of motion at a generic site.In such cases, a simplified diminution function, exp[ ] −πκ 0 f , associated with the path-independent loss of energy can be used to describe the frequency-dependent modifications of the seismic spectrum (Anderson and Hough 1984).Accordingly, in Boore's (2000) computer code, the highfrequency approximated decay parameter κ 0 , which describes the high-frequency roll-off of source spectra propagated to site, is a required input for representing a distance-independent regional average of the attenuation in the shallow crust.Using the recordings collected from the SMART 1 array in Lotung, Taiwan, the kappa values, which depend on site, distance, and source, were determined to be in range of 0.03~0.07sec (Tsai and Chen 2000).Owing to a lack of detailed study at the target area, a kappa value of 0.03 sec is quoted based on the result of Tsai and Chen's (2000), although this may be more suitable for harder sites.This value may represent the regional average of the combined effects for generic site amplification and anelastic attenuation.
Using random vibration theory with the inputs of empirical parameters in the previous description, the estimates of peak ground horizontal accelerations for moment magnitudes M W = 4.0, 5.0, 6.0, and 7.0, respectively, are shown in Fig. 9.The attenuation curves of peak ground acceleration (denoted as g) decay with increasing hypocentral distance, ranging from 5 km to 150 km based on data selected.The result indicates that in practice, for a specified earthquake with a shallow focal depth, a predictive curve can be obtained by giving the moment magnitude of the event.

DISCUSSION AND CONCLUSIONS
The empirical regional attenuation relationship was obtained by regressing a relatively simple functional form containing an average crustal attenuation parameter and a geometrical spreading function.The result was performed by a large number of time histories recorded from the TSMIP network during the period 1993 to 1999.By using a coda normalization technique to eliminate the source term and the site term in equation ( 1), which is done by taking the ratio of the peak amplitude carried by the S or the Lg waves in the filtered time histories to the rms amplitude of the stable seismic coda, the logarithm of band-pass filtered peak amplitude was normalized to zero at an arbitrary reference distance.After the empirical attenuation functions for the sampling frequencies were produced through some averaging procedures, the frequency-dependent crustal quality factor Q(f) was modeled from the functional form represented in equation ( 4) by introducing the geometrical spreading function used in Sokolov et al. (2000) and Roumelioti and Beresnev (2003).
Indeed, the geometrical spreading function herein gives a possible pattern to describe the interference phenomenon between the body waves and the higher-mode surface waves.The effect of propagation can be clearly reflected in the predictive curves, especially in distance larger than 50 km, shown as the heavy lines in Fig. 9.Even with an appropriate amount of deviation existing in different models of the geometrical spreading function, the effect caused by this discrepancy could be compensated with the modifications in the anelastic attenuation function (Atkinson 1989).That was the case in this study.
From calculations of the effective ground-motion duration, a piecewise continuous duration function dependent on the hypocentral distance was derived.The duration is defined herein as the width of the time window that limits the 5%~75% portion of the seismic energy following the S-wave arrival.The source model used here is a typical Brune's one.The high-frequency approximated decay parameter κ 0 = 0.03 sec, which describes the high-frequency roll-off of source spectra propagated to the site, was used to represent a distance-independent regional average of the attenuation in the crust.Referring to related articles, a stress drop of 100 bars for moderate earthquakes was assumed for controlling the high-frequency level of the source spectra.Liu (1999) analyzed about 1000 TSMIP accelerations recorded from 16 shallow earthquakes ( M W = 3.5~6.3) in Chiayi-Tainan area to regress the attenuation relation of peak ground acceleration as the following form: where X represents the hypocentral distance in km.Making the comparison between his result and the prediction obtained in this study shown in Fig. 9, the estimations of this study are smaller than Liu's (1999) results, except in the distance of 60~150 km.Considering focal depth as the major parameter in interpreting ground motion related to earthquake energy, Chang et al. (2000) derived two attenuation laws for shallow crustal earthquakes and subduction earthquakes, respectively, using a two-step stratified regression.
where Dp is the focal depth, and De the epicentral distance in km.In Fig. 9, the dashed curves expressed by equation ( 8) developed by Chang et al. (2000) are for cases of a focal depth of 10 km.From the comparison in the figure, the prediction of PGA in this study is the smallest value among currently discussed attenuation models.A possible reason for underestimation of ground motions is that short-period surface waves with dominated amplitude were inattentively regarded as generic coda waves in the normalization procedure.For example, distinct short-period surface waves were generated by the 1993 Tapu earthquake and made prolonged accelerograms recorded at sedimentary sites of some tens of kilometers away from the source (Chung and Yeh 1997).Such phenomenon was frequently observed at soft-soil stations deployed in the Chia-Nan plain when a very shallow earthquake occurred.A way to avoid this effect is to exclude such sites through comprehensive investigations of site response.Such an endeavor will be pursued in improvement of prediction saved for further investigation.The residual values, which are defined as the natural logarithm of the ratio of the observed value to the predicted one, are illustrated as a function of hypocentral distance (Fig. 10).The 90% of error estimations are confined between the two dashed lines symmetrically above and below the mean value of 0.38, which might be due to the site amplification.In general, the dispersive distribution is similar to the pattern shown in Fig. 3 of Chang et al. (2000).To find the sources introducing local site effect on the observed deviation, the misfits were grouped into four sets based on various site classes B, C, D, and E. As a result, the average residuals for both classes B and C, are very close to zero.Nevertheless, the error means for classes D and E are 0.41 and 0.52, respectively.These evaluations imply that site amplification of a factor of ~1.6 should be considered a correction factor upon predictive acceleration, when the following are satisfied: (1) the station is located on a soft site, and (2) the predictive relationship herein is examined to be suitable only for hard sites (e.g., class B and class C).
Figure 11 shows the comparisons between the observed PGAs and the various attenuation laws for the three moderate-to-large earthquakes that occurred in or near the study area.The predicted curves in this study, for both the October 31 1995 earthquake (M W = 5.03) and the July 17, 1998 Ruey-Li earthquake ( M W = 6.05), reveal being more suitable than the regressions obtained by the other two models.However, for the September 25, 1999 earthquake ( M W = 6.65), one of the major aftershocks of Chi-Chi earthquake, all the predictions are obviously underestimated compared to the data shown in Fig. 11c.An important conceptual statement is that the averaged prediction does not match well with each individual event with specific focal mechanism.Therefore, an understanding of source characteristics based on historical events in the specific region is necessary to build an appropriate corresponding model.
In stochastic simulations, the acceleration spectrum is modeled based on the empirical parameters associated with source characteristics and propagation behavior.Again, to highlight the influence of local site condition on prediction of ground motion, the Fourier acceleration spectra transformed from several accelerograms are shown in Figs. 12, 13, and 14, respectively, for the three discussed earthquakes.For each event, a set of 6~12 horizontal time histories (i.e., 3~6 stations depending on the available data) with almost the same distance (+/-1 km) and with the same category of site class, were analyzed to represent the averaging spectra, for comparison with the model.The length of accelerogram for calculating Fourier spectrum is 30 seconds.Through the truncation in time domain, unusual high amplitudes of spectra at low frequencies (< 0.2 Hz) can be obtained in some cases, which are not shown in these figures.The upper frame of each figure is the results for site class C (sometimes including records from site class B), while the lower frame is for site class E. In the comparisons for the October 31, 1995 earthquake, the shapes and the level of spectra between the prediction and observation are similar (Fig. 12).But for the 1998 Ruey-Li earthquake, the modeled amplitudes in the low frequencies are significantly underestimated for the sites of class E (Fig. 13b).The target distance is 35 km which was enough to generate predominant short-period surface waves due to a moderately large earthquake with a focal depth of only 2.8 km.However, this phenomenon is not shown in the sites of class B and C, even at a distance of 30 km.In such a situation, it can be believed that the soft site conditions play an important role in the site amplification of ground motion.When looking at the spectra for sites with hypocentral distance over 50 km and earthquake magnitude of 6.65 (Figs. 14 a, b), major differences are revealed in the frequency of 0.3~3 Hz, for both hard and soft sites.The amplification around 1 Hz almost reaches an order of 10 comparable with the PGA ratio (Fig. 11c).It would appear that soil condition seems not to be the only factor making such strong amplification.Back to    the essence of the quoted stochastic method, it does not include any phase effects due to propagating rupture and wave propagation routing to the site.Therefore, fault-normal effects, phase differences over horizontal distances, spatial correlation, directivity, etc. are not captured by the simulated motions.Some other efforts should be necessary to solve problems such as these: for example, significant amplification within a frequency band as this is of significant importance in earthquake engineering.
Stress drop also characterizes the level of high-frequency amplitude in the source spectrum and is proportional to the seismic moment of an earthquake.The use of a higher stress drop parameter will lead to larger peak ground acceleration predictions via random vibration theory, since peak acceleration is almost characterized by the high-frequency motion.An apparent limitation is that models based on the stochastic method are fundamentally point-source models.For applications extended to earthquakes with larger magnitude, the effects of a finite-fault averaged over a number of sites distributed around a fault can be captured using the closet distance to faulting as the source-to-site distance and a two-corner source spectrum (Atkinson and Silva 2000).
In conclusion, it is important to check the seismic spectrum to ensure that the assumed source model does not deviate too far from the desired one.Although, stochastic modeling is considered an efficient and powerful tool in the prediction of peak ground motion, it is suggested that the source parameters and the sites chosen should be localized into several subsets as much as possible through analyses of historic earthquakes and strong-motion records.

Fig. 1 .
Fig. 1.Epicenter distribution of analyzed earthquakes (denoted by open circles) and the locations of CHY-stations (indicated by different symbols for the soil classes classified by Lee et al. 2001) of the Taiwan Strong Motion Instrumentation Program (TSMIP) Network used in this study.The radius of the denoted circle is proportional to the local magnitude (M L ) ranging from 4.0 to 7.1.

Fig. 2 .
Fig. 2. Distribution pattern of local magnitude versus hypocentral distance of the 2860 accelerograms recorded in the 158 earthquakes.Most of the data were recorded at the hypocentral distance range of 5~150 km.

Fig. 3 .
Fig. 3. Examples of horizontal acceleration time histories.The S phases as well as the Lg waves are marked.The time windows indicating Lg waves are 5-sec long adopted in this study.R and M represent hypocentral distance and moment magnitude, respectively.The character in parenthesis following the station code represents the soil class based on the classification of Lee et al. (2001).

Fig. 4 .
Fig. 4. Regressive attenuation relationships of the coda-normalized peak acceleration amplitudes carried by S and Lg waves at the sampling frequencies of 0.2, 0.25, 0.33, 0.4, 0.5, 1.0, 2.0, 3.0, 4.0, and 5.0 Hz.The small crosses denote the logarithms of normalized amplitudes, which are normalized to zero at the reference hypocentral distance of 40 km.The solid diamonds are obtained by using a L1 minimization algorithm to form a piecewise empirical attenuation function.
was used herein.That is in the form 1/r b , where b = 1.0 for r < 50 km, b = 0.0 for 50 ≤ r < 150 km.This model was used

Fig. 5 .
Fig. 5. Geometrical spreading function used in this study.It is composed of two linear functions with various slopes for charactering distance-dependent decay of seismic waves.

Fig. 6 .
Fig. 6.Q model derived from the normalized amplitudes of S and Lg waves based on the attenuation functional as equation (4) and the supposed geometrical spreading function shown in Fig. 5.The frequency-dependent Q function (denoted by the solid line) is characterized by equation (5), which is determined by eye from the results for each filtered frequency (shown as the diamonds).Note the unusual high Q value at 0.2 Hz was ignored in modeling of the Q function.

Fig. 8 .
Fig. 8. Empirical distance-dependent duration function.This bilinear function (indicated by solid line) was obtained by matching the averaged estimates of all sampling frequencies at different distances.

Fig. 10 .
Fig. 10.Misfits between the predicted peak acceleration (Ap) and the observed value (Ao).The 90% of the error estimations are confined between the two dashed lines symmetrically above and below the mean value (solid line).The evaluations were made according to various site classes B, C, D, and E, and all sites (left upper frame).

Fig. 11 .
Fig. 11.Comparisons of the different attenuation curves with the observed peak accelerations (indicated by the open squares) for (a) the October 31, 1995 earthquake ( M W = 5.03),(b) the July 17, 1998 Ruey-Li earthquake ( M W = 6.05), and (c) the September 25, 1999 aftershock of Chi-Chi earthquake.The focal depth parameters required in the model proposed by Chang et al. (2000) were determined from Central Weather Bureau, Taiwan.

Fig. 12 .
Fig. 12. Modeled Fourier acceleration spectra (gray heavy dashed lines) and the recorded ones from the sites belonged to (a) class B and class C with the hypocentral distance of 30 km, and (b) class E at a distance of 25 km for the October 31, 1995 earthquake ( M W = 5.03).

Fig. 13 .
Fig. 13.Modeled Fourier acceleration spectra (gray heavy dashed lines) and the recorded ones from the sites belonged to (a) class B and class C with the hypocentral distance of 30 km, and (b) class E at a distance of 35 km for the July 17, 1998 Ruey-Li earthquake ( M W = 6.05).

Fig. 14 .
Fig. 14.Modeled Fourier acceleration spectra (gray heavy dashed lines) and the recorded ones from the sites belonged to (a) class B and class C with the hypocentral distance of 70 km, and (b) class E at a distance of 80 km for the September 25, 1999 aftershock of Chi-Chi earthquake (M W = 6.65).
That is