**1. 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].

**Citation:** Sedghizadeh, M.; Shcherbakov, R. The Analysis of the Aftershock Sequences of the Recent Mainshocks in Alaska. *Appl. Sci.* **2022**, *12*, 1809. https://doi.org/ 10.3390/app12041809

Academic Editor: Amadeo Benavent-Climent

Received: 19 December 2021 Accepted: 5 February 2022 Published: 10 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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–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 magnitudes 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.

**Table 1.** The dates of occurrence, epicentre locations, magnitude and depth of the analyzed mainshocks.


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.

**Figure 1.** *Cont*.

**Figure 1.** Maps of the occurrence of the aftershock sequences generated by the three significant Alaska mainshocks: (**a**) the 2002 Mw 7.9 Denali sequence with *m*0 = 3.0; (**b**) the 2018, Mw 7.9 Kodiak sequence with *m*0 = 3.2; and (**c**) the 2018 Mw 7.1 Anchorage sequence with *m*0 = 2.8. The events during 30 days after each mainshock are plotted. The blue solid circles represent the aftershocks above *m*0. Black points are all events between magnitude 2.5 and *m*0. The focal mechanisms of the studied mainshocks are plotted as beach balls. Quaternary faults are plotted as light brown line segments [33].

#### **2. Materials and Methods**

## *2.1. 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–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].

**Table 2.** The parameters of the elliptical regions used for the identification of the aftershock sequence and the corresponding lower magnitude cutoffs *m*0.


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, [*<sup>T</sup>*0, *Te*], 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, [*<sup>T</sup>*0, *Ts*], and a target time interval, [*Ts*, *Te*]. The seismicity parameters are estimated in the target time interval. A forecasting time interval, [*Te*, *Te* + Δ*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.


**Figure 2.** An illustration of the time intervals used in the analysis.

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, *mi*, and the occurrence time, *ti*, in order to generate a stochastic marked point process during a specific time interval, *S* = {(*ti*, *mi*)} : *i* = 1, 2, . . . , *n*.

#### *2.2. 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:

$$N(m \ge) = 10^{(a-bm)},\tag{1}$$

where *<sup>N</sup>*(*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*):

$$f\_{\theta}(m) = \beta e^{-\beta(m - m\_0)},\tag{2}$$

$$F\_{\theta}(m) = 1 - e^{-\beta(m - m\_0)},\tag{3}$$

where *m*0, is the lower magnitude cut-off that is above the catalogue completeness magnitude *mc* and *θ* = {*β*} is the model parameter which can be obtained from all earthquakes above *m*0 in target time interval [*Ts*, *Te*]. The parameter *β* is related to the *b*-value of the Gutenberg–Richter (GR) scaling:

> *β*

$$
\delta = b \ln(10). \tag{4}
$$

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.

## *2.3. Omori–Utsu Law*

For the first time, Omori [43] introduced a formula for the aftershock sequence decay rate, *<sup>λ</sup>*(*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:

$$
\lambda\_{\omega}(t) = \frac{K}{(t+c)^p},\tag{5}
$$

where *λω*(*t*) is the earthquake rate at a given time *t* with magnitudes above *m*0, and set of parameters *ω* = {*<sup>K</sup>*, *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 *ω* = {*<sup>K</sup>*, *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.

#### *2.4. 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]:

$$\lambda\_{\omega}(t|\mathcal{H}\_t) = \mu + A \sum\_{i:t\_i < t}^{N\_t} \frac{e^{a(m\_i - m\_0)}}{(\frac{t - t\_i}{c} + 1)^p},\tag{6}$$

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*<sup>t</sup>*, during the time interval [*<sup>T</sup>*0, *t*]. *Nt* is the number of the earthquakes with magnitudes above *m*0 in the time interval [*<sup>T</sup>*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:

$$\log L = \sum\_{i:t\_i \le T\_\varepsilon} \lambda\_{\omega}(t\_i | \mathcal{H}\_{t\_i}) - \int\_{T\_\varepsilon}^{T\_\varepsilon} \lambda\_{\omega}(t | \mathcal{H}\_t) \, dt. \tag{7}$$

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 *ti* as

$$
\pi\_i = \int\_0^{t\_i} \lambda\_\omega(t) \, dt. \tag{8}
$$

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.

#### *2.5. 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, [*Te*, *Te* + Δ*T*] can be obtained from the Extreme Value Distribution (EVD) [47]:

$$\Pr\_{EV}\{m\_{cx} \ge m | \theta, \omega, \Delta T\} = 1 - e^{-\left\{\Lambda\_{\omega}(\Delta T) \left[1 - F\_{\theta}(m)\right]\right\}},\tag{9}$$

where *mex* 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 <sup>Λ</sup>*ω*(<sup>Δ</sup>*t*) is a productivity function that is given as:

$$
\Lambda\_{\omega}(\Delta T) = \int\_{T\_{\rm f}}^{T\_{\rm f} + \Delta T} \lambda\_{\omega}(t) \, dt. \tag{10}
$$

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:

$$\begin{split} \Pr\_{EV} \{ m\_{\varepsilon x} \ge m \mid \theta, \omega, \Delta T \} &= \\ 1 - \exp \left\{ - \left[ K \frac{(T\_{\varepsilon} + c)^{1-p} - (T\_{\varepsilon} + \Delta T + c)^{1-p}}{p-1} \right] (\varepsilon^{-\beta(m-m\_0)}) \right\}, \quad \text{(11)} \end{split} $$

for *p* = 1, and the set of parameters {*<sup>θ</sup>*, *ω*} can be obtained during the target time interval [*Ts*, *Te*]. Therefore from Equation (11), the probability of having an earthquake with a magnitude above *m* in a forecast time interval [*Te*, *Te* + Δ*T*] can be obtained, which is the same approach as in Reasenberg and Jones [8].

#### *2.6. 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:

$$\Pr\_B\{m\_{\varepsilon x} \ge m \mid S\_\prime \Delta t\} = \int\_{\Omega} \int\_{\Theta} \Pr\_{EV}(m\_{\varepsilon x} \ge m \mid \theta\_\prime \omega \omega \Delta T) p(\theta\_\prime \omega \mid S) \, d\theta d\omega,\tag{12}$$

where Θ and Ω are the frequency-magnitude distribution and seismicity rate parameter domains, respectively. Pr*EV*(*mex* ≥ *m* | *θ*, *ω*, Δ *T*) is the EVD and *p*(*<sup>θ</sup>*, *ω* | *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 [*Te*, *Te* + Δ *<sup>T</sup>*]. 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.

## *2.7. 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 Λ<sup>2</sup> over the model Λ<sup>1</sup> is defined as *IN*(Λ2, Λ<sup>1</sup>) = *<sup>R</sup>*21/*<sup>N</sup>*obs, where *N*obs is the number of observed earthquakes during the forecasting time interval Δ *T* and *R*<sup>21</sup> = *L*(M|Λ<sup>2</sup>) − *L*(M|Λ<sup>1</sup>) 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].
