Stress Promotion of the 1958 Mw∼7.8 Fairweather Fault Earthquake and Others in Southeast Alaska by Glacial Isostatic Adjustment and Inter‐earthquake Stress Transfer

We assess how recent glacial isostatic adjustment (GIA) and inter‐earthquake stress transfer modulated the state of stress on major faults in Southeast Alaska, and how these effects may have influenced recent large and moderate earthquakes. The Glacier Bay Icefield has lost >3,000 km3 of ice since ca. 1770, with ice thinning up to 1.5 km, and the rest of southeast Alaska has been deglaciating since ca. 1900. The resulting GIA response in the solid earth produces present‐day surface uplift rates of up to ∼4 cm/yr, among the fastest on Earth. This rapid deformation occurs directly alongside and atop the Fairweather Fault, which accommodates 4–5 cm/yr of right‐lateral motion through slip in frequent large earthquakes, recently including a Mw∼7.8 earthquake in 1958. We find that the 1958 earthquake nucleated very close to the location of maximum GIA‐induced Coulomb stress increase on the Fairweather Fault between 1770 and 1958 (estimated at ∼0.2–0.6 MPa), on the only section of fault where GIA‐induced stress changes had involved both right‐lateral shearing and unclamping. The September 10, 1899 Mw∼8.1 Yakutat Bay earthquake also promoted failure on at least the northwest part of the eventual 1958 rupture. In turn, the 1899 and 1958 earthquakes and GIA promoted failure along the St. Elias compressional margin, the site of a Mw∼7.4 event in 1979. While rapid tectonic loading is very likely the main driver of earthquakes in southeast Alaska, we estimate that 23 of 30 instrumentally constrained Mw ≥ 5.0 earthquakes in this region were also promoted by post‐1770 GIA.

