Hybrid Ground Motion Simulation for the 2013 ML 6 . 4 Ruisui , Taiwan Earthquake

The 2013 ML 6.4 Ruisui earthquake struck a seismic gap around the northern part of the Longitudinal Valley where rare earthquakes have occurred for at least two decades. Seismic data with different frequency bands connote a diverse scale of source characteristics. The slip model derived by a previous study using long-period seismic data provided us with the preliminary earthquake source properties in low-frequency energy radiation. However, the high-frequency part of the seismic energy needs to be considered for the strong motions in the comprehensive frequency band. As a case study, we carried out broadband strong motion simulations through the hybrid scheme, which simulated waveforms combining the low-frequency motions (frequency-wavenumber integration method) and high-frequency motions (empirical Green’s function method) to reveal the strong ground motion and source characteristics of the 2013 Ruisui earthquake. The results show that even though the local structure is complicated this hybrid simulation can effectively rebuild the overall broadband strong ground motion characteristics for the 2013 Ruisui earthquake. Additional study issues are required to better understand the nature of localized high peak ground acceleration and the related seismic hazards for future earthquakes.


InTRoducTIon
The NNE-trending 150-km-long Longitudinal Valley fault (LVF), which overlies the suture zone between the Philippine Sea plate and the Eurasian plate, is an active high-angle thrust fault (Yu and Kuo 2001;Hickman et al. 2002).On 31 October 2013, the M L 6.4 Ruisui earthquake struck the northern segment of the LVF (Fig. 1) and caused strong shaking (> 400 gal) recorded in the northern LVF.In the past the Longitudinal Valley has been struck by several large events, e.g., the 1951 Longitudinal Valley earthquake sequence and the 2003 Cheng Kung earthquake (Ching et al. 2007;Chen et al. 2008).Based on the background seismicity from the past two decades, it is noted that the epicen-ter of the Ruisui event is located in a seismic gap (Pu 2013), where there is an unclear segment exists between the locked zone in the north and the creeping fault in the south (Chen et al. 2008).According to the Global Centroid Moment Tensor (CMT), Central Weather Bureau (CWB) CMT, and the Broadband Array in Taiwan for Seismology (BATS) CMT solutions, the 2013 Ruisui event shows the focal mechanism with thrust motion.The aftershock distribution revealed that the rupture extended to the north along a west-dipping fault plane (Chuang et al. 2014;Lee et al. 2014).Furthermore, Lee et al. (2014) indicated that high peak ground acceleration (PGA) north of the 2013 Ruisui mainshock epicenter might respond to the rupture directivity effect.Seven months later the 21 May 2014 M w 5.6 Fanglin earthquake occurred at the northern end of the 2013 Ruisui event aftershock distribution.
The strong motion record in frequency band of 0.5 -2.0 Hz is more important and interesting for most earthquake engineering applications (e.g., Kawase 1998), and the analyses in the shorter-period range are quite effective for better understanding the source characteristics.Nevertheless, low-frequency band data also provided the principal information for the overall faulting process (< 0.66 Hz, Furumura et al. 2002;0.05 -0.2 Hz, Lee et al. 2014).In this study we used small-event records as the empirical Green's functions (EGFs) and investigated the rupture directivity from the source spectral ratios and estimated the strong motion generation area (SMGA) by fitting the observed ground motion records from the synthetic waveforms.To better build up the knowledge of the potential seismic hazard to the Longitudinal Valley area, we explored source characteristics from simulation results on both high-and low-frequency strong ground motions for the 2013 Ruisui earthquake through a hybrid scheme that simulated waveforms by combining the low-frequency motions (frequency-wavenumber integration method) and high-frequency motions (EGF method).

obSERvEd dATA
The deterministic waveform inversion has been broadly and successfully used to obtain the source process and slip model, which contains low-frequency (up to approximately 1 Hz) source information (e.g., Hartzell and Heaton 1983).Thus, using the inverted slip distribution as the input source model is helpful for low-frequency strong motion forward simulation.On the other hand, the EGF method (Irikura 1986) is one of the developed techniques used to analyze the high-frequency properties of the source after cancelling out the propagation path effects and site amplification.The strong motion records are useful in source modeling for estimating strong ground motion as well as investigating the relationship between ground shaking and seismic damage.Thus, we adopted free-field strong motion data maintained by CWB, and a total of 10 stations (Fig. 1) were selected based on the waveform quality for both the mainshock and a small event which was chosen as the EGF of the mainshock.The epicenters and fault plane solutions for the 2013 Ruisui earthquake and the EGF were determined by CWB and BATS, respectively (Table 1).

