Revisiting the 1899 earthquake series using integrative geophysical analysis in Yakutat Bay, Alaska, USA

The Alaska-Aleutian subduction zone undergoes frequent subduction-related earthquakes, including the second-largest earthquake ever recorded, the 1964 Great Alaskan M_W 9.2 earthquake. The 1964 earthquake ruptured an 800 km length of the eastern Aleutian megathrust, generating a trans-Pacific tsunami and up to 10 m of uplift at Montague Island (Plafker, 1969). At its eastern limit, the 1964 earthquake also ruptured part of the Yakutat–North American plate boundary where the Yakutat terrane travels with the Pacific plate (e.g., Eberhart-Philips et al., 2006; Chapman et al., 2014; Finn et al., 2015) and is both subducting beneath and colliding with North America (e.g., Worthington et al., 2012). Paleoseismic evidence indicates that the eastern Aleutian and adjacent Yakutat plate boundaries may have also ruptured simultaneously ~900 and ~1500 yr ago, perhaps generating larger “super ‘64” events (Shennan et al., 2009). This scale of a combined plate-boundary rupture would have a rupture area 15% larger than the 1964 Great Alaskan earthquake (Shennan et al., 2009), leading to increased shaking and tsunami hazards. A rupture of this size would affect the local town of Yakutat (population ~650) and coastal population centers all throughout south-central and southeastern Alaska.

The 1964 earthquake and tsunami have been studied extensively (e.g., Plafker, 1969; Shennan et al., 2009; Parsons et al., 2014; Finn et al., 2015; Haeussler et al., 2015), but another series of large earthquakes that occurred near Yakutat Bay, Alaska, in September 1899 is far less understood and equally critical for characterizing and forecasting large plate-boundary ruptures in Alaska. The largest two events of the 1899 series, M_W 8.1 and M_W 8.2, occurred on 4 and 10 September, respectively (Tarr and Martin, 1912; Abe and Noguchi, 1983; Plafker and Thatcher, 2008). The 1899 events occurred near the complex structural syntaxis between Yakutat terrane underthrusting/collision with North America and the Fairweather strike-slip fault system. The Fairweather fault is seismically active, and its largest recorded earthquake (M_W 7.8) occurred in 1958 near Lituya Bay, causing a 524-m-high tsunami there (Miller, 1960; Doser, 2006; Rollins et al., 2021). The M_W 8.1 event of 4 September 1899 is consistent with underthrusting of the Yakutat terrane because the epicenter was near the Yakutat–North America deformation front, the Pamplona zone (Fig. 1; Coffman et al., 1982; Doser, 2006). It is less clear, however, which faults actually slipped during the larger, more poorly located M_W 8.2 10 September 1899 earthquake (Fig. 1; Plafker and Thatcher, 2008). Eyewitness reports of shaking and an ~6-to ~12-m-high tsunami, as well as coseismic uplift locally measured at over 14 m, suggest that epicentral fault slip during the larger M_W 8.2 main shock of 10 September was located in the region of Yakutat Bay and/or its northern extension, Disenchantment Bay (Fig. 1; Tarr and Martin, 1912). There is evidence that shaking effects of the 10 September event may have stretched well beyond Yakutat Bay along the coast of western North America, with seiching observed as far south as Lake Chelan, Washington (Coffman et al., 1982). Existing teleseismic data cannot reduce the large uncertainty on the 10 September event location or better constrain the focal mechanism (Doser, 2006), so it is possible the 10 September event occurred along the Yakutat terrane underthrusting interface, that slip was limited to shallow, thin-skinned fault systems, and/or that strike-slip faulting along the Fairweather system played a role.

In our study, we aimed to better constrain regional fault geometry and structure in the Yakutat Bay region by compiling various geological and geophysical data, assessing data set reliability, and providing an updated fault model based on our analysis. Available data prior to this study included uplift measurements made in Yakutat Bay in 1905 (Tarr and Martin, 1912) and between 1967 and 2000 (Plafker and Thatcher, 2008), global navigation satellite system (GNSS) data modeling (Elliott et al., 2010, 2013; Elliott and Freymueller, 2019), satellite imagery (landsat.usgs.gov), bathymetric and topographic data (National Oceanic and Atmospheric Administration [NOAA], www.ncei.noaa.gov; Alaska Division of Geological & Geophysical Surveys [DGGS], elevation.alaska.gov; Goff et al., 2012), thermochronology results (Enkelmann et al., 2015; Falkowski et al., 2016; Schartman et al., 2019), and basin-scale seismic reflection data from academic survey EW0408 (Gulick, 2004; Gulick et al., 2007; Elmore et al., 2013). We also examined 2012 high-resolution multichannel seismic (MCS) and compressed high-intensity radar pulse (Chirp) seismic reflection data (Zurbuchen et al., 2015; Gulick, 2016a, 2016b) for active structures and evidence for long-term deformation. Using these combined data sets, we: (1) quantified uncertainties in previous coseismic uplift models, (2) improved active fault maps in Yakutat Bay, (3) created a conceptual tectonic model of faulting on the southeastern side of Yakutat Bay, (4) reevaluated previous hypotheses that the 10 September event was related to Yakutat terrane underthrusting, and (5) highlighted the hazard of a similar, future event in the Yakutat Bay region.

The geology of coastal southeastern Alaska consists of a number of fault-bounded accreted terranes, including the Alexander terrane, the Chugach terrane, and the Yakutat terrane, which is primarily an oceanic plateau with a strip of accretionary prism rocks along its eastern edge (Plafker, 1987; Pavlis et al., 2004; Christeson et al., 2010; Worthington et al., 2012). The Yakutat terrane is currently in the process of subducting beneath and colliding with North America (e.g., Plafker et al., 1994). The terranes are generally composed of highly deformed deep-sea sedimentary and igneous rocks; geologic units observed in outcrop and sampled by drilling include the Cretaceous Yakutat Group, the Eocene Kulthieth Formation, the Oligocene–Miocene Poul Creek Formation, the Miocene–Pleistocene Yakataga Formation, and some volcanic Yakutat terrane rocks (e.g., Plafker, 1987; Plafker et al., 1994; Perry et al., 2009; Chapman et al., 2012; Van Avendonk et al., 2013). Rapid sedimentation due to active orogenesis and glacial erosion has deposited an extensive cover of unconsolidated Pleistocene–Holocene clastic sediments throughout the region (Gulick et al., 2015). In Yakutat and Disenchantment Bays (Fig. 1), the shallowest ~500-m-thick section of sediment was deposited by the Hubbard and Malaspina Glaciers during the Little Ice Age, which had two phases with glaciers culminating at ca. 1300 A.D. and ca. 1700–1791 A.D. (Calkin et al., 2001; Elmore et al., 2013; Zurbuchen et al., 2015). Because of glacial activity, sedimentation rates in Yakutat and Disenchantment Bays are on the order of ~1 m/yr locally (Goff et al., 2012).

Southeastern Alaska tectonics are generally dominated by right-lateral transform motion along the Fairweather–Queen Charlotte fault and flat-slab subduction and underthrusting/collision of the Yakutat terrane. The Yakutat terrane travels with the Pacific plate, separated from it at its southwestern edge by a mostly inactive Transition fault (Gulick et al., 2007, 2013). It is undergoing flat-slab subduction beneath North America due to its thickness and buoyancy (Abers, 2008; Worthington et al., 2012; Gulick et al., 2013). The Pamplona zone is the Yakutat terrane–North American plate deformation front offshore, which links up to the Alaska-Aleutian megathrust (the Pacific–North American plate boundary) to the southwest (Fig. 1; e.g., Pavlis et al., 2004; Gulick et al., 2007, 2013; Worthington et al., 2008). To the northeast, the Pamplona zone links to the St. Elias thrust belt, which includes the Foreland fault zone, the Malaspina fault, and the Chaix Hills fault (Fig. 1; Elliott and Freymueller, 2019; Pavlis et al., 2019). These major foreland thrust fault systems have been mapped in the Chugach–St. Elias Mountains in geologic field studies (Plafker, 1987; Bruhn et al., 2004; Chapman et al., 2012; Pavlis et al., 2012) and by using glacial ice-flow patterns (Bruhn et al., 2012; Cotton et al., 2014). The offshore-onshore transition thus marks the approximate location of a shift from Yakutat terrane flat-slab subduction offshore to underthrusting and collision with North America onshore, which is causing active orogenesis in the Chugach–St. Elias Mountains, the highest coastal mountain range in the world (Mount Logan elevation = 5959 m; Berger and Spotila, 2008; Enkelmann et al., 2009, 2015; Pavlis et al., 2012; Worthington et al., 2012; Gulick et al., 2015; Falkowski et al., 2016; Bootes et al., 2019).

