**3. Results**

In this section the obtained results of the analysis of the recent three significant mainshocks that occurred in Alaska (the 2002, Mw 7.9 Denali; the 2018, Mw 7.9 Kodiak, and the 2018; Mw 7.1 Anchorage earthquakes) are summarized.

#### *3.1. Frequency-Magnitude Statistics Analysis*

The aftershocks of the three mainshocks within elliptical regions as shown in Figure 1 were used to obtain the frequency-magnitude statistics. The fitting of Equation (2) was done to all three sequences and during specific target time intervals. To estimate the parameter *β* from Equation (2) the MLE approach was used [41]. The model parameter uncertainties were estimated using the method of Tinti and Mulargia [42]. In addition, the method of goodness of fit test [50] was utilized to estimate the magnitude of completeness *mc* for the

three sequences. Specifically, *m*0 was selected as the magnitude above which at least 95% of the observed data are modeled by Equation (2). The results are presented in Figure 3. Moreover, to investigate how the earthquake magnitudes evolve over time, they are plotted versus the sequential number in Figure A1. This can be used to inspect the early magnitude incompleteness of aftershock sequences and can be used to justify the use of a chosen magnitude threshold [51].

**Figure 3.** The frequency-magnitude statistics of aftershock sequences for 30 days from the mainshock occurrence. The lines are the Gutenberg–Richter scaling fit, Equation (1). The open symbols represent the cumulative numbers corresponding to each aftershock sequence. The estimated *a* and *b*-value with 95% confidence intervals are given in the legend. The cumulative numbers of aftershocks for the 2002, Denali aftershocks for *m* ≥ 3.0 are plotted as red circles, the 2018, Kodiak aftershocks for *m* ≥ 3.2 cumulative numbers are plotted as purple diamonds, and blue squares are used to depict the 2018, Anchorage aftershocks for *m* ≥ 2.8.

In order to analyze the frequency-magnitude statistics of the 2002, Denali earthquake sequence, *m*0 = 3.0 was considered as a cut-off magnitude, and the analysis was performed during [*Ts*, *Te*]=[0, 30] days after the mainshock occurrence on 3 November 2002 (22:12:41 UTC) for the earthquakes within the elliptical region given in Figure 1a. The total number of aftershocks during the selected time interval was 771 with the maximum magnitude 5.6, respectively. The fit of the GR relation is demonstrated in Figure 3. The estimated *b*-value and *a*-value for the analyzed earthquake sequence are 0.749 ± 0.053 and 5.134 ± 0.176, respectively. The magnitude-frequency statistics analysis of the 2018, Kodiak earthquake sequence was performed during [*Ts*, *Te*]=[0, 30] days after the mainshock that occurred on 23 January 2018 (09:31:4 UTC) for earthquakes with the cut-off magnitude *m*0 = 3.2 within an elliptical region shown in Figure 1b. In total 1207 earthquakes with the magnitude ranging from 3.2 to 5.5 occurred in the analyzed elliptical region during the specified time interval. The estimated *b*-value and *a*-value for selected earthquake sequence are 1.035 ± 0.059 and 6.394 ± 0.214, respectively (Figure 3). For analyzing the frequencymagnitude statistics of the 2018, Anchorage earthquake sequence, *m*0 = 2.8 was considered as a cut-off magnitude, and the analysis was carried out during [*Ts*, *Te*]=[0, 30] days after the mainshock occurrence on 30 November 2018 (17:29:29 UTC) for the earthquakes within an elliptical region of Figure 1c. In total 476 earthquakes within the magnitude range of 2.8–5.2 occurred in the analyzed spatiotemporal window. The obtained *b*-value and *a*-value for Anchorage earthquake sequence are 0.906 ± 0.082 and 5.216 ± 0.246 respectively. The fit of the GR relation is plotted in Figure 3.

#### *3.2. Aftershock Decay Rate Modelling*

To begin with, we obtained the first-order approximation to the background seismicity rate within the presented elliptical regions in Figure 1 for each analyzed aftershock sequence. In this evaluation the background rate was estimated as the ratio of the number of earthquakes to the number of days, during the time interval that started on 1 January 2000, and ended 30 days before each mainshock. The estimated background rates, the corresponding time intervals, and cut-off magnitude for each analyzed mainshock are reported in Table 3.

