A collection of related earthquake events is referred to as an earthquake sequence. When considering the temporal occurrence of the earthquakes, the sequence can be treated as a stochastic marked point process in time. The event times t_i and magnitudes m_i are organized in a set S = (t_i, m_i), i = 1, 2, …, n [38, 39]. If an assumption is made that the magnitudes are not correlated and that t_i can be fully described by the time-dependent rate λ(t), earthquake occurrences can be modeled as a non-homogeneous Poisson process in time [27, 40, 41]. The model parameters describing the aftershock rate can be estimated.

The empirical MOL describes a decay in seismicity rate following a large mainshock [6, 24]. The functional form of the rate model is

where t is the time since the occurrence of the mainshock. The model parameters for the MOL are ω = {K_o, c_o, p_o}. K_o describes the aftershock productivity rate of the sequence with respect to the mainshock. The seismicity rate decays inversely as controlled by the decay parameter p_o. c_o is a characteristic time indicating the transition from the constant to decaying rate.

The ETAS model represents seismicity by means of a trigger model, where each parent event in the sequence can produce its own Omori-like aftershocks as in Equation (1). The functional form of the ETAS model is given by [17]:

where the parameters for the ETAS model are ω = {μ, K, c, p, α} and the rate is dependent on the history H_t of events during the time interval [T₀, t]. μ is a constant background rate during the sequence that can be modeled using a homogeneous Poisson process, and is associated with seismicity from tectonic loading [42]. The model parameters {K, c, p, α} describe the aftershock decay for short-term triggering effects in the ETAS model [43]. c is the time delay for each parent event triggering in an Omori-like aftershock sequence, and p describes the rate of decay for each local subsequence. K and α are productivity parameters associated with each individual event in the sequence. The model parameters describing the aftershock rate can be estimated using a maximum likelihood optimization method [44].

The earthquake rates, Equations (1) and (2), of the two models differ significantly. The MOL model, defined by the deterministic rate (1), constitutes a non-homogeneous Poisson process. Hence, the events occur independently with a decreasing rate. However, the ETAS rate (2) is stochastic in nature with increments for every random event that occur during the evolution of the sequence. This constitutes the self-exciting nature of this doubly stochastic contagious process with probability generating function as in [21]. The ETAS model is not a simple non-homogeneous Poisson process compared to the rate given by Equation (1).

The Gutenberg-Richter law can be represented by the left-truncated exponential distribution, where N is the number of events with magnitude m₀ or larger [45]:

The distribution of magnitudes follows the exponential distribution with the following probability density and cumulative distribution functions

and

respectively [46]. The model parameter for the GR law is described by θ = {β}, where β is related to the b-value by β = ln(10)b.

The model parameters are estimated during the training time period as in [36]. The training time period [T₀, T_s] consists of two components, the preparatory time interval [T₀, T_s] and the target time interval [T_s, T_e]. For the MOL, the preparatory time interval is a short time interval such that the mainshock is not included in the model parameter fitting procedure. For the ETAS model, the preparatory time interval is used to condition the model parameter estimates and allows for the inclusion of foreshocks as additional information. To perform real-time forecasting, the training time interval length is defined relative to the time of the mainshock. The end of the training time interval T_e can be aligned with the number of days following the mainshock, ΔT_m (Figure 1).

The two competing aftershock models can be fitted to the aftershock sequence with increasing ΔT_m using the maximum likelihood estimate (MLE) method to estimate the aftershock rate model parameters ω, for both the ETAS model and MOL, and the GR law parameter θ [47–49]. The MLE produces point parameter estimates based on the training time interval by maximizing the log-likelihood function. The likelihood is the probability of observing the aftershocks within the training time period given the model [44]. For any time dependent point process model with a well-defined intensity function (e.g., Equations 1 and 2), the likelihood function is

where fθ(m)=dFθ(m)dm is the probability density function and Λω(T)=∫TsTeλω(t|Ht)dt is the productivity of the process during the time interval [T_s, T_e] having N_e number of events above a specified magnitude. Its derivation can be found in [44, 50].

