Spatiotemporal seismicity pattern of the Taiwan orogen

We investigate the temporal and spatial seismicity patterns prior to eight M > 6 events nucleating in different regions of Taiwan through a region–time–length algorithm and an analysis of a self-organizing spinodal model. Our results show that the spatiotemporal seismicity variations during the preparation process of impending earthquakes display distinctive patterns corresponding to tectonic settings. Q-type events occur in southern Taiwan and experience a seismic quiescence stage prior to the mainshock. A seismicity decrease of 2.5 < M < 4.5 events occurs around the relatively high b-value southern Central Range, which contributes to the accumulation of tectonic stress for preparing for the occurrence of the Q-type event. On the other hand, A-type events occur in central Taiwan and experience a seismic activation stage prior to the mainshock, which nucleates on the edge of the seismic activation area. We should pay attention when accelerating seismicity of 3 < M < 5 events appears within the low b-value area, which could promote the nucleation process of the A-type event.

1 Introduction

Seismic activity is related to spatiotemporal variations in the stress field and state, and seismicity changes prior to a large earthquake have been widely observed through different techniques, e.g., b-value analysis (Chan et al., 2012; Wyss and Stefansson, 2006), noncritical precursory accelerating seismicity theory (PAST) (Mignan and Giovambattista, 2008), pattern informatics (PI) algorithm (Rundle et al., 2003; Chen et al., 2005), the region–time–length (RTL) algorithm (Chen and Wu, 2006; Wen et al., 2016), and the analysis of self-organizing spinodal (SOS) model (Rundle et al., 2000). Previous studies have mostly focused on a significant earthquake; therefore, it is not easy to understand whether the properties of seismic activation and quiescence patterns respond to regional tectonic stress.

The Taiwan orogenic belt, which is an active and ongoing arc–continent collision zone as a result of the Philippine Sea Plate (PSP) obliquely colliding with the Eurasian Plate (EP), is particularly complex due to the two adjacent subduction zones, the Ryukyu trench and Manila trench to the northeast and south of the island, respectively (Suppe, 1984; Yu et al., 1997). The frequent and significant seismic activities as well as a rapid convergence rate of 85 to 90 mm yr⁻¹ are well observed by the island-wide GPS and seismic networks (Fig. 1). The growth of the Taiwan orogenic belt shows propagation from north to south due to oblique plate convergence and opposing subduction in the southern and northern parts of Taiwan (Suppe, 1984). The central part of Taiwan, which is experiencing rapid to full collision, mainly consists of the Coastal Range, Central Range and Western Foothills (Shyu et al., 2005a, b). A myriad of active and thin-skinned structures are the products of the accretion of the continental sliver to the continental margin. In southern Taiwan, the EP subducting eastward beneath the PSP is in a stage of incipient arc–continent collision (Kao et al., 2000; Shyu et al., 2005a, b). The northwest domain of southern Taiwan, which represents the southern tip of the fold-and-thrust belt in the coastal plain and foothill region and shows very low seismicity, mainly consists of Miocene shallow marine deposits and a Pliocene–Pleistocene foreland basin as well as mudstones.

Figure 1 Horizontal velocities from 2002 to 2017 (Chen et al., 2018) and seismicity between 1991 and 2018. The white star shows the location of the 1999 Chi-Chi earthquake, and the focal mechanisms determined by the Global CMT solution represent the locations of the investigated events in this study. The active faults (thick lines) identified by the Central Geological Survey of Taiwan are also shown.

Over the last 2 decades, several moderate earthquakes have occurred with various seismicity patterns and in GPS velocity field regions. We investigate the temporal and spatial seismicity patterns prior to eight M > 6 events nucleated in different regions of Taiwan through the RTL algorithm and analysis of the SOS model. Our attempt is not to catch the seismic precursor but to focus on the seismicity changes related to the regional tectonics, which might become useful hints for potential seismic hazard assessments. The results show that the temporal and spatial seismicity (2.5 < M < 5) variations during the preparation process of impending earthquakes could display distinctive patterns corresponding to the tectonic setting.

2 RTL algorithm and data

The region–time–length (RTL) algorithm (Sobolev and Tyupkin, 1997, 1999) is a statistical technique to detect the occurrence of seismic quiescence and activation by taking into account the location, occurrence time and magnitude of earthquakes. The RTL value is defined as the product of the three dimensionless factors, R, T and L:

$$\begin{array}{}\text{(1)}& {\displaystyle }R(x,y,z,t)=\left[\sum _{i=\mathrm{1}}^n\exp \left(-{\displaystyle \frac{r_i}{r_{\mathrm{0}}}}\right)\right]-R_{\mathrm{bk}}(x,y,z,t)\text{(2)}& {\displaystyle }T(x,y,z,t)=\left[\sum _{i=\mathrm{1}}^n\exp \left(-{\displaystyle \frac{t-t_i}{t_{\mathrm{0}}}}\right)\right]-T_{\mathrm{bk}}(x,y,z,t)\text{(3)}& {\displaystyle }L(x,y,z,t)=\left[\sum _{i=\mathrm{1}}^n\left({\displaystyle \frac{l_i}{r_i}}\right)\right]-L_{\mathrm{bk}}(x,y,z,t),\end{array}$$

where r_i is the distance between the investigated point (x, y, z) and the ith prior event (with the occurrence time t_i and rupture length l_i). n is the number of prior events that occurred in a defined space–time window with r_i≤2r₀ (r₀, characteristic distance) and $(t-t_i)\le \mathrm{2}{t}_{\mathrm{0}}$ (t₀, characteristic timespan). Rupture length l_i is a function of earthquake magnitude (M_i), $\log l_i=\mathrm{0.5}{M}_i-\mathrm{1.8}$ (Kasahara, 1981). The weighted RTL value reflects the deviation from the background seismicity level (R_bk, T_bk and L_bk) with negative values for seismic quiescence and positive values for activation. r₀ characterizes the decreasing influence of more distant events, and t₀ describes the reducing influence rate of the preceding events as the time of calculation moving on. To diminish the ambiguity in determining the characteristic parameters, we follow the systematic procedure of correlation analysis over pairs of RTL results proposed by Huang and Ding (2012) to obtain the optimal model parameters, $\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}$ and $\stackrel{\mathrm{̃}}{t_{\mathrm{0}}}$, of each event. Details of this technique of correlation analysis are described in Appendix A. We calculate various combinations of r₀ (ranging between 25 and 80 km with a step of 2.5 km) and t₀ (ranging between 0.25 and 2.0 years with a step of 0.05 years). As the correlation coefficient criterion C₀ is set, we can calculate the ratio W (or weight) of the combination with correlation coefficients equal to or larger than C₀ for each model parameter of r_0i (i=1–m​​​​​​​; m=23) and t_0j (j=1–n; n=36).

After testing many criterion sets, the criterion coefficient C₀=0.6 and criterion ratio W₀=0.5 are acceptable for each event, which represents at least 50% of the total combination pairs with correlation coefficient $C\ge C_{\mathrm{0}}=\mathrm{0.6}$. Then, we obtain the average $\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}=\mathrm{49.6}$ km and average $\stackrel{\mathrm{̃}}{t_{\mathrm{0}}}=\mathrm{1.16}$ years. These model parameters are similar to those of previous studies for Taiwan (Chen and Wu, 2006; Wen et al., 2016; Lu, 2017; Wen and Chen, 2017).

For statistical analyses, catalog completeness is an important factor. Since 1991, the Taiwan Telemetered Seismographic Network (TTSN) (Wang, 1989) has merged with the Central Weather Bureau (CWB) seismic network and updated to an integrated earthquake observation system, named the Central Weather Bureau Seismic Network (CWBSN). Wang et al. (1994) pointed out that most shallow earthquakes occurring in Taiwan are distributed at depths less than 35 km. According to previous studies (Wu and Chiao, 2006; Wu et al., 2008; Wen et al., 2016; Hsu et al., 2021), we used the earthquake catalog maintained by the CWB for the entire Taiwan area with M ≥ 2.5 and depth ≤ 35 km between 1991 and 2018 and applied a declustering procedure proposed by Gardner and Knopoff (1974). Considering a sufficient background seismicity and minimizing the influence of the 1999 Chi-Chi earthquake, we only selected the M > 6 inland earthquakes between 2003 and 2016 in Taiwan. Since two events occurring in a close space–time window would show high similarity in RTL function (Lu, 2017), we excluded the event occurring within $\mathrm{2}\stackrel{\mathrm{̃}}{t_{\mathrm{0}}}$ and $\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}$ with respect to the last M > 6 events. For example, two M > 6 events within a distance of 10 km struck the Nantou area on 27 March and 2 June 2013, and we only analyzed the former event. Therefore, we have eight qualified M > 6 events, as listed in Table 1.

