As mentioned above in the PI algorithm section, the parameters including the time parameters (t0, tb, t1, t2), space parameters of the scale of boxes in the studied region and the lower cut-off magnitude Mc from analyzing the FMD which is related to the magnitude earthquake that will occur during the forecast time interval are selected to be the most suitable in a specific area to achieve efficiency and accuracy in the PI anomalies related to a hazardous earthquake. For selecting proper parameters, previous studies are used as guidelines for choosing PI parameters as follows. Holliday et al. (2005) used the change time and the forecast time intervals to be both 10 years to forecast target earthquake M ≥ 5.0 in southern California and central Japan within box size as 0.1° x 0.1° while forecast target earthquake M ≥ 7.0 for worldwide within box size as 1° x 1° by PI method. Chen et al. (2005) applied the PI method to investigate the variations in seismicity pattern changes before the 1999 Chi-Chi, Taiwan, earthquake by selecting the change time interval around 6 years. For large earthquake forecasting in Japan, Nanjo et al. (2006) used the change time interval to be equal to the forecast time interval (10 years) to forecast earthquake event M ≥ 5.0 (Mc = 3.0) by selecting the box size as 0.1° x 0.1° (11 km.). For the retrospective test for the Wenchuan Earthquake, Sichuan-Yunnan region, Jiang and Wu (2010) adapted the variation parameter and then selected the change time and the forecast time intervals to be 5 years within grid size as 0.2° for forecasting target earthquake M ≥ 5.0 (Mc = 3.0). To study the precursory migration of anomalous seismic activity of the 2011 Tohoku earthquake by the PI method, Kawamura et al. (2013) used a grid size of 0.25° x 0.25° and a minimum Mc of 4.0 related to the forecast target magnitude earthquake (≥ 5.0). For the time intervals, 8 years (1 January 2000–1 January 2008) was chosen for the change time while the forecast time was chosen only to cover the time of the case study (11 March 2011). In addition to previous works that specify the parameters for PI analysis, there are also previous works that vary variables to find the most appropriate parameters. Zhang et al. (2013) studied the predictability of the PI method by testing with the M-8.0 Wenchuan and M-7.3 Yutian, varying both of the time intervals (1–10 years) and the grid size (1° x 1° and 2° x 2°). To obtain the suitable parameters, the relative operating characteristics (ROC) diagram and R score are used to verify. These results from the verification test show that the grid size of 2° x 2° and the forecast time interval of 8 years are the best parameters in this area.

Considering these previous studies, in this study, the study region is divided into a grid of square boxes with a size of 0.25° x 0.25°. The lower magnitude cutoff is set as Mw 4.2 and the target earthquake in the forecast time is set to a magnitude larger than 5.0 which is different from the assumption of Holidays et al. (2005), which is set to be greater than Mc + 2 because in this study area, there were not many earthquakes with magnitude larger than 6.2 and if chosen as the case studies, it will not enough. For the time parameters, the initial time was t0 = 1972, the beginning time was starting at 1980 (start of the completeness dataset), and the time step was every 5 years until 2015, while the change time and forecast time intervals were set to be the same interval that varied between 5 or 10 years for the retrospective forecasting of the PI algorithm. Thus, each PI value investigation that depended on time interval conditions was tested iteratively.

For finding the best condition of time interval parameters in this study area, the relative operating characteristics (ROC) diagram was chosen to verify the retrospective test with different time interval conditions between 5 and 10 years. The ROC diagram (Swets 1973; Jolliffe and Stephenson 2003; Mason 2003) is one of the well-known methods for validating the binary forecast, like the PI method. This technique is used to count the hit rate and false alarm rate from real earthquake activity compared with the alarm cells of the forecast region. The forecast is considered successful when determined by maximizing the proportion of earthquakes that occur in alarm cells while decreasing the fraction of alarm cells that occur outside of cells. As a consequence of this test, there are four categories for each cell: a) successful forecast: earthquake occurs in the hotspot cells or within its 8 Moore neighborhood (alarm and earthquake), b) false alarm: no earthquake occurs in a hotspot cell or within its 8 Moore neighborhood (alarm but not earthquake), c) failure to predict: an earthquake occurs in non-hotspot (no alarm but an earthquake) and d) successful forecast of non-occurrence: no earthquake occurs in a white non-hotspot cell (no alarm and no earthquake). In addition, values for the hot spot map consist of values a (Forecast = yes, observed = yes), b (Forecast = yes, Observe = no), c (Forecast = no, Observed = yes), and d (Forecast = no, Observed = no). The fraction of coloured boxes also called the probability of forecast of occurrence, is r = (a + b)/N, where the total number of boxes (N) is N = a + b + c + d. The hit rate (H) and a false alarm rate (F) are used to indicate the success of a forecast as follows: The hit rate (H) is H = a / (a + c); is the fraction of large earthquakes that occur on a hot spot and the false alarm rate (F) is F = b / (b + d); is the proportion of non-observed earthquakes that are incorrectly forecast. From the ROC diagram of both conditions as 5 and 10 years (Fig. 3), it can be seen that the hit rate in 5 years condition (represented in black colour) and in 10 years condition (represented in grey colour) are similar in 4 out of 6 case studies, but two case studies show a noticeable difference in the hit rate. The hit rate of two cases in the 5 years was significantly lower than in the 10 years. Additionally, even though the 10-year condition does not show the highest hit rate, these hit rates are more stable. For these reasons, the 10-year condition is the suitable time interval parameter for the PI investigation.

Figure 3 added here