. (a) Tectonics, glacial isostatic adjustment (GIA), and recent earthquakes in southeast Alaska. Small colored circles: GPS-derived uplift rates from Hu and Freymueller (2019) (HF19) over 1992-2003 (white edges) and 2003-2012 (black edges). Grayscale squares: cumulative 1770-2012 ice loss in HF19 GIA model. Large colored circles and dark blue arrows: HF19-R model predictions of GIA-related vertical and horizontal motions (Section 4), averaged over 1992-2012. Earthquakes are described in Section 4 and supporting information Text S4. The September 4, 1899 M w ∼8.1 earthquake occurred somewhere west of the September 10, 1899 earthquake near the coast, likely off the figure (Doser, 2006). Black fault traces are from Alaska Division of Geological and Geophysical Surveys Quaternary Fault and Fold (QFF) database (R. D. Koehler, 2013); faint gray fault traces are from Plafker et al. (1994); Connector Fault (cf.) is described in Section 4. Black vector: motion of Pacific Plate relative to North America from DeMets et al. (2010). Gray vector: motion of Yakutat block relative to North America from Elliott et al. (2010) Euler pole. GB: Glacier Bay. YI: Yakutat Icefield. (b) Location of study region. Red trace: Alaska-Aleutian subduction trench. ten Brink et al., 2018). In the northeastern Gulf of Alaska, the Pacific Plate is carrying the ∼30-km-thick Yakutat block (Worthington et al., 2012) north-northwestward along the coast of southeast Alaska and wedging it into and beneath southern Alaska. The Yakutat block moves at the same rate as the Pacific Plate, although in a direction ∼9° more westerly (Figure 1, gray arrow), and its motion past southeast Alaska is mostly accommodated by the right-lateral Fairweather Fault (hereafter FWF), the northward continuation of the QCF (e.g., Elliott and Freymueller, 2019;Plafker et al., 1978). (The Chatham Strait Fault, which lies to the east and connects to the eastern Denali Fault, is thought to be inactive at present J. Elliott and Freymueller, 2019;Plafker et al., 1994]) The QCF and FWF therefore dominate the tectonics of southeast Alaska, are often jointly called the Queen Charlotte-Fairweather fault system (hereafter QC-FWF), and frequently rupture in large earthquakes. In 1958, a M w ∼7.8 earthquake ruptured most or all of the ∼350-km onshore length of the FWF (Doser, 2010;Plafker et al., 1978) and triggered a large rockslide that caused a megatsunami in coastal Lituya Bay (e.g., Miller, 1961). Further south, the QC-FWF ruptured in a M w ∼7.0 event in 1927 (Doser & Lomas, 2000), a M w ∼7.9 earthquake offshore British Columbia in 1949 (e.g., Lay et al., 2013), a M w = 7.6 shock offshore Sitka, Alaska in 1972 (Schell & Ruff, 1989), and most recently the 2013 M w = 7.5 Craig, Alaska earthquake (Ding et al., 2015). A portion of the QC-FWF near Sitka may have also ruptured in a M w ∼7.5 event in 1848 (supporting information Text S1).
On its north side, the Yakutat block is colliding with and subducting beneath southern Alaska at ∼4 cm/ yr along the St. Elias fold-and-thrust belt (e.g., J. Elliott et al., 2013;J. Sauber et al., 1997). This rapid compression produces both the high topography of the St. Elias Mountains (the highest coastal mountain range on Earth) and abundant seismicity with pervasive thrust mechanisms ( Figure 1a). The largest recorded earthquake in the St. Elias, a M w ∼7.4 oblique-thrust event in 1979, likely ruptured a subhorizontal decollement at the base of the fold-and-thrust belt (Estabrook et al., 1992) and was followed by a dense aftershock sequence southeast of the epicenter (Figure 1a, beige circles). A M w ∼5.9 oblique-normal earthquake in 1920 also suggests some variability in the local stress state (Doser & Lomas, 2000). tramlines, and slightly modified by HF19. The Glacier Bay ice history is represented by additional load sources that follow individual time histories derived by Larsen et al. (2005).
The Earth model consists of a 55-km-thick elastic layer overlying a 230-km-thick Maxwell viscoelastic asthenosphere with viscosity 3 × 10 19 Pa · s (and two deeper layers with much higher viscosities). In a three-dimensional search over lithosphere thickness, asthenosphere thickness and asthenosphere viscosity, HF19 found that this combination of the three parameters (coupled with the load history) produced the best fit to interannual uplift rates extracted from GPS data in southeast Alaska over 1992-2012. Although southeast Alaska features a complex juxtaposition of North American, Yakutat and Pacific lithosphere, the region of fastest uplift mostly consists of the Alexander Terrane of North America and the eastern Yakutat block (e.g., Colpron et al., 2007), which have similar seismic velocity structures (Morozov et al., 2001;Worthington et al., 2012), and so we think that a 1D Earth model is likely adequate. The uncertainties in Earth structure nonetheless contribute to the uncertainties in our subsequent modeling results, which are difficult to quantify.
HF19 used TABOO (Spada, 2003), which does not compute stress or allow for the addition of earthquake slip sources, so we convert the HF19 model to Relax 2010b), using a uniform shear modulus of 52.5 GPa and a Poisson's ratio of 0.25 (supporting in formation Text S2). We henceforth refer to this as the HF19-R model. Despite the rheological simplifications in HF19-R, it reproduces present-day uplift rates well across our entire study area (Figure 1), and confirms that the rapid present-day uplift around Glacier Bay is mostly driven by the viscoelastic component of GIA (Figure S3), as inferred by Larsen et al. (2005). The HF19-R, HF19 and preceding models (J. L. Elliott et al., 2010;Larsen et al., 2004;Sato et al., 2011) do not incorporate the effects of mass redistribution from sea-level change, but in southeast Alaska that change is likely a sea-level fall at 5-9 mm/yr (Hill et al., 2011), 1-2 orders of magnitude smaller than the ice thickness change (Figures 2a-2d), and therefore likely has a second-order effect on GIA. The HF19-R model provides a good fit to Larsen et al. (2005)  To compute stress changes and model earthquake slip on faults, we construct fault surfaces for the Fairweather-Queen Charlotte, Chatham Strait-Denali, Connector, and Transition faults and St. Elias thrust faults based on constraints from literature (Text S3). In most subsequent results, we use Relax to calculate the shear, normal, and Coulomb stress change on these fault surfaces at every model timestep, assuming an effective friction of 0.4 for Coulomb stress change (e.g., Harris & Simpson et al., 1998;King et al., 1994;Rice, 1992). An accompanying mapview-based approach for calculating stress changes is used and described in Section 4.5. In many of the models, we include stress changes from coseismic and postseismic deformation of the September 10, 1899 M w ∼8.1 Yakutat Bay and 1958 M w ∼7.8 Fairweather Fault earthquakes. For the September 10, 1899 event, we use the multi-segment slip model of Plafker and Thatcher (2008) (Figure 3), with optional along-strike extensions to the southeast and west, as described subsequently. For the 1958 earthquake, we use the slip model of Doser (2010) interpolated onto our modeled FWF surface ( Figure 3), with slip prescribed as uniform from 0-15 km depth and zero below that. To incorporate these earthquakes into the stress change history, we add these coseismic slip models to the HF19-R model at timesteps corresponding to 1899 and 1958 and model their postseismic deformation as viscoelastic relaxation in the same asthenosphere that is undergoing GIA, so the two dynamically interact with each other. Modern-day GPS uplift rates are unaffected by these postseismic signals ( Figure S4), and so the HF19 Earth model (which was derived from GPS uplift rates) remains robust even when accounting for them. Details of other earthquake locations, focal mechanisms and rupture extents are provided in Text S4.   Plafker and Thatcher (2008) and Doser (2010) (Section 4). Gray question marks denote hypothetical southeastward extension of September 10, 1899 rupture (Plafker & Thatcher, 2008), which we explore in models by laterally extending the southeasternmost slip source with 10 m of slip.

