2.1 Central Weather Bureau seismic network catalog

We used data from the seismic network catalog maintained by the Central Weather Bureau (CWB) of Taiwan (R.O.C.) (https://www.cwb.gov.tw/V7e/earthquake/seismic.htm and http://gdms.cwb.gov.tw/index.php, last access: July 2018). The completeness magnitude (M_c) of this catalog is estimated at approximately 2.0 in local magnitude (M_L) (Wu et al., 2008c). In an analysis of focal depth, Wu et al. (2008b) observed that the focal depth of approximately 80 % of earthquakes was shallower than 30 km. Accordingly, we used M_L 2.0 and 30 km as the threshold of magnitude and focal depth, respectively, to select the events to be used in the PI calculation.

2.2 P-alert network

We used the ground motion recordings from the P-alert network to verify the effectiveness of the real-time seismic hazard assessments from our model. The National Taiwan University (NTU) commenced development of the P-alert real-time strong motion network for EEW purposes in 2010 with the support of the Ministry of Science and Technology (MOST) (Wu, 2015). The devices of the P-alert network can record real-time three-component acceleration signals and publish alerts when the peak initial-displacement amplitude (Pd) or the peak ground acceleration (PGA) exceeds predefined thresholds (Wu et al., 2013, 2016b; Wu, 2015). Today, there are more than 600 P-alert stations in Taiwan, most located in elementary schools (Wu et al., 2013; Yang et al., 2018). We mainly adopted the P-alert waveform database maintained by the Taiwan Earthquake Research Center (TEC), and we used the data from the NTU as an auxiliary catalog (data from the P-alert network can be downloaded from the Data Center of the TEC at http://palert.earth.sinica.edu.tw/db/, last access: July 2018, or by contacting Yih-Min Wu at drymwu@ntu.edu.tw for access to the NTU catalog).

The distribution of the P-alert network is still not uniform (see Figs. 2b and 3b), despite the large number of seismic stations covering Taiwan. Obviously, this could cause problems, which will be discussed later.

3.1 Pattern informatics

Phase dynamics is the physical fundamental of the PI method, which describes changes in a system by rotation of the state vector in the Hilbert space (Rundle et al., 2002, 2003). The evolution of the state vector in a dynamic fault system is suggested to be related to stress accumulation and release (Chen et al., 2006). The computational steps we adopted here are a modified version developed by Chang et al. (2016) and Chang (2018) to improve the spatiotemporal resolution of the PI method. The research area (21–26^∘ N, 119–123^∘ E) was divided into boxes of grid size $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$, with each box being indicated by parameter x_i. Because of the M_c and the distribution of the focal depth (mentioned in Sect. 2.1), we selected all the events with M_L≥2.0 and depth ≤30 km. In PI computation, t₁ and t₂ represent the beginning and the end of a change interval, respectively, with the length of a change interval being 4 years. The start time of calculation, t₀, is defined as 12 years before t₂. Then, t_b is a sampling reference time between t₀ and t₁. The t_b starts from t₀ and shifts forward 3 d in each calculation until the length of time between t_b and t₁ is a half change interval. The forecasting interval, t₃, starts after t₂ (Chang, 2018). The seismicity rate in the period t_b to t (t_b to t₁ and t_b to t₂) can be expressed as

We conservatively considered the earthquake number, n, occurring in the x_i and its eight neighboring boxes. The rate change during the change interval can be expressed as

$\mathit{S}(x_i,t_b,t)$ is a vector in a Hilbert space that records present seismic activity, so that ΔS can be interpreted as an angular drift of S (Rundle et al., 2002; Tiampo, 2002). To reduce the time-dependent background seismicity, we used the temporal standard score normalizing ΔS and obtained $\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{S}}$. To compare the high and low levels of seismicity rate change in each grid box at the same t_b, we subsequently used the spatial standard score normalizing $\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{S}}$ and obtained $\mathrm{\Delta }\widehat{\mathit{S}}$. The average of the absolute value at all t_b points in each x_i is

