Characterizing the Foreshock, Main Shock, and Aftershock Sequences of the Recent Major Earthquakes in Southern Alaska, 2016–2018

For each of three major M ≥ 7.0 earthquakes (i.e., the January 24, 2016, M7.1 earthquake 86 km E of Old Iliamna; the January 23, 2018, M7.9 earthquake 280 km SE of Kodiak; and the November 30, 2018, M7.1 earthquake 14 km NNW of Anchorage, Alaska), the study considers characterization of the foreshock and aftershock sequences in terms of their variations and scaling properties, including the behavior of the control parameter η of the unified scaling law for earthquakes (USLE), along with a detailed analysis of the surface wave records for reconstruction of the source in the approximation of the second moments of the stress glut tensor to obtain integral estimation of its length, orientation, and development over time. The three major earthquakes at 600 km around Anchorage are, in fact, very different due to apparent complexity of earthquake flow dynamics in the orogenic corner of the Pacific and North America plate boundary. The USLE generalizes the classic Gutenberg-Richter relationship taking into account the self-similar scaling of the empirical distribution of earthquake epicenters. The study confirms the existence of the long-term periods of regional stability of the USLE control parameter that are interrupted by mid- or even short-term bursts of activity associated with major catastrophic events.


INTRODUCTION
The territories around Anchorage have recently experienced four major earthquakes ( Figure 1) and associated seismic sequences which illustrate the complexity of the tectonic environment at the orogenic corner of the boundary between the Pacific and North America plates where an Eulerian plate model is not expected to be accurate (Bird, 2003). The January 24, 2016 and the November 30, 2018 earthquakes, both of M7.1, occurred under the North American plate at the south and north edges of the Cook Inlet, respectively. The largest of the four, the January 23, 2018, M7.9 earthquakes ruptured the Pacific plate in front of the continental crust of Alaska. The July 22, 2020, M7.8 Alaska Peninsula earthquake (July 22, 2020 06:12:44.719 UTC, 55.0683°N and 158.5543°W, 28 km, M ww 7.83) ruptured a 200-km segment of the Aleutian megathrust fault in the study area after the manuscript submittal. All four appear to rupture the subducting Pacific plate right at the border of or within the extended source region of the March 27, 1964 Great Alaska, M9.3 mega-earthquake (Press and Jackson, 1965;Wyss and Brune, 1967;Kanamori, 1970;Christensen and Beck, 1994). The complexity of the megathrust is characterized with a multiple rupture of several segments of subducting Pacific plate, including the lateral transition faulting along the Yakutat block at the corner of the Pacific-North America plate boundary. The apparent reactivation of this region at the level of the M ≥ 7 earthquakes deserves special attention of seismologists. Therefore, in the following sections, we provide integral characterization of the fore-and aftershock sequences for each of the three earthquakes in terms of their magnitude-space-time distributions and the control parameter of the unified scaling law for earthquakes (Kossobokov and Mazhkenov, 1994;Bak et al., 2002;Kossobokov and Nekrasova, 2019), as well as the average estimates of the rupture extent, duration, and velocity, making use of the low-degree moments of the stress glut rate (Backus, 1977a;Backus, 1977b). Nowadays, when large earthquakes occur in sensible areas, including the Alaska's most populous city of Anchorage, numerous studies provide source tomographies (e.g., Ohta et al., 2006;Wei et al., 2012;Li et al., 2016;Grapenthin et al., 2018;Krabbenhoeft et al., 2018;Lay et al., 2018;Ruppert et al., 2018;Liu et al., 2019;Mann and Abers 2020;Ruppert et al., 2020;West et al., 2020). A large amount of work on the Alaskan earthquakes is based on different types of data and/or different methodological approaches. Therefore, some inconsistency and contradictions in the results obtained by different authors are inevitable, so that rupture process models are often different. That is why simple methods should be used to produce robust constrains to improve resolution of authoritative comprehensive determinations.