Previous Research
It is well established that GIA can influence seismicity in dip-slip settings (e.g., Adams, 1989;Hampel & Hetzel, 2006;Hetzel & Hampel, 2005;Steffen et al., 2014), but evidence is more scant in strike-slip settings. Johnston et al. (1998) showed that GIA following the Last Glacial Maximum (LGM) in Fennoscandia, the British Isles and Iceland would have in theory promoted strike-slip earthquakes outside the areas of the former ice loads. Muir-Wood (2000) showed how stress changes from removal of a disc load could combine with a strike-slip tectonic stress regime to yield alternating quadrants of enhanced and inhibited seismicity, with potential application to Fennoscandia and the British Isles. Physical evidence of such GIA-promoted strike-slip ruptures, however, is debated (Firth & Stewart, 2000), and potential links elsewhere are similarly ambiguous. James and Bent (1994) and Wu and Hasegawa (1996) inferred that post-LGM GIA does not explain modern-day strike-slip seismicity in eastern North America. Grollimund and Zoback (2001) proposed that post-LGM GIA could have triggered Holocene earthquakes in the New Madrid Seismic Zone if the stress changes were locally enhanced by a hypothesized low-viscosity zone at depth; however, those earthquakes have more recently been linked to stress changes from late Pleistocene erosion (Calais et al., 2010). The largest strike-slip earthquake that has been speculated as linked to GIA, to our knowledge, is the 1998 M w ∼8.1 intraplate Antarctic plate event (Tsuboi et al., 2000), but that event's distance offshore made the influence of GIA difficult to prove or disprove (Kreemer & Holt, 2000). In contrast to Southeast Alaska, these other examples are from low-strain-rate regions.
Surface mass changes have also been linked to strike-slip earthquakes in non-GIA contexts. Luttrell et al. (2007) and Brothers et al. (2011) postulated that fluctuations of prehistoric Lake Cahuilla could have modulated the earthquake cycle on the southern San Andreas Fault. Luttrell and Sandwell (2010) and Neves et al. (2015) assessed how post-LGM sea-level rise could have modulated the stress on the San Andreas Fault and strike-slip faults in Portugal, respectively, although those effects would have occurred over timescales much longer than the rapid post-1770 deglaciation in southeast Alaska. Seasonal hydrological loading has also been shown to influence seismicity in strike-slip regimes in central California (Amos et al., 2014;Johnson et al., 2017;Kraner et al., 2018), New Madrid (Craig et al., 2017), central and southwest Japan (Heki, 2003), and possibly in south-central Alaska (Johnson et al., 2020). However, the load changes in these cases are on the order of <1 m equivalent water height, in contrast to the >1 km cumulative ice loss in Glacier Bay.