**Table 3.** Background Seismicity Rate.


Subsequently, the aftershock decay rate was modeled by using the OU law, Equation (5), for the three sequences during the target time interval of [*Ts*, *Te*]=[0.001, 30] days. The obtained parameters of the OU model with 95% confidence interval and the model fits are shown for the Anchorage sequence in Figure 4 and for the Denali and Kodiak sequences in Figure A2.

**Figure 4.** The log-log plot of the aftershock decay rate for the 2018, Mw 7.1 Anchorage aftershock sequence with magnitudes *m* ≥ 2.8 are presented as open squares. The blue solid line is the corresponding fit of the OU law, Equation (5), to the aftershock sequence. The obtained parameters from the OU law, Equation (5), with the 95% confidence intervals are reported in the legend.

The estimated set of parameters of the OU law, Equation (5), with 95% confidence intervals for the 2002, Denali earthquake sequence are *ω* = {*K* = 202.72 ± 48.35, *c* = 0.29 ± 0.13, *p* = 1.22 ± 0.11} (Figure A2a). Similarly for the 2018, Kodiak earthquake sequence the obtained parameters of the OU law, Equation (5) with 95% confidence intervals are *ω* = {*K* = 225.91 ± 46.84, *c* = 0.31 ± 0.18, *p* = 0.88 ± 0.09} (Figure A2b). Finally, for the 2018, Anchorage earthquake sequence the estimated parameters of the OU law are *ω* = {*<sup>K</sup>*, *c*, *p*}, {76.71 ± 11.25, 0.06 ± 0.03, 1.17 ± 0.1}, respectively. The fit of the OU law and the estimated parameters from Equation (5) are demonstrated in Figure 4.

Furthermore, we used the ETAS model, Equation (6), to estimate the aftershock decay rate during the target time interval [*Ts*, *Te*]=[0.06, 30] days for the three aftershock sequences. For comparison we present the OU fit and the ETAS model fit for the 2018, Anchorage earthquake sequence in Figure 5. A similar plot for the other two sequences is

given in Figure A3. Figure 6 illustrates the fit of the ETAS model in transformed time for the 2018, Anchorage earthquake sequence.

**Figure 5.** The aftershock sequence and corresponding earthquake magnitudes during 2018, Mw 7.1 Anchorage sequence with *m* ≥ 2.8. The ETAS model fit, Equation (6), for the target time interval of [*Ts*, *Te*]=[0.06, 30] is plotted as a solid blue line, and the obtained set of parameters are reported with 95% confidence intervals. The OU law fit, Equation (5), is plotted as a black dashed line for comparison.

**Figure 6.** The cumulative number of observed aftershocks in transformed time and corresponding rate of the ETAS model for the 2018, Mw 7.1 Anchorage aftershock sequence with *m* ≥ 2.8. The ETAS model fit, Equation (6), for the target time interval of [*Ts*, *Te*]=[0.06, 30] days is plotted as a solid blue line, and the obtained set of parameters are reported with 95% confidence intervals.

Finally, the statistical properties of the aftershock sequence initiated by the M7.1 Anchorage, Alaska, earthquake are investigated in detail during several additional target time intervals. Specifically, the sequence was analyzed during several target time intervals starting from the occurrence of the mainshock and ending at *Te* = [1, 2, 3, 4, 5, 6, 7, 10, 14, 21, 30] days. The evolution of the estimated parameters with 95% confidence intervals for both models OU and ETAS during 2018, Mw 7.1 Anchorage earthquake sequences are shown in Figure A4. Obtained estimations for the *b*-value of the GR relation, Equation (1) with 95% confidence intervals are demonstrated in Figure A4a. The evolution of the OU model parameters are shown in Figure A4b. We presented the evolution of the estimation of the ETAS model parameters in Figure A4c.

#### *3.3. Retrospective Forecasting of the Largest Expected Aftershocks*

In order to calculate the probability of having the magnitude of the largest expected aftershock to be above a certain magnitude and during a predefined forecasting time interval the EVD, Equation (9), and BPD, Equation (12), are used. In this analysis, the OU law, Equation (5), and the ETAS model, Equation (6), are utilized to calculate the aftershock decay rate, and the frequency-magnitude distribution is estimated from the exponential distribution, Equation (3).