The log-likelihood is maximized using all events within the target time interval above m₀ and above the maximum depth for the sequence. For binned magnitudes, as is the case with earthquake catalog data, θ can be solved using the MLE method for magnitudes exceeding m₀ [47–49].

Using the parameter point estimates from the MLE as prior information, MCMC sampling is applied to the same training time intervals to produce model parameter distributions. The sampling procedure accounts for the uncertainty in the model parameter estimates and uses Bayesian methods to incorporate prior knowledge from the MLE estimates [36, 37, 51].

The MCMC method in this study uses the Metropolis-Hastings algorithm for the MOL and Metropolis-within-Gibbs sampling algorithm for the ETAS model to sample the posterior distributions of the model parameters. The detailed description of the implemented algorithms are given in [37]. The posterior distribution is described by

p(θ, ω|S) provides constraints on the model parameter variability by using prior information π(θ, ω) and information obtained from the training data via the likelihood function L(S|θ, ω).

In this study, we use the Gamma distribution as the prior distribution for the model parameters [37]. This choice of the priors is dictated by the well-defined functional form of the Gamma distribution that allows to limit the parameter values in certain ranges which are inferred from past studies. The Gamma distribution is described by its mean and variance, where θ and ω from the MLE are used as the mean for the Gamma priors. Use of other priors was discussed in [37]. The proposal distributions used in the MCMC sampling are set so that the parameters are chosen to approximate the posterior distribution [37]. These are the multi-variate lognormal distribution for both models. The first 100,000 steps are discarded for “burn-in” and 100,000 steps are sampled for the parameter estimates so the distribution only represents the latter 100,000 steps. Each of the resulting 100,000 model parameter estimates are then used for the forward simulations. The full details of the implementation of the Metropolis-within-Gibbs algorithm can be found in [37].

To properly estimate the aftershock rate and model parameters, the impacts of early aftershock incompleteness need to be considered [17, 52–56]. If many of these smaller events are missing from the catalog when training, this may result in bias during the subsequent aftershock rate analysis. We identify the magnitude completeness of the sequence by plotting the first 3 days of the sequence on a log-log plot using visual inspection. We then select a magnitude cutoff, m₀, where all events ≥ m₀ are included in the analyses. It is also important not to select too high a value for m₀ as doing so reduces the amount of data available for training. In particular, the ETAS model tends to underestimate the number of events during a forecast when selecting higher m₀ and the productivity from small events is not accounted for [57].

The forecasting time interval, [T_e, T_e + ΔT], immediately follows the training time interval. The parameter estimates from the training time interval are used to simulate and forecast events during the forecasting time interval. In this work, the length of the forecasting time interval was set to 7 days, ΔT = 7 days (Figure 1). The simulation of the MOL and ETAS models during the forecasting time interval is performed using the thinning algorithm [58]. For the ETAS model all observed earthquakes prior to the forecasting time intervals are used to properly define the rate during the forecasting time interval.

The point parameter estimates for the MOL and the GR law model parameter from the MLE can be directly used to calculate the probability of the largest events during the forecasting time period by means of the extreme value theory (EVD method) [37]. The probability of a large aftershock during ΔT is P_EV. This method is equivalent to the Reasenberg-Jones method [59]. For the productivity Λω(ΔT)=∫TeTe+ΔTλω(t) dt, P_EV can be calculated as [37]:

where m_ex is the magnitude of the events that we are forecasting for. To incorporate the uncertainty in the model parameter estimates, we use the Bayesian predictive distribution (BPD) [36, 37, 60]:

To compute the BPD, MCMC sampling of the posterior distribution p(θ, ω|S) of the model parameters are generated during the training time interval and are used to forward simulate the sequence in time as an ensemble during the forecasting time interval, ΔT. The simulations are used to compute the probability of large aftershocks during ΔT.

The forecast generates earthquake rates in bins for a range of magnitudes and time intervals [35]. The number of predicted events in each bin is calculated from the rate. The expectations can then be produced by assuming the probability of observing events in each bin follows the Poisson model. The joint log-likelihood is the sum of all the logs of the likelihoods in each bin when comparing the expectation to the observed events in each bin as formulated by [9]. This forms the basis of the likelihood tests.

