Cascading rupture of a megathrust.

Understanding variability in the size and location of large earthquakes along subduction margins is crucial for evaluating seismic and tsunami hazards. We present a coseismic slip model for the 2021 M8.2 Chignik earthquake and investigate the relationship of this earthquake to previous major events in the Alaska Peninsula region and to interseismic coupling. Stress changes from the 2020 M7.8 Simeonof event triggered the Chignik event, and together, the earthquakes partially filled an unruptured section along a 3000-km subduction margin that has experienced a series of ruptures along almost its entire length over the past century. Variations in coupling and structural characteristics may make the region more prone to nucleating M7 to M8 events rather than larger M > 8.5 earthquakes. Stress changes and rupture areas suggest that the two recent earthquakes may be part of an 80-year-long rupture cascade and may have advanced seismic hazard in the region.

INTRODUCTION

During 2020–2021, an earthquake sequence at least partially ruptured two neighboring segments of the Alaska-Aleutian megathrust (where the Pacific plate meets and descends beneath the North American plate) offshore of the Alaska Peninsula. This complex convergent system displays highly variable coupling (the amount of slip that is accumulating toward an earthquake between the downgoing and overriding plates; at 100% coupling, slip is accumulating at the full rate of relative plate motion) () as well as distinct changes in structural properties along the plate interface () and in overriding plate motion. The spatial and temporal proximity of these earthquakes and their locations in relation to previous events offer deeper insights into stress transfer and earthquake triggering, conditions required for rupture propagation, and why some areas may be more prone to generating predominantly magnitude (M) 7 to M8 rather than larger M8.5 to M9+ earthquakes.

The event sequence started on 22 July 2020, when the M7.8 Simeonof megathrust earthquake at least partly filled (–) the Shumagin seismic gap (Fig. 1). This enigmatic segment of the plate boundary last partially ruptured in an M7.4 earthquake in 1917 (). An M7.6 strike-slip event within the downgoing plate struck in October 2020 to the west of the Simeonof event (). The Simeonof event was followed by the 29 July 2021, M8.2 Chignik megathrust earthquake, the largest event in the United States since the 1965 M8.7 Rat Islands earthquake. The Chignik earthquake ruptured a patch of the megathrust directly east of the Simeonof rupture zone. Both earthquakes, separated in time by a year, mostly ruptured the 20- to 40-km depth range within segments of the interface that appear to be partially coupled () and thus did not trigger substantial tsunamis. They occurred in an area that previously generated several M7.5+ earthquakes [1938 M8.3 Alaska Peninsula, 1948 M7.5, and 1788 M8(?)] and is bounded by the rupture zones of much larger earthquakes: the 1957 M8.6 Central Aleutian and 1946 M8.6 Unimak earthquakes to the west and the 1964 M9.2 Prince William Sound earthquake to the east (Fig. 1). The latter three events produced Pacific-wide tsunamis that resulted in substantial damage and casualties in Hawaii and along the west coast of the contiguous United States (, ). With the exception of an ambiguous event in 1788, all of these previous large earthquakes occurred within 30 years of each other, suggesting a megathrust cascade ().

Fig. 1.

Overview of the 2020–2021 Alaska sequence, plate interface coupling, and historic earthquakes.

(A) Location of the main map relative to rest of Alaska. Red outline shows area of main map. (B) 2020 Simeonof earthquake slip model with 1- and 2-m contours () shown in black, 2021 Chignik earthquake slip model with 1- to 5-m contours shown in purple, and outlines for 1- and 2-m contours of the two updated 1938 slip model estimates () shown in light and dark blue. Red outlines show rupture areas for major historic earthquakes (, , ). Segment divisions are derived from historic rupture areas and structural changes along the margin (, , ). The Kodiak segment extends east from the Semidi segment and includes the region that experienced high slip (the Kodiak asperity) during the 1964 M9.2 earthquake (). Yellow stars show epicenters of M ~7+ earthquakes over the past century. Orange star is the epicenter of the M7.4 1917 event. Interseismic coupling estimate is shown in grayscale shading (, ). Earthquake focal mechanisms for the 2020 and 2021 earthquakes are from the USGS. Thick dashed line shows approximate location of a landward-dipping normal fault system imaged in a seismic reflection line just to the west of the Shumagin Islands (). Lateral extent of fault is uncertain. Dotted line marks location of trench.

The Simeonof-Chignik earthquakes are notable for their short temporal and spatial separation and their location along a highly heterogeneous section of the plate interface that evolves from fully creeping to the west to fully coupled in the Kodiak segment to the east (Fig. 1) (, ). This sequence thus presents an opportunity to investigate part of a possible megathrust cascade against the backdrop of stress transfer–driven triggering and coupling-related rupture propagation and arrest. Here, we first resolve the slip distribution for the 2021 M8.2 Chignik earthquake before modeling the stress transfer onto the megathrust due to the 2020–2021 event sequence, which we interpret in the context of previous large tsunamigenic earthquakes, earthquake rupture models, and future scenarios.