Table 1 Earthquake parameters for the investigated events determined by the CWB.

3 Results

3.1 Temporal seismicity variation

The temporal variation in the RTL function represents the different stages of seismicity rate change at the target location with respect to the background level. For consistency, we adopt a 10-year catalog as the background for each investigated event. Figure 2 shows the temporal variation in the RTL functions prior to the investigated events. We can see that before the occurrence of the investigated event, both seismicity changes are observed: the seismic quiescence stage for nos. 1, 2, 5 and 8 (Q-type events hereafter) and the seismic activation stage for nos. 3, 4, 6 and 7 (A-type events hereafter). Q-type events occurred at different locations in southern Taiwan, and most, 3 among 4, of their temporal RTL functions exhibit the seismic quiescence stages during 2002–2004, which was before the occurrence of the 2003 Chengkung earthquake, i.e., event no. 1. The seismicity increase (activation stage) took approximately 2 years following the 2003 Chengkung mainshock (event no. 1). We note that the length of the seismic quiescence stage prior to the Q-type event might correspond to the magnitude.

Figure 2 Temporal variation of the RTL function (blue line) for (a) Q-type events and (b) A-type events. The orange curves and vertical axes on the right represent the enlarged RTL functions of event nos. 3 and 4. The vertical dashed red lines mark the seismic quiescence stage, and the vertical dashed green lines mark the seismic activation stage. The bar chart represents the occurrence time of M ≥ 6.0 events within a distance of 2r₀ from the target event; each number above the bar is the magnitude.

A-type events all occurred in central Taiwan and were located within $\mathrm{2}\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}$ with respect to the 1999 Chi-Chi earthquake. Figure 3 shows the declustered seismicity distribution as a function of time and latitude. Significant seismicity followed the 1999 Chi-Chi earthquake north of 23^∘ N. Since the background seismicity of event nos. 3 and 4 started from 1 January 1999, the RTL functions were obviously affected by the occurrence of the 1999 Chi-Chi earthquake. Therefore, we enlarge the vertical axis to accentuate the seismicity variation prior to event nos. 3 and 4. As shown in Fig. 2, the temporal RTL functions of A-type events mostly show a seismic activation stage between 2004 and 2006, which corresponds to the seismicity increase following the 2003 Chengkung mainshock (event no. 1). However, for the A-type event, we could not see the relationship between the length of the seismic activation stage and the magnitude.

Figure 3 Map view of the earthquake b value and declustered seismicity distribution as a function of time and latitude. The white star indicates the 1999 Chi-Chi earthquake, and the black stars (or focal mechanisms) represent the investigated events in this study. The black arrows indicate the seismicity boundaries. The major geological units in Taiwan are marked by gray curves and labeled from A to G.

3.2 Spatial seismic activation–quiescence distribution

Since Q-type and A-type events are located in southern and central Taiwan, respectively, it would be worth examining the spatial pattern of their abnormal seismicity stages. Wen and Chen (2017) pointed out that various seismic activation or quiescence processes of about 2–4 years were found prior to some events occurring in Taiwan (Chen and Wu, 2006; Wen et al., 2016; Wu et al., 2008). Thus, for consistency, we only consider the last abnormal stage within 4 years prior to the investigated events, as marked by red vertical lines for the quiescence stage of Q-type events and green vertical lines for the activation stage of A-type events. Then, we calculate the summation of the selected period to generate the seismic quiescence–activation distribution. Considering the definition of the weighted RTL function, a sufficient amount of background seismicity should be regarded as a criterion (Wen and Chen, 2017). Using the declustered catalog from 1991 to 2016, we set up two conditions similar to those of Wen and Chen (2017) for each grid to strengthen the reliability: (i) the total number of events within the grid area of 0.1^∘ × 0.1^∘ must be more than 26 (i.e., at least 1 event occurred every year on average), and (ii) the total events within a circle of 2r₀ in radius must be more than 9360 (i.e., at least 30 events occurred every month on average). For each event, we normalize the spatial distribution based on the summed result. The spatial seismic activation–quiescence map provides the information of influence of surrounding seismicity state to the target event during the abnormal stage. Similar to previous studies (e.g., Huang et al., 2001; Huang and Ding, 2012), Fig. 4 shows that Q-type events mostly occurred on the edge of the seismic quiescence area, and seismic activation appeared around the A-type events.