Reference Models for Physical Intuition (Figure 4)
To provide intuition into how ice loss and GIA could modulate the stress on strike-slip faults like the Fairweather Fault, we construct generic models of a 100 × 100 km square ice mass losing 1 m/yr equivalent water height atop or alongside a vertical right-lateral strike-slip fault, within the same Earth structure as the HF19-R model for southeast Alaska. We compute the shear, normal and Coulomb stress changes on the fault through 200 years of constant ice loss and elastic-viscoelastic GIA. Shear stress changes are computed for right-lateral slip and we assume an effective friction of 0.4. These models yield several results that inform our subsequent analysis in southeast Alaska. If the decreasing load is centered over the fault (Figure 4a), GIA unclamps the fault (decreases the normal or confining stress on it) because the lithospheric rebound beneath the load includes an outward radial component of motion (e.g., Figure 1a, horizontal vectors), and so the Coulomb stress on the fault is increased by up to 0.6 MPa (6 bars) beneath the load center. If the fault runs along the edge of the decreasing load (Figure 4b), the unclamping is milder but is accompanied by a shear stress change that is antisymmetric about the load, yielding an asymmetric pattern of Coulomb stress change with a maximum Coulomb stress increase for right-lateral slip of 0.5 MPa (5 bars) at the corner of the load. If the fault runs 50 km from the edge of the decreasing load (Figure 4c), the shear stress changes are milder and the fault is clamped shut because it is in the forebulge, where horizontal motions are radially inward rather than outward; therefore Coulomb stress is decreased over more of the fault than it is increased. Finally, if the fault runs along the opposite edge of the decreasing load (Figure 4d), the unclamping is the same as in Figure 4b but the shear and therefore the Coulomb stress change distributions are reversed. These results are generalizable to other cases of GIA in strike-slip settings, as a complement to the models of Muir-Wood (2000).

GIA-Induced Stress Changes and the 1958 Fairweather Fault Earthquake (Figures 5 and 6)
We then use the HF19-R model to compute the history of stress changes imparted by post-1770 ice loss and GIA to faults in southeast Alaska. This yields our main finding: that the location of maximum GIA-induced Coulomb stress increase on the Fairweather Fault as of 1958 was very close to the epicenter of the 1958 M w ∼7.  decreasing load in the example in Figure 4b. These combined effects also mean that the 1958 epicenter was located near the maximum GIA-induced Coulomb stress increase on the FWF regardless of the effective friction assumed ( Figure S8a).
Assuming an effective friction of 0.4 and that a shear modulus of 52.5 GPa (Section 2) applies along faults, the total 1770-1958 GIA-induced Coulomb stress increase at the 1958 epicenter was ∼0.6 MPa (6 bars) at 2.5 km depth,  loss is stopped in 1900 (and thus confined entirely to Glacier Bay) but viscoelastic relaxation from this ice loss continues after 1900. This model ( Figure S9e) yields a similar Coulomb stress change pattern on the FWF as of 1958 to that of the main model (Figure 5a), showing that most of the GIA-induced stress change on the FWF-including the stress increase at the 1958 epicenter-resulted from the rapid 19th-century ice loss in Glacier Bay and the viscoelastic response to that ice loss. Furthermore, the pre-1958 stress changes on the FWF mostly resulted from viscoelastic deformation rather than the elastic response to short-term interannual ice loss (Figures S10e-S10f). Together, these results indicate that the pre-1958 stress changes specifically mostly resulted from the viscoelastic response to 19th-century ice loss in Glacier Bay.

Comparison of GIA-Induced and Tectonic Stress Changes Along the Fairweather Fault
The estimated GIA-induced ∼0.2-0.6-MPa Coulomb stress increase at the 1958 epicenter is likely equivalent to the stress change from only ∼5 years of shear loading on the fast-slipping QC-FWF (Figures S7c and S7d). As such, GIA likely was a second-order contributor to stress on the FWF, with tectonic loading being the main driver. The 1958 surface rupture also included up to 1 m of southwest-side-up dip slip at the surface (Tocher, 1960); this would have been inhibited by pre-1958 GIA ( Figure S8b), as the lithospheric rebound northeast of the FWF in Glacier Bay would have instead promoted northeast-side-up motion on the fault. The SW-side-up dip slip may be more plausibly explainable as a product of NE-directed tectonic compression (J. Elliott & Freymueller, 2019;McAleer et al., 2009;Schartman et al., 2019). This too suggests that tectonic loading is likely the main driver of earthquake behavior on the FWF, although the timing of failure may be influenced by GIA.

