Paleoseismology of the Sawtooth Fault and Implications for Fault Behavior in the Epicentral Region of the 2020 $Mw$ 6.5 Stanley, Idaho, Earthquake

Regions of distributed, discontinuous, cryptic, or low‐rate faulting present challenges for fault characterizations in probabilistic seismic hazard analysis (PSHA). Complicated fault networks can produce spatial variability in slip modes (Brocher et al., 2017), along‐fault displacement (Manighetti et al., 2015), damage zones (Peacock et al., 2017), and activity rates (Nicol et al., 2009; Spotila and Garvue, 2021), with ruptures influenced by inherited or reactivated structures related to fault growth or previous tectonic episodes (e.g., Improta et al., 2019; Hecker et al., 2021; Nevitt et al., 2023). Earthquakes in these regions can have long‐lasting engineering and socioeconomic impacts but may remain enigmatic owing to a subdued surface‐rupture expression (e.g., Quigley et al., 2016; Koehler et al., 2021) or poorly understood relation to incomplete or simplified source‐fault mapping (e.g., Chiaraluce et al., 2004; DeLong et al., 2016; Hatem et al., 2022). Although difficulties remain in characterizing complex fault networks for PSHA, integrating geologic, geodetic, and seismological data can yield insights into faulting and rupture processes (e.g., Goldberg et al., 2020; Pang et al., 2020).

One such complex fault region is the Centennial Tectonic Belt (CTB) within the northernmost Basin and Range Province (Fig. 1), characterized by moderate‐to‐large‐magnitude seismicity (Stickney and Bartholomew, 1987; Pang et al., 2018), low geodetic strain, and northwest‐trending normal faults that accommodate regional extension and crustal rotation (Payne et al., 2012; McCaffrey et al., 2013). The Eocene Trans‐Challis fault system (TCFS)—a 270 km long extensional zone of northeast‐trending normal faults, grabens, and eruptive centers—forms the northern terminus of the CTB (Bennett, 1986). Along the western border of the CTB, the Sawtooth fault (Fig. 1) bounds a 65 km long east‐facing rangefront (Fig. 1). The Sawtooth fault has evidence of multiple postglacial ruptures (Thackray et al., 2013; Lifton et al. 2023; Shapley et al., 2023); however, their relation to modern seismicity and regional fault kinematics (e.g., Liberty et al., 2020) remains largely unstudied because of the remote nature of the region.

In March 2020, the moment magnitude (⁠$Mw$⁠) 6.5 Stanley, Idaho, earthquake occurred at ∼20 km depth near the northern terminus of the Sawtooth fault (Liberty et al., 2020; Wesnousky, 2021), with abundant aftershocks largely north of the Sawtooth fault (Fig. 1) and along possibly reactivated relict crustal structures (Liberty et al., 2020). The low double‐couple (56%) moment tensor for the main shock (see Data and Resources) and ∼50 km long diffuse zone of aftershock seismicity highlight likely multifault rupture at depth (Yang et al., 2021), which did not propagate to the surface (Liberty et al., 2020). This enigmatic earthquake occurred within a zone of unmapped faults with no apparent surface expression, raising questions about the origin and kinematic history of faults in the intersection zone between the western CTB and TCFS (Liberty et al., 2020; Yang et al., 2021; Luo et al., 2022).

The CTB–TCFS region and 2020 Stanley earthquake offer insights into fault kinematics and seismicity in a region of distributed and intersecting faults, similar to those found in the Basin and Range Province and Walker Lane (Wesnousky, 2005). Here, we seek to better understand the rupture history of the Sawtooth fault and kinematics of the western CTB by integrating results from the Dutch Lake site—the first paleoseismic trench investigation of the Sawtooth fault. We evaluate Sawtooth fault rupture scenarios and complex faulting in the epicentral region of the 2020 earthquake, and explore the role of the relict TCFS in controlling rupture length while accommodating displacement transfer. By unraveling these interconnected fault systems, we improve understanding of recent rupture behavior and current seismic hazard in the CTB.

At the Dutch Lake site on the northern Sawtooth fault (Figs. 1, 2; 44.34381, −115.17643), we excavated a trench across a ∼2 m high, northeast‐facing fault scarp, which traverses a postglacial (≤14–16 ka) alluvial‐fan surface inset between glacial moraines (Fig. 2). The Dutch Lake scarp represents the northern extent of unambiguous faulting along the Sawtooth fault (Lifton et al., 2023). To the north, scarps are less well defined or of nontectonic origin (e.g., sackungen). To the south, a ∼2 km long gap separates the Dutch Lake fault scarp from the remainder of the Sawtooth fault, which continues southeast for ∼60 km (Thackray et al., 2013).

