Time-dependent probabilistic seismic hazard assessment and its application to Hualien City , Taiwan

Here, we propose a time-dependent probabilistic seismic hazard assessment and apply it to Hualien City, Taiwan. A declustering catalog from 1940 to 2005 was used to build up a long-term seismicity rate model using a smoothing Kernel function. We also evaluated short-term seismicity rate perturbations according to the rate-and-state friction model, and the Coulomb stress changes imparted by earthquakes from 2006 to 2010. We assessed both long-term and short-term probabilistic seismic hazards by considering ground motion prediction equations for crustal and subduction earthquakes. The long-term seismic hazard in Hualien City gave a PGA (peak ground acceleration) of 0.46 g for the 2.1 ‰ annual exceedance probability. The result is similar to the levels determined in previous studies. Seismic hazards were significantly elevated following the 2007 ML = 5.8 earthquake that occurred approximately 10 km from Hualien City. This work presents an assessment of a suitable mechanism for time-dependent probabilistic seismic hazard determinations using an updated earthquake catalog. Using minor model assumptions, our approach provides a suitable basis for rapid re-evaluations and will benefit decision-makers and public officials regarding seismic hazard mitigation.


Introduction
Human life and property are threatened by devastating earthquakes.As a result, seismic hazard mitigation is an important issue for seismologists, earthquake engineers, and related scientists.Probabilistic seismic hazard assessments (PSHAs), in terms of the probability of strong ground motion, provide a key reference for the determination of seismic mitigation policies (e.g., building codes and the site selection of pubic structures) (McGuire, 2001).Therefore, scientists in different fields are currently attempting to build reliable systems for PSHAs on different spatial scales.For example, on the national scale, the Headquarters for Earthquake Research Promotion (http://www.jishin.go.jp/main/ index-e.html)proposed a seismic hazard map for Japan.On the regional scale, the subproject JRA2 in the project Network of Research Infrastructures for European Seismology (NERIES, http://www.neries-eu.org/),constructed a probabilistic seismic hazard map for the Euro-Med region.On the worldwide scale, the Global Earthquake Model (GEM, http://www.globalquakemodel.org/)integrated several research institutes, insurance, and re-insurance enterprises in order to evaluate seismic hazards around the world.
Taiwan is located in the plate boundary between the Eurasian and Philippine Sea plates and has a high earthquake activity (Fig. 1).Some of the large earthquakes occurring in this region may lead to a loss of property and human life.Therefore, it is essential to find a means of seismic hazard mitigation.Building a seismic hazard assessment system for Taiwan is one practical approach.
Seismic hazards in Taiwan have been evaluated by many groups.The Global Seismic Hazard Assessment Program (GSHAP, http://www.seismo.ethz.ch/static/GSHAP/)obtained a global probabilistic seismic hazard map that included the Taiwan region.However, they used earthquake catalogs from worldwide seismic networks rather than a detailed seismicity catalog from Taiwan.Cheng (2002) and Cheng et al. (2007) Learning period: 1940-2005Forecasting period: 2006-2010 Interface eq Intraslab eq.
In te rfa ce re gio n In proposed a hazard map by incorporating a catalog from a local network, active fault parameters, and seismogenic zones in Taiwan.In addition, they deaggregated the seismic hazards contributed from different earthquake magnitudes and source distance ranges.Another Taiwan seismic hazard map was proposed by Campbell et al. (2002) using seismic catalogs, active fault parameters, and ground motion prediction equations (GMPEs) for the world, the United States, and Taiwan.Such studies are useful for understanding seismic hazards in Taiwan.However, after the results were published many parameters and the database for seismic hazard assessments, such as distribution of active faults, GMPEs, and earthquake catalogs, were revised and/or updated.Therefore, the evaluation of seismic hazards could be more precise by considering timely information.
The concept of the traditional PSHA has recently been called into question.One of the standard procedures for traditional PSHAs (Cornell, 1968;McGuire, 1976) is the construction of return periods for characteristic earthquakes that are independent of one another.However, it is well recognized that earthquake activity varies with time through periods of seismic-quiescence and -bursts.In other words, earthquakes do not have regularity of occurrence time.Thus, several studies, in respect to physics and statistics, have overruled this hypothesis.Kagan and Knopoff (1978) proposed the time-space Epidemic Type AfterShock model (ETAS).In the ETAS model, each earthquake is regarded as a mainshock that may trigger consequent earthquakes.A physics-based model can explain similar phenomenon.King et al. (1994) concluded that the Coulomb stress change imparted by earthquakes might control the distribution of future earthquakes.Consequent earthquakes generally occur in regions where Coulomb stresses are increased.On the contrary, only a few earthquakes occur in Coulomb stress decreased regions.The case of the 2010 Darfield, New Zealand, earthquake sequence determined disadvantages for the traditional PSHA (Chan et al., 2012a).The M w = 6.3,21 February 2011 Christchurch earthquake is regarded as an aftershock of the M w = 7.1, 4 September 2010 Darfield earthquake.However, the Christchurch earthquake caused severe damage in downtown Christchurch.Besides, not only aftershocks may result in consequent seismic hazards, but also next larger earthquakes can further expand hazards.For example, a M w = 7.4 earthquake took place off the Pacific coast of Tohoku, Japan, on 9 March 2011 (Nettles et al., 2011).Since the epicenter is away from land, the resulting damages are negligible.Fifty-one hours after the earthquake on 11 March, an earthquake with M w = 9.1 took place and resulted in disasters in Japan.Such circumstances point out the importance of an aftershock sequence in seismic hazard evaluations and suggest the re-evaluation of seismic hazards immediately following a large earthquake.
In this study, we propose an approach for the shortterm PSHA, which changes after the occurrence of a large earthquake.Based on a smoothing Kernel, we first evaluated the long-term earthquake occurrence rate.By considering earthquake interactions and seismic characteristics in physical terms, we introduced the rate-and-state friction model (Dieterich, 1994) in order to evaluate seismicity rate changes due to the Coulomb stress changes imparted by previous earthquakes.We applied this approach to Hualien City, which is located in a high seismicity region of the plate boundary between the Eurasian and Philippine Sea plates (Fig. 1).We evaluated the recent seismic hazard evolution of Hualien City and checked the feasibility of this approach.We also discussed the importance and applicability of the timedependent PSHA.

Procedures of the approach
When evaluating time-dependent PSHA, the long-term seismicity density rate, the short-term seismicity rate change, and the GMPE are important factors.According to the PSHA method proposed by Cornell (1968) and McGuire (1976), the total frequency ν T m , with which ground motion parameter A exceeds a specific value z during time period t (usually one year) at a specific site, is estimated as follows: where γ represents the number of earthquakes with magnitudes m greater than selected minimum magnitude M MIN and smaller than selected maximum magnitude M MAX , T m represents the time interval of the complete catalog for the magnitude m, P (A > z|m, x i − x) is the probability that given ground motion parameter z is exceeded a specific value z, conditional on an earthquake m and the epicentral distance from the target site x to the epicenter of the i-th earthquake x i , f M (m) and f R (x i − x|m) are the probability density functions for earthquake magnitude and distance which describe the relative likelihood of different earthquake scenarios, respectively.Below, we outline how we evaluated the long-term seismicity density rate model acquired using a smoothing Kernel function, the short-term seismicity rate change model using the rate-and-state friction model, and the seismic hazard assessment by introducing the path effect and site amplification through GMPEs.

The long-term seismicity rate
In Eq. ( 1), γ dxdm can be regarded as the seismicity rate.We followed Woo (1996) and estimated the long-term seismicity rate as follows: where K (m, x − x i ) represents the Kernel function as a function of the magnitude and the epicentral distance from the target site x to the epicenter of the i-th earthquake x i .In order to acquire the long-term annual seismicity rate, the Kernel function was summed over all of the events in the complete catalog and divided by the corresponding duration.We followed Woo (1996) and described the Kernel function K (m, x − x i ) as follows: where PL denotes 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 can be represented as follows: where c and d are constants.We obtained these constants through a regression of the earthquake catalog.Seismicity rate models based on the magnitude-dependent bandwidth function can illustrate variations in the seismic densities for different magnitudes.

The short-term seismicity rate change
In this study, using the rate-and-state friction model (Dieterich, 1994), we evaluated the short-term seismicity rate change.For the application of this model, we first calculated the Coulomb stress change ( CFS) caused by source events.Based on the constant apparent friction law (Harris, 1998;Cocco and Rice, 2002), the general expression of the CFS can be represented as follows: where τ is the shear stress change computed along the slip direction on the receiver fault, µ is the apparent friction coefficient, and σ n is the normal stress change perpendicular to the receiver fault (note that a positive σ n represents unclamping).
Coulomb model calculations require knowledge of the rupture parameters for source events such as the geometry of the rupturing fault and the slip magnitude.Yen and Ma (2011) determined those relationships using Taiwan data and they were used in this study, in the following equations:  where L is the rupture length in kilometers, M 0 is the seismic moment in Newton-meters, W is the rupture width in kilometers, and AD is the average slip in centimeters.Although the uncertainty terms are proposed in the scaling law, they are not considered in this study for simplicity of calculations.In order to evaluate the seismic moment according to the local magnitude for each earthquake, we followed Wu et al. (2005) and described the relationship between the moment magnitude M W and the local magnitude M L as follows: Another key piece of information for the CFS calculation is the mechanism of the receiver fault.In general, receiver faults can be represented in the following forms: (1) optimally oriented fault planes according to a combination of the regional stress field and the stress change caused by the source event (King et al., 1994); (2) the mechanisms of active faults (Toda et al., 1998); and (3) the fixed focal mechanisms of earthquakes in subregions (Chan and Stein, 2009).Hainzl et al. (2010) investigated the effect of more realistic fault systems in a CFS model and found that considering earthquake nucleation on multiple receiver fault orientations significantly changed the spatial stress change pattern and the percentage of triggered events.In this study, we used the focal mechanisms determined by Wu et al. (2010) as a reference for the receiver faults of the CFS calculation (Fig. 2a).We assumed that a spatial variable receiver fault for each calculation grid consisted of the reference focal mechanism with the shortest epicentral distance (i.e., focal depth is not considered) (Fig. 2b).In other words, we assumed a temporally immovable fault orientation for the study region.The procedure was undertaken in accordance with applications of previous studies (Hainzl et al., 2010;Chan et al., 2012c;Toda and Enescu, 2011).
We estimated CFS on a 0.2 • × 0.2 • grid using spatial variable receiver faults by applying the COULOMB 3.3 program (Toda and Stein, 2002).In order to minimize the depth uncertainty for calculations, we followed the assumptions proposed by Catalli and Chan (2012) and reported the maximum CFS within the seismogenic depth for each calculation grid.
To quantify the impact of CFS on the seismicity rate, we used the rate-and-state friction model (Dieterich, 1994).We followed Chan et al. (2010) and determined the evolution of the seismicity rate change, R (m, x, t), using CFS with the n-th source event, CFS n (x), at the site of interest x as a function of magnitude, m, and time, t, as follows: , where R n−1 (m, x) is the short-term seismicity rate change promptly before the occurrence of the n-th source event, Aσ is a constitutive parameter of the model as described by Dietrich (1994), t n is the occurrence time of the n-th source event, and t na is the aftershock duration.According to this model, the evolution of seismicity rate for different magnitudes is simply controlled by the long-term seismicity rate, i.e., the evolution of the seismicity rate change is magnitude-independent since the influences imparted by the source events ( CFS n (x)) and tectonic loading (Aσ and t na ) are magnitude-independent.Thus, it is difficult to identify the differences of rate changes between small and large events through this procedure.By combining the long-term seismicity rate (Eq.2) and evolution of the seismicity rate change using Eq. ( 10), the short-term seismicity rate at time t can be represented as follows: The relationship describes the short-term seismicity rate evolution by considering a series of source events and represents a generalization of the relationships.The period of the seismicity rate perturbation depends on the magnitude of CFS and aftershock duration t na .

The probabilistic seismic hazard assessment
In addition to a reliable seismicity rate model, another item that should be considered is proper GMPEs for a PSHA (f R (r|m) in Eq. ( 1).Since Lin (2009) pointed out that global GMPEs do not provide a representative description of the local observations in Taiwan, for this work we only considered GMPEs obtained from the regression of ground motion observations in Taiwan.
Over the past decade, many studies have proposed GM-PEs in order to illustrate attenuation behaviors for different types of earthquakes in Taiwan.For crustal events, Wu et al. (2001) selected 60 crustal events with M L ≥ 5.0 and obtained GMPEs in the form of peak ground acceleration (PGA) and peak ground velocity (PGV) as a function of magnitude, as well as the nearest distance to the rupture surface.Liu and Tsai (2005) selected 51 events with M w between 4.0 and 7.1 and established GMPEs in terms of PGA and PGV.In the studies mentioned above, however, attenuation behaviors for different structural periods were not discussed.For this purpose, we expect to consider the GMPEs in the form of an acceleration response spectrum (SA) for more complete presence of attenuation behaviors.Based on SA, attenuation characteristics are represented as response acceleration as a function of response period.Lin (2009) selected 60 earthquakes from 5968 records from the Taiwan Strong Motion Instrumentation Program (TSMIP).He determined groundmotion attenuation relationships for PGA and SA as a function of magnitude m, which is represented as follows: where y is the response acceleration for PGA or SA in g, R is the shortest distance to the rupture surface in kilometers, F NM is 1.0 for the source with normal mechanism, F RV is 1.0 for the source with thrust mechanism, V s 30 is the averaged shear-velocity from the ground surface to the depth of 30 m, 3) for M W > 6.3, C 1 to C 8 and H are constants and their corresponding values for PGA and SA are listed in Table 2. Source effects in the form of different focal mechanisms were considered.The results indicate the largest amplitude for thrust events and the smallest one for normal ones based on the same conditions (i.e., magnitude, distance, and site effects).In addition, these GMPEs analyzed the corresponding standard deviations for PGA and SA (σ ln y in Table 2).It is useful to understand the reliability of these equations.By considering the feasibility of PSHA in terms of SA and its reliability, we used these GMPEs for crustal events in this study.
Many studies have pointed out that ground-motion attenuation behaviors for subduction earthquakes and crustal earthquakes are different (Crouse et al., 1988 and references therein).In the offshore of northeast Taiwan, the Philippine Sea Plate subducts to the north.In southern Taiwan, the Eurasian Plate subducts to the east (Fig. 1).Therefore, it is necessary to consider different GMPEs for subduction earthquakes for both interface and intraslab events for PSHAs.Lin and Lee (2008) proposed the first GMPEs for subduction events in Taiwan and analyzed strong-motion records of subduction earthquakes obtained by TSMIP arrays and the Strong Motion Array in Taiwan Phase I (SMARTI) to establish the GMPEs for PGA and SA.The GMPEs are functions of magnitude m and their form is represented as follows: where y is the response acceleration for PGA or SA in g, R is the hypocentral distance in kilometers, H is the focal depth in kilometers, Z t is the subduction zone earthquake type (Z t = 0 for interface earthquakes, whereas Z t = 1 for intraslab earthquakes), C 1 to C 7 are constants and their corresponding values for PGA and SA are presented in Table 3.The corresponding standard deviations are also analyzed (σ ln y in Table 3).For this study, we used the GMPEs of Lin and Lee (2008) for the Taiwan subduction system.

Magnitude completeness and the declustering process
By analyzing the earthquake catalog, the smoothing Kernel function was used to evaluate the long-term seismicity rate.For reliability of the results, we considered the catalog during the period when the network recorded all of the earthquakes with a certain magnitude threshold (known as magnitude completeness, M c ).In the Taiwan region, a better seismic network that spanned the entire island of Taiwan was constructed after 1935.To check the temporal evolution of M c , the maximum curvature approach (Wiemer and Wyss, 2000) was employed.We determined that M c was smaller than 5.0 after 1940.For construction of a short-term rate change model and retrospective forecast, the catalog with a period is considered to test the feasibility of models (defined as "forecasting period").Since the rate perturbation duration for a M ≤ 7.0 earthquake could be years (Burkhard and Grünthal, 2009), the time window from 2006 to 2010 is assumed to be the "forecasting period".Earthquakes with M L ≥ 5.0 during the "forecasting period" (Fig. 1b) are introduced.On the other hand, we considered the time window from 1940 to 2005 to be the "learning period" and introduced earthquakes with M L ≥ 5.0 during the "learning period" (Fig. 1a) to establish the long-term rate model.
Since we modeled aftershock behavior with the rate-andstate friction model, the catalog in the time window from 1940 to 2005 used for the evaluation of the long-term seismicity rate needed to be declustered.We followed the approach of Burkhard and Grünthal (2009) for declustering.According to this approach, earthquakes are considered dependent when their distance and time are within the derived distance and time window, respectively.Distance window, dR(M w ), as a function of moment magnitude, M w , is represented as follows: The unit of dR(M w ) is in kilometers.The earthquakes are considered as foreshocks if the time intervals with mainshock are smaller than time window dTf (M w ), which is represented as follows: The unit of dTf (M w ) is in days.The earthquakes are considered as aftershocks if the time intervals with mainshock are smaller than time window dT a(M w ), which is represented as follows: dT a(M w ) = e (6.44+0.06•Mw ) for M w ≥ 6.6.
The unit of dT a(M w ) is in days.Although this approach was originally developed for conditions in central Europe, it has been proven to be well applicable to the Italian region (Chan et al., 2010), and even for a subduction environment (Suckale and Grünthal, 2009).We used this approach and below discuss the feasibility of its application to Taiwan.The catalog from 1940 to 2005 for evaluating the long-term earthquake occurrence rate is a declustered dataset, whereas the catalog used to calculate the short-term seismicity rate change (Table 1) is a raw one.The rate perturbations imparted by foreshocks, mainshocks, and aftershocks can be calculated using those procedures mentioned above.Thus, spatial-temporal evolution seismicity rate for a sequence can be modeled.3 Results

The long-term seismicity rate
For determining the long-term seismicity rate we first applied the smoothing Kernel function (Eq.2).Values for the bandwidth function could be determined using a linear regression of the catalog during the learning period (Fig. 3).The bin interval is 0.5 units for calculating the mean distance of the nearest event.We found that the c and d values of the bandwidth function were 2.18 and 0.41, respectively.With a small standard deviation of the regression (dashed lines in Fig. 3), it suggested that the regression results fit the observation well and a reliable forecasting result could be expected.
Here, we do not discuss the impacts of different magnitude bin intervals for regression since the purpose of this study focuses on seismic hazard assessment.Only 192 earthquakes with M L ≥ 5.0 in the declustered catalog (gray dots and gray stars in Fig. 1a) are considered for determining the long-term seismicity rate.A relatively small number of earthquakes were used in this study.It is difficult to discuss the impacts of different magnitude intervals for the seismicity rate model.Molina et al. (2001) recommended that the power law index, PL, in Eq. ( 2) should be between 1.5 and 2.0, corresponding to the cubic or quadratic decay of seismic activity within a hypocentral distance.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.
The long-term seismicity rate as a function of magnitude was calculated using Eq.(2) according to the catalog during the learning period (Fig. 4).Generally, higher seismicity rates are determined for smaller magnitude ranges (e.g., Fig. 4a) than for larger ones (e.g., Fig. 4d).This phenomenon corresponds to the Gutenberg-Richter law (Gutenberg and Richter, 1954).However, one of the highest seismicity rates was evaluated in southwestern Taiwan, corresponding to  frequent seismicity during the learning period (Fig. 1a), and can be associated with a high deformation rate according to the GPS observations (Yu et al., 1997).Another region with a high seismicity rate was found on the southeastern and eastern offshore regions.It corresponds to the boundary between the Eurasian and Philippine Sea plates (Fig. 1).
According to the magnitude-dependent bandwidth function (Fig. 3), the seismicity rate distributions between magnitude bins are slightly different from each other.The highest rates for magnitudes between 5.0 and 5.5 are evaluated in southwestern Taiwan (Fig. 4a).Those for the larger magnitudes (M L ≥ 5.5), by contrast, are on the southeastern and eastern offshore regions (Fig. 4b-d).

The short-term seismicity rate change
To model the distribution of aftershock sequences or triggered earthquakes, we used the CFS incorporated into the rate-and-state friction model.We assumed that small earthquakes or those that occurred far in the past did not significantly influence the current seismicity rate (Catalli et al., 2008;Chan et al., 2010).We considered earthquakes of M L ≥ 5.0 that occurred between 2006 and 2010 (Table 1) as source events for the rate change calculation.The focal mechanism and the depth of each event were determined using the moment tensor inversion listed on the website for the Broadband Array in Taiwan for Seismology (BATS, http://bats.earth.sinica.edu.tw/).For the CFS calculation, we considered an intermediate value of µ = 0.4.The value is in good agreement with a reasonable range of µ between 0.2 and 0.5, as inferred from the study of earthquake focal mechanisms in Taiwan (Hsu et al., 2010).Previous studies (Toda and Stein, 2003;Toda et al., 2005;Catalli et al., 2008) have suggested that the physically reasonable range for Aσ is between 0.1 and 0.4 bars, for application of the rate-and-state friction model.In this study we assumed a fixed Aσ of 0.2 bars.For aftershock duration, t a , we first assume it to be the function of magnitude proposed by Burkhard and Grünthal (2009) and below discuss the influence of the PSHA results.
The calculated seismicity rate changes for different moments are presented in Fig. 5.The results indicate that a significant seismicity rate increases near source events that have just occurred and then retreats to the long-term rate gradually with time, corresponding to aftershock decays succeeded by large earthquakes.It should be mentioned that the 2007 M L = 5.8 earthquake (event 6 in Table 1) took place approximately 10 km away from Hualien City and caused a seismicity rate increase of 153 % near the epicenter, which is approximately 10 km away from Hualien City.By comparing with the distribution of consequent events (both white and gray circles in Fig. 5), it is found that most of events are located in the region with seismicity rate increase (the area in red or yellow).For the rate drop region (the area in blue), by contrast, there are relatively few events occurred.

Both the long-term and short-term probabilistic seismic hazard assessment in Hualien City
To test the feasibility of the probabilistic seismic hazard assessment, we applied this approach to Hualien City, which is located in one of high seismicity regions in Taiwan (Fig. 1).By considering the magnitude of completeness for the catalog, we introduced the seismicity with M L ≥ 5.0 during the learning period .We introduced the GMPEs by Lin (2009) (Eq.12) and Lin and Lee (2008) (Eq.13) for crustal and subduction earthquakes, respectively.Based on these GMPEs, the averaged shear-velocity from the ground surface to the depth of 30 m (V s 30 ) is a requirement.We introduced the V s 30 in the TSMIP sites of HWA007-HWA019 (Lee and Tsai, 2008), which are in Hualien City.We obtained the averaged V s 30 of 555 m s −1 .Based on the information, seismic hazards in terms of probability of strong ground motion, spectral acceleration, and hazard curve are assessed.

The temporal evolution of PSHA
We calculated the temporal evolution of the seismic hazard in Hualien City from 2006 until 2010, corresponding to the forecasting period for construction of the short-term rate change model (black line in Fig. 6).We evaluated the seismic hazard for the 2.1 ‰ annual exceedance probability, corresponding to a 10 % probability of exceedance in 50 yr of exposure time (typical design lifetime of structures).Based on the long-term seismicity rate for M L ≥ 5.0 (Fig. 4) and corresponding GMPEs, a long-term hazard of 0.46 g was determined (black line in the beginning of 2006 in Fig. 6).Afterward, seismic hazard perturbations that followed the occurrence of source events (Table 1) were determined.Since the 2007 M L = 5.8 earthquake (event 6 in Table 1) was close to Hualien City and caused a significant rise in the seismicity rate in the vicinity (Fig. 5), such a significant rate increase will also raise seismic hazards in the surrounding regions following this earthquake.Such results could be associated with a possible aftershock sequence or even a consequent larger event following this event.During this period the largest earthquake occurred at the end of 2006 (event 4 in Table 1, which is as known as "the Pintung earthquake") (Ma and Liang, 2008).Since the epicenter was far from Hualien City (farther than 200 km) its influence on the seismic hazard for Hualien City was trivial (Fig. 6).Fig. 7.The acceleration response spectra in Hualien City for the probability of 2.1 ‰ as a function of the structural period of the two moments.The gray dashed line represents the long-term hazard curve (the short-term rate perturbation is not considered).The black solid line represents the hazard curve one day following the 2007 M L = 5.8 earthquake (event 6 in Table 1).

The spectral acceleration
Seismic hazards in the form of acceleration response spectra for Hualien City for two moments were evaluated (Fig. 7).As inferred from the long-term seismicity rate (Fig. 4) without considering earthquake perturbations, the long-term seismic hazard was relatively low (the dashed curve in Fig. 7).A significantly higher hazard of ≥1.0 m s −2 was obtained for the response period between 0.15 and 0.6 s.We obtained a significantly higher hazard one day following the 2007 M L = 5.8 earthquake (event 6) that was approximately twice as high as the long-term hazard for the entire spectrum of response periods.

The hazard curve
To show the probability of annual exceedance as a function of ground motion level for the response period of 0.2 s (SA 0.2 s), in Fig. 8 we present hazard curves.We compared hazard curves for the two moments and found the influence of the short-term seismic rate perturbation.Following the occurrence of event 6, more than three times a higher probability was determined when the ground motion level was lower than 1.0 g.In contrast, the contribution of the short-term perturbation was less significant for higher ground levels.This phenomenon can be attributed to the relatively small magnitudes of the earthquakes in the vicinity of Hualien city (include the M L = 5.8 event 6), which contribute lower probability for high ground motion levels.The temporal evolution of PSHA, spectral acceleration, and hazard curves are scenario results and will be difficult to compare with observations.
Long-term seismic hazard Short-term seismic hazard (1 day after event 6) Hazard curves for the response period of 0.2 seconds (in g) The gray dashed line represents the long-term hazard curve (i.e., the short-term rate perturbation was not considered).The black solid line represents the hazard curve one day after the 2007 M L = 5.8 earthquake (event 6 in Table 1).

The spatial distribution of seismicity in Taiwan
In this study, we evaluated the long-term seismicity density rate according to the distribution of the declustered catalog during the learning period.Our results indicated the highest seismicity rate in southwestern Taiwan for small events (Fig. 4a).However, there were only a few M L ≥ 6.0 events (the stars in Fig. 1a) that could result in high b values for southwestern Taiwan, as reported by Chan et al. (2012b).On the southeastern and eastern offshore regions, higher rates are expected for large earthquakes (Fig. 4b-d).
According to the non-declustering catalog (the transparent circles in Fig. 1a), more M L ≥ 5.0 earthquakes were observed on the southeastern and eastern offshore regions.The observations could be associated with more consequent events triggered by mainshocks.Such clustering behavior fulfills the observation of Zhuang et al. (2005), who calculated the ratio of the clustering seismicity rate to the total seismicity rate and found a high ratio for this region due to related high aftershock activities triggered by mainshocks.Such clustering behavior may imply the importance of earthquake sequences in seismic hazard evaluations for this region.

A comparison to the traditional PSHA approach
The application of traditional PSHA approaches (Cornell, 1968;McGuire, 1976) requires knowledge of the properties of the seismic source zones, which are determined by subjective judgments that may be different in various studies.In this study, the earthquake catalog was the only input for the long-term PSHA.Despite that fact that our approach only considered a few factors, we concluded similar hazards as for evaluations using traditional approaches.We obtained a background seismic hazard of 0.46 g for the peak ground acceleration for a 2.1 ‰ annual exceedance probability (black line in the beginning of 2006 shown in Fig. 6).The result is similar to those obtained from previous studies (Cheng, 2002;Cheng et al., 2007;Campbell et al., 2002).

The feasibility of seismicity rate models
An accurate model for the seismicity density rate plays a key role in obtaining a reliable PSHA.We evaluated the longterm seismicity density rate and the short-term rate change based on a smoothing Kernel, and the rate-and-state friction model, respectively; and proposed a retrospective forecast in order to validate their feasibility.We compared models with the distribution of non-declustering earthquakes during the forecasting period (denoted as "forecasting earthquakes" in Fig. 1b) using the Molchan diagram (Molchan, 1990(Molchan, , 1991)).The diagram was designed for an evaluation of earthquakeforecasting ability and is presented as the fraction of space occupied by the expectation versus the fraction of failure to predict.We present the "fraction of space occupied by alarm" as the proportion of the study area having a forecasted seismicity rate equal to or higher than the threshold, defined as "alarm".The "fraction of failure to predict" indicates the proportion of forecasted earthquakes that had a lower forecasted seismicity rate than the alarm.In other words, when data points were distributed along the diagonal line, the distribution of forecasting earthquakes was uniform or independent of the forecasted rate.When a convex was presented it suggests that the majority of forecasting earthquakes occurred within regions with a lower forecasted rate as compared to the entire area.On the other hand, when a concave was presented it suggests that the majority of the forecasting earthquakes occurred in an area with a higher forecasted rate.An optimistic forecast result represents a condition of having the least space occupied by alarms and the lowest percentage of forecasting earthquakes with a failure to predict.When we compared the seismicity rate obtained using the smoothing Kernel function with the locations of forecasting earthquakes (as shown in Fig. 4) in the Molchan diagram (the open circles in Fig. 9), a positive correlation between each was evident.All of the forecasting earthquakes were located within 61 % of the study area having the highest seismicity rate.When we compared the forecasted seismicity rate evaluated from a combination of the smoothing Kernel function and the rate-and-state friction model (the crosses in Fig. 9) we found somewhat better forecasting ability in comparison to that from using solely the smoothing Kernel function.The results indicate the importance of considering a short-term rate perturbation in order to improve the seismicity rate model.Fraction of space occupied by alarm Fig. 9.A Molchan diagram that investigates the correlation between seismicity rate models and non-declustering seismicity during the forecasting period (2006)(2007)(2008)(2009)(2010).Gray circles represent the correlation with the long-term rate model using the smoothing Kernel function.Black crosses represent the correlation with the short-term rate model according to a combination of the smoothing Kernel function and the rate-and-state friction model.

The feasibility of the declustering approach
We modeled seismicity rate evolution following the occurrence of large events according to the rate-and-state friction model.To ward off the reduplication of earthquake interactions in the long-term rate model, we considered a declustering catalog according to the approach proposed by Burkhard and Grünthal (2009).In order to test the feasibility of the declustering approach applied to Taiwan, we compared seismicity rate models with the distribution of declustering earthquakes during the forecasting period (as shown in Fig. 1b).Therefore, it is expected that when consequent events are triggered by previous earthquake(s), the model considering short-term rate evolution represents better forecasting ability.However, the model was found to have better forecasting ability without the consideration of a short-term rate change in the Molchan diagram (Fig. 10).The result suggests that earthquakes in the declustering catalog are independent of one another, and indicates that the declustering approach by Burkhard and Grünthal (2009) is applicable for the Taiwan catalog.

Impacts of the aftershock duration
We evaluated the temporal evolution of the seismic hazard in Hualien City from 2006 until 2010 by assuming aftershock  duration, t a , as a function of the magnitude proposed by Burkhard and Grünthal (black line in Fig. 6).To discuss the impact of t a , we assumed this parameter to be 10 yr and 10 days and evaluated corresponding evolution of seismic hazards (dashed and light-gray lines, respectively, in Fig. 6).When a short t a is assumed, the seismic hazard perturbations are transient and trivial following source events; and vice versa the assumption of a long t a results in longlasting and significant perturbations.According to the statistical test through the Molchan diagram (Fig. 10), it suggests that the aftershock duration proposed by Burkhard and Grünthal (2009) is applicable for the Taiwan catalog.However, this parameter could be further constrained through different aspects, such as focal mechanisms of earthquake sequences, loading rates of source faults/regions (Stein and Liu, 2009), etc.

The importance of the short-term probabilistic seismic hazard assessment
According to a traditional PSHA (Cornell, 1968;McGuire, 1976), the seismicity rates generally represent long-term behavior of seismic sources and do not include short-term variations in their rates due to the time since the last large earthquake occurrence.However, several experiences have suggested that aftershock durations can be longer than the designed lifetime of structures (typically for 50 yr) and demonstrated the importance of fault interactions for the evaluation of seismic rate and seismic hazard.For example, after the 2010 Darfield, New Zealand, mainshock, the M w = 6.3, 2011 Christchurch earthquake took place 5 months after the mainshock and resulted in severe damage and fatality in downtown Christchurch (Chan et al., 2012a).Another example for longer-period perturbation is the 4 March 2010 M w = 6.4 Jiashian, Taiwan, earthquake sequence (Chan and Wu, 2012).A significantly high seismicity rate is found in the southern Taiwan region, where an increase of 170 % in seismicity rate is modeled in the vicinity of the 26 February 2012 M w = 6.1 Wuatai epicenter at the end of 2012.Furthermore, Stein and Liu (2009) concluded that durations of aftershock sequences could be sustained for from days to centuries.So, it is important to consider short-term rate perturbation for PSHA.
One possible approach is by consideration of deterministic scenarios caused by aftershock sequence.The behaviors of aftershock sequence can be integrated into seismic hazard assessments by considering models such as the modified Omori's law (Utsu, 1961;Utsu et al., 1995) and Båth's law (Båth, 1965).In these models, however, it is difficult to evaluate the probability of consequent earthquakes with larger magnitudes, such as the foreshock-mainshock sequence off the Pacific coast of Tohoku, Japan, in 2011 (Nettles et al., 2011).
In our approach, short-term rate change is controlled by the rate-and-state friction model and the period of seismicity rate perturbation depends on aftershock duration, t na , and the magnitude of CFS.It is expected that a higher hazard would sustain for a longer period if longer aftershock duration, t na , is assumed (dashed line in Fig. 6) or/and a larger event takes place.Additionally, based on this approach, the foreshock-mainshock behaviors can also be modeled, i.e., a source event can trigger consequent events with larger magnitudes.The feasibility of this approach has been tested through the applications to the 2010 Jiashian, Taiwan, earthquake sequence (Chan and Wu, 2012), the 2010 Darfield, New Zealand, earthquake sequence (Chan et al., 2012a), and the 1906 Meishan, Taiwan, earthquake sequence (Chan et al., 2013).

Possible applications of the short-term probabilistic seismic hazard assessment
Based on the ability of updating information, the short-term probabilistic seismic hazard assessment could be useful for decision-makers and engineers.One possible application is to provide probabilistic seismic hazard/risk information after the occurrence of a devastating earthquake.Our results show that we can assess the probability as a function of time and assuming ground motion levels (Fig. 11).probabilities for PGA, SA 1.0 s and 0.2 s are 0.1, 0.5, and 0.9 %, respectively.After the occurrence of event 6 (Fig. 11), the corresponding probabilities become 0.9, 3.6, and 7.1 %, respectively.About six to eight times of probability is higher than estimation before the influence.Since this approach can be widely applied for each administrative region, when the information of vulnerability (fragility curve) and exposure (distribution of structures and population) is further considered, a short-term probabilistic seismic risk map can be assessed.Such information is important to provide emergency responses regarding victim sheltering and relocation.Another possible application is the combination with deterministic earthquake scenarios.Scenarios are also useful for hazard and risk mitigation purposes, in terms of building codes and the site selection of pubic structures (McGuire, 2001).Probability of occurrence and characteristics of a single (scenario) earthquake can be re-evaluated after a large earthquake.In a formal way, such estimations may be done using deaggregation and determination of dominant earthquakes.Our approach provides different insights for seismic hazard.It would be useful as an alternative credible model in the logic-tree approach for PSHA.

Conclusions
Traditional PSHA approaches (Cornell, 1968;McGuire, 1978) are often applied, and are required in order to obtain some of the assumptions and information for seismic source zones, which are highly subjective and time consuming.Therefore, understanding the seismic hazards that result from a large earthquake cannot be rapidly evaluated.However, using the approach proposed in this study employing a smoothing Kernel and the rate-and-state friction model, PSHAs could be evaluated very quickly following the occurrence of a large earthquake.We demonstrated the method's feasibility by applying the results to Hualien City, Taiwan.Using minor model assumptions, our approach provides a suitable basis for rapid evaluations.The concept of a near real-time and updating seismic hazard assessment benefit decision-makers and public officials in seismic hazard mitigation.
evaluated seismic hazards for Taiwan and C.-H. Chan et al.: Time-dependent probabilistic seismic hazard assessment

Fig. 1 .
Fig. 1.The seismicity during the following: (a) the learning period (1940-2005), and (b) the forecasting period (2006-2010).The declustering and non-declustering earthquakes are shown with light gray and transparent colors, respectively.Earthquakes with M L ≥ 5.0 and M L ≥ 6.0 are shown with circle and star shapes, respectively.Red lines represent the Longitudinal Valley fault, the boundary between the Eurasian and Philippine Sea plates.Dashed lines represent the shapes of the subduction zones.The earthquake parameters are determined by the Central Weather Bureau Seismic Network (CWBSN).

Fig. 2 .
Fig. 2. (a) The reference focal mechanisms collected by Wu et al. (2010).(b) The focal mechanisms as spatial variable receiver faults for the CFS calculation.Note that the actual calculation grids are denser than the spacing presented in the figure.

Fig. 3 .
Fig. 3.The bandwidth functions (solid line) and the corresponding standard deviation (dashed lines) obtained by the linear regression of observations (black dots).

Fig. 6 .
Fig.6.Assumed aftershock durations, t a , as theBurkhard and Grünthal (2009) relation (black line), 10 yr (dashed line), and 10 days (light-gray line); the temporal evolution of the seismic hazard for the 2.1 ‰ annual exceedance probability in Hualien City from 2006 until 2010.Variations in the seismic hazard can be attributed to the seismicity rate change imparted by source events (Fig.5).

Fig. 8 .
Fig.8.The probability of annual exceedance as a function of hazard curves for the response period of 0.2 seconds at two moments.The gray dashed line represents the long-term hazard curve (i.e., the short-term rate perturbation was not considered).The black solid line represents the hazard curve one day after the 2007 M L = 5.8 earthquake (event 6 in Table1).

Fig. 10 .
Fig. 10.A Molchan diagram that investigates the correlation between seismicity rate models and declustering seismicity during the forecasting period (2006-2010).Circles represent the correlation of the rate model using the smoothing Kernel function.Crosses represent the correlation of the short-term rate model according to a combination of the smoothing Kernel function and the rate-and-state friction model.

Fig. 11 .
Fig. 11.Assumed ground shaking of 0.6 g for PGA (black line), SA1.0 s (gray line), and 0.2 s (light-gray line), the temporal evolution of the annual exceedance probability in Hualien City from 2006 until 2010.

Table 1 .
Source parameters for source events used to calculate the seismicity rate change with the rate-and-state friction model.The earthquakes of M L ≥ 5.0 that occurred from 2006 to 2010 were considered.The parameters of focal mechanism are determined by the Broadband Array in Taiwan for Seismology (BATS, http://bats.earth.sinica.edu.tw/).

Table 2 .
Lin (2009)ing parameters and standard deviations (σ ln y ) of the GMPEs for PGA and SA of various response periods for crustal events (Eq.12) byLin (2009).

Table 3 .
Lin and Lee (2008)meters and standard deviations (σ ln y ) of the GMPEs for PGA and SA of various response periods for subduction events (Eq.13) byLin and Lee (2008).
Here we represent three examples of 0.6 g ground shaking for PGA, SA 1.0 s, and 0.2 s.The long-term (i.e., at the beginning of 2006)