In this work, the forecasts tests are consider the number test (N-test) and magnitude test (M-test) to evaluate these aspects individually. Two comparative tests, the ratio test (R-test) and T-test are applied to assess the relative performance of the forecasts [10, 34]. These tests are based on the assumption that the events are Poisson distributed and use a simple Poisson based likelihood function. This differs from the definition of the likelihood used in the formulation of the ETAS model, Equation (6) [19]. As a result, this may introduce confounding effects when applied to forecasts produced by the ETAS model. In addition, the Bayesian p-test is used to evaluate the performance of the BPD [36]. The first four were adapted from tests used by the Collaboratory for the Study of Earthquake Predictability (CSEP) and follow the same assumptions where the magnitudes are binned and the number of earthquakes in each forecast bin are Poisson-distributed and independent of each other [9, 35, 61]. The effective significance level of 5% for all the tests are used following [62]. Quantile scores exceeding the significance level are considered a passing score on the test. However, the applicability of the CSEP framework to the ETAS model can be considered approximate, because the ETAS model deviates from the Poisson assumption placed on the occurrence of events in the CSEP framework [19–23].

The N-test compares the number of events simulated during the forecasted time period to the number of observed events during the same time period [9]. There is an associated quantile score considered for the N-test, δ. The quantile score indicates if the generated sequences produced forecasted event numbers N_i, i = 1, ..., N_sim, above or below the observed values N_obs. In this work, the estimation of the quantile is performed through simulations of the respective point processes [9, 35]

where N_sim is the total number of simulations and I(x) is the indicator function. A value of δ close to 1 indicates that the forecast underpredicts the observed number of events, and a small δ indicates that the forecast overpredicts the observed number of events. Thus, if the probability of δ is smaller than the half of the significance level α_q/2 or larger than 1 − α_q/2, the forecast can be rejected [9]. The significance level for the respective N-tests is set to 5% (α_q = 0.05) [9].

The M-test evaluates the distribution of the magnitudes of the forecasted events during the forecasting time period compared to the true magnitude distribution of the observed events [9, 35]. For this test, the magnitude distribution is conditioned on the number of forecasted events by considering the forecasts such that the number of simulated events N_i = N_obs. The testing is done by computing the joint log-likelihood score of the simulated magnitude events. The M-test statistic is indicated by the κ as shown in

where M_obs is the joint log-likelihood of the observed events and M_iis the joint log-likelihood computed for every simulation i = 1, ..., N_sim [9].

The R-test compares two different forecasting models relative to each other by assuming one of the models is true and stands in as the null hypothesis, which is referred to as the reference hypothesis H₂. The competing model hypothesis is H₁. Each of the models has a likelihood score from the L-test based on the modified observed events assuming the null hypothesis is true (H₂ corresponds with likelihood score L₂ and H₁ corresponds with likelihood score L₁). The L-test formulation and details can be found in [35, 36]. The R-test provides the ratio of the likelihoods for two models R₂₁, where R₂₁ = L₂ − L₁ [35]. If R₂₁ > 0, then the reference model H₂ performs better. Similarly, a larger quantile score for the R-test, the quantile score α, indicates that H₂ performs better than H₁. The definition of the quantile score α is given in [10]. This test is reversible, and the null hypotheses can be set as the other model. If both models when set as the reference hypotheses result in a quantile score α larger than the effective significance level, then neither model can be rejected using the R-test [10].

To provide additional context for the R-test and to reduce computational time, the T-test, inspired by Student's t-test, was introduced as a proxy for the R-test [10]. The results of the T-test indicate whether a competing hypothesis has performed better than or is inconsistent with the reference hypothesis by evaluating the significance of the information gain. The T-test score results in the sample information gain of one model over another, where the sample information gain from the reference hypothesis (H₂) over the alternate hypothesis (H₁) is shown as IG(H₂, H₁) for N_obs, the number of observed events during the forecasting time interval ΔT [10]:

Positive values demonstrate information gain, where information gain IG > 1 can be considered significant. Negative values on the T-test indicate that the alternative hypothesis performs better than the reference hypothesis. As with the R-test, the T-test can also be reversed by changing the reference and alternative hypotheses to investigate the information gain in the other direction.

For forecasts that use the BPD method to generate probabilities, the Bayesian p-test (herein p-test) can be conducted. The value p_B gives the probability that the largest event in the simulations for the forecasted sequence will be more extreme than the observed largest event during the forecasting time interval. p_B is described as [36]:

The p-test has a test quantity (t, θ, ω) for the observed variable y and the simulated quantity ŷ. In this case, maxy is the largest event in the simulation, and p_B is the proportion of the test quantities from the simulated maximum events that are greater than or equal to the observed largest event during the forecasting time interval. Extreme values of p_B ≈ 0 indicate that the features are not well demonstrated by the model [36].

The evaluated sequences include the 2009 L'Aquila (LAQ), Italy; 2010 Darfield (DFL), New Zealand; 2016 Amatrice-Norcia (AMA), Italy; 2016 Kumamoto (KUM), Japan; and 2020 Monte Cristo Range (MCR), United States of America. The mainshock magnitudes and dates are presented in Table 1. The spatial distributions of the epicenters of earthquakes for each sequence are given in Supplementary Figures 1–5 as indicated by the dashed ellipses. A brief description of the sequence conditions are as follows. The 2009 L'Aquila (LAQ) sequence began in the central Apennines next to large and normal faults [64]. In the months leading up to the LAQ sequence, there was increased background seismicity [65]. The 2016 Amatrice-Norcia (AMA) sequence took place in central Italy, which activated a 60 km normal fault system [66]. This sequence was evaluated independently of the LAQ sequence, though it is worth noting the aftershocks from the LAQ sequence appeared to migrate toward the AMA sequence region. The catalogs for the LAQ and AMA sequence were retrieved from the Italian Seismological Instrumental and Parametric Database.

The 2016 Kumamoto (KUM) sequence began with a foreshock sequence before the mainshock of Mw 7.3. The events of the sequence indicated a right lateral strike-slip fault along the Futagawa fault segment [67]. The catalog for this sequence was retrieved from the Japan Meteorological Agency catalog. The Monte-Cristo (MCR) sequence took place along a previously unmapped, steeply dipping fault near Walker Lane, Nevada. The aftershocks took place with a mix of left-lateral, right-lateral, and normal fault motions [68]. The catalog for the MCR sequence was retrieved from the United States Geological Survey Advanced National Seismic System Comprehensive Earthquake Catalog. The 2010 Darfield (DFL) sequence took place on previously unmapped faults, in association with strike-slip movement with a reverse component for the mainshock [69, 70]. The catalog for the DFL sequence was retrieved from GeoNet.

The forecasts that are being evaluated are generated using the BPD method after an initial model parameter estimation using the MLE, which are used as prior distributions for the model parameters. Here we use the LAQ sequence as an example. The model parameter estimates using the MLE are provided in Supplementary Figures 10–14 for all of the five sequences and training time intervals. Reviewing the MLE estimates, it can be observed that the model parameters tend to stabilize and converge with increased ΔT_m, generally stabilizing after ΔT_m = 3 days for all of the sequences and aftershock models (Supplementary Figures 10–14).

The forecast can be produced for both models using the BPD method by MCMC sampling of the posterior distributions and forward simulating the models during the forecasting time interval. Plots indicating the probability of the largest expected event ≥ m_ex during the forecasting time interval ΔT = 7 days are shown for each training time interval in Figure 2 for the 2009 L'Aquila sequence and in Supplementary Figures 6–9 for the rest of the sequences. The end of each training time interval T_e corresponds with the number of days following the mainshock ΔT_m. We considered ΔT_m = [1, 2, 3, 4, 5, 6, 7, 10, 14, 21, 30] days for the training time intervals. Overall, the probability for the largest expected events decreases as the training time interval increases. Slight increases in probability are aligned with changes in the model parameter estimates for rate and decay parameters. For example, in the LAQ sequence, the model parameter K for the ETAS model jumps from 0.74 to 16.05 and p decreases from 3.2 to 2.2. There is a clear increase in probability for this training time interval and the increased rate outweighs the decrease in productivity α. The probability for an event with magnitude ≥ m_ex increases following ΔT_m = 3 days. This increase is accompanied by a decrease in decay rate (increase in p) and increase in productivity parameter α. The observed seismicity rate can be seen in the density of the black open diamonds in Figure 2. The spike in the MOL forecast for the LAQ sequence follows a slight delay that is reflected in the adjustment of the decay rate p_o.