The Dutch Lake trench (Fig. 3) revealed three main sedimentary and pedogenic units: (1) alluvial‐fan gravel, (2) a paleosol within the upper part of unit 1, and (3) scarp‐derived colluvium, with a modern soil A horizon at the surface (unit 3A). Gray sand and gravel in unit 1 is >1 m thick, clast supported, locally imbricated, and includes indistinct to well‐defined bedding and pedogenic carbonate. We interpret unit 1 as alluvial‐fan sediments transported in a postglacial, sediment‐rich environment. Unit 2 includes gray–brown silt, sand, and gravel and is similar to unit 1, but with a greater concentration of fines that yield a more massive and matrix‐supported texture. Although soil structure in unit 2 is weak, we infer a pedogenic origin based on the unit’s fine‐grained and cohesive texture, lateral continuity, and constant 0.2–0.3 m thickness. Colluvium in unit 3 consists of tan‐brown silt, sand, and gravel with laterally variable texture and discontinuous slope‐parallel bedding. We interpret a scarp‐derived origin for the colluvium based on its wedge geometry, which tapers from 0.8–1.0 m beneath the scarp to <0.2 m at the trench ends. Unit 3 includes macrofloral remains, is locally bioturbated, and is overprinted by a ≤0.2 m thick dark brown soil A horizon (unit 3A) near the surface.

Units 1 and 2 are faulted by vertical to steeply west‐dipping zones of shearing. We identified three primary fault zones in the trench, which together span ∼3–4 m of the ∼8 m wide scarp face (Fig. 3). The vertical displacement of the top of unit 2 is ≥1.0–1.8 m in the south wall and ≥1.5–1.9 m in the north wall. These displacement observations are minima due to the limited exposure of the postglacial fan units in the fault footwall. Units 1 and 2 are unconformably overlain by unit 3, which is the thickest in the hanging wall of the easternmost shear zone (south wall; Fig. 3). We interpret units 1 and 2 as correlative across the fault zone based on lateral consistency in their texture and color. Although the vertical to steep west dip and distributed nature of faulting supports a lateral component of slip, the contribution is likely minor owing to the lack of laterally offset features apparent in the site geomorphology.

Radiocarbon (⁠$C14$⁠) and feldspar infrared‐stimulated luminescence (IRSL) ages provide temporal constraints for the sedimentary deposits (DuRoss et al., 2023; Fig. 3). Thirteen $C14$ ages are for Pinus charcoal sampled ∼0.3–0.5 m below the surface within parts of unit 3 proximal to the fault zones. These calendar‐calibrated ages constrain unit 3 to ∼2.4–4.2 ka. An additional age (R13, ∼0.2 ka) for a single 0.3 mg twig from unit 1 is unreliable (inverted with all ages for units 2 and 3) and not used in subsequent age analyses. Eight IRSL samples include sediment from units 1–3. We apply the minimum age model (MAM) to the multialiquot luminescence data considering the coarse and locally poorly sorted nature of the sediments (i.e., short transport distance and poor bleaching potential). MAM ages for the uppermost part of unit 1 (∼18.6 ka) and unit 2 (∼11.7–14.6 ka) are most consistent with the postglacial (≤14–16 ka) interpretation of the alluvial‐fan setting. We exclude two IRSL samples (L1 and L2) from further interpretation, as their resulting 36–69 ka ages are substantially older than overlying ages for units 1 and 2 (11.7–18.6 ka), likely as a result of partial bleaching during sediment transport. IRSL ages for L4, L5, and L7 constrain unit 3 to ∼6.6–8.6 ka.

We interpret a single postglacial surface‐faulting earthquake at the Dutch Lake site. This interpretation is based on the (1) vertical displacement of units 1–2; (2) multiple fault terminations at the units 2–3 contact; (3) unconformable relation between units 1–2 and unfaulted unit 3; (4) lack of paleosols within unit 3; (5) agreement between the estimates of vertical displacement (≥1.0–1.9 m), the cumulative thickness of unit 3 across the shear zones (∼2.0–2.3 m), and the vertical separation across the fault scarp (2.0–2.2 m; Fig. 2); and (6) lateral consistency in the height of the postglacial fault scarp at the site and on glacial surfaces to the north and south (Fig. 2).