Then, the mean squared change in probability

was computed (Chen et al., 2005; Chang et al., 2016). We further divided the magnitude range of earthquakes into several segments to separately calculate the relative probabilities ΔP(x_i). The divided magnitude range is from magnitude 2.0 with a window length of 0.5 magnitude, and it shifts forward by 0.2 each time. Then, we calculated the relative probability each time, such as ΔP(x_i)_2.0–2.5, ΔP(x_i)_2.2–2.7. Finally, we multiplied all the relative probabilities.

ΔP_M to forecast the occurrence of earthquakes is referred to as the modified pattern informatics method (Chang, 2018). According to Chang et al. (2016), the forecasting interval of the PI method reaches 90 d. Finally, the PI method produced a forecasting probability distribution of seismic sources for M_L≥5.0 within the forecasting interval.

3.2 Real-time PSHA

In the traditional PSHA framework (Cornell, 1968; Wang et al., 2016), the probability of an earthquake occurrence follows the Poisson process and the average recurrence interval for an annual frequency of exceedance can be expressed as

where $f_{M_i}\left(m\right)$ and $f_{R_i}\left(r\right)$ are the probability density functions of magnitude and distance, respectively; $P\left(Z>z\mathrm{|}m,r\right)$ is the conditional probability of ground motion Z exceeding a specified value z for a specific magnitude m and distance r. $\dot{N_i}$ is the annual occurrence rate of earthquakes and is described by the truncated exponential model (Cosentino et al., 1977) and the characteristic earthquake model (Schwartz and Coppersmith, 1984). Finally, to consider all scenarios, the total probability of N_s earthquakes is summarized in a given region.

In real-time PSHA, the occurrence rate of earthquakes used in the traditional PSHA framework is replaced by seismic forecasting probability to achieve spatiotemporal variability in the hazard assessment. Then, considering the gridded space, real-time PSHA can be expressed as

where $P_{M_i,{\text{Loc}}_i}(m,\text{loc})$, the forecasting probability distribution, is a function of magnitude and location. It specifies an occurrence probability for specific magnitude, M_i, at each spatial location, Loc_i. The summations are to consider the whole of the contribution from any possible magnitude, M_s, and location, Loc_s. We adopted the forecasting probability from the PI method as $P_{M,\text{Loc}}(m,\text{loc})$. Loc refers to x_i in the PI method. The forecasting probability of the PI method presents a distribution of cumulative forecasting probability for M_L≥5.0. We referred to the average character of the Gutenberg–Richter law in Taiwan (Gutenberg and Richter, 1944; Wang et al., 2015) to convert it into the probability density function (PDF). It can correspond to the specific magnitude conditions for $P(Z>z\mathrm{|}m,\text{loc})$. To evaluate the ground motion, we used the GMPE published by Lin et al. (2012),

which was also adopted for the Taiwan PSHA in Lee et al. (2017). In Eq. (7), C₁ to C₈ and H are the regression coefficients (Table 2); R is the closest distance (km); F_NM and F_RV represent the earthquake type, namely F_NM=1 and F_RV=0 for a normal fault earthquake and F_NM=0 and F_RV=1 for a reverse fault earthquake. In this GMPE, earthquake type is regarded as an important parameter. However, the division of seismic sources in the PI method is no longer based on the geological classification but on the grid box, x_i. Considering that most faults in Taiwan are reverse faults (Shyu et al., 2016), we adopted the reverse fault parameter setting for the entire research area. V_s30 is eclectically assigned V_s30=760. Using the conversion equation from Lin and Lee (2008), which was adopted in Lin (2012), turns M_L into M_w. Finally, the forecasting maximum PGA from real-time PSHA is transferred to seismic intensity according to the seismic intensity scale of the CWB listed in Table 3 (Wu et al., 2003; Central Weather Bureau, 2020). This implies that the seismic intensity forecasting map presents the maximum seismic intensity that every site will encounter over the following 90 d.