STRonG MoTIon SIMuLATIon
The broadband strong motions are simulated through a hybrid scheme that combines the low-frequency motions (0.08 -0.5 Hz) from the frequency-wavenumber (F-K) integration method (Zhu and Rivera 2002) and high-frequency motions (0.5 -10 Hz) from the EGF method.

Simulation of Low-Frequency Motion
The 2013 Ruisui earthquake was located near the east coast of Taiwan; therefore, most near-field stations were to the east of the mainshock (Fig. 1). Lee et al. (2014) used teleseismic waveforms, geodetic data, and strong motion records to build a detailed source rupture model of the 2013 Ruisui event.We adopted their inverted slip distribution as the input source model for generating low-frequency strong motion.In this study our goal is not to obtain the best fit model but to investigate the main source properties and reproduce the overall characteristics of the observed ground motions.Thus, following the fault setting of Lee et al. (2014), we used the F-K method (Zhu and Rivera 2002) with a 1D velocity model modified from a 3D velocity model (Kuo-Chen et al. 2012) to compute the synthetic Green's functions.
Figure 2 shows a comparison of the observed and synthetic acceleration waveforms, both with bandpass filters at 0.08 -0.5 Hz.It is clear that the amplitudes of the synthetics are underestimated at these near-field stations and this underestimation might be caused by the absence of soft surface layers (Kamae et al. 1998).Because CWB measured the logging P-and S-wave velocities in the shallowest 300 m at the downhole seismometer station HWA (Shin et al. 2013;Wen et al. 2014a), we adopted the near-surface velocity structure  2).The synthetics can reasonably explain the main features and amplitudes of the observations, as shown in Fig. 2. We noticed that the synthetics of the station between two faults (HWA054) showed much higher amplitudes than the observed records, which might have been caused by the simplification of the velocity structure.This improvement suggests that the adopted near-surface velocity structure is necessary when calculating the near-field Green's functions.

Simulation of High-Frequency Motion
The EGF method has been developed and applied to investigate seismic source properties in the high-frequency range up to approximately 10 Hz (e.g., Irikura 1986;Irikura and Kamae 1994;Miyake et al. 2001).Due to the available data and waveform quality, we selected a small event that occurred approximately two hours after the mainshock, with M L 4.6 as the EGF (Fig. 1).Following the source spectral ratio fitting approach proposed by Miyake et al. (1999), the observed source spectral ratio between the mainshock and EGF event was fitted with a theoretical ratio.The theoretical source spectral ratio function, which follows the omegasquared source model of Brune (1970Brune ( , 1971)), can be repre-sented as M 0 /m 0 indicates the seismic moment ratio between the mainshock and EGF event at the lowest frequency, and f cm and f ca are the corner frequencies of the mainshock and the EGF event, respectively.Using the formulas from Irikura (1986) and Miyake et al. (2003), we can derive the constant levels of the acceleration and displacement spectra of the mainshock and EGF event and then estimate the scaling parameters N and C, which are the ratios of the fault dimensions and stress drop, respectively, between the mainshock and the EGF event.Here, U 0 /u 0 indicates the ratio between the mainshock and EGF event for the flat levels of the displacement spectra.The above parameters can be estimated through the weighted least-squares method of Miyake et al. (1999Miyake et al. ( , 2003)).We selected four stations (open triangles in Fig. 1) to obtain the source spectral ratios of the mainshock to the EGF event.The frequency range of strong  motion record is 0.5 -10 Hz, and the three component amplitude spectra vector summation was calculated for a 30-s waveform duration to include the entire S-wave.Figure 3 shows the mean observed source spectral ratio of the mainshock to the EGF event and the parameters determined from the source spectral ratio fitting.
The SMGA, which can be defined as the characteristic area with a uniform high slip velocity in the total rupture area and divided into N × N subfaults with sizes equivalent to the EGF event rupture area, reproduces the near-source strong ground motions in the frequency range up to approximately 10 Hz (Irikura 1986;Irikura and Kamae 1994;Miyake et al. 2001).The parameters related to SMGA were obtained after minimizing the residuals between the observed and synthetics for the displacements and the envelopes of the acceleration (Miyake et al. 1999), including the initiation position (rupture starting point), the size of the SMGA, the rupture velocity (V r ), and the rise time ( x), as listed in Table 3. Comparisons of the observed and synthetic waveforms at four stations using source spectral ratio fitting are shown in Fig. 4. Inaccuracy in both the focal mechanism and attenuation effect might introduce errors into the synthetic waveforms, and differences in the radiation patterns between the mainshock and EGF are also an important factor.Nevertheless, most features of the observed velocity and displacement records could be well reproduced.

