Chiang Mai Journal of Science

Reconciling Conflicting Hazard Maps: A Methodological Reappraisal of Poisson and McGuire-Style PSHA Models in Urban Taiwan

1. INTRODUCTION

     Seismic hazard analysis plays a pivotal role in disaster prevention, seismic-resistant design, and risk-informed urban planning, particularly in tectonically active regions such as Taiwan. Positioned at the convergent boundary between the Philippine Sea Plate and the Eurasian Plate, Taiwan experiences frequent seismic activity. Significant events, such as the 1999 Chi-Chi earthquake (Mw 7.6) and the 1906 Meishan earthquake (Mw ~7.1), underscore the importance of robust seismic hazard frameworks in supporting engineering and public safety. The foundational theory of Probabilistic Seismic Hazard Assessment (PSHA) was introduced by Cornell [1], providing a quantitative method for estimating the likelihood of exceeding various levels of ground shaking at a specific site. This classical formulation treats earthquake occurrences as a Poisson process, which assumes time independence and a constant annual occurrence rate. The method has been widely adopted globally and underpinned the first PSHA maps of Taiwan [2]. Its mathematical simplicity and computational tractability have made it the dominant approach in seismic hazard zoning.
     However, mounting evidence has highlighted the limitations of the time-independent Poisson assumption in seismically complex regions. Earthquake clustering, fault interactions, and non-stationary seismicity can lead to hazard underestimation or temporal misrepresentation ([3-4]). For instance, Pailoplee and Charusiri utilized Poisson-based PSHA to evaluate seismic hazard in northern Thailand [5], while Triyoso et al. [6] extended this framework to assess tsunami scenarios in the Sumatran Subduction Zone. While effective for broad hazard zoning, these applications do not account for time-varying fault activity or inter-event dependencies.
     In response, more advanced modeling strategies have been developed. Notably, McGuire’s ([7-8]) extensions introduced time-dependent PSHA frameworks, capable of incorporating fault recurrence, stress interactions, and scenario-based simulations. This McGuire-style hazard modeling approach emphasizes forward-looking estimations, incorporating a continuum of possible events across magnitudes and distances, rather than relying solely on past discrete occurrences. Such frameworks have been adopted in landmark studies, such as [9] and the US National Hazard Map updates [10-11], and are increasingly advocated for in regions with complex fault systems. For example, Mosca et al. and Lindholm et al. employ the PSHA in modifying seismic hazard maps for the United Kingdom, mainland Norway, Svalbard, and the adjacent offshore regions [12-13].
     Recent studies, such as [14-15], have modernized Taiwan’s PSHA by integrating updated fault geometries, tectonic segmentation, and region-specific GMPEs. Lu and Chan invest in hazard maps around the Longitudinal Valley, combining pure characteristic earthquake and seismic hazard, as well as earthquake rates in fault systems. This fact reveals noticeable differences around the Milun Fault and the southern tip of the Longitudinal Valley [16]. However, these models often rely on a single probability framework. Recently, Xu et al. [17] have indicated that earthquake sequences in Taiwan exhibit non-random temporal patterns, with smaller events (Mw < 7.0) showing evidence of being influenced by the occurrence of larger magnitude earthquakes (Mw ≥ 7.0). The temporal occurrence of moderate-to-strong earthquakes (Mw < 7.0) is well characterized by the Gamma distribution, whereas the exponential distribution more accurately models large earthquakes (Mw≥ 7.0). This result challenges the applicability of a uniform Poisson model and highlights the need for a comparative evaluation using alternative probabilistic frameworks.
     Despite this, few studies have systematically compared how model choice—specifically, Poisson versus McGuire-style— affects hazard estimates across Taiwan. This gap is particularly significant for urban regions such as Taipei, Taichung, and Tainan, where PSHA results inform building codes and emergency response plans. Addressing this, the present study conducts a comparative sensitivity analysis using updated earthquake catalogs from 1973 to April 2025, fault recurrence intervals, and Taiwan-specific GMPEs. We evaluate the exceedance probabilities of Peak Ground Acceleration (PGA) under both models for six major cities and 28 seismic zones. Our framework includes seismic deaggregation to reveal dominant magnitude–distance contributions and examines model performance in regions with high seismicity or significant fault activity. Furthermore, we test a hybrid approach that applies BPT to known faults and Poisson to background sources, reflecting best practices in current international research [18-19].
     The contributions of this research are threefold:
1. It provides a critical methodological comparison between discrete historical modeling (Poisson) and forward-forecasting approaches (McGuire-style) within a single regional context.
2. It highlights regional hazard discrepancies, particularly in PGA zones exceeding 0.3g, where model outputs diverge significantly.
3. It proposes a pathway to integrate fault recurrence information using hybrid probability models, enhancing hazard realism without overcomplicating implementation. Ultimately, this comparative framework not only supports ongoing efforts to refine Taiwan’s national hazard model but also informs the global PSHA discourse on balancing simplicity, accuracy, and applicability across tectonic settings.
     This study presents one of the first comprehensive comparisons between Poisson-based and McGuire-style PSHA models applied to real-world seismic zones in Taiwan. The novelty lies not only in the side-by-side quantitative evaluation of model outputs—including hazard curves, deaggregation behavior, and spatial patterns—but also in the investigation of model suitability under different tectonic and recurrence conditions. Previous PSHA studies often apply a single model, but few question the implications of the underlying statistical assumptions. Our work fills this gap and provides a methodological reappraisal that can influence hazard forecasting frameworks globally. Furthermore, this study includes new hazard estimates reflecting the destructive earthquakes that occurred after 2019 in Taiwan, addressing the need for updated risk assessments.
     Our PSHA framework, although tested in Taiwan, is fully generalizable and adaptable to Southeast Asian settings such as Thailand, where limited historical or paleoseismic data are available. The inclusion of fault-specific recurrence models and probabilistic deaggregation makes this methodology a valuable tool for national hazard mapping and infrastructure planning across Southeast Asia. This result is particularly relevant following the April 2025 Myanmar earthquake, as the methods presented here can inform regional disaster preparedness and seismic risk mitigation in Thailand and neighboring countries.

2. MATERIALS AND METHODS

2.1 Methodological Framework
     This study employs a comparative approach to PSHA, evaluating both the classical Poisson-based model and the time-integrative McGuire-style model across Taiwan. The methodology consists of four key components: seismic catalog processing, source characterization, ground motion modeling, and hazard integration.

