Scenario for a Short-Term Probabilistic Seismic Hazard Assessment ( PSHA ) in Chiayi , Taiwan

Using seismic activity and the Meishan earthquake sequence that occurred from 1904 to 1906, a scenario for short-term probabilistic seismic hazards in the Chiayi region of Taiwan is assessed. The long-term earthquake occurrence rate in Taiwan was evaluated using a smoothing kernel. The highest seismicity rate was calculated around the Chiayi region. To consider earthquake interactions, the rate-and-state friction model was introduced to estimate the seismicity rate evolution due to the Coulomb stress change. As imparted by the 1904 Touliu earthquake, stress changes near the 1906 Meishan and Yangshuigang epicenters was higher than the magnitude of tidal triggering. With regard to the impact of the Meishan earthquake, the region close to the Yangshuigang earthquake epicenter had a +0.75 bar stress increase. The results indicated significant interaction between the three damage events. Considering the path and site effect using ground motion prediction equations, a probabilistic seismic hazard in the form of a hazard evolution and a hazard map was assessed. A significant elevation in hazards following the three earthquakes in the sequence was determined. The results illustrate a possible scenario for seismic hazards in the Chiayi region which may take place repeatly in the future. Such scenario provides essential information on earthquake preparation, devastation estimations, emergency sheltering, utility restoration, and structure reconstruction.


InTroduCTIon
Due to the interaction between the Eurasian and Philippine Sea Plates, many devastating earthquakes have taken place in Taiwan (Wang and Shin 1998;Hsu 2003;Wu et al. 2008;Cheng et al. 2012).The Chiayi region, located in southwestern Taiwan with a high deformation rate (Yu et al. 1997), has, in the past, been struck by several destructive earthquakes.In 1792, a devastating earthquake struck the Chiayi region which resulted in more than seven hundred fatalities and eight thousand collapsed buildings (Cheng et al. 2012).This earthquake can be attributed to the rupture of the Meishan Fault (Fig. 1).Additionally, the Meishan earthquake sequence (Fig. 1) which took place between 1904 and 1906 caused thousands of fatalities and thousands of building collapses (Hsu 2003).The earthquake sequence included three major events, the November 5 th , 1904 M w 6.1 Touliu earthquake, the March 16 th , 1906 M w 6.9 Meishan earthquake (the largest event), and the April 13 th , 1906 M w 6.4 Yanshuigang earthquake (Prof.Shih-Nan Cheng, personal communication).Since the 1906 Meishan earthquake, the population has greatly increased and no M w > 6.5 event has occurred.Therefore, it is necessary to understand seismic hazards in this region.
For the purpose outlined above, it is important to apply a PSHA in terms of the probability of strong ground motion for Taiwan, especially within the Chiayi region.In 1999, the Global Seismic Hazard Assessment Program (GSHAP) published a global probabilistic seismic hazard map that included the Taiwan region (Shedlock et al. 2000).However, the GSHAP considered earthquake catalogs from global seismic networks rather than a detailed seismicity catalog from the Taiwan regional seismic network.Also, a detailed assessment of special areas, such as those in the Chiayi region, was not provided.Another hazards assessment for Taiwan was proposed by Campbell et al. (2002) who considered seismic catalogs and active fault parameters for the Taiwan region.However, the ground motion prediction equations (GMPEs) they considered were obtained from the world rather than from Taiwan.Cheng et al. (2007) evaluated seismic hazards for Taiwan and proposed a hazard map by incorporating a catalog from a local network, active fault parameters, seismogenic zones, and local GMPEs in Taiwan.Such studies are useful for understanding seismic hazards in Taiwan.However, many parameters and the database for seismic hazard assessments have been revised and/or updated.For example, Lin (2009) and Lin and Lee (2008) obtained GMPEs for crust and subduction events in the Taiwan region, respectively.Hopefully, these GMPEs can better represent ground motion attenuation behaviors in Taiwan.As a result, it is necessary to use these new datasets for an analysis.
The basic concept of the traditional PSHA has recently been called into debate.Based on traditional PSHAs (e.g., Cornell 1968;McGuire 1976), seismic hazards can be calculated using the assumption that earthquakes occur with a known time-independent rate.However, the case of the 2010 Darfield, New Zealand earthquake sequence has shed light on the limitations of this assumption (Chan et al. 2012a).The M w 6.3, February 21 st , 2011 Christchurch, New Zealand earthquake is regarded as an aftershock of the M w 7.1, September 4 th , 2010 Darfield mainshock.However, the Christchurch earthquake caused severe damage in downtown Christchurch.Obviously, the Darfield earthquake raised the seismic hazard probability for Christchurch.Earthquakes that contribute to seismic hazards may not be totally independent of each other.For such cases, earthquake interactions and time-dependent seismic hazards should be analyzed.The Meishan earthquake sequence could also increase the awareness of the unsuitability of traditional PSHAs for the Chiayi region.According to the time and space relationship, the 1904 Touliu and 1906 Yanshuigang earthquakes can be regarded as a foreshock and an aftershock of the 1906 Meishan earthquake, respectively.However, all three earthquakes caused severe disasters (Hsu 2003).Such a circumstance points to the importance of seismic hazards imparted by an earthquake sequence.Therefore, a near real-time updating or a time-dependent seismic hazards evaluation is necessary.
In this study, we applied an assessment for short-term probabilistic seismic hazards (Fig. 2) to a scenario in the Chiayi region that included the Meishan earthquake sequence.We first evaluated the long-term seismicity rate based on a smoothing kernel (Woo 1996) in surroundings wherein earthquakes were independent of the Meishan earthquake sequence.When considering earthquake interactions and seismic characteristics on physical terms, the rateand-state friction model (Dieterich 1994) was introduced for evaluating the seismicity rate change due to the Coulomb stress change (Harris 1998;Cocco and Rice 2002) imparted by the earthquake sequence.Considering the source, path, and site effect using GMPEs, probabilistic seismic hazards in terms of strong ground motion could be assessed.In this work, we discuss the spatial and temporal distribution of seismic hazards during the Meishan earthquake sequence.