Hybrid Simulation of broadband Ground Motion
We obtained broadband synthetic ground motions through a pair of matched filters and linearly summing lowfrequency (0.08 -0.5 Hz) synthetics from the F-K method and high-frequency (0.5 -10 Hz) synthetics from the EGF method.Figure 5 shows a comparison of the synthetic and observed broadband ground motion waveforms, including the stations not used in the SMGA modeling (solid triangles in Fig. 1).Due to limited information on velocity structure, uncertainty in the slip model, bias in the focal mechanism between the mainshock and the EGF event, and so on (as discussed in previous two sections), we cannot expect a perfect match for all waveform details between the observation and simulation.Except for some northern stations (e.g., HWA001, HWA032, and HWA033), the agreement is satisfactory for most of the stations.The results revealed here preliminarily manifest the ability of this hybrid scheme in broadband ground motion simulation.

Application of broadband Waveforms
We simulated both high-and low-frequency ground motion waveforms in this study.In the high-frequency part, the EGF method is a useful technique for retrieving the primary rupture characteristics of an earthquake by introducing minimal a priori assumptions for the available seismic data, and Fig. 4 indeed shows that the synthetics can reasonably explain the observations.As long as a suitable small event is available, by applying the EGF method, we can retrieve    important source information on the main causes affecting the strong ground motion behavior.
Most previous studies adopted 3D simulation approaches to estimate long-period ground motions (e.g., Kamae et al. 1998;Furumura et al. 2002;Graves and Pitarka 2010).In addition, Lee et al. (2014) showed that their model explained long-period ground motions well, so we also tried to estimate the long-period ground motions following the source rupture model derived by Lee et al. (2014).However, the 3D synthetic ground motions also revealed the problem that without the near-surface velocity structure, the amplitudes of synthetics were also insufficient to fit the observations, as discussed by Kamae et al. (1998).The model of Huang et al. (2014) indicates that the velocity structure beneath those stations for our study has a relatively small 3D effect.The effects of velocity structure and site might be corrected for each station to obtain better waveform fitting, however, our goal is to reproduce the main characteristics of the broadband strong motion but not to match all details of the data records.In addition, the 1D F-K method with nearsurface velocity structure requires less computation time and fits our requirements in this study.Lee et al. (2014) indicated that the 2013 Ruisui earthquake displayed rupture directivity effects toward the northeast.In Fig. 5, we noticed that some northern stations (e.g., HWA001, HWA032, and HWA033) display clearly smaller amplitudes of later pluses for synthetics.The recent 2013 Nantou earthquakes were both thrust events showing different rupture directivity effects, which could be observed through the differences in the source spectral characteristics between stations located in the directions of the forward and backward ruptures (Wen et al. 2014b).Figure 3 shows the source spectral ratios of the 2013 Ruisui mainshock to the EGF event.Although station HWA032 exhibits a much higher low-frequency spectral ratio, it is not easy to identify that the southern stations (e.g., HWA021 or HWA037) were located in the direction of backward rupture, which should show lower corner frequency and gentle spectral decay (Boore and Joyner 1978;Miyake et al. 2001).In addition, the rupture starting point in the SMGA (Fig. 1) did not suggest clear rupture directivity propagating to the north.Aagaard et al. (2004) mentioned that the propagation direction perpendicular to the slip direction during the faulting would show the minimal directivity, and this is closer to the rupture pattern of the 2013 Ruisui event.Nevertheless, the slip model of Lee et al. (2014) showed that there was a small asperity in the northern part, which was also located north of station HWA032.This small asperity might be the factor causing localized high PGA for these more northern stations.The low-frequency simulations (which considered the full slip model) better explain the observations for the northern stations.As mentioned by Kamae et al. (1998), a more detailed 3D velocity structure, including the nearsurface part beneath each station, would better improve the waveform fitting.On the other hand, because the spectral ratio fitting was applied to the averaged spectral ratio and because we considered the main SMGA, this process might lead to underestimation of the effect of the other asperity.Most earthquakes occurred in Taiwan are moderate-sized events, so we usually could assume one asperity to represent the main faulting feature (Yen and Ma 2011).Since we used the 2013 Ruisui earthquake as a case study to investigate the earthquake source characteristics using broadband strong motion simulation through the hybrid scheme, analysis with multiple SMGAs could be considered for more comprehensive investigations in the future.