2.1.1 Dataset and zone zoning
     The seismic catalog used spans from 1973 to April 2025, compiled from the Central Weather Bureau Seismic Network (CWBSN) and other regional datasets. The occurrence model of earthquakes is assumed to follow the Poisson procedure when applying the PSHA approach by Cornell [1]. To remove aftershocks and foreshocks, we employed the declustering approach [20] to obtain the earthquakes recorded in the catalog, fitting the Poisson model. According to the suggestion offered by [21], the shallow regional sources are divided into 28 seismic zones. The Gutenberg-Richter recurrence relationship [22] with the relation between magnitude and frequency is described as

$$\log N(> ML) = a - bML                                           (1)$$

where $N(>ML)$ expresses the frequency of earthquakes with magnitudes larger than or equal to $ML$, $a$ is a constant related to the level of seismic activity, and $b$ is the slope of the power law. The $a$ and $b$ values of each seismic zone are estimated using least squares fitting. The zones are clearly classified into three groups (A, B, C) based on their $a$ and $b$ values, which represent different levels of seismic productivity and trends in magnitude distribution.

2.2 Methods
2.2.1 Fault source modeling
     In addition to background seismicity, Central Weather Bureau (CWB) integrated 38 active faults across Taiwan into the hazard model. The detailed information on these 38 active faults is listed in Appendix A1. For faults with known historical rupture dates (e.g., Chelungpu, Meishan, Milun), we applied the Brownian Passage Time (BPT) model [23-24] to account for time-dependent rupture probability. The time-dependent rupture probability within a forecasting window T can be estimated using the BPT model as below

$$P(t_0, T) = \Phi \left( \frac{t_0 + T - \mu}{\alpha\mu} \right) - \Phi \left( \frac{t_0 - \mu}{\alpha\mu} \right)                    (2)$$

where $\Phi$ denotes the standard normal cumulative distribution function, $\mu$ is the Mean recurrence interval of the earthquake (years), $t_0$ is the elapsed time since the last characteristic earthquake event (years), and $\alpha$ is the aperiodicity coefficient (i.e., the coefficient of variation of the recurrence interval), typically ranging from 0.3 to 0.6 for active crustal faults. For faults lacking temporal constraints, [25] adapted a Poisson model assuming a constant rupture rate. The Poisson model is represented as

$$P = 1 - e^{-\lambda(RT)}                                           (3)$$

     The recurrence intervals $RT$ were parsed from published literature [14]. In modeling the recurrence behavior of active faults, selecting between a time-independent (Poisson) and time-dependent (renewal-based) formulation has significant implications for seismic hazard forecasts. The Poisson model assumes memoryless behavior and constant annual rupture probability, suitable when recurrence intervals are poorly constrained or when uniform hazard mapping across many faults is required. By contrast, the BPT model incorporates temporal memory and aperiodicity into the rupture cycle, making it particularly relevant for faults with well-documented historical ruptures. The BPT model, by incorporating aperiodicity, reflects more realistic rupture dynamics and is thus preferable for detailed risk assessments of critical infrastructure. In this research suggests that for faults with known last-event dates and recurrence intervals, the BPT model yields significantly higher conditional rupture probabilities when the elapsed time approaches or exceeds the mean recurrence interval. In contrast, the Poisson model underestimates near-term probabilities but may be more stable when data is limited. Selecting an appropriate fault rupture model is crucial for reliable seismic hazard assessment. The choice depends on the availability of paleoseismic records, the required spatial resolution, and the criticality of the assets being assessed. Based on our findings and current best practices, we recommend the following:
     First, the Poisson model is suitable when the rupture history of a fault is unknown, incomplete, or unreliable. Its simplicity and computational efficiency make it particularly useful for regional hazard mapping where rapid processing is necessary. It also serves as a stable baseline model when data constraints prevent the use of time-dependent approaches. Second, the BPT model is preferred when paleoseismic or historical earthquake data provide a well-constrained mean recurrence interval and a known time since the last rupture. The BPT model accounts for the temporal clustering of events and reflects the increasing conditional probability as time since the last event grows, which the Poisson model cannot capture. Third, for critical infrastructures located near active faults such as nuclear power plants, high-speed railways, or urban centers, the BPT model is strongly recommended. It enables a more realistic time-dependent probability estimate, essential for assessing near-term risk under long-overdue seismic conditions. Ultimately, in the context of seismic zonation with limited data resolution, the Poisson model remains a suitable approximation. It enables the incorporation of regionalized Gutenberg–Richter parameters while maintaining manageable computational complexity.
     These guidelines support a hybrid modeling strategy, wherein Poisson-based estimates serve broad-scale planning, and BPT-based time-dependent models are used for site-specific or high-priority infrastructure risk analysis.

2.2.2 Ground motion prediction equations (GMPEs)
     We employed the Taiwan-specific crustal GMPEs [26-27] to estimate median peak ground acceleration (PGA) as a function of magnitude, distance, site condition ($VS_{30}$), and fault mechanism. A standard deviation ($\sigma = 0.4$) was assumed for logarithmic PGA distribution. Site effects were incorporated using local $VS_{30}$ maps from [28], while distances were computed using the haversine formula and geometric fault centers.

2.2.3 Hazard computation
     For the Poisson model, annual exceedance probabilities were computed per magnitude-distance bin and then aggregated. In the McGuire-style model, a Monte Carlo integration over magnitude and distance distributions was employed, convolving recurrence rates, GMPE-based exceedance probabilities, and distance PDFs. We constructed hazard curves for each city by combining background and fault contributions, with exceedance probabilities computed for PGA values ranging from 0.01 g to 1.0 g over a 30-year exposure period. Hazard deaggregation was also performed for selected cities to identify dominant magnitude–distance contributions to risk, offering further insight into regional hazard drivers under each model.

2.2.4 Poisson-based PSHA implementation
     The Poisson-based model follows the classical Cornell [1] framework, where earthquakes are treated as statistically independent events in time. For each recorded earthquake $i$, the annual occurrence rate is computed as:

$$\lambda_i = \frac{10^{a-bM_i}}{T}                               (4)$$

where $a$ and $b$ are the Gutenberg–Richter parameters for the assigned seismic zone, $M_i$ is the magnitude of the event, and $T$ is the catalog duration. The probability that the PGA at a site exceeds a threshold $y^*$ is given by:

$$P_i = 1 - \Phi \left( \frac{\log_{10}(y^*) - \mu(M_i, R_i)}{\sigma} \right)                                  (5)$$

where $\mu$ and $\sigma$ are the mean and standard deviation from the GMPE, and $R_i$ is the source-to-site distance. The total background exceedance probability during period time $\Delta t$ (year) which the exposure time window converts annual rates into cumulative risk is calculated via the exponential probability of non-exceedance:

$$P_{bg} = 1 - \exp(-\sum_i \lambda_i \cdot P_i \cdot \Delta t)                        (6)$$

     Fault-based seismic hazard is calculated independently based on the rupture probability of each fault using a hybrid time-dependent model and then integrated into the total hazard via:

$$P_{total} = 1 - (1 - P_{bg})(1 - P_{fault})                                                  (7)$$

where the term $P_{fault}$ means the probabilities of faults rupture described in section 2.2.

2.2.5 McGuire-style PSHA implementation
      The McGuire-style model generalizes hazard by integrating over a continuum of magnitudes $M$ and distances $R$, bypassing the limitations of catalog dependence. The resulting value $P_Y(y^*)$ is the total probability that the ground motion exceeds the threshold $y^*$ at site $Y$, typically within a specified time frame ([7-8]):

\[ P_Y(y^*) = \int_{M_{\min}}^{M_{\max}} \int_{R_{\min}}^{R_{\max}} M\, f_M(M)\, f_R(r)\, \left[ 1 - \Phi\!\left(\frac{\log_{10} y^* - \mu(M,R)}{\sigma}\right) \right] \, dr\, dM \tag{8} \]

where:
• M_min and M_max denote the minimum and maximum considered earthquake magnitudes, respectively, which define the magnitude range of interest for hazard analysis.
• R_min and R_max represent the minimum and maximum source-to-site distances, defining the spatial range within seismic sources can contribute to the ground motion hazard at the site.
• f_M(M) is an unnormalized magnitude-frequency function based on the Gutenberg–Richter relation, describing the relative contribution of different magnitudes to seismic activity. The commonly used expression is \( f_M(M)=b\ln(10)\cdot10^{-bM} \). This functional form is commonly seen in PSHA-related hazard integration, and arises from the differentiation of the Gutenberg–Richter cumulative law combined with a base-10 to natural logarithm conversion [8]. From a mathematical perspective, within the hazard integral, this function acts not as a conventional PDF but rather as a seismicity weighting function, scaling the contribution of each magnitude to the overall exceedance probability.
• f_R(r) is a distance weighting function derived from the spatial geometry of the seismic source zones, representing the relative likelihood or contribution of rupture distances from the source region to the target site. While often treated analogously to a probability density function, it is generally used as a weighting term in the PSHA hazard integral and may not be strictly normalized. The distance-weighting function \( f_R(r) \), which determines the geometric contribution of each source zone is obtained numerically using a Monte Carlo sampling scheme. Within each seismic-source polygon D, N random points are uniformly generated, and their great-circle (Haversine) distances r to the site are computed. The distance range is discretized into bins [r_k, r_k+1] of width Δr_k, and the number of sampled points falling in each bin n_k is counted. The resulting Monte Carlo estimator

\( \hat f_R(r_k)=\dfrac{n_k}{N\Delta r_k} \)         (9)

represents the empirical probability density of source-to-site distances, where n_k is the number of randomly sampled points within the k^th distance bin of width Δr_k, and r_k is the bin midpoint. As N → ∞ and Δr_k → 0, \( \hat f_R(r) \) converges to the true f_R(r), which depends explicitly on the geometry and relative position of the seismic-source zone with respect to the site. The integral of f_R(r) over all distances equals unity:

\( \int_{0}^{R_{\max}} f_R(r)\,dr = 1 \)                  (10)

     The function f_R(r) denotes the probability density of source–to–site distances for all potential hypocenters uniformly distributed over a seismic source polygon D (area A(D)). Because it measures how much of the site-centered iso-distance curve at radius r lies inside D, f_R(r) is inherently geometry-dependent, reflecting the shape, area, orientation, and relative position of the site. Let L_D(r) be the length of the intersection between D and the circle of radius r centered at the site, in a planar approximation, f_R(r) can be written as:

\( f_R(r)=\dfrac{L_D(r)}{A(D)} \)                            (11)

Alternatively, if Θ(r) ⊂ [0,2π] is the set of polar angles for which the circle of radius r lies inside D, then

\( f_R(r)=\dfrac{r}{A(D)}\times|\Theta(r)| \)               (12)

     These formulas make explicit that a longer intersection (i.e., larger L_D(r) or |Θ(r)|) yields higher density at that distance, indicating that f_R(r) varies systematically with the geometry of D. According to equation (9), different geometries of D alter the counts {n_k} in a consistent manner—producing narrower or broader peaks for compact or circular zones, sharp shoulders or skewed distributions for offshore or remote polygons. Consequently, wide zones produce broader f_R(r) curves, elongated or asymmetric zones shift the distribution toward preferred distances, and sites located near the zone margin yield highly skewed f_R(r) functions. This formulation ensures that complex source shapes and spatial configurations are faithfully incorporated into the hazard integral of Equation (8), while the Monte Carlo approach preserves all geometric effects without requiring analytical simplifications.
• μ(M,R) is the median predicted logarithmic ground motion for a given magnitude M and distance R, calculated using a GMPE.
• σ is the standard deviation of the logarithmic ground motion, reflecting aleatory uncertainty (variability) in ground motion prediction.
• Φ(·) is the cumulative distribution function (CDF) of the standard normal distribution, used to compute the probability that the ground motion is less than y*. Subtracting this value from 1 gives the exceedance probability.
     The annual exceedance rate v(y*) is integrated over all M–R combinations to yield:

\( v(y^*)=\int_M\int_R \lambda(M)f_R(R)\left[1-\Phi\!\left(\dfrac{\log_{10}y^*-\mu(M,R)}{\sigma}\right)\right]\,dR\,dM \)                   (13)

The total hazard is then:

\( P_{bg}=1-\exp(-v(y^*)\cdot\Delta t)=1-\exp[-\lambda_{tot}P_y(y^*)\Delta t] \)                    (14)

where λ_tot is the total annual occurrence rate which is defined as

\( \lambda_{tot}=\int_{M_{min}}^{M_{max}} \lambda(M)\,dM \)                          (15)

     Accordingly, Equations (9) and (10) guarantee that complex source geometries and spatial configurations are rigorously incorporated into the hazard integral of Equation (8). In this formulation, P(y*) in Equation (8) constitutes a purely spatial–magnitude integration, whereas the temporal dependence of the site-specific exceedance probability is governed by the external rate term v(y*) and the exposure time Δt in Equation (14). Fault hazard is incorporated using the same rupture probability and GMPE formulations as in the Poisson or BPT models. Using Eq. (8), total probabilities of the McGuire-style model can be derived. The combined model framework enables both broad-scale hazard planning, utilizing Poisson assumptions, and site-specific risk assessments that incorporate BPT-based time dependence. This hybrid approach improves hazard realism and supports resilient infrastructure planning under varying data constraints.

3. RESULTS AND DISCUSSION