Table 2Coefficients in the GMPE.

Download Print Version | Download XLSX

Table 3Seismic intensity scale of CWB.

Download Print Version | Download XLSX

3.3 Performance verification

3.3.1 Receiver operating characteristic curve

The receiver operating characteristic (ROC) diagram is a binary classification model used widely as a tool to quantify the performance of earthquake prediction (Holliday et al., 2006; Nanjo et al., 2006; Wu et al., 2016a). We used the ROC diagram as an objective quantitative indicator to evaluate the performance of the seismic forecasting probability computed by the PI method. For each box, x_i, there are four situations (parameters) when comparing forecasting hot spots and target earthquakes. Namely, a means any target earthquake in a hot spot, b means no target earthquake in a hot spot, c means no hot spot but at least one target earthquake, d means no target earthquake and no hot spot. The true positive rate (TPR) is defined as $a/(a+c)$ and the false positive rate (FPR) is defined as $b/(b+d)$. The values of a, b, c, and d change with the threshold of forecasting probability and, therefore, TPR and FPR change as well. The value of the area under the ROC curve (AUC) varies between 0 and 1. AUC=1 is a perfect prediction and AUC=0.5 is a random guess. For each PI forecasting map, we generated 1000 random test maps by re-distributing the hot spots randomly over the research area to examine the possibility that a specific distribution of hot spots could be generated by chance. In Fig. 1c and d, the blue line is the 95 % confidence interval based on 2 standard deviations. The standard deviation is calculated by the random test results in each bin of the x axis. The 95 % confidence interval helps to differentiate the distributing range of random tests and the significance of the forecasting probability.

Figure 1Panels (a) and (b) show the forecasting probability maps of the Meinong earthquake and the Hualien earthquake, respectively. Panels (c) and (d) are the ROC curves of (a) and (b), respectively. Red, gray, blue, and black curves represent the forecasting probability map, random tests, 95 % confidence interval, and the average of random tests, respectively.

4.1 Forecasting earthquake occurrences

Figure 1a and b show the forecasting probability maps computed with the PI method, and Fig. 1c and d show the corresponding forecasting performance verified by the ROC tests. In the case of the 2016 Meinong earthquake, t₀, t₁, and t₂ are 31 January 2004, 31 January 2012, and 31 January 2016. In the case of the 2018 Hualien earthquake, t₀, t₁, and t₂ are 31 January 2006, 31 January 2014, and 31 January 2018. The forecasting intervals of both cases are 90 d after t₂. The cyan stars in Fig. 1a and b indicate the main shock of the 2016 Meinong and 2018 Hualien earthquakes and the largest earthquake in the forecasting interval. The gray circles in Fig. 1a and b are the earthquakes with magnitude M_L≥5.0 in the forecasting interval, with more-detailed information about these earthquakes presented in Table 1. Notably, both main shocks and most large earthquakes are located in or in close proximity to hot spots. Overall, simply from visual inspection, the performance of the PI forecasting probabilities appeared satisfactory.

In Fig. 1c and d, the red curves are located far above the blue curves (95 % confidence interval). The AUCs of the red curves are 0.91 and 0.94 and are apparently larger than the AUCs of the blue curves, which are 0.73 and 0.70. The ROC tests quantitatively verified that the performance of the PI forecasting probability was significant and that these patterns were not generated by chance by the random distribution of hot spots. Both distributions of hot spots were found to be physically meaningful. In view of the above, we were able to use these probability maps as the function of earthquake occurrence rate in subsequent calculations for real-time PSHA.

4.2 Real-time PSHA

