Research on the forecasting strategy of early aftershocks in North China

To systematically investigate short-term aftershock forecasts following moderate-to-strong earthquakes of North China so as to develop "operational" aftershock forecasting strategies tailored to regional seismic characteristics, we adopt the widely used Reseanberg-Jones (R-J) model and the Omi-R-J model to explore aftershock forecasting strategies of 24 earthquake sequences of North China, and use the N-test and T-test to evaluate the forecasting effectiveness. Early forecast results after mainshock show that the R-J model and the Omi-R-J model have an average effectiveness rate of 77.0% and 87.9% for the selected sequences, respectively. The R-J model has a lower ratio of forecasting "too low" numbers of earthquakes while the Omi-R-J has a rather low overall "failure rate". With the rapid development of monitoring network after 2008, the efficacy of earthquake sequence forecasting has improved significantly, with monitoring capability being an important factor constraining forecasting effectiveness. The possible scientific strategies for the aftershocks forecasting in North China include strengthening the construction of seismic networks and applying better algorithms for detecting or simulating aftershocks, giving priority to using the Omi-R-J model within a short period of time (within 2 hours) after the mainshock, and weigh the range of the actual number of future aftershocks by appling the R-J model and the Omi-R-J model simultaneously.


Introduction
It is crucial to release the reliable and authoritative aftershock forecast information in time for earthquake tracking, relief, emergency management, and order maintenance in disaster areas after a moderate-to-strong earthquake [Reasenberg and Jones, 1989;Gerstenberger et al., 2005;Marzocchi and Lombardi, 2009;Woessner et al., 2011;Nanjo et al., 2012;Ogata et al., 2013]. Moderate-to-strong earthquakes can cause great losses to life and property, which is in no small portion due to secondary disasters including aftershocks especially within three days after the mainshock, while the occurrence time of larger aftershocks varies from a few days to several months after the mainshock, and may even be as close as within one day after [Japan Meteorological Agency (JMA), 2009]. In the early stage after the mainshock, the drowning effect of the large surface wave amplitude of the mainshock and the superposition of signals from a large number of simultaneous aftershocks will reduce the seismic monitoring capacity in the region and even globally within several hours after the mainshock, presenting significant challenges for aftershock probability forecasting [Ogata, 1983;Utsu et al., 1995;Enescu et al., 2007;Iwata, 2008]. Therefore, the research on the scientific strategy of the aftershock forecasting in the early stage after mainshock is a realistic problem in a general sense.
In the current Collaboratory for the Study of Earthquake Predictability (CSEP) program, the short-term aftershock forecasting model has seen rapid development [Schorlemmer et al., 2018]. More than 400 earthquake forecast models have been submitted to the four international testing centers of the CSEP program for "retrospective" and "prospective" research, [Jiang et al., 2015[Jiang et al., , 2018Omi et al., 2015;Bi and Jiang, 2017;Han et al., 2017;Taroni et al., 2018;Ogata et al., 2018]. The Operational Aftershock Forecasting (OAF) conducted by the US Geological Survey (USGS) and the Global Earthquake Model (GEM) project represents an important application in earthquake forecasting strategies [Helmstetter et al., 2006;Console et al., 2010].
City clusters of North China are developing rapidly, and the coordinated development of the Beijing-Tianjin-Hebei region and the construction of Xiong'an New Area are well underway. To ensure the smooth implementation of national strategies, it is necessary to upgrade security measures, such as seismic hazard assessment, which includes aftershock probability forecasting, as well as risk reduction. It is important to rely on modern technology at home and abroad, and mature regional applications and innovative research to effectively improve the technical capabilities of aftershock forecasting and risk management decisions. Technically, by testing in North China the internationally accepted and comparable R-J model [Reasenberg and Jones, 1989] and the Omi-R-J model [Omi et al., 2013[Omi et al., , 2015[Omi et al., , 2016 capable of using incomplete seismic records, we hope to derive knowledge about the suitability of these two aftershock forecasting models, in a bid to inform aftershock forecasting strategies of North China for the benefit of actual earthquake disaster prevention and relief work.

Data and Method
Distributed across the North China region (34~43 o N, 108~125 o E) are complex seismic belts such as the Zhangjiakou-Bohai, the Hebei plain, the Shanxi, and the Tanlu seismic belts, featuring strong tectonic activity, the development of faults and folds, various types of earthquake sequences, and frequent occurrence of seismic events, causing great losses to people's lives and property. The region is prone to strong earthquakes, with historical data indicating more than 20 earthquakes of magnitude 7 or above. The 1976 Tangshan Earthquake in particular is the deadliest earthquake in the 20 th century, destroying the entire city and left 242,000 people dead and 164,000 seriously injured.
For the study on the aftershock sequences of moderate-to-strong earthquakes in North China, we adopted the "National Unified Official Catalogue" 1 provided by the China Earthquake Networks Center. We use Gardner-knopoff method [Gardner and Knopoff, 1974] to separate 10145 earthquake catalogues with completeness magnitude of 3.0 and above in North China. After deleting the aftershocks, we obtained 2629 earthquake catalogues, 159 of which are above 4.5. In order to select the earthquake sequence, we use the method of the "natural boundary method" Jiang, 2017, 2019] based on the combination of latitude-time plot, longitude-time plot and epicenter distribution map to select desired sequences. In addition, according to the calculation requirements of R-J model, a certain amount of aftershock events are needed in the earthquake sequence. Therefore, in the process of preliminary selection of earthquake sequence, we need not less than 60 earthquake events, and not less than 30 earthquake events above the completeness magnitude for parameter fitting.
In the fitting process of the R-J model, to ensure the completeness of earthquake catalogs and the inclusion of a sufficient number of aftershocks in the calculation, we also use the "Magnitude-Rank" method [Huang, 2006;Jiang and Wu, 2011;Zhuang et al., 2012] to choose completeness magnitude of different sequence. The "Magnitude-Rank" method is based on the sequence of earthquake occurrence time, and then the completeness magnitude is determined according to the distribution of magnitude and rank. In addition, in the process of parameter fitting, it is necessary to set the starting time of fitting C 0 , which is related to the cutoff magnitude M c : the larger M c , the smaller C 0 . By adjusting the relationship between the two parameters, we find out the values of M c and C 0 suitable for each seismic sequence, so as to ensure that enough earthquakes participate in the fitting and calculation.
According to the above selection rules, a total of 24 earthquake sequences in North China are selected, and their spatial distribution map is shown in Figure 1, and their mainshock parameters are shown in Table 1, together with the values of M c and C 0 . Taking the M S 5.6 Fengzhen earthquake in Inner Mongolia on August 13, 1981 as an example, Figure 2 gives a schematic diagram for the selection of this sequence. Iterating each seismic event step by step, the Omi-R-J method can obtain seismic detection rate changes over time during the study period. Figure 3 shows the distribution of the 50% seismic detection rate ( ) of the M S 5.6 Fengzhen earthquake in Inner Mongolia according to Ogata and Katsura (1993)'s model. To investigate the stability of the fitting process, detection rate distributions in 5 periods, namely, 0~0.10 days, 0~1.00 days, 0~3.00 days, 0~5.00 days, and 0~30.00 days after the mainshock, were separately investigated. Reasenberg and Jones [1989] developed the Reseanberg-Jones (R-J) model based on the "Omori-Utsu" formula [Omori, 1894;Utsu, 1961] as the intensity constraint and the G-R law [Gutenberg and Richter, 1944;Aki, 1965] as the frequency constraint for short-term forecasts of aftershocks. Therefore, the function of aftershock intensity whose magnitude is at or above M at time t in the earthquake sequence can be written as: In order to overcome the shortcomings of not fully utilizing a large number of small-magnitude earthquake events below the M c , and inability to quickly and reliably estimate model parameters in the early stage of aftershock sequences, Omi et al. [2013] developed the "Omi-R-J model" which is based on the R-J model but incorporates some new technologies. Omi et al. [2013] used the expression of the detection rate function ( ) provided by Ogata andKatsura [1993, 2006] to describe detection degrees of the incomplete parts of seismic records.
The actually recorded probability density function can be written as: ( 2) After obtaining parameters , and and the dynamic changesof detection rate by the "state-space" model, where V is a hyper-parameter that controls the degree of smoothness of ( ). The parameters , and from the Omori-Utsu (O-U) formula for aftershock decay and the parameter from the Gutenberg-Richter (G-R) formula for the magnitude-frequency relation ( = bln10) are constants to be estimated. Thereafter, the number of forecasting aftershocks in the magnitude range > and any time interval [t 2 , t e ] can be calculated using equation (1): It is essential to obtain stable model sequence parameters before investigating corresponding strategies for aftershock probability forecasting [Jiang et al., , 2018Bi and Jiang, 2017]. For the parameters , and of the R-J model, we can use the maximum likelihood method to fit the Omori-Utsu formula, and use the Fisher information matrix to estimate the standard deviation of these parameters [Ogata, 1983]. The parameter b-value and its standard deviation can be calculated using the maximum likelihood method [Aki, 1965] and the methods given by Shi and Bolt [1982], respectively.
In order to study the continuous change of the R-J model parameters, we used a series of time windows with fixed starting points but increasing length to select earthquake sequences and perform model parameter fitting. The setting of the starting time C 0 and the ending time (sequence duration time or starting time of forecasting) ₂ of these time windows is shown in Tables 1 and 2, respectively. Among them, ₂ generally increases gradually from 0.05 day after the mainshock to 5.00 days in steps of 0.05 day, but in order to ensure that the time windows contain a sufficient number of earthquakes, the forecast starting point of t 2 for each earthquake sequence is slightly different, see Table 2. Taking the Zhangbei M S 6.3 earthquake in Hebei Province on January 10, 1998 as an example, Figure   4 The parameters of the Omi-R-J model were fitted using Bayesian estimates given by Omi et al. [2013]. Specifically, the Newton's iteration and the expectation maximum (EM) algorithm were used to optimize the hyperparameters , , V and (t), respectively, while the maximum likelihood method is used to estimate the parameters , , and . Different from the traditional R-J method, the Omi-R-J method utilizes all aftershock events recorded after the earthquake, and performs parameter fitting immediately (within 2 hours) after the earthquake, requiring only a small number of aftershock events to obtain seismic sequence parameters. Taking the Zhangbei M S 6.3 earthquake in Hebei Province on January 10, 1998 as an example, we fitted the model parameters, and selected three time periods with relatively sharp changes in the early stage, namely, 0~0.10 days, 0~1.00 days, 0~2.00 days ( Figure 5). As Figure 5 shows