The 1958 Rupture and Stress
Changes From the September 10, 1899 M w ∼8.1 Earthquake (Figures 7 and 8) The 1958 earthquake ruptured unilaterally northwestward from the epicenter (Doser, 2010). In addition to increasing the shear and Coulomb stress on the FWF near the 1958 epicenter, pre-1958 GIA likely unclamped most or all of the 1958 rupture ( Figure 5c). However, it would also have decreased the shear and thus the Coulomb stress along the rupture northwest of the epicenter (Figures 5a and 5b), even given a larger effective friction coefficient ( Figure S8a Figures 4b and 4c.) In the Section 5, we discuss two factors that could have potentially counteracted this shear stress decrease (dynamic rupture propagation and time since the last rupture). A third factor that may have promoted northwestward propagation in the 1958 earthquake was the stress change imparted by the 1899 M w ∼8.1 Yakutat Bay earthquakes, particularly the September 10, 1899 M w ∼8.1 event, which likely ruptured multiple thrust or oblique-thrust faults seaward of the FWF (Figure 3). Using the Plafker and Thatcher (2008) dislocation model for the September 10, 1899 earthquake, we calculate that that event increased both the shear and normal stress on the northwest end of the future 1958 rupture (Figure 7). The September 10, 1899 event also may have also ruptured further southeast down the coast (Plafker & Thatcher, 2008), and if so it would have unclamped more of the 1958 rupture ( Figure S11), although the total Coulomb stress effect on the FWF depends on whether slip on the southeast extension was closer to pure reverse ( Figure S9) or right-reverse ( Figure S12). Still, despite uncertainties in the rupture geometry and slip distribution, we conclude that the September 10, 1899 event likely promoted failure on at least the northwest section of the 1958 rupture (Figures 7, S11 and S12). (The September 4, 1899 M w ∼8.1 event was further west and likely had a lesser stress effect, Figure S13). The combination of GIA and the September 10, 1899 event would therefore have promoted failure along much of the 1958 rupture (Figure 8).
The time evolution of Coulomb stress at three points on the FWF in this combined model shows that although the near-field stress change in the September 10, 1899 event may have been very large (Figure 8b, yellow line), stress was rising faster at the 1958 epicenter (red line) than along any other part of the FWF before 1958 (even accounting for postseismic deformation of the September 10, 1899 event ( Figure S13), which is included in the model in Figure 8 as described in Section 4); therefore the epicenter was approaching the Coulomb failure threshold faster than anywhere else on the FWF. This is relevant because: (1) the absolute stress at any point on the fault is unknown, so each of the lines in Figure 8b could be shifted up or down; and (2) the magnitude and distribution of slip in the September 10, 1899 event is highly uncertain. Therefore it is not known how close any point on the FWF was to failure in an absolute sense.

