Slip and Strain Accumulation Along the Sadie Creek Fault, Olympic Peninsula, Washington

Upper‐plate faulting in the Olympic Peninsula of Washington State results from relative motion of crustal blocks within the Cascadia forearc and earthquake cycle processes along the Cascadia subduction zone. We reconstruct fault slip rates since ∼14 ka on the Sadie Creek fault (SCF), north of the Olympic Mountains, using airborne lidar and field mapping of surficial deposits and landforms and optically stimulated luminescence and radiocarbon dating. The SCF is a ≥14 km‐long northwest‐striking, subvertical, dextral strike‐slip fault with a subordinate dip‐slip component. Laterally, offset debris flow channels cut into Late‐Pleistocene and younger surfaces show dextral slip of 4.0–24.5 m and dip slip of 0.7–6.5 m. Re‐evaluation of fault slip on the adjacent Lake Creek Boundary Creek fault (LCBCF) shows similar dextral (4.5–29.7 m) and dip slip (0.8–4.6 m). A deglacial age of 14 ka paired with the largest—and presumably oldest—slip measurements produce a minimum dextral slip rate of 1.3–2.3 mm/yr and dip‐slip rate of 0.05–0.5 mm/yr. Similarities in geometry, kinematics, slip rate, and earthquake timing between the SCF and LCBCF suggest these faults represent one continuous geologic structure, the North Olympic fault zone. Geodetically constrained boundary element method models considering the effects of coseismic subduction zone stresses on upper plate structures produce comparable kinematics to those measured on the SCF and LCBCF, suggesting that the North Olympic fault zone acts as the main strain‐accommodating structure in the northern Olympic Peninsula and may be modulated by stress transferred from subduction zone earthquakes.

geologic context for the tectonic development of the NOFZ, and lend support to an overall kinematic and seismologic connection between the SCF and LCBCF as a continuous upper plate fault (the NOFZ). Finally, we compare our slip rate data with the predictions of a boundary element method model that estimates the slip on the NOFZ (Figure 2a) necessary to relieve stresses imposed by various earthquake cycle processes on the underlying Cascadia subduction zone. This comparison sheds light on potential interactions between the megathrust earthquake cycle and slip on upper plate faults.

Geologic and Tectonic Setting
The Olympic Mountains of western Washington (Figure 1b) represent the structural and topographic high of the Cascadia forearc, and the only location where the Cascade accretionary wedge has been exposed subaerially (Brandon & Vance, 1992;Brandon et al., 1998). The accretionary complex consists of marine sedimentary rocks underthrust beneath Eocene marine to nonmarine basalts of the Crescent Formation and Eocene-Miocene marginal marine sedimentary rocks (Brandon & Vance, 1992;Brandon et al., 1998;Tabor & Cady, 1978a). Within the northern Olympic Peninsula, Eocene-Miocene marine sedimentary rocks dip moderately north to northeast and are cut by northeast-dipping reverse faults (Tabor & Cady, 1978a).  show similar geometries to bedrock faults in the Olympic Peninsula, possibly reactivating these older structures (Barnett et al., 2015;Nelson et al., 2017). Active faults across this region accommodate north-south shortening through some combination of reverse slip and conjugate strike slip faulting (Figure 1b). Generally, slip on northwest-striking faults is dextral-oblique (e.g., the southern Whidbey Island fault; Johnson et al., 1996;Sherrod et al., 2008), on northeast-striking faults is sinistral-oblique (e.g., the Saddle Mountain fault; Barnett et al., 2015;Witter et al., 2008), and on west-striking faults is primarily reverse (e.g., the Seattle fault; Blakely et al., 2002;Johnson et al., 1999; Figure 1b).

Late Pleistocene Glaciation
The Juan de Fuca lobe of the Cordilleran ice sheet, which split and flowed westward from the Puget Lobe at ∼48° latitude (Figure 1), reached elevations of up to 1 km in the northern Olympic Peninsula (Polenz et al., 2004;Porter & Swanson, 1998;Thorson, 1980). As such, glacial landforms and deposits from the Juan de Fuca lobe define the morphology and surficial geology of the region (Mosher & Hewitt, 2004;Tabor & Cady, 1978a). Alpine glaciers were extensive within the core of the Olympic Mountains during the last glacial maximum, however, deposits around the NOFZ were solely derived from the Juan de Fuca lobe, except for possible input of alpine-derived glacial outwash into the Elwha River drainage (Figure 2a; Polenz et al., 2004;Schasse, 2003;Tabor & Cady, 1978a). Existing radiocarbon ages suggest that the Juan de Fuca lobe of the Cordilleran ice sheet advanced to its maximum extent after ca.16 ka and had retreated by DUCKWORTH ET AL.   (Nelson et al., 2017) and unstudied scarps to the west of the West Twin River. White lines show the extent of lidar coverage and are labeled by the year collected. The lidar released in 2019 was not made available in time for detailed analysis in this study. (b) Generalized surficial geologic map from Duckworth (2019) showing the locations of samples (triangles) and lateral offset measurements (circles). Yellow and purple markers/text show the location, names, and ages of OSL and radiocarbon samples and black and blue markers/text show the location, names, and preferred dextral offsets of the offset features. All locations where separate offset measurements were derived from the channel and the interfluve are labeled with a subscript "c" and "i," respectively. Sample names which begin with "L" represent luminescence samples. The extents of Figures 3 and 4a are shown by dashed white boxes.
ca.14 ka. The timing of Juan de Fuca lobe advance is supported by radiocarbon dating of wood samples (15,200 calibrated years before present [Cal yrs. BP]) from Lake Ozette, WA (Figure 1b) originally published by Heusser (1973) and later reinterpreted as pre-glacial wood in lodgment till by Haugerud and Hendy (2016). Radiocarbon dating of post-glacial wood (14,095 Cal yrs. BP) from near Bellingham, WA (Kovanen & Easterbrook, 2001) and bone fragments (13,010-13,380 Cal yrs. BP) from the Manis Mastodon site in Sequim, WA (Waters et al., 2011) constrains the timing of Juan de Fuca lobe retreat (Figure 1b). In this study, we contribute new radiocarbon and luminescence ages that bear on the timing of glacial retreat along the SCF and northern Olympic Peninsula.