RESULTS

M8.2 Chignik finite fault model

Using teleseismic and high-rate Global Navigation Satellite System (GNSS) waveforms, static GNSS offsets, and interferometric synthetic aperture radar (InSAR) line-of-sight (LOS) offsets for separate patches of coherent pixels (Fig. 2) that are either tied to static GNSS offsets or “floating” with an ambiguity parameter estimated during the inversion (), we infer the coseismic slip model in Fig. 3 (details in Materials and Methods). We adopted the epicenter from the U.S. Geological Survey (USGS) estimate and pinned its depth (26.2 km) to the slab depth as estimated by Slab2.0 (). The Chignik event originated at the northeast edge of the Simeonof rupture zone () and propagated slightly to the southwest but not substantially past the 1-m Simeonof rupture contour (Fig. 1). The bulk of the energy release occurred near the hypocenter during the first 20 s of the Chignik event, with a maximum slip amplitude of about 6 m between a depth of 20 and 30 km (Fig. 3). From 20 to 60 s after the nucleation time, slip propagated northeastward toward Chirikof Island, releasing between 1 and 1.5 m of slip in the eastern part of the Semidi segment, with a slight migration of higher slip toward shallower depths (Fig. 3B). No substantial slip (>1.5 m) occurred above a depth of 15 km along the entire rupture zone. Our estimated rupture velocity suggests rather slow but generally homogeneous slip release up-dip of the hypocenter (1 to 1.5 km/s) and a much faster rupture down-dip (up to 4 km/s; Fig. 3C). Our slip model yields a moment rate function that is different from the USGS estimate in that the bulk of energy is released slightly earlier, peaking at 23 s, and resembles other teleseismic estimates including Institut de Physique du Globe de Paris (IPGP)/SCARDEC (). The slightly lower peak energy release of our model seems to be compensated by a longer tail of moment release past 110 s that ends at 130 s. Differences between the moment rate functions may be a result of the frequency bands used in the inversions. As the USGS model uses long-period surface waves, a low-frequency band (0.01 to 0.125 Hz) was applied in that model. In our model, which uses static geodetic data that provide strong constraints on the spatial slip distribution, a higher-frequency band (0.01 to 1 Hz) was used to allow recovery of more details of the slip evolution.

Fig. 2.

Overview of near-field geodetic data.

GNSS stations are marked and show static horizontal (red) and vertical (blue) coseismic displacements. Traces of high-rate GNSS time series are plotted next to AB13, AC21, AC13, and AC45 and include east (black), north (blue), and vertical (magenta) components (see AB13 plot for component labels and time and displacement scales). Wrapped Sentinel-1 InSAR data used in the slip inversion are included with insets for Chirikof Island (path 102, frame 407, descending) and the Semidi Islands (path 7, frame 180, ascending); other areas were heavily affected by atmospheric noise. OT, origin time.

Fig. 3.

Slip model for the 2021 M8.2 Chignik earthquake.

(A) Slip amplitude (colors) and rake (arrows) along strike and dip (depth) for each rectangular subfault. Rake shows motion of upper plate relative to the downgoing slab. The red star marks the hypocenter location. Contours indicate location of the rupture front in seconds after earthquake nucleation. (B) Slip model from (A) projected into map view with Slab2.0 () depth contours in dashed white lines. The red star is the Chignik epicenter, and the white star is the relocated Simeonof epicenter (). (C) Rupture velocity for each subfault with a total slip larger than 80 cm. (D) Moment rate function (MRF) estimated from our preferred model (black line and gray shaded) and the USGS (pink line).

Details of the model fit to the data are shown in figs. S1 to S3. Slight misfits in the near-field static GNSS data (fig. S1) fall within the data uncertainties. The most substantial differences exist for farther-field vertical static GNSS offsets (fig. S1D), where several stations have misfits of ~1 cm. Adjusting the fault dip angle did not substantially improve the overall model fit, although our model fault had a constant dip. On the basis of seismicity, the slab begins to experience a steeper dip beneath the Alaska Peninsula. A more complicated model that included varying dip with depth may be required to improve fit in this region.

The high-rate GNSS waveforms are very well reproduced in phase and amplitude (fig. S2). We see notable misfit in amplitude at the near-field site AC13, which experienced some of the largest dynamic motions. This mismatch is likely due to model regularization or weighting decisions made. As the static offset at this site is well matched, we do not consider this problematic. Higher-frequency oscillations in the data are noise, which is generally larger in the vertical than the horizontal components due to satellite-station geometry. Fits to teleseismic waveforms are equally good in phase and amplitude for different seismic wave phases (fig. S3).