GIA-and Earthquake-Induced Stress Changes on Thrust Faults in the St. Elias (Figure 9)
We find that GIA also increased the Coulomb stress on north-dipping thrust faults of the St. Elias margin (Figures 5a and 9a), in agreement with Sauber and Molnia (2004) and Sauber and Ruppert (2008). In contrast to the FWF, this GIA-related stress change predominantly resulted from post-1900 ice loss overlying these faults (Figures 2b, 2c, S7d, and S7f), which would have unclamped them (Figure 5c (Figure 9c). This depends on its location and extent; the Doser (2006) relocated epicenter is 59.43°N, 143.05°W, and the observation of an uplifted terrace at Cape Yakataga (Plafker & Thatcher, 2008) suggests that the September 4, 1899 event may have been a shallow thrust earthquake like the September 10, 1899 event six days later, but this remains speculative.

GIA's Overall Effect on Earthquakes in Southeast Alaska (Figure 10)
To assess the overall effect of post-1770 GIA on recorded earthquakes in southeast Alaska, we compile a unified earthquake catalog for the region from 1920 to present (supporting information Text S4) and use the HF19-R model to calculate the GIA-induced stress changes at the sites and times of these earthquakes. Within an area extending from 57° to 61° N and from 142° to 132° N, the unified catalog includes ROLLINS ET AL.   To show this graphically and to complement Figures 5-9, we make mapview plots of GIA-related stress changes on a grid of representative fault orientations throughout Southeast Alaska. We derive 2D grids of local fault orientations ( Figure S15) and slip senses ( Figure S16) by spatially interpolating from our fault model (Section 2; we add the Plafker and Thatcher [2008] model of the September 10, 1899 earthquake) and from the strikes, dips and rakes of M w ≥ 5.0 events in the unified catalog, using natural-neighbor interpolation in MATLAB. We then calculate the history of the stress change tensor at 10 km depth and resolve the stress changes on the gridded fault planes and rakes. Figures 10a-10d show the distributions of Coulomb, shear and normal stress changes as of 1958. ( Figure S6 shows an example for a uniform receiver fault orientation based on the strike, dip and rake of the 1958 FWF earthquake; Figure S17 shows the result if the fault orientation grid is interpolated only from major fault surfaces, with focal mechanisms excluded.) When incorporating the coseismic and postseismic deformation of the September 10, 1899 and 1958 earthquakes ( Figure S18), 24 of 30 M w ≥ 5.0 events would have been promoted and six inhibited.
As one test of whether the number of GIA-promoted M w ≥ 5.0 events is meaningful, we use these maps to assess how many of the 30 recorded M w ≥ 5.0 events would have been promoted if they had occurred at random locations within our map area. We run 100 models in which we compute the Coulomb stress change at sets of 30 random points within the map area as interpolated from the distribution in Figure 10a. In these 100 models, the mean and median number of promoted locations are 11.26 and 11 out of 30, respectively, and none of the 100 models has more than 17 promoted locations out of 30, in contrast to the 23 in the real set. (If one interpolates the stress changes from the distribution in Figure S6a, in which the fault orientation grid is based only on major fault surfaces, the mean, median and maximum in 100 are 9/30, 9/30, and 15/30, respectively.)

Factors Promoting Failure and Rupture Propagation in the 1958 Fairweather Fault Earthquake
Our key finding is that post-1770 GIA, in particular the viscoelastic response to 19th-century ice loss in Glacier Bay, likely promoted Coulomb failure at the epicenter of the 1958 M w ∼7.8 Fairweather Fault earthquake. The 1958 epicenter coincides with the largest GIA-related Coulomb stress increase on the FWF ( Figure 5; Section 4.1). Pre-1958 GIA also unclamped most or all of the 1958 rupture zone ( Figure 5c); however, it also decreased the shear and Coulomb stress on the FWF northwest of the epicenter (Figures 5a  and 5b), which in principle should have inhibited the northwestward rupture in the earthquake. However, the September 10, 1899 M w ∼8.1 Yakutat Bay earthquake likely increased the Coulomb stress on at least the northwest end of the FWF, possibly counteracting the effect of GIA there (Figures 7 and 8; Section 7).
Two other factors may have additionally promoted northwestward rupture in the 1958 earthquake. The first is stress transfer during rupture propagation. During a strike-slip rupture, the slip that has occurred should quasi-statically promote shear failure beyond the rupture tip, which could overcome the GIA-induced shear stress decrease. As the 1958 earthquake was not a supershear event (Doser, 2010), shaking from the slip that had occurred at a given time would have also dynamically modulated the shear stress state beyond the propagating rupture tip. In addition, slip on a bimaterial interface (which the FWF likely is) should induce a dynamic reduction in the normal stress (Ben-Zion, 2001;ten Brink et al., 2018;Weertman, 1980), which would have compounded the GIA-induced unclamping ( Figure 5c) and further promoted rupture.
The second factor is that the FWF northwest of the 1958 epicenter may have been particularly "due" for rupture. In the relocations of Doser and Lomas (2000) and the historical record of Sykes et al. (1981), the only two pre-1958 events attributed to the QC-FWF are the 1927 M w ∼7.0 event to the southeast (Figure 5a) and the 1848 M w ∼7.5 event near Sitka (then Novo-Arkhangelsk) (Figure 2a), which may have been similar to the 1972 M w = 7.6 Sitka earthquake (supporting information Text S1). Therefore, it is possible that the 1958 rupture zone had not ruptured since at least 1848 (Sykes et al., 1981) and perhaps in even longer if the 1848 event was further south, whereas the QC-FWF to the southeast had likely ruptured in 1927 and perhaps in 1848. Plafker et al. (1978) and Koehler and Carver (2018) estimate the recurrence interval for earthquakes on the FWF to be 60-110 years and 70-155 years, respectively, consistent with the 110-year interval between 1848 and 1958). However, historical earthquake records in this region may be incomplete, particularly along the remote coast spanned by the 1958 rupture. Previous large waves (local tsunamis) likely occurred in Lituya Bay in ca. 1853-1854, ca. 1874and ca. 1899(Miller, 1961, any of which could have in theory been triggered by a nearby earthquake. (The ca. 1899 wave may have been triggered by the September 10, 1899 M w ∼8.1 Yakutat Bay shock or by a nearby event it triggered [Miller, 1961].) An uncatalogued earthquake was also felt in Lituya Bay on 30 August 1933, and Jim Huscroft, the sole permanent resident of Lituya Bay between 1915 and 1917 and 1939, described earthquakes as "a frequent occurrence in the bay" (Fradkin, 2001). Nevertheless, the instrumental evidence for the M w ∼7.0 shock to the southeast in 1927 (Doser & Lomas, 2000) and the lack of instrumental evidence of any other pre-1958 events on the FWF suggests that the northwestward rupture in 1958 may have been favored by a lack of large recent ruptures on that section of fault, with GIA-induced stress changes perhaps favoring the 1958 epicenter as the first point of failure.
A follow-on question is whether the 1927 M w ∼7.0 event itself could have promoted failure in 1958. Using the Doser and Lomas (2000) rupture extent and modeling coseismic and postseismic deformation with the HF19-R earth model, we infer that the 1927 event likely occurred too far south to have an appreciable stress effect on the 1958 rupture ( Figure S8e), although this depends on the extent of rupture in 1927. If there were indeed a gap between the 1927 and 1958 rupture zones ( Figure S8e), it could also have been filled by a smaller, uncatalogued event (presumably early in the 20th century or before that), aseismic slip, or some combination thereof, any of which could have promoted failure in 1958. We calculate that the 1949 M w ∼7.9 QCF earthquake likely had no stress effect at the 1958 epicenter ( Figure S8e). Poroelastic relaxation of the elastic stress changes during GIA also likely had only a small effect on fault stress over the 1770-1958 period compared to the effect of GIA itself ( Figure S19). GIA also may have promoted failure in the 1927 earthquake ( Figure 5a) and possibly the 1848 event ( Figure S8b; supporting information Text S1), depending on their rupture extents.
Finally, we note that the 1958 earthquake occurred in July, the warmest month of the year in Southeast Alaska (http://oldclimate.gi.alaska.edu/Climate/Location/TimeSeries/Juneau.html; http://oldclimate. gi.alaska.edu/Climate/Location/TimeSeries/Yakutat.html), when seasonal ice loss would have been at its most rapid. 1958 was also the warmest year in the region since 1944 (see same references), which would presumably have further enhanced short-term ice loss. While the magnitude and distribution of this ice loss is unknown, the major topographic feature of the 1958 epicentral region is the high Fairweather Range, which lies north of the epicenter and on the northeast side of the FWF, contrasting with lower, less glaciated topography and the shore on the southwest. If the largest summer snow and ice loss occurs within the Fairweather Range, it would also increase the Coulomb stress at the 1958 epicenter (Figures 4b and 4c). Therefore, the high temperature of 1958 and the month of July could have also promoted failure at the epicenter. (The same would not apply to the 1979 M w ∼7.4 St. Elias earthquake, which occurred in late February, although it could perhaps apply to the September 1899 M w ∼8.1 earthquakes [J. Sauber & Ruppert, 2008].)

Impact of GIA-Induced Stress Changes on Earthquake Probabilities on the Fairweather Fault
In the context of the rapid tectonic loading on the FWF, we can quantitatively assess the significance of the GIA-induced stress change at the 1958 epicenter by estimating its impact on earthquake probabilities. We assume that the GIA-induced stress change has the effect of reducing the FWF's earthquake recurrence interval (e.g., Working Group on California Earthquake Probabilities, 1990), as if the fault were being loaded faster (which may be viable to first order as the stress change occurred gradually). Assuming that the stress change accrued linearly between 1770 and 1958 and was in total equivalent to a 5-year reduction in recurrence interval over that time (Section 4.2), it translates to a 2.7% reduction in recurrence interval. We then assume that FWF earthquakes follow a stationary Poisson process in which the probability at least one event in timespan t given recurrence interval R is 1 -exp (-t/R). Assuming a background recurrence interval of 70 years (R. Koehler & Carver, 2018), the 50-year probability of at least one earthquake on the FWF without GIA would be 1 -exp (-50/70) = 51.1%, whereas with the GIA-induced stress change it would be 1 -exp (-50/(70-0.027 × 70)) = 52.0%, a probability gain of 0.9%. If the background recurrence interval were instead 155 years (R. Koehler & Carver, 2018), GIA would change the 50-year probability of at least one event from 27.6% to 28.2%. These subtle changes underline that the contribution of GIA is likely second-order in the context of the rapid tectonic loading in southeast Alaska. This parallels the finding that aftershock activity is controlled more by the preexisting stress state than by the stress changes in earthquakes (Segou & Parsons, 2020).

Promotion of Failure in Other M w ≥ 5.0 Events by GIA
We find that post-1770 GIA likely promoted Coulomb failure in at least 23 of the 30 M w ≥ 5.0 earthquakes recorded in southeast Alaska since 1920 ( Figure 10), a correlation well above what randomized models would predict (Section 4.5). (Note that GIA's impact on the 1899 M w ∼8.1 earthquakes is poorly constrained as regionwide ice loss is estimated to have begun ca. 1900). This suggests that although rapid tectonic loading is the first-order driver of stress and earthquake activity in this region, GIA may nevertheless play an accompanying role (perhaps promoting earthquake occurrence in time), as inferred by Sauber and Molnia (2004) and Sauber and Ruppert (2008) in the St. Elias. For comparison, the estimated ∼0.2-0.6-MPa Coulomb stress increase at the 1958 Fairweather Fault earthquake's epicenter is well above the estimated ∼0.01-0.02-MPa threshold value needed for static stress changes to impact seismicity (Hardebeck et al., 1998;Rydelek & Sacks, 1999). Although GIA has been frequently linked to dip-slip earthquakes (Section 3), the southeast Alaska example shows that it may also affect strike-slip earthquakes. This may be applicable (Figure 4) to other cases of rapid GIA in proximity to strike-slip faults, as for example in northern Patagonia (Dietrich et al., 2010;Wang et al., 2007) and possibly along the Alpine Fault in New Zealand (e.g., Chinn et al., 2012).

Stress Interactions Between the Fairweather Fault, St. Elias Margin and Connector Fault
The September 10, 1899 M w ∼8.1 thrust earthquake likely promoted failure on the northwest Fairweather Fault (Figure 7), and in turn the 1958 FWF earthquake likely promoted failure on nearby thrust faults of the St. Elias (Figure 9b). Positive feedback loops like these may be characteristic of earthquake interactions in this tectonic corner over the long term and, if so, may assist in the transfer of stress between the FWF and St. Elias. Bufe and Boyd (2008) found a similar feedback loop between the 1958 rupture, the great 1964 M w = 9.2 subduction earthquake to the west, and the 2002 Denali earthquake to the northwest. We note also that the 1958 earthquake likely inhibited Coulomb failure on the postulated Connector Fault ( Figure S7f). While the September 10, 1899 earthquake likely promoted failure on the southernmost Connector Fault (the Art Lewis Glacier Fault), it would have also inhibited failure on it north of ∼60° N (Figures 7, S11-S13). This may have contributed to the Connector Fault's relative seismic quiescence over the instrumental period (Figure 1a), and if the 1958 and 1899 earthquakes are characteristic of typical FWF and great St. Elias events, these stress shadow effects would complicate the transfer of tectonic strain from the FWF to the Connector Fault over the long term.

Conclusion
We find that GIA driven by rapid 19th-century ice loss in Glacier Bay promoted Coulomb failure at the eventual epicenter of the 1958 M w ∼7.8 Fairweather Fault earthquake ( Figure 5, Section 4.1), favoring the epicenter for failure and perhaps affecting the earthquake's timing and/or characteristics. Regionwide GIA driven by recent ice loss promoted Coulomb failure at the sites of at least 23 of 30 known M w ≥ 5.0 earthquakes in southeast Alaska (Figure 10). While rapid tectonic loading is very likely the main driver of fault stress and earthquakes in this region, our findings suggest that GIA may also impact them. This is consistent with the findings of Sauber and Molnia (2004) and Sauber and Ruppert (2008) that ice load can influence seismicity on seasonal and interannual timescales in the St. Elias. The possible influence of GIA on strike-slip earthquakes on the Queen Charlotte-Fairweather Fault complements the cases of GIA's influence on dip-slip earthquakes in the literature. The 1958 earthquake and other events may have also been favored by static stress transfer from preceding earthquakes (Figures 7-9). These two processes may thus play nonnegligible roles in earthquake occurrence in southeast Alaska and perhaps in similar settings elsewhere.

Data Availability Statement
GPS data are provided in Hu and Freymueller (2019).