Four Bayesian time‐stratigraphic models (Fig. 4) constrain the timing of Sawtooth fault surface rupture at the Dutch Lake site. With the exception of R13, all $C14$ ages and IRSL sample L6 (∼18.6 ka) are used in all models. Model variations 1–4 (Fig. 4a,b) allow for different combinations of limiting IRSL ages. Model 1, which includes all IRSL ages, constrains the Dutch Lake earthquake rupture to 10.0 ± 1.0 ka (⁠$1σ$⁠), with L3 (∼11.7 ka; unit 2) and L4 (∼8.6 ka; unit 3) providing the closest temporal limits. Model 2, which relies on L8 (∼14.6 ka) as the maximum constraint and $C14$ ages (≤4.2 ka) as the minimum, yields a broad earthquake time range of 9.5 ± 2.9 ka (⁠$1σ$⁠). In model 3, the event is constrained to 9.1 ± 1.5 ka (⁠$1σ$⁠) using L3 as the maximum and L7 (∼7.2 ka) as the minimum; this model assumes that the youngest (L7) of the three IRSL ages for unit 3 is the most accurate. Model 4 is similar to model 2, but relies on L3 as the closest maximum age, and generates a rupture time of 8.0 ± 2.1 ka (⁠$1σ$⁠).

We equally weight and sum the four earthquake‐timing probability density functions (PDFs) generated in the Bayesian models, after DuRoss et al. (2018). The resulting PDF for the Dutch Lake earthquake (herein, DL1) is asymmetric, with a negative skew and a mean age of 9.1 ± 2.1 ka (⁠$1σ$⁠) (Fig. 4b). The summed PDF indicates that DL1 occurred between 4.8 and 13.6 ka at 95% confidence, with the most likely (modal) time at ∼9–11 ka. The occurrence of DL1 from ∼4.4–6.2 ka to ∼12.5–14.5 ka is possible, but less likely. We consider it highly unlikely that DL1 occurred <4 ka or >16 ka, that is, outside of the limits of the DL1 PDF distribution. An early Holocene time range for DL1 is also consistent with the weak paleosol (unit 2), which likely formed and was subsequently faulted and buried <5 ky following deglaciation.

Lacustrine paleoseismic data from glacial–moraine impounded lakes along the central part of the Sawtooth fault facilitate a comparison of earthquake timing along the fault (Fig. 5a). Shapley et al. (2023) used disturbance horizons within sedimentary cores from Redfish, Little Redfish, and Pettit Lakes to interpret two to three Holocene ruptures of the fault. The youngest earthquake near Redfish Lake (RL1) occurred at 4.1–4.6 ka (Redfish Lake) to 4.0–4.4 ka (Little Redfish Lake) based on $C14$ ages that bracket a homogenite to mud–clast conglomerate. In Pettit Lake, 16 km to the south, a ∼3 cm thick normally graded bed provides evidence for an earthquake at 4.8–5.3 ka (PL1). Shapley et al. (2023) concluded that PL1 and RL1 likely represent separate shaking events based on the high‐fidelity lacustrine strata and the minimal overlap (0.2 ky at 95% confidence) in the event PDFs. Disturbed Redfish Lake sediments provide the basis for an older earthquake (RL2), which occurred between deglaciation (<14 ka) and deposition of Mazama tephra at ∼7.6 ka (Shapley et al., 2023).

Earthquake DL1 overlaps temporally with Redfish Lake events RL1 at 4.0–4.6 ka and RL2 at 7.6–14.0 ka, as well as Pettit Lake event PL1 at 4.8–5.3 ka. However, the most substantial area of overlap is between the modal time for DL1 at 9–11 ka and RL2 at 7.6–14.0 ka. Although DL1 could have occurred near the extreme young end of its PDF distribution (∼4 ka), we consider this scenario unlikely as hundreds of years likely elapsed between surface rupture, scarp degradation, colluvial wedge deposition, and eventual incorporation of macrofloral remains (≤4.2 ka) into the colluvial deposits. Further, the weak paleosol (unit 2) beneath the colluvium indicates an early Holocene DL1 time. Thus, we consider it most likely that DL1 is correlative with the RL2 (7.6–14.0 ka) rather than the RL1 or PL1 (Fig. 5a).

Here, we evaluate the kinematics and rupture of the Sawtooth fault in the context of the 2020 Stanley earthquake epicentral region.

Terrestrial and lacustrine paleoseismic data permit rupture of the majority of the mapped Sawtooth fault in the early Holocene (Fig. 5a). This rupture could have involved discontinuous scarps along two subparallel and ∼13 km long fault strands that form the northern end of the Sawtooth fault (Lifton et al., 2023; Fig. 2). The southern extent of this rupture remains unconstrained. Scenarios including rupture of the entire Sawtooth fault yield 64 km (end‐to‐end) or 69 km (curving trace) of rupture length and an $Mw$ 7.2–7.4 earthquake (using $Mw$–area scaling relations; Shaw, 2023).

