Taiwan Earthquake Occurrence Probability Estimation from Regional Source Model Since 1900

Taiwan is located on the boundary between the Eurasia Plate and Philippine Sea Plate, which is a very high seismicity rate area. We begin calculating the earthquake occurrence probability using four recurrence models to mitigate seismic disasters. We focus on estimating the occurrence probabilities for regional earthquake sources based on the catalog released by the Central Weather Bureau over the period from 1900 to 2011. According to the tectonic and seismicity characteristics areas in and around Taiwan are divided into several zones for shallow and deep earthquakes. We utilize four recurrence models to estimate the earthquake occurrence probabilities over the next 30 and 50 years, respectively. In addition, the grid-based probabilities in 0.1° × 0.1° spatial size are calculated using the micro-zoning method. The results obtained from four recurrence models show that areas with high values over the next 30 and 50 years are correlated with two subduction zones and a suture zone. High probabilities in the western foothills appear highly active tectonic. Moreover, the high values appear at in Eastern Taiwan, offshore Hualien County. There are discrepancies between the results from the four models. The highest grid-based probability is about 3.0, 3.5, 2.5 and 3.5% for the Lognormal, Gamma, Exponential, and Weibull models, respectively. The inland probabilities are below 0.5% for the results from Lognormal, Gamma, and Weibull models. Even so, the results from the Exponential model are upmost in the range from 0.5 1.0%.


INTRODUCTION
Taiwan is located on an arc and continent collision zone, which is the boundary between the Eurasian plate (EP) and the Philippine Sea plate (PSP) (Fig. 1).The seismicity rate is very high in and around the Taiwan region.Several disastrous earthquakes occurred in the recent past (Tsai 1985;Cheng and Yeh 1989).After a disastrous event on September 21 st , 1999, the issue regarding earthquake prediction became popular.However, there are still difficult problems that must be conquered because the earthquake mechanism occurrence is not clear enough.Calculating the earthquake occurrence probabilities is important work for disaster reduction.Hence, statistical methods are utilized to analyze the earthquake occurrence probabilities.The results can be applied to long-term seismic hazard mitigation.
The National Science and Technology Center for Disaster Reduction (NCDR) began a campaign against the disasters brought by large-scale earthquakes in recent years.Due to limited time and cost resources, it is suggested that regions prone to devastation by large-scale earthquakes be considered as the first step.Those characterized with high earthquake disaster potential due to political and economic importance should be assigned as priority regions where mitigation resources should be promptly applied.Once the priority regions are assigned, the seismic scenario simulation technology can be applied to evaluate the potential damage and losses, i.e., buildings collapsed and casualties.Accordingly, a corresponding policy and strategy to cope with disasters can be developed to improve the seismic safety of buildings and guarantee government functionality and business operations.The hazard mitigation plan for largescale earthquakes is shown in Fig. 2.
The priority region designation mechanism is to organize a working group that stresses the tasks of surveying the fault parameters, estimating the earthquake occurrence probabilities and evaluating the seismic risks for Taiwan.Results from the working group will conclude the priority regions for seismic hazard mitigation.To achieve this goal the working group is composed of three forums (Fig. 3).There are duty institutes to invite relevant scholars and experts to participate in each of the forums.The first one hosted by the Central Geological Survey (CGS) is the Fault Parameter Forum which proposes applicable fault parameters.The second one hosted by NCDR and the Central Weather Bureau (CWB) is the Earthquake Potential Forum which calculates the occurrence probabilities for regional source earthquakes and predicts the ground shaking of a specific earthquake.The third one hosted by the National Center for Research on Earthquake Engineering (NCREE) is the Risk Evaluation Forum which assesses the possible damage and losses induced by earthquakes in Taiwan.According to the results from these three forums, an assembly of expertise will be held to discuss and conclude the priority regions.In the following text we focus on the results from the second forum.
In general, earthquake sources can be classified into two groups depending on whether fault ruptures can be ob-served or investigated in the field.The first one is earthquakes associated with the faults that rupture to the surface, also called fault sources (Ma et al. 1999;Chen et al. 2001).The second one belongs to earthquakes regarding undefined faults, namely blind faults (also called regional sources).Due to insufficient information on fault sources at present we focus on evaluating the occurrence probability for regional sources under the assumption that all earthquakes occur independently, known as the Poisson process (Utsu  2002).In fact, the status of regional tectonic tends to be stability stress after a main shock has occurred, leading to aftershocks.In order to follow the Poisson process assumption, the aftershock sequences must be removed from the catalogs (Gardner and Knopoff 1974;Knopoff 2000;Kagan 2002).After removing the aftershock sequence the returned periods for each region will be calculated and used as inputs in the earthquake recurrence models.Based on the catalogs over the period from 1900 to 2011 released by the CWB, we then estimate the earthquake occurrence probabilities over the next 30 and 50 years, respectively.We also utilize the micro-zoning method to elaborately calculate the occurrence probabilities for each grid with a size of 0.1° × 0.1° covering the Taiwan area over the next 50 years to meet seismic hazard mitigation needs.

