More than one-third of the global population experiences earthquakes, with frequent slight or greater damage to lives and property from such events (Marti et al., 2019). However, at present, earthquake forecasting remains a major unsolved scientific problem and is not sufficiently reliable to rapidly predict earthquakes. Seismologists can also contribute to society by performing long-term probabilistic seismic hazard assessments. As a basic input for building codes, such assessments are crucial to ensure good building practices to save lives (Hiemer et al., 2013).

Mapped faults and past earthquakes have historically been the two main types of data used by seismic hazard models to estimate the likelihood of earthquake occurrence (Rhoades et al., 2017). Some short-term models have been built using only cataloged earthquakes, including the short-term earthquake probability (STEP) (Gerstenberger et al., 2005) and epidemic-type aftershock sequence (ETAS) (Ogata, 1988; Ogata, 1998; Gerstenberger and Rhoades, 2010) models, which are based on the decay of aftershock rates and the Omori–Utsu law (Utsu, 1961). Other useful medium- or long-term models, such as the proximity-to-past-earthquakes (PPE) (Kagan and Jackson, 1994; Jackson and Kagan, 1999) and the every-earthquake-a-precursor-according-to-scale (EEPAS) (Rhoades and Evison, 2004; Rhoades and Evison, 2005; Rhoades, 2007) models, rely only on past earthquakes. Meanwhile, long-term models such as the proximity-to-mapped-faults (PMF) model (Rhoades and Stirling, 2012) often use mapped faults with their slip rates. In addition, the increasing availability of crustal movement observational data, such as Global Positioning System (GPS) network data and derived models of crustal strain rate (Bird et al., 2015; Rhoades et al., 2017) are important to estimate the earthquake occurrence rate (Rhoades et al., 2015).

A parallel development is that of hybrid methods, which combine two or more different models to gain predictive skills from a single idea or data source. The results of retrospective and pseudo-prospective analyses have shown that some multiplicative ensembles provide statistically better forecasts than their constituent forecasting models (Rhoades and Stirling, 2012; Rhoades 2013; Gerstenberger et al., 2014; Rhoades et al., 2015; Rhoades et al., 2016). Therefore, the international Collaboratory for the Study of Earthquake Predictability (CSEP) supports methods to evaluate combinations of two or more individual models or to assimilate new gridded covariates into existing models (Schorlemmer et al., 2018).

The robustness of various earthquake forecast models cannot be fully described without first being assessed in diverse regions with distinct seismic features. The Sichuan-Yunnan region provides a natural earthquake forecast experiment because, following the M _S 8.0 Wenchuan earthquake, it became more seismically active. For example, two earthquakes of magnitude ≥7.0, the M _S 7.0 Lushan and M _S 7.0 Jiuzhaigou earthquakes, have occurred in this region in the last 10 years since the M _S 8.0 Wenchuan earthquake; however, no earthquakes had occurred over the same time period before the M _S 8.0 Wenchuan earthquake. Moreover, substantial infrastructure, such as reservoirs, railway networks, and electricity grids, is spread throughout the Sichuan-Yunnan region. While earthquakes of great intensity occurring due to natural processes cannot be avoided, forewarning can help minimize their often catastrophic and damaging impacts (Herrera et al., 2022). Much work is required regarding probabilistic seismic hazard assessments in this area. Therefore, we introduce individual probabilistic models based on past earthquakes or mapped faults and a hybrid model based on both to estimate the probability of earthquake occurrence; we then take a statistical approach to evaluate the efficiency of earthquake forecasting in the Sichuan-Yunnan region.

Section 2 includes a full description of the mapped faults and cataloged earthquakes that were used in this study. The two individual models and the hybrid model are introduced in Section 3. Section 4 and Section 5 present, successively, the results of the estimated seismicity rate for each model, along with an evaluation of their predictive effectiveness.

We chose the Sichuan-Yunnan region for testing because of the good monitoring capability in this region with respect to seismicity and fault movement. In the test area of the Sichuan-Yunnan region, which is marked in Figure 1, 1,234 earthquakes with magnitudes ≥4.0 and depth ≤50 km were recorded in the China Earthquake Networks Center (CENC) catalog from January 1, 1970, to September 10, 2022, including the M _S 8.0 Wenchuan earthquake, the largest earthquake in the test region. Of these, 14 target earthquakes in the test region with M ≥ 6.0 and depth ≤ 50 km that occurred during the three periods of May 20, 2008–May 20, 2011; January 1, 2013–January 1, 2016; and September 10, 2019–September 10, 2022, listed in Table 1, were used to evaluate the forecast efficiency of the three different models. Figure 1 shows the epicenters of the earthquakes that occurred in the test region.