To illustrate the applicability of the methods, one particular example is illustrated in case of the 2018 Anchorage sequence. The training time interval was set to [*Ts*, *Te*]=[0.06, 14] days and the forecasting time interval of Δ *T* = 7 days was considered. The lower magnitude cut-off *m*0 = 2.8 was used. The computed distributions using the EVD, Equation (9), and BPD, Equation (12) are plotted in Figure 7. For the BPD analysis, total of 20,000 MCMC sampling of the posterior distribution was performed. The first 10,000 iterations were discarded as "burn-in" and the remaining 10,000 samples were utilized to perform stochastic simulations of the ETAS or OU processes. The resulting distributions of the OU and ETAS model parameters estimated from the MCMC chains are reported in Table A1 and plotted in Figures A5 and A6, respectively.

Moreover, two cases are considered for computing the probabilities for the occurrence of the largest expected aftershocks above a certain magnitude during the evolution of the three sequences. For the first case, we considered a constant forecasting time interval Δ *T* = 7 days. As for the target time intervals, we considered the following ending times *Te* = [1, 2, 3, 4, 5, 6, 7, 10, 14, 21, 30] days with the lower magnitude thresholds of 3.0, 3.2, and 2.8 for the 2002, Denali, 2018, Kodiak, and 2018, Anchorage sequences, respectively. In this analysis, the BPD, Equation (12), with the exponential distribution, Equation (3), for the frequency magnitude statistics, and the ETAS model, Equation (6), for the occurrence rate of the earthquakes were utilized. The obtained results for the probabilities of the largest expected earthquakes greater than *m*ex ≥ 5.0, 6.0, 7.0 are shown in Figure 8. Furthermore, the probabilities of having the largest expected aftershock to be above magnitude 6.0 were computed utilizing the EVD, Equation (8), combined with the OU law, Equation (5), and using the BPD, Equation (12), combined with the OU law, Equation (5), or the ETAS model, Equation (6). The obtained results for analyzed mainshocks are presented in Figure 9.

**Figure 7.** The EVD and BPD for the 2018, Mw 7.1 Anchorage aftershock sequence during the 7-day forecast time interval after the training time interval of 14 days. The blue solid line represents the BPD using the ETAS model with 10000 MCMC sampling steps using the Gamma prior. The orange line represents the obtained BPD using the OU model and the yellow line is the plot of the EVD with the OU law, Equation (11).

**Figure 8.** The probabilities to have the largest expected aftershocks to be larger than *m*ex ≥ 5.0, 6.0, 7.0 using the BPD, Equation (12), during a constant forecasting time interval Δ*T* = 7 days and for the varying target time intervals. (**a**) The 3 November 2002, Mw 7.9 Denali sequence with *m* ≥ 3.0. (**b**) The 23 January 2018, Mw 7.9 Kodiak sequence with *m* ≥ 3.2. (**c**) The 30 November 2018, Mw 7.1 Anchorage sequence with *m* ≥ 2.8.

**Figure 9.** The comparison of the probabilities to have the largest expected aftershock during the forecasting time interval Δ*T* = 7 days for the three sequences: (**a**) the 3 November 2002, Mw 7.9 Denali for *m* ≥ 3.0; (**b**) the 23 January 2018, Mw 7.9 Kodiak for *m* ≥ 3.2; (**c**) the 30 November 2018, Mw 7.1 Anchorage for *m* ≥ 2.8. The blue triangles are computed using the BPD, Equation (12), with an earthquake decay rate given by the ETAS model, Equation (6). The purple squares are computed using BPD, Equation (12), with an earthquake decay rate given by OU law, Equation (5). The green circles give probabilities computed using the EVD, Equation (11).

In the second case, a constant target time interval [*Ts*, *Te*]=[0.06, 2] days was considered. However, the forecasting time interval was varied as Δ *T* = [1, 2, 5, 7, 10, 14] days to compute the probabilities of the occurrence of the largest expected aftershocks. The computed probabilities for the largest anticipated earthquakes *m*ex ≥ 5.0, 6.0, 7.0 are illustrated in Figure A7. In addition, the comparison of the two approaches to compute the probabilities (EVD versus BPD) combined with either the OU law or the ETAS model are shown in Figure A8.