Active Faults in the Northern Olympic Mountains
The SCF, originally identified through analysis of gravity data (MacLeod et al., 1977) and later named by Joyner (2016), was first mapped by Nelson et al. (2017) in a paleoseismic study of the LCBCF after the release of a high resolution lidar data set in 2015 (Puget Sound Lidar Consortium, 2015). The LCBCF, SCF, and newly discovered scarps west of the West Twin River (U.S. Geological Survey, 2018; Figure 2a) are termed the Northern Olympic fault zone by Schermer et al. (2020). Paleoseismic trenching (Nelson et al., 2017;Schermer et al, 2020;Figure 2a) and coring and seismic imaging of megaturbidites within Lake Crescent (Leithold et al., 2019) document three to five surface-rupturing earthquakes since retreat of the Juan de Fuca lobe. Schermer et al. (2020) correlated three to four of these earthquakes across most of the length of the NOFZ, suggesting that the NOFZ has at times ruptured in M7.0-7.5 earthquakes. Nelson et al. (2017) measured dextral slip and vertical separation on the LCBC and SCF as ∼11-28 m and 1-2 m respectively and, using deglaciation at ∼13 ka as a maximum age, calculated a slip rate of ∼1-2 mm/ yr. Nelson et al. (2017) included some analysis of channels offset by the SCF, although they did not perform detailed geologic mapping or dating of landforms. The data presented here builds on their analysis by providing the necessary geologic context for understanding the incremental and cumulative displacement along the SCF and by modeling the dynamic processes that drive deformation on the entire NOFZ.

Fault Morphology and Surficial Deposits
New surficial geomorphic mapping based on both fieldwork and lidar interpretation allows for interpretation of the history of the geomorphic and tectonic processes that have sculpted the landscape around the SCF (Figures 2b and 3). Glacially scoured bedrock overlain by till (Qgt), ice contact deposits (Qgoi), moraines, and outwash terraces (Qgo) reflect the imprint of the continental ice sheets that flowed westward across the area during the last glacial maximum. Glacial units are reworked by postglacial processes to form colluvial slopes (Qc), alluvial fans (Qf), landslide deposits (Qls), and alluvial gravel/terrace deposits (Qal/ Qt) (Figures 2b and 3).
We used lidar-derived hillshade and slope images computed using ArcGIS software, 3-D visualization in Quick Terrain Modeler, and scarp-perpendicular topographic profiles to map the surface trace of the SCF from Lake Crescent westward to the northwestern extent of the 2015 lidar (Puget Sound Lidar Consortium, 2015;Figure 2a). The surface trace of the SCF is roughly linear, cuts glacial and post-glacial deposits and landforms (e.g., Figures 2b, 3, and 4), and is laterally continuous for ≥14 km westward from Lake Crescent ( Figure 2). Lidar data released in late 2019 (U.S. Geological Survey, 2018) reveals well defined fault scarps that continue northwestward for another >25 km that were not analyzed in this study (Figure 1b;Schermer et al, 2020). The SCF scarp switches from north-facing east of the Lyre River to south-facing west of the Lyre River and predominantly strikes northwest-southeast, except between Sadie Creek and the East Twin River where it strikes east-west ( Figure 2).
The SCF dextrally offsets 11 debris flow channels (e.g., Figure 4 and Table 1, supporting information S1) between Lake Crescent and Susie Creek. We did not identify any laterally offset features west of Susie Creek where the SCF is defined by scarps at lower elevations and cuts the valley floor (e.g., Figure 3) because there are no linear features capable of preserving lateral displacement in that region. All laterally offset debris flow channels are incised into till deposits on steep hillslopes aside from one channel ("J Channel") which is incised into outwash (Table 1, supporting information S1). Changes in scarp facing direction along with a relatively linear overall fault trace and dextral offset of debris flow channels support the interpretation of DUCKWORTH ET AL.   Tabor and Cady (1978b). Samples locations and ages are shown by yellow (OSL samples) and purple (radiocarbon and OSL) triangles, with the last three digits of their sample name. Nelson et al. (2017) that the SCF is a steeply dipping, dominantly right lateral strike-slip fault with a subordinate dip slip component.

Fault Slip Analysis
We combined field and lidar measurements to map and analyze linear features displaced by the SCF. This study uses the channel analysis program, "3D_Fault_Offsets," of Stewart et al. (2017) to measure the strikeslip component of the total slip vector (herein referred to as "strike slip" or "strike-slip displacement") from a lidar-derived digital elevation model (DEM). We observe that in these laterally offset channels, vertical scarps are rarely preserved, suggesting that either the dip slip component is either removed or modified by post-earthquake debris flow erosion or deposition. As such, we measure the total dip-slip component of displacement (herein referred to as "dip slip" or "dip-slip displacement") on the intervening geomorphic surfaces or interfluves using fault-perpendicular topographic profiles and the methods of Thompson et al. (2002). Importantly, our dip slip measurements likely represent cumulative offset since formation of the deglacial or post-glacial surfaces, while our strike slip measurements are more likely to capture incremental offset from earthquakes intervening with channel formation.