3.1 Results
3.1.1 Seismicity model validation and zonal clustering analysis
     This section provides a comprehensive validation of the seismicity model used in this study, followed by a clustering analysis of regional seismic behavior. In this study, local magnitude from the CWBSN catalog is referred to as ML, while moment magnitude derived from fault parameters or major earthquakes is called Mw. This distinction aligns with international seismological standards [14]. All figures, tables, and text have been updated accordingly. The analysis compares the Poisson and McGuire-style PSHA results across Taiwan, focusing on six urban centers: Taipei, Taichung, Tainan, Kaohsiung, Hualien, and Taitung. Since Taiwan lies at the junction of the Philippine Sea Plate and the Eurasian Plate, it experiences high seismicity. Figure 1(a) shows the distribution of earthquake hypocenters from January 1973 to April 2025.
     To meet the independence assumption of Poisson-based PSHA, we applied a space–time windowing algorithm to decluster the catalog [20]. Figure 1(b) compares the frequency–magnitude histograms before and after declustering. After declustering, the total number of events decreases by about two-thirds across all magnitudes, with proportionally greater reductions in highly active zones where clustering is most intense. Only ~35% of the original events were retained, yielding a more suitable dataset for hazard modeling. Declustering improved the estimation of the completeness magnitude (Mc). Using the maximum curvature method, a stable Mc ≈ 2.5 was determined across most zones and applied as the lower bound in Gutenberg–Richter (G–R) recurrence model fitting. This approach prevents artificial inflation of b-values due to aftershocks, preserving the physical reliability of long-term seismicity trends.
     Figure 1(c) presents the regional zonation framework, consisting of 28 seismic zones (S01–S21) [21], overlaid with 38 active faults [14]. The locations of six major cities are also marked, highlighting the intersection between urban regions and seismogenic zones, forming the foundation for hybrid PSHA modeling. Because no higher-level groupings are defined, we introduce a visual classification for readability: zones are organized into six groups based on tectonic province and geographic continuity, following a centroid-in-province rule (ties resolved by maximum overlap). This grouping is used solely for figure presentation and does not influence the hazard calculations. Figure 1(d) illustrates the administrative boundaries of Taiwan's counties, serving as a spatial reference framework for the subsequent analysis of seismic probability distributions.
     G–R relationships were fitted for each zone, with varying a and b-values reflecting differences in earthquake productivity and magnitude scaling. For instance:
1. Zone S04 exhibits a high b-value (1.03), indicative of mature or creeping faults.
2. Zone S06 exhibits a low b-value (0.77) and nonlinear deviation for ML > 5.5, indicating the dominance of significant thrust events.
3. The offshore zone S15 yields the highest a-value (6.61), reflecting its high productivity.
4. Zone S21 has a low cumulative event count and a b-value of 0.79, consistent with a low-seismicity regime.
     To quantify uncertainties, bootstrap resampling was used to estimate confidence intervals for b-values across all zones (Figure 2(a)). Most zones fall within the 0.8–1.1 range, consistent with global observations [28]. Outliers such as S03, S13, and S06 exhibit lower b-values with narrow bounds, indicating a higher likelihood of significant events. Conversely, S04, S05A, and S14B show higher b-values, indicating predominance of small events. Figure 2(b) plots the relationship between a and b-values. While no strong linear correlation exists, distinct seismic patterns emerge, justifying a zonal rather than uniform recurrence model in PSHA. Based on a–b characteristics, the zones are preliminarily classified into the following clusters:
     Group A (High a, b ≈ 1): e.g., S15, S14B, S5A
• High microseismic productivity with balanced magnitude scaling. It likely associated with active compressional stress or shallow fault rupture environments. Group B (Moderate a, Low b): e.g., S06, S08A, S03
• Characterized by moderate activity and a larger proportion of strong events. These zones warrant long-term monitoring due to their capacity for infrequent but destructive earthquakes. Group C (Low a, b ≈ 0.9): e.g., S13, S20, S01
• Likely represent stable or recently relaxed regions with minimal seismicity, serving as comparative baselines.
• Notable outlier behavior further enriches this framework:
• S5A (High a, High b): High small-event frequency, possibly due to geothermal or hydrothermal systems. • S03 (Low a, Low b): May indicate a large locked fault segment with long recurrence intervals and high rupture potential.

     In summary, this zonation and clustering analysis not only validates the seismic model through empirical a–b relationships and Mc estimation, but also provides interpretive insight for regional PSHA refinement. It emphasizes the need to distinguish between background and fault-based sources in Taiwan's tectonically complex setting.

3.1.2 Comparisons of exceeding probabilities (McGuire vs. Poisson Models)
     The uses of both Poisson-based and McGuire-style PSHA frameworks are grounded in well-established seismic hazard literature ([1], [8-9]), ensuring theoretical rigor and methodological validity. Here, we initially generated a series of exceedance probability maps under the McGuire-style and Poisson PSHA framework for four PGA thresholds: 0.10g, 0.30g, 0.50g, and 0.70g. The results are displayed in Figure 3(a) and (b), respectively, to visualize the spatial variation in seismic hazard across Taiwan. These maps illustrate the probability that ground shaking at a site will exceed the specified PGA value within 30 years. The ID, names, and corresponding data of the faults are listed in Table 1. These models allow for a systematic comparison of hazard intensity and regional vulnerability.

    Figure 3(a) displays the McGuire-based maps, revealing pronounced regional disparities in the probability of ground motion exceedance. Four key patterns emerge:
1. Concentrated Hazard in Central and Eastern Taiwan Regions along the Central Mountain Range and Longitudinal Valley, such as near the Longitudinal Valley Fault (ID 33), maintain moderate exceedance probabilities (>0.3) even at 0.70g, indicating sustained fault activity and elevated seismic risk. These zones demand priority in infrastructure strengthening and urban resilience planning.
2. Sharp Hazard Decline in Western Plains While the western plains (e.g., Yunlin, Chiayi, Tainan) show high hazard (>0.9) at 0.10g, their exceedance probabilities drop sharply below 0.1 at 0.50g, suggesting more stable crustal conditions and reduced potential for high-intensity shaking.
3. Regional Stratification and Planning Implications Southern areas like Tainan, Kaohsiung, and Pingtung retain notable hazard levels (e.g., ~0.3 at ≥0.30g, ~0.1 at ≥0.50g), with Pingtung exceeding 0.1 even at 0.70g. In contrast, central-western regions (Nantou, Taichung, Changhua Plain) exhibit rapid attenuation (<0.05 at 0.50g), underscoring the need for regional hazard-based planning and tiered seismic design.
4. Strong Correlation with Active Faults High-hazard zones align closely with active faults, such as the Milun (ID 32), Hengchun (ID 30), Luyeh (ID 35), and Chaochou (ID 29) faults. These faults exhibit high slip rates (>6 mm/yr) and short recurrence intervals (66–400 years), reinforcing the persistent threat posed by fault systems. Central Taiwan faults (e.g., Chelungpu, Meishan) contribute less to long-term hazard due to longer recurrence intervals and more localized influence.

     These spatial insights emphasize the McGuire model's utility for integrating fault-specific recurrence into resilience planning. The resulting maps support data-informed design thresholds for critical facilities and targeted public safety initiatives. Further integration with population density and infrastructure vulnerability layers can enhance disaster risk mapping.
     Figure 3(b) presents the Poisson-based exceedance maps, which reflect a temporally independent hazard structure. Although simpler, several key features emerge:
1. Rapid Hazard Attenuation with Increasing PGA As the PGA threshold rises, exceedance probabilities fall below 0.3 across most regions. The model's assumption of time independence limits its ability to account for accumulated tectonic stress, making it conservative in forecasting strong shaking in previously quiet areas.
2. Localized Peaks in Southern Taiwan At PGA ≥ 0.5g, southern Taiwan retains several hazard peaks aligned with faults such as the Chiayi Frontal (ID 21) and Chaochou (ID 29), indicating the model's sensitivity to spatial seismicity rates despite limited recurrence modeling.
3. 0.30g as a Zoning Benchmark The 0.30g contour effectively separates moderate-hazard zones in the east and south from low-hazard areas in the west and north. This result establishes a practical threshold for seismic zoning and prioritization of structural retrofitting and policy interventions.
4. Consistently Low Hazard in Northern Taiwan Taipei, New Taipei, and Keelung exhibit low hazard across all PGA thresholds, consistent with their lower fault density and historically moderate seismicity. These findings underscore the importance of incorporating geologic context when interpreting hazard results.
     The McGuire-style model, which incorporates time-dependent fault activity, offers enhanced hazard resolution in seismically active zones and is better suited for performance-based design and critical infrastructure planning. In contrast, the Poisson model provides a conservative, smoothed hazard landscape, which is helpful for generalized zoning and public education. A hybrid strategy that leverages the strengths of both models is recommended to optimize seismic resilience across Taiwan.

3.1.3 Comparisons of uniform-probabilities ground-motion maps
     Uniform-probability ground motion maps (Figure 4) offer insight into the expected PGA values at fixed exceedance probabilities (5%, 10%, 30%, and 50%) over 30 years. These maps, generated using both McGuire-style and Poisson-based PSHA frameworks, reveal critical differences in hazard magnitude, spatial pattern, and policy implications. The McGuire-style model (Figure 4(a)), which incorporates fault-specific recurrence intervals and slip rates, shows progressively decreasing PGA values with increasing exceedance probability. At the 5% level, central and southern Taiwan experience the highest PGA values (>0.6–0.8g), especially near major active faults such as Changhua (ID 16), Chaochou (ID 29), Hengchun (ID 30), and the Longitudinal Valley Fault (ID 33). These areas align with high slip rates and short recurrence intervals, indicating the potential for rare but catastrophic shaking.
     At 10–30% probabilities, PGA values decline to ~0.3–0.4 g, with remaining hotspots located near the Chaochou Fault and the Ilan region (ID 38)—these persistent zones of elevated hazard support regional prioritization for seismic retrofitting and emergency preparedness. By the 50% level, PGA values fall below 0.2g across most of Taiwan, defining a baseline for routine infrastructure design. This gradient illustrates the inverse relationship between ground motion intensity and exceedance probability in PSHA, confirming the McGuire model’s effectiveness in capturing both rare, high-impact and frequent, low-impact shaking scenarios.

    The Poisson model (Figure 4(b)) assumes temporally independent seismicity and produces smoother hazard maps. The Poisson maps are spatially smoother, quantified by the variance of the exceedance probability over land (σ²). For the same PGA threshold, σ²(Poisson) ≈ 0.05 vs σ²(McGuire) ≈ 0.22; smaller σ² indicates lower spatial variability, confirming a markedly more heterogeneous pattern for the McGuire-style maps. While it replicates the inverse PGA-probability relationship, the magnitude and spatial concentration of the hazard are lower than those of the McGuire model, particularly at high-intensity thresholds. PGA values at 5% and 10% probabilities peak at ~0.45–0.6g in eastern and southern Taiwan but attenuate rapidly across the island.
     Despite these limitations, the Poisson model effectively captures key high-hazard areas, such as the Longitudinal Valley and the Tainan–Pingtung region, supporting its use in general hazard stratification and public education. However, its conservative hazard estimates at high PGA levels suggest potential underestimation in fault-dense zones due to the lack of time-dependent recurrence input.
     Table 2 summarizes the comparative performance. The McGuire-style model demonstrates greater resolution in fault-dense areas, especially at higher PGA thresholds. It is more suitable for critical infrastructure design—such as hospitals, transportation hubs, and emergency centers—where preparedness for rare events is crucial. Conversely, the Poisson model provides a generalized hazard landscape, ideal for large-scale zoning or planning in areas with limited paleoseismic data.
     These findings align with the results of [9] and [29], leading to the following recommendations:
1. Dual-Model Strategy: Use McGuire-style PSHA for site-specific engineering in high-risk zones and Poisson-based mapping for baseline hazard and public awareness.
2. Hazard Tiering for Seismic Zoning: Classify PGA levels into low (<0.2g), moderate (0.2–0.5g), and high (>0.5g) hazard zones. Eastern and southern Taiwan require higher code standards due to persistent hazards at all thresholds.
3. Targeted Investment: Allocate resources to regions that show consistently high PGA values across models, including Hualien, Taitung, Chiayi, and Nantou.
4. Integrated Risk Mapping: Future research should overlay hazard data with population, infrastructure, and soil amplification layers to support multi-hazard mitigation planning.

     In summary, the uniform-probability ground motion maps reveal a clear spatial stratification of seismic risk in Taiwan. The McGuire model’s time-dependent framework better identifies persistent hazard zones, while the Poisson model offers a simplified alternative for regional-scale risk communication. A hybrid approach strengthens resilience planning and supports sustainable seismic risk governance. To complement the spatial and probabilistic comparisons presented above, we further performed a comprehensive statistical validation to examine whether the McGuire-style framework better represents the temporal behavior of historical seismicity in Taiwan than the classical Poisson model. We analyzed inter-event times from 1973–2025 in 28 zones (using both original and declustered catalogs) by fitting Exponential (Poisson surrogate), Gamma, and Lognormal models via maximum likelihood and comparing them using the Akaike Information Criterion (AIC = 2k − 2 log L). The Exponential model was never selected as best in any zone; after declustering, Lognormal dominated (≈25/28), with Gamma preferred in a few cases. Whole-island analyses likewise favored Lognormal by decisive ΔAIC margins. The Kolmogorov–Smirnov (KS) statistic (D) provided consistent quantitative evidence. In the original catalog, KS D values averaged about 0.58 for the Gamma model and ≈ 0.07 for the Lognormal model, whereas after declustering, they decreased sharply to ≈ 0.16 and ≈ 0.06, respectively. The Poisson (Exponential) model showed a similar reduction from roughly 0.5–0.6 before declustering to about 0.08–0.10 afterward. However, this apparent improvement mainly reflects the statistical randomization effect of declustering rather than a genuine enhancement of model performance. Indeed, AIC comparisons show that the Poisson model remained decisively inferior to both the Gamma and Lognormal distributions, and it was never the best-fitting distribution in any of the zones. Dispersion diagnostics also reject the Poisson assumption: coefficients of variation (CV) and annual Fano factors were greater than 1 in most zones (original) and remained greater than one after declustering, indicating temporal clustering and over-dispersion. Taken together, these independent indicators (AIC ranking, KS D, CV/Fano) demonstrate that seismicity in Taiwan is non-Poissonian and is better captured by McGuire-style, time-dependent models—particularly the Lognormal distribution—consistent with stress diffusion, memory effects, and rate evolution. Overall, these combined statistical indicators demonstrate that the McGuire-style, time-dependent model offers a significantly better fit to Taiwan's historical seismicity than the classical Poisson process. Nonetheless, analyses incorporating more extended temporal coverage or regionally refined datasets could further consolidate this conclusion.