The Chinese Earthquake Administration (CEA) and other groups launched the Project Crustal Movement Observation Network of China (CMONC) to monitor the crust dynamics using various techniques including GPS. Based on the GPS measurements and the fault locations in the test region, the Second Monitoring and Application Center of CEA calculated the fault rates in the Sichuan-Yunnan region and provided the fault-model datasets for this study (Figure 2). Figure 2 shows the higher slip rates in the Ganzi-Yushu fault (F5), the fault intersection segment of the Xianshuihe fault (F13), the Longmenshan fault (F10), the Jinpingshan fault (F15), the Anninghe fault (F14), and the Daliangshan fault (F12). The Qinghai-Tibet Plateau Block moves eastward continuously, blocked by the relatively stable Sichuan Basin, and the South China Block in the Sichuan-Yunnan region moves from east to south along these faults at a relatively high slip rate. In Figure 2, four aftershocks of the M _S 8.0 Wenchuan earthquake with magnitudes of M _S 6+ occurred near the northern segment (latitude >32°N) of the Longmenshan fault (F14), however, the segment of the fault does not have the expected relatively high slip rate (Shen et al., 2009; Xiong et al., 2021).

A hybrid model may underperform in a separate testing period if the seismicity rates that contribute to it are not well correlated with the earthquake locations or if the model is overfitted to a training period in which the seismicity rates are more closely correlated with the earthquake locations than those in the testing period (Rastin et al., 2022). Unaffected by any model fitting, an error diagram is a valuable tool to illustrate the relationship between the seismicity rate and the earthquake occurrence during various time periods and is also commonly used to evaluate earthquake forecast strategies (Molchan, 1991; Zechar and Jordan, 2008).

In this study, we defined the following earthquake forecast strategy: there will likely be earthquakes of M _S ≥ 6.0 occurring in the subsequent 3-year period at locations with relatively high seismicity rates based on the PPE, PMF, and PKS models.

The testing region was divided into multiple cells 0.5° × 0.5° in size. The seismicity rate of each cell was assigned the average of all the seismicity rates in the cells. A threshold was then specified and an alarm was generated if the seismicity rate of the cells was higher than the threshold. The fraction of the space-time occupied by the alarm (τ) of the Poisson is defined such that τ = ∑ i = 1 , n I S i > s l e v e l A i ∑ i = 1 , n A i (4a)

The fraction of earthquakes missed (ν) is defined such that ν = ∑ i = 1 , n I S i > s l e v e l E i ∑ i = 1 , n E i (5) where S(i) is the seismicity rate of the ith cell, slevel is the threshold value of the seismicity rate, A(i) is the area of the ith cell, and E(i) is the number of the targeted earthquakes at the ith alarm. We used the area skill score (ASS) and the probability gain to summarize the potential performance of the earthquake forecast strategy. The region above the error diagram inside the unit square, called ASS, serves as a summary of the prospective performance of the seismicity rate in alarm-based earthquake occurrence forecasting (Zechar and Jordan, 2008; Zechar and Jordan, 2010). The diagonal line indicates ASS = 0.5. A positive connection between the seismicity rate and the frequency of earthquakes is denoted by ASS >0.5, while ASS <0.5 denotes a negative correlation (Rhoades et al., 2017). The probability gain (Gain) is defined such that G a i n = 1 − ν τ (6)

The sample value of Gain corresponds to the slope of the line connecting (0,1) to (τ, ν), as shown in Figure 6. The diagonal line corresponds to a completely random guess, and Gain = 1. A higher Gain value indicates that fewer cells are needed for the same “hits,” which means better performance of the earthquake forecast strategy.

In Section 4, the spatial variation of the estimated seismicity rate under the PPE, PMF, and PKS models was calculated at three dates, May 20, 2008; January 1, 2013; and September 10, 2019. By adjusting the threshold of the seismicity rate from the highest to the lowest, the fraction of the space-time occupied by alarm and alarmed targeted earthquakes was obtained and their correlations are illustrated in a Molchan error diagram, as seen in Figure 6. All 14 targeted earthquakes during the three periods of May 20, 2008–May 20, 2011, January 1, 2013–January 1, 2016, and September 10, 2019–September 10, 2022, are included.