It can be observed that the increased probability in the ETAS model aligns well with the local cluster of high density observed events following ΔT_m = 3 days, whereas the MOL peak probability for large events occurs after this cluster. This is expected, as MOL does not account for potential clustering during the forecast and incorporates the increased seismicity only after it takes place. In addition, when considering all of the sequences and their forecasts, the probability of large aftershock increases more for the ETAS model than the MOL following a large aftershock that is associated with an increase in the seismicity rate during the training time interval. The probability of large aftershocks decreases following the mainshock in general.

The N-test applied to the sequences indicates that the number of events is typically overestimated for all training time intervals for the ETAS model as indicated by the low δ scores for the DFL and AMA, and MCR sequences (Figure 3). The number of events for the MOL forecasts alternates between over- and underestimation of the observed number of events. The plots of the forecasted and observed number of events for all five sequences are given in Supplementary Figures 15–19. From the N-test plots, the effective significance score is not consistently passed on both sides for most training time intervals for both aftershock rate models. The overestimation of the MOL can be directly related with the model parameter K_o, which is larger when the rate is overestimated, or with the p_o estimate, when it is small and the rate of aftershock decay is slow (Supplementary Figures 10–14). The ETAS model overestimation cannot be directly related to a single model parameter. For most of the sequences, the background seismicity was very low and could not be the primary contributor to the overestimation from the ETAS model with the exception of the KUM sequence which had relatively high background seismicity estimates.

The MOL aftershock rate is not affected by the GR law. However, for the ETAS model, it is important to consider whether the events are impacted by the magnitude distribution from which the model draws from. For the M-test, it can be observed that when the number of events forecasted is scaled to the observed events, the magnitude distribution of the forecast closely represents the observed magnitude distribution (Figure 4). This is not surprising as this is isolated from the aftershock rate models. The impact of the GR law can be interpreted to be minimal on the number of forecasted events for the ETAS model. The MOL aftershock rate and number of events during the forecast are not impacted by the GR law. From the M-test scores, the ETAS model performance can be considered reliable and independent of the impact of the GR law fitting.

It can be observed that the performance of the MOL on the M-test typically results in a larger quantile score than the ETAS model. However, it does mimic the same trends as observed for the ETAS model. For the MOL, the training time interval is slightly longer. This suggests that the magnitude distribution immediately following the mainshock strongly influences the performance of the M-test and does not represent the magnitude distribution of events well. It follows that the b-value may be better estimated when the effect of the mainshock and in the immediate surrounding are minimized.

Both models perform well on the Bayesian p-test, except for the DFL sequence, suggesting that both models sufficiently captured the probability of the largest observed event in the forecast, so neither was rejected (Figure 5). This confirms the procedure is suitable for providing a forecast.

When comparing the performance of the two aftershock models directly, the R-test demonstrates that it is possible for both sides of the test to be passed or failed regardless of the reference model [10]. However, we find that the MOL is more likely to be rejected as an alternative hypothesis when the ETAS model is considered the reference model (Figure 6). This suggests that the MOL model tends to perform better than the ETAS model for the selected training time intervals and sequences. This is supported by the results of the T-test. The information gain of the MOL when the ETAS model is the reference model is typically negative, indicating a loss of information (Figure 7).