DATA
Seismicity of the region around Anchorage from January 01, 2006 through May 2019 is considered within 54°-64°N and 140°-160°W. The online search of the U.S. Geological Survey Advanced National Seismic System (ANSS) database provides a reasonably complete record of magnitude 2.5 or above earthquakes in the study area ( Figure 2); in particular, the catalog is complete in the circles around the two M7.1 epicenters, while has an apparent deficiency of earthquakes below magnitude 3.0 in the larger circle around the M7.9 epicenter in front of the continental Alaska. The graphs of the monthly number of the M ≥ 2.5 earthquakes confirm the stability of hypocenters' determinations in the ANSS catalog in advance of the major events.
For each of the total 20,803 earthquakes considered, the catalog reports the ANSS preferred magnitude M. The Gutenberg-Richter plot (Gutenberg and Richter, 1944) of the cumulative number of earthquakes in magnitude range from 2.5 to 7.9 ( Figure 2) follows the exponential best-fit trend line with the b-value of 0.868 (R 2 0.993). The plot is below its trend line in the magnitude range from below 5 up to 6.8, and then rises above due to the three major earthquakes on January 24, 2016, January 23, 2018, and November 30, 2018, whose parameters are listed in FIGURE 1 | Epicenters of the M ≥ 2.5 earthquakes (ANSS, 2006(ANSS, -2019 small blue circles) and the three major earthquakes (black stars) in Southern Alaska. Note: The two major earthquakes of M W 7.1 and the one of M W 7.9 under study are encircled with R 1°and 2.5°circles, respectively. Red line marks the boundary of the North American and the Pacific plates. The epicenters of the 1964 Great Alaskan earthquake (big red star) and its first aftershocks (red circles) are given on top the subsurface rupture zone (shaded pink). The epicenters of the 22 July 2020, M W 7.8 (red outline star) and its first day M ≥ 2.5 aftershocks (small yellow circles) at 105 km SSE of Perryville, Alaska and the 19 October 2020, M W 7.6 earthquake at 97 km SSE of Sand Point, Alaska (blue outline star) and its first day aftershocks (small green circles) occurred after the manuscript submittal and acceptance, correspondingly. The directions and rates (in cm/year) of Pacific plate convergence are indicated by the arrows. Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 584659 Table 1. Figure 3 displays the 4-D distribution of earthquakes in plots of magnitude, latitude, and longitude vs. time. The determination of depth is naturally biased by preset constants; therefore, to avoid potentially misleading conclusions, the depth vs. time distribution is not shown. Earthquakes with data points on pale yellow are not considered in the study. Figure 4 provides the map of the coefficients of unified scaling law for earthquakes, USLE, which generalizes the Gutenberg-Richter relationship as follows (Nekrasova and Kossobokov, 2005;Kossobokov, 2020): where N(M, L) is the number of earthquakes of a certain magnitude M expected in a year within an earthquake-prone area of diameter L; A and B are constants characterizing the annual rate of magnitude 5 events and the magnitude exponent analogous to a-and b-values of the Gutenberg-Richter relationship; and C estimates the fractal dimension of the epicenter loci at the given site. As evident from Figure 4C, the values of fractal dimension about 1.4 or larger highlight the highly fractured zone where subduction interacts with slippage along the Queen Charlotte transform fault. One can see that the highest values of C spread around the city of Anchorage accompanied with relatively low values of the logarithm of the seismic rate (A about −0.6). In contrast the area of Alaska-Aleutian trench is characterized by the highest level of earthquake rate (A about 0) and the lowest estimates of seismic locus dimension (C down to about 1). A, B, and C values at the epicenters of the three major earthquakes are listed in Table 1.

METHODOLOGY Characterizing Earthquake Sequences
As clearly noted by Bak et al. (2002), in essence, the USLE states that the distribution of inter-event times between earthquakes depends only on the value of the variable control parameter where τ is the time between the two successive earthquakes, M is the magnitude of the second one, and L is the distance between the two. Recent studies of seismic sequences associated with strong M ≥ 6 earthquakes in Central Italy (Kossobokov and Nekrasova, 2017) and major M ≥ 7 earthquakes in New Zealand confirmed the existence of steady levels of the USLE control parameter η, called "periods of stability," which, apparently, may change at the shorter times of critical transitions of a seismic regime, including those associated with regional or local catastrophic events, that is, major or strong earthquakes. To characterize seismic dynamics around each of the three major earthquakes, we use the same choices of the space-time limits as in Kossobokov and Nekrasova (2019). In particular, we apply the 50-point moving averages to control parameter η evaluated at an angular distance of 1°(or 2.5°) from the epicenter of the M ≥ 7 (or M ≥ 7.5) main shock for 128 days before and 128 days after its origin time ( Figure 5). (The choice of time interval allows for a sevenfold doubling of the 1-day period which is appropriate in analyzing either acceleration or deceleration of a daily time series.)

Characterizing an Earthquake Source
We analyze surface wave amplitude spectra at periods much longer than the earthquake duration to determine its source parameters derived from the moments of the stress glut rate. A detailed description of this methodology is given in Bukchin (1995). These techniques have been successfully applied in a number of studies of earthquake rupture processes (e.g., Aoudia et al., 2000;Lasserre et al., 2001;Clévédé et al., 2004). Considering the source of an earthquake in approximation of an instant point, we determine the source depth and focal mechanism (in terms of strike, dip, and slip angles) by systematic exploration of 4D parametric space, as well as seismic moment by minimizing the residual between the observed and calculated long period surface wave amplitude spectra. To improve the resolution of surface wave inversion, we performed joint inversion of surface waves and first arrival polarities. We calculate amplitude spectra using the model of weak lateral inhomogeneity for the Earth structure (Babich et al., 1976;Woodhouse, 1974). As is well known, the focal mechanism  FIGURE 4 | The Unified Scaling Law for Earthquakes (USLE) coefficients in Southern Alaska (Nekrasova and Kossobokov, 2002;Nekrasova and Kossobokov, 2019). Note: stars mark the epicenters of the three major earthquakes.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 584659 4 cannot be uniquely determined from the surface waves' amplitude spectra. We use P-wave first arrival polarities to select the optimal solution, which in detail description and constrains are given in Lasserre et al. (2001). Before inversion, we apply a smoothing procedure to the observed polarities, then use the frequency-time analysis program (FTAN) (Lander 1989a;Lander 1989b;Levshin et al., 1989) to isolate the Love and Rayleigh fundamental modes recorded by selected stations of the IRIS, GEOSCOPE, and GEOFON networks (Institut de physique du globe de paris, & ecole et observatoire des sciences de La terre de strasbourg, 1982;Scripps Institution of Oceanography, 1986;GEOFON, 1993). We simulated separately Love and Rayleigh fundamental modes, and compared them with observed records located by FTAN. The observed surface wave records were corrected for attenuation using PREM model and instrumental response. Analyzing the surface wave spectra at shorter periods, we describe the source as an approximation of the stress glut second moments (Backus, 1977a;Backus, 1977b). Using these moments, we determine integral estimates of the length of the source major and minor axes, l max and l min , its duration Δt, modulus of the axis average value of the instant centroid velocity |v|, and two angles: the angle between the strike axis and major axis, φ l , and the angle between the strike axis and instant centroid velocity axis, φ v .
Source dimension and spatial extent depend on the distribution of the stress glut rate in time and along the fault. In assumption of Gaussian distributions and 99% confidence levels, duration and length are set 2.5 and 3 times larger than their integral estimations, respectively.
To measure the effect of directivity, we consider a 1-D bilateral fault model with uniform slip from the nucleation point in two directions and define the modulus of uniform rupture velocity |v| and the ratio α β/λ, where β is the shorter length to the source edge from the nucleation point and λ is the total source length. These two parameters can be expressed through the second moments of the stress glut rate. The vector v is directed along the vector of the average instant centroid velocity.
Identification of the fault plane based on the source description in terms of the second moment of stress glut rate approximation is impossible, if the major axis of the source is much larger than the minor one and rupture propagation is directed along the line of nodal planes' intersection. In the case of a frequently encountered type of earthquake in which the major axis and the average instantaneous centroid velocity are directed along the strike axis, the fault plane can be reliably identified for a strike-slip event (λ 0°or λ 180°), and cannot be reliably identified for a dip-slip event (λ 90°or λ 270°). If the earthquake considered is not of the type mentioned above, then its fault plane identification must be resolved with other available evidence (Bukchin, 2017). Figure 5 shows epicenters of the M ≥ 2.5 earthquakes at angular distance of 1°and 2.5°from the epicenter of each of the three major earthquakes around Anchorage sampled in the time interval from 10 years before to 3 years after the origin time of the main shock or to the end of the available catalog. We specify those events occurring 128 days before (yellow circles) and after the main shock (red circles) and report a few related characteristics in the last rows of Table 1. The numbers of fore-and aftershocks differ dramatically; the ratio between them is 1:8, 1:262, and 1:21 for M ≥ 4.0 events around the January 2016, January 2018, and November 2018 main shocks, respectively. The delay time from the last foreshock differs from less than 7 h to about 5 days for the M ≥ 2.5 foreshocks and rises to a few months when M ≥ 4.0 foreshocks are considered. It is notable that the corresponding distances are much shorter for M ≥ 4.0 foreshocks (in particular, for both 2018 major earthquakes).

RESULTS
The sizes and shapes of aftershocks' time series are also very different in Figure 5: for the 2016, M7.1 event, we observe a nearly straight line extending for about 50 km; for the January 20518, M7.9 event, the pattern of epicenters looks like a butterfly,  Figure 6. The single period of stability of <η> for the entire interval of time (within a decimal order of their values 10 4 ≤ <η> < 10 5 ) was interrupted just for a few months after the January 24, 2016, M7.1 earthquake. In case of earthquake series around the January 23, 2018, M7.9 epicenter, the period of stability with 2 × 10 5 ≤ <η> < 2 × 10 6 was violated by a burst of activity from the end of March to the middle of May 2009 when seven earthquakes in magnitude range M5.0-5.9 occurred and then after the major shock. By the end of May 2019, the values of <η> increased to about 10 4 and may remain at this level as characterized by more intensive seismic activity in the area (confirmed by 4 × 10 3 ≤ <η> < 4 × 10 4 in June 2019−July 2020). Among the three series, the period in advance of the November 30, 2018, M7.1 Anchorage earthquake, characterized with 5 × 10 4 ≤ <η> < 2 × 10 5 , is the most stable in regard to the variance of the USLE control parameter. Moreover, by the end of May 2019, the values of <η> have increased above 2 × 10 4 , which indicates reaching the same decimal order of values as those in advance of the main major shock (confirmed by 5 × 10 4 ≤ <η> < 10 5 in June 2019-July 2020).  Let us zoom in to the 128-day intervals in advance and after the three major earthquakes. In the expanded view of <η> shown in Figure 7, we notice an apparent increase of <η> by a factor of 2 in advance of the January 24, 2016, M7.1 event, as well as a drop of <η> by a factor of 2 in advance of the other two major earthquakes, which does not appear to offer a method for short-term earthquake forecasting. The decay of the aftershock series is also rather diverse: the best fit of the 50-point moving average <η> in the three aftershock series of the January 2016, January 2018, and November 2018, as a function of time after the main shock, t, is the power law equal to 80.8 × t 1.18 (R 2 0.988), 22.6 × t 1.18 (R 2 0.949), and 19.9 × t 0.88 (R 2 0.932), respectively.
To determine the source parameters of the three major earthquakes at 600 km around Anchorage as derived from the moments of stress glut rate (Bukchin, 1995), we selected the following spectral bands (spectral bands in brackets correspond to the second moment approximation): 70-250 s (40-60 s) for the two M7.1 major earthquakes in 2016 and 2018 and 150-250 s (40-60 s) for the M7.9 earthquake in 2018. Table 2 shows for each of the three major earthquakes, the two nodal planes (first row) and, for each nodal plane, the six integral source parameters derived from the second stress glut moments (rows 2-7) along with the values of total residuals in the optimal solution (row 8). The six parameters characterize the earthquake model size, orientation, duration, and rupture propagation. The characteristics were obtained from making use of 36, 42, and 40 fundamental Love and Rayleigh mode records of the January 24, 2016, January 23, 2018, and November 30, 2018 earthquakes,  respectively. Selecting the records, we did our best to satisfy the following conditions: i) signal quality, characterized by the value of the polarization anomaly, should be less than 15°, and ii) the azimuthal distribution of stations should be rather homogeneous. Figure 8 summarizes modeling of the January 24, 2016, M W 7.1 earthquake source by displaying the location of stations used in the determination a), the residual function which minimum defines the depth of the instant point source b), the double couple beach ball solution c), the focal mechanism modeled by three angles d), and the 1-D bilateral model of faulting e). The solution is not a pure dip slip, which complicates identification the fault plane. The residuals for the nodal planes of this earthquake ( Table 2, last row) are about the same and do not allow for a confident choice of the fault plane; however, the abovementioned shape of the distribution of its Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 584659 8 aftershocks in Figure 5 favors the mechanism presented in Figure 8D.
By analogy, Figure 9 summarizes modeling of the January 23, 2018, M W 7.9 earthquake source. In this case, the focal mechanism is practically pure strike-slip motion. The residuals, that differ by about 17%, favor the slippage along the near meridian strike of the narrow "butterfly body" of the aftershock distribution in Figure 5. As evident from Figures 3, 5, 10, the aftershocks of this earthquake split into three clusters with the narrow middle one associated with the largest number of aftershocks ( Figure 10). Figure 11 is analogous to Figures 8, 9. The November 30, 2019, M W 7.1 earthquake is modeled as a pure dip slip. In this case, there is no evidence for a confident choice of the fault plane, so that the optimal choice based on nearly the same residuals (which differ by a negligible 3% of their value, Table 2) remains rather uncertain; an option of both nodal planes involved in the normal faulting at a depth of about 70 km ( Figure 11B) cannot be dismissed (the apparently bimodal shape of the graph may explain different determinations reported by other agencies).
It is worth mentioning that there are no significant differences in focal mechanisms from the parameters provided by USGS, Global CMT, and ESMC. However, the depths of the three major earthquakes provided by these agencies differ from our determination (e.g., compare 100, 20, and 70 km with 129, 14.7, and 46.7 km by USGS, respectively). A significant difference is the choice of a nodal plane representing the fault plane for the January 23, 2018 event; the USGS determination favor E-W plane, while our analysis of the complex aftershock data (Figures 3, 5, 10) strongly suggests the orthogonal one as the correct representative of the fault plane.

