The Analysis of the Aftershock Sequences of the Recent Mainshocks in Alaska

: The forecasting of the evolution of natural hazards is an important and critical problem in natural sciences and engineering. Earthquake forecasting is one such example and is a difﬁcult task due to the complexity of the occurrence of earthquakes. Since earthquake forecasting is typically based on the seismic history of a given region, the analysis of the past seismicity plays a critical role in modern statistical seismology. In this respect, the recent three signiﬁcant 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) presented an opportunity to analyze these sequences in detail. This included the modelling of the frequency-magnitude statistics of the corresponding aftershock sequences. In addition, the aftershock occurrence rates were modelled using the Omori–Utsu (OU) law and the Epidemic Type Aftershock Sequence (ETAS) model. For each sequence, the calculation of the probability to have the largest expected aftershock during a given forecasting time interval was performed using both the extreme value theory and the Bayesian predictive framework. For the Bayesian approach, the Markov Chain Monte Carlo (MCMC) sampling of the posterior distribution was performed to generate the chains of the model parameters. These MCMC chains were used to simulate the models forward in time to compute the predictive distributions. The calculation of the probabilities to have the largest expected aftershock to be above a certain magnitude after a mainshock using the Bayesian predictive framework fully takes into account the uncertainties of the model parameters. Moreover, in order to investigate the credibility of the obtained forecasts, several statistical tests were conducted to compare the performance of the earthquake rate models based on the OU formula and the ETAS model. The results indicate that the Bayesian approach combined with the ETAS model produced more robust results than the standard approach based on the extreme value distribution and the OU law.