MeTHodS And ProCedureS
The long-term seismicity density rate, the short-term seismicity rate evolution, and the GMPE are important factors for short-term PSHAs.In this study, the long-term seismicity density rate was estimated according to a smoothing kernel function; the short-term rate change was calculated using the rate-and-state friction model; and, the probabilistic seismic hazard was assessed by introducing the source, the path effect, and the site amplification with GMPEs (Fig. 2).The procedures of the assessment are provided in subsequent sections.

The Long-Term Seismicity rate
Based on the study by Woo (1996), the long-term seis-micity rate λ(M, x) at the target site (x) can be described as a function of the magnitude (M) as follows: where K(M, xx i ) is the kernel function as a function of the epicentral distance between the target site (x) and the epicenter of the i'th earthquake (x i ); T M is the complete time period for M in a catalog; N M is the number of earthquakes with a magnitude of M in a catalog.In order to acquire the long-term annual seismicity rate, the kernel function was summed over all of the events in a catalog and was divided by the complete period of the catalog.In other words, the behavior of seismic activity was assumed to be temporally stationary when rate was not perturbed by an earthquake.
According to the definition of Woo (1996), the kernel function K(M, xx i ) can be described as follows: where PL is the power law index.The bandwidth function, H(M), is defined as the distance between two events, as an exponential function of the magnitude (M), and is represented as follows: (3) where c and d are constants that can be obtained by regression in an earthquake catalog.Note that variations in the seismic densities for different magnitudes can be illustrated by the magnitude-dependent bandwidth function.