The fit to the InSAR data in both predicted phase and LOS deformation is very good, suggesting that our model captures the geometry of and the total slip on the fault well (figs. S4 to S8). Slight misfits to the phase observations are likely due to residual atmosphere errors (e.g., Sutwik Island; fig. S7C).

Checkerboard tests (fig. S11) show that major features are resolved given the available data distribution. Resolution of slip patches near the down-dip end of the model fault plane is good because of the nearby presence of data along the Alaska Peninsula. At shallower depths, only a very smooth version of the input slip distribution is recovered because of the lack of data coverage except along the southwestern edge of the model, where GNSS and InSAR data on the Shumagin Islands improve resolution.

Stress change

We investigate coseismic static Coulomb stress change (see Materials and Methods) to assess the impact of the stress changes generated by the 2020 M7.8 Simeonof earthquake (), particularly up-dip of the rupture and surrounding the 2021 M8.2 Chignik hypocenter, the 2020 M7.6 Sand Point strike-slip event on the area of the Chignik hypocenter, and the M8.2 Chignik earthquake up-dip of the rupture and along strike to the east, particularly in the eastern Semidi and Kodiak segments.

Static coseismic Coulomb stress change due to slip during the 2020 M7.8 Simeonof earthquake projected onto the megathrust is shown in Fig. 4A. An area of substantial positive stress increase of >2 bars occurs just south of the Shumagin Islands. In the shallow region directly inboard of the trench, stress increases reach 0.5 to 1 bar. Limited areas of very high stress loading (>5 bars) exist to the southwest of the Simeonof hypocenter and along the down-dip edge of the model fault plane. In terms of shear and normal stress changes (fig. S12), the area directly up-dip of the Simeonof hypocenter experienced a decrease in shear stress and positive normal stress change (fault unclamping). Up-dip and to the west of the hypocenter, shear stress increased, while normal stress underwent a negative change, indicating fault clamping. The stress change around the 2021 M8.2 Chignik hypocenter is slightly more complex with regions of stress loading and stress drop. A vertical cross section (Fig. 4B) from the surface to a 50-km depth shows the Chignik hypocenter clearly embedded in an area of stress loading between 0.5 and 1.5 bars ranging in depth from about 15 to 35 km on the megathrust surface. Here, shear stresses increased on the megathrust, while normal stresses showed positive change, suggesting shear loading and fault unclamping (fig. S12B). This result is consistent with the peak slip region of the Chignik event being advanced toward failure.

Fig. 4.

Coulomb stress change due to the 2020 Simeonof earthquake.

Receiver fault is the model fault plane for the 2021 Chignik earthquake. (A) Stress change projected onto the megathrust as estimated by Slab2.0 (). Black dot with red border marks the hypocenter of Simeonof event. Red dot marks the hypocenter of Chignik event, and arrow points to the focal mechanism for that earthquake. Grid shows the model fault plane for Simeonof event (). Two coefficients of friction (0.1 and 0.4) are shown. Limits of scale were chosen to avoid saturation. Gray dashed lines show depth contours of slab. (B) Vertical profiles of stress change. Profile lines are shown in (A). White star marks the hypocenter of 2021 Chignik earthquake. Arrows indicate the sense of fault motion.

While it is beyond the scope of this study to develop a postseismic model for the Simeonof earthquake, observed postseismic deformation (fig. S13) displays similar patterns in terms of spatial extent and direction as the coseismic deformation. This suggests that afterslip may be occurring in areas immediately surrounding (and perhaps within) the coseismic rupture zone and that the overall stress changes due to afterslip will be similar to the coseismic, albeit of significantly smaller magnitude.

To assess the impact of the stress changes due to the 2020 M7.6 Sand Point strike-slip earthquake, we estimate a preliminary slip model (Materials and Methods and fig. S14) using the same approach as for the Chignik event but using only static GNSS and teleseismic waveform data as constraints. Our model provides a very good fit to the geodetic (fig. S15) and seismic (fig. S16) data without requiring any significant slip on the megathrust even when it is included in the model. Tsunami waveforms may require additional sources than a strike-slip fault, with megathrust slip or underwater landslides being potential contributors (). Coulomb stress changes from the Sand Point (fig. S17) event resolved onto our model fault plane for the Chignik earthquake suggest minimal impact (≤0.01 bar of stress loading) from the former on the Chignik hypocenter. Given this result, impact from stress changes due to postseismic slip from the Sand Point earthquake is expected to be negligible.