SEISMICITY OF TAIWAN
Prior to 1897 the historical records on earthquakes appear primarily in either local government documents or personal diaries.Certainly, the records are definitely incomplete.The seismograph installed in Taiwan since 1897 begins the instrumental seismicity observations.The early observation instruments were few in number and had a low magnification, so only earthquakes strong enough to be felt by people could be detected.Only those earthquakes with magnitude greater than 5.5 that caused substantial damage were recorded.
The records for smaller earthquakes with magnitudes less than 5.5 are incomplete prior to 1935.It is reasonable to believe that the record is complete for earthquakes with magnitudes greater than or equal to 6.0 during this period.A great emphasis on seismological observation was given following the destructive event of April 21, 1935 in Central Taiwan.Since then the seismic network has been improved to increase the detection and location capability for earthquakes.However, the catalogs for earthquakes with magnitudes less than 5.0 during this period are still incomplete.Since 1972 the Institute of Earth Sciences, Academia Sincia has developed the Taiwan Telemetered Seismograhic Network (TTSN).The earthquake location and magnitude resolution has been enhanced since then.The magnitude threshold was lowered to 2.0.
Different magnitude scales were used in the catalogs reported by Lee et al. (1976), Hsu (1980), andTsai (1985).Since the local magnitude, M L is more pertinent for earthquake engineering purposes, Dr. Shih-Nan Cheng (2012 unpublished result) compiled an earthquake catalog of instrumental observations in the M L scale for consistency.

ESTIMATE OF OCCURRENCE PROBABILITY
As mentioned above earthquakes can be classified into regional and fault sources.We used regional sources in this study to estimate earthquake occurrence probabilities based on four general statistical models used frequently for earthquake recurrence (Utsu 2002).Figure 4 shows a flowchart of earthquake occurrence probability evaluations.We will discuss the processes in detail in this section and also show an example to describe the process and methodology.The example assumes the total Taiwan area as a single region.

Regional Source Occurrence Probability
Because shallow earthquake sources are different from deep sources we divided the earthquake catalog into two groups according to the focal depth.In traditional seismic hazard analysis it is suggested that the earthquake focal depths 35 km or less are shallow source.In order to reduce the influence on the location precision by the Taiwan Central Weather Bureau Seismic Network (CWBSN), we defined the source depth ≤ 40 km as shallow source and source depth ≥ 35 km as deep source.There are 5-km focal depth buffers to distinguish between shallow and deep sources.

Earthquake Catalog
We used the CWB unified earthquake catalog, for which the location and magnitude of earthquakes are revised (Shin 1993;Yeh et al. 1995;Cheng 2012 unpublished result).The earthquake magnitude is the local magnitude, M L .The catalog period 1900 to 2011 is used in this study.