Figure 4 (a) The summed and normalized seismic quiescence map for the selected time window of the temporal RTL function of Q-type events, and (b) the summed and normalized seismic activation map for the selected time window of the temporal RTL function of A-type events. Stars represent the locations of the investigated events. The active faults (thick lines) identified by the Central Geological Survey of Taiwan are also shown.

4 Discussion

4.1 Spatiotemporal characteristics of seismicity changes

The RTL analysis accounts for the background seismicity prior to the investigated event. Therefore, the RTL analyses account for almost the same background period for event nos. 3 and 4 (1999–2009) and for event nos. 6 and 7 (2003–2013), respectively. As the temporal RTL functions show the seismic activation stage prior to the mainshocks during a similar period, we could expect similar seismic activation maps for event nos. 3 versus 4 and event nos. 6 versus 7, as shown in Fig. 4. Furthermore, the seismic quiescence stage of event no. 5 occurred in a similar period as the seismic activation stage of event no. 3 (Fig. 2), and the seismic quiescence area of event no. 5 complements the seismic activation area of event no. 3 (Fig. 4). In contrast, although event nos. 3 and 7 occurred at close locations, the difference in the 10-year background period affects the weighting of the deviation. For example, the seismic quiescence stage during 2007–2009 shown in the temporal RTL function of event no. 7 (Fig. 2) is evaluated as the background seismicity level (RTL value is equal to zero) in the temporal RTL function with respect to event no. 3. On the other hand, Wen and Chen (2017) pointed out that an abnormal seismic stage derived with various background periods cannot be produced by chance. The temporal RTL functions of five events (nos. 1–5 in Fig. 2) accounting for different background periods all exhibit the seismic quiescence stage before the occurrence of event no. 1. This phenomenon is consistent with the seismic quiescence map of event no. 1 (Fig. 4) and the Z-value map of Wu et al. (2008) in which the seismic activity decreased during 2002–2003 for a large area in Taiwan. In addition, the widespread seismic activation distribution of nos. 6 and 7 (Fig. 4) also responded to the seismic activity increase during 2011–2012 (nos. 6–8 in Fig. 2). Wen et al. (2016) suggested that, after the 2010 Jiashian earthquake (event no. 5), the 2-year seismicity increase is caused by the increase in Coulomb stress change, which is consistent with the seismic activation period in the temporal RTL function of event no. 8.

Rundle et al. (2000) proposed the self-organizing spinodal (SOS) model for characteristic earthquakes and suggested that small earthquakes occurred uniformly at all times, while the occurrence rate of intermediate-sized earthquakes varied during the earthquake cycle. Chen (2003) investigated the SOS behavior of the 1999 Chi-Chi earthquake and proposed the seismic activation of moderate-size (5 < M < 6) events prior to the mainshock. Here, we also calculate the cumulative frequency–magnitude distributions for these eight events using the same catalog periods of the RTL analysis. For each investigated event, we only compared the distribution diagrams of the long-term (background period) and abnormal seismic stages marked by dashed lines in Fig. 2, within a radius of 25 km with respect to the epicenter. As shown in Fig. 5, cumulative frequency–magnitude distributions of long-term seismicity (red dots) generally exhibit linear power law distributions. For the Q-type events, the cumulative frequency distributions of the seismic quiescence stage (black dots) appear to lack 2.5 < M < 4.5 events (Fig. 5a), and the lack of a level corresponds to the seismic quiescence distribution near the epicenter (Fig. 4). This indicates that within the seismic quiescence stage before the occurrence of the Q-type event, the quiescence of 2.5 < M < 4.5 activity contributes to the accumulation of tectonic stress. On the other hand, the cumulative frequency distributions of the seismic activation stage of the A-type events (black dots in Fig. 5b) show that the seismic activation of 3 < M < 5 events within the seismic activation stage before the occurrence of the A-type earthquake can be found, which is similar to the results of the 1999 Chi-Chi earthquake (Chen, 2003). Event nos. 6 and 7, which are located very close to the high seismic activation area (Fig. 4), display the more obvious increase in the number of 4 < M < 5 events during the seismic activation stage (Fig. 5b).

Figure 5 The cumulative frequency–magnitude distributions prior to the investigated events. Red and black dots represent the long-term and abnormal seismic stage marked in Fig. 2, respectively.