We use our slip model for the 2021 M8.2 Chignik event to assess the coseismic static stress changes that it induced on the megathrust (Fig. 5). In addition to broad stress decreases in the region of the rupture itself, our results clearly show substantial stress increases exceeding 2 bars up-dip of the entire rupture (Fig. 5). Past the eastern end of the modeled slip region, near Chirikof Island, the stress change becomes more complex as an area of increased stress is interrupted by a narrow region of stress drop occurring at a depth of 10 to 40 km. Further east, toward Kodiak Island, stress increases of 0.5 to 1.5 bars occur along the entire down-dip width of the megathrust. This increased stress band, while narrow, is a notable stress increase ~200 km from the hypocenter along the western edge of the Kodiak asperity of the 1964 rupture.

Fig. 5.

Coulomb stress change due to the 2021 Chignik earthquake.

Receiver fault is the Kodiak segment of the megathrust, with contours (gray dashed lines) defined by Slab2.0 (). Stress change is projected onto megathrust. Grid shows the model fault for Chignik event, and white star marks the hypocenter. (A) Stress change for coefficient of friction of 0.1. (B) Stress change for coefficient of friction of 0.4.

DISCUSSION

The Simeonof-Chignik sequence occurred within a section of the Alaska-Aleutian megathrust that exhibits highly heterogeneous interseismic coupling. Geodetic estimates of the coupling show a transition from very strong coupling values within the Kodiak segment to moderate-to-low values through the Semidi and Shumagin segments to creeping in the Sanak segment (Fig. 1) (). Lateral variations in the coupling estimates are robustly determined, while along-dip variations have greater uncertainties due to the lack of offshore data constraints and are dependent on the chosen model regularization (). It is important to note that geodetic estimates of coupling are based on the data collected within the past few decades and may not fully represent coupling distributions over longer geologic periods. In some areas, lateral changes in coupling align with megathrust segment boundaries defined by previous ruptures and structural differences (Fig. 1) such as along the Sanak-Shumagin and Shumagin-Semidi boundaries. This alignment may represent more long-lived coupling distributions. In other places, the present-day coupling does not neatly align with the segment boundaries. This is the case for the area of the eastern Semidi and western Kodiak segments, where paleoseismological and paleotsunami data suggest that a nonpersistent boundary may exist ().

Present-day lateral coupling boundaries may influence earthquake rupture nucleation and arrest as well as slip release, especially in regions with long-lived coupling boundaries. The Simeonof earthquake nucleated just to the west of the Shumagin-Semidi boundary and ruptured westward, while the Chignik earthquake nucleated just to the east of that boundary and ruptured eastward (Fig. 1). Both the Simeonof and Chignik events ruptured predominantly deeper (15 to 40 km, with most slip released between 20 and 30 km), moderately coupled portions of the Shumagin and Semidi segments. Almost all of the substantial slip (>1 m) during the Simeonof earthquake stayed within the coupling-defined bounds of the eastern Shumagin segment. All of the substantial slip from the Chignik earthquake was contained within the Semidi segment, with high-slip (>3 m) areas limited to the moderately coupled western half of the Semidi segment. Although the exact rupture area is uncertain, the 1917 event appears to have nucleated near the Shumagin-Semidi boundary and ruptured westward, while the 1938 event nucleated to the east of the boundary and propagated eastward (). Earlier large events in 1788 and 1847 may also have been bilateral or paired earthquakes centered around the Shumagin-Semidi boundary (, ).

Expanded sets of detailed global observations have revealed that subduction zones can exhibit a wide variety of behavior, with rupture modes and extent varying spatially and temporally within a single system () instead of following the simple model of major earthquakes repeating in the same area with a recurrence time based on plate motion rates (). Many aspects of this variability, including the occurrence of complementary ruptures or megathrust rupture cascades (temporally clustered ruptures that fill in a megathrust like puzzle pieces without large areas of overlap), are influenced by changes in stress loading of the megathrust fault system. Stress triggering and cascades have been observed in a number of other subduction systems. A prime example is Sumatra, along which the 2004 M9.1 event began a southward cascade of earthquakes that included the M8.6 Nias event in 2005 and M8.5 and M7.9 events in 2007 (). Similarly, over a 20-year period during the 20th century, the Kurile-Hokkaido margin ruptured almost its entire span in a series of M8+ earthquakes ().

The Alaska-Aleutian margin exhibited a spectacular cascading rupture sequence when five M8+ events broke almost the entire margin from the far western Aleutians to South Central Alaska in a period of less than 30 years, the last of which—the 1965 M8.7 Rat Islands earthquake—occurred 55 years ago. A feature of cascading ruptures is that each rupture usually heightens the hazard in neighboring asperities (), although several factors can complicate this scenario (see below). Therefore, a major question is whether the Simeonof-Chignik sequence, with two M7.5+ events occurring a year apart on adjacent segments of the megathrust with distinctly different properties, is part of an emerging cascade or whether these events are the last components of the previous one, as they filled in regions that did not rupture during the main 30-year cascade.