Removing the Effect of Aftershock
All earthquakes occur independently under the Poisson process assumption.We adopted a method to remove the aftershock effect.The aftershock spatiotemporal distribution depends on the main shock magnitude (Gardner and Knopoff 1974;Knopoff 2000;Kagan 2002).We choose the method proposed by Gardner and Knopoff (1974) in form of Eqs. ( 1) and (2).
where M is the magnitude, T is time in day units and L is distance in km units.Gardner and Knopoff (1974) utilized the Southern California catalog to conclude the a 1 , b 1 , a 2 , b 2 as -2.007, 0.725, -0.567, 0.347 values, respectively.Figure 5a shows the relationship between the magnitude and distance in the form of a solid line with circles.Figure 5b displays the relationship between the magnitude and time in the form of a solid line with circles, in which the time exceeds 1000 days while the magnitude is at 7.0.The dependence on magnitude parameters are not suitable for the Taiwan area.Therefore, we adopted the equations modified by Dr. Tzay-Chyn Shin (personal communication), who is a committee member of the Earthquake Potential Forum hosted by NCDR.The equations can be expressed as follows: Therefore, the solid line with squares in Fig. 5b displays the relationship between the magnitude and time in these equations.According to Eq. ( 4) we have a time span over 500 days while the magnitude is equal to 7.An example shown in Fig. 6a is a sequence of earthquakes including aftershocks.According to Eqs. ( 3) and ( 4), lots of aftershocks are filtrated from the sequence in Fig. 6b.

Zoning
According to the seismicity in and around the Taiwan area there are strong correlations between the earthquake locations and tectonic process.As mentioned above, the EP and PSP collision results in active seismic areas, such as two subduction zones, suture zone and the accretion wedges.Tsai et al. (1977) stated that the boundary between EP and PSP in Northern Taiwan lies parallel with the 121.5°E meridian.The dipping angle of the subduction zone is about 45° -50°.Therefore, lots of earthquakes with deep focal depth occur in the northeastern part of Taiwan in the area east of 121.5°E and north of 24°N.The seismicity occurs Fig. 4. Procedures for calculating the occurrence probabilities for regional sources.mostly in the western foothill province where the fold-andthrust of imbricate structures and thin-skinned tectonics control the earthquake behaviors (Suppe 1987).
Based on the geological, tectonic structures, subduction models and seismological information, the seismogenic zones were selected to a 12 and 4 for shallow sources and deep ones, respectively, as shown in Fig. 7.The zoning schemes are similar to that of Loh et al. (2004), which were designed to evaluate the seismic hazards for the Taiwan Power Company nuclear power plant.In Fig. 7a the gray dots denote the shallow earthquake distribution with M L ≥ 5 over a period from 1900 to 2011, while the solid lines depict the 12 regions with codes of S01, S02, etc.The deep source occurrence probabilities for 4 regions are evaluated.In Fig. 7b the dots denote the deep earthquake distribution with M L ≥ 5 over a period from 1900 to 2011, while the solid lines depict the 4 regions with codes of D01, D02, etc.

Magnitude of Maximum Potential Earthquake
From the earthquake catalogues we can estimate the maximum earthquake magnitude in the whole area as well as in each subzone (Iida 1976;Makropoulos and Burton 1983).Iida (1976) proposed the relation between the magnitude and strain energy, E, to evaluate the upper bound magnitude in form of Eq. ( 5): The maximum magnitude evaluation, M u, c , can be obtained through two enveloped lines for which the released strain energy is cumulated with time.The slopes of the two lines assume the energy is released constantly and repeatedly in this region providing long enough observations to average out large short-term fluctuations.The M u, c value can be determined under the assumption of entirely released strain energy between two enveloped lines.As shown in Fig. 8 the M u, c value is the maximum magnitude from Eq. ( 5).The M u, r value is the maximum magnitude from earthquake catalogues over the period from 1900 to 2011.Tables 1 and 2 show the M u, c and M u, r values in each zone for shallow sources and deep ones, respectively.Comparing the estimated maximum magnitude to the recorded maximum magnitude for each subzone, the estimated maximum magnitude appears to be a reasonable value.The seismic hazard analysis should take greater consideration into this evaluation.Hence, we choose the larger one between M u, c and M u, r as the potential magnitude for the specific zone, in this case the potential magnitude is 7.7.