Event no. 4 occurred only one month later than event no. 3; however, the seismic activation stage of event no. 4 was much longer than that of event no. 3. Furthermore, the cumulative frequency distributions of the seismic activation stage of event no. 4 display a lower intercept (Fig. 5b), which represents the overall decreasing seismicity within this seismic activation stage. Here, we further divide the seismic activation stage of event no. 4 into three periods for discussion:

• i.

  P1 encompasses February 2008–March 2009 before the seismic activation stage of event no. 3.

• ii.

  P2 encompasses April–September 2009 matching the seismic activation stage of event no. 3.

• iii.

  P3 encompasses October 2009 between the occurrences of event nos. 3 and 4.

The seismic activation distributions in Fig. 6 are all normalized with respect to the maximum RTL value of the seismic activation distribution of event no. 4 through periods P1–P3. We can see that before the seismic activation stage of event no. 3 during February 2008–March 2009 (P1), the location of event no. 3 indeed shows no seismic activation, as exhibited in the temporal RTL function (Fig. 2b). On the other hand, for the location of event no. 4, the seismic activation remains through all three periods P1–P3. Combined with the overall decreasing seismicity indicated by the lower intercept in Fig. 5b, these results suggest that this seismic activation prior to event no. 4 was mainly contributed by the relatively accelerating activity of 3.5 < M < 4 events.

Figure 6 The summed seismic activation map for different periods of the seismic activation stage prior to event no. 4; all maps are normalized based on the summed results of P1–P3. Stars represent the locations of event nos. 3 and 4. The active faults (thick lines) identified by the Central Geological Survey of Taiwan are also shown.

4.2 Implication for the tectonic setting

Several major active faults in southwestern Taiwan have been identified, and most of them have been dominated by thrust movement. Some strike-slip structures, e.g., the Zuochen and Hsinhua faults, acted as the transfer structures between these thrust faults (Ching et al., 2011; Deffontaines et al., 1994, 1997; Rau et al., 2012). These transfer structures develop at around 23^∘ N, which is the northern limit of the Wadati–Benioff zone (Kao et al., 2000) and close to the seismicity boundary indicated in Fig. 3. Geodetic data displayed various rates and orientations of horizontal shortening with rapid uplift rates in southern Taiwan (Fig. 1), which might be caused by underplating beneath the Central Range sustaining crustal thickening and exhumation (Simoes et al., 2007). The seismic b value, which is the relative earthquake size distribution, can be derived from the Gutenberg–Richter relation (Gutenberg and Richter, 1944): $\log N=a-bM$, where constant a is related to seismicity and N is the number of earthquakes with magnitudes greater than M. In general, a high b value indicates a larger proportion of small events, and a low b value suggests that large earthquakes dominate over small ones. Using the same declustered catalog from 1991 to 2018, we search for events within a radius of 25 km with respect to the center of each grid (0.1^∘ × 0.1^∘). Only for the grids with more than 30 events, we calculate the b value using the weighted least-squares fitting method (Shi and Bolt, 1982) and the spatial distribution of b values, as shown in Fig. 3. The seismicity in the southern Central Range is active but shows significant heterogeneity in faulting types (Fig. 1; Chen et al., 2017; Wu et al., 2018), and relatively high b values suggest the predominance of small earthquakes in this region (Fig. 3 and red dots in Fig. 5a; Wu et al., 2018). Wen et al. (2016) found the decreased seismicity and increased Coulomb stress change in the southern Central Range prior to the 2010 Jiashian earthquake (i.e., event no. 5) and suggested both variations in Coulomb stress and seismicity rate play important roles in contributing to the nucleation process of impending earthquakes. The seismicity rate change can be considered a proxy for the stress state change (Dieterich, 1994; Dieterich et al., 2000), and this implies that the quiescence of seismicity contributes to the accumulation of tectonic stress. Since this relatively high b-value region in the southern Central Range has been observed to have a seismicity decrease (2.5 < M < 4.5 events) before the occurrence of Q-type events, it can be an indicator of stress change.