3.1.4 Comparisons of seismic hazard curves in six major urban cities
     Taiwan’s six major metropolitan areas—Taipei, New Taipei, Taoyuan, Taichung, Tainan, and Kaohsiung—are highly populated, economically vital, and seismically vulnerable. Accurate assessment of probabilistic seismic hazard in these cities is critical for disaster mitigation, infrastructure resilience, and public safety.
     We compare Seismic hazard curves derived from the McGuire-style with Poisson-based PSHA models across the six cities (Figure 5(a–b)). Instead of relying on peak values from single nodes, the curves were based on spatially averaged exceedance probabilities within each city’s boundary, smoothing out local anomalies caused by site effects or nearby faults.
     The McGuire-style curves exhibit smoother decay with increasing PGA, reflecting a wide integration of plausible earthquake scenarios and time-dependent recurrence along known faults. In contrast, Poisson-based curves are jagged and reflect sensitivity to the limited number of historical events in the declustered catalog (1973–2025). The divergence is especially evident at PGA ≥ 0.2g: in Tainan and Kaohsiung cities, McGuire-style models yield exceedance probabilities of 0.4–0.5, while Poisson-based values remain around 0.1–0.2 due to the absence of recent significant events in the historical record. Even in northern cities like Taipei and New Taipei, the McGuire model predicts higher seismic risk by incorporating regional subduction events. These differences highlight the McGuire model’s advantage in capturing epistemic uncertainty and rare-event potential, providing a more conservative and comprehensive hazard outlook that is essential for infrastructure design and insurance modeling. Table 3 summarizes the significant discrepancies of the two Seismic hazard curves derived from the two models.
     Deaggregation analyses further clarify the dominant seismic sources contributing to the hazard in each city (Figure 5(c–d)). Under the McGuire-style PSHA, exceedances at design-level ground motions (~0.3g) are mainly driven by moderate-to-large magnitude earthquakes (Mw 6.0–6.8) at near-field distances (13–33 km). For example, the cities of Tainan and Kaohsiung exhibit sharply localized hazard contributions from the Hsinhua Fault (ID 24) and Chaochou Fault (ID 29). At the same time, Taichung displays more diffuse contributions from Chelungpu (ID 17) and Tamaopu–Shuangtung (ID 18), reflecting a more complex fault interaction.
     In contrast, Poisson-based deaggregation attributes hazard to moderate events (Mw 6.0–6.5) at broader distance ranges. While it captures contributions from known faults such as Shanchiao (ID 1) in the north and Chelungpu in central Taiwan, the hazard is more diffusely distributed and closely tied to past cataloged events. For cities like Taoyuan and Tainan, the spatial dispersion of sources may stem from incomplete records or unresolved background sources. A fundamental limitation of the Poisson model is its inability to account for unrecorded or infrequent significant events due to its time-independent formulation. As a result, it tends to underestimate tail risks, especially in tectonically complex regions such as Taiwan, where clustered faults and short observation windows prevail. By contrast, the McGuire-style model incorporates fault recurrence, magnitude–frequency distributions, and geometries, yielding sharper and more realistic Mw–R peaks.

     Table 4 summarizes key differences. The McGuire-style model offers a physically meaningful decomposition of hazard, identifies critical magnitude–distance pairs, and is thus highly suited for scenario planning, early warning, and risk-informed infrastructure design. While the Poisson model adequately captures historical seismic trends, it may underestimate hazard levels in regions characterized by sparse data or high uncertainty. To assist readers, the principal abbreviations and notations used in this study are compiled in Table 5, providing the full terms and concise explanations of all acronyms relevant to the PSHA framework.
     These findings reaffirm the importance of selecting appropriate PSHA models based on the intended application. For performance-based design and high-risk urban planning, the McGuire-style model provides superior insight. For general zoning and public education, the Poisson model remains a valuable tool for analysis and understanding. Together, these models support a dual-framework strategy to enhance seismic safety and resilience across Taiwan’s urban landscape.