The Short-Term Seismicity rate evolution
In this study, the rate-and-state friction model (Dieterich 1994) was considered in order to evaluate short-term seismicity rate evolution.For the application of this model, the Coulomb stress change (ΔCFS) by earthquakes was calculated.According to the constant apparent friction law (Harris 1998;Cocco and Rice 2002), the general expression for the Coulomb stress change, ΔCFS, can be represented as follows: where Δτ is the shear stress change computed along the slip direction on the receiver fault; nl is the apparent friction coefficient; and Δσ n is the normal stress change perpendicular to the receiver fault.Previous studies (e.g., Toda et al. 2011;Chan et al. 2012b;Sevilgen et al. 2012) have found that a positive stress change enhances subsequent events.A negative stress change, in contrast, inhibits the occurrence of future seismicity.
One key piece of information for the ΔCFS calculation is the mechanism of the receiver fault.In this study, the focal mechanisms determined by Wu et al. (2010) were used as a reference for receiver faults.We assumed that a temporal-stationary and spatial-variable receiver fault for each calculation grid consisted of the reference focal mechanism with the shortest epicentral distance, i.e., a temporally immovable fault orientation is assumed.The procedure is in accordance with those of previous studies (Hainzl et al. 2010;Toda and Enescu 2011;Chan et al. 2012a, b).
We estimated ΔCFS on each 0.01° × 0.01° grid by solving spatial variable receiver faults using the COU-LOMB 3.3 program (Toda and Stein 2002).To minimize the depth uncertainty for the ΔCFS calculation, we followed the procedure proposed by Catalli and Chan (2012) and reported the maximum ΔCFS among seismogenic depths for each calculation grid.

Probabilistic Seismic Hazard Assessment
Proper GMPEs are important for PSHA applications.For crustal events, we used the GMPEs obtained by Lin (2009).These GMPEs are determined using the regression of 5968 seismic records from 60 earthquakes recorded by the seismic network of the Taiwan Strong Motion Instrumentation Program (TSMIP).The results indicated magnitude scaling, distance scaling, and site amplification effects for ground-motion data in Taiwan.Additionally, the GMPEs considered the focal mechanism and determined that thrust events have the largest ground-motion, whereas normal mechanism events display the smallest ground-motion under the same conditions (i.e., magnitude, distance, and site effects).Additionally, the GMPEs provided the corresponding standard deviations.Understanding the variability within these equations is useful for PSHA.
Previous studies (Crouse et al. 1988, and references therein) have suggested that ground-motion attenuation behaviors for subduction earthquakes are distinct from those for crustal earthquakes.In northeastern Taiwan, the Philippine Sea Plate subducts to the north beneath the Eurasian Plate; in southern Taiwan, the Eurasian Plate subducts to the east beneath the Philippine Sea Plate (Fig. 3).Therefore, it is necessary to consider specific GMPEs for subduction earthquakes (i.e., for both interface and intraslab events).Lin and Lee (2008) proposed the first GMPEs for Taiwan subduction events and analyzed strong-motion records for subduction earthquakes obtained using a TSMIP network and the Strong Motion Array Taiwan Phase I (SMARTI).The GMPESs proposed by Lin and Lee (2008) were used in this study since they are the only GMPEs for subduction events in Taiwan.In order to consider site amplification, the averaged shear-velocity down to 30 m (known as Vs30), obtained by Lee and Tsai (2008), was utilized.