Our results show that the stress changes induced by the Simeonof earthquake rupture promoted the Chignik earthquake (Fig. 4). The relationship of previous earthquakes to the Simeonof earthquake is more complicated. The exact rupture area and hypocenter of the 1917 earthquake are uncertain, but it may have ruptured a similar area to the 2020 Simeonof event. It was likely a smaller event with slip of 0.5 to 1.5 m () and, assuming constant present-day coupling since 1917 (), 1.6 m of slip has accumulated since 1917. A recent study () investigated how the coseismic static stress changes from several M7+ earthquakes may have influenced the 2020 Simeonof earthquake and suggested that the static stress changes from the coseismic slip alone did not promote the 2020 earthquake. However, their model used a slip model for the 1938 earthquake that included substantial coseismic slip in the vicinity of the Simeonof hypocenter. A recent reevaluation of the 1938 event () suggests that most of the slip may have occurred closer to the trench. A comparison of the rupture of the 2021 event with the revised 1938 estimates (Fig. 1) shows that the high-slip regions (>1.5 m) of the events do not substantially overlap but appear to be complementary ruptures, suggesting that shallower slip during the earlier event may be a more accurate representation. Depending on the version of the updated slip model used, positive Coulomb stress changes were induced either in the vicinity of the hypocenter or in the area of major energy release of the Simeonof earthquake, or both. Viscoelastic postseismic relaxation from pre-2020 events (), which was not accounted for in the previous Simeonof Coulomb model (), could also load the Simeonof-Chignik region. Postseismic relaxation, especially due to large events such as the 1964 M9.2 megathrust event, has promoted failure on faults, even in the far field (, ). As mentioned previously, the postseismic deformation following the Simeonof earthquake (fig. S13) likely enhanced the positive stress change around the Chignik hypocenter generated by the mainshock.

The Alaska Peninsula presents a unique opportunity to investigate what may govern rupture size. The Shumagin segment has nucleated four major earthquakes over the past century, all with M < 8 (Fig. 1). The Semidi segment has generated several M7+ events and two M8.2 to M8.3 events directly adjacent to the boundary with the Shumagin segment. Why does this part of a prolific, great earthquake–generating margin rupture in multiple M7 to M8 events rather than in fewer, larger earthquakes, especially in the instance of the recent sequence? Various mechanisms have been proposed to explain the sequences of partial megathrust ruptures. Fracture mechanics suggests that a series of partial ruptures may be promoted in a scenario where creeping (no coupling) sections of a fault are adjacent to more highly coupled sections if thresholds related to stress loading, mechanical properties of the fault, and the migration of creep into the coupled zone are exceeded (). In this scenario, ruptures nucleate near the boundary between the creeping and coupled zones, and temporal clustering of events is related to the rate of afterslip. Ultimately, sequences of partial ruptures are driven by variations in the stress field and differences in local frictional conditions. Investigations of earthquakes along the Sumatra margin suggest that sequences may occur in areas where only part of the accumulated slip deficit is released during an earthquake and that this partial release may be related to nonpermanent rupture barriers caused by stress changes from past earthquakes (). Along the Chile margin, it has been proposed that frictional variations on the interface, especially those related to pore fluid pressure, may result in more frequent, moderate-sized earthquakes in the deeper portions of the megathrust, while the shallower megathrust ruptures in less frequent great earthquakes ().

Multiple mechanisms, in particular abrupt changes in coupling and frictional differences caused by variations in the downgoing plate fabric, may have contributed to partial megathrust ruptures along the Alaska Peninsula. In the Shumagin segment, the megathrust is weakly to moderately coupled (<40%), while the adjacent western half of the Semidi segment is more highly coupled (maximum <70%), setting up a scenario where partial ruptures may be promoted as described above (). The 2021 M8.2 Chignik earthquake nucleated just to the east of the Shumagin-Semidi boundary and propagated eastward, almost fully contained within the Semidi segment. The 1938 M8.3 earthquake also nucleated very close to this boundary, with most of the slip occurring to the east within the strongly coupled eastern Semidi segment (). Downgoing plate fabric, structure, and hydration may play a major role in encouraging partial megathrust ruptures. Basement faulting and thin sediment cover on the downgoing plate in the Shumagin segment leads to roughness along the megathrust that promotes the occurrence of series of moderate-sized earthquakes rather than a smaller number of great earthquakes (). The irregular surface of the plate interface in the Shumagin segment may also be directly related to the low-to-moderate coupling observed there as rough surfaces on the downgoing plate may lead to heterogeneous stresses that lead to creeping behavior (, ). The idea that the Shumagin-Semidi region is more prone to nucleating M7 to M8 earthquakes does not preclude larger earthquakes nucleating in adjacent regions and propagating through the area. This may have occurred during events in 1788 and 1847, although these larger earthquakes may have been closely spaced series of smaller events rather than single, great earthquakes (, , ).