Source Scaling
Figure 6 shows the source scaling characteristic analysis result, and both the SMGA and rise time are in good  agreement with the scale.Here, we also estimated the stress drop for the SMGA and obtained a value of 10.1 MPa, which is consistent with most earthquakes occurred in Taiwan (Yen and Ma 2011).Wen et al. (2014b) suggests that the high stress drop of the two 2013 Nantou earthquakes might be related to immature blind fault.Chuang et al. (2014) and Lee et al. (2014) both concluded that the most rupture of 2013 Ruisui earthquake is below 5 km, and this might be the reason for the high stress drop for the 2013 Ruisui event.Therefore, it is important to put more attention to the events originating from the buried faults for the mitigation of any future seismic hazards.

concLuSIon
We carried out a hybrid ground motion simulation to investigate the source and strong ground motion characteristics for the 2013 Ruisui earthquake.The results suggest that for the long-period strong motion simulation, the near-surface velocity structure is an important factor for controlling the waveform amplitude.However, a smaller asperity could also cause localized high PGA for the stations more distant from the main asperity.This case study helps us to further understand earthquake source properties through the different viewpoints of various frequency bands, and it is important and integral for future seismic hazard assessment.

Fig. 1 .
Fig. 1.Epicenters (black stars) and focal mechanisms of the 2013 Ruisui earthquake and the EGF event listed in Table 1.Gray star shows the location of 2014 M w 5.6 Fanglin earthquake.Triangles indicate free-field strong motion stations.Open triangles show stations used in the source spectral ratio analysis and the EGF simulation.The active faults (thick gray lines) identified by the Central Geological Survey of Taiwan are also shown.The upper inset shows the rupture starting point (star) on the strong motion generation area (SMGA).(Color online only)

Fig. 2 .
Fig. 2. Comparison of low-frequency observed (black lines) and synthetic (gray lines) waveforms."obs" shows the observed record, "1D-syn" shows the 1D simulation without near-surface velocity structure, and "1DS-syn" shows 1D simulation with near-surface velocity structure.The numbers above the waveforms indicate the maximum amplitudes of the observed records, in cm s -2 .

Fig. 3 .
Fig. 3. Source spectral ratios of the mainshock to the EGF event and average observed source spectral ratio (thick gray line).The parameter values determined from the source spectral ratio fitting are listed.

Fig. 4 .
Fig. 4. Comparison of observed (black lines) and synthetic (gray lines) waveforms at strong motion stations used for source modeling in the EGF method, with the upper number indicating the maximum amplitudes of the observed records and the lower number indicating the cross-correlation coefficient of the synthetic and observed data.

Fig. 6 .
Fig. 6.The scaling relationship for seismic moment versus (a) the SMGA and (b) the rise time from Miyake et al. (2003).Squares indicate the two 2013 Nantou earthquakes estimated by Wen et al. (2014b).Inverted triangle shows our result for the 2013 Ruisui earthquake superimposed on the figures.

Table 1 .
Earthquake Parameters for the 2013 Ruisui Earthquake and the EGF event.The epicenters and fault plane solutions were determined by CWB and BATS, respectively.