Introduction
The Pacific Ring of Fire is one of the most seismically active regions of the world. Alaska and western Canada are a part of this ring and are prone to the occurrence of significant earthquakes. This geographic region is characterized by high seismic activity and is capable of producing megathrust earthquakes. These earthquakes can pose significant hazard and are also capable of triggering tsunamis or intense ground shaking [1] and subsidiary hazards such as liquefaction, landslides and aftershocks [2]. While tsunamis pose a serious threat to coastal areas, ground shaking can cause damage to infrastructure and endanger human life. Therefore, it is important to perform a comprehensive statistical analysis of the aftershock sequences in the Aleutian subduction zone and central Alaska. Moreover, the occurrence of large aftershocks poses a significant risk to the infrastructure that has been affected by a mainshock. Therefore, estimating the probabilities for the occurrence of the largest expected aftershocks plays an important role in post-mainshock decision-making [3,4].
One of the earliest empirical studies of the difference between the magnitude of the mainshock and its largest aftershock was conducted by Båth [5], who postulated that the largest aftershock is on average 1.2 magnitudes lower than the mainshock regardless of the magnitude of the mainshock. Vere-Jones [6,7] proposed that the magnitude difference between the mainshock and the largest aftershock was independent of the number of events. Reasenberg and Jones [8] were one of the first in developing an aftershock forecasting model. They introduced a parametric model that was capable of computing the probabilities of aftershocks in a certain time window after a mainshock for California. Michael et al. [9] proposed the methodology which aforementioned model parameters can be estimated with Bayesian updating from both the ongoing aftershock sequence and from historic aftershock sequences.
An important step in the calculation of the probability of having an earthquake above a certain magnitude is the estimation of the model parameters that describe the seismicity rate and the frequency-magnitude distribution. Those parameters are highly dependent on the lower cut-off magnitude, m 0 . The correct estimation of the cut-off magnitude plays a crucial role in earthquake forecasting and modelling. Mignan and Woessner [10] emphasized that a high-value of the cut-off magnitude can result in under-sampling of useful data and a low-value of the cut-off magnitude can result in uncertainty and bias of the estimated seismicity parameters and forecasting model, respectively.
The other issue in aftershock forecast modelling is the catalogue incompleteness right after the occurrence of strong mainshocks [11,12]. This early catalogue incompleteness can affect significantly the estimation of the parameters of the earthquake decay rate. The uncertainties in the estimation of the parameters of the aftershock decay rate can result in significant miscalculation of the probabilities for the occurrence of largest events. The empirical prior probability distribution was presented by Omi et al. [13] to reduce the uncertainty of the parameter estimation of the ETAS model regarding the incompleteness of the earthquake catalogues.
Utilizing generic parameters to create an aftershock forecast model for the early days after the mainshocks is one of the possible ways to control the catalogue incompleteness. Page et al. [14] introduced a method for generic parameter estimation by using tectonic zoning of García et al. [15] to improve the spatial distribution of forecasted events. In this approach, Bayes' rule and aftershock records are used to update the generic parameters. In addition, the distribution of the regional generic parameters can be considered as a prior and the aftershock data can be used to calculate the posterior distribution. Michael et al. [9] applied this approach to the 2018 Anchorage aftershock sequence. They reported that the use of the generic parameters for the forecast model leads to the overestimation of the seismic activity.
One of the critical tasks in statistical seismology is the ability to accurately and reliably forecast the evolution of earthquake sequences. A consistent approach for earthquake forecast testing has been implemented in the Collaboratory for the Study of Earthquake Predictability (CSEP) [16][17][18][19]. In this framework, the gridded rate forecast is used in which the selected geographic area is separated into zones then the number of earthquakes in each zone is estimated [19]. In addition, the number of earthquakes in each forecast bin is considered to be independent of the other bins and follows the Poisson distribution. Several statistical tests were developed as part of the CSEP framework to examine earthquake forecasts. As a result of these developments, it is possible to determine if a particular forecasting scheme is able to accurately replicate locations, magnitudes, and the observed numbers of earthquakes [19,20]. Various forecasting algorithms can also be compared using the aforementioned likelihood-based tests. For example, retrospective aftershock forecasting of the 2011 Tohoku, Japan; 2010 Canterbury, New Zealand; 2016 Kaikoura, New Zealand; and 2019 Ridgecrest, California earthquakes were tested by using this approach [4,[21][22][23][24].
In this study, the analysis of three major earthquake sequences that occurred in Alaska in the past 20 years was conducted to test retrospectively the ability to forecast the magni-tudes of the largest expected aftershocks. Specifically, the 23 January 2018 Mw 7.9 Kodiak earthquake occurred in the Gulf of Alaska near Kodiak Island at 09: 31:40.89 UTC at a depth of 14 km [25,26]. There was no significant damage reported. The earthquake woke residents in Anchorage which was located 560 km northeast from the epicenter. It was also felt in parts of British Columbia, Canada. The 2018 Mw 7.1 Anchorage earthquake happened approximately 15 km north from Anchorage, Alaska on 30 November of 2018 at 17:29:29.33 UTC at a depth of 46.7 km [9,27]. A few minutes later, a magnitude 5.8 aftershock shook the region. Significant damage has been reported to infrastructure, buildings, and airports [9]. Moreover, we investigated the characteristics of the 3 November 2002, Mw 7.9 Denali earthquake that occurred in central Alaska along a shallow strike-slip fault on the Denali-Totschunda fault system [28,29]. The details of the selected mainshocks are listed in Table 1 and the spatial distributions of the mainshocks with the corresponding aftershock sequences during the first 14 days are shown in Figure 1. In this study, the left truncated exponential distribution was utilized to model the magnitude frequency statistics [30]. Moreover, the modified Omori-Utsu (OU) law [31] and Epidemic Type Aftershock Sequence (ETAS) model [32] were used to approximate the rate of the aftershocks. In addition, two statistical approaches including the extreme value distribution and Bayesian predictive distribution were utilized to compute the probabilities of having the largest expected aftershocks to be above a certain magnitude during the evolution of each sequence.
The paper is structured as follows: Section 2 begins with the specification of the earthquake catalogues and follows by defining the statistical methods to analyze the aftershock sequences. The results of the statistical analysis are provided in Section 3. In Section 4, discussion of the results and concluding remarks are given.

Earthquake Catalogue
In order to analyze the 2002, Mw 7.9 Denali; the 2018, Mw 7.9 Kodiak, and the 2018, Mw 7.1 Anchorage earthquake sequences, the United States Geological Survey (USGS) earthquake catalogue https://earthquake.usgs.gov/earthquakes/search/ (accessed on 18 December 2021) was used. The spatial distribution of aftershocks during 30 days after each mainshock are shown in Figure 1. The focal mechanisms of the mainshocks were obtained from the USGS website [34][35][36].
The 2002 Denali, Alaska, earthquake sequence occurred along the Denali-Totschunda faults which is a right-lateral strike-slip fault system. The Mw 7.9 mainshock nucleated on the Susitna Glacier thrust fault and propagated further along the Denali fault and continued along the Totschunda fault [29]. The parameters of the elliptical region for this aftershock sequence, are given in Table 2 and the sequence for 30-day is depicted in Figure 1a. The plotted fault plane solution for this mainshock in Figure 1a was taken from the the USGS website [34].
The 2018 Mw 7.9 Kodiak, Alaska, earthquake took place in the Gulf of Alaska southeast of Kodiak Island. The mainshock location and the focal mechanism reflect a strike-slip faulting system within the shallow lithosphere of the Pacific plate near the subduction zone [37]. In Table 2 the details of the studied sequence for this mainshock are reported. In addition, the earthquake sequence during 30-day after the mainshock occurrence is demonstrated in Figure 1b. The focal mechanism of the 2018 Mw 7.9 Kodiak, Alaska, earthquake suggests a steeply dipping fault either as a right-lateral system that strikes the north-northwest or as a left-lateral fault that strikes west-southwest. The fault plane solution for the mainshock in Figure 1b was obtained from the USGS website [35].
The 2018 Mw 7.1 Anchorage, Alaska, earthquake happened as the result of a normal faulting rupturing north from Anchorage, Alaska. The location and mechanism of the focal mechanism reflect a moderately dipping north-south fault system fault within the subducting Pacific slab [38,39]. Details of the analyzed earthquake sequence are presented in Table 2 and the sequence of 30-day is illustrated in Figure 1c. The indicated moment tensor for this mainshock in Figure 1c was acquired from USGS [36]. For the statistical analysis of seismicity, several time intervals were utilized to estimate properly the parameters of the models describing the evolution and the statistics of the aftershock sequences. For the estimation of the model parameters, the training time interval, [T 0 , T e ], is considered. In order to properly account for the impact of preceding earthquakes on the earthquake rate, the training time interval is divided into an initial time interval, [T 0 , T s ], and a target time interval, [T s , T e ]. The seismicity parameters are estimated in the target time interval. A forecasting time interval, [T e , T e + ∆T] is also considered to analyze the evolution and the statistics of the seismicity. The schematic illustration of the time intervals for the analysis of the aftershock sequences is shown in Figure 2.

Mainshock
Training time interval Forecasting time interval Earthquakes occur due to sudden energy release associated with the slippage of faults and are characterized by finite rupture areas. However, for statistical analysis of seismicity, the point assumption is utilized to characterize each earthquake. On time scales larger than the propagation of rupture along the fault, earthquakes can be treated as points in time and space. This idealization helps to describe the earthquake process by point process models. The point process becomes a marked point process by assigning magnitudes to each event. Therefore, each earthquake can be characterized by the magnitude, m i , and the occurrence time, t i , in order to generate a stochastic marked point process during a specific time interval, S = {(t i , m i )} : i = 1, 2, . . . , n.

Gutenberg-Richter Scaling and the Exponential Distribution
Gutenberg and Richter [40] proposed the relationship that describes the frequencymagnitude statistics of earthquakes. This relationship between the frequency of event occurrences and the event magnitudes is one of the most commonly used empirical laws in statistical seismology. They suggested the following equation: where N(m ≥) is the total number of earthquakes above magnitude m and N(0 ≥) = 10 a and b is the value of the slope of the fitted line to the N on a logarithmic scale. Vere-Jones [30] emphasized that the distribution of earthquake magnitudes is described by the exponential distribution for m ≥ m 0 with the probability density, f θ (m), and cumulative distribution function, F θ (m): where m 0 , is the lower magnitude cut-off that is above the catalogue completeness magnitude m c and θ = {β} is the model parameter which can be obtained from all earthquakes above m 0 in target time interval [T s , T e ]. The parameter β is related to the b-value of the Gutenberg-Richter (GR) scaling: The Maximum Likelihood Estimation (MLE) is the most common approach to estimate the b-value or parameter β. Bender [41] suggested an estimator for β by taking into account the binning of the magnitude. Tinti and Mulargia [42] proposed an approach to calculate the uncertainties of the parameter β at a given confidence level.

Omori-Utsu Law
For the first time, Omori [43] introduced a formula for the aftershock sequence decay rate, λ(t), that is inversely related to the elapsed time after the mainshock. Utsu [31] proposed a modification of the Omori law which is known as the Omori-Utsu (OU) law. Utsu [31] modified the original intensity to the following form: where λ ω (t) is the earthquake rate at a given time t with magnitudes above m 0 , and set of parameters ω = {K, c, p}, and t is the time elapsed since the occurrence of the mainshock at T 0 = 0. The parameter K is the productivity of the sequence, c is the characteristic time, and p is the rate of the decay in time. By considering the non-homogeneous Poisson process for the occurrence of earthquakes, the parameters ω = {K, c, p} can be determined by using the MLE approach [44,45]. In addition, in this model, the parameter uncertainties can be estimated from the inverse of the Fisher information matrix that is computed from the likelihood function.

The Epidemic Type Aftershock Sequence (ETAS) Model
A more realistic approximation of the earthquake rate was proposed by Ogata [46], where he suggested that each earthquake could be considered as a trigger for the next events in the sequence. The conditional intensity of the temporal ETAS model, λ ω (t|H t ), at time t is defined as [46]: where ω = {µ, A, c, p, α} is the set of parameters of temporal conditional intensity with a reference magnitude m 0 and the occurrence history of earthquakes, H t , during the time interval [T 0 , t]. N t is the number of the earthquakes with magnitudes above m 0 in the time interval [T 0 , t]. In the ETAS model, µ specifies the average rate of background events that transpire independently of any other events. c is the temporal characteristic time, p governs the rate of decay of triggered events as a power law, and A controls the event productivity. The parameter α determines the degree of aftershock clustering. Larger values of α correspond to more pronounced aftershock sequences with stronger variability in earthquake magnitudes. In contrast, the impact of event's magnitude on aftershock generation is reduced by smaller α values. The estimation of the ETAS model parameters is achieved by maximizing the log-likelihood function: In general, the consistency of the ETAS model is measured on the basis of a transformed time. The transformed time τ i for a given event is computed by using the cumulative conditional intensity at time t i as If the fit of the model is accurate, the sequence of earthquakes should obey a stationary Poisson process in the transformed time. Furthermore, the cumulative number of observed earthquakes in transformed time can be close to a straight line [13]. The deviation of the cumulative number of observed events from the straight line indicates that the model does not fit well the earthquake sequence.

Extreme Value Distribution
By considering a non-homogeneous Poisson sequence of earthquakes, the probability of having an extreme earthquake with a magnitude above m in the forecasting time interval, [T e , T e + ∆T] can be obtained from the Extreme Value Distribution (EVD) [47]: where m ex is the magnitude of the largest expected event, ω is the set of parameters of seismicity rate λ ω (t), F θ (m) is the cumulative distribution function of the events' magnitude with the set of parameters θ, and Λ ω (∆t) is a productivity function that is given as: By considering the exponential model, Equation (3), for describing the magnitude distribution and the UO model, Equation (5), for the intensity of the productivity function, Equation (9) can be rewritten as: for p = 1, and the set of parameters {θ, ω} can be obtained during the target time interval [T s , T e ]. Therefore from Equation (11), the probability of having an earthquake with a magnitude above m in a forecast time interval [T e , T e + ∆T] can be obtained, which is the same approach as in Reasenberg and Jones [8].

Bayesian Predictive Distribution
The obtained parameters of the aftershock sequence model during the training time interval play a crucial role in calculating the EVD. The uncertainty of the parameters have a significant impact on the calculation of the corresponding probabilities. Shcherbakov et al. [3,48] incorporated the model uncertainties into the computation of the probabilities for the occurrence of the largest expected earthquakes by applying the Bayesian predictive distribution (BPD) approach, in which the BPD can be defined as: where Θ and Ω are the frequency-magnitude distribution and seismicity rate parameter domains, respectively. Pr EV (m ex ≥ m | θ, ω, ∆T) is the EVD and p(θ, ω | S) is the posterior distribution function, which quantifies the uncertainties of the model parameters. Since the ETAS model deviates from a non-homogeneous Poisson process the EVD for the largest magnitudes is not given by Equation (9). Shcherbakov et al. [3] suggested to use the stochastic simulations to approximate the extreme value distribution and ultimately the BPD. In this approach, the Metropolis-within-Gibbs algorithm is used to sample from the conditional posterior distribution to generate the chain of the model parameters using the Markov Chain Monte Carlo (MCMC), then the model parameter chain is used to simulate the ETAS model during the forecasting time interval [T e , T e + ∆T]. At the end, the maximum magnitude is taken from each sequence of events to construct a distribution that approximates the BPD.
When performing MCMC sampling a certain initial part of the parameter chain is discarded as "burn-in". The Gamma distribution was considered for the prior distribution of the model parameters. As burn-in the first 50% of Markov chains were discarded and the second half was utilized for calculation of the BPD.

Forecast Validation
To evaluate the number of forecasted earthquakes by a specific model in the forecasting time interval, the N-test can be used [4,17,19,49]. It tests the distribution range of the number of the forecasted events versus the number of observed earthquakes. In addition, in order to test the magnitude distribution of the forecasted earthquakes the M-test can be applied [4,17,19,49]. The N and M-tests examine the consistency of the forecasts with respect to observations, and the R-test can be used to compare the performance of different forecasting models [17]. In addition, to evaluate statistical forecast the T-test can be applied [49]. In formulating the T-test, the sample information gain per earthquake of the model Λ 2 over the model Λ 1 is defined as I N (Λ 2 , Λ 1 ) = R 21 /N obs , where N obs is the number of observed earthquakes during the forecasting time interval ∆T and R 21 = L(M|Λ 2 ) − L(M|Λ 1 ) is the log-likelihood ratio of the two models. The detailed explanation and implementation of these tests applied to the time-dependent models such as the ETAS and OU rates can be found in Shcherbakov [4].

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.

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 m c 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].

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.  Figure 4 and for the Denali and Kodiak sequences in Figure A2.   Figure 4.
Furthermore, we used the ETAS model, Equation (6), to estimate the aftershock decay rate during the target time interval [T s , T e ] = [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.  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 T e = [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

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 [T s , T e ] = [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 T e = [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.       (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 [T s , T e ] = [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.

Testing the Model Forecasts
Several tests were conducted to evaluate the forecast during the time interval [T e , T e + ∆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 T e = [1,2,3,4,5,6,7,10,14,21,30] days are given in Figure 10.  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. 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  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.   (c)

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.  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 [T e , T e + ∆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. 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][35][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.

Abbreviations
The following abbreviations are used in this manuscript:         (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). The aftershock magnitudes are modelled using Equation (3).