The 2009 L'Aquila sequence had a large aftershock of M = 5.4 2 days following the mainshock. The ETAS model reflected the increase in seismicity and corresponding large aftershock probability before the occurrence of the large aftershock (Figure 2). The ETAS model starts to overestimate the number of events after ΔT_m = 7 days (Supplementary Figure 15B). The MOL underestimates the number of forecasted events earlier in the sequence (Supplementary Figure 15A). This corresponds with a decline on the M-test performance for both aftershock rate models, which suggests that the ETAS underestimation was influenced by the magnitude distribution, while the MOL independently underestimated the number of forecasted events. In contrast to [71], the ETAS model performed comparably to the MOL based on the T-test for most of the training time intervals.

The 2010 Darfield sequence was less seismically active than would be expected based on the New Zealand aftershock decay model [72]. Previous work has shown that the ETAS model performs better than the regional models using the 3 month and 1 day forecast of the sequence [73]. On the N-test, the number of events forecasted during the simulations is overestimated by the ETAS model for the various training time intervals (Supplementary Figure 16B). The MOL fluctuates between over- and underestimation of the observed events during the forecast simulations, though the number of forecasted events is closer to the observed events for longer training time intervals ΔT_m ≥ 10 days (Supplementary Figure 16A). The performance on this test also impacts the results of the p-test. For extremely high values on the p-test, where p_B ≈ 1, we can observe that the number of events on the N-test was overestimated at these time intervals. On the M-test, the magnitude distributions of the two aftershock models paralleled each other and performed well, with quantile scores above the effective significance level (Figure 4).

The 2016 Kumamoto sequence is an example of a foreshock sequence with a typical aftershock sequence following the mainshock. The ETAS model was fitted in such a way to include the foreshock sequence during the model parameter estimations. The MOL was limited to data following the mainshock. Both models stabilized around the same training time interval length. However, the MOL underestimated the number of forecasted events during the first four training time intervals following the mainshock (Supplementary Figure 17A). This is unsurprising as decay of the foreshock sequence may have continued to contribute to the seismicity and could not be explicitly accounted for using the MOL. The ETAS model tended to slightly overestimate the number of events, which partly corresponds to the high background seismicity estimates (Supplementary Figure 17B). The M-test was passed by both models for all training time intervals. Comparison of the two models on the T-test indicate that the MOL model consistently demonstrates small information gain over the ETAS. We note that the model parameters for the KUM sequence show a distinct shift in the estimated background seismicity around ΔT_m = 7 days, which may contribute to the observed features (Supplementary Figure 12). Despite the estimation of a high background rate, the performance on the N-test indicated that the ETAS model alternated between over- and under-estimation of the number of events during the forecasting time interval. When evaluating the forecasting scores for ΔT_m ≥ 10 days, then the ETAS model performs slightly better than the MOL.

The 2016 Norcia sequence is an example of a classic foreshock–aftershock sequence. The background seismicity of the sequence was considered elevated due to the occurrence of the Mw 6.0 Amatrice mainshock on August 24, 2016, and set to μ = 1.0. It can be observed that the MOL model parameters are stable and estimated consistently during all of the training time intervals (Supplementary Figure 13A). The ETAS model parameters stabilize following ΔT_m = 4 days (Supplementary Figure 13B). The ETAS model forecast tended to overestimate the number of events in the forecasting time interval, though not greatly (Supplementary Figure 18B). The MOL tended to underestimate the number of events during the forecasting time interval (Supplementary Figure 18A). In the comparison tests, the MOL model consistently demonstrated information gain over the ETAS and the R-test was failed by the MOL in several instances which agrees with [66].

The 2020 Monte Cristo Range sequence demonstrates an example of multiple ETAS model parameters being constrained to produce consistent estimates. The MOL parameters converged after ΔT_m = 3 days, with the K_o parameter remaining consistent (Supplementary Figure 14). The performance on the N- and M-test are similar, though the MOL performs better on the N-test by more frequently scoring a pass on both sides of the test. The ETAS model typically overestimates the number of forecasted events (Supplementary Figure 19B). The MOL model demonstrates marginal information gain over the ETAS model. For this sequence, while the MOL model appears to perform better than the ETAS model, it is not fully clear if the improvement is significant.