Magnitude and Frequency Relationship
The seismicity statistics following the Gutenberg-Richter relationship (G-R relationship) (Gutenberg and Richter 1944) describes the relation between magnitude and frequency in the form of Eq. ( 6): where N(≥ M L ) expresses the frequency of earthquakes with magnitudes larger and equal to M L .Under the seismology and earthquake engineering consideration and to avoid catalog incompleteness, we take the m o value as a lower limit for the hazard estimation and the m u value as an upper limit for a potential earthquake in a specific region, respectively.In this study the m o value is set to 4. The coefficient of a and b can be calculated from the least square fit of Eq. ( 6).Because the earthquake catalogs in Taiwan are not complete a two-step method is used to estimate the a and b values.In the first step we selected a short term earthquake catalog (over the period of from 1973 to 2011) to estimate the b value using According to the Probability Seismic Hazard Analysis approach (PSHA) (Cornell 1968) we can develop the relationship between magnitude and the return period, T r .However, the G-R relationship only resolves the mean value of T r , but not the coefficient of covariance, COV(σ T /T r ), where the σ T value is the standard deviation of interevent time obtained using statistical techniques.If the earthquake events are incomplete the COV(σ T /T r ) values will be unreliable.To reduce the consequence of catalog incompleteness we did not evaluate the occurrence probabilities for shallow earthquakes for M L < 5 and those for deep earthquakes for M L < 3, respectively.

Recurrence Models
According to the cumulative energy concept the stress energy accumulates with time elapsed if no event occurs.The occurrence probability for the next earthquake increases with the time elapsed.The earthquake occurrence probability is evaluated using a recurrence model in form of Eq. ( 7) as ( , ) ( ) where P is the occurrence probability for an earthquake of magnitude (E m ) in next T p year; and T e is the elapsed time.
A cartoon shown in Fig. 10 displays how to evaluate the probability in the future in which T c and T 0 are the present time and last earthquake occurrence time, respectively.The value of T e can be obtained by subtracting T 0 from T c .As shown in Fig. 11 the vertical lines denote the occurrence time corresponding with the earthquake magnitude over the period from 1900 to 2011.The thick line displays an example to obtain the elapsed time for a specific magnitude.In this case finding out the elapsed time of M L = 4.0 we start at the present time and track backward in time to meet the next event with M L ≥ 4.0.The time interval can be obtained by subtracting between every meeting event.