The Long-Term Seismicity rate
Employing the smoothing kernel function, we evaluated the long-term seismicity rate by analyzing an earthquake catalog.Since the aftershock pattern was evaluated using the rate-and-state friction model, the catalog for the evaluation of the long-term seismicity rate had to be declustered.In other words, foreshocks and aftershocks had to be removed from the catalog.In this study, the approach proposed by Burkhard and Grünthal (2009) was applied for the de-clustering process.The approach was originally developed in central Europe; however, it also performs well for the Taiwan region (Chan et al. 2012b), as well as for subduction environments (Suckale and Grünthal 2009).The seismic catalog used in this study is obtained from Cheng et al. (2010).In order to improve reliability, we considered the catalog for a period when the seismic network recorded all earthquakes within a certain magnitude threshold, known as the magnitude of completeness, M c .This catalog obtained a better record quality since 1935 due to improvement of a local network that spans the entire island of Taiwan.The maximum curvature approach (Wiemer and Wyss 2000) was employed to check the temporal evolution of the M c , and it was determined that the M c was smaller than 5.0 after 1940.Therefore, earthquakes with M ≥ 5.0 between 1940 and 2010 (Fig. 3) were considered in order to establish the long-term rate model.Note that the period for the catalog we used is after the occurrence of the Meishan sequence.In other words, we assumed a temporally immovable longterm seismicity rate when stress does not perturbed by large earthquakes.
To make use of the smoothing kernel function [Eq.( 2)], the bandwidth function in the form of c and d values [Eq.( 3)] for the catalog is needed.Using a linear regression for the catalog, we obtained c and d values that were 1.643 and 0.451, respectively.
Recommended values for the power law index, PL, in [Eq.( 2)] are between 1.5 and 2.0.The boundaries correspond to cubic or quadratic decays with a hypocentral distance for seismic activity (Molina et al. 2001).Chan et al.
(2010) found insignificant differences between the results when the power law index was assumed to be between the recommended values.Therefore, we assumed an intermediate recommended value of 1.75 for calculation.
Based on the parameters of the smoothing kernel function acquired above, the long-term seismicity rate was obtained (Fig. 4).The highest seismicity rate was determined in southwestern Taiwan near the Chiayi region corresponding to frequent seismicity (Fig. 3).Other regions with high seismicity rates were found within the southeastern and eastern offshore regions of Taiwan.They corresponded to the boundary between the Eurasian and Philippine Sea Plates.The active seismicity can be associated with a high deformation rate observed by GPS stations (Yu et al. 1997).
Using the kernel function, long-term seismicity rate could be evaluated.Based on the characteristics of this approach, the spatial distribution of the rate can be presented  With a verified applicability for seismicity forecasts, the approach will be useful for subsequent seismic hazard assessments.In addition, their result confirmed a temporally immovable longterm seismicity rate in Taiwan.