Statistical test of forecast results
We used the N-test [Kagan and Jackson, 1995;Schorlemmer et al., 2007;Zechar, 2010] as adopted in CSEP for earthquake number verification to determine the deviation of the forecasting number of earthquakes by the R-J model and the Omi-R-J model from the actual number of earthquakes. We use to represent the number of forecasted aftershocks and to represent the actual number of aftershocks. In the simplified processing of zechar (2010), the "too low" and "too high" in the actual number of aftershocks is tested by the scores 1 and 2 , respectively.
The effective significance level = 0.025 was used to unilaterally test the N-test results. 1 < indicates "too low" numbers of forecasted aftershocks, while 2 < indicates "too high" numbers of forecasted aftershocks.
As an example, Figure 5(d), (e) and (f) show the forecasts at three "turning points" by the Omi-R-J model for 0.10~1.10 days, 1.00~2.00 days, and 2.00~3.00 days after the earthquake. As can be seen, the number of actual earthquakes in the three time periods fall within the 95% confidence interval of the forecasts, indicating relatively high forecasting efficacy. As an example, we provide the test results of the forecasting efficacy of the R-J and Omi-R-J models in predicting the future 1-day aftershock probability of Zhangbei M S 6.3 earthquake in Hebei Province on January 10, 1998 with sequence duration time t 2 =2.00 ( Figure 6). Figure 6 Tables 2 and 3. The corresponding N-test results of the two models are shown in Figure 7, where the red color blocks represent the results when 1 <0.025 and 2 <0.025, the slash time period represents the time period that cannot participate in the calculation because of insufficient data, and the blank area indicates that no earthquake event of or above the cut-off magnitude has occurred in the forecasting period. As the N-test cannot be performed if no actual events have taken place, to ensure objective evaluations, this paper does not include situations with zero actual occurrence despite a forecasted probability of earthquakes.
For complex earthquake sequences, forecasting "failures" were observed even with the well-accepted R-J model and Omi-R-J model. As Table 2, Table 3 and Figure 7 suggest, in all time periods involved in the calculation, the R-J model is more effective than the Omi-R-J model with respect to having "too low" numbers of forecasted events, while the Omi-R-J model considerably outperforms the R-J model in respect to having "too high" numbers of forecasted events. The Omi-R-J model also shows considerably better outcomes in terms of overall failure rates.
As the fitting and calculation of the R-J model is impossible in the early stage of certain sequences, and given relatively poor results of early stage forecasts in statistical test, the actual forecasting efficacy of the R-J model may Jinmeng Bi et al. The starting times of forecasts is marked by a black vertical dashed lines, and the ending times is the right edge of the picture. Three plots in each line refer to the different test periods (0.1,1.1), (1,2) and (2,3), respectively.
be overestimated. To further compare forecast performances between the two models, we analyzed 10 sequences in the same forecasting periods of the two models. By comparison, the R-J model is more effective than the Omi-R-J model with respect to having "too low" numbers of forecasted aftershocks, while the Omi-R-J model has significantly higher overall efficacy than the R-J model. Therefore, utilizing the R-J and Omi-R-J models by focusing on the bottom line thinking of "no more than" and "no less than" the corresponding forecasting earthquake number respectively may present practical value to post-earthquake relief in North China. Detailed evaluations of the efficacy for each earthquake sequence are given in Tables 2 and 3.
To further explore contributing factors for forecasting failures, we analyzed changes in the seismic network's monitoring capacity before and after 2008. The average failure rate of the R-J model turned from 27.0% before 2008 to 7.7% after, while that of the Omi-R-J model turned from 13.3% before 2008 to 8.9% after. We conducted fitting to better show the relationship between occurrence time and failure rate ( Figure 8). As Figure 8 shows, with increased monitoring capacity over time, both models have witnessed improved forecasting capacity, which is more pronounced for the R-J model than the Omi-R-J model, also indicating smaller impact of the incomplete records of early aftershocks on the latter. Similar conclusions were made by Jiang et al. [2018] by randomly "deleting" small magnitudes. Thus, at the early stage after the mainshock when numerous aftershocks are unobserved due to dense overlapping of seismic waves, the Omi-R-J model has more practical usefulness. The development of monitoring technologies for improved monitoring capacity may be a key factor to raising aftershock forecasting efficacy.