#### *3.4. Testing the Model Forecasts*

Several tests were conducted to evaluate the forecast during the time interval [*Te*, *Te* + Δ *T*] by comparing the simulated results with the observed seismicity. To check the performance of forecasts for the number of aftershocks and magnitude distribution, the N and M-tests were performed, respectively. The details of the implementation of the tests can be found in Shcherbakov [4].

For the three aftershock sequences the number of forecasted aftershocks in the forecasting time interval Δ *T* = 7 days and using the following target time intervals *Te* = [1, 2, 3, 4, 5, 6, 7, 10, 14, 21, 30] days are given in Figure 10. For comparison, in the same figure the observed number of earthquakes are shown as blue circles for each prediction time interval Δ *T*. In addition, for the constant target time interval [*Ts*, *Te*]=[0.06, 2] days and varying forecasting time intervals Δ *T* = [1, 2, 5, 7, 10, 14] days the number of forecasted and observed earthquakes are shown in Figure A9. To investigate the effect of the magnitude cutoff, we also performed the same analysis for earthquakes above magnitude 3.5. This is reported in Figure A10.

**Figure 10.** *Cont*.

**Figure 10.** The number of forecasted and observed aftershocks in the forecasting time interval Δ*T* = 7 days for the three sequences: (**a**) the 3 November 2002, Mw 7.9 Denali sequence for *m* ≥ 3.0; (**b**) the 23 January 2018, Mw 7.9 Kodiak sequence for *m* ≥ 3.2; (**c**) the 30 November 2018, Mw 7.1 Anchorage sequence for *m* ≥ 2.8. The red triangles show the average number of forecasted earthquakes using the ETAS model and the black squares illustrate the average number of forecasted earthquakes using OU law. The shading bands represent 95% confidence intervals. The blue circles represent the observed number of earthquakes in each forecasting time interval.

In addition, M-test was performed to assess the consistency of the distribution of the magnitudes of the forecasted events. The results of the performance of the OU law and ETAS model are reported by computing the quantile score, *κ* [4,19]. *κ* is defined as the proportion of the forecasted magnitudes compared to the observed magnitudes in each magnitude bin. The obtained quantile scores for the constant forecasting time interval Δ*T* = 7 days are given in Figure 11. In addition, in Figure A11 the outcomes for the M-test for the constant target time interval [*Ts*, *Te*]=[0.06, 2] days are plotted for the three aftershock sequences.

**Figure 11.** *Cont*.

**Figure 11.** The obtained quantile scores from the M-test for the constant forecasting time interval Δ*T* = 7 days for the three sequences: (**a**) the 3 November 2002, Mw 7.9 Denali for *m* ≥ 3.0; (**b**) the 23 January 2018, Mw 7.9 Kodiak for *m* ≥ 3.2; (**c**) the 30 November 2018, Mw 7.1 Anchorage for *m* ≥ 2.8. The red triangles demonstrate the obtained quantile scores from the ETAS model and the black squares illustrate the quantile scores of OU law. The blue dashed lines represent the 0.025th and 0.05th quantiles.

In order to evaluate and compare the models, the R-test and T-test were applied for both cases by considering the ETAS model versus the OU law. In the R-test the quantile score, *α*, was calculated. *α* is the proportion of the simulated likelihood ratios, over the observed likelihood ratios [17]. The values of *α* that are greater than a specific level of significance support the model that was chosen as a base model, in this case it is the ETAS model. The obtained result for the *α* from the OU law versus the ETAS model is shown in Figure 12. Furthermore, in Figure 13, the ratio of the likelihood score of the ETAS model and the OU law over the number of the observed events in the forecasting time interval is given. The T-test is used to assess whether the sample information gain is statistically different from zero. This is used to select the preferred model [49]. Lastly, the Bayesian *p*-value of the BPD analysis using either the ETAS model or the OU law are illustrated for both cases in Figure 14.