Lacustrine paleoseismic data and light detection and ranging (lidar) fault mapping (Lifton et al., 2023) support partial ruptures of the Sawtooth fault. For example, separate rupture of the southern (PL1) and central–northern (RL1) parts of the fault agrees well with an along‐strike change in geomorphic expression between these lakes (Fig. 5a). The central part of the fault consists of a 27 km long, east‐convex zone of mostly continuous fault scarps. In contrast, the southern ∼30 km of the fault is west‐convex with highly discontinuous scarps (≤7 km apart). The diminished geomorphic expression to the south could be the result of more degraded scarps, smaller displacement, and/or more distributed faulting at the surface. However, coarse glacial debris, rough surfaces, and incomplete lidar coverage could also help explain the variable geomorphic expression. Thackray et al. (2013) and Shapley et al. (2023) used this change in scarp expression and differences in the RL1 and PL1 event times to infer a point of rupture termination between the central and southern portions of the fault (Fig. 5a). Ruptures of the central or southern parts of the fault would generate $Mw$ 6.8–7.0 earthquakes (above the $Mw$ ∼6.5 threshold for surface rupture; Wells and Coppersmith, 1994), whereas combinations of the central, northern, or southern portions could range from 38 to 56 km and $Mw$ 7.0–7.2.

The Stanley earthquake occurred in an area of distributed Quaternary faulting (Lifton et al., 2023), with aftershocks aligned both north–south and east–west (Fig. 5b). In a model that integrates aftershock seismicity, teleseismic records, and interferometric synthetic aperture radar, Yang et al. (2021) showed that subsurface slip in the Stanley earthquake occurred on two north–south striking strike‐slip faults (F1 and F2, Fig. 5b) that have opposing dips and bound the northernmost Sawtooth fault (e.g., Wilbur, 2022). These faults lack previously recognized surficial expression (Liberty et al., 2020). The northern fault, F1, experienced peak slip of ∼0.7 m at 11–12 km depth (Yang et al., 2021), likely dips 72°–74° west (Luo et al., 2022), and agrees with the kinematic model of Pollitz et al. (2020). F2 experienced about ∼0.2 m of slip (Yang et al., 2021) and is possibly vertical to ∼60° east‐dipping (Yang et al., 2021; Luo et al. 2022). Yang et al. (2021) noted that compression in the F1–F2 overlap area could generate orthogonal reverse faulting, subparallel to the TCFS.

Lidar mapping in the epicentral region of the 2020 earthquake shows a broad zone of disconnected faults both within and north of the TCFS (Fig. 5b,c), including the 8 km long Cape Horn fault (Lifton et al., 2023). These surface traces spatially agree with faults F1 and F2 of Yang et al. (2021); however, inferred fault dips for the Cape Horn fault (west) and F2 (east) disagree. Although part of F2 is in the footwall of the Sawtooth fault, no obvious kinematic link exists between these faults. In general, the discontinuous surface geometry of faults in this region argues for spatially limited rupture. Considering the lack of connecting (e.g., transfer) structures at the surface between the epicentral region faults and the Sawtooth fault, we infer that the Sawtooth fault and the 2020 sequence and causative faults are not directly linked.

Ruptures in the complex intersection of the northernmost Sawtooth fault, TCFS, and Cape Horn fault are likely limited spatially (and in $Mw$⁠) by the distributed, disconnected, and possibly incipient nature of the faults accommodating regional extension (Fig. 5c). Ruptures in this region, especially within or north of the TCFS, may have lengths ≤25–35 km (based on fault mapping or Stanley earthquake models), steep dips (≥60°), and $Mw≤6.9$⁠. However, the cryptic nature of faulting in this region limits our ability to infer fault connectivity in this system.

We hypothesize that the northeast‐trending TCFS, which includes a 7–12 km wide multistranded zone of extension without evidence of Quaternary slip, acts as a barrier to ruptures along and between the Sawtooth fault and the distributed fault system to the north (Fig. 5c). The TCFS may have acted as a termination point for the early Holocene and the previous ruptures of the Sawtooth fault, as no clear geologic evidence of Holocene faulting is present north of the Dutch Lake fault scarp. The TCFS also appears to exert control on the lateral extent of faults such as the Cape Horn fault (Fig. 5c) in the 2020 epicentral region. Although 2020 seismicity spans the TCFS, most seismic moment release occurred within this structure (e.g., 76% on fault F1; Yang et al., 2021) and distributed aftershocks illuminate complex faulting possibly influenced by inherited structures (Wilbur, 2022). Considering the potential for intersecting faults at depth (Fig. 5c), the TCFS has likely acted as a barrier to rupture propagation (e.g., King and Nábělek, 1985).