Occurrence Probability
There are several statistical models frequently utilized in evaluating earthquake recurrence (Anagnos and Kiremi-djian 1988;Utsu 2002).The recurrence models are selected given the probability density function (PDF) to evaluate the occurrence probabilities.For comparison, we selected four models, such as log-normal, exponential, Weibull and  gamma distribution to evaluate the occurrence probabilities, respectively.The adopted recurrence model follows the log-normal distribution given the PDF, f(x, μ, σ), and can be expressed as follows: where x is the interevent time, m and v are the mean value and the standard deviation of f(x, μ, σ), respectively.We assumed the return period (T r ), it can be found from the G-R relationship and COV(σ T /T r ) obtained from interevent time of earthquake using statistical methods.
The gamma distribution has two parameters, α and β, which are derived from the m and v values.These equations can be expressed as follows: where Γ(x) represents the incomplete gamma function.
We selected an exponential distribution for testing.Its  PDF has only one parameter, β, which can be obtained using the mean value m.The model can be expressed as follows: The fourth model is selected as the Weibull distribution in which there are two parameters, λ and k in its PDF.The parameters, λ and k, are determined by m and v through the gamma function.The equations can be expressed as follows: ( , , ( ) As shown in Fig. 12 future earthquake occurrence probability results for the entire Taiwan area are shown based on the four recurrence models.The diagrams (a), (b), (c), and (d) show the results for the log-normal, Gamma, Exponential and Weibull models, respectively.The lines with squares, circles, diamonds, inverted triangles and stars represent the occurrence probabilities over the next 10, 20, 30, 40, and 50 years from 2011, respectively.The curves vary decreasingly with the magnitude increasing.It is understandable that larger magnitude earthquakes occur less frequently compared to small earthquakes.

Occurrence Probability of Micro-Zoning Method
The regional source occurrence probability in the aforementioned context shows the seismic potential for the whole region and each subzone, which do not meet practical application needs for the seismic hazard assessment.We applied the micro-zoning method to calculate the occurrence probabilities for each grid.The entire region is divided into smaller sized grids.The occurrence probability evaluation process is based on the micro-zoning method, refer to   13.This is similar to the regional source process in Fig. 4, but some steps are different from the regional source case.In step 3 the grid size is 0.1° by 0.1° (Fig. 14) and the grid center is moved an interval of 0.1° at each time (Fig. 15).Due to the location precision problem for the early period seismic observation system, we considered the events covered using a gird size of 0.2° by 0.2°, adjusting the completeness.In step 5 we implemented the relation between the magnitude and frequency for each grid.The b value of each grid comes from the regional source results (grid within the region), but the a value is estimated from the grid catalog.

RESULTS
The occurrence probabilities over the next 30 and 50 Fig. 13.Assessment processes of evaluating the occurrence probabilities for the micro-zones.years for earthquakes with M L ≥ 6.0, ≥ 6.7, and ≥ 7.0 are calculated, respectively.

Occurrence Probability from Region Sources
Table 3 shows the a and b values of the G-R relationship using the two-step-regression method in each zone for the shallow sources.The relation between the magnitude and return period can be derived from the a and b values.As mentioned above, the catalog completeness is adequate for magnitudes greater than 2.0 after 1973.We derived the reliable b values in the first step using the short term data over the period from 1973 to 2011 with M L ≥ 2. In the second step we corrected the b value to obtain the a value with long term data with M L ≥ 4. Given the occurrence time for the last event the occurrence probabilities over the next 30 and 50 years can be evaluated based on the recurrence model.We defined three earthquake magnitudes, 6.0, 6.7, and 7.0 as the potential magnitude.The occurrence probabilities with M L ≥ 6.0, ≥ 6.7, and ≥ 7.0 over the next 30 and 50 years are listed in Table 4.The asterisks in Table 4 denote the estimated M u, c values less than the potential magnitude.For earthquakes with M L ≥ 7.0, the highest probability is the S04 region resulting in 45.2, 43.9, 49.1, and 43.9% over the next 30 years from the four statistical models.The S10 region is an inland with high population, leading to more hazardous earthquakes with M L ≥ 7.0.The probability exceeds 16% over the next 50 years.
Conversely, the S01 and S12 regions belong to areas with lower seismic hazard.
Table 5 shows the a and b values of G-R relationship in each zone for the deep sources.The b values are all less than 1.The earthquake occurrence probabilities with M L ≥ 6.0, ≥ 6.7, and ≥ 7.0 over the next 30 and 50 years are listed in Table 6.The highest earthquake occurrence probability with M L ≥ 7.0 appears in the D01 region.The probability rate exceeds a value of about 50% over the next 50 years.The Okinawa trough stretches its length to region D01, which belongs to the PSP intra-slab.The D04 region belongs to an area with lower deep hazard seismic sources.

Occurrence Probability from Micro-Zones
The intention behind this study to estimate the occurrence probability utilizing the micro-zoning technique is to focus on locations where the highest probability occurs.The occurrence probabilities over the next 50 years are calculated based on the four probability models for the shallow sources.As shown in Fig. 16 diagrams (a), (b), (c), and (d) represent the spatial probability for earthquakes with M L ≥ 6.0 based on the Log-normal, Gamma, Exponential, and Weibull models, respectively.Obviously, the probability exceeding 50% is located at Eastern Taiwan offshore.There are high occurrence probabilities areas located on the boundary between the Eurasia Plate and the Philippines Sea Plate. Figure 17 shows the earthquake occurrence probability with M L ≥ 7.0 based on the four probability models.similar to that from the Weibull model.There is a discrepancy in the results from the four models.The highest probability is about 3.0, 3.5, 2.5 and 3.5% for the Lognormal, Gamma, Exponential, and Weibull models, respectively.The inland probabilities are below 0.5% for the Lognormal, Gamma, and Weibull model results.However, the inland probabilities from the Exponential model are upmost in the range of from 0.5 -1.0%.The occurrence probabilities over the next 50 years are calculated based on the four probability models for the deep sources.As shown in Fig. 18 diagrams from left to right (a), (b), (c), and (d) represent the spatial probability for earthquakes with M L ≥ 6.0 based on the Log-normal, Gamma, Exponential, and Weibull models, respectively.The high probabilities appear at northeastern Taiwan, offshore Ilan County.There are high occurrence probabilities areas located on the boundary, where the Philippines Sea Plate subducts beneath the Eurasia Plate. Figure 19 shows the earthquake occurrence probability with M L ≥ 7.0 based on the four probability models.The high probability values are located at Eastern Taiwan, offshore Ilan County.There is a discrepancy in the results from the four models.The highest   probability is below 0.5% for the Lognormal, Gamma, and Weibull models.Only the results from the Exponential model can reach a value of 1.0%.The inland probabilities are almost 0.0% for the results from the four models.It is reasonable that the seismicity concentrates in the shallow zone associated with the fold-and-thrust of imbricate structure and thin-skinned tectonics.

CONCLUDING REMARKS
This study calculated the regional earthquake source oc-currence probabilities classified into two groups, i.e., shallow and deep earthquakes utilizing four statistical models over the next 30 and 50 years.The probabilities in scale of zones and micro-zones are revealed for M L ≥ 6.0, ≥ 6.7, and ≥ 7.0.The high occurrence probabilities are associated with active tectonics, such as the fold-and-thrust of imbricate structure and thin-skinned tectonics for the inland areas and subduction between PSP and EP for the offshore areas.The regions or areas with high occurrence probabilities are identified and provided for seismic hazard mitigation plans in Taiwan.Estimating the seismic potential for Taiwan is a preliminary task.It is suggested that the mitigation work should utilize new technologies and data to revise the seismic potential map every 3 to 5 years.

Fig. 2 .
Fig.2.The hazard mitigation planning scheme for large-scale earthquake disasters.

Fig. 5 .Fig. 6 .
Fig. 5.The spatiotemporal aftershock dependence on main shock magnitude: (a) The dependent distance in various magnitudes; (b) The dependent duration in various magnitudes.

Fig. 7 .
Fig. 7.The distribution of earthquakes with M L ≥ 5 over a period of from 1900 to 2011.The gray dots denote the distribution of earthquakes and the solid lines depict the regions, respectively, for (a) the shallow sources and (b) the deep sources.

Fig. 8 .
Fig. 8.The strain energy accumulated with time over a period of from 1900 to 2011.

Fig. 9 .
Fig. 9.The relationship between magnitude and frequency.The diamonds and squares denote the data from short term and long term catalog, respectively.The solid lines are the best-fitting results.

Fig. 10 .
Fig. 10.A cartoon illustrates the assessment of the occurrence probability.

Fig.
Fig. 13.This is similar to the regional source process in Fig.4, but some steps are different from the regional source case.In step 3 the grid size is 0.1° by 0.1° (Fig.14) and the grid center is moved an interval of 0.1° at each time (Fig.15).Due to the location precision problem for the early period seismic observation system, we considered the events covered using a gird size of 0.2° by 0.2°, adjusting the completeness.In step 5 we implemented the relation

Fig. 16 .
Fig. 16.The occurrence probabilities (%), shown by the color bar, of shallow zone (depth ≤ 40 km) with M L ≥ 6.0 for the four recurrence models over the next 50 years.

Fig. 17
Fig. 17.The occurrence probabilities (%), shown by the color bar, of shallow zone (depth ≤ 40 km) with M L ≥ 7.0 for the four recurrence models over the next 50 years.
Fig. 17.The occurrence probabilities (%), shown by the color bar, of shallow zone (depth ≤ 40 km) with M L ≥ 7.0 for the four recurrence models over the next 50 years.
Fig. 18.The occurrence probabilities (%), shown by the color bar, of deep zone (depth ≥ 35 km) with M L ≥ 6.0 for the four recurrence models over the next 50 years.

Table 1 .
Maximum magnitude of potential earthquake for each shallow zone.

Table 3 .
The Coefficients a and b, of G-R relationship for each shallow zone.

Table 4a .
Probabilities of shallow zone for M L ≥ 6.0, M L ≥ 6.7, and M L ≥ 7.0 over the next 30 years.

Table 4b .
Probabilities of shallow zone for M L ≥ 6.0, M L ≥ 6.7, and M L ≥ 7.0 over the next 50 years.
Note: log: Log-normal model, gam: Gamma model, exp: Exponential model, wbl: Weibull model.The asterisks denote M u, c < M u, r .

Table 6a .
The high probability values are located at Eastern Taiwan, offshore Hualien County.The probability pattern from the gamma model is Probabilities of deep zone for M L ≥ 6.0, M L ≥ 6.7, and M L ≥ 7.0 over the next 30 years.

Table 6b .
Probabilities of deep zone for M L ≥ 6.0, M L ≥ 6.7, and M L ≥ 7.0 over the next 50 years.