**Figure 12.** The obtained quantile scores from the R-test to compare the forecast based on the OU law versus the ETAS model for the three sequences: (**a**) the 3 November 2002, Mw 7.9 Denali for *m* ≥ 3.0; (**b**) the 23 January 2018, Mw 7.9 Kodiak for *m* ≥ 3.2; (**c**) the 30 November 2018, Mw 7.1 Anchorage for *m* ≥ 2.8. The black squares illustrate the quantile scores in the case with the constant forecasting time interval of Δ*T* = 7 days. The red triangles show the obtained quantile scores from the second case with the constant target time interval of [*Ts*, *Te*]=[0.06, 2] days. The blue dashed lines represent the 0.025th and 0.05th quantiles.

**Figure 13.** The sample information gain of the ETAS model versus the OU law over the number of observed data for the three sequences: (**a**) the 3 November 2002, Mw 7.9 Denali for *m* ≥ 3.0; (**b**) the 23 January 2018, Mw 7.9 Kodiak for *m* ≥ 3.2; (**c**) the 30 November 2018, Mw 7.1 Anchorage for *m* ≥ 2.8. The black solid squares illustrate the sample information gain for the first case with the fixed the forecasting time interval of Δ*T* = 7 days. The red triangles demonstrate the sample information gain for the second case for the constant target time interval [*Ts*, *Te*]=[0.06, 2] days.

**Figure 14.** For both cases the Bayesian *p*-value of the BPD from the ETAS model and the OU law are illustrated for the three sequences: (**a**) the 3 November 2002, Mw 7.9 Denali for *m* ≥ 3.0; (**b**) the 23 January 2018, Mw 7.9 Kodiak for *m* ≥ 3.2; (**c**) the 30 November 2018, Mw 7.1 Anchorage for *m* ≥ 2.8. The green squares and the red triangles illustrate the obtained *p*-value from the first case, with the fixed forecasting time interval Δ*T* = 7 days. The yellow diamonds and purple circles demonstrate the *p*-value for the second case with the constant target time interval of [*Ts*, *Te*]=[0.06, 2] days. The blue dashed lines represent the 0.025th and 0.05th quantiles.

#### **4. Discussion and Conclusions**

To describe the three aftershock sequences which occurred in the Alaska region, statistical models were used in this study. To be more precise, the frequency-magnitude statistical analysis was performed using the GR relation, and the occurrence rates of the aftershock sequence were estimated by the OU law and the ETAS model. The EVD and BPD approaches were used to calculate the probability of having the largest expected aftershock above a certain magnitude during evolution of each sequence for various training and forecasting time intervals.

The frequency-magnitude distributions and estimated GR parameters for analyzed sequences are shown in Figure 3. The frequency-magnitude distribution of the Denali aftershock sequence, was characterized by a broad distribution of event magnitudes which led to a relatively low *b*-value (0.749 ± 0.053) (Figure 3). The 2002, Denali mainshock was followed by several large aftershocks with the largest being 5.6 magnitude which occurred in the first 24 h after the mainshock. The frequency-magnitude statistical analysis of the aftershock sequence of the Kodiak 7.9 magnitude earthquake with the cut-off magnitude 3.2 for the target time interval [*Ts*, *Te*]=[0, 30] days was performed (Figure 3). The frequencymagnitude distribution of the 2018, Kodiak sequence indicates a typical GR fit with *b*-value (1.035 ± 0.059). The mainshock was followed by several large aftershocks with the largest being 5.5 magnitude event. It should be noted the epicenter of the 2018, Kodiak earthquake was located in a remote area in the North Pacific. The analysis of the 2018 Anchorage sequence produced the *b*-value of 0.906 ± 0.082 (Figure 3). The largest aftershock with a magnitude 5.8, was close to the expected magnitude from Båth's law [5] which states that the largest aftershock is on average 1.2 magnitudes lower than the mainshock.

In order to analyze the occurrence rate of the aftershock sequence of the selected mainshocks, the Omori–Utsu law, Equation (5), and the ETAS model, Equation (6), were utilized. The obtained results from the analysis of the decay rate of the aftershock sequences show that the parameter *p* is comparable for both models (the OU law and the ETAS model) except for the 2018 Kodiak sequence.

Computing the probability of having a largest expected aftershock with a magnitude above a given value during different forecasting time intervals after the mainshock was one of the main objectives of this work. The EVD and BPD approaches were used to accomplish this objective.