The onshore-offshore Fairweather–Queen Charlotte fault is the 900-km-long Pacific-North American transform plate boundary, extending along the continental margin of western Canada and southeastern Alaska. The offshore Queen Charlotte fault, which accommodates ~50–57 mm/yr of right-lateral slip along a single narrow fault zone (Brothers et al., 2020), links up with its onshore, more transpressive counterpart, the Fairweather fault, near Icy Point, Alaska. The Fairweather fault marks the northeastern boundary of the Yakutat terrane, separating it from North America (e.g., Chase and Tiffin, 1972; Carlson et al., 1988; Fig. 1). The Fairweather fault accommodates ~36–45 mm/yr of right-lateral motion, which increases from north to south (Elliott and Freymueller, 2019). The Fairweather fault continues northward past Yakutat Bay, where it meets faults related to Yakutat terrane underthrusting and collision in a structural syntaxis. There, just northwest of Yakutat Bay, the Fairweather fault effectively ends, leading to several splay faults associated with the termination of a major strike-slip fault (e.g., Bruhn et al., 2012). Transpression at the termination of the Fairweather fault is accommodated along the subparallel Boundary, Yakutat, and previously proposed Otmeloi faults (Figs. 1 and 2; Plafker and Thatcher, 2008; Elliott et al., 2010; Elliott and Freymueller, 2019).

The 1899 earthquakes occurred in a complexly faulted syntaxis that includes both underthrusting/collision and transpression near Yakutat Bay, Alaska. Convergent stress related to Yakutat terrane flat-slab subduction, which reaches rates of up to ~46 mm/yr (Elliott and Freymueller, 2019), is accommodated via underthrusting in northwestern Yakutat Bay along the Esker Creek, Bancas Point, and Chaix Hills faults (Figs. 1 and 2; e.g., Chapman et al., 2012; Cotton et al., 2014). The 1899 sequence and the 1958 Lituya Bay earthquake along the Fairweather fault are the largest recorded earthquakes in the region. Other historical seismicity has clustered along the Yakutat terrane–North American plate deformation front faults, with a few events located along the Fairweather fault; the Yakutat Bay area is otherwise largely devoid of seismicity (Fig. 1).

Proposed fault models in the area have been largely based on structural mapping north and west of Yakutat Bay and modeling of the 1979 M_S 7.2 St. Elias earthquake (Figs. 1 and 2; e.g., Plafker, 1987; Estabrook et al., 1992; Chapman et al., 2012; Pavlis et al., 2019), along with mapping and geophysical modeling more local to Yakutat Bay by Plafker and Thatcher (2008) and Bruhn et al. (2012). The Plafker and Thatcher (2008) fault model (Fig. 2) focused on 1899 earthquake-related deformation and included inferred faults crossing Yakutat Bay (e.g., Shennan et al., 2009; Elliott et al., 2010; Enkelmann et al., 2010). In our marine geophysical analysis, we aimed to verify and update this fault model by providing improved constraints on the nearshore Esker Creek, Bancas Point, Yakutat, Boundary, and Otmeloi fault systems, which are critical for understanding the 1899 events and their relationship to Yakutat terrane underthrusting.

In this study, we utilized several bathymetric data sets, including a set of 2018 Yakutat Bay multibeam surveys (H13069–H13075) collected by the NOAA ship Fairweather, available through the NOAA National Centers for Environmental Information (NCEI; www.ncei.noaa.gov). We resampled the variable 32-to 52-m-resolution NOAA grids to a uniform 50 m grid size using a bilinear algorithm and then mosaicked these seven surveys with blended edges. Where we lacked coverage in Yakutat Bay from the 2018 NOAA ship Fairweather data, we utilized a bathymetric compilation from Goff et al. (2012), which incorporated several sources of high-resolution bathymetry data in Yakutat Bay, including multibeam bathymetry data obtained during survey EW0408, which were collected aboard the R/V Maurice Ewing in 2004, from Newport, Oregon, to Kodiak, Alaska. The Goff et al. (2012) data were previously analyzed by Elmore et al. (2013) and Zurbuchen et al. (2015). We used bilinear resampling of the Goff et al. (2012) 50 m grid to project it into the same spatial reference frame as the NOAA data, and we merged the NOAA data and Goff et al. (2012) data together, superseding the Goff et al. (2012) data set where the newer NOAA data existed to generate a combined Yakutat Bay bathymetric data set with 50 m grid spacing (Fig. 2).

Bathymetric sounding data from 1892 (survey H02159) sparsely cover upper Yakutat Bay and lower Disenchantment Bay and are available through NCEI (www.ncei.noaa.gov); a data report accompanies the data set but provides little information about how the navigation or depth data were acquired. Positioning for the 1892 sounding data was likely done via triangulation to local landmarks, providing a potentially significant source of spatial uncertainty. Depth value uncertainties, in addition to those from the tidal range, may also partly derive from the difficulty that older sounding systems had in avoiding sideswipe from adjacent steep slopes (where features out of plane are recorded); in general, if the seafloor slopes are low, depth errors are minimized. The 1892 bathymetric sounding data were too coarse (~2 km) to provide meaningful gridded images, so we differenced the individual 1892 data points with corresponding bathymetric depths from our merged 50 m gridded multibeam map. Because of the 1892 data-acquisition methods, depth uncertainties are likely high (approximately tens of meters) and variable; thus, it was difficult to quantify uncertainties for the 1892 data or their difference with the merged 50 m grid.

High-resolution onshore topographic data sets included processed interferometric synthetic aperture radar (IfSAR) digital elevation data available through the Alaska Division of Geological and Geophysical Surveys (elevation.alaska.gov). The 5-m-resolution IfSAR data were collected by Fugro EarthData, Inc., between the years 2013 and 2016 and are available as standard digital surface models (DSM) and as bare-earth digital terrain models (DTM); here, we utilized the DSMs in order to observe the most realistic representation of the geologic surface. We also use Landsat color satellite imagery (landsat.usgs.gov) with 30 m resolution to inform our onshore interpretations (e.g., Fig. 2).