DISCUSSION AND CONCLUSION
The determination of the spatial and temporal distribution of slip over the fault area of large earthquakes is essential for better comprehension of the mechanics of the seismic rupture and its environmental impact (e.g., Yagi and Kikuchi, 2000;Bouchon et al., 2001;Clevede et al., 2012). We use the integral source parameters given by the stress glut rate moments of total degree 2 for constraining the models of the three major earthquakes at 600 km around Anchorage, Southern Alaska. Beyond point source moment tensor and hypocenter determinations, we provide average estimates of the rupture length and width, duration, and velocity, making use of the low degree moments of stress glut rate (Backus, 1977a;Backus, 1977b). To our knowledge, a few earthquakes have been studied from this point of view (Gusev and Pavlov, 1988;Bukchin, 1995;Gomez et al., 1997a;Gomez et al., 1997b;Aoudia et al., 2000;McGuire et al., 2000;McGuire et al., 2001;Clévédé et al., 2004).
Our analysis of the source models for the three earthquakes yields similar results to other moment tensor solutions. However, for the 2018 Kodiak earthquake, we infer that the conjugate plane in the final USGS NEIC solutions is the fault plane. This is similar to the findings presented in the comprehensive study by Lay et al. (2018) based on seismic, GPS, and tsunami data. They conclude a multiple fault rupture dominated by right-lateral slip on an SSE trending, Specifically, their "primary faulting with ϕ 155°and δ 72°, with slip of up to 15.6 m in the crust and uppermost mantle," is confirmed and complimented by our robust estimates based on the low-degree moments of the stress glut rate (Table 2; Figure 9).
In general, major earthquakes of Mw ≥ 7 are complex. The geometry, timing, and slip distribution of the complex nonlinear system of the lithospheric blocks and faults (Keilis-Borok, 1990) involved are usually not well resolved. Their rupture may propagate both upward into the crust and downward into the mantle (Liu et al., 2019;Ruppert et al., 2020). As a result, the expected ground shaking could be under-or over-predicted by an order of magnitude (Grapenthin et al., 2018;Mann and Abers 2020), causing unexpected societal impacts and errors in seismic hazard assessment (Wyss et al., 2012;West et al., 2020).
Seismic energy accumulation, release, and redistribution in the lithosphere are not yet well studied. It appears that seismic diversity is a cornerstone of understanding complex dynamics of a system of blocks and faults in advance and after catastrophic events, the main contributors to energy release (Ben-Zion, 2008;Zaliapin et al., 2008). Characterization of the fore-and aftershock sequences of the recent major earthquakes in Southern Alaska confirms the existence of the long-term periods of seismic stability defined by the averages <η> of the USLE control parameter that are interrupted by mid-or even short-term bursts of activity associated with catastrophic events. However, at this time, neither of the two M7.1 events considered in this study showed a change in the level of <η> observed in advance of their origin times; the January 23, 2018, M7.9 earthquake may become an exception, if the level of <η> reached by the end of May 2019 persists. It seems premature to discuss if the variability of <η> can be useful in operational earthquake forecasting (Kossobokov et al., 2015) of seismic catastrophes, due to the yet rather small number of case studies: five cases in Central Italy (Kossobokov and Nekrasova, 2017), nine in New Zealand , and three in Southern Alaska. Nevertheless, the apparent major seismic activity, in particular, the one associated with the M7.9 strike-slip earthquake 280 km SE of Kodiak and the most recent M7.8 Alaska Peninsula thrust faulting, right at the borders of the rupture zone of the 1964 Great Alaska, M9.3 mega-earthquake (Figure 1), deserves special attention, with continued further studies and monitoring of this region.

AUTHOR CONTRIBUTIONS
BB and AF contributed the solutions and interpretation of the three major earthquake sources; VK and AN contributed presentation of seismic flow dynamics in advance and after the three major earthquakes.

FUNDING
This study was carried out as part of the RF state task of Scientific Research Works on "Development of methods of seismic data analysis in order to study the source, environment, seismic hazard" (No. 0143-2019-0007). No other funding received.