The Short-Term Seismicity rate Change
The ΔCFS calculation requires knowledge of the rupture parameters for source events, such as the dimension of the rupturing fault and the slip magnitude.We obtained these parameters from the scaling laws proposed by Yen and Ma (2011), as follows: where L is the rupture length in kilometers, M 0 is the seismic moment; W is the rupture width in kilometers; and AD is the average slip in centimeters.Three of the events with M w ≥ 6.0 in the Meishan earthquake sequence (i.e., the 1904Touliu 1906Meishan, and 1906 Yangshuigang earthquakes) (Fig. 1) were regarded as source events for ΔCFS calculations.The corresponding source parameters (i.e., occurrence time, hypocenter locations, moment magnitudes, and focal mechanisms) and the rupture parameters (i.e., rupture length, rupture width, and averaged displacements) are presented in Table 1.
For the ΔCFS calculation, an intermediate value of nl = 0.4 was considered.The value corresponds to a reasonable range of nl between 0.2 and 0.5, as inferred from the characteristics of earthquake focal mechanisms in Taiwan (Hsu et al. 2010).
The ΔCFSs imparted by the Meishan earthquake sequence were evaluated in the Chiayi region (Fig. 5).The size of the study region is determined by the influence of the Meishan earthquake sequence since the ΔCFSs in the regions far away from the sequence were not significantly perturbed.The results indicated remarkable increases in the vicinity of each event.The M w 6.9 Meishan earthquake (Fig. 5b) caused a stress increase with a wider range and a higher magnitude.In contrast, the M w 6.1 Touliu earthquake (Fig. 5a) caused less significant stress perturbations.Such a discrepancy can be attributed to the magnitude difference between earthquakes.As imparted by the Touliu earthquake, stress changes close to the Meishan and Yangshuigang epicenters were +0.07 and +0.06 bars, respectively (Fig. 5a).The values are higher than the magnitude of tidal triggering (+0.01 bars) (Stein 2004).Considering the stress change of the Meishan earthquake, the epicenter of the Yangshuigang earthquake had a +0.75 bar stress increase.The result confirms the mainshock-aftershock relationship between these two events.According to the results, a ΔCFS ≥ 0.01 bar can trigger consequent earthquakes and can be a threshold for earthquake warnings, as well as the monitoring of the Coulomb stress evolution.
For modeling the spatial and temporal distribution of aftershock sequences or triggered earthquakes, ΔCFS incorporated into the rate-and-state friction model (Dieterich 1994) was proposed in order to evaluate the short-term seismicity rate change.For application of the rate-and-state friction model, previous studies (e.g., Toda and Stein 2003;Toda et al. 2005;and Catalli et al. 2008) suggested that the physically reasonable range for Aσ was between 0.1 and 0.4 bars.In this study, a fixed Aσ of 0.2 bars was assumed.The aftershock duration, t a , was assumed to be a function of the moment magnitude, M w , as proposed by Burkhard and Grünthal (2009)  According to the corresponding magnitudes, the t a for the 1904 M w 6.1 Touliu, the 1906 M w 6.9 Meishan, and the 1906 M w 6.4 Yanshuigang earthquakes were 256, 490, and 328 days, respectively.
Seismicity rate changes within the Chiayi region at various time spots were calculated using the rate-and-state friction model (Fig. 6).Remarkable seismicity rate increases were determined in the region close to the epicenter of the events when they just occurred.Rate decays with time are illustrated with the rate changes at various time spots (i.e., from Figs. 6a to b; from Figs. 6c to d; and from Figs. 6e to f).Prior to the occurrence of the 1906 Meishan earthquake (Fig. 6b), the change of the seismicity rate at the epicenter was insignificant.The result can be attributed to the relatively small magnitude of the 1904 Touliu earthquake and long time interval between the two earthquakes.Following the occurrence of the Meishan earthquake (Figs.6c and d), a significant seismicity rate increase was determined in the Chiayi region.The results indicated that a large earthquake near the epicentral region of the Yanshuigang earthquake could be expected.It is worth mentioning that some regions with the highest Coulomb stress increases did not occur consequent to damaging earthquakes (Fig. 5).Such results can be ascribed to the low level of cumulative stress, far from next rupture, or to corresponding tectonics not detected with seismogenics.
Table 1.Rupture parameters of the source events for calculating the Coulomb stress change.The source parameters (i.e., occurrence time, hypocenter locations, moment magnitudes, and focal mechanisms) were obtained by Cheng et al. (personal communication).The rupture parameters (i.e., rupture length, rupture width, and averaged displacement) were determined by scaling laws (Yen and Ma 2011)

The Temporal evolution of Seismic Hazards
The temporal evolution of seismic hazards for some cities within the Chiayi region from 1904 until 1908 was evaluated (Fig. 7).Prior to the 1904 Touliu earthquake, the annual exceedance probabilities for peak ground acceleration (PGA) = 0.6 g were low (between 1.9 and 2.6‰).Following the occurrence of each earthquake, a significant increase in the probability was obtained.Note that the highest hazards were expected for the city with the shortest epicenter distance (Fig. 1).Following the 1904 Touliu earthquake, the highest hazard was determined in Xingang City; following the 1906 Meishan earthquake, the highest hazard was determined in Meishan City; and following the 1906 Yanshuigang earthquake, the highest hazard was determined in Chiayi City.The results might be associated with the distribution of the damage and fatalities for each event (Prof.Shih-Nan Cheng, personal communication).At the end of 1908, the probability for strong ground motion was expected to be similar to background levels.

Seismic Hazard Maps
The temporal evolution of seismic hazards in the three cities surrounding the Chiayi region is presented (Fig. 7).The results indicated a significant higher probability for strong ground motion following the three earthquakes of the Meishan sequence.In this section, the spatial distributions of seismic hazards are evaluated further.Seismic hazard maps at four time spots were obtained (Fig. 8).A relatively low hazard was determined along the western coast region prior to earthquake perturbation (Fig. 8a).On the other hand, hazards were higher in the east within the study region.The pattern corresponded to the results of previous studies (Campbell et al. 2002;Cheng et al. 2007).Following the occurrence of each earthquake (Figs.7b to d), significantly higher hazards were evaluated in the vicinity of the epicenter.A lower hazard was expected following the 1904 M w 6.1 Touliu earthquake than after the 1906 M w 6.9 Meishan earthquake.The discrepancy could be attributed to the magnitude difference of the ΔCFS imparted by the two earthquakes (Figs.5a and b).Such results could provide essential information for emergency response for victim sheltering and relocation immediately following a devastating earthquake.