The TCFS as a barrier to rupture is supported in the framework of regional crustal extension and rotation (Fig. 5d). Brocher et al. (2017) explain CTB extension in the context of clockwise crustal motions about a geologic pole in eastern Washington. Extension is accommodated by normal faults radial to the pole, such as the Sawtooth, Lost River, Lemhi, and Beaverhead faults. In contrast, zones of strike slip and distributed displacement along small circles around the pole may serve to transfer slip between radial fault systems (Brocher et al., 2017). We infer that the Eocene TCFS, in the 2020 epicentral region, also serves as a conduit for the transfer of slip between the Sawtooth fault and multifault networks to the southwest (e.g., the normal oblique Deadwood fault; Fig. 5d) or northeast along the northern margin of the CTB, enabling displacement transfer to faults such as the Bitterroot fault in western Montana. In general, the TCFS is misoriented by ∼20° with small circle paths (Fig. 5d), and, thus, we assume that displacement transfer occurs either as moderate seismicity such as the Stanley or 2014–2017 Challis (northern Lost River fault zone; Pang et al., 2018) sequences (Fig. 1) or oblique slip along distributed, cryptic, and spatially limited faults rather than along a continuous, reactivated shear zone. In the northern CTB, increasing the density of the regional seismic network and spatial coverage of high‐resolution topographic data would improve understanding of the mechanics and hazard implications of low rate, distributed slip along the TCFS.

Our results emphasize that integrating geological, geodetic, and seismic data can facilitate the exploration of complex faulting at the intersection of differing fault systems. Such an approach could expand our understanding of the kinematics and hazard of complex systems globally that exhibit slip transfer through complicated networks (e.g., central Walker Lane, Wesnousky, 2005; central Apennines, Galli et al., 2010), host distributed surface rupture (e.g., Koehler et al., 2021), and contrast with displacement transfer along more organized systems (e.g., Garlock fault; Davis and Burchfiel, 1973).

We draw the following conclusions from our comparison of paleoseismic and modern seismicity in western CTB:

1. Dutch Lake paleoseismic data confirm early Holocene rupture of the northern Sawtooth fault at 9.1 ± 2.1 ka. Terrestrial and lacustrine data support varying Sawtooth fault rupture lengths and magnitudes (⁠$Mw$ 6.8–7.4).

2. Complex faulting in the 2020 Stanley epicentral region includes distributed, discontinuous, and opposing‐dip faults. Oblique extension in this distributed system is likely to generate earthquakes of moderate length (≤35 km) and magnitude (⁠$Mw≤6.9$⁠).

3. The TCFS likely plays a role in controlling fault length and rupture continuity for adjacent active faults while also facilitating displacement transfer between normal faults that accommodate regional extension and clockwise crustal rotation.

4. Relict structures have the potential to influence fault growth, rupture extent, and magnitude and may help inform fault characterizations for PSHA.

The 2020 Stanley mainshock is available at https://earthquake.usgs.gov/earthquakes/eventpage/us70008jr5/executive; aftershocks are available at https://earthquake.usgs.gov/earthquakes/map/. Quaternary faults are available at doi: 10.5066/P9BCVRCK; Trans‐Challis fault system (TCFS) mapping: Lewis et al. (2012). National Elevation Dataset (NED) and 0.5 m light detection and ranging (lidar) data are available at https://apps.nationalmap.gov/downloader/. All the sites were last accessed in October 2023. Trench photomosaics (Fig. S1), code for Bayesian models (Code S1), and $Mw$ calculations (Table S1) are included in the supplemental material.

The authors acknowledge that there are no conflicts of interest recorded.

This work was supported by the U.S. Geological Survey (USGS) Earthquake Hazards Program and Geologic Hazards Science Center. The authors thank Ellyson Long and the U.S. Department of Agriculture (USDA) Forest Service for assistance with trench permitting. Thanks to Dan O’Connell and Bob Creed for conducting a nodal geophysical survey prior to trenching, and Claudio Berti and Lee Liberty for discussions of the Dutch Lake trench. Charles Trexler, Will Yeck, two anonymous reviewers, and Editor‐in‐Chief Keith Koper provided constructive comments that improved this article.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.