Many devastating earthquakes with surface ruptures have occurred in central Taiwan, including the 1935 M 7.1 Hsinchu–Taichung earthquake, the 1951 Longitudinal Valley earthquake sequence and the 1999 Chi-Chi earthquake (Lee et al., 2007; Chen et al., 2008; Lin et al., 2013). Hsu et al. (2009) derived the consistent orientations of principal strain rate and crust stress axes in central Taiwan, which implies that faulting style corresponds to stress buildup accumulating from interseismic loading (Fig. 1). They also pointed out that, for central Taiwan, small events tend to surround the locked fault zone, where major earthquakes might occur, during the interseismic period. The 1999 Chi-Chi earthquake ruptured the area near the end of the décollement with a high contraction rate (Dominguez et al., 2003; Hsu et al., 2003, 2009). In addition, similar to the 1999 Chi-Chi earthquake, the A-type events occurred in the low b-value area surrounded by small and active events. Chen and Wu (2006) derived the temporal RTL function of the 1999 Chi-Chi earthquake, showing a pattern similar to that of A-type events with the activation stage prior to the mainshock. Furthermore, Wu (2006) calculated the seismic activation map of the 1999 Chi-Chi event and found that the 1999 Chi-Chi mainshock occurred on the edge of the seismic activation area, which is a low b-value region. This is similar to the seismic activation maps of A-type events, which display the hot-spot pattern contracting within the low b-value area (Figs. 3 and 4). The nucleation of the A-type mainshock can be attributed to the perturbation of background seismicity (3 < M < 5 events) by the stress state change (Dieterich, 1994; Dieterich et al., 2000).

The cumulative frequency distributions of long-term seismicity in Fig. 5 show a b value of 0.8–1.0 around these eight events, which is consistent with the pattern shown in Fig. 3. However, the cumulative frequency distributions of long-term seismicity exhibit different trends of magnitudes larger than 4.5 for the two types of events. The seismicity for M > 4.5 events is lower in the area around the Q-type event but higher in the area around the A-type event. Event nos. 1, 2, 3 and 7 occurred in eastern Taiwan with an average GPS velocity of about 60 mm yr⁻¹ (Fig. 1), and the cumulative frequency distributions of long-term seismicity display a high intercept (Fig. 5). This rapid convergence rate generally remains in the western part of southern Taiwan, which indicates that only a little shortening is consumed from east to the west in southern Taiwan. This corresponds to the active seismicity of small earthquakes, as indicated by the high intercept of the cumulative frequency distributions of long-term seismicity for event nos. 1, 2, 5 and 8 (Fig. 5). Therefore, for the pre-collisional rapid and distributed convergence in southern Taiwan (Shyu et al., 2005a), the quiescence of 2.5 < M < 4.5 activity contributes to the accumulation of tectonic stress for preparing for the occurrence of the Q-type event. On the other hand, the shortening rate is obviously consumed in the mountainous area of central Taiwan. Therefore, the lowest intercept of the cumulative frequency distributions of long-term seismicity for event no. 4 (Fig. 5) reflects the slow GPS velocity and low seismicity in the western part of central Taiwan (Fig. 1). For central Taiwan, small events tend to surround the locked fault zone of the potential major events during the interseismic period, and the 1999 Chi-Chi earthquake is the case affected by the accelerating seismicity of moderate-size events and ruptured the area near the end of the décollement with a high contraction rate (Chen, 2003; Dominguez et al., 2003; Hsu et al., 2003, 2009). Tectonic stress accumulating from the interseismic loading with the perturbation of the accelerating activity of 3 < M < 5 events could promote the nucleation process of the A-type event. The mechanisms causing these different phenomena are not clear, and further study is still needed. For example, some studies using machine-learning-based earthquake detectors and template-matching techniques will be helpful to build a more complete earthquake catalog in Taiwan (Liao et al., 2021; Zhai et al., 2021) and to get more useful data on small earthquakes with a magnitude below 2.5.

5 Conclusion

Through statistical analyses of recent large earthquakes that occurred in Taiwan, we summarize various temporal and spatial seismicity patterns prior to the earthquakes that nucleated in different regions of Taiwan:

• Q-type events occurred in southern Taiwan, with the northern boundary of 23.2^∘ N, and experienced a seismic quiescence stage prior to the mainshock. A seismicity decrease of 2.5 < M < 4.5 events in the relatively high b-value southern Central Range could be an indicator of stress change related to the preparation process of such events.

• A-type events occurred in central Taiwan and experienced a seismic activation stage prior to the mainshock, which nucleated on the edge of the seismic activation area. We should consider when accelerating seismicity of 3 < M < 5 events appears within the low b-value area.