Strike-Slip Displacement
We measured 11 laterally offset debris flow channels (e.g., Figure 4) on the SCF and two channels on the LCBCF not previously measured by Nelson et al. (2017) to calculate the amount of strike-slip displacement (Table 1; supporting information S4). In each case, points along morphologic markers of the channel (e.g., channel bottom, riser midpoint, etc.) were identified within a lidar DEM based on relative elevations, steepness, and changes in concavity using 3D_Fault_Offsets. Markers were then regressed and projected onto the fault trace in order to measure the strike-slip displacement (Stewart et al., 2017; Figure 4d). For each channel, we visually assess the quality of marker identification (Table 1), removing measurements where automatically identified points do not appear to accurately follow their respective geomorphic feature, and interpret the distribution of measurements in terms of the geomorphic development of the channel (supporting information S3 and S7). For example, a lower magnitude of strike-slip displacement measured from the channel bottom than the midpoint of the riser may suggest that channel erosion has minimized the record of strike-slip displacement (e.g., Figure S3a; Reitman et al., 2019). Thus, the measurement made from the channel riser midpoint represents a more accurate representation of strike-slip displacement. After all measurements are assessed, each slip measurement is assigned a probability distribution based on the spread of the individual marker measurements and the interpretation of channel development. We assign gaussian distributions if we believe the lateral slip record is well preserved (e.g., multiple markers show overlapping measurements). However, if the record of lateral slip does not appear well preserved and there is no geomorphic reason to favor the measurement from one marker over another, we assign a uniform distribution to equally account for the uncertainty across all markers. Finally, we evaluate the measurement by back-slipping the offset channel (Haddon et al., 2016;Stewart et al., 2017;Zielke, 2009) along the fault trace to the preferred value, lower bound, and upper bound of the reported offset. If the backslip does not visually show a satisfactory pre-fault reconstruction, we adjust the probability distributions and the resultant errors (see supporting information S3 for further discussion of lateral slip probability distributions). Comparison with our field measurements suggests overall agreement with our lidar analysis ( Figure S3b).
Strike-slip measurements from offset features on the LCBCF presented by Nelson et al. (2017) were verified through backslipping and adjusted at several locations (Table 1). For each of these measurements we assigned triangular probability distributions to reflect the maximum, minimum, and mostly likely value of the range of each reported offset. Additionally, we measured strike-slip displacement at three locations ("NOL03 channel," "LCBC-CD1," and "LCBC-CD2," Table 1) on the LCBCF not reported by Nelson et al. (2017) using the same channel analysis methods outlined above.
Strike-slip displacement for both the SCF and LCBCF ranges from 4.0 to 29.7 m ( Figure 5, Table 1, Figure  S3d). The largest strike-slip displacements on the SCF occur on the eastern section, between the Lyre River and Lake Crescent ( Figure 5 between ∼12 and 15 km) and at three locations on the LCBCF ( Figure 5). The smallest strike-slip displacement was measured at "J Channel" (Table 1) a location where a scarp is preserved within the channel and the channel is offset 4.0 ± 0.8 m, significantly less than the interfluve ("J Interfluve") which is offset 8.3 ± 1.8 m (Schermer et al., 2020). Additionally, channels that preserve varying magnitudes of strike-slip displacement are often located in close proximity (e.g., Channels B, C, and D; Figure 4a; Table 1) suggesting a range of channel ages that record different numbers of surface-rupturing earthquakes.

Dip-Slip Displacement
We analyzed 48 scarp profiles on the SCF to calculate the amount of dip-slip displacement of post-glacial deposits and landforms. Scarp profiles were extracted from the ground returns of the lidar point cloud, in three-meter-wide swaths orthogonal to the scarp face (e.g., Figure 4a). In cases where geomorphic piercing lines, such as the top edge of a channel riser were traceable across the fault, we extracted and combined profiles from multiple line segments (e.g., Figure 4, profiles NOL8E and NOL8W). We assume generally steep fault dips ranging from 60° to 90° and follow the procedures and parameters of Thompson et al. (2002) to calculate the amount of vertical separation, the component of dip slip on the fault plane, and the associated errors for each profile from a Monte Carlo simulation. Given that the resulting dip slip and vertical separation distributions from each Monte Carlo simulation resemble normal distributions, we report the preferred dip slip as the median of the distribution produced by Monte Carlo simulation and the uncertainty for each measurement as the 95% confidence interval around the median (e.g., Zechar & Frankel, 2009). Because Nelson et al. (2017) reported vertical separation (without errors) as opposed to dip slip, we calculated dip slip at their scarp profile locations using the same scarp profiling methods used in this study; supporting information S5).
Dip-slip displacement ranges from 0.5 to 6.5 m along the SCF and 0.7 to 4.6 m along the LCBCF ( Figure 5, Tables S5a and S5b). Dip-slip displacement does not appear to have any distinct groupings or spatial correlations along the SCF with the exception of three measurements on the west-striking portion of the SCF that have relatively large dip-slip displacement (∼4.1-6.5 m; Figure 5, between 3 km and 5 km), and measurements of faulted postglacial deposits that have relatively less dip-slip displacement (1.0-3.7 m; Figure 5, blue triangles between 5 km and 8 km).
The consistency in the geometry and range of post-glacial strike slip and dip slip between the SCF (4.0-24.5 m strike slip; 0.5-6.5 m dip slip) and LCBCF (4.5-29.7 m strike slip; 0.7-4.6 m dip slip) lends support to the inference that SCF and LCBCF should be considered as one connected fault system (the NOFZ; Schermer et al., 2020). As such, we consider the measurements made on the SCF and LCBCF together for the calculation of a single post-glacial slip rate (Section 6) and as a single fault system in models of forearc deformation (Section 7).