3.2 Discussion
     This study presents a comprehensive comparison between Poisson-based and McGuire-style PSHA models across Taiwan, highlighting significant differences in peak ground acceleration (PGA) estimates, spatial hazard patterns, and source attributions. The McGuire-style model, which incorporates fault-specific recurrence intervals, slip rates, and magnitude–distance distributions, consistently yields higher PGA values—often by 20–40%—particularly in tectonically active regions such as eastern and southern Taiwan. These elevated estimates are especially pronounced in low-probability, high-consequence scenarios, underscoring the importance of accounting for significant, infrequent events in hazard planning. In contrast, the Poisson model, constrained by historical earthquake catalogs and the assumption of time invariance, tends to underestimate hazard in regions characterized by long recurrence intervals or limited seismic records. Spatial analysis reveals that McGuire-style outputs align closely with known fault structures, producing hazard contours that accurately reflect geological reality.
     In contrast, Poisson-based models diffuse risk more uniformly, potentially obscuring site-specific threats. Deaggregation analyses further highlight the superiority of the McGuire model in identifying near-field fault contributions—information critical for seismic design and emergency response planning. These differences have direct implications for seismic hazard zoning and land-use regulation. For instance, cities like Tainan, Hualien, and Taitung exhibit exceedance probabilities greater than 0.5 for 0.3g and 0.5g PGA thresholds in the McGuire model, while these same zones appear muted or diffuse in Poisson-based maps. Such discrepancies could lead to an underestimation of hazard exposure if official zoning is based solely on historical catalogs.

     Furthermore, the source-specific resolution of the McGuire model enables more precise identification of high-risk corridors, such as the Longitudinal Valley and Hengchun Peninsula (ID 30-31) supporting adaptive infrastructure planning and prioritization of seismic retrofitting. Although the McGuire-style model has been critiqued for its complexity or presumed neglect of historical data, these concerns are unfounded. It integrates catalog information via statistically calibrated recurrence parameters, and its outputs—hazard curves and deaggregation plots—remain interpretable and suitable for regulatory and engineering applications. The probabilistic framework also explicitly addresses epistemic uncertainty by combining fault-based and catalog-based data, thereby enhancing rather than undermining resilience-oriented planning. Based on these insights, we advocate for an integrated PSHA framework that leverages both models: the Poisson approach for regional baseline estimation and data-limited areas, and the McGuire model for fault-specific refinement in high-risk zones. A proposed “hazard sensitivity” metric—quantifying the variability of hazard estimates across model types—further adds diagnostic value for policy and academic use. Ultimately, such a hybrid approach enhances the precision, interpretability, and decision relevance of seismic hazard assessments, ensuring that both historical seismicity and future fault behavior are incorporated into risk-sensitive planning for Taiwan’s seismic resilience.

4. CONCLUSIONS

     This study presents a comprehensive comparison between McGuire-style and Poisson-based PSHA models, applied across a broad range of urban centers and seismically active regions throughout Taiwan. Rather than limiting the scope to only six major metropolitan areas, the analysis incorporates a wide grid-based assessment that captures both high-population zones and areas in proximity to active faults. By integrating deaggregation results, hazard curve comparisons, and PGA-based hazard maps, the study demonstrates that the McGuire-style model provides a more robust, physically consistent, and forward-looking evaluation of seismic hazard. Notably, the McGuire-style PSHA reveals localized zones of elevated hazard that closely align with known fault structures, particularly in eastern and southern Taiwan. Regions such as Tainan, Kaohsiung, Hualien, and Taitung frequently exhibit exceedance probabilities greater than 0.5 for PGA thresholds of 0.3g and 0.5g—thresholds that are highly relevant for seismic design.
     In contrast, the Poisson-based model tends to show rapid hazard attenuation with increasing PGA levels and weaker spatial correlation with fault traces. Its reliance on historical catalogs and time-independent assumptions may result in an underestimation of hazard in tectonically complex regions. Hazard curves reinforce this divergence: McGuire-style curves are smoother and integrate magnitude–distance distributions with fault recurrence statistics, while Poisson-based curves remain jagged and catalog-driven.
     Deaggregation analyses further confirm these differences. The McGuire-style model attributes dominant contributions to large-magnitude (Mw 6.6–6.8), near-field events occurring within 13–38 km of a given site, often linked to specific fault systems. In contrast, Poisson-based deaggregation highlights smaller, more distant events, reflecting a retrospective, less fault-sensitive perspective on hazard. From an engineering and disaster policy standpoint, these findings underscore the value of adopting McGuire-style PSHA frameworks in regions like Taiwan, where fault density, recurrence potential, and urban exposure are high, yet long-term seismic records remain incomplete. The McGuire-style model offers significant advantages for performance-based design, critical infrastructure planning, and evidence-based resilience strategies. While Poisson-based PSHA can still serve as a conservative baseline in data-sparse areas, it lacks the predictive resolution required for modern seismic safety planning.
     This research directly contributes to the goals of Sustainable Development Goals 11 (Sustainable Cities and Communities) and 13 (Climate Action) by supporting risk-informed urban development, refining seismic codes, and promoting forward-looking adaptation to earthquake hazards. It also aligns with the Sendai Framework for Disaster Risk Reduction, particularly in advancing the understanding of disaster risk and strengthening risk governance. The comparative findings of this study not only highlight key modeling differences between Poisson-based and McGuire-style PSHA frameworks but also underscore the importance of region-specific hazard models. As recent seismic swarms in Southeast Asia, such as the 2012 Phuket event [30], have demonstrated the presence of active crustal processes beyond primary fault zones, this study’s comparative PSHA framework offers a timely reference model for regional seismic resilience planning in Thailand and neighboring regions. In light of recent seismic events in Taiwan and Myanmar, this work provides a valuable foundation for extending hazard assessments to Thailand and other countries in Southeast Asia, thereby supporting data-informed disaster risk reduction policies.

ACKNOWLEDGEMENTS

     This research was supported by funding from Nantong Institute of Technology. The authors also thank the Central Weather Bureau and the Central Geological Survey of Taiwan for providing earthquake catalog and active fault data. We are grateful for the helpful discussions with colleagues and the valuable suggestions from anonymous reviewers that improved the quality of this manuscript.

     This research was supported by Nantong Institute of Technology. Dr. Liao was supported by the Natural Science Foundation Project of Nantong City (2023) at Nantong Institute of Technology under Contract No. JC2023109.

AUTHOR CONTRIBUTIONS

Boi-Yee Liao: Conceptualization, Methodology, Resources, Software, Validation, Formal analysis, Data curation, Writing-Original draft preparation, Writing-Reviewing and Editing. Xu Wu: Supervision, Funding acquisition, Project administration.

CONFLICT OF INTEREST STATEMENT

     The authors declare that they hold no competing interests.

DECLARATION OF GENERATIVE AI IN PREPARATION OF MANUSCRIPT

     During the preparation of this work, the authors used ChatGPT (GPT-4, OpenAI) to improve the readability and language of the manuscript. After using this tool, the authors reviewed and edited the content as necessary and take full responsibility for the publication's content.

FUNDING

     This research was financially supported by Nantong Institute of Technology Fund (Grant Number JC2023109)

REFERENCES

[1] Cornell C.A., Engineering seismic risk analysis. Bulletin of the Seismological Society of America, 1968; 58(5): 1583–1606. DOI 10.1785/BSSA0580051583.

[2] Cheng C.T., Chiou S.J., Lee C.T. and Tsai Y.B., Study on probabilistic seismic hazard maps of Taiwan after Chi-Chi earthquake. Journal of GeoEngineering, 2007; 2: 19-28. DOI 10.6310/jog.2007.2(1).3.

[3] Malakar S., Rai A.K., Kannaujiya V.K. and Gupta A.K., Probability and recurrence interval of large earthquakes (Mw $\ge$ 6) in the Himalayan seismicity zone. Geomatics, Natural Hazards and Risk, 2024; 15(1): 2393667. DOI 10.1080/19475705.2024.2393667.