Given that the 2020 and 2021 earthquakes occurred along the deeper portion of the interface only, our modeled Coulomb stress changes suggest that the shallower up-dip regions from both ruptures could now be closer to failure. The 2021 Chignik earthquake resulted in a positive Coulomb stress change at shallow depths with particularly large (up to 2 bars) increases along the southern edge of the rupture area (Fig. 5). This coincides with the northern edge of the high-slip (2+ m) region of the 1938 rupture as recently estimated (Fig. 1) (). Assuming that the megathrust gained present-day coupling within a few years of that event, up to 5 (based on easternmost coseismic estimate, which is in a highly coupled region; Fig. 1) or 2.5 m of slip (based on westernmost coseismic estimate, which is in lower-coupled region) had accumulated by 2020, meeting or exceeding the estimated slip during the 1938 earthquake. A major caveat to this potential increased seismic hazard is that the 1938 earthquake occurred at shallower depths where we do not have direct geodetic constraints and coupling estimates are strongly influenced by model regularization (). Another reason to question the seismic hazard around the region of the 1938 earthquake is the extremely rapid afterslip observed at AC13 on Chirikof Island (fig. S18). The magnitude of postseismic deformation matched the coseismic deformation within 2 months of the Chignik earthquake, suggesting that (at least temporarily) the shallow megathrust in this region may be creeping. A thick layer of overpressured sediment on the downgoing plate likely promotes creep along the shallow portion of the megathrust along the Semidi segment, while sediment dewatering and lithification would allow higher coupling and earthquake nucleation at greater depths (). This change in behavior with depth along the Semidi segment may explain why the Chignik earthquake ruptured a deeper portion of the megathrust and did not propagate to the trench. Explanations for the possible contradiction between this occurring in the same area as a historic great earthquake include a more compact 1938 rupture zone located east of Chirikof Island (in which case the slip during that earthquake would be higher) or much more heterogeneity along the megathrust than we can resolve with the current geodetic and offshore seismic datasets. As noted above, the areas of high slip for the Chignik and 1938 ruptures do not appear to substantially overlap, which suggests that they may be complementary ruptures that broke different parts of the megathrust rather than repeated ruptures of the same fault patch (Fig. 1). The 1938 earthquake and the Chignik earthquake may be part of a related but semi-independent rupture cycle in which the shallow interface ruptures during less frequent, larger earthquakes than the intermediate depth megathrust. These superimposed cycles have been observed at other subduction zones including Sumatra, Nankai, and Tohoku-Oki ().

Another region of major concern is the shallow interface outboard of the 2020 Simeonof earthquake, as no historic earthquakes (other than an enigmatic earthquake in 1788) are known to have ruptured that region, and it shares characteristics such as a heterogeneous plate interface, a small frontal prism, and strongly faulted crust with margins known to host tsunamigenic earthquakes (). The Simeonof earthquake induced positive stress change values exceeding 0.5 bars across a broad portion of the shallow interface with a limited area of up to 2 bars of stress increase directly outboard of the hypocenter (Fig. 4). Although this suggests that the shallow interface should be closer to failure, the region may lack the frictional properties needed to rupture in a major tsunamigenic earthquake. Geologic investigations have not been able to detect earthquake-induced uplift or high tsunami deposits in the Shumagin Islands (). For that to be the case, an earthquake nucleating in the area would either have to be deeper (like the 2020 Simeonof earthquake) or rupture only the shallowest interface in an M < 8 (). Seismic imaging shows that, within the Shumagin segment, the subducting oceanic crust is deeply faulted and has a fairly thin sediment cover, a setting that may promote rupture in multiple, smaller earthquakes (, ). This is in contrast to the oceanic crust of the Semidi segment, which is covered with a substantial layer of sediment and may be more prone to generating great earthquakes as a result (, ). Offshore seismic data have also revealed that a large, landward dipping normal fault cuts through the upper slope inboard of the trench and intersects the megathrust at a depth of ~35 km (). The normal fault marks the seaward limit of abundant seismicity and, interestingly, of the 2020 Simeonof rupture (Fig. 1), suggesting that it may influence slip behavior on the megathrust and limit rupture extent. Investigations of normal faults in subduction systems, including Tohoku-Oki, suggest that lower effective friction exists along the interface up-dip of the faults compared with the down-dip region (, ). While these observations suggest that the shallow region of the Shumagin segment may not be expected to generate great tsunamigenic earthquakes, it should be noted that great earthquakes may have propagated into the Shumagin segment from neighboring regions and may still pose substantial hazard ().