In order to map offshore fault structures, a high-resolution seismic reflection survey took place aboard the U.S. Geological Survey (USGS) R/V Alaskan Gyre in August 2012, funded by the USGS Earthquake Hazards Program (Fig. 2; here we refer to the survey as AG1208, though it is the same as USGS identifier 2012–602-FA). The survey targeted subsurface nearshore areas in Yakutat and Disenchantment Bays that had not been previously imaged, and where faulting was proposed based on geologic uplift measurements (Tarr and Martin, 1912; Plafker and Thatcher, 2008). The 2012 MCS reflection data included ~153 line km of two-dimensional (2-D) seismic profiles shot using a Sercel mini-GI gun (15/15 cubic inches, or 246/246 cm³), which has a peak acoustic frequency of ~250 Hz, recorded on a 24 channel, ~100 m receiver streamer with group spacing of 3.125 m. Vertical resolution of the data is generally 1–3 m, and the total record length is 2 s. Processing was completed to poststack time migration. Further processing details can be found in Zurbuchen et al. (2015), who first published the data as part of an analysis of Hubbard Glacier advance-retreat history and glacigenic sedimentary sequences. Processed AG1208 MCS data are available through the Lamont-Doherty Observatory Marine Geoscience Data System (MGDS; Gulick, 2016a, https://doi.org/10.1594/IEDA/500076). Coincident Chirp data were collected along most MCS lines during the AG1208 survey using a pole-mounted Knudsen 3.5 kHz Chirp unit and along additional lines in Yakutat Bay when ice coverage inhibited deploying MCS gear (Fig. 2). AG1208 Chirp data are also available through MGDS (Gulick, 2016b, https://doi.org/10.1594/IEDA/500077).

Several lower-resolution, deeper-penetrating seismic reflection data sets also exist in or near Yakutat Bay, including a 2004 R/V Maurice Ewing dual 45/45 cubic inches (737/737 cm³) GI-gun survey (EW0408; for survey details, see Gulick, 2004; for processing, see Gulick et al., 2007; see also Fig. 2), a 2008 R/V Marcus G. Langseth 36 air-gun 6600 in³ (108,000 cm³) survey (MGL0814; for processing and survey details, see Worthington et al., 2010), and a 1979 Western Geophysical survey (W1279EG; for survey details, see Elmore et al., 2013). These lower-resolution MCS surveys, particularly line 2101 from the EW0408 survey, were primarily used for planning the AG1208 survey and to provide regional subsurface context for mapped structures and sedimentary deposits. For this study, we also utilized 2010 Chirp envelope data collected aboard the R/V Quest in Russell Fiord (Goff et al., 2012). The EW0408 MCS data and the 2010 Chirp data shown here are both available via MGDS (https://doi.org/10.1594/IEDA/500105 and https://www.marine-geo.org/tools/search/entry.php?id=HUBCHIRP2010, respectively).

Two sets of historical geologic uplift measurements have been taken along the coastlines of Yakutat and Disenchantment Bays. The first set was measured by Tarr and Martin in 1905 and published in 1912 (hereafter, the “TM12” data set; Tarr and Martin, 1912). Their regionally extensive set of measurements covered the northwest and southeast sides of Yakutat Bay, where they noted and photographed areas of uplift and subsidence. Their detailed report also included remarkable interviews with local inhabitants and eyewitness accounts of earthquake shaking and the 10 September 1899 tsunami. In their field work, Tarr and Martin (1912) utilized elevation markers including exposed barnacles and mussels, submerged trees, and uplifted beaches, sea caves, and sea cliffs (e.g., Evelpidou and Pirazzoli, 2015). As their field work was completed shortly after the 1899 events (in 1905), measured uplift was interpreted to be coseismic with the 1899 events, although it may have included a small amount of afterslip. The other set of uplift measurements was taken by Plafker and Thatcher (2008) (hereafter, the “PT08” data set) in several field seasons between 1967 and 2000, in which they revisited many of the TM12 sites and attempted to corroborate the uplift values. The majority of the PT08 field data were collected in 1973 and 1980; many of their reported values are comparable to or the same as the TM12 data set. Direct comparison of the measurement methodologies between PT08 and TM12 is challenging due the differing and occasionally sparse level of detail about their acquisition in the source studies.

Geodetic and thermochronology data are available for the Yakutat Bay region. GNSS data and tectonic models have been published by Larsen et al. (2005), Elliott et al. (2010, 2013), and Elliott and Freymueller (2019). The most recent Alaska tectonic block model is constrained by a total of 321 GNSS stations (Elliott and Freymueller, 2019). Tectonic and glacial isostatic adjustment (GIA) models fit to GNSS data by Elliott et al. (2010, 2013) informed our understanding of the regional fault geometry and kinematics throughout southeastern Alaska. Additionally, low-temperature thermochronology, including apatite and zircon ages, provided exhumation rates along fault systems around Yakutat Bay, which could be used as a proxy for active uplift and relative fault motion (Enkelmann et al., 2015; Falkowski et al., 2016; Schartman et al., 2019).

Seismic reflection, bathymetry, topography, and satellite imagery data were plotted using ESRI's ArcGIS Desktop software (www.esri.com) alongside contoured GNSS GIA uplift data from Elliott et al. (2010, 2013), onshore geologic uplift measurements (Tarr and Martin, 1912; Plafker and Thatcher, 2008), and previous fault models (Plafker and Thatcher, 2008; Bruhn et al., 2012) to assist with remapping fault geometry. Coulomb modeling was completed using Coulomb 3.3 software available through the USGS (Toda et al., 2011). AG1208 MCS data were processed using Paradigm Geophysical Echos seismic processing software and imported into Landmark DecisionSpace Desktop (www.halliburton.com) for visualization and interpretation alongside other seismic reflection data sets.

We analyzed new and existing geophysical data, including gridded multibeam bathymetry and high-resolution seismic reflection data, for evidence of active faulting and deformation in Yakutat Bay. We also reexamined coseismic uplift measurements from Tarr and Martin (1912) and Plafker and Thatcher (2008) (TM12 and PT08) and quantified uncertainty associated with each data set. Finally, we placed our observations in the context of other recently published data, including geodetic models and thermochronology.

The merged 50 m Yakutat Bay bathymetric grid (Figs. 2 and 3) highlights glacial bedforms on the seabed of Yakutat Bay in both topographic (Fig. 2) and slope (Fig. 3) map representations. The modern glacial geomorphology (Fig. 3) is a result of Little Ice Age activity (culminating in phases between ca. 1300 and ca. 1791 A.D. in this area) from the tidewater Hubbard Glacier, which is still actively calving at the northernmost extent of Disenchantment Bay, and the piedmont Malaspina Glacier (Zurbuchen et al., 2015). Glacial bedforms in Yakutat Bay include a prominent glacial trough located within Disenchantment Bay (Hubbard Channel) and along northeastern Yakutat Bay (Knight Island Channel), each flanked by lateral moraines (Fig. 3). The southern extent of Hubbard Channel has a terminal moraine (the Blizhni Point moraine) located at the mouth of Disenchantment Bay at Blizhni Point (Fig. 3), which has been associated with the most recent Little Ice Age advance of Hubbard Glacier ca. 1700–1791 A.D. (Zurbuchen et al., 2015). Toward the mouth of Yakutat Bay, moraines associated with a past advance of Malaspina Glacier have also been preserved (Fig. 3; Zurbuchen et al., 2015). Several previous studies have gone into further detail about the glacial bedforms and the stratigraphy of Yakutat and Disenchantment Bays (Goff et al., 2012; Elmore et al., 2013; Zurbuchen et al., 2015).

The gridded bathymetric data were analyzed for evidence of faulting in Yakutat Bay, including fault scarps, offset bedforms, and piercing points (Figs. 2 and 3); however, neither the Goff et al. (2012) data set nor the 2018 Fairweather data set shows evidence for these features (Figs. 2 and 3). We acknowledge that the 50 m resolution of the gridded seafloor data may be insufficient to resolve finer-scale features, and that interpreting structural features within the glacial geomorphology, particularly moraines, can also be challenging. However, we corroborated our seafloor observations with subsurface data (see “Subsurface Data Analysis” section below).

The 1892 bathymetric sounding data (Fig. 4) provided a unique opportunity to examine the offshore bathymetric expression prior to the historical 1899 earthquake events. The difference between the 1892 sounding data points and the 2018 grid has the potential to show signals of sediment accumulation, the tectonic response to the 1899 events, and glacial isostatic uplift in the intervening years. The difference between individual 1892 soundings and the same locations on the merged 50 m bathymetry grid (Fig. 4) revealed values ranging from −183 to 137 m. At face value, negative values indicate subsidence or erosion, and positive values indicate uplift or deposition. However, positioning errors in the 1892 data may be the cause of some of the depth differences. We have no way to quantify the errors in the 1892 data but note that locations near steep slopes are more likely to be affected by positioning and depth uncertainties than areas of flat seafloor. The wide range (−183 to 137 m) in differencing depth results suggests that the 1892 sounding depth uncertainty is likely as high as tens of meters, especially near slopes. The uncertainty of the sounding depths eclipses any uncertainty that might be due to GIA (~1–2 m), so we did not apply the effects of GIA to the offshore sounding data.

The comparison of the difference between the 1892 sounding values and modern bathymetry revealed some regional patterns. There is a cluster of low values within the glacially carved Knight Island Channel, a cluster of high values near Knight Island, and generally higher values in Disenchantment Bay closer to the Hubbard Glacier's calving front (Fig. 4). Most positive values in the center of Yakutat Bay are consistent with sedimentation rates in the bays over the last century, as we expect up to ~100 m of deposition since 1892 in the glacial channel near Hubbard Glacier (assuming ~1 m/yr; Goff et al., 2012). Anomalously low values in the channel (e.g., −180.9 m in Knight Island Channel; see Fig. 4) might be artifacts due to steep slopes along the channel flanks, and high values near Knight Island (Fig. 4) are likely to be similarly anomalous due to the presence of land and/or lateral mismatch between data sets due to location uncertainties in the 1892 data.

We observed evidence of several known faults onshore using recent topographic data (5 m IfSAR data and 30 m Landsat satellite imagery) and prior fault interpretations (e.g., Bruhn et al., 2004; Plafker and Thatcher, 2008; Bruhn et al., 2012). On the west side of Yakutat Bay, we interpreted the surficial trace of the Esker Creek fault to be located west of Bancas Point along the range front of the St. Elias Mountains where the change in slope was highest (Fig. 3). A deflection in a tongue of glacial ice near the head of Malaspina Glacier observed by Cotton et al. (2014) supports the position of the fault at or near the range front west of Yakutat Bay. East of Yakutat Bay, NW-SE–trending linear valleys indicate the presence of the Fairweather fault and the Boundary fault (Fig. 3); the Boundary fault goes through Russell Fiord, and the Fairweather fault crosses beneath the terminus of Hubbard Glacier. Eastward of Yakutat Bay, the Fairweather fault continues along a valley through the Fairweather Range (the southernmost range of the St. Elias Mountains; Fig. 1). For the Boundary and Fairweather faults, we only slightly modified previous fault maps (Plafker and Thatcher, 2008; Bruhn et al., 2012) where there was evidence for surficial fault scarps in our topographic data compilation (Fig. 3). Other fault systems east of Yakutat Bay include the Yakutat fault, which we mapped along the base of another range front trending NW-SE (Fig. 3). Our mapping of the Yakutat fault more closely follows the topographic inflection than previous maps of the Yakutat fault (Fig. 3; Plafker and Thatcher, 2008; Bruhn et al., 2012). Plafker and Thatcher (2008) used their uplift measurements to propose an additional thrust fault south of the Yakutat fault, the Otmeloi fault, but we found no topographic or subsurface (see following section) evidence for an additional fault here (Figs. 2 and 3).

Plafker and Thatcher (2008) inferred at least one thrust fault (the Yakutat fault) crossing Yakutat Bay (Figs. 2 and 3) to explain their observations and interpretations of 1899 coseismic uplift patterns. To test this model, we collected and analyzed MCS data throughout Yakutat Bay, Disenchantment Bay, and Russell Fiord in order to characterize any offshore, bay-crossing fault structures associated with the Esker Creek and/or Bancas Point faults, the Boundary fault, the Yakutat fault, and the proposed Otmeloi fault. Our objective was to examine these data for any evidence of offset stratigraphy or deformed regions that might represent a through-going fault structure. Line 2101 from survey EW0408, which goes down the length of Disenchantment Bay and through Knight Island Channel (Fig. 2), became the basis for the later AG1208 survey design because it did not reveal any evidence of faulting or tectonic deformation (Fig. 5A); AG1208 was designed to further test this initial result.

Here, we report little to no evidence of active faults or folds (i.e., offsets, growth folds, deformed stratigraphy) within the upper ~300 m of the subsurface in the AG1208 seismic reflection and Chirp data in Yakutat or Disenchantment Bays (Figs. 5B, 5C, 6A, and 6B). The seismostratigraphic facies here, described in detail by Zurbuchen et al. (2015), generally consist of acoustically chaotic glacial moraine facies flanking flat-lying glacial retreat facies beneath the glacial channels. We observed some NE-dipping subsurface seismic reflections near the proposed location of the Yakutat fault (Fig. 5B), potentially indicative of a nearshore fault or high-impedance glacial moraine deposits. One potential near-seafloor offset was observed near Bancas Point (Fig. 6A), although this offset could be related to the steep moraine deposits and/or local slope failures, and no offsets were visible on adjacent MCS lines (Fig. 6B). We observed flat-lying, undeformed glacial retreat strata throughout Yakutat and Disenchantment Bays, which consistently onlapped glacial moraine deposits, without evidence of active faulting and deformation (Figs. 5 and 6). We emphasize that the AG1208 survey was designed to cross fault traces proposed by Plafker and Thatcher (2008), with 14 MCS crossings and 32 Chirp crossings of the proposed Esker Creek/Bancas Point, Yakutat, and Otmeloi fault locations (Fig. 2), yet we did not image any continuous faults, demonstrably active faults, or deformation uniquely attributable to tectonics (e.g., folding, tilting) in the Yakutat Bay data.

Conversely, we did observe potential fault activity in Russell Fiord along the Boundary fault, where a Chirp line from the HUBCHIRP2010 survey crossed the mapped fault orthogonally (Fig. 7). This line, while difficult to interpret due to seismic artifacts, exhibited a notched seafloor and some possible offset reflections and growth strata (wedge-shaped layers that thicken toward the fjord edges) in the shallow sediments, potentially consistent with an active fault (Fig. 7). Because of the difficulty of accessing Russell Fiord due to ice mélange and large-magnitude currents, no other geophysical data (including multibeam bathymetric data) exist in Russell Fiord.

In order to understand the subsurface observations, we must consider the glacial history and the age of the sediments in Yakutat Bay. Recent deposits in Yakutat Bay have been attributed to advances of the Hubbard and Malaspina Glaciers during the early phase of the Little Ice Age (ca. 1300 A.D.), and the terminal moraine and retreat sequence in Disenchantment Bay were deposited by the Hubbard Glacier in the later phase around 1700–1791 A.D. (Zurbuchen et al., 2015). Due to the relatively short time frame (~700 yr) of the shallowest glacial stratigraphy, it is possible that the imaged stratigraphic record does not well exhibit long-term fault deformation. The 1899 events may be the only earthquakes that could have deformed the upper ~200–300 m of the post–1791 A.D. glacial sediments in Disenchantment Bay, and if sedimentation rates are high enough, 1899 earthquake fault offsets may be seismically invisible due to burial by post–1899 A.D. sediments. Seismically chaotic moraine deposits can also make it difficult in some areas to interpret potential tectonic features.

Despite some of the challenges involved with imaging tectonic offsets in young glacial stratigraphy, we contend that the resolution of our data set should allow for imaging of relatively small amounts of offset or uplift (i.e., the 4–14 m observed by Plafker and Thatcher [2008] along shorelines), particularly within the flat-lying glacial retreat sequences. The AG1208 MCS data have a vertical resolution of 1–3 m, the Chirp data have a vertical resolution of 10–20 cm, and sedimentation rates in the bays are very high (up to and occasionally exceeding ~1 m/yr) due to the presence of the tidewater Hubbard Glacier (Goff et al., 2012). Rapid sediment accumulation provides a high-resolution sedimentary succession, which, when coupled with high-resolution subsurface data, provides a subsurface record capable of resolving seismic reflectors with nearly annual temporal resolution in some locations.

Finally, Zurbuchen et al. (2015) concluded that the Little Ice Age advances likely reoccupied previously carved glacial channels and did not erode all previous sediment. We were able to image the upper ~300 m of sediment (Figs. 5 and 6), which includes pre–Little Ice Age stratigraphy. If this record encompasses a long enough time period, we would expect to see offshore evidence for deep-seated and long-term deformation manifested in structures such as fold-and-thrust systems, growth strata, and continuous, regionally mappable, linear or curvilinear faults, none of which we observed. Similarly, with a high enough slip rate on bay-crossing faults, we would expect to see evidence of deformation even in post–Little Ice Age sediment, which was likewise not observed.

The previous fault model in Yakutat Bay was based largely on interpreted coseismic uplift from onshore measurements made in Yakutat Bay following the 1899 events (Tarr and Martin, 1912; Plafker and Thatcher, 2008). Lack of evidence for faults crossing Yakutat and Disenchantment Bays in the offshore seismic reflection data led us to revisit these onshore geologic uplift measurements and assess possible sources of uncertainty, as neither data set quantified uncertainty or the potential influence of GIA on their values. We started by comparing, point-by-point, published PT08 values with digitized TM12 data (Table 1). We also quantified uncertainty on the uplift measurements, notably including the effects of GIA, one of the most important topographic forcing mechanisms in the region. Larsen et al. (2005) and Elliott et al. (2010) modeled the effects of GIA due to ice retreat since the cessation of the Little Ice Age using GNSS data in the region just southeast of Yakutat Bay, where GIA is particularly significant (their models began at 1770 A.D.). Uplift due to GIA reached rates of 3 cm/yr in locations of former heavy ice coverage, particularly the Yakutat Icefield (Larsen et al., 2005). These uplift rates may have significantly affected measurement and/or interpretation of coseismic tectonic uplift in the region, and they are important to consider when evaluating the reliability of geologic uplift observations.

In order to compare the TM12 and PT08 data, we needed to consider the effects of GIA since the 1899 events. We assumed that the TM12 uplift data were collected in 1905 and that the PT08 data were collected in 1980 (to match 2008 publication discussion and air photo dates). For each PT08 data point, we extracted GIA uplift in that location from Elliott et al. (2010) and calculated a cumulative expected GIA value from 1905 to 1980, assuming constant rates (Table 1). Of course, the rate has not been constant over the twentieth century; Larsen et al. (2003) showed some decadal variability. We attest, however, that this method likely underestimates total GIA because GIA rates beginning ca. 1770 A.D. were initially higher than the values we used due to rapid deglaciation immediately following the Little Ice Age (Larsen et al., 2005). Seventy-five years of cumulative GIA (1905–1980) should be close to the difference between the PT08 and TM12 measurements, so we included cumulative GIA as a source of uncertainty for the PT08 data. We utilized the following equation to estimate uncertainty due to GIA in the PT08 data set (column 5 of Table 1):

where U is uncertainty due to GIA in meters, U₁ is uplift measured by Tarr and Martin (1912) in meters at a given location, GIA is the glacial isostatic adjustment rate in m/yr at the given location, 75 (1980–1905) is a constant in years, and U₂ is the measurement by Plafker and Thatcher (2008) at the given location in meters. Note that one Bancas Point site measured in 1905 was not remeasured in 1980, so uncertainty estimates are not provided for this site (Table 1).

Total uncertainty estimates for the PT08 data (right-hand column of Table 1) were calculated using a standard square root of the sum of the squares error propagation calculation combining the above-estimated GIA uncertainty in meters, an assumed 1 m measurement uncertainty for PT08 values, and cumulative geodetic uplift uncertainties in meters from Elliott et al. (2010). For this aggregate calculation, we assumed that total PT08 measurement uncertainty, including field method, tools used, rugged terrane, and relative sea level, was on the order of ±1 m, which is likely a low estimate. Tarr and Martin (1912) estimated their measurement error to be 0.3–0.6 m, which we included in our data table (Table 1) and discussion. The uncertainty estimates of Tarr and Martin (1912) may also be underestimated, though it seems likely that their paleoshoreline measurements were more accurate than those of Plafker and Thatcher (2008) due to fresher exposure in 1905. The GNSS uncertainty is a maximum of ~20 cm when extrapolated across 75 yr, and, though comparatively small, it was included in our uncertainty quantification because it could be determined relatively accurately. The effects of sea-level rise are also relatively small, on the order of ~10 cm total from 1905 to 1980 (Church and White, 2011). This effect is essentially negligible when discussing uplift on the scale of meters and GIA of ~2–3 cm/yr. Similarly, we might also consider that the PT08 data may include more afterslip than the TM12 measurements, but we contend that the overall similarity between the PT08 and TM12 measurements suggests minimal afterslip after the TM12 data were collected in 1905. Due to the difficulty of separating the coseismic and postseismic signals, we assumed that both the PT08 and TM12 data measured the same net coseismic and postseismic motion.

Given the decades between the two studies, the PT08 measurements should be ~1–2 m higher than the TM12 measurements due to GIA (Elliott et al., 2010; Table 1). In comparing the two onshore data sets, several PT08 values (measured ca. 1980 A.D.) nearly exactly matched the TM12 values (measured in 1905 A.D.), where others were ~1–3 m higher (Table 1). At Bancas Point, the site of interpreted 14.4 m of uplift during the 10 September 1899 event, both Tarr and Martin (1912) and Plafker and Thatcher (2008) measured 14.4 m of uplift (Figs. 8A–8C). However, evidence for the 14.4 m of uplift was only measured at a single site in the TM12 data set, and other TM12 Bancas Point uplift values ranged from 10.3 m to 13.5 m (see Tarr and Martin, 1912). The PT08 data included several measurements of 14.4 m at Bancas Point sites (Fig. 8A), and these measurements should all be expected to be ~2 m higher than the TM12 data due to GIA effects (Table 1). Some PT08 values were indeed ~2 m higher than TM12 values, suggesting actual coseismic uplift at those locations closer to the 12–13 m reported by TM12 (Table 1).

We also compared the TM12 and PT08 data along the east side of Yakutat Bay, where the Yakutat fault, Boundary fault, and Otmeloi fault were mapped (Plafker and Thatcher, 2008). Along Logan Beach, a 20 km stretch of coastline along southeastern Yakutat Bay, Tarr and Martin (1912) measured little to no uplift in 1905 except at Point Latouche, and they even measured subsidence in some areas (Fig. 8A; Table 1). Upon revisiting the TM12 sites, Plafker and Thatcher (2008) observed ~3–5 m of uplift along the length of Logan Beach, reclassifying TM12 subsidence as nontectonic slumping perhaps induced by earthquake shaking. Plafker and Thatcher (2008) also observed a linear nearshore “scarp” visible along Logan Beach from the air, now overgrown by trees (Fig. 8D). Below the “scarp” (and closest to the shore), all trees were dated as post–1899 A.D., suggesting rapid uplift of the beach and subsequent growth of the trees along the uplifted shore (Plafker and Thatcher, 2008). The TM12 data showed minor uplift (0.5–2.7 m) in northern Russell Fiord, with the PT08 data set exhibiting similar values (Fig. 8A); higher values in both data sets were consistently measured on the eastern side of Russell Fiord and the mapped Boundary fault (Fig. 8A). In southern Russell Fiord, both the TM12 and PT08 data sets showed generally higher uplift values in the ~2–3 m range (Fig. 8A). In the Discussion, we consider the PT08 measurements that may have been impacted by GIA uplift.

To investigate potential causes of uplift during the 10 September 1899 earthquake, we generated a simple Coulomb stress change model of the 4 September 1899 event to visualize potential stress loading onto nearby fault planes resulting from a rupture on the Yakutat plate boundary interface (Fig. 9). We modeled the 4 September 1899 event as shallow slip along the Yakutat–North American décollement beginning along the Pamplona zone deformation front, consistent with the event location from Doser (2006). We modeled a deformation front along a simplified Foreland fault zone, a décollement dip of 5°, and rupture extending to 10 km depth, following parameters used by Elliott et al. (2013). We opted for this fault geometry instead of the 1979 M_S 7.2 St. Elias earthquake fault geometry (Estabrook et al., 1992; Pavlis et al., 2019) due to the 4 September 1899 epicenter location (Doser, 2006) and the larger 1899 M_W 8.1 magnitude, necessitating a significantly larger rupture patch than the 1979 M_S 7.2 event (Wells and Coppersmith, 1994; Hanks and Bakun, 2002). Although the area, geometry, and slip during the actual 4 September event are poorly constrained, the preferred model we present here requires 6 m of slip to generate a M_W 8.1 earthquake, modeled on a source fault area of ~70 km by ~115 km extending to a depth of 10 km (Fig. 9).

The receiver faults in our model are the simplified Esker Creek/Bancas Point faults (modeled as one structure), Yakutat, Boundary, and Fairweather faults (geometries compiled from this study; see Fig. 9 and Discussion for details). After evaluating several fault geometries, all models consistently indicated stress loading at the edges of the 4 September rupture patch at a depth of 7.5 km (Fig. 9). The largest stress change in the model occurred on the northeast side of the rupture, downdip of the Esker Creek system (Fig. 9). Note that the model we show here used a 15° dip on the Esker Creek fault based on the dip of the Malaspina fault from Elliott et al. (2013), and that this is somewhat shallower than dip estimates of 30°–60° from geologic studies (Plafker and Thatcher, 2008; Bruhn et al., 2012; Chapman et al., 2012; Cotton et al., 2014; Enkelmann et al., 2015). Our model results are consistent with Coulomb models described in Cotton et al. (2014).

Our new observations of offshore deformation and integration with other geophysical data sets were used to create an updated map of fault geometry in Yakutat Bay (Fig. 10), a conceptual structural model (Fig. 10), uncertainty estimates on previous coseismic uplift measurements (Table 1), and a refined understanding of the tectonic framework and coseismic slip involved in the 1899 events. Our fault model (Fig. 10) represents our best understanding of the fault geometry local to Yakutat Bay and is consistent with recent studies (Elliott et al., 2013; Cotton et al., 2014; Enkelmann et al., 2015; Falkowski et al., 2016; Elliott and Freymueller, 2019; Pavlis et al., 2019; Schartman et al., 2019). In our discussion below, we begin by outlining our interpretations of the fault geometry and structure in the context of new geophysical data, geologic uplift data (along with newly quantified uncertainties), and previously published literature. We then discuss our proposed fault model in the context of the 1899 events and outline the modern tectonic framework in terms of potential future geohazards, including recurrence of an 1899-type event.

Due to a lack of new offshore fault observations, our revised structural interpretations rely on onshore fault mapping and our assessment of previously unreported uncertainty in 1899 coseismic uplift measurements. Our calculated uncertainties for PT08 uplift data are significant, particularly in less-uplifted areas such as Logan Beach and Russell Fiord. Assuming the GIA models are accurate for this area (Elliott et al., 2010), we should see values in the PT08 data set consistently ~2 m higher than the TM12 data set due to effects of GIA, but there is no appreciable geographic pattern indicating this systematic increase in the PT08 data (Table 1). Overall, we suggest that the TM12 measurements might generally be considered a truer indicator of 1899 coseismic motion due to the short time interval between the events and the measurements, which shortly followed the 1899 events (~6 yr), and the well-recorded uplift evidence, including raised shoreline, cliffs, sea caves, and dead, uplifted barnacles located at 80% of TM12 sites. However, we are mindful of the difficulty of measuring uplifts exactly, particularly paleoseismic uplifts based on the morphology of raised shorelines and the use of barnacles and mussels as sea-level indicators (Evelpidou and Pirazzoli, 2015). Tarr and Martin (1912) may therefore have overestimated uplifts in some cases, and we point out specific examples of this in our analysis below. Additionally, a vertical reference datum was not reported by either Tarr and Martin (1912) or Plafker and Thatcher (2008), which may be another significant source of uncertainty in both data sets, but one that we are not able to quantify here.

We found no evidence for a thrust fault or contractional deformation crossing Yakutat Bay or Disenchantment Bay in the offshore MCS or Chirp data (Figs. 5 and 6). We acknowledge the possibility of deeper blind thrusts, faults beneath glaciers, or faults within seismically reflective glacial moraine deposits, which are difficult to image in our data set. However, all the proposed offshore faults (Plafker and Thatcher, 2008) would cut across the flat-lying glacial retreat sequences in the center of Yakutat and Disenchantment Bays, which would have been at least partially deposited by 1899. Seismic reflection data should have sufficient resolution to image faulting-related deformation in flat-lying sequences, especially in high-resolution images of sediments present prior to the 1899 earthquakes (Figs. 5 and 6). This is a primary observation from our study. Additionally, there is no bathymetric evidence of coseismic uplift within Yakutat Bay that can be uniquely attributed to the 1899 earthquakes. There are no lineaments or ridges in the 2018 multibeam bathymetry that might indicate the presence of faults (e.g., Fig. 4). Lastly, while the 1892 sounding data have large uncertainties and are thus unlikely to show coseismic uplift on the scale of several meters, we also did not observe any linear patterns in these data and interpreted only expected patterns of glacial sedimentation (Fig. 4); however, the large uncertainties prevented us from using these data for robust tectonic analysis.

All lines of evidence suggest that local fault systems on the northwestern side of Yakutat Bay, namely, the Esker Creek fault and the Bancas Point fault, are onshore active thrust or reverse faults accommodating contractional strain at rates similar to those experienced by other fault strands associated with Yakutat terrane underthrusting (Table 2). Our results are consistent with numerous previous studies that have found dips between 30° and 60° on these faults (Plafker and Thatcher, 2008; Bruhn et al., 2012; Chapman et al., 2012; Cotton et al., 2014; Enkelmann et al., 2015). We infer that the Esker Creek fault lies along the edge of the range front west of Yakutat Bay (Figs. 3 and 10), which supports substantial north-side-up thrust motion across it, with convergent slip rates perhaps similar to the values of ~8–22 mm/yr inferred along other faults in the St. Elias thrust system (Foreland fault zone, Malaspina fault, Yakataga–Chaix Hills fault; Elliott and Freymueller, 2019).

Both TM12 and PT08 geologic uplift data along the Esker Creek and Bancas Point faults exhibit compelling evidence for coseismic uplift during the 10 September 1899 event (Figs. 8B and 8C; Table 1). Based on our comparisons between the PT08 and TM12 data sets, we interpret that Bancas Point was almost certainly uplifted at least 12 m during the 10 September 1899 event. The TM12 data set has 14.4 m of uplift in only one location, suggesting that 14.4 m of uplift was not ubiquitous across the Bancas Point area and may be an anomalously high measurement. Although the PT08 data set includes several 14.4 m values at Bancas Point, post-1899 GIA of ~2 m occurred by the time of the PT08 measurements, making actual coseismic uplift closer to 12–13 m, which is consistent with TM12 measured uplift in most Bancas Point locations (Table 1).

In our revised fault model (Fig. 10), we interpreted that the onshore Esker Creek and Bancas Point faults connect to each other and are likely linked with the St. Elias thrust system (specifically, the Foreland fault zone and/or Malaspina fault) westward via a deeper décollement (e.g., Pavlis et al., 2019). Our interpretation—that the Esker Creek and Bancas Point faults are active onshore faults accommodating Yakutat terrane tectonic motion—is consistent with extensive geologic mapping by Chapman et al. (2012), which supported interpretation of the Esker Creek fault and the neighboring Chaix Hills fault as major onshore thrust systems. We did not modify the mapping results of Chapman et al. (2012) aside from contributing the new observation that the Esker Creek, Bancas Point, and Chaix Hills systems likely exist exclusively onshore or nearshore.

Due to the large magnitude (M_W 8.2) of the 10 September 1899 event, it seems likely that the Esker Creek and Bancas Point faults are related to Yakutat terrane underthrusting, though it is possible that they may instead be thin-skinned faults associated with Yakutat terrane collision (Chapman et al., 2012). Both of these interpretations are consistent with thermochronology data, which indicate rapid exhumation (3–5 mm/yr) along the Esker Creek and Chaix Hills faults since ca. 5 Ma (Enkelmann et al., 2015; Falkowski et al., 2016). We suggest that the Esker Creek and Bancas Point thrust faults are explicitly linked and that they likely ruptured together in 1899. The Bancas Point fault may be a lateral ramp of the Esker Creek fault, a structure common at the end points of thrust systems and capable of producing large uplift in an earthquake event (e.g., the Pitas Point fault system in southern California; Rockwell et al., 2016). The lack of obvious rupture evidence across Yakutat Bay and the 1899 coseismic uplift pattern support a connection between the Bancas and Esker Creek faults near the shoreline.

Our data, combined with past observations, suggest that the Yakutat fault is an onshore or nearshore, steeply dipping (~60°) reverse fault (Table 2). The Yakutat fault was previously mapped as a northwest-striking, shallowly dipping thrust fault located along Logan Beach and connecting across the bay to the Esker Creek fault (Figs. 2 and 3; Plafker and Thatcher, 2008). We prefer a dip of ~60° for the Yakutat fault, which is consistent with recent GNSS (Elliott et al., 2010) and thermochronology (Enkelmann et al., 2015; Schartman et al., 2019) models. Our preferred dip of 60° is steeper than the 30° dip modeled by Plafker and Thatcher (2008) using uplift data, for which we highlight uncertainties (Table 1); however, our analysis does not preclude a shallower 30° dip. The Yakutat fault has been modeled as a vertical plane in GNSS models (modeled as the Foothills fault; Elliott et al., 2010) and in a zone of relatively high transpression, with the latest models predicting an average of 5.4 ± 0.7 mm/yr convergence and 3.7 ± 0.7 mm/yr right-lateral shear across the Yakutat fault (modeled as the Foothills fault; Elliott and Freymueller, 2019). Thermochronology data show relatively rapid recent exhumation (2–3 mm/yr) along the Yakutat fault (Enkelmann et al., 2015; Schartman et al., 2019), again supporting significant reverse motion across this fault. The geometry of the Fairweather fault also supports transpression across the Yakutat fault. As the Fairweather fault steps onshore at Icy Point, it takes an ~20° convergent left bend, and recent results suggest a total of ~20 mm/yr convergence across the Fairweather fault zone here due to its oblique geometry (Brothers et al., 2020). Our results are consistent with previous studies (e.g., Bruhn et al., 2004), suggesting that this convergence is largely accommodated by thrusts and reverse faults such as the Yakutat fault in the St. Elias Mountains.

Although the Bancas Point fault featured the largest observed uplift during the 10 September 1899 event, smaller amounts of uplift (~1–3 m) on the southeastern side of Yakutat Bay also likely occurred. Observation of a raised 1899-attributed paleoshoreline by Plafker and Thatcher (2008) along Logan Beach (Fig. 8D) strongly indicates at least some uplift along Logan Beach associated with the 1899 events. The biggest discrepancy between the TM12 and PT08 uplift data sets is along central Logan Beach, where the TM12 data set shows zero or near-zero values, and the PT08 data set has significantly higher values of 4–5 m (Table 1). The sandy Logan Beach shoreline lacks sea-level reference features, in contrast to the rocky coastlines of Bancas Point (with exposed barnacles and mussels; e.g., Evelpidou and Pirazzoli, 2015), perhaps affecting the TM12 measurements here. Some of the discrepancy between the two data sets can also be explained by GIA. The expected ~2 m of GIA-related uplift could account for nearly half of the measured uplift at PT08 sites along Logan Beach (Table 1). We therefore interpret that 2–3 m of actual coseismic (or postseismic) uplift occurred along Logan Beach due to the 1899 earthquakes. This interpretation is again consistent with either a shallow (~30° dip) thrust or a steeply dipping (~60°) Yakutat reverse fault located along Logan Beach (Fig. 10).

The Otmeloi fault was proposed at a location south of, and subparallel to, the Yakutat fault in southeastern Yakutat Bay by Plafker and Thatcher (2008) (Figs. 2 and 3). However, there is no clear topographic or bathymetric signature of the fault (Figs. 3 and 10), and the only uplift data are within the uncertainty due to GIA (Table 1). We found no evidence for the fault offshore in the bathymetric data, 2 MCS crossings, or 3 Chirp crossings of the original proposed fault trace (Plafker and Thatcher, 2008; Fig. 2A). The Otmeloi fault is also not required to fit geodetic models (Elliott et al., 2010; Elliott and Freymueller, 2019). Therefore, we suggest that the Otmeloi fault does not exist (Table 2), though it may have small (<2 m) cumulative offsets below the resolution of our geophysical data and calculated uncertainties (Table 1).

The Boundary fault has been mapped along the center of Russell Fiord subparallel to the Fairweather fault (e.g., Plafker and Thatcher, 2008; Figs. 2 and 3), a geometry that is consistent with our data. Chirp data from 2010 potentially support the presence of the Boundary fault within Russell Fiord (Fig. 7), but the data do not constrain the dip or sense of motion across the fault. Imaging of the fault shows seafloor notching with some possible offset thin layers, which we interpret to be due to lesser accumulated sediment and lower sedimentation rate in Russell Fiord, which may allow for better preservation of long-term tectonic deformation. We favor the interpretation of the Boundary fault as a steeply dipping (~85°) transpressional fault, which differs from the interpretation of Plafker and Thatcher (2008), in which the Boundary fault is a shallowly dipping (~10°) thrust. Our interpretation is, however, consistent with asymmetric uplift data collected by Plafker and Thatcher (2008), in which the northeastern side of the Boundary fault has generally higher uplift values (Fig. 8), potentially supporting a northeast-side-up oblique fault. The latest geodetic models fit the northern Boundary fault as a near-vertical fault (Elliott et al., 2010) with 2.8 ± 0.9 mm/yr of right-lateral shear and 1.7 ± 0./8 mm/yr of convergence across the fault, consistent with oblique transpression (Elliott and Freymueller, 2019). Moving southward, geodetic models show that the Boundary fault transitions to become transtensional (Elliott and Freymueller, 2019), further supporting a steep dip (~85°) and a larger share of right-lateral as opposed to dip-slip motion. The geodetic model is also consistent with structural mapping of Bruhn et al. (2004), who showed that the Boundary fault has a steep (90°) dip and dextral slip sense where it emerges from Russell Fiord in the south. Thermochronology data from Enkelmann et al. (2015) suggest <1.5 mm/yr of exhumation along the Boundary fault, or around half that observed along the Yakutat fault; thermochronology data collected by Schartman et al. (2019) also support a transpressive interpretation for the Boundary fault with both dip-slip and strike-slip motion across the fault.

The Schartman et al. (2019) study proposed an additional transpressive fault between, and subparallel to, the Boundary fault and the Yakutat fault called the Calahonda fault (Fig. 10). The Calahonda fault goes through a mélange (Hudson et al., 1977), which may have made it difficult to detect previously, so its existence is based primarily on an offset in apatite and zircon cooling age data (Schartman et al., 2019) as well as a NW-SE linear topographic trend in that location (fault dashed in our Figs. 3 and 10). The presence of an additional fault between the Yakutat and Boundary faults could certainly be consistent with our interpretations, though we do not include the fault in our block model (Fig. 10) or discussion because little is known about the relative motion across the fault.

Russell Fiord likely experienced variable 1899 coseismic or postseismic uplift of ~1–2 m, generally smaller than that at Logan Beach. At four Russell Fiord sites included in both the TM12 and PT08 data, two are located in southern Russell Fiord (east of Logan Beach), and two are located in northern Russell Fiord (close to Disenchantment Bay; Fig. 8A; Table 1). At the southern Russell Fiord sites, both the TM12 and PT08 data show uplift (2.2–3.2 m). After considering uncertainties in southern Russell Fiord (Table 1), it is probable there was 1–3 m of coseismic (or postseismic) 1899 uplift there, perhaps caused by uplift of the hanging wall along the Yakutat fault. In northern Russell Fiord, 0.5–2.7 m values were measured (Table 1), and in both the TM12 and PT08 data sets, lower values (0–0.5 m) are clustered on the western bank, and higher values (2–3 m) are clustered on the eastern bank (Fig. 8A). This spatial distribution of values is consistent with a Boundary fault located in Russell Fiord having a steep (~85°) NE dip and its hanging wall on the northeast side of Russell Fiord. Our uncertainty calculations on the 1899 uplift data (Table 1), alongside interpretations of the fault in geodetic and thermochronology data (Elliott et al., 2010; Enkelmann et al., 2015; Elliott and Freymueller, 2019; Schartman et al., 2019; Table 2), lead us to interpret small amounts of 1899 uplift (~1–2 m) along an obliquely transpressive Boundary fault.

In summary, geological and geophysical data support the existence of strike-slip motion along the Fairweather fault, transpression along the Boundary fault, and transpression with a larger component of convergence across the Yakutat fault. Given our combined observations, we suggest that the effective termination of the Fairweather fault just northwest of Yakutat Bay (e.g., Bruhn et al., 2012) is accompanied by regional right-transpression, similar to a horsetail structure (Fig. 10; e.g., Woodcock and Fischer, 1986; Cunningham and Mann, 2007) and the one-sided flower structure described by Enkelmann et al. (2015). Decreasing northward slip on the Fairweather fault in geodetic models also potentially supports its termination near Yakutat Bay, with some right-lateral strain accommodated northward of the termination via diffuse deformation along the newly forming Totschunda-Fairweather connector fault (Fig. 1; Kalbas et al., 2008; Doser, 2014; Marechal et al., 2015; Elliott and Freymueller, 2019). Strike-slip systems often exhibit complex splay behavior at their end points (e.g., Woodcock and Fischer, 1986; Cunningham and Mann, 2007). In our conceptual model, the Boundary and Yakutat faults are splays of the Fairweather fault accommodating transpression via slip partitioning. Thus, we suggest a horsetail-style slip accommodation pattern wherein the Fairweather fault is a vertical strike-slip fault accommodating purely (or near-purely) right-lateral motion, the Boundary fault is a northeast-dipping, near-vertical (~85°) transpressional splay accommodating both right-lateral and convergent motion, and the Yakutat fault is a northeast-dipping (~60°) reverse fault accommodating primarily convergence (Fig. 10). There is no compelling evidence to suggest the existence of the Otmeloi fault.

Across Yakutat Bay to the northwest, fault orientations shift, and the tectonic regime becomes convergent (Fig. 10). On the northwestern side of Yakutat Bay, the tectonics are dominated by underthrusting and collision of the Yakutat terrane, with convergence accommodated along the Esker Creek and Bancas Point reverse faults. During (or shortly following) the 1899 events, uplift of at least 12 m on the Esker Creek/Bancas Point fault system suggests total coseismic slip of at least 12 m, and this is consistent with a Yakutat terrane underthrusting event; on the other side of the bay, the Yakutat fault slipped 2–3 m and the Boundary fault slipped 1–2 m, consistent with 1899 coseismic or postseismic triggering of the terminus of the Fairweather fault system.

Alternatives to our conceptual fault model exist. We did not image a linking, bay-crossing fault structure in our data set (Fig. 5), but the high-resolution MCS data may not be able to image a deep blind structure crossing Yakutat Bay or structures coincident with glacial moraines. If the Esker Creek and Yakutat faults connect across Yakutat Bay beneath glacial moraines (leading to a geometry something like that proposed by Plafker and Thatcher, 2008), multiple faults may have ruptured coseismically or via dynamic triggering in 1899. If structures link across the bay via a décollement, rupture along the deeper interface in 1899 may have led to coseismic slip along linked splays (like the splay faults in the 1964 Alaska earthquake; Plafker, 1969). Complex coseismic uplift patterns have been observed along subduction zones elsewhere (e.g., New Zealand; Clark et al., 2019). However, we note that GNSS models do not support a low-angle interface beneath Yakutat Bay, instead favoring locked upper-crustal thrust and reverse faults west of the bay (that may or may not sole out into a deeper décollement; Elliott et al., 2013) and higher-angle transpressional faults east of the bay, consistent with our proposed model (Elliott and Freymueller, 2019).

Our Coulomb model of the 4 September 1899 event (Fig. 9) supports loading of fault systems west of Yakutat Bay (Esker Creek/Bancas Point), indicating that the 10 September 1899 event either (1) ruptured part of the underthrusting interface, with slip propagating to the surface locally along the Esker Creek and/or Bancas Point faults, or (2) was caused by stress loading and subsequent rupture along these local upper-crustal faults.

If the 4 September 1899 event was related to Yakutat terrane underthrusting, then that event may have loaded an adjacent Yakutat terrane underthrusting-related fault downdip of the Esker Creek fault (Fig. 9), implying the 10 September 1899 rupture occurred on a deep-seated thrust or an extension of the Yakutat–North American décollement. Plafker and Thatcher (2008) modeled the 10 September 1899 event as a multifault rupture with three distinct slip patches on both sides of Yakutat Bay. High uncertainties on their uplift values on the southeastern side of Yakutat Bay (Table 1) suggest that coseismic slip there may have been more limited than previously modeled and that the majority of slip during the 10 September 1899 earthquake occurred to the northwest along the Esker Creek/Bancas Point system, although new rupture modeling is required to verify this hypothesis. As Plafker and Thatcher (2008) proposed, however, it seems likely that the 10 September 1899 event was a multifault rupture to some degree. Recalculated coseismic uplift of 1–3 m on the southeast side of Yakutat Bay along the Yakutat and Boundary faults (Table 1) likely occurred during (i.e., dynamically triggered) or following the 10 September event. If the Yakutat and Boundary faults slipped postseismically after the 10 September event, coseismic movement on the Esker Creek/Bancas Point faults may have loaded these faults; new rupture modeling of the 10 September 1899 event is also needed to verify this hypothesis.

Because our 4 September 1899 Coulomb stress change model indicates elevated stress at the northeastern edge of the rupture patch closest to the Esker Creek system (Fig. 9), which seems to be the case with varying source thrust fault geometries, the 4 September event may have loaded onshore upper-crustal faults (like the Esker Creek and Bancas Point faults) themselves, instead of, or in addition to, loading the underthrusting interface. The localized, large-magnitude slip potentially supports a shallow crustal rupture on thin-skinned faults related to Yakutat terrane collision with North America rather than a rupture connected to a deeper underthrusting interface. Our recalculated value of ~12 m of coseismic uplift is one of the largest ever measured; Rockwell et al. (2016) also measured up to 12 m of paleoseismic uplift at Pitas Point, California, and by comparison, the 1964 Great Alaskan M_W 9.2 megathrust subduction earthquake produced ~10 m of coseismic uplift at Montague Island (Plafker, 1969).

Regardless of whether the 1899 events occurred due to loading of the deeper underthrusting interface or the upper-crustal faults, there is significant risk of another 1899-type earthquake. If we assume constant convergence of ~46 mm/yr across the St. Elias thrust belt deformation front faults from geodetic models (Elliott and Freymueller, 2019), elastic buildup from over 5 m of shortening has accumulated in that zone since 1899, which would be enough to cause a M_W 7+ event (Wells and Coppersmith, 1994; Hanks and Bakun, 2002). More than 5 m of shortening may have built up if the 1899 events did not fully relieve the accumulated stress following the previous Yakutat terrane underthrusting event ~900 yr ago (Shennan et al., 2009).

Apart from earthquake recurrence, there are other remaining questions related to geohazards in Yakutat Bay, especially landslide and tsunami hazards. Movement along the onshore/blind thrust faults mapped in our study does not readily explain an ~6–12 m tsunami observed by eyewitnesses on 10 September 1899 near the mouth of Russell Fiord (Tarr and Martin, 1912). Tsunami runup of ~12 m on Logan Beach was interpreted as likely to be from the same eyewitness-observed wave (Tarr and Martin, 1912). While it is possible that the ~12-m-high 1899 tsunami was generated by uplifting the coast at Bancas Point or elsewhere, we suggest it is more likely the 1899 tsunami was generated due to slope failure. The tsunami was local to Disenchantment and Yakutat Bays, with variable runups and limited evidence for a regional tsunami (Tarr and Martin, 1912), suggesting a local, nontectonic tsunami source. Several other lines of evidence point to ongoing landslide activity in Yakutat Bay, including slumping observed by Plafker and Thatcher (2008), initially interpreted as tectonic subsidence by Tarr and Martin (1912), and evidence of submarine slope failure from our study (Fig. 6A).

Subaerial landslides transported into water are particularly proficient at generating localized large-runup tsunami waves (e.g., 530 m at Lituya Bay in 1958, 193 m in Taan Fiord in 2015; Miller, 1960; Dufresne et al., 2017; Haeussler et al., 2018; Higman et al., 2018), although a subaerial-to-submarine event is perhaps too efficient to match the comparatively small ~12 m 1899 tsunami wave (Tarr and Martin, 1912). Additionally, no terrestrial landslide scars were observed by Tarr and Martin (1912) or Plafker and Thatcher (2008). Therefore, we speculate that submarine slumping near the glacier termini, where glacial moraines and steep fjord-edge slopes exist (e.g., Fig. 6A), might have provided sufficient source material for the observed ~6 to ~12 m tsunami wave in 1899. However, we did not map a submarine landslide event in our data set near the head of Disenchantment Bay, as the ongoing advance of Hubbard Glacier and glacial outburst floods from Russell Fiord would have likely erased geomorphic evidence of 1899 deposition (Cowan et al., 1996; Trabant et al., 2003; Willems et al., 2011; Goff et al., 2012; Zurbuchen et al., 2015). Tsunami modeling is needed to test the possibility of a submarine source.

Our study highlights the fact that more data and modeling are required to adequately map and understand local fault systems and the secondary impacts (e.g., landslides) of earthquakes near Yakutat Bay. An obvious target is to examine the onshore-offshore connections between the offshore Malaspina fault and Foreland fault zone and the onshore Esker Creek fault, which lie beneath the Malaspina Glacier. Understanding the geometry, character, and recent offsets along the Malaspina fault and Foreland fault zone will assist with linking offshore Yakutat–North American underthrusting structures to onshore collisional tectonics in and around Yakutat Bay and Icy Bay, thus improving understanding of the recurrence risk of an 1899-type event. Additional modeling that considers our calculated uplift uncertainties, specifically earthquake slip and Coulomb models of the 10 September 1899 earthquake, would also be valuable to better understand stress, strain, and slip patterns from that event.

Our study synthesized existing geological and geophysical data to examine fault geometries and fault-related geohazards in Yakutat Bay, Alaska, in the context of the 1899 A.D. earthquakes. Below, we summarize our primary results and conclusions.

1. Based on our new observations of comprehensive MCS and bathymetric data alongside previously published data sets, we found no evidence for major active fault systems crossing Yakutat Bay. This finding suggests that the Esker Creek, Bancas Point, and Yakutat faults exist entirely onshore or nearshore. Using these results, we created an updated map of fault geometry and a conceptual structural model of fault systems in Yakutat Bay (Fig. 10).

2. Dextral transpression likely dominates southeastern Yakutat Bay with strain partitioning through a horsetail-type termination of the Fairweather strike-slip fault. Convergent structures dominate in northwestern Yakutat Bay, where the majority of the 10 September 1899 slip occurred. This finding is consistent with the presence of a structural syntaxis in the area, representing a transition between Yakutat terrane underthrusting/collision and Fairweather strike-slip motion/transpression.

3. We quantified uncertainties in interpreted 1899 coseismic uplift measurements for the Plafker and Thatcher (2008) data set, estimating between 1 and 3 m of uncertainty in their values and finding that uplift related to GIA may have influenced some of their interpretations.

4. Considering coseismic uplift uncertainties due to GIA, our analysis suggests that slip along the relatively steep Esker Creek and Bancas Point reverse faults generated as much as ~12 m of uplift at the surface during the 10 September 1899 M_W 8.2 earthquake. Uplift during or following the same event on the southeastern side of Yakutat Bay along the Yakutat and Boundary faults was in the ~1–3 m range.

5. The 4 September 1899 event stress-loaded fault systems west of Yakutat Bay, leading to a 10 September event sourced either (a) along the Yakutat terrane underthrusting interface, with slip propagating to the surface along local faults like the Esker Creek and Bancas Point faults, or (b) as shallow slip on these upper-crustal faults themselves.

6. The 10 September 1899 event may have been a multifault rupture with dynamically triggered coseismic slip along horsetail-splay fault structures on the southeastern side of the bay. Alternatively, uplift on the southeastern side of the bay may have been postseismic, caused by stress loading from the 10 September 1899 event with a potential contribution from glacial unloading.

7. Recurrence risk of an 1899-type event of at least M_W 7.0 in the Yakutat Bay region is likely high, as a minimum of 5 m of unrelieved elastic strain has likely accumulated along the fault systems on the northwestern side of the bay since 1899.

Many thanks go to the captain and science party of the R/V Alaskan Gyre during the August 2012 cruise, particularly ship captain Greg Snedgen and seismic technician Steffen Saustrup. Thanks go to Julie Zurbuchen for processing the 2012 seismic data. For helpful scientific discussion and feedback, we also thank Julie Elliott, John Goff, Brian Horton, David Mohrig, Nick Hayman, Jake Walter, Bobby Reece, John Swartz, Sebastian Ramirez, Kaustubh Thirumalai, and David Walton. Additionally, we appreciate the U.S. Geological Survey for support of the 2012 cruise and providing much of the legacy data, and the computing facilities and systems administrators at University of Texas Institute for Geophysics (UTIG). This project was funded in part by the National Earthquake Hazards Reduction Program (USGS G12-AP20009) and the UTIG Ewing/Worzel Fellowship. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. government. This is UTIG contribution #3871.