For the retrospective test with the PI method, the spatial distribution map of PI value was evaluated and mapped into the hotspot areas (colour codes) that represent areas where ∆P is positive and is associated with a relative probability increase for strong earthquakes. For the colour code that represents the seismicity changes region, it is explained as follows. Warm tones represent the areas with large changes in seismicity (seismic quiescence and activation) during the change time interval, which indicates that these areas have a high chance of earthquake occurrence after the change interval, with the red colour indicating the highest probability of earthquake occurrence. The cold tones show the areas with small changes in seismicity, which indicates that these areas have a low chance of earthquake occurrence after the change interval (Fig. 4).

From the PI investigation in six case studies, these forecast results of earthquakes with magnitude ≥ 5.0 (shown by stars; red stars represent the epicenter of earthquakes that occur in or near the anomalous area and grey stars represent those that occur not close to the anomalous area.) with different time intervals from 1990 until 2020, respectively, are shown in Fig. 4a-f and explained as follows. From the Fig. 4c, the forecast hotspots for the occurrence of Mw ≥ 5.0 from 2000 to 2010, defined as t1 = 1 January 1990 and t2 = 1 January 2000 shows that 12 (red stars) out of 15 (grey stars) earthquakes occurred in or near the anomalous activity area, which is the Myitkyina and surrounding areas, Naypyidaw, and the Andaman Sea. In particular, the Myitkyina and Naypyidaw are high PI value zones, which are in the range of 0.8–0.9 and 0.4–0.6, respectively, and the epicenter of more than 6 events was located in a cluster of hotspots in these areas with 5 events (composed of Mw 6.0 (2008), Mw 5.4 (2000), Mw 5.2 (2007), Mw 5.1 (2005), and Mw 5.0 (2005)) in the Myitkyina and 3 events (composed of Mw 6.6 (2003), Mw 6.4 (2007), and Mw 5.0 (2009) in the Naypyidaw, as shown in red stars. Likewise, for the forecast interval from 2005 to 2015 (Fig. 4d), defined as t1 = 1 January 1995 and t2 = 1 January 2005, it is clear that the anomalous zones with a high PI value remain located in the same areas as Myitkyina (PI = 0.7–0.8), Naypyidaw (PI = 0.3–0.5). In addition to the earthquakes that occurred in 2000–2010, it was also found that in these anomalous zones, there were three new earthquakes during the periods 2010–2015, which were the magnitude of 5.5, 5.1, and 5.9 in 2011, 2012, and 2014, respectively. While the anomaly in the Andaman Sea remains in the same zones with a low PI value of approximately 0.1–0.2, it is considered that there is a low probability of an earthquake occurring in this area and is in accordance with a forecast that earthquakes will occur only in the surrounding of the anomalous zones, which is Mw 5.2 in 2008 and Mw 5.3 in 2009. Another anomalous zone that has increased is the Bago. Although no earthquakes occurred near the epicenter, earthquakes fall into the surrounding anomalous area as well.

Thereafter, select t1 = 1 January 1990 and t2 = 1 January 2000 to forecast the possible locations of the occurrence of earthquakes during the period 2010–2020, the anomalous zone around Bago still shows high PI values (0.2–0.5), and Mw 5.1 earthquake fell into this area in 2017 as seen in the red star in Fig. 4e. Moreover, it can be seen that earthquakes (Mw 5.0 in 2016 and Mw 5.1 in 2018) will occur in the Andaman Sea, which is close to the anomalously high PI zones, with a value of approximately 0.7–0.8. The vicinity of Myitkyina is still active and earthquakes as well. When moving the period to 5 years (t1 = 1 January 2005 and t2 = 1 January 2015), the forecast result for anomalous zones in 2015–2025 is still exposed near the same areas as in the period 2010–2020, namely the Myitkyina, Naypyidaw, Bago and Andaman Sea. Moreover, it can be seen that the anomalous zones around Bago and Andaman Sea have expanded and have high values of 0.4–0.6 and 0.7–0.9, respectively, which makes the forecast more accurate, as indicated by the Mw 5.1 Bago earthquake in 2017 and the Mw 5.0 Andaman Sea earthquake in 2016 as seen in Fig. 4f.

According to the spatial distribution of anomalous PI value during the forecast period 1990–2000, defined as t1 = 1 January 1980 and t2 = 1 January 1990, the anomalous high PI values were evident in the Andaman Sea, the northwest of the Myitkyina, and the Naypyidaw. Figure 4a shows that the epicenters of the 1991 earthquake located in an anomalous area with a high PI value in the Andaman Sea (PI = 0.6–0.7) that is fairly conforming to the high probability of occurrence. Meanwhile, the Mw 5.2 (1999) and Mw 5.1 (1995) earthquakes are located around the anomalous areas of Naypyidaw and the northwest of Myitkyina, with low PI values of 0.2–0.3, which is also fairly conforming to the low probability of occurrence. In addition to the case of forecast period 1990–2000, the case of forecast period 1995–2005, defined as t1 = 1 January 1985 and t2 = 1 January 1995, it was shown that although the epicenter of all earthquakes (represented by the stars) was not located in the anomalous areas, 4 earthquakes occurred in the surrounding of areas with anomalous PI values. Figure 4b shows that the magnitudes of 6.6 earthquakes in 2003 and 5.2 earthquakes in 1999 occurred close to the anomalous areas in Naypyidaw, while the magnitudes of 5.1 earthquakes in 1995 and 5.0 earthquakes in 2003 occurred close to the anomalous areas in the northwest of Myitkyina.

According to the retrospective test and the ROC diagram verification technique, selecting the characteristic parameters that include the time intervals equal to 10 years, a spatial grid size of 0.25°, and a target earthquake magnitude ≥ 5.0 (Mc = 4.2) will result in high PI values which indicates the seismogenic region reliably with reasonable accuracy to evaluate the prospective areas of upcoming major earthquake along the SFZ in the next session.

Figure 4 added here