Our Coulomb stress change modeling suggests that the Chignik earthquake produced positive stress changes in the vicinity of the western edge of the Kodiak segment that last ruptured during the M9.2 1964 earthquake. That region experienced a substantial coseismic stress drop during 1964, but the effects of viscoelastic postseismic deformation, which is ongoing in the present day, and strain accumulation above the relocked asperity have likely partially reloaded the system. Assuming that the Kodiak segment gained present-day coupling values within a few years of the 1964 earthquake, ~3 m of slip have accumulated () or roughly 15 to 20% of the slip released in that area during the M9.2 event (). Geologic evidence suggests that the southwestern end of the 1964 rupture contains a nonpersistent boundary, leading to variable-sized ruptures (). It is possible that the Kodiak segment could rupture on its own, and several such ruptures have occurred in the recent geologic past in addition to the multisegment ruptures such as in 1964 (, ). The 2020–2021 sequence, and particularly the Chignik earthquake, may have brought parts of this segment of the megathrust closer to failure.

The question remains as to whether the 2020–2021 earthquakes are part of the 20th-century rupture cascade or whether they represent the beginning of a new series. Given the calculated and implied stress changes and the spatiotemporal relationships between the 2020 and 2021 ruptures and past earthquakes, we propose that the Simeonof and Chignik earthquakes are related to the cascade of megathrust earthquakes that began along the Alaska-Aleutian megathrust during the last century. This would suggest a cascade lasting more than 80 years, a span also observed in Sumatra and Japan (), and may be consistent with the variable coupling behavior. While surrounding regions of the megathrust were brought closer to failure by the sequence, the paucity of offshore data and earthquake records make the pre-2020 stress state unclear, especially for the shallowest portion of the megathrust. A more complete understanding of the region will require substantially more offshore data, particularly seafloor geodetic data, making this area a prime target for future additional instrumentation and focused studies.

MATERIALS AND METHODS

InSAR analysis

Our analysis closely follows that of Xiao et al. (). We acquired Sentinel-1 pre- and postevent single-look complex data from the Alaska Satellite Facility on several ascending and descending paths. Our analysis uses data from ascending path 7, frame 180 with an acquisition from 19 July 2021 as base image, and a repeat on 31 July 2021. We also use acquisitions on the same dates on descending path 102, frames 402 and 407. We perform two-pass processing with GMTSAR (). Atmospheric conditions are mixed; some regions show excellent phase coherence and clearly earthquake-induced phase gradients, whereas others are obviously affected by a thick cloud cover. Because the rupture area is mostly covered by ocean with only a few islands and the Alaska Peninsula in the near field, we unwrapped islands with favorable phase signatures individually using SNAPHU (). The resulting LOS displacement fields are floating in that they are not connected to an area of zero motion; hence, we cannot know the absolute displacements of these patches, unless we have static GNSS offsets that can be used as ties to the InSAR LOS displacements. As neither Kodiak Island nor the nearby Trinity Islands showed significant static deformation in the continuous data, we could only tie Chirikof Island and a part of the Alaska Peninsula to continuous GNSS sites (AC13 and AC40, respectively). The other islands remain floating, and we estimate an ambiguity factor in the inversion, effectively allowing us to fit the phase gradient at these locations. We uniformly downsampled the LOS observations, retaining a total of 2258 InSAR data points. A complete list of islands that we used in the analysis is given in table S1.

Static GNSS analysis

To estimate daily station positions in ITRF2014, we use the GIPSY-OASIS goa-6.5 software () developed by NASA’s Jet Propulsion Laboratory (JPL). We use JPL’s nonfiducial (i.e., loosely constrained) orbit and clock products. Antenna phase models are provided by the International GNSS Service (IGS14) (). We apply tropospheric corrections using VMF1GRID (, ). Ocean tidal loading corrections are based on the TPXO 7.2 and ATLAS model () as implemented in SPOTL (). Once all data from stations in and around Alaska are combined, we transform them into the ITRF2014 reference frame through a seven-parameter transformation. The calculated coseismic offsets are listed in table S2.

High-rate kinematic GNSS analysis

We use the GipsyX analysis software developed at JPL () to estimate 1-Hz kinematic position time series for all near-field GNSS stations. We use IGS antenna phase models (IGS14) () and JPL orbit and clock products in the position analysis. To correct tropospheric delays, we apply the Global Temperature and Pressure (GPT2) model (). Ocean tidal loading corrections are estimated using the SPOTL software (), within which we use the TPXO 7.2 and ATLAS model (). In a last step, we rotate the resulting high-rate precise point positions from an ITRF2014 reference frame into a local east-north–vertical coordinate system. The GNSS waveforms are shown in fig. S2.