In Figs. 2 and 3, panel (a) shows the map forecasting the maximum seismic intensity estimated by real-time PSHA for the forecasting interval, and panel (b) shows the map indicating the maximum seismic intensity recorded by the P-alert network during the forecasting interval. To ensure that it was the absolute maximum intensity during the forecasting interval, we used only the stations that had recorded all the target events (M_L≥5.0) in the forecasting interval. Although there are over 600 P-alert stations distributed widely in Taiwan, some boxes do not contain any station, e.g., the Central Mountain Range (see Fig. 5a and b). Therefore, we estimated the intensities in these boxes by interpolating. However, clearly, our strategy generated an artificial effect, which will be shown later.

Figure 2The 2016 Meinong earthquake. (a) Map of forecasted maximum seismic intensity by real-time PSHA. The forecasting interval of seismic intensity is 90 d. (b) Map of maximum seismic intensity recorded by the P-alert network. Black and white triangles indicate the P-alert stations that we used and did not use, respectively, in the verification. (c) Difference in seismic intensity between the forecast and the record. The cyan star represents the Meinong earthquake, and the gray circles represent the earthquakes with magnitude M_L≥5 in this forecasting interval.

Figure 3The 2018 Hualien earthquake. (a) Map of forecasting maximum seismic intensity. (b) Map of maximum seismic intensity recorded by the P-alert network. (c) Difference in seismic intensity between the forecast and the record. The cyan star represents the Hualien earthquake.

Comparing Fig. 2a and b, we suggest that both seismic intensity distributions are remarkably similar. An apparent deviation between the forecasted seismic intensities and the recorded values occurs in southwestern Taiwan, particularly the area closer to the 2016 Meinong main shock. Figure 2c shows the difference in seismic intensity between Fig. 2a and b, with the blue and red colors indicating that the forecasting value in a box was underestimated or overestimated, respectively. Most boxes have an intensity difference in the range −1 to 1, but some boxes in southwestern Taiwan are underestimated, with the differences being mostly 2 or even up to 3.

Comparing Fig. 3a and b, we suggest that both seismic intensity distributions are still extremely similar. In this instance, an apparent deviation between the forecasted seismic intensities and the recorded values occurs in southern Taiwan and a part of the southwestern area. Figure 3c shows that most boxes in southern Taiwan have a smaller recorded intensity, with the recorded intensities in a part of southwestern Taiwan being larger than the forecasting values.

Figure 4 shows the verifications generated by the APHR to quantitatively evaluate the performance of forecasting the seismic intensity. We considered the denominator of two classifications in Eq. (8), i.e., the total number of P-alert stations and the total number of boxes in the research area. The results are indicated by “P-alert” and “Map”, respectively, in Fig. 4. When comparing forecasting intensity with recorded value, both cases “forecasting = recorded” and “forecasting = recorded + 1” indicate “successful forecasting”. However, defining the tolerance range, which depends on the perspectives and allowances of different users, is debatable (Hsu et al., 2018). In this study, we tolerated an overestimation of 1 intensity rather than underestimation, as, in relation to the prevention or mitigation of earthquake disasters, “overestimation” was considered preferable to “underestimation”.

Figure 4Performance test of APHR. The red line indicates the forecasts of real-time PSHA, the gray circle indicates the result of a random test by randomly redistributing seismic intensities, and the blue error bar indicates the interval with 2 standard deviations over all random tests.

Download

All the red lines are above the maximum hit rate of the random tests and higher than 0.5, not to mention the random guesses of the eight choices of the seismic intensity scale on each station or box. This implies that the forecasting ability of the generated seismic intensity maps is significantly effective and that this satisfactory performance could not be ascribed to chance. Furthermore, both hit rates of the “P-alert” cases are higher than the rates of the “map” cases. However, this result could be attributed to the influence of the artificial effect generated by the interpolation of seismic intensity from the P-alert data of nonuniform distribution. Finally, it should be emphasized that we focused only on earthquakes with M_L≥5 and we cannot deny the possibility that a M_L<5 earthquake could cause large seismic intensity in the near field.