Our results show that the spatiotemporal seismicity variations during the preparation process of impending earthquakes could display a distinctive pattern corresponding to the tectonic setting.

Appendix A

In the systematic correlation analysis for searching the optimal model parameters, we calculate various combinations of r₀ (ranging between 25 and 80 km with a step of 2.5 km) and t₀ (ranging between 0.25 and 2.0 years with a step of 0.05 years). As the correlation coefficient criterion C₀ is set, we can calculate the ratio W (or weight) of the combination with correlation coefficients equal to or larger than C₀ for each model parameter of r_0i (i=1–m​​​​​​​; m=23) and t_0j (j=1–n​​​​​​​; n=36). Then, the contour map for the ratio W is generated, as shown in Fig. A1.

$$\begin{array}{}\text{(A1)}& W_{ij}={\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{m}I\left(C_{ik}\ge C_{\mathrm{0}}\right)+{\sum }_{l=\mathrm{1}}^{n}I\left(C_{jl}\ge C_{\mathrm{0}}\right)}{m+n}},\end{array}$$

where the logical function I(Φ) is defined as

$$\begin{array}{}\text{(A2)}& I\left(\mathrm{\Phi }\right)=\left\{\begin{array}{ll}\mathrm{1},& \mathrm{\Phi }\mathrm{is}\mathrm{true}\\ \mathrm{0},& \mathrm{otherwise}\end{array}\right..\end{array}$$

As the criterion ratio W₀ is set, the optimal model parameters, $\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}$ and $\stackrel{\mathrm{̃}}{t_{\mathrm{0}}}$, can be obtained by the following formulas:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^n\sum _{i=\mathrm{1}}^m{W}_{ij}I\left(W_{ij}\ge W_{\mathrm{0}}\right)r_{\mathrm{0}i}}{{\sum }_{j=\mathrm{1}}^n\sum _{i=\mathrm{1}}^m{W}_{ij}I\left(W_{ij}\ge W_{\mathrm{0}}\right)}}\text{(A4)}& {\displaystyle }\stackrel{\mathrm{̃}}{t_{\mathrm{0}}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m\sum _{j=\mathrm{1}}^n{W}_{ij}I\left(W_{ij}\ge W_{\mathrm{0}}\right)t_{\mathrm{0}j}}{{\sum }_{i=\mathrm{1}}^m\sum _{j=\mathrm{1}}^n{W}_{ij}I\left(W_{ij}\ge W_{\mathrm{0}}\right)}}.\end{array}$$

Using event no. 6 as an example, we considered criterion coefficient C₀ = 0.6 and criterion ratio W₀ = 0.5, which indicates that at least 50 % of the total combination pairs had a correlation coefficient $C\ge C_{\mathrm{0}}=\mathrm{0.6}$. Then, we obtained $\stackrel{\mathrm{̃}}{r_{\mathrm{0}}}$ = 50.0 km and $\stackrel{\mathrm{̃}}{t_{\mathrm{0}}}$ = 1.14 years (diamond in Fig. A1) by averaging the parameter values that passed the criterion.

In addition, Nagao et al. (2011) proposed the RTM algorithm to reduce the dual effect of the distance (r_i) by introducing the new factor

$$\begin{array}{}\text{(A5)}& M(x,y,z,t)=\left[\sum _{i=\mathrm{1}}^n\left(M_i\right)\right]-M_{\mathrm{bk}}(x,y,z,t),\end{array}$$

where M_i is the earthquake magnitude of the ith prior event. Here, we also calculate the RTM function of each investigated event with the same characteristic parameter set of the RTL model, and both functions display very similar trends with minor differences, as shown in Fig. A2. The reason for this could be that, for these eight events, no large earthquakes occurred in the vicinity of the epicenter. The bar chart in Fig. A2, which represents the occurrence time of M ≥ 6.0 events within a distance of 2r₀ from the target event, also supports this explanation.

Figure A1 Contour map of ratio W for various combinations of model parameters of r₀ and t₀, with C₀ = 0.6 for event no. 6. The diamond shows the optimal model parameters as selecting criterion ratio W₀ = 0.5.

Figure A2 Temporal variation of the RTL (solid line) and RTM (dotted line) functions for (a) Q-type events and (b) A-type events. The bar chart represents the occurrence time of M ≥ 6.0 events within a distance of 2r₀ from the target event; each number above the bar is the magnitude.