Total Fault Slip and Kinematics
Characterization of total fault slip and kinematics requires measurements of dip-slip and strike-slip displacement which, ideally, are made on the same feature. However, the only record of strike-slip displacement preserved on the SCF comes from debris flow channels that do not consistently preserve a scarp within the channel. As such, we calculate total slip, strike slip to dip slip ratios, and rake using strike-slip measurements from dextrally offset debris flow channels and dip-slip measurements from scarp profiles in the immediate vicinity (<15 m) of the offset channel (Table S5a) shows the correlation of scarp profiles with offset channels). In order to combine dip-slip and strike-slip measurements in this manner, we assume that: (1) the initial offset of the channel and formation of the fault scarp are the same age; (2) laterally offset markers within the channel are protected from geomorphic modification during the interseismic period; and (3) scarp degradation outside of the channels is minimal during the interseismic period. Given that one of these assumptions is likely incorrect for some channels, we assert that the measured total slip is less than the true total cumulative fault slip. Additionally, the strike slip to dip slip ratio will most likely appear DUCKWORTH ET AL.
10.1029/2020JB020276 9 of 22 less than the true strike slip to dip slip ratio. Interseismic geomorphic modification is likely to be greatest in the channels, thus minimizing the amount of strike slip recorded by channel markers relative to dip slip preserved by fault scarps located on the interfluve (e.g., Figure S3a; Reitman et al., 2019). Calculated dip-slip is likely a minimum owing to scarp degradation, but this process is likely less efficient than erosion within the channels. Taken together, all calculated strike slip to dip slip ratios, total slip, and rakes are minima (rake is skewed towards ± 90°).
We fit kernel density estimates to the distributions of dip-slip displacements of all profiles considered at each offset channel on the SCF and LCBCF to capture the variability of dip-slip measurements within their associated errors (Table 1). In some instances, we omitted a profile from the fit of density functions if there was substantial geologic evidence that the dip-slip displacement measurement was of low quality. For example, we omitted profiles if there was uncertainty in the type of deposit on either side of the fault trace or evidence of geomorphic modification at the scarp or such as SP02 (Figure 4c). We then conducted a Monte Carlo simulation using the kernel density distributions of the dip-slip measurements and the assigned probability distributions of the strike-slip measurements to calculate the strike slip to dip slip ratio, rake, and the total slip at each slip site. Given that the resulting distribution often are skewed, we report strike slip to dip slip ratios, rakes, and total slip amounts and their uncertainties as the peak value (mode) of the kernel density fit to the result of the Monte Carlo simulation and the 95% confidence interval around the median, respectively (Thompson et al., 2002;Zechar & Frankel, 2009; Table 1).
Our calculations show that the rake of the inferred slip vector is generally >150° (>2:1 strike slip to dip slip ratio) along strike of the SCF and LCBCF but most measurements show rakes of ∼170° (∼6:1 strike slip to dip slip ratio) (Table 1, Figure S3c). We also note a slightly more oblique rake (lower strike slip to dip slip ratio) for the SCF (164° average rake) than the LCBCF (172° average rake), due to the greater maximum dipslip displacement on the SCF ( Figure 5, Table 1). We infer that the largest lateral offsets are recorded by the oldest channels (see Section 5), and thus the rake of the slip vector from these channels is likely to be most representative of the cumulative slip recorded by the interfluve scarps. The preferred rakes for channels that exhibit >15 m lateral slip ranges from 162 to 175 on the SCF and 172 to 179 on the LCBCF (Table 1), consistent with the range of rakes recorded by all channels with their nearby scarp profiles, and consistent with the interpretation that the lateral to vertical ratio of the SCF is smaller than that on the LCBCF.