ConCLuSIonS
In this study, we assessed short-term probabilistic seismic hazards during the period of the Meishan earthquake sequence.Our results indicated the migration of seismic hazards during the sequence.According to the Coulomb stress changes imparted by the three events of the Meishan earthquake sequence (Fig. 5), stress interactions confirmed inter-dependency.The seismicity rate evolution was illustrated using the rate-and-state friction model (Fig. 6), and displayed a remarkable rate increase in the region close to the epicenters immediately following each earthquake.Short-term seismic hazards can be estimated by considering seismicity rate evolution and ground motion prediction Fig. 7.The temporal evolution of the annual exceedance probability in three study cities (the squares in Fig. 1) for PGA = 0.6 g.Variations in the seismic hazards are attributed to the seismicity rate change imparted by the three study events (Fig. 6).
equations.The results suggested that a significant increase in hazards followed the three earthquakes of the sequence, especially in cities adjacent to earthquakes (Fig. 7).
According to the results of the seismic hazard evolution (Fig. 7), the seismic hazard in cities of interest can be forecast in near real time mode.Spatial distributions of the hazards are represented according to seismic hazard maps (Fig. 8) and could be of benefit not only to decision-makers but also to public officials for seismic hazard mitigation.For example, results from the scenario illustrate the possible evolution of seismic activity and seismic hazards during an earthquake sequence within the Chiayi region.They can provide essential information for the long-term preparation for seismic hazard responses prior to the occurrence of a catastrophic earthquake.Additionally, the rapid re-evaluation for short-term seismic hazards immediately following devastating earthquakes could be beneficial for devastation estimations, emergency sheltering, utility restoration, and/ or structural reconstruction subsequent to a devastating earthquake.

Fig. 1 .
Fig. 1.Distribution of three earthquakes with M w ≥ 6.0 for the Meishan earthquake sequence during 1904 and 1906.Stars indicate the epicenters of the earthquakes.Active faults are shown as red lines.Three cities of interest (Xingang, Chiayi, and Meishan) used for seismic hazards evaluations (Fig. 7) are shown as solid squares.The location of the study region is shown in Fig. 3.

Fig. 2 .
Fig. 2. The flow chart for the procedures of the approach used for this study.

Fig. 3 .
Fig. 3.The distribution of de-clustered earthquakes with M ≥ 5.0 from 1940 to 2010.Crustal, interface, and intraslab earthquakes are shown as white, gray, and black circles, respectively.The three study events for the Meishan earthquake sequence are shown as stars.The solid rectangle represents the study region.Dashed polygons illustrate the shapes of the subduction zones.The inset represents the tectonic setting surrounding Taiwan.Note that some interface events in northeastern Taiwan take place east of the longitude of 122° are not shown in the figure.
as a function of the magnitude.Chan et al. (2012b) have proved the forecasting ability of this approach with an application to the Taiwan region by constructing a forecasting model by considering the earthquake catalog from 1973 -2007.To compare seismicity during 2008 -2009, they found a positive correlation between the two.

Fig. 4 .
Fig.4.The estimated long-term seismicity rate for M ≥ 5.0.Three study events with M w ≥ 6.0 are shown as white stars.The solid rectangle represents the study region.Note that the magnitude-dependent seismicity density was obtained.However, for simplicity the seismicity rate for the M ≥ 5.0 threshold is presented instead of that for each magnitude bin.
is represented as follows:

Fig. 8 .
Fig. 8.The seismic hazard maps (a) before Touliu, (b) immediately following Touliu, (c) immediately following Meishan, and (d) immediately following Yanshuigang.The stars represent the epicenters of the three earthquakes.The solid squares represent the cities of interest, used for seismic hazards evaluations in Fig. 7.