The obtained result of this analysis indicates that the BPD method using the ETAS model is more conservative than BPD using the OU law and the EVD approach. In addition, for the aforementioned approaches, the probabilities of having an earthquake with magnitude 6 and above were calculated for both cases (Figures 9 and A8).

Moreover, in order to compare characteristics of analyzed sequences the probabilities to have the largest expected aftershocks to be larger than *m*ex ≥ 5.0, 6.0, 7.0 were estimated by using the BPD, Equation (12), for both cases (Figures 8 and A7). The results of this analysis show that the Anchorage sequence had a lower potential to generate aftershocks with *m*ex ≥ 5.0, 6.0, 7.0 compared to other analyzed sequences. These statistical results can be explained directly by the number of events in the aftershock sequence and the magnitude of the mainshock. This increases the probability of occurrence of an aftershock with a certain magnitude in a predefined time interval after the mainshock. In the present implementation of the EVD and BPD analysis we used the unbounded GR distribution. It was suggested that more realistic truncated magnitude distribution can be more appropriate for forecasting [52,53]. This can be easily incorporated in the analysis as well.

The N-tests, M-test, R-test, and T-test were performed to evaluate the goodness of the models' results in the forecasting time interval [*Te*, *Te* + Δ *T*] by comparison of the simulated results and observed seismicity. The number of the forecasted earthquakes was evaluated by the N-test and the results are shown in Figures 10 and A9 for both cases. In both cases, for the Anchorage sequence, a more accurate forecast for the number of earthquakes was accomplished by the ETAS model, while for two other sequences the OU law performed better. It should be noted as a result of the branching nature of the ETAS model, it shows a

wider confidence interval range compared to the OU law. The obtained results from the M-test demonstrate higher consistency in generating the distribution of the magnitudes for the ETAS model compared to the OU law for both cases, Figures 11 and A11. For the model comparison, the R-test was performed (Figure 12). The *α* quantile score is higher than the thresholds representing the rejection of the OU hypothesis in favor of the ETAS model. The T-test results are given in Figure 13. The ETAS model performed better in case of the Denali and Kodiak sequences, however, the OU model was more accurate in estimating the rate and the corresponding forecasting performance in case of the Anchorage sequence in its early days. This is also evident when plotting the information gain both for the fixed forecasting time interval Δ *T* = 7 days with varying training time intervals and in case of the fixed training time interval with varying forecasting time intervals (Figure 13c). The posterior predictive *p*-value test was performed to assess the fit of the posterior distribution of Bayesian models by comparison of the posterior predictive distribution and the observed data. In Figure 14 the results of Bayesian p-value analysis are given. They indicate that the forecasts based both on the ETAS model and OU formula are consistent in reproducing the maximum event during each corresponding forecasting time interval.

The obtained results indicate that for the sequences analyzed the forecasting based on the ETAS model and the OU formula produce comparable results for shorter time intervals after the mainshocks. However, the EAST model is more realistic in terms of reproducing the seismicity on longer time scales. Moreover, the ETAS model performs better when the mainshock sequence is preceded by a well defined foreshock sequence [3].

The ETAS model typically performs better with increased number of events in the sequence. However, this is limited by the current earthquake catalogues which typically have relatively high level of completeness that results in fewer events.

**Author Contributions:** Conceptualization, M.S. and R.S.; methodology, M.S. and R.S.; software, M.S. and R.S.; validation, M.S. and R.S.; writing—original draft preparation, M.S.; writing—review and editing, R.S.; visualization, M.S. and R.S.; supervision, R.S.; project administration, R.S.; funding acquisition, R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** M.S. and R.S. would like to acknowledge the support from the UWO IDI grant. R.S. was partially supported by NSERC Discovery grant.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The United States Geological Survey (USGS) search engine was used to extract the earthquake catalogue https://earthquake.usgs.gov/earthquakes/search/ (accessed on 18 December 2021). The United States Geological Survey (USGS), was used to obtain the Quaternary fault and fold database for the United States [33]. The mainshock focal mechanism was obtained from the USGS Moment Tensor catalogue [34–36]. The data analysis was performed using a computer code written in Matlab and can be requested from the author.

**Acknowledgments:** The authors would like to thank Matteo Taroni and two other anonymous reviewers for their useful and constructive comments that helped to improve the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.