Jinmeng Bi et al.
10   In addition, the T-test (Student's) method in CSEP was used to explore the relative merits of the two models.
The T-test method can be expressed as the calculation of the average "information gain per earthquake" (IGPE) of model A relative to the model B in the confidence interval [Imoto, 2007]: where is the total number of "target earthquakes", and ln and ln are the likelihood functions of model and model , respectively.
For calculations of the space grid, we adopted the same simplified method as Jiang et al. [2017], that is, directly using forecasting numbers of earthquakes in corresponding time windows. Based on the IGPE calculation method and the T-test, this study evaluated the forecasting efficacy of the Omi-R-J model relative to the R-J model for the sequences of 17 identical forecasting periods in North China according to their respective magnitudes. Test results at 95% confidence interval are shown in Figure 9. T-test results suggest that the Omi-R-J model significantly outperforms the R-J model for 88.2% (15/17) of the sequences.
12 Figure 9. T-test results of the Omi-R-J model against the R-J model for aftershock forecasting of the earthquake sequence in North China, with different "target magnitudes". The black dots and horizontal lines represent the information gain per earthquake and the 95% confidence interval respectively, the number of observed earthquakes is shown above the T-test results.

Conclusion and Discussion
To explore forecasting strategies for early aftershocks in North China, we have conducted a systematic comparative study of the widely used R-J model and the Omi-R-J model, which is capable of fully using large numbers of incomplete small-magnitude events in the early stage of the aftershock sequence, for 24 earthquake sequences in North China. We have arrived at the following conclusions by performing continuous sliding and fitting of multiple time windows, aftershock occurrence rate forecasting, and efficacy evaluation by the N-test and T-test: 1. Both the R-J model and Omi-R-J model prove rather effective for the early stage of sequences in the North China region, although the Omi-R-J model is generally better than R-J model, including the IGPE results given by T-test and the lower failure rate shown by N-test. Thus, it may be useful for post-disaster relief in North China to develop aftershock forecasting strategies that capitalize on the relative advantages of the R-J model and the Omi-R-J model in terms of relatively "low ratio of having 'too few' forecasted numbers of earthquakes" and "low total failure rate".
2. We analyzed changes in the seismic network's monitoring capacity before and after 2008 and we found that the average failure rate of the R-J model dropped from 27.0% before 2008 to 7.7% after; while that of the Omi-R-J model dropped from 13.3% before 2008 to 8.9% after. This decline suggests considerable improvement in forecasting efficacy for earthquake sequences with the rapid growth of the seismic network following 2008, and that monitoring capacity is a key factor constraining forecasting efficacy. Considering that seismic network in China has improved significantly around 2008, especially the increase in the number of seismic stations and the quality of instruments, further enhancement of seismic monitoring capabilities in the future is still an important factor to improve the short-term aftershock forecasting capabilities in North China.
Immediate access to relatively stable model parameters is essential to improving forecasting capacity within a short time following strong earthquakes, while the missing aftershock sequences resulting from "submerged" small earthquake waveforms in the early stage may be one major factor affecting the stability of model parameters. In replenished early aftershock sequences datasets, tentative practices have been adopted to gain access to stable parameters after earthquakes, such as using "matched filtering" to pick up more repetitive aftershock events [Peng et al., 2006], a temporal point process with time-independent marks [Zhuang and Wang, 2016;Zhuang et al., 2017], etc. The results of this paper reveal the substantial impact of seismic monitoring capabilities on aftershock forecasting performance, therefore, the development of these new technologies provides another way to increase the potential of aftershock forecasting efficiency compared to investing in more seismic stations. The Omi-R-J model, which is formed by technical improvements based on the traditional R-J model, allows a large number of incomplete smallmagnitude earthquakes to participate in aftershock forecasting, further exerting the stability of model parameters fitting and the reliability of aftershock forecasting results. The application of these technologies has practical significance for formulating more scientific and reliable aftershock forecasting strategies in North China.