With respect to the individual models, the PPE model slightly outperforms the PMF model in the efficiency of earthquake forecast because the ASS value of the PPE model is higher than that of the PMF model. However, the PMF model has a more informative probability gain than the PPE model because when τ < 0.3 and the number of “alarmed” earthquakes is < 6, which indicates that less than six targeted earthquakes occur near the mapped faults having a relatively high slip rate, as seen in Figure 4. Conversely, when τ ≥ 0.3, the efficiency of PPE is better than that of PMF. The target catalogs include surprises; i.e., the M _S 8.0 Wenchuan aftershocks (Figure 4A) and the M _S 6.1 Panzhihua (Figure 4A) and M _S 6.6 Jinggu (Figure 4B) earthquakes, which did not correspond to previously mapped faults.

As shown in Figure 6, ASS_PMF < ASS_PPE; therefore, the PPE model is more informative than the PMF model. Meanwhile, ASS_PMF < ASS_PPE ≈ ASS_PKS ≈ ASS_PKSW ≈ ASS_PKS(μ=0.1), which implies that, when we specify a low value of μ close to zero, the hybrid PKS model outperforms the PMF model in any combination but does not obtain a better forecast efficiency than the PPE model.

In addition, we obtained the performance for the two methods of combining the PMF model with the PPE model. A comparison of the results reveals that, according to the ASS value, there is no obvious difference between the two ensemble methods. However, the combination method in which μ is manually specified has a more informative gain when τ is in an appropriate range, as shown in Figure 6.

Predicting the time, location, and magnitude of future seismic events is possible using information from earthquakes in historical and instrumental catalogs and geologically mapped active faults, as proven in multiple studies (Rhoades and Stirling, 2012; Rhoades et al., 2015). In the present study, we applied the PPE and PMF models to fit the seismicity in the Sichuan-Yunnan region in China. Using the optimized parameters, we estimated the seismicity rates of the PPE and PMF models at three different times: May 20, 2008; January 1, 2013; and September 10, 2019. Retrospective synthetic testing has consistently shown that both independent and hybrid models can provide informative forecasts (Rhoades and Gerstenberger, 2009; Bayona et al., 2022). Next, we introduced the PKS model proposed by Rhoades and Stirling (2012) to fit the seismicity and to calculate the seismicity rate for the same dates using the PPE and PMF models. Finally, we collected the earthquakes with M _S 6.0+ in the test region during three time periods, May 20, 2008–May 20, 2011, January 1, 2013–January 1, 2016, and September 10, 2019–September 10, 2022, as target events for prospective testing. The Molchan error diagram method was used to evaluate the model performances with respect to earthquake forecasts. Based on our results, our main conclusions are as follows.

(1) The PPE, PMF, and PKS models show better forecast efficiency with respect to moderate earthquakes compared to a homogeneous Poisson distribution in the Sichuan-Yunnan test region. While the PPE model is straightforward, its earthquake forecast ability is significant. In some cases, the PMF model has a more informative probability gain than the PPE model when τ is below 0.3, which indicates that, in the Sichuan-Yunnan region, the earthquake forecast is more effective when assigning a relatively higher threshold. We obtained a poor performance with the PMF model when τ was high.

The retrospective tests in our study also had limitations; for instance, the earthquake forecast strategy was originally designed for a specified 3-year interval. The performance may change somewhat when applied to other periods because of the degree of time dependence in the data used to generate the models. In addition, we did not account for location uncertainties when estimating the seismicity rates of the three models. Because we calculated their seismicity rates on a rather coarse grid of 0.1° × 0.1° cells, location errors can naturally be accommodated without significantly affecting the correlation with the earthquake locations (Rastin et al., 2022).

The Sichuan-Yunnan region is one of the most active seismic regions on the Chinese mainland, and various spatial parameters carrying different types of information have been used to analyze the seismic activity. The PPE model, which is straightforward, has good forecast efficiency and can be used as a base model to evaluate the effectiveness of various earthquake forecast methods. Moreover, the results of our study indicate that an approach to earthquake forecasting that is model-driven and hyper-parameter controlled could be a promising means to implement operational earthquake forecasting.