Finite fault inversion

The method to estimate time evolving slip on a fine fault from seismic and geodetic data is described by Xiao et al. () who expand a nonlinear wavelet-based inversion routine () to estimate parameters that capture missing phase cycles of the floating InSAR patches.

The Chignik model fault geometry is a local average of Slab2.0 () and thus does not include any curvature. Consequently, the dip angle is depth-averaged from 15 to 40 km, and the strike is not exactly parallel to the trench. We use a modified USGS hypocenter that has the depth fixed to 26.2 km to match the Slab2 depth at that location. Similar to the inversion used in the previous Simeonof inversion (), we align a fault plane of 320-km length and 160-km width, subdivided into 10 km by 10 km patches, with the hypocenter using Slab 2.0 strike (235.2°) and dip (15.1°) angles. We calculate Green’s functions for each subfault using the one-dimensional layered regional velocity structure from LITHO1.0 ().

To infer the kinematic rupture model, we simultaneously estimate rupture initiation time, slip amplitude, rake angle, and an asymmetric rise time function (, , ) for each subfault. We allow a total slip amplitude on each subfault between 0 and 8 m, rake angles range from 60° to 120° in increments of 3°, and the rupture velocity is constrained to between 1.0 and 4.0 km/s. Because we constrain the start and end times of the asymmetric rise time function from 0 to 10 s in 1-s intervals, the slip duration on each subfault is limited to 20 s.

To determine the relative weights of each data type against each other, we tested different weighting schemes using only the indicated data types in the inversion and calculated the model misfit relative to the data for GNSS and InSAR data (fig. S9A), teleseismic waveforms (fig. S9B), and static geodetic data (GNSS and InSAR) against waveform data (high-rate GNSS and teleseismic waveforms; fig. S9C). This is identical to the approach of our Simeonof inversion (). As the misfit for one data type decreased, it increased for other data types, and we picked an optimal weighting to remain within the data uncertainty of each data type.

To avoid unphysical behavior in slip amplitude and rupture propagation, we impose spatial and temporal L2 regularization. We pick the optimal regularization parameters where the L-curve (fig. S10), representing the trade-off between model misfit and model roughness, achieves maximum curvature.

The model for the 2020 M7.6 Sand Point earthquake followed the same methodology, except that only static GNSS and teleseismic waveforms were used. The model fault had a strike of 350° and dipped 49° to the east. The fault had a length of 120 km and a width of 70 km. The model fault plane had 12 subfaults along strike and seven subfaults in the dip direction. Each subfault had dimensions of 10 km by 10 km.

To evaluate model resolution for our Chignik earthquake estimate, we created two checkerboard input models with alternating patches of 0- and 3-m slip for 60 km by 60 km patches (test 1) and 40 km by 40 km patches (test 2). Synthetic data were generated for all the datasets with the rake angle set to 90° and the rupture velocity set up 2.5 km/s. A symmetrical rise time function with a length of 6 s was applied for the simulated seismic and high-rate GNSS waveforms. The resulting symmetrical moment is Mw (moment magnitude) 8.3, very similar to the moment for the 2021 mainshock. To investigate our preferred model’s resolution in the high-slip (≥2 m) area, we performed a test with an input model that approximated the shape of the high-slip area (test III). The input slip was set to 3 m; all other parameters for the simulation were the same as in the checkerboard tests. Results of these tests are shown in fig. S11.

Coulomb stress change modeling

We use the Coulomb 3 Deformation and Stress-Change software (–) to compute stress changes following the earthquakes. We assume an elastic half-space with uniform isotropic elastic properties (). Finite-fault models from the previous Simeonof model () and this work serve as slip and geometry inputs for the computation and drive stress changes on neighboring receiver faults. For the Simeonof and Sand Point stress change models, the receiver fault was the Chignik model fault plane (strike = 235.2, dip = 15.1, and rake = 95). For the Chignik stress change model, the receiver fault was the eastern Kodiak segment of the megathrust (strike = 235.2, dip = 14, and rake = 83). The receiver faults are reverse faults, and the associated fault parameters are chosen on the basis of both geometry of the subducting slab (derived from Slab2) () and fault parameters of previous earthquakes in the same region (e.g., the 1964 M9.2 earthquake). Two coefficients of friction for all scenarios were tested: 0.1 and 0.4. For the vertical cross sections, stress changes were calculated for depth intervals of 1 km and then interpolated within the Coulomb 3 routine. To project stress changes onto the megathrust, stress changes were calculated for depths of 1 to 50 km at intervals of 1 km, and values spatially closest to the Slab2 model interface were extracted and plotted. The python interpolation package KDTree was used to project the depth slices calculated by Coulomb 3 onto the model slab surface.