[4] Zhuang J., Ogata Y. and Vere-Jones D., Analyzing earthquake clustering features by using stochastic reconstruction. Journal of Geophysical Research: Solid Earth, 2004; 109(B05): B05301. DOI 10.1029/2003JB002879.

[5] Pailoplee S. and Charusiri P., Probabilistic analysis of the seismic activity and hazard in northern Thailand. Geosciences Journal, 2015; 19: 731–740. DOI 10.1007/s12303-015-0011-8.

[6] Triyoso W., Kongko W. and Prasetya G.S., Understanding earthquake potential for future hazard mitigation. Frontiers in Built Environment, 2025; 10:1499848. DOI 10.3389/fbuil.2024.1499848.

[7] McGuire R.K., FORTRAN computer program for seismic risk analysis (Report No. 76-67). USGS Open-File Report. 1976. DOI 10.3133/ofr7667.

[8] McGuire R. K., Seismic hazard and risk analysis. Earthquake Engineering Research Institute, MNO-10, 2004.

[9] Field E. H., Arrowsmith J. R., Biasi G. P., Bird P., Dawson T. E., Felzer K. R., et al., Uniform California Earthquake Rupture Forecast, Version 3 (UCERF3)—The time-independent model. Earthquake Spectra, 2013; 29(3): 1021–1050. DOI: 10.3133/ofr20131165.

[10] Petersen M.D., Moschetti M.P., Powers P.M., Mueller C.S., Haller K.M., Frankel A.D., et al., Documentation for the 2014 update of the United States national seismic hazard map (Report No. 2014-1091). USGS Open-File Report. 2014. DOI 10.3133/ofr20141091.

[11] Woessner J., Danciu L., Giardini D., Crowley H., Cotton F., Grünthal G., et al., The 2013 European seismic hazard model: Key components and results. Bulletin of Earthquake Engineering, 2015; 13: 3553–3596. DOI 10.1007/s10518-015-9795-1.

[12] Mosca I., Sargeant S., Baptie B., Musson R.M.W. and Pharaoh T.C., The 2020 national seismic hazard model for the United Kingdom. Bulletin of Earthquake Engineering, 2022; 20: 633–675. DOI 10.1007/s10518-021-01281-z.

[13] Lindholm C., Bungum H., Ghione F. Meslem A., Huang C. and Oye V., Earthquakes and seismic hazard for Norway and Svalbard. Journal of Seismology, 2025; 2: 107–126. DOI: 10.1007/s10950-024-10270-z.

[14] Wang Y.J., Chan C.H., Lee Y.T., Ma K.F., Shyu J.B.H., Rau R.J., et al., Probabilistic seismic hazard assessment for Taiwan: TEM PSHA 2020. Terrestrial, Atmospheric and Oceanic Sciences, 2016; 27: 325-340. DOI 10.1177/8755293020951587.

[15] Pagani M., Hao K.X., Fujiwara H., Gerstenberger M. and Ma K.F., Appraising the PSHA earthquake source models of Japan, New Zealand, and Taiwan. Seismological Research Letters, 2016; 87(6): 1240–1253. DOI 10.1785/0220160101.

[16] Lu T.Y. and Chan C.H., Probabilistic seismic hazard assessment in complex fault systems: Exploring the longitudinal valley of Taiwan. Geoscience Letters, 2025; 12: 13. DOI 10.1186/s40562-025-00385-7.

[17] Xu W., Li X. and Gao M., Temporal distribution characteristics of earthquakes in Taiwan, China. Frontiers in Earth Science, 2023; 10: 930468. DOI: 10.3389/feart.2022.930468.

[18] Pagani M., Garcia-Pelaez J., Gee R., Johnson K., Poggi V., Styron R., et al., The 2018 version of the global earthquake model: Hazard component. Earthquake Spectra, 2020; 36: 226-251. DOI 10.1177/8755293020931866.

[19] Gouache C., Tinard P. and Bonneau F., Stochastic generator of earthquakes for mainland France. Applied Sciences, 2022; 12: 571. DOI 10.3390/app12020571. [20] Gardner J.K. and Knopoff L., Is the sequence of earthquakes in southern California, with aftershocks removed, Poissonian? Bulletin of the Seismological Society of America, 1974; 64: 1363–1367. DOI 10.1785/BSSA0640051363.

[21] Cheng C.T., Hsieh P.S., Lin P.S., Yen Y.T. and Chan C.H., Probability Seismic Hazard Mapping of Taiwan; in; Beer M., Kougioumtzoglou I.A., Patelli E. and Au S.K., eds., Encyclopedia of Earthquake Engineering. Springer, Berlin, Heidelberg, 2015: 1-25. DOI 10.1007/978-3-642-36197-5_100-1.

[22] Gutenberg B. and Richter C.F., Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 1944, 34: 185-188. DOI 10.1785/BSSA0340040185.

[23] Kagan Y.Y. and Knopoff L., Random stress and earthquake statistics: Time dependence. Geophysical Journal International, 1987; 88: 723-731. DOI 10.1111/j.1365-246X.1987.tb01653.x.

[24] Ellsworth W.L., Matthews M.V., Nadeau R.M., Nishenko S.P., Reasenberg P.A. and Simpson R.W., A physically-based earthquake recurrence model for estimation of long-term earthquake probabilities. Open-File Report 99-522, U.S. Geological Survey, Reston, VA, 22. Available at https://pubs.usgs.gov/of/1999/0522/ . (Accessed on 22 July 2025).

[25] Anagnos h. and Kiremidjian A.S., A review of earthquake occurrence models for seismic hazard analysis. Probabilistic Engineering Mechanics, 1988; 1: 3-11. DOI 10.1016/0266-8920(88)90002-1.

[26] Lin P.S. and Lee C.T., Ground-motion attenuation relationships for subduction-zone earthquakes in northeastern Taiwan. Bulletin of the Seismological Society of America, 2008; 98: 220-240. DOI 10.1785/0120060002.

[27] Lin P.S., Ground-Motion Attenuation Relationship and Path-Effect Study Using Taiwan Data Set. PhD Thesis, Institute of Geophysics, Central University, Taiwan, 2009.

[28] Utsu T., Representation and analysis of the earthquake size distribution: A historical review and some new approaches. Pure and Applied Geophysics, 1999; 155: 509-535. DOI 10.1007/S000240050276.

[29] Parsons T., Significance of stress transfer in time-dependent earthquake probability calculations, Journal of Geophysical Research, 2005; 110: B05S02. DOI 10.1023/2004JB003190.

[30] Nazaruddin D.A. and Duerrast H., 2012 Earthquake swarm in Phuket, southern Thailand. Chiang Mai Journal of Science, 2021; 48(2): 648-663.