OSL and Radiocarbon Dating
We collected samples from various glacial (till, outwash, ice contact, etc.) and post-glacial (terraces, debris flows, alluvial fans, etc.) deposits along the SCF to establish the ages of offset deposits and features. We used the single aliquot, regenerative-dose procedure for OSL dating of quartz and feldspar sand grains (Murray and Wintle, 2000, 2003 and radiocarbon dating of detrital charcoal. The ages provide additional constraints on the timing of late-Pleistocene deglaciation in the northern Olympics and the necessary geochronologic context for our interpretation of the late-Pleistocene to Holocene slip rate along the SCF.

OSL Methods and Results
In this study, we collected eight OSL samples from deposits of till, outwash, ice contact, alluvial terraces, and alluvial fans (Figures 2b and 5a, Table 2) to determine the time of deposition of the units given the assumption of complete bleaching of quartz and feldspar grains (Duller, 2004;Rhodes, 2011). Five OSL samples have ages between ∼12.5 and 17.9 ka (Figure 6a, Table 2) and likely represent late glacial or recessional sedimentation (see discussion in Section 5.3).

Radiocarbon Methods and Results
All radiocarbon samples were analyzed using accelerator mass spectrometry by DirectAMS (www.directams. com) and were converted to Cal yrs. B.P. using the program OxCal v 4.3 (Ramsey, 2008(Ramsey, , 2009) and the INT-CAL13 calibration data of Reimer et al. (2013). Radiocarbon ages of detrital charcoal represent a maximum age of the deposit as the charcoal must already be present in the landscape before being mobilized and buried and additionally contain an "inbuilt" age from the age of the wood that was burned (e.g. Gavin, 2001) We dated 11 detrital charcoal samples collected from different deposits along the trace of the SCF (Figure 2b, Table 3). Seven samples were collected from different debris flow deposits related to deposition within the offset channels. Of these, six samples were collected from inside channels and one sample came from a debris flow fan at the base of the channel (Figure 6b and Table 3). Samples were collected as deep in the exposure as possible in attempt to date the oldest deposits within that channel, but were limited by the availability of charcoal, the thickness of deposits, and a shallow water table (which filled many sampling pits). The other four samples were collected from till, an alluvial terrace, and a fan deformed by the SCF. Two of the samples, GB32 and BL02, were sampled from the paleoseismic trenches of Schermer et al., 2020. Most calibrated dates from charcoal samples collected within channel-associated debris flow deposits fall between ∼3,300 and ∼1,400 Cal yrs. BP (Table 3 and Figure 6b). Two older debris flow samples, SCF18-10-06 (5,751-5,982 Cal yrs. BP) and SCF18-11B-01 (10,707-11,104 Cal yrs. BP) may represent recycled charcoal or less active channels which have not eroded away older debris flow deposits. We interpret samples from till (BL02; 31,180-30,780 Cal yrs. BP) and an alluvial fan (SCF18-CD6; 47,812-45,500 Cal yrs. BP) to have been recycled because OSL ages from these deposits are approximately 10,000 and 30,000 years younger, respectively (Figure 6a, Tables 2 and 3).

Geochronologic Interpretation
The OSL and charcoal ages reported in this study are consistent with a 14 ± 1.0 ka deglacial age for the Juan de Fuca Lobe of the Cordilleran ice sheet in the northern Olympic Mountains (Haugerud & Hendy, 2016;Waters et al., 2011). The deglacial age of the Juan de Fuca lobe (14,940-14,095 Cal yrs. BP) suggested by Kovanen and Easterbrook (2001) is similar to OSL ages from outwash, and alluvial fans (12.5-14.5 ka), and to the radiocarbon age of charcoal (14,794 Cal yrs. BP) sampled from coarse grained fluvial deposits (sample GB32, Table 3). This estimate is also consistent with a range of 13.9-12.2 ka from four other radiocarbon samples of a correlative post-glacial fluvial unit to GB32 (Schermer et al., 2020), suggesting that ice had retreated from the study area by 14 ka. OSL samples from ice contact and alluvial fan deposits that date to ∼17.5 ka may represent deposition during periods of minor ice recession during late glacial time or possibly staggered retreat on parts of the Juan de Fuca lobe. The abundance of younger charcoal (∼11 ka-1.5 ka) within debris flow deposits and fluvial terraces (  of landscape evolution in paraglacial settings. For example, sediment supply exhaustion and re-entry of vegetation into the landscape both restrict geomorphic activity to major drainages, small channels on steep hillslopes, and landslides (Ballantyne, 2017).

Fault Slip Rates and Seismic Hazard
In order to estimate the slip rate on the NOFZ, we compare the magnitude of laterally offset channels with radiocarbon ages obtained from channel deposits. If charcoal within channel deposits recorded the age of channel incision (before faulting), then the magnitude of strike-slip displacement would increase with channel age. However, comparison of channel offset with the radiocarbon dates shows no correlation DUCKWORTH ET AL.

10.1029/2020JB020276
12 of 22 between the magnitude of strike-slip displacement and channel deposit age (Figure 6b). Instead, most ages from debris flow deposits within channels are young (between ∼3,300 and ∼1,400 Cal yrs. BP) with two older ages (∼5,800 and ∼11,000 Cal yrs. BP; Figure 6a). We interpret these ages to suggest that channels have continuously deposited and mixed radiocarbon material and debris flow sediment throughout post-glacial time. While geomorphic activity in channels has removed the record of dip-slip displacement within the channels in most cases, these channels still appear to be robust recorders of strike-slip displacement because a wide range of strike-slip magnitudes (including large amounts of slip) are preserved-often in close proximity. For example, Channels B, C, and D are located on the same hillslope, within ∼100 m of each other (Figure 4a), and show preferred offsets of 15.6, 6.5, and 26.0 m, respectively (Table 1). We are confident that these closely spaced channels (and others along the NOFZ) preserve an accurate record of cumulative fault slip because they are spaced farther apart than the total amount of offset, the variance between measurements is >30% of the preferred offset value, we used multiple geomorphic markers to calculate each offset value, and we carefully considered the interseismic geomorphic modification during the measuring processes (see Section 4 and supporting information S3 and S4; Reitman et al., 2019). Therefore, strike-slip displacement may correlate better with channel incision age than the deposit age dated through detrital charcoal. As such, the channels with most strike-slip displacement (∼25 ± 5 m) are inferred to be the oldest channels in the study area.
Channel incision on steep, drift-mantled hillslopes occurs within a few hundred years after glacial retreat due to the large sediment supply and lack of vegetation (Ballantyne, 2017). Therefore, older channels, which record the most strike-slip displacement, likely formed soon after deglaciation, and are the best piercing lines in construction of a post-glacial fault slip rate for the NOFZ. The combination of a deglacial age of 14 ± 1 ka, constrained through existing radiocarbon ages (Haugerud & Hendy, 2016;Kovanen & Easterbrook, 2001) and new radiocarbon and OSL dating presented here (Figure 6a, Table 2) and in Schermer et al. (2020), with the total slip recorded by channels with the most offset (25 ± 5 m) produces a minimum post-glacial strike-slip rate of 1.3-2.3 mm/yr. Similarly, combination of the deglacial age with the range of dip-slip displacements on glacially derived deposits along the NOFZ (0.7-6.5 m) yield a dip-slip rate of 0.05-0.5 mm/yr. If we subtract the amount of slip accumulated in the last earthquake (4 ± 1 m) and instead only consider the amount of slip accumulated over the timespan of the dated earthquake cycles from 12.2 to 10.0 ka to 1.2-0.9 ka (11,900 to 8,800 years; Schermer et al., 2020), the strike slip rate is 1.  The magnitude of earthquakes along the NOFZ is constrained by the rupture length and the amount of displacement in a single surface rupturing event. Schermer et al. (2020) documented four to five earthquakes on the SCF and correlated three to four of these across the SCF, Lake Crescent and LCBCF suggesting that the NOFZ has ruptured along much of the ∼80 km length. Empirical scaling relationships based on fault length (Stirling et al., 2002(Stirling et al., , 2013Wesnousky, 2008) applied to the entire 80 km-long fault system produces a magnitude estimate of M w 7.0-7.5. The smallest dextral offsets measured on the NOFZ ("J Channel," "H Channel," and "OL03 Channel"; Figures 2b and 5, and Table 1) show preferred strike-slip displacements of 4.0 ± 0.8 m, 4.1 ± 1 m, and 4.5 ± 2.4 m. Schermer et al. (2020) interpreted these displacements, particularly "J Channel," to represent the amount of dextral slip that occurred in a single earthquake. Estimates of single-event displacement (4 ± 1 m; Schermer et al., 2020) are consistent with the maximum dextral displacement on the SCF (25 ± 5 m) in four to five surface rupturing earthquakes. Additionally, empirical scaling relations based on single-event displacements (e.g., Hemphill-Haley & Weldon, 1999;Wells & Coppersmith, 1994) yield estimates of Mw 7.2-7.5, consistent with estimates based on the ∼80 km rupture length (M w 7.0-7.5). A ∼80 km long rupture producing ∼4 m of coseismic slip would suggest that the NOFZ is structurally mature (i.e., has a long and developed slip history; Manighetti et al., 2007) consistent with the conclusions of Nelson et al. (2017) that the NOFZ documents a slip reversal from left-lateral in the Miocene to right-lateral in the Holocene. Altogether, our study shows that the NOFZ is the fastest slipping known crustal fault in western Washington and represents a significant seismic hazard to communities and infrastructure in the northern Olympic Peninsula.

Models of Forearc Deformation
To investigate the relationship between slip on crustal faults and subduction zone earthquake cycle processes, we compare fault slip rates and kinematics determined from paleoseismic trenches and geomorphic mapping with numerical models of forearc deformation constrained by GPS velocities. GPS velocity fields within the Cascadia forearc show the influence of crustal block rotation and interseismic coupling along the Cascadia subduction zone McCaffrey et al, 2007McCaffrey et al, , 2013Wells & Simpson, 2001;R. Wells et al., 2014). Quantification of tectonic processes through GPS-constrained models presents the opportunity to understand the influence of subduction zone processes on upper plate faults as done on the Canyon River fault (CRF; Figure 1b) by Delano et al. (2017). These authors compared the late Quaternary slip rate of the CRF with GPS-constrained models of Cascadia subduction zone interseismic coupling to suggest that interseismic stress from a locked subduction zone may drive slip on crustal faults. However, these authors focused on the dip-slip rate of the CRF, despite evidence of left-lateral displacement in paleoseismic trenches (Bennett et al., 2017;Walsh & Logan, 2007), and did not consider the effects of a coseismic stress transfer during a full length subduction zone megathrust event.
We build on the work of Delano et al. (2017) by implementing a boundary element method model (Crouch & Starfield, 1983;Thomas, 1993) to estimate the stress imposed onto the NOFZ by (1) interseismic coupling on the Cascadia subduction zone and (2) a full length Cascadia subduction zone coseismic rupture. Though these models make some simplifying assumptions about rheology, earthquake recurrence, and fault mechanics, we use them to provide a first-order estimate of the slip rate on the crustal faults required to relieve the imposed stress. We assume that the triangular elements representing the crustal faults are shear traction-free surfaces (e.g., Cooke & Dair, 2011) and that the slip completely relieves the accumulated shear stress imposed by the tectonic process being modeled. The resulting slip rate distributions represent the coseismic slip distribution on the crustal faults, normalized by the recurrence interval. The rate and distribution of slip calculated from the modeling are then compared to the measurements of long-term (∼14 ka) fault slip rate on the NOFZ determined in this study and the slip rate on the CRF determined by Delano et al. (2017). Finally, in Section 8, we compare the slip rates and kinematics of the NOFZ and the CRF to existing GPS-constrained elastic block models (McCaffrey et al., 2013;Wells et al., 1998) to evaluate the role these faults play in deformation partitioning across the forearc.

Geologic Constraints on NOFZ and CRF Slip
In this study, we characterize the NOFZ as a dominantly dextral, strike-slip fault with a lateral slip rate of 1.3-2.3 mm/yr and dip slip rate of 0.05-0.5 mm/yr. Bedrock mapping of the northern Olympic Peninsula suggests that the NOFZ has reverse motion (Polenz et al., 2004;Schasse, 2003;Tabor & Cady, 1978a). However, the geomorphic expression of the NOFZ (straight crested scarps, changes in scarp facing direction, etc.) and paleoseismic trenching (Schermer et al., 2020) only suggest that the fault dip is steep and find evidence of both normal and reverse motion. Delano et al. (2017) suggested a south-side up dip-slip rate of ∼0.1-0.3 mm/yr over the past ∼12-37 ka on the CRF based on geomorphic mapping, measurements of scarp profiles, and OSL dating of fluvial terraces. Paleoseismic trenching of south side-up scarps along the CRF show both south dipping, reverse faults (Walsh & Logan, 2007) and north dipping, normal faults (Bennett et al., 2017) suggesting that the dip slip kinematics of the CRF are not well constrained. Additionally, these paleoseismic studies (Bennett et al., 2017;Walsh & Logan, 2007) observed sinistral-oblique slickenlines indicating that the CRF has a left-lateral component that is not well preserved in the geomorphic record .

Boundary Element Method Modeling Process and Parameters
The boundary element method model represents the 3-D geometry of the Cascadia subduction zone, NOFZ, and the CRF using triangular dislocation elements embedded in an elastic half-space characterized by a shear modulus of 3 × 10 10 Pa and Poisson's ratio of 0.25. Within the model, our mapped surface trace of the NOFZ is projected to a depth of 10 km and rotated to a ∼80° dip to the north, consistent with bedrock mapping (Polenz et al., 2004;Schasse, 2003;Tabor & Cady, 1978a). We maintain the same boundary element method model parameters and geometry for the CRF as described by Delano et al. (2017). The geometry of the Cascadia subduction zone is based on slab depth contours derived from seismicity data (McCrory et al., 2012). We prescribe the traction rate or slip rate in the strike and dip direction of each element. On all elements, we prescribe zero slip rate conditions in the element-normal direction to prevent opening or interpenetration across elements. On subduction zone elements, we prescribe the geodetically estimated slip (or slip rate) distribution from the modeled tectonic process in the strike and dip directions. On crustal fault elements, we assign zero shear traction in the strike and dip directions, allowing them to slip freely in response to the tectonic loading applied by the modeled subduction zone process as well as mechanical interactions among the crustal fault segments (Cooke & Dair, 2011). By assuming that crustal fault slip completely relieves the traction imposed by subduction zone processes, we are estimating maximum slip rates.

Crustal Slip Rates Due to Cascadia Subduction Zone Interseismic Coupling
To model the effects of interseismic Cascadia subduction zone coupling on the NOFZ and SCF, we use a slip deficit rate distribution on the subduction interface from a GPS-constrained elastic block model , with microplate boundaries defined by Holocene-active faults (U.S. Geological Survey, 2006) and gradients in GPS velocities (Evans et al., 2015;Loveless & Meade, 2011). In order to isolate the effects of interseismic subduction zone coupling on the NOFZ and CRF, the block model does not include either of these fault systems as block bounding segments when computing the subduction zone slip deficit distribution. The subduction zone slip deficit distribution (Figure 7a) is then input in the boundary element method model along with the NOFZ and SCF geometries and the boundary conditions described above to estimate the slip on the crustal faults necessary to relieve the stress imposed by interseismic coupling ( Figure S8a). We use the same Cascadia subduction zone slip deficit distribution used by Delano et al. (2017) to maintain consistency when comparing model results between the CRF and NOFZ.
The slip rate distribution estimated by the model for the NOFZ and CRF is consistent with dip-slip observations but inconsistent with strike-slip observations from paleoseismic trenches and geomorphic mapping. The model predicts up to 1.6 mm/yr of left lateral slip and up to 0.5 mm/yr of reverse slip on the NOFZ (Figure 8a, Table 4) and up to 0.7 mm/yr of right-lateral slip and 0.1-0.5 mm/yr of reverse slip on the CRF ; Figure 8c, Table 4). The rate of reverse slip on the NOFZ estimated by interseismic coupling is consistent with the magnitude and direction of the post-glacial dip-slip rate, however, the strikeslip component of modeled slip produces the opposite kinematics than shown by the geomorphic record of the NOFZ (Table 4). Similarly, the rate of dip slip on the CRF estimated by the model is consistent with the rate of south-side up slip (0.1-0.3 mm/yr, Table 4) measured by Delano et al. (2017) but shows right-lateral kinematics rather than the left-lateral slickenlines found in paleoseismic trenches (Bennett et al, 2017;Walsh & Logan, 2007).

Crustal Slip Rates Due to Cascadia Subduction Zone Coseismic Stress
We used 33 distributions of coseismic slip originally presented by Frankel et al. (2018) and Wirth et al. (2018) to model the stress imposed by a Cascadia subduction zone megathrust earthquake on the NOFZ and CRF. Each slip distribution is composed of several high stress-drop M w 8.0 subevents superimposed on large, shallower background slip such that the total seismic moment of each model corresponds to M w 9.0 . Coseismic slip distributions are randomly seeded across three prescribed down-dip rupture extents based on the 2014 National Seismic Hazard Maps (Frankel et al., 2015;Petersen et al., 2002). These limits are: (1) the top of the nonvolcanic tremor zone (∼12-25 km depth), (2) the 1 cm/yr locking contour derived from inversion of GPS and uplift data (∼15-30 km), and (3) the midpoint of the edge of the locked zone in thermal modeling and the 1 cm/yr locking location (∼27-35 km). The up-dip rupture limit is specified at a depth of 5 km for all models . We convert each slip distribution to a slip rate distribution assuming a recurrence interval of 500 years (Atwater & Hemphill-Haley, 1997;McCaffrey & Goldfinger, 1995;Priest et al., 2017) and interpolate onto the triangular dislocation elements representing the geometry of the Cascadia subduction zone (McCrory et al., 2012). Two iterations of the model are run for each slip distribution, with variations in subduction zone slip rake (90°, 115°) to account for uncertainty in the rake of coseismic slip (McCaffrey et al., 2013). These slip rate distributions (e.g., Figure 7b  The resulting modeled slip rate distributions on the NOFZ and CRF are consistent with strike-slip observations from paleoseismic trenches and geomorphic mapping, with ∼0.3-1.5 mm/yr of right-lateral and 0.1-0.3 mm/yr of normal slip estimated on the NOFZ (e.g., Table 4, Figure 8b). Accounting for the range in slip rate and motion sense predicted by the boundary element method model for the CRF, we summarize the modeled CRF slip rates as 0.1 mm/yr of right lateral slip to 1.5 mm/yr of left lateral slip and 0.2 mm/yr of reverse slip to 0.8 mm/yr of normal slip (Table 4, Figure 8d). The variation in the magnitude of crustal fault slip rates between model runs is related to the depth, location of greatest slip, and rake of the input subduction zone slip distribution. Generally, the models with deeper down-dip rupture limits, rakes of 90° and/or randomly seeded large magnitudes of slip in the northern Cascadia subduction zone predicted faster crustal slip rates (both dip-slip and strike-slip) relative to models without these attributes. In the scenarios that produce faster crustal strike-slip rates, the magnitude and direction of lateral slip are consistent with the dextral slip we measure on the NOFZ during post-glacial time, and modeled left-lateral motion on the CRF is consistent with paleoseismic trenching (Bennett et al., 2017;Walsh & Logan, 2007). These models, however, predict normal slip on the crustal faults, inconsistent with the interpretation that Olympic Peninsula faults have reverse motion (Polenz et al., 2004;Schasse, 2003;Tabor & Cady, 1978a).

Interpretation of Forearc Deformation Patterns
Comparison of long-term (10 3 yr) with model-predicted short term rates (10 1 yr) and kinematics of deformation on the Olympic Peninsula is complicated by ambiguity in long-term crustal fault kinematics. While bedrock maps indicate that the dip-slip component of the NOFZ is reverse (Polenz et al., 2004;Schasse, 2003;Tabor & Cady, 1978a), we here characterize the NOFZ to be a dominantly dextral strike slip fault, as we did not find any direct geomorphic evidence to constrain a primary dip direction or the kinematics of the dip-slip component over post-glacial time. Model-predicted strike slip on the NOFZ in response to stress imparted by coseismic slip on the underlying subduction zone is consistent in sense and magnitude with our geologic constraints.
The consistency between measured slip rates and kinematics on the NOFZ and those estimated by models of stress transfer from coseismic slip on the subduction zone does not necessarily indicate that a megathrust earthquake will trigger an upper plate earthquake (e.g., Gomberg & Sherrod, 2014), as the boundary element method model only considers the magnitude of stress changes due to subduction zone earthquakes, not details in the timing of crustal earthquake occurrence. We suggest that crustal faults will accumulate stress imposed by regional tectonic processes, occurring over multiple megathrust earthquake cycles, and will slip once their strength is exceeded by the cumulative stress level, which may or may not be coincident with a megathrust earthquake (e.g., Loveless et al., 2010).
An understanding of the stresses that drive slip on upper plate faults requires consideration of stresses within the Cascadia forearc apart from those of the subduction zone earthquake cycle.  the earthquake cycle and microplate motion have commonly been modeled using GPS-constrained elastic block models. Block models of the Pacific Northwest predict shortening across block-bounding faults that cross the Olympic Peninsula as a result of clockwise block rotation relative to stable North America (McCaffrey et al., 2013;Savage & Wells, 2015;Wells et al., 1998). McCaffrey et al. (2013) present one such blocktheir "Olym" block-which accommodates shortening through ∼0-1 mm/yr of reverse and ∼1-2 mm/yr of right-lateral slip along its northern boundary, which is similar in location and strike to the mapped trace of the NOFZ. Similarly, the "Olym" block is bounded on the southern side of the Olympic Peninsula by a left-lateral (∼0-1 mm/yr), reverse (∼1-2 mm/yr) fault near the CRF. Though other mapped crustal faults on the Olympic Peninsula could be considered as alternative or additional block boundaries, we suggest that the proximity, similar fault geometry, and consistent slip rates/kinematics between geologic and geomorphic observations of the NOFZ and CRF with those modeled for the northern and southern boundaries of the "Olym" block of McCaffrey et al. (2013) imply that the NOFZ and CRF may be treated as microplate boundaries within elastic block models of the Cascadia forearc and that long-term slip patterns along the NOFZ are consistent with the modern, nominally interseismic GPS velocity field.
Interpreted together, results from boundary element method modeling and published block models suggest that right-lateral slip on the NOFZ requires subduction zone coseismic stress transfer and/or that this fault system acts as microplate boundary among a broader deformation field that extends beyond the subduction zone. While elastic block models are constrained by nominally interseismic GPS velocities, the fault slip rates that are estimated can be taken to represent long-term relative block motion across a complete seismic cycle (McCaffrey, 2002;). Block models do not, however, make assumptions about how an episode of slip on a given fault may be triggered by stress transfer from other structures. We therefore, interpret the boundary element and block model results as consistent: the NOFZ fault slip estimated in the block model could be achieved, at least in part, by stress imposed on it by subduction zone earthquakes. Interaction between the subduction zone and the network of crustal fault zones-of which the NOFZ and CRF are but two-indicates a complex partitioning of fore-arc deformation. During any single phase of the subduction zone earthquake cycle, crustal faults may be activated, introducing complexity into the simple concept of a perfectly elastic earthquake cycle.

Conclusions
Based on geomorphic mapping, fault slip analysis, and dating of offset geomorphic features and deposits along the SCF, we characterize a post-glacial chronology of geomorphic activity and constrain the long term (10 3 -10 4 years) rate and kinematics of faulting. Radiocarbon and OSL results are consistent with previous work and suggest that the northern Olympic Mountains were deglaciated by ∼14 ka. Immediately following the retreat of the Juan de Fuca lobe of the Cordilleran ice sheet, geomorphic instability in the landscape mobilized sediment, deposited fans and terraces, and carved debris flow channels into the steep till-mantled hillslopes. Subsequent stabilization of the hillslopes restricted Holocene geomorphic activity to channel incision, debris flow deposition, and landslides. Measurements of slip, derived from laterally offset channels and scarp profiles, suggest that the SCF and LCBCF should be considered as one continuous fault system, the NOFZ, characterized by dominantly dextral motion (∼6:1 strike slip to dip slip ratio) along a steeply dipping fault plane, and a slip rate of 1.3-2.3 mm/yr.
The similarities between long-term and short-term deformation from geodetically constrained elastic block models and models of Cascadia subduction zone coseismic stress transfer suggest that the NOFZ plays an important role in the accumulation of permanent strain within the forearc. Slip rates arising from forearc block motions are consistent with slip rates and kinematics on the NOFZ, and boundary element method models suggest that distinct slip events on these crustal faults may be modulated by stresses related to subduction zone earthquakes.